00001 #include "fluxforco2inj.h"
00002 #include "numericmethods.h"
00003 #include "exception.h"
00004 #include <algorithm>
00005 #include "sfunctions.h"
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00020 FluxForCO2Inj::FluxForCO2Inj(OrthoMesh &mesh,const VecDouble &cK,double v1,double v2,double sr1,double sr2,double pc,double pw,double g)
00021 :FaceFluxFunction(2),
00022 m_mesh(mesh),
00023 m_fMob(v1,v2,sr1,sr2),
00024 m_fMobGrav(v1,v2,sr1,sr2),
00025 m_DfMob(v1,v2,sr1,sr2),
00026 m_DfMobGrav(v1,v2,sr1,sr2),
00027 m_Vn(mesh.numFaces()),
00028 m_cK(cK)
00029 {
00030 m_maxK = m_cK.linfty_norm();
00031 m_pc=pc;
00032 m_pw=pw;
00033 m_g = g;
00034 m_G = (m_pc-m_pw)*m_g;
00035
00036 }
00037
00038 FluxForCO2Inj::~FluxForCO2Inj()
00039 {
00040 }
00041
00042
00053 void FluxForCO2Inj::fluxAtFace(VecDouble &vFlux, const FaceInfo &face,const VecDouble &Q1,const VecDouble &Q2)
00054 {
00055 assert( (unsigned) face.index < m_Vn.size());
00056 assert( (unsigned) face.cell1 < m_cK.size());
00057 assert( (unsigned) face.cell2 < m_cK.size());
00058 assert(Q1.size() == 2);
00059 assert(Q2.size() == 2);
00060
00061 double Sc1 = Q1(_Sc);
00062 double Sc2 = Q2(_Sc);
00063
00064 double vel = m_Vn(face.index);
00065 double C1 = getVolFractionOfCO2InWater(Sc1,Q1(_Sd));
00066 double C2 = getVolFractionOfCO2InWater(Sc2,Q2(_Sd));
00067
00068 vFlux(_Sc) = (m_fMob(Sc1)+m_fMob(Sc2))/2.0*vel;
00069 vFlux(_Sd) = ( (1.0-m_fMob(Sc1))*C1 + (1.0-m_fMob(Sc2))*C2 )/2.0*vel;
00070
00071
00072
00073 if (face.normal == Z && face.at_boundary )
00074 {
00075 const VecDouble &cK = getPermeability();
00076 double K1 = cK(face.cell1);
00077 double K2 = cK(face.cell2);
00078 vFlux(_Sd) += ( C1*m_fMobGrav(Sc1)*(1-C1) + C2*m_fMobGrav(Sc2)*(1-C2) ) * K1*K2/(K1+K2) * m_G;
00079
00080 vFlux(_Sc) -= ( m_fMobGrav(Sc1)*(1-C1) + m_fMobGrav(Sc2)*(1-C2) ) * K1*K2/(K1+K2) * m_G;
00081 }
00082
00083
00084
00085 }
00086
00087
00088
00089 double FluxForCO2Inj::maxLocalCharVelocity(const FaceInfo &face,const VecDouble &Q1, const VecDouble &Q2)
00090 {
00091 double grav = std::max(m_cK(face.cell1),m_cK(face.cell2))*m_G;
00092 double Sc1 = Q1(_Sc);
00093 double Sc2 = Q2(_Sc);
00094 if (Sc2 < Sc1)
00095 std::swap(Sc1,Sc2);
00096 double mobMaxNormGrad=m_DfMob.getMaxNorm(Sc1,Sc2);
00097 if (mobMaxNormGrad==0)
00098 mobMaxNormGrad=1;
00099 if (face.normal == OrthoMesh::NORMAL_Z)
00100 return fabs(m_Vn(face.index)*mobMaxNormGrad) + fabs(m_DfMobGrav.getMaxNorm(Sc1,Sc2)*grav);
00101 else
00102 return fabs(m_Vn(face.index)*mobMaxNormGrad);
00103 }
00104
00105
00106
00107
00108
00109 double FluxForCO2Inj::getAquousDensity(double pc,double pw,VecDouble& vQ)
00110 {
00111 double Sa = (1.0-vQ(_Sc));
00112 if (fabs(Sa) < 1.0E-10)
00113 return pw;
00114
00115 double C = vQ(_Sd)/Sa;
00116 return C*pc + (1-C)*pw;
00117
00118 }
00119
00120
00121 double FluxForCO2Inj::getVolFractionOfCO2InWater(double Sc,double Cd)
00122 {
00123 if (Sc < 1.0)
00124 return Cd/(1.0-Sc);
00125 else
00126 return 0.;
00127 }
00128
00129
00130 void FluxForCO2Inj::updateDynamicData(DynamicBase &dynMod)
00131 {
00132 dynMod.getNormalVelocityAtFaces(m_Vn);
00133 }
00134
00135 void FluxForCO2Inj::maxGlobalCharVelocity(double vel[3])
00136 {
00137 VecDouble maxAbsVel(3);
00138
00139 m_mesh.getMaxAbsVelocitiesInFacesByComponents(maxAbsVel,m_Vn);
00140
00141 vel[X]= maxAbsVel(X)*m_DfMob.getMaxNorm(0,1);
00142 vel[Y]= maxAbsVel(Y)*m_DfMob.getMaxNorm(0,1);
00143 vel[Z]= maxAbsVel(Z)*m_DfMob.getMaxNorm(0,1) + m_maxK*fabs(m_DfMobGrav.getMaxNorm(0,1)*m_G);
00144
00145 }
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
