OGS
ProcessLib::RichardsMechanics Namespace Reference

Classes

struct  IntegrationPointData
 
struct  LocalAssemblerInterface
 
struct  MicroPorosityParameters
 
struct  MicroPorosityStateSpace
 
class  RichardsMechanicsLocalAssembler
 
class  RichardsMechanicsProcess
 
struct  RichardsMechanicsProcessData
 
struct  SecondaryData
 

Functions

template<int DisplacementDim>
std::ostream & operator<< (std::ostream &os, MicroPorosityStateSpace< DisplacementDim > const &state)
 
template<int DisplacementDim>
MicroPorosityStateSpace< DisplacementDim > computeMicroPorosity (MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &I_2_C_el_inverse, double const rho_LR_m, double const mu_LR, MicroPorosityParameters const &micro_porosity_parameters, double const alpha_B, double const phi, double const p_L, double const p_L_m_prev, MaterialPropertyLib::VariableArray const &, double const S_L_m_prev, double const phi_m_prev, ParameterLib::SpatialPosition const pos, double const t, double const dt, MaterialPropertyLib::Property const &saturation_micro, MaterialPropertyLib::Property const &swelling_stress_rate)
 
void checkMPLProperties (std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media)
 
template<int DisplacementDim>
std::unique_ptr< ProcesscreateRichardsMechanicsProcess (std::string const &name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< ProcessVariable > const &variables, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, std::optional< ParameterLib::CoordinateSystem > const &local_coordinate_system, unsigned const integration_order, BaseLib::ConfigTree const &config, std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media)
 
template std::unique_ptr< ProcesscreateRichardsMechanicsProcess< 2 > (std::string const &name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< ProcessVariable > const &variables, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, std::optional< ParameterLib::CoordinateSystem > const &local_coordinate_system, unsigned const integration_order, BaseLib::ConfigTree const &config, std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media)
 
template std::unique_ptr< ProcesscreateRichardsMechanicsProcess< 3 > (std::string const &name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< ProcessVariable > const &variables, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, std::optional< ParameterLib::CoordinateSystem > const &local_coordinate_system, unsigned const integration_order, BaseLib::ConfigTree const &config, std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const &media)
 
template<int DisplacementDim, typename IPData >
void updateSwellingStressAndVolumetricStrain (IPData &ip_data, MaterialPropertyLib::Medium const &medium, MaterialPropertyLib::Phase const &solid_phase, MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > const &C_el, double const rho_LR, double const mu, std::optional< MicroPorosityParameters > micro_porosity_parameters, double const alpha, double const phi, double const p_cap_ip, MPL::VariableArray &variables, MPL::VariableArray &variables_prev, ParameterLib::SpatialPosition const &x_position, double const t, double const dt)
 

Function Documentation

◆ checkMPLProperties()

void ProcessLib::RichardsMechanics::checkMPLProperties ( std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const & media)

Definition at line 31 of file CreateRichardsMechanicsProcess.cpp.

