22 template <
int DisplacementDim>
42 template <
int DisplacementDim>
46 return os <<
"phi_m: " << state.
phi_m <<
", "
47 <<
"e_sw: " << state.
e_sw <<
", "
48 <<
"p_L_m: " << state.
p_L_m <<
", "
49 <<
"sigma_sw: " << state.
sigma_sw.transpose();
62 template <
int DisplacementDim>
66 double const rho_LR_m,
69 double const alpha_B,
double const phi,
double const p_L,
70 double const p_L_m_prev,
72 double const S_L_m_prev,
double const phi_m_prev,
78 static constexpr
int kelvin_vector_size =
81 static constexpr
int nls_size = 1 + 1 + 1 + kelvin_vector_size;
83 static constexpr
int i_phi_m = 0;
84 static constexpr
int i_e_sw = 1;
85 static constexpr
int i_p_L_m = 2;
86 static constexpr
int i_sigma_sw = 3;
88 using ResidualVectorType = Eigen::Matrix<double, nls_size, 1>;
89 using JacobianMatrix =
90 Eigen::Matrix<double, nls_size, nls_size, Eigen::RowMajor>;
92 Eigen::FullPivLU<Eigen::Matrix<double, nls_size, nls_size, Eigen::RowMajor>>
95 JacobianMatrix jacobian;
98 ResidualVectorType solution = ResidualVectorType::Zero();
100 if (p_L >= 0 && p_L_m_prev >= 0)
102 return {solution[i_phi_m], solution[i_e_sw], solution[i_p_L_m],
103 solution.template segment<kelvin_vector_size>(i_sigma_sw)};
106 double const alpha_bar =
109 auto const update_residual = [&](ResidualVectorType& residual)
111 double const delta_phi_m = solution[i_phi_m];
112 double const delta_e_sw = solution[i_e_sw];
113 auto const& delta_sigma_sw =
114 solution.template segment<kelvin_vector_size>(i_sigma_sw);
115 double const delta_p_L_m = solution[i_p_L_m];
117 double const phi_m = phi_m_prev + delta_phi_m;
118 double const p_L_m = p_L_m_prev + delta_p_L_m;
122 variables_prev[
static_cast<int>(MPL::Variable::liquid_saturation)] =
130 variables[
static_cast<int>(MPL::Variable::liquid_saturation)] = S_L_m;
131 double const delta_S_L_m = S_L_m - S_L_m_prev;
133 auto const sigma_sw_dot =
134 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
136 variables, variables_prev, pos, t, dt));
138 residual[i_phi_m] = delta_phi_m - (alpha_B - phi) * delta_e_sw;
139 residual[i_e_sw] = delta_e_sw + I_2_C_el_inverse.dot(delta_sigma_sw);
140 residual.template segment<kelvin_vector_size>(i_sigma_sw).noalias() =
141 delta_sigma_sw - sigma_sw_dot * dt;
145 (phi_m * delta_S_L_m - (alpha_B - phi) * S_L_m * delta_e_sw) +
146 phi_m * S_L_m * rho_LR_m * delta_e_sw -
151 auto const update_jacobian = [&](JacobianMatrix& jacobian)
153 jacobian = JacobianMatrix::Identity();
155 double const delta_phi_m = solution[i_phi_m];
156 double const delta_e_sw = solution[i_e_sw];
157 double const delta_p_L_m = solution[i_p_L_m];
159 double const phi_m = phi_m_prev + delta_phi_m;
160 double const p_L_m = p_L_m_prev + delta_p_L_m;
169 variables_prev[
static_cast<int>(MPL::Variable::liquid_saturation)] =
171 variables[
static_cast<int>(MPL::Variable::liquid_saturation)] = S_L_m;
172 double const delta_S_L_m = S_L_m - S_L_m_prev;
176 auto const dsigma_sw_dS_L_m =
177 MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
179 variables, variables_prev, MPL::Variable::liquid_saturation,
182 jacobian(i_phi_m, i_e_sw) = -(alpha_B - phi);
184 jacobian.template block<1, kelvin_vector_size>(i_e_sw, i_sigma_sw) =
185 I_2_C_el_inverse.transpose();
187 jacobian.template block<kelvin_vector_size, 1>(i_sigma_sw, i_p_L_m) =
188 -dsigma_sw_dS_L_m * dS_L_m_dp_cap_m;
190 jacobian(i_p_L_m, i_phi_m) =
191 rho_LR_m * (delta_S_L_m + S_L_m * delta_e_sw);
193 jacobian(i_p_L_m, i_e_sw) = -rho_LR_m * S_L_m * (alpha_B - phi - phi_m);
195 jacobian(i_p_L_m, i_p_L_m) =
196 alpha_bar / mu_LR * dt -
197 rho_LR_m * (phi_m - (alpha_B - phi - phi_m) * delta_e_sw) *
201 auto const update_solution = [&](ResidualVectorType
const& increment)
202 { solution += increment; };
206 decltype(update_jacobian), ResidualVectorType,
207 decltype(update_residual),
208 decltype(update_solution)>(
209 linear_solver, update_jacobian, update_residual, update_solution,
212 auto const success_iterations = newton_solver.
solve(jacobian);
214 if (!success_iterations)
217 "Could not find solution for local double structure nonlinear "
221 return {solution[i_phi_m], solution[i_e_sw], solution[i_p_L_m],
222 solution.template segment<kelvin_vector_size>(i_sigma_sw)};
std::optional< int > solve(JacobianMatrix &jacobian) const
std::array< VariableType, static_cast< int >(Variable::number_of_variables)> VariableArray
@ saturation_micro
capillary pressure saturation relationship for microstructure.
Eigen::Matrix< double, kelvin_vector_dimensions(DisplacementDim), 1, Eigen::ColMajor > KelvinVectorType
constexpr int kelvin_vector_dimensions(int const displacement_dim)
Kelvin vector dimensions for given displacement dimension.
MicroPorosityStateSpace< DisplacementDim > computeMicroPorosity(MathLib::KelvinVector::KelvinVectorType< DisplacementDim > const &I_2_C_el_inverse, double const rho_LR_m, double const mu_LR, MicroPorosityParameters const µ_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)
std::ostream & operator<<(std::ostream &os, MicroPorosityStateSpace< DisplacementDim > const &state)
double mass_exchange_coefficient
NumLib::NewtonRaphsonSolverParameters nonlinear_solver_parameters
MicroPorosityStateSpace & operator+=(MicroPorosityStateSpace const &state)
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
MathLib::KelvinVector::KelvinVectorType< DisplacementDim > sigma_sw