Skip to content

Commit

Permalink
Merge pull request #114 from wenqing/fixing
Browse files Browse the repository at this point in the history 
Several improvements
  • Loading branch information
@wenqing
wenqing authored Nov 2, 2018
2 parents 9f957a3 + 22b7655 commit 5c4b28c
Show file tree
Hide file tree
Showing 12 changed files with 224 additions and 208 deletions.
8 changes: 8 additions & 0 deletions FEM/FEMEnums.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -283,6 +283,14 @@ enum SolidReactiveSystem
SolidReactiveSystem convertSolidReactiveSystem(const std::string& reactive_string);
std::string convertSolidReactiveSystemToString(SolidReactiveSystem reactive_system);

enum InitDataReadWriteType
{
NO_IO = 0,
READ,
WRITE,
READ_WRITE,
};

} // end namespace FiniteElement

struct TimType
Expand Down
6 changes: 4 additions & 2 deletions FEM/fem_ele.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -611,7 +611,8 @@ void CElement::computeJacobian(const int gp, const int order, const bool need_in
DetJac = _Jacobian[0] * _Jacobian[3] - _Jacobian[1] * _Jacobian[2];
if (fabs(DetJac) < MKleinsteZahl)
{
std::cout << "\n*** Jacobian: Det == 0 " << DetJac << "\n";
std::cout << "\n*** Jacobian: Det == 0 " << DetJac
<< " in element " << MeshElement->GetIndex() << "\n";
abort();
}
invJacobian[0] = _Jacobian[3];
Expand Down Expand Up @@ -662,7 +663,8 @@ void CElement::computeJacobian(const int gp, const int order, const bool need_in

if (fabs(DetJac) < MKleinsteZahl)
{
std::cout << "\n*** Jacobian: DetJac == 0 " << DetJac << "\n";
std::cout << "\n*** Jacobian: DetJac == 0 " << DetJac
<< " in element " << MeshElement->GetIndex() << "\n";
abort();
}
invJacobian[0] = _Jacobian[4] * _Jacobian[8] - _Jacobian[7] * _Jacobian[5];
Expand Down
219 changes: 87 additions & 132 deletions FEM/fem_ele_std.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1522,100 +1522,48 @@ void CFiniteElementStd::CalNodalEnthalpy()
**************************************************************************/
double CFiniteElementStd::CalCoefMass()
{
int Index = MeshElement->GetIndex();
const int Index = MeshElement->GetIndex();
double val = 0.0;
double humi = 1.0;
double rhov = 0.0;
double arg[2];
double biot_val, poro_val = 0.0, rho_val = 0.0, Se;
int tr_phase = 0; // SB, BG
double saturation = 0.0; // SB, BG
CompProperties* m_cp = NULL;
double drho_dp_rho = 0.0; // drho/dp * 1/rho

switch (PcsType)
{
default:
std::cout << "Fatal error in CalCoefMass: No valid PCS type"
<< "\n";
exit(EXIT_FAILURE);
break;
case EPT_LIQUID_FLOW: // Liquid flow
// Is this really needed?
val = MediaProp->StorageFunction(Index, unit, pcs->m_num->ls_theta);

// get drho/dp/rho from material model or direct input
if (FluidProp->compressibility_model_pressure > 0)
{
rho_val = FluidProp->Density();
arg[0] = interpolate(NodalVal1); // p
arg[1] = interpolate(NodalValC1); // T
drho_dp_rho = FluidProp->drhodP(arg) / rho_val;
}
else
drho_dp_rho = FluidProp->drho_dp;

// JT 2010, needed storage term and fluid compressibility...
// We derive here the storage at constant strain, or the inverse of Biot's "M" coefficient
// Assumptions are the most general possible:: Invarience under "pi" (Detournay & Cheng) loading.
// Se = 1/M = poro/Kf + (alpha-poro)/Ks :: Cf = 1/Kf = 1/rho * drho/dp :: alpha = 1 - K/Ks
// Second term (of Se) below vanishes for incompressible grains
// WW if(D_Flag > 0 && rho_val > MKleinsteZahl)
if (dm_pcs && MediaProp->storage_model == 7) // Add MediaProp->storage_model. 29.09.2011. WW
{
biot_val = SolidProp->biot_const;
poro_val = MediaProp->Porosity(Index, pcs->m_num->ls_theta);
val = 0.; // WX:04.2013

// WW if(SolidProp->K == 0) //WX: if HM Partitioned, K still 0 here
if (fabs(SolidProp->K) < DBL_MIN) // WW 29.09.2011
val = MediaProp->StorageFunction(Index, unit, pcs->m_num->ls_theta);
double drho_dp_rho = 0.0;
// get drho/dp/rho from material model or direct input
if (FluidProp->compressibility_model_pressure > 0)
{
if (SolidProp->Youngs_mode < 10 || SolidProp->Youngs_mode > 13) // JM,WX: 2013
SolidProp->K = SolidProp->E / 3 / (1 - 2 * SolidProp->PoissonRatio);
else
{
double E_av; // average Youngs modulus
double nu_av; // average Poisson ratio
double nu_ai; // Poisson ratio perpendicular to the plane of isotropie, due to strain in the
// plane of isotropie
double nu_ia; // Poisson ratio in the plane of isotropie, due to strain perpendicular to the
// plane of isotropie
double nu_i; // Poisson ratio in the plane of isotropy

E_av = 2. / 3. * (*SolidProp->data_Youngs)(0) + 1. / 3. * (*SolidProp->data_Youngs)(1);

nu_ia = (*SolidProp->data_Youngs)(2);
nu_ai = nu_ia * (*SolidProp->data_Youngs)(1)
/ (*SolidProp->data_Youngs)(0); // nu_ai=nu_ia*Ea/Ei

nu_i = SolidProp->Poisson_Ratio();
// 12 13 21 23 31 32
// ai ai ia ii ia ii
nu_av = 1. / 3. * (nu_ai + nu_ia + nu_i);

SolidProp->K = E_av / 3 / (1 - 2 * nu_av);
}
const double rho_val = FluidProp->Density();
double arg[2];
arg[0] = interpolate(NodalVal1); // p
arg[1] = interpolate(NodalValC1); // T
drho_dp_rho = FluidProp->drhodP(arg) / rho_val;
}
else
{
drho_dp_rho = FluidProp->drho_dp;
}
val += poro_val * drho_dp_rho + (biot_val - poro_val) * (1.0 - biot_val) / SolidProp->K;

// Will handle the dual porosity version later...
}
else
{
poro_val = MediaProp->Porosity(Index, pcs->m_num->ls_theta);
const double poro_val = MediaProp->Porosity(Index, pcs->m_num->ls_theta);
val += poro_val * drho_dp_rho;
}

// AS,WX: 08.2012 storage function eff stress
if (MediaProp->storage_effstress_model > 0)
{
double storage_effstress = 1.;
CFiniteElementStd* h_fem;
h_fem = this;
storage_effstress = MediaProp->StorageFunctionEffStress(Index, nnodes, h_fem);
val *= storage_effstress;
}
// AS,WX: 08.2012 storage function eff stress
if (MediaProp->storage_effstress_model > 0)
{
double storage_effstress = 1.;
CFiniteElementStd* h_fem;
h_fem = this;
storage_effstress = MediaProp->StorageFunctionEffStress(Index, nnodes, h_fem);
val *= storage_effstress;
}

//val /= time_unit_factor;
//val /= time_unit_factor;
}
break;
case EPT_UNCONFINED_FLOW: // Unconfined flow
break;
Expand All @@ -1640,12 +1588,12 @@ double CFiniteElementStd::CalCoefMass()

// JT 2010, generalized poroelastic storage. See single phase version in case "L".
// Se = 1/M = poro/Kf + (alpha-poro)/Ks
rho_val = FluidProp->Density();
const double rho_val = FluidProp->Density();
if (D_Flag > 0 && rho_val > MKleinsteZahl)
{
biot_val = SolidProp->biot_const;
poro_val = MediaProp->Porosity(Index, pcs->m_num->ls_theta);
Se = poro_val * FluidProp->drho_dp + (biot_val - poro_val) * (1.0 - biot_val) / SolidProp->K;
const double biot_val = SolidProp->biot_const;
const double poro_val = MediaProp->Porosity(Index, pcs->m_num->ls_theta);
const double Se = poro_val * FluidProp->drho_dp + (biot_val - poro_val) * (1.0 - biot_val) / SolidProp->K;
// The poroelastic portion
val += Se * MMax(0., Sw);
}
Expand All @@ -1660,6 +1608,7 @@ double CFiniteElementStd::CalCoefMass()
// dSedPc always return positive numbers for default case
// However, the value should be negative analytically.
double dSwdPc = m_mmp->PressureSaturationDependency(Sw, true);
const double poro_val = MediaProp->Porosity(Index, pcs->m_num->ls_theta);
val += poro_val * dSwdPc;
}
}
Expand All @@ -1682,65 +1631,71 @@ double CFiniteElementStd::CalCoefMass()
break;
//....................................................................
case EPT_MASS_TRANSPORT: // Mass transport //SB4200
// Porosity
val = MediaProp->Porosity(Index, pcs->m_num->ls_theta);
// SB Transport in both phases
tr_phase = cp_vec[this->pcs->pcs_component_number]->transport_phase;
// Multi phase transport of components
saturation = PCSGetEleMeanNodeSecondary_2(Index, pcs->flow_pcs_type, "SATURATION1", 1);
if (tr_phase == 0) // Water phase
val *= saturation;
else if (tr_phase == 10) // non wetting phase
val *= (1.0 - saturation);
m_cp = cp_vec[pcs->pcs_component_number];
// Retardation Factor
val *= m_cp->CalcElementRetardationFactorNew(Index, unit, pcs);
{
// Porosity
val = MediaProp->Porosity(Index, pcs->m_num->ls_theta);
// SB Transport in both phases
const double tr_phase = cp_vec[this->pcs->pcs_component_number]->transport_phase;
// Multi phase transport of components
const double saturation = PCSGetEleMeanNodeSecondary_2(Index, pcs->flow_pcs_type, "SATURATION1", 1);
if (tr_phase == 0) // Water phase
val *= saturation;
else if (tr_phase == 10) // non wetting phase
val *= (1.0 - saturation);
CompProperties* const m_cp = cp_vec[pcs->pcs_component_number];
// Retardation Factor
val *= m_cp->CalcElementRetardationFactorNew(Index, unit, pcs);
}
break;
case EPT_OVERLAND_FLOW: // Liquid flow
val = 1.0;
break;
case EPT_RICHARDS_FLOW: // Richards
Sw = 1.0;
dSdp = 0.;
PG = interpolate(NodalVal1); // 12.02.2007. Important! WW

if (PG < 0.0) // JM skip to call these two functions in saturated case
{
Sw = MediaProp->SaturationCapillaryPressureFunction(-PG);
dSdp = -MediaProp->PressureSaturationDependency(Sw, true); // JT: dSdp now returns actual sign (i.e. <0)
}
Sw = 1.0;
dSdp = 0.;
PG = interpolate(NodalVal1); // 12.02.2007. Important! WW

poro = MediaProp->Porosity(Index, pcs->m_num->ls_theta);
rhow = FluidProp->Density();
if (PG < 0.0) // JM skip to call these two functions in saturated case
{
Sw = MediaProp->SaturationCapillaryPressureFunction(-PG);
dSdp = -MediaProp->PressureSaturationDependency(Sw, true); // JT: dSdp now returns actual sign (i.e. <0)
}

if (FluidProp->compressibility_model_pressure > 0)
{ // drho_dp from rho-p-T relation
arg[0] = PG;
arg[1] = interpolate(NodalValC1); // T
drho_dp_rho = FluidProp->drhodP(arg) / rhow;
}
else
drho_dp_rho = FluidProp->drho_dp;
poro = MediaProp->Porosity(Index, pcs->m_num->ls_theta);
rhow = FluidProp->Density();

// Storativity
val = MediaProp->StorageFunction(Index, unit, pcs->m_num->ls_theta) * Sw;
double drho_dp_rho = 0.0;
if (FluidProp->compressibility_model_pressure > 0)
{ // drho_dp from rho-p-T relation
double arg[2];
arg[0] = PG;
arg[1] = interpolate(NodalValC1); // T
drho_dp_rho = FluidProp->drhodP(arg) / rhow;
}
else
drho_dp_rho = FluidProp->drho_dp;

// Fluid compressibility
if (rhow > 0.0)
val += poro * Sw * drho_dp_rho;
// Capillarity
if (PG < 0.0) // dSdp gives always a value>0, even if p>0!
val += poro * dSdp;
// WW
if (MediaProp->heat_diffusion_model == 1)
{
// PG = fabs(interpolate(NodalVal1));
TG = interpolate(NodalValC) + PhysicalConstant::CelsiusZeroInKelvin;
humi = exp(PG / (SpecificGasConstant::WaterVapour * TG * rhow));
rhov = humi * FluidProp->vaporDensity(TG);
//
val -= poro * rhov * dSdp / rhow;
val += (1.0 - Sw) * poro * rhov / (rhow * rhow * SpecificGasConstant::WaterVapour * TG);
// Storativity
val = MediaProp->StorageFunction(Index, unit, pcs->m_num->ls_theta) * Sw;

// Fluid compressibility
if (rhow > 0.0)
val += poro * Sw * drho_dp_rho;
// Capillarity
if (PG < 0.0) // dSdp gives always a value>0, even if p>0!
val += poro * dSdp;
// WW
if (MediaProp->heat_diffusion_model == 1)
{
// PG = fabs(interpolate(NodalVal1));
TG = interpolate(NodalValC) + PhysicalConstant::CelsiusZeroInKelvin;
const double humi = exp(PG / (SpecificGasConstant::WaterVapour * TG * rhow));
const double rhov = humi * FluidProp->vaporDensity(TG);
//
val -= poro * rhov * dSdp / rhow;
val += (1.0 - Sw) * poro * rhov / (rhow * rhow * SpecificGasConstant::WaterVapour * TG);
}
}
break;
case EPT_FLUID_MOMENTUM: // Fluid Momentum
Expand Down
Loading

0 comments on commit 5c4b28c

Please sign in to comment.