33{
34 std::array const required_medium_properties = {
41 std::array const required_liquid_properties = {
43 std::array const required_solid_properties = {MaterialPropertyLib::density};
44
45 for (auto const& m : media)
46 {
47 checkRequiredProperties(*m.second, required_medium_properties);
48 checkRequiredProperties(m.second->phase("AqueousLiquid"),
49 required_liquid_properties);
50 checkRequiredProperties(m.second->phase("Solid"),
51 required_solid_properties);
52 }
53}
void checkRequiredProperties(Component const &c, std::span< PropertyType const > const required_properties)
Definition Component.cpp:60

References MaterialPropertyLib::biot_coefficient, MaterialPropertyLib::bishops_effective_stress, MaterialPropertyLib::density, MaterialPropertyLib::porosity, MaterialPropertyLib::reference_temperature, MaterialPropertyLib::relative_permeability, MaterialPropertyLib::saturation, and MaterialPropertyLib::viscosity.

Referenced by createRichardsMechanicsProcess().

◆ computeMicroPorosity()

template<int DisplacementDim>
MicroPorosityStateSpace< DisplacementDim > ProcessLib::RichardsMechanics::computeMicroPorosity ( MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const & I_2_C_el_inverse,
double const rho_LR_m,
double const mu_LR,
MicroPorosityParameters const & micro_porosity_parameters,
double const alpha_B,
double const phi,
double const p_L,
double const p_L_m_prev,
MaterialPropertyLib::VariableArray const & ,
double const S_L_m_prev,
double const phi_m_prev,
ParameterLib::SpatialPosition const pos,
double const t,
double const dt,
MaterialPropertyLib::Property const & saturation_micro,
MaterialPropertyLib::Property const & swelling_stress_rate )

Definition at line 64 of file ComputeMicroPorosity.h.

77{
78 namespace MPL = MaterialPropertyLib;
79 static constexpr int kelvin_vector_size =
81 // phi_m, e_sw, p_L_m, sigma_sw
82 static constexpr int nls_size = 1 + 1 + 1 + kelvin_vector_size;
83
84 static constexpr int i_phi_m = 0;
85 static constexpr int i_e_sw = 1;
86 static constexpr int i_p_L_m = 2;
87 static constexpr int i_sigma_sw = 3;
88
89 using ResidualVectorType = Eigen::Matrix<double, nls_size, 1>;
90 using JacobianMatrix =
91 Eigen::Matrix<double, nls_size, nls_size, Eigen::RowMajor>;
92
93 Eigen::FullPivLU<Eigen::Matrix<double, nls_size, nls_size, Eigen::RowMajor>>
94 linear_solver;
95
96 JacobianMatrix jacobian;
97
98 // Agglomerated solution vector construction.
99 ResidualVectorType solution = ResidualVectorType::Zero();
100
101 if (p_L >= 0 && p_L_m_prev >= 0)
102 {
103 return {solution[i_phi_m], solution[i_e_sw], solution[i_p_L_m],
104 solution.template segment<kelvin_vector_size>(i_sigma_sw)};
105 }
106
107 double const alpha_bar =
108 micro_porosity_parameters.mass_exchange_coefficient;
109
110 auto const update_residual = [&](ResidualVectorType& residual)
111 {
112 double const delta_phi_m = solution[i_phi_m];
113 double const delta_e_sw = solution[i_e_sw];
114 auto const& delta_sigma_sw =
115 solution.template segment<kelvin_vector_size>(i_sigma_sw);
116 double const delta_p_L_m = solution[i_p_L_m];
117
118 double const phi_m = phi_m_prev + delta_phi_m;
119 double const p_L_m = p_L_m_prev + delta_p_L_m;
120 MPL::VariableArray variables_prev;
121 variables_prev.capillary_pressure = -p_L_m_prev;
122 variables_prev.liquid_saturation = S_L_m_prev;
123
124 MPL::VariableArray variables;
125 variables.capillary_pressure = -p_L_m;
126
127 double const S_L_m =
128 saturation_micro.template value<double>(variables, pos, t, dt);
129 variables.liquid_saturation = S_L_m;
130 double const delta_S_L_m = S_L_m - S_L_m_prev;
131
132 auto const sigma_sw_dot =
133 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
134 swelling_stress_rate.template value<Eigen::Matrix3d>(
135 variables, variables_prev, pos, t, dt));
136
137 residual[i_phi_m] = delta_phi_m - (alpha_B - phi) * delta_e_sw;
138 residual[i_e_sw] = delta_e_sw + I_2_C_el_inverse.dot(delta_sigma_sw);
139 residual.template segment<kelvin_vector_size>(i_sigma_sw).noalias() =
140 delta_sigma_sw - sigma_sw_dot * dt;
141
142 residual[i_p_L_m] =
143 rho_LR_m *
144 (phi_m * delta_S_L_m - (alpha_B - phi) * S_L_m * delta_e_sw) +
145 phi_m * S_L_m * rho_LR_m * delta_e_sw -
146 micro_porosity_parameters.mass_exchange_coefficient / mu_LR *
147 (p_L - p_L_m) * dt;
148 };
149
150 auto const update_jacobian = [&](JacobianMatrix& jacobian)
151 {
152 jacobian = JacobianMatrix::Identity();
153
154 double const delta_phi_m = solution[i_phi_m];
155 double const delta_e_sw = solution[i_e_sw];
156 double const delta_p_L_m = solution[i_p_L_m];
157
158 double const phi_m = phi_m_prev + delta_phi_m;
159 double const p_L_m = p_L_m_prev + delta_p_L_m;
160 MPL::VariableArray variables_prev;
161 variables_prev.capillary_pressure = -p_L_m_prev;
162 MPL::VariableArray variables;
163 variables.capillary_pressure = -p_L_m;
164
165 double const S_L_m =
166 saturation_micro.template value<double>(variables, pos, t, dt);
167 variables_prev.liquid_saturation = S_L_m_prev;
168 variables.liquid_saturation = S_L_m;
169 double const delta_S_L_m = S_L_m - S_L_m_prev;
170
171 double const dS_L_m_dp_cap_m = saturation_micro.template dValue<double>(
172 variables, MPL::Variable::capillary_pressure, pos, t, dt);
173 auto const dsigma_sw_dS_L_m =
174 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
175 swelling_stress_rate.template dValue<Eigen::Matrix3d>(
176 variables, variables_prev, MPL::Variable::liquid_saturation,
177 pos, t, dt));
178
179 jacobian(i_phi_m, i_e_sw) = -(alpha_B - phi);
180
181 jacobian.template block<1, kelvin_vector_size>(i_e_sw, i_sigma_sw) =
182 I_2_C_el_inverse.transpose();
183
184 jacobian.template block<kelvin_vector_size, 1>(i_sigma_sw, i_p_L_m) =
185 -dsigma_sw_dS_L_m * dS_L_m_dp_cap_m;
186
187 jacobian(i_p_L_m, i_phi_m) =
188 rho_LR_m * (delta_S_L_m + S_L_m * delta_e_sw);
189
190 jacobian(i_p_L_m, i_e_sw) = -rho_LR_m * S_L_m * (alpha_B - phi - phi_m);
191
192 jacobian(i_p_L_m, i_p_L_m) =
193 alpha_bar / mu_LR * dt -
194 rho_LR_m * (phi_m - (alpha_B - phi - phi_m) * delta_e_sw) *
195 dS_L_m_dp_cap_m;
196 };
197
198 auto const update_solution = [&](ResidualVectorType const& increment)
199 { solution += increment; };
200
201 auto newton_solver =
202 NumLib::NewtonRaphson<decltype(linear_solver), JacobianMatrix,
203 decltype(update_jacobian), ResidualVectorType,
204 decltype(update_residual),
205 decltype(update_solution)>(
206 linear_solver, update_jacobian, update_residual, update_solution,
207 micro_porosity_parameters.nonlinear_solver_parameters);
208
209 auto const success_iterations = newton_solver.solve(jacobian);
210
211 if (!success_iterations)
212 {
213 OGS_FATAL(
214 "Could not find solution for local double structure nonlinear "
215 "problem.");
216 }
217
218 return {solution[i_phi_m], solution[i_e_sw], solution[i_p_L_m],
219 solution.template segment<kelvin_vector_size>(i_sigma_sw)};
220}
#define OGS_FATAL(...)
Definition Error.h:26
std::optional< int > solve(JacobianMatrix &jacobian) const
@ saturation_micro
capillary pressure saturation relationship for microstructure.
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.

