Skip to content

Commit

Permalink
Fixed storage mode 7
Browse files Browse the repository at this point in the history
  • Loading branch information
@wenqing
wenqing committed Oct 2, 2018
1 parent 4b4647b commit 8af6253
Show file tree
Hide file tree
Showing 4 changed files with 79 additions and 88 deletions.
99 changes: 25 additions & 74 deletions FEM/fem_ele_std.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1583,87 +1583,38 @@ double CFiniteElementStd::CalCoefMass()
<< "\n";
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
#ifdef USE_FREESTEAM
if (FluidProp->compressibility_model_pressure > 0 || FluidProp->density_model == 8)
#else
if (FluidProp->compressibility_model_pressure > 0)
#endif
{
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 == 8) // 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

SolidProp->Calculate_Lame_Constant();
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;
}
// Ks = K / (1-alpha_B)
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 Down
30 changes: 18 additions & 12 deletions FEM/rf_mmp_new.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -812,16 +812,13 @@ std::ios::pos_type CMediumProperties::Read(std::ifstream* mmp_file)
break;
case 7: // RW/WW
{
in >> storage_model_values[0]; // Biot's alpha
in >> storage_model_values[1]; // Skempton's B coefficient
in >> storage_model_values[2]; // macroscopic drained bulk modulus
double val_l = storage_model_values[0] * (1. - storage_model_values[0] * storage_model_values[1])
/ storage_model_values[1] / storage_model_values[2];
storage_model_values[1] = val_l;
if(!PCSGet("DEFORMATION"))
{
ScreenMessage("Error: Porosity model 7 must be combined with deformation process");
exit(EXIT_FAILURE);
}
break;
}
case 8:
break;
default:
cout << "Error in MMPRead: no valid storativity model"
<< "\n";
Expand Down Expand Up @@ -7455,10 +7452,19 @@ double CMediumProperties::StorageFunction(long index, double* gp, double theta)
*/
break;
case 7: // poroelasticity RW
storage = storage_model_values[1];
break;
case 8:
return 0.0;
{
// Moved the following comment from double CFiniteElementStd::CalCoefMass() to here by WW
// 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

SolidProp::CSolidProperties* solid_prop = msp_vector[Fem_Ele_Std->GetMeshElement()->GetPatchIndex()];
const double biots_constant = solid_prop->getBiotsConstant();
const double porosity = Porosity(index, theta);
return (biots_constant - porosity) * (1.0 - biots_constant) / solid_prop->getBulkModulus();
}
case 10:
if (permeability_saturation_model[0] == 10) // MW
storage = porosity_model_values[0] / (gravity_constant * gravity_constant * mfp_vector[0]->Density());
Expand Down
32 changes: 32 additions & 0 deletions FEM/rf_msp_new.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8207,6 +8207,38 @@ double CSolidProperties::getBishopCoefficient(const double effectiveS, const dou
}
return p;
}

double CSolidProperties::getBulkModulus() const
{
if (K > DBL_MIN)
return K;

if (Youngs_mode < 10 || Youngs_mode > 13)
return E / 3 / (1 - 2 * PoissonRatio);

// average Youngs modulus
double const E_av = 2. / 3. * (*data_Youngs)(0) + 1. / 3. * (*data_Youngs)(1);

// Poisson ratio perpendicular to the plane of isotropie, due to strain in the
// plane of isotropie
const double nu_ia = (*data_Youngs)(2);

// Poisson ratio in the plane of isotropie, due to strain perpendicular to the
// plane of isotropie
const double nu_ai = nu_ia * (*data_Youngs)(1)
/ (*data_Youngs)(0); // nu_ai=nu_ia*Ea/Ei

// Poisson ratio in the plane of isotropy
const double nu_i = Poisson_Ratio();

// average Poisson ratio
// 12 13 21 23 31 32
// ai ai ia ii ia ii
const double nu_av = 1. / 3. * (nu_ai + nu_ia + nu_i);

return E_av / 3 / (1 - 2 * nu_av);
}

} // end namespace

/////////////////////////////////////////////////////////////////////////////
Expand Down
6 changes: 4 additions & 2 deletions FEM/rf_msp_new.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,8 +93,6 @@ class CSolidProperties
int GetConductModel() const { return Conductivity_mode; }
double Thermal_Expansion() const { return ThermalExpansion; }
double Poisson_Ratio() const { return PoissonRatio; }
double getBulkModulus() const { return K; }
double getBiotsConstant() const { return biot_const; }

// Elasticity
void Calculate_Lame_Constant();
Expand Down Expand Up @@ -166,6 +164,10 @@ class CSolidProperties
std::vector<double>& eps_K_curr, std::vector<double>& eps_M_curr,
std::vector<double>& eps_pl_curr, double& e_pl_v, double& e_pl_eff, double& lam,
Math_Group::Matrix& Consistent_Tangent, double Temperature, double& local_res);

double getBulkModulus() const;
double getBiotsConstant() const { return biot_const; }

private:
// CMCD
FiniteElement::CFiniteElementStd* Fem_Ele_Std;
Expand Down

0 comments on commit 8af6253

Please sign in to comment.