00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include "WaterPropsIAPWS.h"
00017 #include "ctexceptions.h"
00018 #include "stringUtils.h"
00019 #include <cmath>
00020 #include <cstdio>
00021 #include <cstdlib>
00022
00023 namespace Cantera {
00024
00025
00026
00027
00028 const doublereal T_c = 647.096;
00029
00030 static const doublereal P_c = 22.064E6;
00031
00032 const doublereal Rho_c = 322.;
00033
00034 static const doublereal M_water = 18.015268;
00035
00036
00037
00038
00039
00040
00041
00042
00043 static const doublereal Rgas = 8.314371E3;
00044
00045 #ifndef MAX
00046 # define MAX(x,y) (( (x) > (y) ) ? (x) : (y))
00047 #endif
00048
00049 #ifndef MIN
00050 # define MIN(x,y) (( (x) < (y) ) ? (x) : (y))
00051 #endif
00052
00053
00054
00055 WaterPropsIAPWS:: WaterPropsIAPWS() :
00056 m_phi(0),
00057 tau(-1.0),
00058 delta(-1.0),
00059 iState(-30000)
00060 {
00061 m_phi = new WaterPropsIAPWSphi();
00062 }
00063
00064
00065
00066
00067
00068 WaterPropsIAPWS::WaterPropsIAPWS(const WaterPropsIAPWS &b) :
00069 m_phi(0),
00070 tau(b.tau),
00071 delta(b.delta),
00072 iState(b.iState)
00073 {
00074 m_phi = new WaterPropsIAPWSphi();
00075 m_phi->tdpolycalc(tau, delta);
00076 }
00077
00078
00079
00080
00081
00082 WaterPropsIAPWS & WaterPropsIAPWS::operator=(const WaterPropsIAPWS &b) {
00083 if (this == &b) return *this;
00084 tau = b.tau;
00085 delta = b.delta;
00086 iState = b.iState;
00087 m_phi->tdpolycalc(tau, delta);
00088 return *this;
00089 }
00090
00091
00092 WaterPropsIAPWS::~WaterPropsIAPWS() {
00093 delete (m_phi);
00094 m_phi = 0;
00095 }
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105 void WaterPropsIAPWS::calcDim(doublereal temperature, doublereal rho) {
00106 tau = T_c / temperature;
00107 delta = rho / Rho_c;
00108
00109
00110
00111 if (temperature > T_c) {
00112 iState = WATER_SUPERCRIT;
00113 } else {
00114 if (delta < 1.0) {
00115 iState = WATER_GAS;
00116 } else {
00117 iState = WATER_LIQUID;
00118 }
00119 }
00120 }
00121
00122
00123
00124 doublereal WaterPropsIAPWS::helmholtzFE() const {
00125 doublereal retn = m_phi->phi(tau, delta);
00126 doublereal temperature = T_c/tau;
00127 doublereal RT = Rgas * temperature;
00128 return (retn * RT);
00129 }
00130
00131
00132
00133
00134
00135
00136
00137 doublereal WaterPropsIAPWS::pressure() const {
00138 doublereal retn = m_phi->pressureM_rhoRT(tau, delta);
00139 doublereal rho = delta * Rho_c;
00140 doublereal temperature = T_c / tau;
00141 return (retn * rho * Rgas * temperature/M_water);
00142 }
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 doublereal WaterPropsIAPWS::density(doublereal temperature, doublereal pressure,
00160 int phase, doublereal rhoguess) {
00161
00162 doublereal deltaGuess = 0.0;
00163 if (rhoguess == -1.0) {
00164 if (phase != -1) {
00165 if (temperature > T_c) {
00166 rhoguess = pressure * M_water / (Rgas * temperature);
00167 } else {
00168 if (phase == WATER_GAS || phase == WATER_SUPERCRIT) {
00169 rhoguess = pressure * M_water / (Rgas * temperature);
00170 } else if (phase == WATER_LIQUID) {
00171
00172
00173
00174
00175 rhoguess = 1000.;
00176 } else if (phase == WATER_UNSTABLELIQUID || phase == WATER_UNSTABLEGAS) {
00177 throw Cantera::CanteraError("WaterPropsIAPWS::density",
00178 "Unstable Branch finder is untested");
00179 } else {
00180 throw Cantera::CanteraError("WaterPropsIAPWS::density",
00181 "unknown state: " + Cantera::int2str(phase));
00182 }
00183 }
00184 } else {
00185
00186
00187
00188
00189 rhoguess = pressure * M_water / (Rgas * temperature);
00190 }
00191
00192 }
00193 doublereal p_red = pressure * M_water / (Rgas * temperature * Rho_c);
00194 deltaGuess = rhoguess / Rho_c;
00195 setState_TR(temperature, rhoguess);
00196 doublereal delta_retn = m_phi->dfind(p_red, tau, deltaGuess);
00197 doublereal density_retn;
00198 if (delta_retn >0.0) {
00199 delta = delta_retn;
00200
00201
00202
00203
00204 density_retn = delta_retn * Rho_c;
00205
00206
00207
00208
00209 setState_TR(temperature, density_retn);
00210
00211
00212 } else {
00213 density_retn = -1.0;
00214 }
00215 return density_retn;
00216 }
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244 doublereal WaterPropsIAPWS::density_const(doublereal pressure,
00245 int phase, doublereal rhoguess) const {
00246 doublereal temperature = T_c / tau;
00247 doublereal deltaGuess = 0.0;
00248 doublereal deltaSave = delta;
00249 if (rhoguess == -1.0) {
00250 if (phase != -1) {
00251 if (temperature > T_c) {
00252 rhoguess = pressure * M_water / (Rgas * temperature);
00253 } else {
00254 if (phase == WATER_GAS || phase == WATER_SUPERCRIT) {
00255 rhoguess = pressure * M_water / (Rgas * temperature);
00256 } else if (phase == WATER_LIQUID) {
00257
00258
00259
00260
00261 rhoguess = 1000.;
00262 } else if (phase == WATER_UNSTABLELIQUID || phase == WATER_UNSTABLEGAS) {
00263 throw Cantera::CanteraError("WaterPropsIAPWS::density",
00264 "Unstable Branch finder is untested");
00265 } else {
00266 throw Cantera::CanteraError("WaterPropsIAPWS::density",
00267 "unknown state: " + Cantera::int2str(phase));
00268 }
00269 }
00270 } else {
00271
00272
00273
00274
00275 rhoguess = pressure * M_water / (Rgas * temperature);
00276 }
00277
00278 }
00279 doublereal p_red = pressure * M_water / (Rgas * temperature * Rho_c);
00280 deltaGuess = rhoguess / Rho_c;
00281
00282 delta = deltaGuess;
00283 m_phi->tdpolycalc(tau, delta);
00284
00285
00286 doublereal delta_retn = m_phi->dfind(p_red, tau, deltaGuess);
00287 doublereal density_retn;
00288 if (delta_retn > 0.0) {
00289 delta = delta_retn;
00290
00291
00292
00293
00294 density_retn = delta_retn * Rho_c;
00295
00296 } else {
00297 density_retn = -1.0;
00298 }
00299
00300 delta = deltaSave;
00301 m_phi->tdpolycalc(tau, delta);
00302 return density_retn;
00303 }
00304
00305
00306
00307
00308
00309
00310
00311 doublereal WaterPropsIAPWS::density() const {
00312 return (delta * Rho_c);
00313 }
00314
00315
00316
00317
00318
00319 doublereal WaterPropsIAPWS::temperature() const {
00320 return (T_c / tau);
00321 }
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334 doublereal WaterPropsIAPWS::psat_est(doublereal temperature) const {
00335
00336 static const doublereal A[8] = {
00337 -7.8889166E0,
00338 2.5514255E0,
00339 -6.716169E0,
00340 33.2239495E0,
00341 -105.38479E0,
00342 174.35319E0,
00343 -148.39348E0,
00344 48.631602E0
00345 };
00346 doublereal ps;
00347 if (temperature < 314.) {
00348 doublereal pl = 6.3573118E0 - 8858.843E0 / temperature
00349 + 607.56335E0 * pow(temperature, -0.6);
00350 ps = 0.1 * exp(pl);
00351 } else {
00352 doublereal v = temperature / 647.25;
00353 doublereal w = fabs(1.0-v);
00354 doublereal b = 0.0;
00355 for (int i = 0; i < 8; i++) {
00356 doublereal z = i + 1;
00357 b += A[i] * pow(w, ((z+1.0)/2.0));
00358 }
00359 doublereal q = b / v;
00360 ps = 22.093*exp(q);
00361 }
00362
00363
00364
00365 ps *= 1.0E6;
00366 return ps;
00367 }
00368
00369
00370
00371
00372
00373
00374 doublereal WaterPropsIAPWS::isothermalCompressibility() const {
00375 doublereal dpdrho_val = dpdrho();
00376 doublereal dens = delta * Rho_c;
00377 return (1.0 / (dens * dpdrho_val));
00378 }
00379
00380
00381
00382
00383
00384
00385
00386
00387 doublereal WaterPropsIAPWS::dpdrho() const {
00388 doublereal retn = m_phi->dimdpdrho(tau, delta);
00389 doublereal temperature = T_c/tau;
00390 doublereal val = retn * Rgas * temperature / M_water;
00391 return val;
00392 }
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403 doublereal WaterPropsIAPWS:: coeffPresExp() const {
00404 doublereal retn = m_phi->dimdpdT(tau, delta);
00405 return (retn);
00406 }
00407
00408
00409
00410
00411
00412
00413
00414 doublereal WaterPropsIAPWS:: coeffThermExp() const {
00415 doublereal kappa = isothermalCompressibility();
00416 doublereal beta = coeffPresExp();
00417 doublereal dens = delta * Rho_c;
00418 return (kappa * dens * Rgas * beta / M_water);
00419 }
00420
00421
00422
00423 doublereal WaterPropsIAPWS::Gibbs() const {
00424 doublereal gRT = m_phi->gibbs_RT();
00425 doublereal temperature = T_c/tau;
00426 return (gRT * Rgas * temperature);
00427 }
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443 void WaterPropsIAPWS::
00444 corr(doublereal temperature, doublereal pressure, doublereal &densLiq,
00445 doublereal &densGas, doublereal &delGRT) {
00446
00447 densLiq = density(temperature, pressure, WATER_LIQUID, densLiq);
00448 if (densLiq <= 0.0) {
00449 throw Cantera::CanteraError("WaterPropsIAPWS::corr",
00450 "Error occurred trying to find liquid density at (T,P) = "
00451 + Cantera::fp2str(temperature) + " " + Cantera::fp2str(pressure));
00452 }
00453 setState_TR(temperature, densLiq);
00454 doublereal gibbsLiqRT = m_phi->gibbs_RT();
00455
00456 densGas = density(temperature, pressure, WATER_GAS, densGas);
00457 if (densGas <= 0.0) {
00458 throw Cantera::CanteraError("WaterPropsIAPWS::corr",
00459 "Error occurred trying to find gas density at (T,P) = "
00460 + Cantera::fp2str(temperature) + " " + Cantera::fp2str(pressure));
00461 }
00462 setState_TR(temperature, densGas);
00463 doublereal gibbsGasRT = m_phi->gibbs_RT();
00464
00465 delGRT = gibbsLiqRT - gibbsGasRT;
00466 }
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478 void WaterPropsIAPWS::
00479 corr1(doublereal temperature, doublereal pressure, doublereal &densLiq,
00480 doublereal &densGas, doublereal &pcorr) {
00481
00482 densLiq = density(temperature, pressure, WATER_LIQUID, densLiq);
00483 if (densLiq <= 0.0) {
00484 throw Cantera::CanteraError("WaterPropsIAPWS::corr1",
00485 "Error occurred trying to find liquid density at (T,P) = "
00486 + Cantera::fp2str(temperature) + " " + Cantera::fp2str(pressure));
00487 }
00488 setState_TR(temperature, densLiq);
00489 doublereal prL = m_phi->phiR();
00490
00491 densGas = density(temperature, pressure, WATER_GAS, densGas);
00492 if (densGas <= 0.0) {
00493 throw Cantera::CanteraError("WaterPropsIAPWS::corr1",
00494 "Error occurred trying to find gas density at (T,P) = "
00495 + Cantera::fp2str(temperature) + " " + Cantera::fp2str(pressure));
00496 }
00497 setState_TR(temperature, densGas);
00498 doublereal prG = m_phi->phiR();
00499
00500 doublereal rhs = (prL - prG) + log(densLiq/densGas);
00501 rhs /= (1.0/densGas - 1.0/densLiq);
00502
00503 pcorr = rhs * Rgas * temperature / M_water;
00504 }
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525 doublereal WaterPropsIAPWS::psat(doublereal temperature, int waterState) {
00526 static int method = 1;
00527 doublereal densLiq = -1.0, densGas = -1.0, delGRT = 0.0;
00528 doublereal dp, pcorr;
00529 if (temperature >= T_c) {
00530 densGas = density(temperature, P_c, WATER_SUPERCRIT);
00531 setState_TR(temperature, densGas);
00532 return P_c;
00533 }
00534 doublereal p = psat_est(temperature);
00535 bool conv = false;
00536 for (int i = 0; i < 30; i++) {
00537 if (method == 1) {
00538 corr(temperature, p, densLiq, densGas, delGRT);
00539 doublereal delV = M_water * (1.0/densLiq - 1.0/densGas);
00540 dp = - delGRT * Rgas * temperature / delV;
00541 } else {
00542 corr1(temperature, p, densLiq, densGas, pcorr);
00543 dp = pcorr - p;
00544 }
00545 p += dp;
00546
00547 if ((method == 1) && delGRT < 1.0E-8) {
00548 conv = true;
00549 break;
00550 } else {
00551 if (fabs(dp/p) < 1.0E-9) {
00552 conv = true;
00553 break;
00554 }
00555 }
00556 }
00557
00558 if (waterState == WATER_LIQUID) {
00559 setState_TR(temperature, densLiq);
00560 } else if (waterState == WATER_GAS) {
00561 setState_TR(temperature, densGas);
00562 } else {
00563 throw Cantera::CanteraError("WaterPropsIAPWS::psat",
00564 "unknown water state input: " + Cantera::int2str(waterState));
00565 }
00566 return p;
00567 }
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579 int WaterPropsIAPWS::phaseState(bool checkState) const {
00580 if (checkState) {
00581 if (tau <= 1.0) {
00582 iState = WATER_SUPERCRIT;
00583 } else {
00584 doublereal T = T_c / tau;
00585 doublereal rho = delta * Rho_c;
00586
00587 doublereal rhoMidAtm = 0.5 * (1.01E5 * M_water / (8314.472 * 373.15) + 1.0E3);
00588 doublereal rhoMid = Rho_c + (T - T_c) * (Rho_c - rhoMidAtm) / (T_c - 373.15);
00589 int iStateGuess = WATER_LIQUID;
00590 if (rho < rhoMid) {
00591 iStateGuess = WATER_GAS;
00592 }
00593 doublereal kappa = isothermalCompressibility();
00594 if (kappa >= 0.0) {
00595 iState = iStateGuess;
00596 } else {
00597
00598 doublereal rhoDel = rho * 1.000001;
00599
00600
00601 doublereal deltaSave = delta;
00602 doublereal deltaDel = rhoDel / Rho_c;
00603 delta = deltaDel;
00604 m_phi->tdpolycalc(tau, deltaDel);
00605
00606 doublereal kappaDel = isothermalCompressibility();
00607 doublereal d2rhodp2 = (rhoDel * kappaDel - rho * kappa) / (rhoDel - rho);
00608 if (d2rhodp2 > 0.0) {
00609 iState = WATER_UNSTABLELIQUID;
00610 } else {
00611 iState = WATER_UNSTABLEGAS;
00612 }
00613
00614 delta = deltaSave;
00615
00616 m_phi->tdpolycalc(tau, delta);
00617 }
00618 }
00619 }
00620 return iState;
00621 }
00622
00623
00624
00625
00626
00627
00628 doublereal WaterPropsIAPWS::densSpinodalWater() const {
00629 doublereal temperature = T_c/tau;
00630 doublereal delta_save = delta;
00631
00632
00633
00634 if (temperature >= T_c - 0.001) {
00635 return Rho_c;
00636 }
00637 doublereal p = psat_est(temperature);
00638 doublereal rho_low = 0.0;
00639 doublereal rho_high = 1000;
00640
00641 doublereal densSatLiq = density_const(p, WATER_LIQUID);
00642 doublereal dens_old = densSatLiq;
00643 delta = dens_old / Rho_c;
00644 m_phi->tdpolycalc(tau, delta);
00645 doublereal dpdrho_old = dpdrho();
00646 if (dpdrho_old > 0.0) {
00647 rho_high = MIN(dens_old, rho_high);
00648 } else {
00649 rho_low = MAX(rho_low, dens_old);
00650 }
00651 doublereal dens_new = densSatLiq* (1.0001);
00652 delta = dens_new / Rho_c;
00653 m_phi->tdpolycalc(tau, delta);
00654 doublereal dpdrho_new = dpdrho();
00655 if (dpdrho_new > 0.0) {
00656 rho_high = MIN(dens_new, rho_high);
00657 } else {
00658 rho_low = MAX(rho_low, dens_new);
00659 }
00660 bool conv = false;
00661
00662 for (int it = 0; it < 50; it++) {
00663 doublereal slope = (dpdrho_new - dpdrho_old)/(dens_new - dens_old);
00664 if (slope >= 0.0) {
00665 slope = MAX(slope, dpdrho_new *5.0/ dens_new);
00666 } else {
00667 slope = -dpdrho_new;
00668
00669
00670 }
00671 doublereal delta_rho = - dpdrho_new / slope;
00672 if (delta_rho > 0.0) {
00673 delta_rho = MIN(delta_rho, dens_new * 0.1);
00674 } else {
00675 delta_rho = MAX(delta_rho, - dens_new * 0.1);
00676 }
00677 doublereal dens_est = dens_new + delta_rho;
00678 if (dens_est < rho_low) {
00679 dens_est = 0.5 * (rho_low + dens_new);
00680 }
00681 if (dens_est > rho_high) {
00682 dens_est = 0.5 * (rho_high + dens_new);
00683 }
00684
00685
00686 dens_old = dens_new;
00687 dpdrho_old = dpdrho_new;
00688 dens_new = dens_est;
00689
00690 delta = dens_new / Rho_c;
00691 m_phi->tdpolycalc(tau, delta);
00692 dpdrho_new = dpdrho();
00693 if (dpdrho_new > 0.0) {
00694 rho_high = MIN(dens_new, rho_high);
00695 } else if (dpdrho_new < 0.0) {
00696 rho_low = MAX(rho_low, dens_new);
00697 } else {
00698 conv = true;
00699 break;
00700 }
00701
00702 if (fabs(dpdrho_new) < 1.0E-5) {
00703 conv = true;
00704 break;
00705 }
00706 }
00707
00708 if (!conv) {
00709 throw Cantera::CanteraError(" WaterPropsIAPWS::densSpinodalWater()",
00710 " convergence failure");
00711 }
00712
00713 delta = delta_save;
00714 m_phi->tdpolycalc(tau, delta);
00715
00716 return dens_new;
00717 }
00718
00719
00720
00721
00722
00723
00724 doublereal WaterPropsIAPWS::densSpinodalSteam() const {
00725 doublereal temperature = T_c/tau;
00726 doublereal delta_save = delta;
00727
00728
00729
00730 if (temperature >= T_c - 0.001) {
00731 return Rho_c;
00732 }
00733 doublereal p = psat_est(temperature);
00734 doublereal rho_low = 0.0;
00735 doublereal rho_high = 1000;
00736
00737 doublereal densSatGas = density_const(p, WATER_GAS);
00738 doublereal dens_old = densSatGas;
00739 delta = dens_old / Rho_c;
00740 m_phi->tdpolycalc(tau, delta);
00741 doublereal dpdrho_old = dpdrho();
00742 if (dpdrho_old < 0.0) {
00743 rho_high = MIN(dens_old, rho_high);
00744 } else {
00745 rho_low = MAX(rho_low, dens_old);
00746 }
00747 doublereal dens_new = densSatGas * (0.99);
00748 delta = dens_new / Rho_c;
00749 m_phi->tdpolycalc(tau, delta);
00750 doublereal dpdrho_new = dpdrho();
00751 if (dpdrho_new < 0.0) {
00752 rho_high = MIN(dens_new, rho_high);
00753 } else {
00754 rho_low = MAX(rho_low, dens_new);
00755 }
00756 bool conv = false;
00757
00758 for (int it = 0; it < 50; it++) {
00759 doublereal slope = (dpdrho_new - dpdrho_old)/(dens_new - dens_old);
00760 if (slope >= 0.0) {
00761 slope = dpdrho_new;
00762
00763
00764 } else {
00765
00766 slope = MIN(slope, dpdrho_new *5.0 / dens_new);
00767
00768 }
00769 doublereal delta_rho = - dpdrho_new / slope;
00770 if (delta_rho > 0.0) {
00771 delta_rho = MIN(delta_rho, dens_new * 0.1);
00772 } else {
00773 delta_rho = MAX(delta_rho, - dens_new * 0.1);
00774 }
00775 doublereal dens_est = dens_new + delta_rho;
00776 if (dens_est < rho_low) {
00777 dens_est = 0.5 * (rho_low + dens_new);
00778 }
00779 if (dens_est > rho_high) {
00780 dens_est = 0.5 * (rho_high + dens_new);
00781 }
00782
00783
00784 dens_old = dens_new;
00785 dpdrho_old = dpdrho_new;
00786 dens_new = dens_est;
00787
00788 delta = dens_new / Rho_c;
00789 m_phi->tdpolycalc(tau, delta);
00790 dpdrho_new = dpdrho();
00791 if (dpdrho_new < 0.0) {
00792 rho_high = MIN(dens_new, rho_high);
00793 } else if (dpdrho_new > 0.0) {
00794 rho_low = MAX(rho_low, dens_new);
00795 } else {
00796 conv = true;
00797 break;
00798 }
00799
00800 if (fabs(dpdrho_new) < 1.0E-5) {
00801 conv = true;
00802 break;
00803 }
00804 }
00805
00806 if (!conv) {
00807 throw Cantera::CanteraError(" WaterPropsIAPWS::densSpinodalSteam()",
00808 " convergence failure");
00809 }
00810
00811 delta = delta_save;
00812 m_phi->tdpolycalc(tau, delta);
00813
00814 return dens_new;
00815 }
00816
00817
00818
00819
00820
00821 void WaterPropsIAPWS::setState_TR(doublereal temperature, doublereal rho) {
00822 calcDim(temperature, rho);
00823 m_phi->tdpolycalc(tau, delta);
00824 }
00825
00826
00827
00828
00829
00830 doublereal WaterPropsIAPWS::enthalpy() const {
00831 doublereal temperature = T_c/tau;
00832 doublereal hRT = m_phi->enthalpy_RT();
00833 return (hRT * Rgas * temperature);
00834 }
00835
00836
00837
00838
00839
00840 doublereal WaterPropsIAPWS::intEnergy() const {
00841 doublereal temperature = T_c / tau;
00842 doublereal uRT = m_phi->intEnergy_RT();
00843 return (uRT * Rgas * temperature);
00844 }
00845
00846
00847
00848
00849
00850 doublereal WaterPropsIAPWS::entropy() const {
00851 doublereal sR = m_phi->entropy_R();
00852 return (sR * Rgas);
00853 }
00854
00855
00856
00857
00858
00859 doublereal WaterPropsIAPWS::cv() const {
00860 doublereal cvR = m_phi->cv_R();
00861 return (cvR * Rgas);
00862 }
00863
00864
00865
00866 doublereal WaterPropsIAPWS::cp() const {
00867 doublereal cpR = m_phi->cp_R();
00868 return (cpR * Rgas);
00869 }
00870
00871
00872
00873 doublereal WaterPropsIAPWS::molarVolume() const {
00874 doublereal rho = delta * Rho_c;
00875 return (M_water / rho);
00876 }
00877
00878 }