References MaterialPropertyLib::VariableArray::capillary_pressure, MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::VariableArray::liquid_saturation, ProcessLib::RichardsMechanics::MicroPorosityParameters::mass_exchange_coefficient, ProcessLib::RichardsMechanics::MicroPorosityParameters::nonlinear_solver_parameters, OGS_FATAL, and NumLib::NewtonRaphson< LinearSolver, JacobianMatrix, JacobianMatrixUpdate, ResidualVector, ResidualUpdate, SolutionUpdate >::solve().

◆ createRichardsMechanicsProcess()

template<int DisplacementDim>
std::unique_ptr< Process > ProcessLib::RichardsMechanics::createRichardsMechanicsProcess ( std::string const & name,
MeshLib::Mesh & mesh,
std::unique_ptr< ProcessLib::AbstractJacobianAssembler > && jacobian_assembler,
std::vector< ProcessVariable > const & variables,
std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const & parameters,
std::optional< ParameterLib::CoordinateSystem > const & local_coordinate_system,
unsigned const integration_order,
BaseLib::ConfigTree const & config,
std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const & media )
Input File Parameter
prj__processes__process__type
Input File Parameter
prj__processes__process__RICHARDS_MECHANICS__coupling_scheme

Process Variables

Input File Parameter
prj__processes__process__RICHARDS_MECHANICS__process_variables

Primary process variables as they appear in the global component vector:

Input File Parameter
prj__processes__process__RICHARDS_MECHANICS__process_variables__pressure
Input File Parameter
prj__processes__process__RICHARDS_MECHANICS__process_variables__displacement

Process Parameters

