1060{
1062
1063 auto const [p_L, u] =
localDOF(local_x);
1064 auto const [p_L_prev, u_prev] =
localDOF(local_x_prev);
1065
1067 typename ShapeMatricesTypeDisplacement::template MatrixType<
1072
1074 typename ShapeMatricesTypeDisplacement::template VectorType<
1077
1080 DisplacementDim)>::identity2;
1081
1083 ShapeMatricesTypePressure::NodalMatrixType::Zero(
pressure_size,
1085
1087 ShapeMatricesTypePressure::NodalMatrixType::Zero(
pressure_size,
1089
1091 ShapeMatricesTypePressure::NodalMatrixType::Zero(
pressure_size,
1093
1095 ShapeMatricesTypePressure::NodalMatrixType::Zero(
pressure_size,
1097
1098 typename ShapeMatricesTypeDisplacement::template MatrixType<
1100 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
1103
1104 typename ShapeMatricesTypeDisplacement::template MatrixType<
1106 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
1109
1110 auto const& medium =
1112 auto const& liquid_phase = medium->phase("AqueousLiquid");
1113 auto const& solid_phase = medium->phase("Solid");
1116
1119
1120 unsigned const n_integration_points =
1122 for (unsigned ip = 0; ip < n_integration_points; ip++)
1123 {
1129
1130 auto const& w =
_ip_data[ip].integration_weight;
1131
1132 auto const& N_u =
_ip_data[ip].N_u;
1133 auto const& dNdx_u =
_ip_data[ip].dNdx_u;
1134
1135 auto const& N_p =
_ip_data[ip].N_p;
1136 auto const& dNdx_p =
_ip_data[ip].dNdx_p;
1137
1138 auto const x_coord =
1142 auto const B =
1144 ShapeFunctionDisplacement::NPOINTS,
1147
1148 double p_cap_ip;
1150
1151 double p_cap_prev_ip;
1153
1156
1157
1159
1161 medium->property(MPL::PropertyType::reference_temperature)
1162 .template value<double>(variables, x_position, t, dt);
1164
1165 std::get<StrainData<DisplacementDim>>(SD).eps.noalias() = B * u;
1166
1168 t, dt, x_position,
_ip_data[ip], variables, variables_prev, medium,
1170 CapillaryPressureData<DisplacementDim>{
1171 p_cap_ip, p_cap_prev_ip,
1172 Eigen::Vector<double, DisplacementDim>::Zero()},
1173 CD, SD, SD_prev, this->
process_data_.micro_porosity_parameters,
1175
1176 {
1177 auto const&
C = *std::get<StiffnessTensor<DisplacementDim>>(CD);
1178 local_Jac
1179 .template block<displacement_size, displacement_size>(
1181 .noalias() += B.transpose() *
C * B * w;
1182 }
1183
1185
1186 {
1187 auto const& sigma_eff =
1190 .sigma_eff;
1191 double const rho = *std::get<Density>(CD);
1193 .noalias() -= (B.transpose() * sigma_eff -
1195 w;
1196 }
1197
1198
1199
1200
1201
1202 double const alpha =
1203 *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD);
1204 double const dS_L_dp_cap =
1205 std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(
1206 CD)
1207 .dS_L_dp_cap;
1208
1209 {
1210 double const chi_S_L =
1211 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
1212 .chi_S_L;
1213 Kup.noalias() +=
1214 B.transpose() *
alpha * chi_S_L * identity2 * N_p * w;
1215 double const dchi_dS_L =
1216 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
1217 .dchi_dS_L;
1218
1219 local_Jac
1220 .template block<displacement_size, pressure_size>(
1222 .noalias() -= B.transpose() *
alpha *
1223 (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
1224 identity2 * N_p * w;
1225 }
1226
1227 double const phi =
1228 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi;
1229 double const rho_LR = *std::get<LiquidDensity>(CD);
1230 local_Jac
1231 .template block<displacement_size, pressure_size>(
1233 .noalias() +=
1234 N_u_op(N_u).transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245 if (!medium->hasProperty(MPL::PropertyType::saturation_micro) &&
1246 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
1247 {
1248 using DimMatrix = Eigen::Matrix<double, 3, 3>;
1249 auto const dsigma_sw_dS_L =
1251 solid_phase
1252 .property(MPL::PropertyType::swelling_stress_rate)
1253 .template dValue<DimMatrix>(
1254 variables, variables_prev,
1255 MPL::Variable::liquid_saturation, x_position, t,
1256 dt));
1257 local_Jac
1258 .template block<displacement_size, pressure_size>(
1260 .noalias() +=
1261 B.transpose() * dsigma_sw_dS_L * dS_L_dp_cap * N_p * w;
1262 }
1263
1264
1265
1266 double const S_L =
1267 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
1269 .S_L;
1270 if (this->
process_data_.explicit_hm_coupling_in_unsaturated_zone)
1271 {
1272 double const chi_S_L_prev = std::get<PrevState<
1274 ->chi_S_L;
1275 Kpu.noalias() += N_p.transpose() * chi_S_L_prev * rho_LR *
alpha *
1276 identity2.transpose() * B * w;
1277 }
1278 else
1279 {
1280 Kpu.noalias() += N_p.transpose() * S_L * rho_LR *
alpha *
1281 identity2.transpose() * B * w;
1282 }
1283
1284
1285
1286
1287
1288 double const k_rel =
1290 DisplacementDim>>(CD)
1291 .k_rel;
1292 auto const& K_intrinsic =
1294 DisplacementDim>>(CD)
1295 .Ki;
1296 double const mu =
1297 *std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(
1298 CD);
1299
1301
1302 laplace_p.noalias() +=
1303 dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
1304
1305 auto const beta_LR =
1306 1 / rho_LR *
1307 liquid_phase.property(MPL::PropertyType::density)
1308 .template dValue<double>(variables,
1309 MPL::Variable::liquid_phase_pressure,
1310 x_position, t, dt);
1311
1312 double const beta_SR =
1313 std::get<
1315 CD)
1316 .beta_SR;
1317 double const a0 = (
alpha - phi) * beta_SR;
1318 double const specific_storage_a_p = S_L * (phi * beta_LR + S_L * a0);
1319 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
1320
1321 double const dspecific_storage_a_p_dp_cap =
1322 dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0);
1323 double const dspecific_storage_a_S_dp_cap =
1324 -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
1325
1326 storage_p_a_p.noalias() +=
1327 N_p.transpose() * rho_LR * specific_storage_a_p * N_p * w;
1328
1329 double const DeltaS_L_Deltap_cap =
1330 std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap;
1331 storage_p_a_S.noalias() -= N_p.transpose() * rho_LR *
1332 specific_storage_a_S * DeltaS_L_Deltap_cap *
1333 N_p * w;
1334
1335 local_Jac
1338 .noalias() += N_p.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
1339 rho_LR * dspecific_storage_a_p_dp_cap * N_p * w;
1340
1341 double const S_L_prev =
1342 std::get<
1343 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
1345 ->S_L;
1346 storage_p_a_S_Jpp.noalias() -=
1347 N_p.transpose() * rho_LR *
1348 ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
1349 specific_storage_a_S * dS_L_dp_cap) /
1350 dt * N_p * w;
1351
1352 if (!this->
process_data_.explicit_hm_coupling_in_unsaturated_zone)
1353 {
1354 local_Jac
1357 .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap *
alpha *
1358 identity2.transpose() * B * (u - u_prev) / dt *
1359 N_p * w;
1360 }
1361
1362 double const dk_rel_dS_l =
1363 medium->property(MPL::PropertyType::relative_permeability)
1364 .template dValue<double>(variables,
1365 MPL::Variable::liquid_saturation,
1366 x_position, t, dt);
1368 grad_p_cap = -dNdx_p * p_L;
1369 local_Jac
1372 .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
1373 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1374
1375 local_Jac
1378 .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
1379 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1380
1381 local_rhs.template segment<pressure_size>(
pressure_index).noalias() +=
1382 dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
1383
1384 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
1385 {
1386 double const alpha_bar =
1388 ->mass_exchange_coefficient;
1389 auto const p_L_m =
1392 .noalias() -=
1393 N_p.transpose() * alpha_bar / mu * (-p_cap_ip - p_L_m) * w;
1394
1395 local_Jac
1398 .noalias() += N_p.transpose() * alpha_bar / mu * N_p * w;
1399 if (p_cap_ip != p_cap_prev_ip)
1400 {
1401 auto const p_L_m_prev = **std::get<PrevState<MicroPressure>>(
1403 local_Jac
1404 .template block<pressure_size, pressure_size>(
1406 .noalias() += N_p.transpose() * alpha_bar / mu *
1407 (p_L_m - p_L_m_prev) /
1408 (p_cap_ip - p_cap_prev_ip) * N_p * w;
1409 }
1410 }
1411 }
1412
1414 {
1415 storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
1416 storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
1417 storage_p_a_S_Jpp =
1418 storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
1419 }
1420
1421
1422 local_Jac
1425 .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
1426
1427
1428 local_Jac
1429 .template block<pressure_size, displacement_size>(
pressure_index,
1431 .noalias() = Kpu / dt;
1432
1433
1434 local_rhs.template segment<pressure_size>(
pressure_index).noalias() -=
1435 laplace_p * p_L +
1436 (storage_p_a_p + storage_p_a_S) * (p_L - p_L_prev) / dt +
1437 Kpu * (u - u_prev) / dt;
1438
1439
1441 .noalias() += Kup * p_L;
1442}
static void assembleWithJacobianEvalConstitutiveSetting(double const t, double const dt, ParameterLib::SpatialPosition const &x_position, IpData &ip_data, MPL::VariableArray &variables, MPL::VariableArray &variables_prev, MPL::Medium const *const medium, TemperatureData const T_data, CapillaryPressureData< DisplacementDim > const &p_cap_data, ConstitutiveData< DisplacementDim > &CD, StatefulData< DisplacementDim > &SD, StatefulDataPrev< DisplacementDim > const &SD_prev, std::optional< MicroPorosityParameters > const µ_porosity_parameters, MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material, ProcessLib::ThermoRichardsMechanics::MaterialStateData< DisplacementDim > &material_state_data)
BaseLib::StrongType< double, struct TemperatureDataTag > TemperatureData
ConstitutiveModels< DisplacementDim > createConstitutiveModels(TRMProcessData const &process_data, MaterialLib::Solids::MechanicsBase< DisplacementDim > const &solid_material)
std::tuple< StiffnessTensor< DisplacementDim >, ProcessLib::ThermoRichardsMechanics::PorosityData, Density, LiquidDensity, ProcessLib::ThermoRichardsMechanics::BiotData, ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv, ProcessLib::ThermoRichardsMechanics::LiquidViscosityData, ProcessLib::ThermoRichardsMechanics::SolidCompressibilityData, ProcessLib::ThermoRichardsMechanics::BishopsData, PrevState< ProcessLib::ThermoRichardsMechanics::BishopsData >, ProcessLib::ThermoRichardsMechanics::PermeabilityData< DisplacementDim >, SaturationSecantDerivative > ConstitutiveData
Data that is needed for the equation system assembly.
MatrixType< ShapeFunction::NPOINTS, ShapeFunction::NPOINTS > NodalMatrixType
VectorType< GlobalDim > GlobalDimVectorType
Public Member Functions inherited from ProcessLib::RichardsMechanics::LocalAssemblerInterface< DisplacementDim >
1.11.0