00001 #include "flashco2brine.h"
00002 #include "flash.h"
00003 #include <fstream>
00004 #include <algorithm>
00005 #include <assert.h>
00006 #include "units.h"
00007 #include "hdf5orthowriter.h"
00008 #include "numericmethods.h"
00009
00010
00011 #define R Flash::R
00012
00013
00014
00015
00016 FlashCO2Brine::FlashCO2Brine(OrthoMesh &mesh, double referenceT, double v1, double v2)
00017 :FlashCompositional(2,3),m_T(referenceT),m_mesh(mesh),m_v1(v1),m_v2(v2)
00018 {
00019 m_volCell = mesh.cellVolume();
00020
00021
00022
00023 m_A=0.38384020e-3*m_T*m_T -0.55953850*m_T + 0.30429268e+3 -0.72044305e+5/m_T + 0.63003388e+7/(m_T*m_T);
00024 m_B=-0.57709332e-5*m_T*m_T +0.82764653e-2*m_T -0.43813556e+1 +0.10144907e+4/m_T - 0.86777045e+5/(m_T*m_T);
00025
00026
00027 m_components_molar_mass.reinit(3);
00028 m_components_molar_mass(WATER)=18.01528e-3;
00029 m_components_molar_mass(CO2 )=44.01e-3;
00030 m_components_molar_mass(SALT )=58.443e-3;
00031
00032 }
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 void FlashCO2Brine::flash(double P, const VecDouble &compTotalMoles2, FlashData &data)
00043 {
00044
00045 static VecDouble compTotalMoles(3);
00046 compTotalMoles(0)=std::max(compTotalMoles2(0),0.0);
00047 compTotalMoles(1)=std::max(compTotalMoles2(1),0.0);
00048 compTotalMoles(2)=std::max(compTotalMoles2(2),0.0);
00049
00050
00051 assert(compTotalMoles.size()==3);
00052 assert(data.getNPhases()==2);
00053 assert(data.getNComponents()==3);
00054
00055
00056 assert(compTotalMoles(0)>=0.0);
00057 assert(compTotalMoles(1)>=0.0);
00058 assert(compTotalMoles(2)>=0.0);
00059
00060
00061
00062 if (compTotalMoles(WATER) == 0.0)
00063 {
00064 data.zeroEntries();
00065 data.setMoles(CO2_RICH,CO2,compTotalMoles(CO2));
00066 return;
00067 }
00068
00069
00070 Flash::BrineCO2Data brineco2Data;
00071 double nT=compTotalMoles(WATER)+compTotalMoles(CO2)+compTotalMoles(SALT);
00072
00073 double zH2O =compTotalMoles(WATER)/nT;
00074 double zCO2 =compTotalMoles(CO2) /nT;
00075 double zSalt=compTotalMoles(SALT) /nT;
00076
00077
00078 Flash::ProgressiveFlashBrineCO2(brineco2Data, m_T, Units::NewtonToBars(P), zSalt, zH2O, zCO2);
00079
00080 data.setMoles(CO2_RICH,WATER,nT*brineco2Data.betaC*brineco2Data.yH2O);
00081 data.setMoles(CO2_RICH,CO2 ,nT*brineco2Data.betaC*brineco2Data.yCO2);
00082 data.setMoles(CO2_RICH,SALT ,0);
00083 data.setMoles(AQUEOUS ,WATER,nT*brineco2Data.betaA*brineco2Data.xH2O);
00084 data.setMoles(AQUEOUS ,CO2 ,nT*brineco2Data.betaA*brineco2Data.xCO2);
00085 data.setMoles(AQUEOUS ,SALT ,nT*brineco2Data.betaA*brineco2Data.xNaCl);
00086
00087 if ( brineco2Data.betaS != 0)
00088 printf("Warning, Salt Production\n");
00089 if ( brineco2Data.stateCO2 == Flash::LIQUID)
00090 printf("Warning, Liquid CO2 produced P=%g Bars\n",Units::NewtonToBars(P));
00091 }
00092
00093
00094
00095
00101 double FlashCO2Brine::getFluidCompressibility(double P, FlashData &data)
00102 {
00103
00104 Flash::BrineCO2Data brineco2Data;
00105 double h=P*1.e-07;
00106 double nT;
00107 static VecDouble compC(3);
00108 data.getMixtureComposition(compC,nT);
00109
00110 double Ph=Units::NewtonToBars(P+h);
00111 if (compC(WATER) == 0.0)
00112 {
00113 brineco2Data.xCO2 = 0.0;
00114 brineco2Data.xH2O = 0.0;
00115 brineco2Data.xNaCl= 0.0;
00116 brineco2Data.yCO2 = 1.0;
00117 brineco2Data.yH2O = 0;
00118 brineco2Data.betaC = compC(CO2)/(compC(CO2)+compC(SALT));
00119 brineco2Data.betaS = 1.-brineco2Data.betaC;
00120 brineco2Data.stateCO2=Flash::GAS;
00121 }
00122 else
00123 Flash::ProgressiveFlashBrineCO2(brineco2Data, m_T, Ph, compC(SALT), compC(WATER), compC(CO2));
00124 double Vpos=Units::cm3Tom3(getTotalMolarVolume(brineco2Data,Ph,m_T));
00125
00126
00127 Ph=Units::NewtonToBars(P-h);
00128 if (compC(WATER) != 0.0)
00129 {
00130 Flash::ProgressiveFlashBrineCO2(brineco2Data, m_T, Ph, compC(SALT), compC(WATER), compC(CO2));
00131 }
00132 double Vneg=Units::cm3Tom3(getTotalMolarVolume(brineco2Data,Ph,m_T));
00133
00134
00135 return -nT*(Vpos-Vneg)/(2*h);
00136 }
00137
00138
00139
00140 void FlashCO2Brine::getPhasesVolume(double P, const FlashData &data,VecDouble &phasesVol)
00141 {
00142 assert(phasesVol.size()==2);
00143 assert(data.getNPhases()==2);
00144 assert(data.getNComponents()==3);
00145 double PBars=Units::NewtonToBars(P);
00146 phasesVol(AQUEOUS)=Units::cm3Tom3(data.getMoles(AQUEOUS,WATER)*LIQ_WATER_MOLAR_VOLUME + data.getMoles(AQUEOUS,SALT)*DISSOLVED_SALT_MOLAR_VOLUME + data.getMoles(AQUEOUS,CO2)*getDissolvedCO2ApparentVolume(PBars));
00147 phasesVol(CO2_RICH)=Units::cm3Tom3((data.getMoles(CO2_RICH,WATER)+data.getMoles(CO2_RICH,CO2))*Flash::ComputeVolumeCO2(m_T,PBars));
00148
00149 }
00150
00151
00152
00153
00154 const VecDouble& FlashCO2Brine::getComponentsMolarMass() const
00155 {
00156 return m_components_molar_mass;
00157 }
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00183 void FlashCO2Brine::printFlashPropertiesTable(double P0,double P1,unsigned nP, double T,double zSalt)
00184 {
00185 std::string solid;
00186 double dP=(P1-P0)/nP;
00187 double dZ=(1.0-zSalt)/nP;
00188 Flash::BrineCO2Data brineco2Data;
00189 static VecDouble m_aq(3),m_co(3),m_t(3),w(3);
00190 double P=P0;
00191 Matrix GradVaq(3,3),GradVco(3,3);
00192 printf("#Table of properties of the flash computations\n");
00193 printf("#Fixed Parameters: T=%g zSalt=%g\n",T,zSalt);
00194 printf("#PRESSURE\tzCO2\tMolar Volume(cm^3)\tFluid Compressibility\tdSg/dm_g*m_aq");
00195 for (unsigned i=0;i<=nP;P+=dP,i++)
00196 {
00197 double zCO2=0;
00198 for (unsigned j=0;j<nP;zCO2+=dZ,j++)
00199 {
00200 double zH2O=(1-zCO2-zSalt);
00201 Flash::ProgressiveFlashBrineCO2(brineco2Data, T, P, zSalt, zH2O, zCO2);
00202
00203 m_t(SALT )=zSalt;
00204 m_t(WATER)=zH2O;
00205 m_t(CO2 )=zCO2;
00206
00207 getGradient(GradVaq,GradVco,m_t,P,T,0.01);
00208 getMolesAqueousPhase(m_aq,brineco2Data,1);
00209 getMolesCO2RichPhase(m_co,brineco2Data,1);
00210 GradVaq.vmult(w,m_co);
00211
00212 if (brineco2Data.betaS != 0.0)
00213 solid="S";
00214 else
00215 solid="";
00216
00217
00218 printf("%g\t%g\t%g\t<%g,%g,%g> %s\n",P,zCO2,getTotalMolarVolume(brineco2Data,P,T),w(0),w(1),w(2),solid.c_str());
00219 }
00220 printf("\n");
00221 }
00222 }
00223
00224
00225 void FlashCO2Brine::printCO2RichMolarVolume(double P0,double P1,unsigned nP, double T)
00226 {
00227 double dP=(P1-P0)/nP;
00228 double P=P0;
00229 std::ofstream file1("co2GasMolarVolume.dat");
00230 std::ofstream file2("co2CriticMolarVolume.dat");
00231 char str[1500];
00232 sprintf(str,"#Molar Volume of Gas CO2\n#Parameters: T=%gK\n#PRESSURE (BARS)\tVolume (cm^3)\n",T);
00233 file1 << str;
00234 sprintf(str,"#Molar Volume of Super Critic CO2\n#Parameters: T=%gK\n#PRESSURE (BARS)\tVolume (cm^3)\n",T);
00235 file2 << str;
00236
00237 for (unsigned i=0;i<nP;P+=dP,i++)
00238 {
00239 file1 << P << "\t" << Flash::ComputeVolumeCO2(T,P,Flash::GAS)<<std::endl;
00240 file2 << P << "\t" << Flash::ComputeVolumeCO2(T,P,Flash::SUPERCRITICAL)<<std::endl;
00241 }
00242 }
00243
00244
00245
00246
00247
00248
00256 double FlashCO2Brine::getTotalMolarVolume(Flash::BrineCO2Data &data,double P, double T)
00257 {
00258 double tV;
00259
00260 tV=getCO2MolarVolume(data,P,T);
00261
00262 tV+=getAqueousMolarVolume(data,P,T);
00263
00264
00265
00266
00267 return tV;
00268 }
00269
00270
00271
00272
00273 double FlashCO2Brine::getAqueousMolarVolume(Flash::BrineCO2Data &data,double P,double T)
00274 {
00275 return data.betaA*(data.xH2O*LIQ_WATER_MOLAR_VOLUME + data.xNaCl*DISSOLVED_SALT_MOLAR_VOLUME + data.xCO2*getDissolvedCO2ApparentVolume(P));
00276 }
00277
00278 double FlashCO2Brine::getCO2MolarVolume(Flash::BrineCO2Data &data,double P, double T)
00279 {
00280 return data.betaC*Flash::ComputeVolumeCO2(T,P,data.stateCO2);
00281 }
00282
00288 double FlashCO2Brine::getPureCO2MolarVolume(double P, double T)
00289 {
00290 return Units::cm3Tom3(Flash::ComputeVolumeCO2(T,Units::NewtonToBars(P)));
00291 }
00292
00293
00294
00295
00301 double FlashCO2Brine::getDissolvedCO2ApparentVolume(double P)
00302 {
00303
00304 return LIQ_WATER_MOLAR_VOLUME*(1+m_A+P*m_B);
00305 }
00306
00307
00308
00309
00310 void FlashCO2Brine::getVolGradient(VecDouble &GradVaq,VecDouble &GradVco,double P,double T,double zSalt,double zH2O,double zCO2,double h)
00311 {
00312 assert(GradVaq.size()==3);
00313 assert(GradVco.size()==3);
00314 Flash::BrineCO2Data brineco2Data;
00315 double Hpos,Hneg;
00316 double VposAq,VnegAq,VposCO,VnegCO;
00317 double dx;
00318
00319
00320
00321
00322 dx=2*h;
00323 Hpos=zH2O+h;
00324 Hneg=zH2O-h;
00325 if (Hpos>1.0)
00326 {
00327 Hpos=zH2O;
00328 dx=h;
00329 }
00330 if (Hneg<0.0)
00331 {
00332 Hneg=zH2O;
00333 dx=h;
00334 }
00335 Flash::ProgressiveFlashBrineCO2(brineco2Data, T, P, zSalt, Hpos, zCO2);
00336 VposAq=getAqueousMolarVolume(brineco2Data,P,T);
00337 VposCO=getCO2MolarVolume(brineco2Data,P,T);
00338 Flash::ProgressiveFlashBrineCO2(brineco2Data, T, P, zSalt, Hneg, zCO2);
00339 VnegAq=getAqueousMolarVolume(brineco2Data,P,T);
00340 VnegCO=getCO2MolarVolume(brineco2Data,P,T);
00341 GradVaq(FlashCO2Brine::WATER)=(VposAq-VnegAq)/(dx);
00342 GradVco(FlashCO2Brine::WATER)=(VposCO-VnegCO)/(dx);
00343
00344
00345
00346 dx=2*h;
00347 Hpos=zCO2+h;
00348 Hneg=zCO2-h;
00349 if (Hpos>1.0)
00350 {
00351 Hpos=zCO2;
00352 dx=h;
00353 }
00354 if (Hneg<0.0)
00355 {
00356 Hneg=zCO2;
00357 dx=h;
00358 }
00359 Flash::ProgressiveFlashBrineCO2(brineco2Data, T, P, zSalt, zH2O, Hpos);
00360 VposAq=getAqueousMolarVolume(brineco2Data,P,T);
00361 VposCO=getCO2MolarVolume(brineco2Data,P,T);
00362 Flash::ProgressiveFlashBrineCO2(brineco2Data, T, P, zSalt, zH2O, Hneg);
00363 VnegAq=getAqueousMolarVolume(brineco2Data,P,T);
00364 VnegCO=getCO2MolarVolume(brineco2Data,P,T);
00365 GradVaq(FlashCO2Brine::CO2)=(VposAq-VnegAq)/(dx);
00366 GradVco(FlashCO2Brine::CO2)=(VposCO-VnegCO)/(dx);
00367
00368
00369
00370 dx=2*h;
00371 Hpos=zSalt+h;
00372 Hneg=zSalt-h;
00373 if (Hpos>1.0)
00374 {
00375 Hpos=zSalt;
00376 dx=h;
00377 }
00378 if (Hneg<0.0)
00379 {
00380 Hneg=zSalt;
00381 dx=h;
00382 }
00383 Flash::ProgressiveFlashBrineCO2(brineco2Data, T, P, Hpos, zH2O, zCO2);
00384 VposAq=getAqueousMolarVolume(brineco2Data,P,T);
00385 VposCO=getCO2MolarVolume(brineco2Data,P,T);
00386 Flash::ProgressiveFlashBrineCO2(brineco2Data, T, P, Hneg, zH2O, zCO2);
00387 VnegAq=getAqueousMolarVolume(brineco2Data,P,T);
00388 VnegCO=getCO2MolarVolume(brineco2Data,P,T);
00389 GradVaq(FlashCO2Brine::SALT)=(VposAq-VnegAq)/(dx);
00390 GradVco(FlashCO2Brine::SALT)=(VposCO-VnegCO)/(dx);
00391
00392
00393
00394
00395 }
00396
00397
00403 void FlashCO2Brine::getMolesAqueousPhase(VecDouble &M_aq,const Flash::BrineCO2Data &data,double nT)
00404 {
00405 assert(M_aq.size()==3);
00406 M_aq(FlashCO2Brine::WATER)=nT*data.betaA*data.xH2O;
00407 M_aq(FlashCO2Brine::CO2)= nT*data.betaA*data.xCO2;
00408 M_aq(FlashCO2Brine::SALT)= nT*data.betaA*data.xNaCl;
00409
00410 }
00411
00412
00413
00419 void FlashCO2Brine::getMolesCO2RichPhase(VecDouble &M_g,const Flash::BrineCO2Data &data,double nT)
00420 {
00421 assert(M_g.size()==3);
00422 M_g(FlashCO2Brine::WATER)=nT*data.betaC*data.yH2O;
00423 M_g(FlashCO2Brine::CO2)= nT*data.betaC*data.yCO2;
00424 M_g(FlashCO2Brine::SALT)= 0.0;
00425
00426 }
00427
00436 void FlashCO2Brine::getGradient(Matrix &gradMaq,Matrix &gradMg,const VecDouble &m_T,double P,double T, double h)
00437 {
00438 assert(WATER==0 && CO2==1 && SALT==2);
00439 static VecDouble maqPos(3),maqNeg(3);
00440 static VecDouble mgPos(3),mgNeg(3);
00441 static VecDouble m_TPos(3),m_TNeg(3);
00442
00443 double dx=2*h;
00444
00445 for (unsigned i=0;i<3;i++)
00446 {
00447 m_TPos=m_T;
00448 m_TNeg=m_T;
00449 m_TPos(i)+=h;
00450 m_TNeg(i)-=h;
00451 flash(maqPos,mgPos,m_TPos,P,T);
00452 flash(maqNeg,mgNeg,m_TNeg,P,T);
00453
00454 gradMaq(WATER,i)=(maqPos(WATER) - maqNeg(WATER))/dx;
00455 gradMaq(CO2 ,i)=(maqPos(CO2 ) - maqNeg(CO2 ))/dx;
00456 gradMaq(SALT ,i)=(maqPos(SALT ) - maqNeg(SALT ))/dx;
00457 gradMg (WATER,i)=(mgPos (WATER) - mgNeg(WATER ))/dx;
00458 gradMg (CO2 ,i)=(mgPos (CO2 ) - mgNeg(CO2 ))/dx;
00459 gradMg (SALT ,i)=(mgPos (SALT ) - mgNeg(SALT ))/dx;
00460 }
00461
00462 }
00463
00464 void FlashCO2Brine::flash(VecDouble &m_aq,VecDouble &m_g,const VecDouble &m_T,double P,double T)
00465 {
00466 assert(m_T.size()==3);
00467 assert(m_aq.size()==3);
00468 assert(m_g.size()==3);
00469
00470
00471 Flash::BrineCO2Data brineco2Data;
00472 double nT=m_T(WATER)+m_T(CO2)+m_T(SALT);
00473
00474 double zH2O =m_T(WATER)/nT;
00475 double zCO2 =m_T(CO2) /nT;
00476 double zSalt=m_T(SALT) /nT;
00477
00478 Flash::ProgressiveFlashBrineCO2(brineco2Data, T, P, zSalt, zH2O, zCO2);
00479 getMolesCO2RichPhase(m_g ,brineco2Data,nT);
00480 getMolesAqueousPhase(m_aq,brineco2Data,nT);
00481
00482
00483 }
00484
00485
00486
00487
00488
00489
00490
00491
00492 void FlashCO2Brine::getPhasesViscosities(double P,const FlashData &data,VecDouble &visc)
00493 {
00494 visc(AQUEOUS)= m_v1;
00495 visc(CO2_RICH)=m_v2;
00496 }
00497
00498
00499 void FlashCO2Brine::printOutput()
00500 {
00501 HDF5OrthoWriter &hdf5 = HDF5OrthoWriter::getHDF5OrthoWriter();
00502 unsigned nCells = m_mesh.numCells();
00503 VecDouble vc(nCells);
00504 VecDouble vcSatC(nCells);
00505 FlashData data(numPhases(),numComponents(),NULL);
00506 VecDouble vv(m_mesh.numVertices());
00507 VecDouble phasesVol(2);
00508 for (unsigned i=0;i<nCells;i++)
00509 {
00510 FlashCompositional::flash(i,data);
00511 getPhasesVolume(getDynamicModule().getPressureAtCells()(i),data,phasesVol);
00512 vcSatC(i)=phasesVol(CO2_RICH);
00513 }
00514 m_mesh.projectCentralValuesAtVertices(vcSatC,vv);
00515 hdf5.writeScalarField(vv,"SatC");
00516
00517 getFlashDataArray().getVectorData(CO2_RICH,WATER,vc);
00518 m_mesh.projectCentralValuesAtVertices(vc,vv);
00519 hdf5.writeScalarField(vv,"ngw");
00520
00521
00522 getFlashDataArray().getVectorData(CO2_RICH,CO2,vc);
00523 m_mesh.projectCentralValuesAtVertices(vc,vv);
00524 hdf5.writeScalarField(vv,"ngc");
00525
00526
00527 getFlashDataArray().getVectorData(AQUEOUS,WATER,vc);
00528 m_mesh.projectCentralValuesAtVertices(vc,vv);
00529 hdf5.writeScalarField(vv,"naqw");
00530
00531 getFlashDataArray().getVectorData(AQUEOUS,CO2,vc);
00532 m_mesh.projectCentralValuesAtVertices(vc,vv);
00533 hdf5.writeScalarField(vv,"naqc");
00534
00535 {
00536 VecDouble vc2(nCells);
00537 getFlashDataArray().getVectorData(AQUEOUS,WATER,vc2);
00538 NumericMethods::divideEachElement(vc,vc,vc2);
00539 m_mesh.projectCentralValuesAtVertices(vc,vv);
00540 hdf5.writeScalarField(vv,"Cco2");
00541
00542
00543 }
00544
00545
00546
00547 getFlashDataArray().getVectorData(AQUEOUS,SALT,vc);
00548 m_mesh.projectCentralValuesAtVertices(vc,vv);
00549 hdf5.writeScalarField(vv,"naqs");
00550
00551
00552 m_mesh.projectCentralValuesAtVertices(vcSatC,vv);
00553 hdf5.writeScalarField(vv,"SatC");
00554
00555
00556 }
00557
00558
00559 void FlashCO2Brine::getTotalVolumeDerivatives(double P,FlashData &data,VecDouble &dv_dm)
00560 {
00561 static VecDouble m(3);
00562 data.getTotalComponentsMoles(m);
00563 getNumericTotalVolumeDerivative(P,m,dv_dm,1e-04);
00564
00565
00566
00567 }