Input File Parameter
prj__processes__process__RICHARDS_MECHANICS__specific_body_force
Input File Parameter
prj__processes__process__RICHARDS_MECHANICS__initial_stress
Input File Parameter
prj__processes__process__RICHARDS_MECHANICS__micro_porosity
Input File Parameter
prj__processes__process__RICHARDS_MECHANICS__micro_porosity__nonlinear_solver
Input File Parameter
prj__processes__process__RICHARDS_MECHANICS__micro_porosity__mass_exchange_coefficient
Input File Parameter
prj__processes__process__RICHARDS_MECHANICS__mass_lumping
Input File Parameter
prj__processes__process__RICHARDS_MECHANICS__explicit_hm_coupling_in_unsaturated_zone

Definition at line 56 of file CreateRichardsMechanicsProcess.cpp.

67{
69 config.checkConfigParameter("type", "RICHARDS_MECHANICS");
70 DBUG("Create RichardsMechanicsProcess.");
71
72 auto const coupling_scheme =
74 config.getConfigParameterOptional<std::string>("coupling_scheme");
75 const bool use_monolithic_scheme =
76 !(coupling_scheme && (*coupling_scheme == "staggered"));
77
79
81 auto const pv_config = config.getConfigSubtree("process_variables");
82
83 ProcessVariable* variable_p;
84 ProcessVariable* variable_u;
85 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>
86 process_variables;
87 if (use_monolithic_scheme) // monolithic scheme.
88 {
91 auto per_process_variables = findProcessVariables(
92 variables, pv_config,
93 {
94 "pressure",
96 "displacement"});
97 variable_p = &per_process_variables[0].get();
98 variable_u = &per_process_variables[1].get();
99 process_variables.push_back(std::move(per_process_variables));
100 }
101 else // staggered scheme.
102 {
103 using namespace std::string_literals;
104 for (auto const& variable_name : {"pressure"s, "displacement"s})
105 {
106 auto per_process_variables =
107 findProcessVariables(variables, pv_config, {variable_name});
108 process_variables.push_back(std::move(per_process_variables));
109 }
110 variable_p = &process_variables[0][0].get();
111 variable_u = &process_variables[1][0].get();
112 }
113
114 DBUG("Associate displacement with process variable '{:s}'.",
115 variable_u->getName());
116
117 if (variable_u->getNumberOfGlobalComponents() != DisplacementDim)
118 {
119 OGS_FATAL(
120 "Number of components of the process variable '{:s}' is different "
121 "from the displacement dimension: got {:d}, expected {:d}",
122 variable_u->getName(),
123 variable_u->getNumberOfGlobalComponents(),
124 DisplacementDim);
125 }
126
127 DBUG("Associate pressure with process variable '{:s}'.",
128 variable_p->getName());
129 if (variable_p->getNumberOfGlobalComponents() != 1)
130 {
131 OGS_FATAL(
132 "Pressure process variable '{:s}' is not a scalar variable but has "
133 "{:d} components.",
134 variable_p->getName(),
135 variable_p->getNumberOfGlobalComponents());
136 }
137
138 auto solid_constitutive_relations =
139 MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
140 parameters, local_coordinate_system, config);
141
143 // Specific body force
144 Eigen::Matrix<double, DisplacementDim, 1> specific_body_force;
145 {
146 std::vector<double> const b =
148 config.getConfigParameter<std::vector<double>>(
149 "specific_body_force");
150 if (b.size() != DisplacementDim)
151 {
152 OGS_FATAL(
153 "The size of the specific body force vector does not match the "
154 "displacement dimension. Vector size is {:d}, displacement "
155 "dimension is {:d}",
156 b.size(), DisplacementDim);
157 }
158
159 std::copy_n(b.data(), b.size(), specific_body_force.data());
160 }
161
162 auto media_map =
164 DBUG("Check the media properties of RichardsMechanics process ...");
165 checkMPLProperties(media);
166 DBUG("Media properties verified.");
167
168 // Initial stress conditions
169 auto const initial_stress = ParameterLib::findOptionalTagParameter<double>(
171 config, "initial_stress", parameters,
172 // Symmetric tensor size, 4 or 6, not a Kelvin vector.
174 &mesh);
175
176 std::optional<MicroPorosityParameters> micro_porosity_parameters;
177 if (auto const micro_porosity_config =
179 config.getConfigSubtreeOptional("micro_porosity"))
180 {
181 micro_porosity_parameters = MicroPorosityParameters{
184 micro_porosity_config->getConfigSubtree("nonlinear_solver")),
186 micro_porosity_config->getConfigParameter<double>(
187 "mass_exchange_coefficient")};
188 }
189
190 auto const mass_lumping =
192 config.getConfigParameter<bool>("mass_lumping", false);
193
194 auto const explicit_hm_coupling_in_unsaturated_zone =
196 config.getConfigParameter<bool>(
197 "explicit_hm_coupling_in_unsaturated_zone", false);
198
199 RichardsMechanicsProcessData<DisplacementDim> process_data{
200 materialIDs(mesh),
201 std::move(media_map),
202 std::move(solid_constitutive_relations),
203 initial_stress,
204 specific_body_force,
205 micro_porosity_parameters,
206 mass_lumping,
207 explicit_hm_coupling_in_unsaturated_zone};
208
209 SecondaryVariableCollection secondary_variables;
210
211 ProcessLib::createSecondaryVariables(config, secondary_variables);
212
213 return std::make_unique<RichardsMechanicsProcess<DisplacementDim>>(
214 std::move(name), mesh, std::move(jacobian_assembler), parameters,
215 integration_order, std::move(process_variables),
216 std::move(process_data), std::move(secondary_variables),
217 use_monolithic_scheme);
218}
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
std::string const & getName() const
int getNumberOfGlobalComponents() const
Returns the number of components of the process variable.
MaterialSpatialDistributionMap createMaterialSpatialDistributionMap(std::map< int, std::shared_ptr< Medium > > const &media, MeshLib::Mesh const &mesh)
PropertyVector< int > const * materialIDs(Mesh const &mesh)
Definition Mesh.cpp:268
NewtonRaphsonSolverParameters createNewtonRaphsonSolverParameters(BaseLib::ConfigTree const &config)
void checkMPLProperties(MeshLib::Mesh const &mesh, MaterialPropertyLib::MaterialSpatialDistributionMap const &media_map)
std::vector< std::reference_wrapper< ProcessVariable > > findProcessVariables(std::vector< ProcessVariable > const &variables, BaseLib::ConfigTree const &pv_config, std::initializer_list< std::string > tags)
void createSecondaryVariables(BaseLib::ConfigTree const &config, SecondaryVariableCollection &secondary_variables)

