From 8af62535f21a353947e7e96e290f2d2c179f2cf1 Mon Sep 17 00:00:00 2001 From: Wenqing Wang Date: Tue, 2 Oct 2018 12:17:46 +0200 Subject: [PATCH] Fixed storage mode 7 --- FEM/fem_ele_std.cpp | 99 ++++++++++++--------------------------------- FEM/rf_mmp_new.cpp | 30 ++++++++------ FEM/rf_msp_new.cpp | 32 +++++++++++++++ FEM/rf_msp_new.h | 6 ++- 4 files changed, 79 insertions(+), 88 deletions(-) diff --git a/FEM/fem_ele_std.cpp b/FEM/fem_ele_std.cpp index 4e538bdc4..7d228a9ec 100644 --- a/FEM/fem_ele_std.cpp +++ b/FEM/fem_ele_std.cpp @@ -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; diff --git a/FEM/rf_mmp_new.cpp b/FEM/rf_mmp_new.cpp index 8f774ff8a..fc914616f 100644 --- a/FEM/rf_mmp_new.cpp +++ b/FEM/rf_mmp_new.cpp @@ -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"; @@ -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()); diff --git a/FEM/rf_msp_new.cpp b/FEM/rf_msp_new.cpp index e1ef780c8..1532d8c5c 100644 --- a/FEM/rf_msp_new.cpp +++ b/FEM/rf_msp_new.cpp @@ -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 ///////////////////////////////////////////////////////////////////////////// diff --git a/FEM/rf_msp_new.h b/FEM/rf_msp_new.h index ca1ce8ec0..5f234b4e8 100644 --- a/FEM/rf_msp_new.h +++ b/FEM/rf_msp_new.h @@ -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(); @@ -166,6 +164,10 @@ class CSolidProperties std::vector& eps_K_curr, std::vector& eps_M_curr, std::vector& 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;