freesteam: solver2.h Source File

00001 /* 00002 00003 freesteam - IAPWS-IF97 steam tables library 00004 Copyright (C) 2004-2005 John Pye 00005 00006 This program is free software; you can redistribute it and/or 00007 modify it under the terms of the GNU General Public License 00008 as published by the Free Software Foundation; either version 2 00009 of the License, or (at your option) any later version. 00010 00011 This program is distributed in the hope that it will be useful, 00012 but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 GNU General Public License for more details. 00015 00016 You should have received a copy of the GNU General Public License 00017 along with this program; if not, write to the Free Software 00018 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 00019 00020 */ 00021 00022 #ifndef SOLVER2_H 00023 #define SOLVER2_H 00024 00025 #define SOLVER2_DEBUG 00026 #define SOLVER2BASE_DEBUG Solver2Base<FirstProp,SecondProp,FirstPropAlt,SecondPropAlt>::debug 00027 00028 #include "steamproperty.h" 00029 #include "units.h" 00030 #include "satcurve.h" 00031 #include "exception.h" 00032 #include "convergencetest.h" 00033 #include "b23curve.h" 00034 #include "boundaries.h" 00035 00037 00045 template<class FirstProp,class SecondProp,int FirstPropAlt=0, int SecondPropAlt=0> 00046 class Solver2Base{ 00047 00048 protected: 00049 00050 #ifdef SOLVER2_DEBUG 00051 Solver2Base(const bool debug=false){ 00052 this->debug=debug; 00053 } 00054 #else 00055 Solver2Base(){} 00056 #endif 00057 00058 virtual ~Solver2Base(){} 00059 00060 virtual int whichRegion(const FirstProp &fp,const SecondProp &sp) = 0; 00061 00062 virtual SteamCalculator solve(const FirstProp &fp, const SecondProp &sp) = 0; 00063 00064 FirstProp getFirstProp(SteamCalculator &S){ 00065 return SteamProperty<FirstProp,FirstPropAlt>::get(S); 00066 } 00067 00068 SecondProp getSecondProp(SteamCalculator &S){ 00069 return SteamProperty<SecondProp,SecondPropAlt>::get(S); 00070 } 00071 00072 #ifdef SOLVER2_DEBUG 00073 bool debug; 00074 #endif 00075 00076 }; 00077 00079 00106 template<class FirstProp,class SecondProp,int FirstPropAlt=0, int SecondPropAlt=0> 00107 class Solver2 00108 : public Solver2Base<FirstProp,SecondProp,FirstPropAlt,SecondPropAlt>{ 00109 00110 private: 00111 static const int maxiter=30; 00112 00114 00117 SteamCalculator makeRegion1Guess(const FirstProp& f, const SecondProp &s){ 00118 SteamCalculator S; 00119 S.set_pT(10.0 * bar, fromcelsius(100)); 00120 return S; 00121 } 00122 00123 00125 00130 SteamCalculator makeRegion2Guess(const FirstProp& f, const SecondProp &s){ 00131 SteamCalculator S; 00132 S.set_pT(6.0 * MPa, fromcelsius(400)); 00133 return S; 00134 } 00135 00137 00142 SteamCalculator makeRegion3Guess(const FirstProp& f, const SecondProp &s){ 00143 SteamCalculator S; 00144 S.setRegion3_rhoT(1 / 0.00317 * kg_m3, fromcelsius(400)); 00145 return S; 00146 } 00147 00149 00154 SteamCalculator makeRegion4Guess(const FirstProp& f, const SecondProp &s){ 00155 SteamCalculator S; 00156 S.setRegion4_Tx(fromcelsius(263.9),0.5); 00157 return S; 00158 } 00159 00160 public: 00161 00162 #ifdef SOLVER2_DEBUG 00163 Solver2(const bool debug=false) 00164 : Solver2Base<FirstProp,SecondProp,FirstPropAlt,SecondPropAlt>(debug){ 00165 #else 00166 Solver2() 00167 : Solver2Base<FirstProp,SecondProp,FirstPropAlt,SecondPropAlt>(){ 00168 #endif 00169 //cerr << endl <<"Solver2<" << SteamProperty<FirstProp,FirstPropAlt>::name() << "," << SteamProperty<SecondProp,SecondPropAlt>::name() << ">::Solver2()"; 00170 } 00171 00173 00186 int whichRegion(const FirstProp &fp, const SecondProp &sp); 00187 00189 00195 SteamCalculator solve(const FirstProp &fp, const SecondProp &sp){ 00196 int region = 0; 00197 SteamCalculator S,S2; 00198 00199 try{ 00200 //cerr << endl << "Solver2: solving for " << SteamProperty<FirstProp,FirstPropAlt>::name() << " = " << fp << ", " << SteamProperty<SecondProp,SecondPropAlt>::name() << " = " << sp; 00201 00202 region = whichRegion(fp,sp); 00203 00204 }catch(Exception *E){ 00205 stringstream ss; 00206 ss << "Solver2<" << SteamProperty<FirstProp,FirstPropAlt>::name() << "," << SteamProperty<SecondProp,SecondPropAlt>::name() << ">::solve (no first guess): " << E->what(); 00207 delete E; 00208 throw new Exception(ss.str()); 00209 } 00210 00211 try{ 00212 //cerr << endl << "Solver2: solving in region " << region; 00213 switch(region){ 00214 case 1: 00215 S2 = solveRegion1(fp,sp,makeRegion1Guess(fp,sp)); 00216 break; 00217 case 2: 00218 S2 = solveRegion2(fp,sp,makeRegion2Guess(fp,sp)); 00219 break; 00220 case 3: 00221 S2 = solveRegion3(fp,sp,makeRegion3Guess(fp,sp)); 00222 break; 00223 case 4: 00224 S2 = solveRegion4(fp,sp,makeRegion4Guess(fp,sp)); 00225 break; 00226 default: 00227 stringstream ss; 00228 ss << "Invalid return from whichRegion(" << fp << "," << sp << "): region = " << region; 00229 throw new Exception(ss.str()); 00230 } 00231 }catch(Exception *E){ 00232 stringstream s; 00233 s << "Solver2<" << SteamProperty<FirstProp,FirstPropAlt>::name() << "," << SteamProperty<SecondProp,SecondPropAlt>::name() << ">::solve (no first guess; found region=" << region << "): " << E->what(); 00234 delete E; 00235 throw new Exception(s.str()); 00236 } 00237 00238 return S2; 00239 } 00240 00242 00249 SteamCalculator solve(const FirstProp &fp, const SecondProp &sp, const SteamCalculator &firstguess){ 00250 int region=0; 00251 try{ 00252 region = whichRegion(fp,sp); 00253 00254 if(!firstguess.isSet() || firstguess.whichRegion()!=region){ 00255 return solve(fp,sp); 00256 } 00257 00258 switch(region){ 00259 case 1: 00260 return solveRegion1(fp,sp,firstguess); 00261 case 2: 00262 return solveRegion2(fp,sp,firstguess); 00263 case 3: 00264 return solveRegion3(fp,sp,firstguess); 00265 case 4: 00266 return solveRegion4(fp,sp,firstguess); 00267 } 00268 00269 }catch(Exception *E){ 00270 stringstream ss; 00271 00272 ss << "Solver2<" << SteamProperty<FirstProp,FirstPropAlt>::name() << "," << SteamProperty<SecondProp,SecondPropAlt>::name() << ">::solve (" << SteamProperty<FirstProp,FirstPropAlt>::name() << "=" << fp << ", " << SteamProperty<SecondProp,SecondPropAlt>::name() << "=" << sp; 00273 00274 ss << "; found region="<< region; 00275 00276 if(firstguess.isSet()){ 00277 ss << "; given first guess with " << SteamProperty<FirstProp,FirstPropAlt>::name() << "=" << SteamProperty<FirstProp,FirstPropAlt>::get(firstguess) << ", " << SteamProperty<SecondProp,SecondPropAlt>::name() << "=" << SteamProperty<SecondProp,SecondPropAlt>::get(firstguess); 00278 }else{ 00279 ss << "; given uninitialised first guess"; 00280 } 00281 00282 ss << "): " << E->what(); 00283 00284 delete E; 00285 throw new Exception(ss.str()); 00286 } 00287 } 00288 00289 private: 00290 00292 SteamCalculator solveRegion3(const FirstProp &f1, const SecondProp &s1,const SteamCalculator &firstguess){ 00293 try{ 00294 // Iterate on rho, T 00295 00296 // Just like region 1, except for it's T,x 00297 00298 //cerr << endl << "Solver2<" << SteamProperty<FirstProp,FirstPropAlt>::name() << "," << SteamProperty<SecondProp,SecondPropAlt>::name() << ">::solveRegion3:"; 00299 00300 SteamCalculator guess = firstguess; 00301 00302 //cerr << endl << "Solver2::solveRegion4: firstguess copied to guess OK"; 00303 00304 if(firstguess.whichRegion()!=3){ 00305 throw new Exception("Solver2::solveRegion3: First guess is not region3"); 00306 } 00307 00308 SteamCalculator petT, petrho; 00309 Temperature T,dT; 00310 Density rho,drho; 00311 FirstProp f,Df, DfT, Dfrho; 00312 SecondProp s,Ds,DsT, Dsrho; 00313 Pressure p; 00314 00315 // cerr << endl << "Solver2::solveRegion3: Starting iterations"; 00316 00317 00318 int niter=0; 00319 while(1){ 00320 00321 T = guess.temp(); 00322 rho = guess.dens(); 00323 p = guess.pres(); 00324 00325 //cerr << endl << "Iter " << niter << ": T = " << T << ", rho = " << rho; 00326 00327 f = SteamProperty<FirstProp,FirstPropAlt>::get(guess); 00328 s = SteamProperty<SecondProp,SecondPropAlt>::get(guess); 00329 00330 //cerr << ": " << SteamProperty<FirstProp,FirstPropAlt>::name() << " = " << f; 00331 //cerr << ", " << SteamProperty<SecondProp,SecondPropAlt>::name() << " = " << s; 00332 00333 Df = f1 - f; 00334 Ds = s1 - s; 00335 00336 // In this template it's hard to know the units of the convergence test, so make it as a new, separate, template: 00337 if( 00338 ConvergenceTest<FirstProp,FirstPropAlt>::test(Df,p,T) 00339 && ConvergenceTest<SecondProp,SecondPropAlt>::test(Ds,p,T) 00340 ){ 00341 // cerr << endl << " ... SOLUTION OK (" << niter << " iterations)" << endl; 00342 return guess; 00343 } 00344 00345 dT = T * 0.001; 00346 00347 if(T + dT > T_MAX){ 00348 dT = -dT; 00349 } 00350 00351 petT.setRegion3_rhoT(rho, T + dT); 00352 ASSERT(petT.whichRegion()==3); 00353 00354 // ensure we don't go over rho limits 00355 drho = rho * 0.001; 00356 if(rho + drho > 50.0 * kg_m3){ 00357 drho = -drho; 00358 } 00359 00360 petrho.setRegion3_rhoT(rho + drho, T); 00361 ASSERT(petrho.whichRegion()==3); 00362 00363 DfT = SteamProperty<FirstProp,FirstPropAlt>::get(petT) - f; 00364 DsT = SteamProperty<SecondProp,SecondPropAlt>::get(petT) - s; 00365 00366 Dfrho = SteamProperty<FirstProp,FirstPropAlt>::get(petrho) - f; 00367 Dsrho = SteamProperty<SecondProp,SecondPropAlt>::get(petrho) - s; 00368 00369 ASSERT(DfT * Dsrho != DsT * Dfrho); 00370 00371 // Solve new peturbations 00372 dT = dT * (Df*Dsrho - Ds*Dfrho) / (DfT*Dsrho - DsT*Dfrho); 00373 drho = drho * (DfT*Ds - DsT*Df) / (DfT*Dsrho - DsT*Dfrho); 00374 00375 ASSERT(!isnan(dT)); 00376 ASSERT(!isnan(drho)); 00377 00378 // Regulate maximum change in temperature 00379 00380 Temperature dTMax = 0.1 * T; 00381 Temperature dTAbs = fabs(dT); 00382 if(dTAbs > dTMax){ 00383 //cerr << endl << " ... limiting dT, too great"; 00384 dT = dT * dTMax/dTAbs; 00385 } 00386 00387 // Regulate max change in pressure 00388 00389 Density drhoMax = 0.4 * rho; 00390 Density drhoAbs = fabs(drho); 00391 if(drhoAbs > drhoMax){ 00392 //cerr << endl << " ... limiting drho, too great"; 00393 drho = drho * drhoMax/drhoAbs; 00394 } 00395 00396 //cerr << endl << " ... calculated dT = " << dT << ", drho = " << drho; 00397 T = T + dT; 00398 rho = rho + drho; 00399 00400 if(T > T_MAX){ 00401 // Strange, this test never seems to be used... 00402 00403 //cerr << endl << " .... Applying T_MAX limit"; 00404 T = T_MAX; 00405 } 00406 00407 if(T < T_REG1_REG3){ 00408 //cerr << endl << " .... Applying T_REG1_REG3 limit"; 00409 T = T_REG1_REG3; 00410 } 00411 00412 Density rho_b134; 00413 SteamCalculator S; 00414 S.setSatWater_T(T_REG1_REG3); 00415 rho_b134 = S.dens(); 00416 if(rho < rho_b134 && T < T_CRIT){ 00417 if(rho < RHO_CRIT){ 00418 Density rho_g = Boundaries::getSatDensSteam_T(T); 00419 if(rho > rho_g){ 00420 //cerr << endl << " ... Applying rho_g limit"; 00421 rho = rho_g; 00422 } 00423 }else if(rho > RHO_CRIT){ 00424 Density rho_f = Boundaries::getSatDensWater_T(T); 00425 if(rho < rho_f){ 00426 //cerr << endl << " ... Applying rho_f limit"; 00427 rho = rho_f; 00428 } 00429 } 00430 00431 } 00432 00433 B23Curve<Density,Temperature> B23; 00434 Density rho_b23 = B23.solve(T); 00435 if(rho < rho_b23){ 00436 rho = rho_b23; 00437 } 00438 00439 S.setRegion3_rhoT(rho,T); 00440 int pmaxiter = 0; 00441 if(S.pres() > P_MAX){ 00442 do{ 00443 //cerr << endl << " ... found p = " << S.pres() / MPa << " MPa"; 00444 00445 drho *= 0.5001; 00446 //dT *= 0.5001; 00447 rho -= drho; 00448 //T -= dT; 00449 S.setRegion3_rhoT(rho,T); 00450 //cerr << endl << " ... Applying P_MAX limit at T = " << T << ": rho = " << rho; 00451 if(++pmaxiter > 20){ 00452 throw new Exception("Solver2::solveRegion3: Failed to find rho of P_MAX"); 00453 } 00454 00455 }while(S.pres() > P_MAX); 00456 //cerr << endl << " ... pressure capped to " << S.pres(); 00457 } 00458 00459 guess.setRegion3_rhoT(rho,T); 00460 00461 niter++; 00462 00463 if(niter > maxiter){ 00464 throw new Exception("Solver2::solveRegion3: Exceeded maximum iterations"); 00465 } 00466 } 00467 00468 }catch(Exception *E){ 00469 stringstream s; 00470 s << "Solver2<" << SteamProperty<FirstProp,FirstPropAlt>::name() << "," << SteamProperty<SecondProp,SecondPropAlt>::name() << ">:solveRegion3: " << E->what(); 00471 delete E; 00472 throw new Exception(s.str()); 00473 } 00474 } 00475 00477 SteamCalculator solveRegion4(const FirstProp &f1, const SecondProp &s1,const SteamCalculator &firstguess){ 00478 00479 try{ 00480 // Just like region 1, except for it's T,x 00481 00482 //cerr << endl << "Solver2::solveRegion4: "; 00483 00484 SteamCalculator guess = firstguess; 00485 00486 //cerr << endl << "Solver2::solveRegion4: firstguess copied to guess OK"; 00487 00488 if(firstguess.whichRegion()!=4){ 00489 throw new Exception("Solver2::solveRegion4: First guess is not region 4"); 00490 } 00491 00492 SteamCalculator petT, petx; 00493 Temperature T,dT; 00494 Num x,dx; 00495 FirstProp f,Df, DfT, Dfx; 00496 SecondProp s,Ds,DsT, Dsx; 00497 Pressure p; 00498 00499 //cerr << endl << "Solver2::solveRegion4: Starting iterations"; 00500 00501 int niter=0; 00502 while(1){ 00503 00504 T = guess.temp(); 00505 x = guess.quality(); 00506 p = guess.pres(); 00507 00508 //cerr << endl << "Iter " << niter << ": x = " << x << ", T = " << T; 00509 00510 f = SteamProperty<FirstProp,FirstPropAlt>::get(guess); 00511 s = SteamProperty<SecondProp,SecondPropAlt>::get(guess); 00512 00513 //cerr << ": " << SteamProperty<FirstProp,FirstPropAlt>::name() << " = " << f; 00514 //cerr << ", " << SteamProperty<SecondProp,SecondPropAlt>::name() << " = " << s; 00515 00516 Df = f1 - f; 00517 Ds = s1 - s; 00518 00519 // In this template it's hard to know the units of the convergence test, so make it as a new, separate, template: 00520 if( 00521 ConvergenceTest<FirstProp,FirstPropAlt>::test(Df,p,T) 00522 && ConvergenceTest<SecondProp,SecondPropAlt>::test(Ds,p,T) 00523 ){ 00524 //cerr << endl << " ... SOLUTION OK" << endl; 00525 return guess; 00526 } 00527 00528 dT = -T * 0.001; 00529 // ensure we don't go over the T limits with our peturbation 00530 if(T + dT < T_TRIPLE){ 00531 dT = -dT; 00532 } 00533 00534 Boundaries::isSat_Tx(T+dT,x, true); 00535 00536 petT.setRegion4_Tx(T + dT, x); 00537 ASSERT(petT.whichRegion()==4); 00538 00539 // ensure we don't go over x limits 00540 dx = 0.001; 00541 if(x + dx > 1.0){ 00542 dx = -dx; 00543 } 00544 00545 Boundaries::isSat_Tx(T, x+dx, true); 00546 00547 petx.setRegion4_Tx(T, x + dx); 00548 ASSERT(petx.whichRegion()==4); 00549 00550 DfT = SteamProperty<FirstProp,FirstPropAlt>::get(petT) - f; 00551 DsT = SteamProperty<SecondProp,SecondPropAlt>::get(petT) - s; 00552 00553 Dfx = SteamProperty<FirstProp,FirstPropAlt>::get(petx) - f; 00554 Dsx = SteamProperty<SecondProp,SecondPropAlt>::get(petx) - s; 00555 00556 //ASSERT(!(DfT == 0.0 * FirstProp() && Dfx == 0.0 * FirstProp())); 00557 //ASSERT(!(DsT == 0.0 * SecondProp() && Dsx == 0.0 * SecondProp())); 00558 //ASSERT(!(DsT == 0.0 * SecondProp() && DfT == 0.0 * FirstProp())); 00559 //ASSERT(!(Dsx == 0.0 * SecondProp() && Dfx == 0.0 * FirstProp())); 00560 ASSERT(DfT * Dsx != DsT * Dfx); 00561 00562 //ASSERT(0.0 * FirstProp() * SecondProp() != DfT * Dsx - DsT *Dfx); 00563 00564 //ASSERT(!((DfT==0.0 * FirstProp() || Dsx==0.0 * SecondProp()) && (Dfx==0.0 * FirstProp() || DsT==0.0 * SecondProp()))); 00565 00566 00567 // Solve new peturbations 00568 dT = dT * (Df*Dsx - Ds*Dfx) / (DfT*Dsx - DsT*Dfx); 00569 dx = dx * (DfT*Ds - DsT*Df) / (DfT*Dsx - DsT*Dfx); 00570 00571 ASSERT(!isnan(dT)); 00572 ASSERT(!isnan(dx)); 00573 00574 // Regulate maximum change in temperature 00575 00576 Temperature dTMax = 0.1 * T; 00577 Temperature dTAbs = fabs(dT); 00578 if(dTAbs > dTMax){ 00579 //cerr << endl << " ... limiting dT, too great"; 00580 dT = dT * dTMax/dTAbs; 00581 } 00582 00583 // Regulate max change in pressure 00584 00585 Num dxMax = 0.2; 00586 Num dxAbs = fabs(dx); 00587 if(dxAbs > dxMax){ 00588 //cerr << endl << " ... limiting dx, too great"; 00589 dx = dx * dxMax/dxAbs; 00590 } 00591 00592 //cerr << endl << " ... calculated dT = " << dT << ", dx = " << dx; 00593 T = T + dT; 00594 x = x + dx; 00595 00596 // Keep new temperature in limits 00597 if(T > T_CRIT){ 00598 //cerr << endl << " ... applying T_CRIT limit"; 00599 T = T_CRIT; 00600 }else if(T < T_TRIPLE){ 00601 //cerr << endl << " ... applying T_TRIPLE limit"; 00602 T = T_TRIPLE; 00603 } 00604 00605 // Keep pressure in limits 00606 if(x > 1){ 00607 //cerr << endl << " ... applying upper x limit"; 00608 x = 1; 00609 }else if(x < 0){ 00610 //cerr << endl << " ... applying lower x limit"; 00611 x = 0; 00612 } 00613 00614 Boundaries::isSat_Tx(T, x, true); 00615 guess.setRegion4_Tx(T,x); 00616 00617 niter++; 00618 00619 if(niter > maxiter){ 00620 throw new Exception("Solver2::solveRegion4: Exceeded maximum iterations"); 00621 } 00622 } 00623 }catch(Exception *E){ 00624 stringstream s; 00625 s<< "Solver2::solveRegion4: " << E->what(); 00626 delete E; 00627 throw new Exception(s.str()); 00628 } 00629 } 00630 00632 SteamCalculator solveRegion1(const FirstProp &f1, const SecondProp &s1,const SteamCalculator &firstguess){ 00633 00634 SteamCalculator guess = firstguess; 00635 00636 if(firstguess.whichRegion()!=1){ 00637 throw new Exception("Solver2::solveRegion1: First guess is not region 1"); 00638 } 00639 00640 SteamCalculator petT, petp; 00641 Temperature T,dT; 00642 Pressure p,dp; 00643 FirstProp f,Df, DfT, Dfp; 00644 SecondProp s,Ds,DsT, Dsp; 00645 00646 int niter=0; 00647 // If we are in region 1, then we will be iterating with pressure and temperature 00648 00649 while(1){ 00650 00651 T = guess.temp(); 00652 p = guess.pres(); 00653 00654 //cerr << endl << "Iter " << niter << ": p = " << p << ", T = " << T; 00655 00656 f = SteamProperty<FirstProp,FirstPropAlt>::get(guess); 00657 s = SteamProperty<SecondProp,SecondPropAlt>::get(guess); 00658 00659 //cerr << ": " << SteamProperty<FirstProp,FirstPropAlt>::name() << " = " << f; 00660 //cerr << ", " << SteamProperty<SecondProp,SecondPropAlt>::name() << " = " << s; 00661 00662 Df = f1 - f; 00663 Ds = s1 - s; 00664 00665 // In this template it's hard to know the units of the convergence test, so make it as a new, separate, template: 00666 if( 00667 ConvergenceTest<FirstProp,FirstPropAlt>::test(Df,p,T) 00668 && ConvergenceTest<SecondProp,SecondPropAlt>::test(Ds,p,T) 00669 ){ 00670 //cerr << endl << " ... SOLUTION OK" << endl; 00671 return guess; 00672 } 00673 00674 dT = -T * 0.001; 00675 // ensure we don't go over the T limits with our peturbation 00676 if(T + dT < T_TRIPLE){ 00677 dT = -dT; 00678 } 00679 00680 petT.set_pT(p,T + dT); 00681 ASSERT(petT.whichRegion()==1); 00682 00683 // ensure we don't go over p limits 00684 dp = p * 0.001; 00685 if(p > 99.0 * MPa){ 00686 dp = -dp; 00687 } 00688 00689 petp.set_pT(p + dp, T); 00690 ASSERT(petp.whichRegion()==1); 00691 00692 DfT = SteamProperty<FirstProp,FirstPropAlt>::get(petT) - f; 00693 DsT = SteamProperty<SecondProp,SecondPropAlt>::get(petT) - s; 00694 00695 Dfp = SteamProperty<FirstProp,FirstPropAlt>::get(petp) - f; 00696 Dsp = SteamProperty<SecondProp,SecondPropAlt>::get(petp) - s; 00697 00698 // Solve new peturbations 00699 dT = dT * (Df*Dsp - Ds*Dfp) / (DfT*Dsp - DsT*Dfp); 00700 dp = dp * (DfT*Ds - DsT*Df) / (DfT*Dsp - DsT*Dfp); 00701 00702 // Regulate maximum change in temperature 00703 00704 Temperature dTMax = 0.1 * T; 00705 Temperature dTAbs = fabs(dT); 00706 if(dTAbs > dTMax){ 00707 //cerr << endl << " ... limiting dT, too great"; 00708 dT = dT * dTMax/dTAbs; 00709 } 00710 00711 // Regulate max change in pressure 00712 00713 Pressure dpMax = 0.4 * p; 00714 Pressure dpAbs = fabs(dp); 00715 if(dpAbs > dpMax){ 00716 //cerr << endl << " ... limiting dp, too great"; 00717 dp = dp * dpMax/dpAbs; 00718 } 00719 00720 //cerr << endl << " ... calculated dT = " << dT << ", dp = " << dp; 00721 T = T + dT; 00722 p = p + dp; 00723 00724 // Keep new temperature in limits 00725 if(T > T_REG1_REG3){ 00726 //cerr << endl << " ... applying upper T limit"; 00727 T = T_REG1_REG3; 00728 }else if(T < T_TRIPLE){ 00729 //cerr << endl << " ... applying lower T limit"; 00730 T = T_TRIPLE; 00731 } 00732 00733 // Keep pressure in limits 00734 if(p > P_MAX){ 00735 //cerr << endl << " ... applying upper p limit"; 00736 p = P_MAX; 00737 }else if(p < P_TRIPLE){ 00738 //cerr << endl << " ... applying lower p_triple limit"; 00739 p = P_TRIPLE; 00740 }else if(p < PB_LOW){ 00741 Pressure p_sat = Boundaries::getSatPres_T(T); 00742 if(p < p_sat){ 00743 //cerr << endl << " ... applying saturation limit to p"; 00744 p = p_sat; 00745 } 00746 } 00747 00748 ASSERT(guess.whichRegion()==1); 00749 ASSERT(Boundaries::isRegion1_pT(p,T)); 00750 00751 guess.setRegion1_pT(p,T); 00752 00753 niter++; 00754 00755 if(niter > maxiter){ 00756 throw new Exception("Solver2::solveRegion1: Exceeded maximum iterations"); 00757 } 00758 } 00759 } 00760 00762 SteamCalculator solveRegion2(const FirstProp &f1, const SecondProp &s1,const SteamCalculator &firstguess){ 00763 00764 // Most of this is the same as for solveRegion1: 00765 00766 #ifdef SOLVER2_DEBUG 00767 if(SOLVER2BASE_DEBUG){ 00768 cerr << endl << "---------------------------------" << endl << "SOLVE REGION 2"; 00769 } 00770 #endif 00771 00772 SteamCalculator guess = firstguess; 00773 00774 if(firstguess.whichRegion()!=2){ 00775 throw new Exception("Solver2::solveRegion2: First guess is not region 2"); 00776 } 00777 00778 SteamCalculator petT, petp; 00779 Temperature T,dT; 00780 Pressure p,dp; 00781 FirstProp f,Df, DfT, Dfp; 00782 SecondProp s,Ds,DsT, Dsp; 00783 00784 int niter=0; 00785 // If we are in region 1, then we will be iterating with pressure and temperature 00786 00787 while(1){ 00788 00789 T = guess.temp(); 00790 p = guess.pres(); 00791 00792 #ifdef SOLVER2_DEBUG 00793 if(SOLVER2BASE_DEBUG){ 00794 cerr << endl << "Iter " << niter << ": p = " << p / MPa << " MPa, T = " << T << " (" << tocelsius(T) << "°C)"; 00795 } 00796 #endif 00797 00798 f = SteamProperty<FirstProp,FirstPropAlt>::get(guess); 00799 s = SteamProperty<SecondProp,SecondPropAlt>::get(guess); 00800 00801 #ifdef SOLVER2_DEBUG 00802 if(SOLVER2BASE_DEBUG){ 00803 cerr << ": " << SteamProperty<FirstProp,FirstPropAlt>::name() << " = " << f; 00804 cerr << ", " << SteamProperty<SecondProp,SecondPropAlt>::name() << " = " << s; 00805 } 00806 #endif 00807 00808 Df = f1 - f; 00809 Ds = s1 - s; 00810 00811 if( 00812 ConvergenceTest<FirstProp,FirstPropAlt>::test(Df,p,T) 00813 && ConvergenceTest<SecondProp,SecondPropAlt>::test(Ds,p,T) 00814 ){ 00815 //cerr << endl << " ... SOLUTION OK" << endl; 00816 return guess; 00817 } 00818 00819 // Peturb T but keep inside T_MAX 00820 dT = T * 0.001; 00821 if(T + dT > T_MAX){ 00822 dT = -dT; 00823 } 00824 00825 petT.set_pT(p,T+ dT); 00826 ASSERT(petT.whichRegion()==2); 00827 00828 // Peturb p but keep above P_TRIPLE 00829 dp = -p * 0.001; 00830 if(p + dp < P_TRIPLE){ 00831 dp = -dp; 00832 } 00833 petp.set_pT(p + dp, T); 00834 ASSERT(petp.whichRegion()==2); 00835 00836 DfT = SteamProperty<FirstProp,FirstPropAlt>::get(petT) - f; 00837 DsT = SteamProperty<SecondProp,SecondPropAlt>::get(petT) - s; 00838 00839 Dfp = SteamProperty<FirstProp,FirstPropAlt>::get(petp) - f; 00840 Dsp = SteamProperty<SecondProp,SecondPropAlt>::get(petp) - s; 00841 00842 // Solve new peturbations 00843 dT = dT * (Df*Dsp - Ds*Dfp) / (DfT*Dsp - DsT*Dfp); 00844 dp = dp * (DfT*Ds - DsT*Df) / (DfT*Dsp - DsT*Dfp); 00845 00846 // Regulate maximum change in temperature 00847 00848 Temperature dTMax = 0.2 * T; 00849 Temperature dTAbs = fabs(dT); 00850 if(dTAbs > dTMax){ 00851 //cerr << endl << " ... limiting dT, too great"; 00852 dT = dT * dTMax/dTAbs; 00853 } 00854 00855 // Regulate max change in pressure 00856 00857 Pressure dpMax = 0.5 * p; 00858 Pressure dpAbs = fabs(dp); 00859 if(dpAbs > dpMax){ 00860 //cerr << endl << " ... limiting dp, too great"; 00861 dp = dp * dpMax/dpAbs; 00862 } 00863 00864 //cerr << endl << " ... calculated dT = " << dT << ", dp = " << dp; 00865 T = T + dT; 00866 p = p + dp; 00867 00868 // Keep new temperature in limits 00869 if(T > T_MAX){ 00870 //cerr << endl << " ... applying T_MAXlimit"; 00871 T = T_MAX; 00872 }else if(T < T_TRIPLE){ 00873 //cerr << endl << " ... applying T_TRIPLE limit"; 00874 T = T_TRIPLE; 00875 } 00876 00877 // Keep pressure in limits 00878 if(p < P_TRIPLE){ 00879 //cerr << endl << " ... applying P_TRIPLE limit"; 00880 p = P_TRIPLE; 00881 }else if(T > TB_LOW){ 00882 Pressure pbound = Boundaries::getpbound_T(T); 00883 if(p > pbound){ 00884 //cerr << endl << " ... applying pbound limit"; 00885 p = pbound; 00886 } 00887 }else{ 00888 Pressure psat = Boundaries::getSatPres_T(T); 00889 if(p > psat){ 00890 //cerr << endl << " ... applying psat limit"; 00891 p = psat; 00892 } 00893 } 00894 00895 ASSERT(guess.whichRegion()==2); 00896 ASSERT(Boundaries::isRegion2_pT(p,T)); 00897 00898 guess.setRegion2_pT(p,T); 00899 00900 niter++; 00901 00902 if(niter > maxiter){ 00903 throw new Exception("Solver2::solveRegion2: Exceeded maximum iterations"); 00904 } 00905 } 00906 00907 } 00908 }; 00909 00910 template<> 00911 SteamCalculator 00912 Solver2<Pressure,Temperature,0,0>::solveRegion3( 00913 const Pressure &p, const Temperature &T, const SteamCalculator &firstguess 00914 ); 00915 00916 #endif