References BaseLib::ConfigTree::checkConfigParameter(), checkMPLProperties(), MaterialPropertyLib::createMaterialSpatialDistributionMap(), NumLib::createNewtonRaphsonSolverParameters(), ProcessLib::createSecondaryVariables(), DBUG(), ProcessLib::findProcessVariables(), BaseLib::ConfigTree::getConfigParameter(), BaseLib::ConfigTree::getConfigParameterOptional(), BaseLib::ConfigTree::getConfigSubtree(), BaseLib::ConfigTree::getConfigSubtreeOptional(), ProcessLib::ProcessVariable::getName(), ProcessLib::ProcessVariable::getNumberOfGlobalComponents(), MathLib::KelvinVector::kelvin_vector_dimensions(), and OGS_FATAL.

◆ createRichardsMechanicsProcess< 2 >()

template std::unique_ptr< Process > ProcessLib::RichardsMechanics::createRichardsMechanicsProcess< 2 > ( std::string const & name,
MeshLib::Mesh & mesh,
std::unique_ptr< ProcessLib::AbstractJacobianAssembler > && jacobian_assembler,
std::vector< ProcessVariable > const & variables,
std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const & parameters,
std::optional< ParameterLib::CoordinateSystem > const & local_coordinate_system,
unsigned const integration_order,
BaseLib::ConfigTree const & config,
std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const & media )

◆ createRichardsMechanicsProcess< 3 >()

template std::unique_ptr< Process > ProcessLib::RichardsMechanics::createRichardsMechanicsProcess< 3 > ( std::string const & name,
MeshLib::Mesh & mesh,
std::unique_ptr< ProcessLib::AbstractJacobianAssembler > && jacobian_assembler,
std::vector< ProcessVariable > const & variables,
std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const & parameters,
std::optional< ParameterLib::CoordinateSystem > const & local_coordinate_system,
unsigned const integration_order,
BaseLib::ConfigTree const & config,
std::map< int, std::shared_ptr< MaterialPropertyLib::Medium > > const & media )

◆ operator<<()

template<int DisplacementDim>
std::ostream & ProcessLib::RichardsMechanics::operator<< ( std::ostream & os,
MicroPorosityStateSpace< DisplacementDim > const & state )

◆ updateSwellingStressAndVolumetricStrain()

template<int DisplacementDim, typename IPData >
void ProcessLib::RichardsMechanics::updateSwellingStressAndVolumetricStrain ( IPData & ip_data,
MaterialPropertyLib::Medium const & medium,
MaterialPropertyLib::Phase const & solid_phase,
MathLib::KelvinVector::KelvinMatrixType< DisplacementDim > const & C_el,
double const rho_LR,
double const mu,
std::optional< MicroPorosityParameters > micro_porosity_parameters,
double const alpha,
double const phi,
double const p_cap_ip,
MPL::VariableArray & variables,
MPL::VariableArray & variables_prev,
ParameterLib::SpatialPosition const & x_position,
double const t,
double const dt )

Definition at line 33 of file RichardsMechanicsFEM-impl.h.

43{
44 auto const& identity2 = MathLib::KelvinVector::Invariants<
46 DisplacementDim)>::identity2;
47
48 auto& sigma_sw = ip_data.sigma_sw;
49 auto const& sigma_sw_prev = ip_data.sigma_sw_prev;
50 if (!medium.hasProperty(MPL::PropertyType::saturation_micro))
51 {
52 // If there is swelling, compute it. Update volumetric strain rate,
53 // s.t. it corresponds to the mechanical part only.
54 sigma_sw = sigma_sw_prev;
55 if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
56 {
57 auto const sigma_sw_dot =
58 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
60 solid_phase[MPL::PropertyType::swelling_stress_rate]
61 .value(variables, variables_prev, x_position, t,
62 dt)));
63 sigma_sw += sigma_sw_dot * dt;
64
65 // !!! Misusing volumetric strain for mechanical volumetric
66 // strain just to update the transport porosity !!!
67 variables.volumetric_strain +=
68 identity2.transpose() * C_el.inverse() * sigma_sw;
69 variables_prev.volumetric_strain +=
70 identity2.transpose() * C_el.inverse() * sigma_sw_prev;
71 }
72 }
73
74 // TODO (naumov) saturation_micro must be always defined together with
75 // the micro_porosity_parameters.
76 if (medium.hasProperty(MPL::PropertyType::saturation_micro))
77 {
78 double const phi_M_prev = ip_data.transport_porosity_prev;
79 double const phi_prev = ip_data.porosity_prev;
80 double const phi_m_prev = phi_prev - phi_M_prev;
81 double const p_L_m_prev = ip_data.liquid_pressure_m_prev;
82
83 auto const S_L_m_prev = ip_data.saturation_m_prev;
84
85 auto const [delta_phi_m, delta_e_sw, delta_p_L_m, delta_sigma_sw] =
86 computeMicroPorosity<DisplacementDim>(
87 identity2.transpose() * C_el.inverse(), rho_LR, mu,
88 *micro_porosity_parameters, alpha, phi, -p_cap_ip, p_L_m_prev,
89 variables_prev, S_L_m_prev, phi_m_prev, x_position, t, dt,
90 medium.property(MPL::PropertyType::saturation_micro),
91 solid_phase.property(MPL::PropertyType::swelling_stress_rate));
92
93 auto& phi_M = ip_data.transport_porosity;
94 phi_M = phi - (phi_m_prev + delta_phi_m);
95 variables_prev.transport_porosity = phi_M_prev;
96 variables.transport_porosity = phi_M;
97
98 auto& p_L_m = ip_data.liquid_pressure_m;
99 p_L_m = p_L_m_prev + delta_p_L_m;
100 { // Update micro saturation.
101 MPL::VariableArray variables_prev;
102 variables_prev.capillary_pressure = -p_L_m_prev;
103 MPL::VariableArray variables;
104 variables.capillary_pressure = -p_L_m;
105
106 ip_data.saturation_m =
107 medium.property(MPL::PropertyType::saturation_micro)
108 .template value<double>(variables, x_position, t, dt);
109 }
110 sigma_sw = sigma_sw_prev + delta_sigma_sw;
111 }
112}
template Eigen::Matrix< double, 3, 3 > formEigenTensor< 3 >(MaterialPropertyLib::PropertyDataType const &values)

References MaterialPropertyLib::VariableArray::capillary_pressure, MaterialPropertyLib::formEigenTensor< 3 >(), MaterialPropertyLib::Medium::hasProperty(), MaterialPropertyLib::Phase::hasProperty(), MathLib::KelvinVector::kelvin_vector_dimensions(), MaterialPropertyLib::Medium::property(), MaterialPropertyLib::Phase::property(), MaterialPropertyLib::VariableArray::transport_porosity, and MaterialPropertyLib::VariableArray::volumetric_strain.