1101{
1103
1104 auto const [p_L, u] =
localDOF(local_x);
1105 auto const [p_L_prev, u_prev] =
localDOF(local_x_prev);
1106
1108 typename ShapeMatricesTypeDisplacement::template MatrixType<
1113
1115 typename ShapeMatricesTypeDisplacement::template VectorType<
1118
1121 DisplacementDim)>::identity2;
1122
1124 ShapeMatricesTypePressure::NodalMatrixType::Zero(
pressure_size,
1126
1128 ShapeMatricesTypePressure::NodalMatrixType::Zero(
pressure_size,
1130
1132 ShapeMatricesTypePressure::NodalMatrixType::Zero(
pressure_size,
1134
1136 ShapeMatricesTypePressure::NodalMatrixType::Zero(
pressure_size,
1138
1139 typename ShapeMatricesTypeDisplacement::template MatrixType<
1141 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
1144
1145 typename ShapeMatricesTypeDisplacement::template MatrixType<
1147 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
1150
1151 auto const& medium =
1153 auto const& liquid_phase = medium->phase("AqueousLiquid");
1154 auto const& solid_phase = medium->phase("Solid");
1157
1158 unsigned const n_integration_points =
1160 for (unsigned ip = 0; ip < n_integration_points; ip++)
1161 {
1167
1168 auto const& w =
ip_data_[ip].integration_weight;
1169
1170 auto const& N_u =
ip_data_[ip].N_u;
1171 auto const& dNdx_u =
ip_data_[ip].dNdx_u;
1172
1173 auto const& N_p =
ip_data_[ip].N_p;
1174 auto const& dNdx_p =
ip_data_[ip].dNdx_p;
1175
1183
1184 auto const B =
1186 ShapeFunctionDisplacement::NPOINTS,
1189
1190 double p_cap_ip;
1192
1193 double p_cap_prev_ip;
1195
1198
1199
1201
1203 medium->property(MPL::PropertyType::reference_temperature)
1204 .template value<double>(variables, x_position, t, dt);
1206
1207 std::get<StrainData<DisplacementDim>>(SD).eps.noalias() = B * u;
1208
1210 t, dt, x_position,
ip_data_[ip], variables, variables_prev, medium,
1212 CapillaryPressureData<DisplacementDim>{
1213 p_cap_ip, p_cap_prev_ip,
1214 Eigen::Vector<double, DisplacementDim>::Zero()},
1215 CD, SD, SD_prev, this->
process_data_.micro_porosity_parameters,
1217
1218 {
1219 auto const& C = *std::get<StiffnessTensor<DisplacementDim>>(CD);
1220 local_Jac
1221 .template block<displacement_size, displacement_size>(
1223 .noalias() += B.transpose() * C * B * w;
1224 }
1225
1227
1228 {
1229 auto const& sigma_eff =
1232 .sigma_eff;
1233 double const rho = *std::get<Density>(CD);
1235 .noalias() -= (B.transpose() * sigma_eff -
1236 N_u_op(N_u).transpose() * rho * b) *
1237 w;
1238 }
1239
1240
1241
1242
1243
1244 double const alpha =
1245 *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD);
1246 double const dS_L_dp_cap =
1247 std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(
1248 CD)
1249 .dS_L_dp_cap;
1250
1251 {
1252 double const chi_S_L =
1253 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
1254 .chi_S_L;
1255 Kup.noalias() +=
1256 B.transpose() *
alpha * chi_S_L * identity2 * N_p * w;
1257 double const dchi_dS_L =
1258 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
1259 .dchi_dS_L;
1260
1261 local_Jac
1262 .template block<displacement_size, pressure_size>(
1264 .noalias() -= B.transpose() *
alpha *
1265 (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
1266 identity2 * N_p * w;
1267 }
1268
1269 double const phi =
1270 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi;
1271 double const rho_LR = *std::get<LiquidDensity>(CD);
1272 local_Jac
1273 .template block<displacement_size, pressure_size>(
1275 .noalias() +=
1276 N_u_op(N_u).transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287 if (!medium->hasProperty(MPL::PropertyType::saturation_micro) &&
1288 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
1289 {
1290 using DimMatrix = Eigen::Matrix<double, 3, 3>;
1291 auto const dsigma_sw_dS_L =
1293 solid_phase
1294 .property(MPL::PropertyType::swelling_stress_rate)
1295 .template dValue<DimMatrix>(
1296 variables, variables_prev,
1297 MPL::Variable::liquid_saturation, x_position, t,
1298 dt));
1299 local_Jac
1300 .template block<displacement_size, pressure_size>(
1302 .noalias() +=
1303 B.transpose() * dsigma_sw_dS_L * dS_L_dp_cap * N_p * w;
1304 }
1305
1306
1307
1308 double const S_L =
1309 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
1311 .S_L;
1312 if (this->
process_data_.explicit_hm_coupling_in_unsaturated_zone)
1313 {
1314 double const chi_S_L_prev = std::get<PrevState<
1316 ->chi_S_L;
1317 Kpu.noalias() += N_p.transpose() * chi_S_L_prev * rho_LR *
alpha *
1318 identity2.transpose() * B * w;
1319 }
1320 else
1321 {
1322 Kpu.noalias() += N_p.transpose() * S_L * rho_LR *
alpha *
1323 identity2.transpose() * B * w;
1324 }
1325
1326
1327
1328
1329
1330 double const k_rel =
1332 DisplacementDim>>(CD)
1333 .k_rel;
1334 auto const& K_intrinsic =
1336 DisplacementDim>>(CD)
1337 .Ki;
1338 double const mu =
1339 *std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(
1340 CD);
1341
1343
1344 laplace_p.noalias() +=
1345 dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
1346
1347 auto const beta_LR =
1348 1 / rho_LR *
1349 liquid_phase.property(MPL::PropertyType::density)
1350 .template dValue<double>(variables,
1351 MPL::Variable::liquid_phase_pressure,
1352 x_position, t, dt);
1353
1354 double const beta_SR =
1355 std::get<
1357 CD)
1358 .beta_SR;
1359 double const a0 = (
alpha - phi) * beta_SR;
1360 double const specific_storage_a_p = S_L * (phi * beta_LR + S_L * a0);
1361 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
1362
1363 double const dspecific_storage_a_p_dp_cap =
1364 dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0);
1365 double const dspecific_storage_a_S_dp_cap =
1366 -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
1367
1368 storage_p_a_p.noalias() +=
1369 N_p.transpose() * rho_LR * specific_storage_a_p * N_p * w;
1370
1371 double const DeltaS_L_Deltap_cap =
1372 std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap;
1373 storage_p_a_S.noalias() -= N_p.transpose() * rho_LR *
1374 specific_storage_a_S * DeltaS_L_Deltap_cap *
1375 N_p * w;
1376
1377 local_Jac
1380 .noalias() += N_p.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
1381 rho_LR * dspecific_storage_a_p_dp_cap * N_p * w;
1382
1383 double const S_L_prev =
1384 std::get<
1385 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
1387 ->S_L;
1388 storage_p_a_S_Jpp.noalias() -=
1389 N_p.transpose() * rho_LR *
1390 ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
1391 specific_storage_a_S * dS_L_dp_cap) /
1392 dt * N_p * w;
1393
1394 if (!this->
process_data_.explicit_hm_coupling_in_unsaturated_zone)
1395 {
1396 local_Jac
1399 .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap *
alpha *
1400 identity2.transpose() * B * (u - u_prev) / dt *
1401 N_p * w;
1402 }
1403
1404 double const dk_rel_dS_l =
1405 medium->property(MPL::PropertyType::relative_permeability)
1406 .template dValue<double>(variables,
1407 MPL::Variable::liquid_saturation,
1408 x_position, t, dt);
1410 grad_p_cap = -dNdx_p * p_L;
1411 local_Jac
1414 .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
1415 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1416
1417 local_Jac
1420 .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
1421 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1422
1423 local_rhs.template segment<pressure_size>(
pressure_index).noalias() +=
1424 dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
1425
1426 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
1427 {
1428 double const alpha_bar =
1430 ->mass_exchange_coefficient;
1431 auto const p_L_m =
1434 .noalias() -=
1435 N_p.transpose() * alpha_bar / mu * (-p_cap_ip - p_L_m) * w;
1436
1437 local_Jac
1440 .noalias() += N_p.transpose() * alpha_bar / mu * N_p * w;
1441 if (p_cap_ip != p_cap_prev_ip)
1442 {
1443 auto const p_L_m_prev = **std::get<PrevState<MicroPressure>>(
1445 local_Jac
1446 .template block<pressure_size, pressure_size>(
1448 .noalias() += N_p.transpose() * alpha_bar / mu *
1449 (p_L_m - p_L_m_prev) /
1450 (p_cap_ip - p_cap_prev_ip) * N_p * w;
1451 }
1452 }
1453 }
1454
1456 {
1457 storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
1458 storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
1459 storage_p_a_S_Jpp =
1460 storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
1461 }
1462
1463
1464 local_Jac
1467 .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
1468
1469
1470 local_Jac
1471 .template block<pressure_size, displacement_size>(
pressure_index,
1473 .noalias() = Kpu / dt;
1474
1475
1476 local_rhs.template segment<pressure_size>(
pressure_index).noalias() -=
1477 laplace_p * p_L +
1478 (storage_p_a_p + storage_p_a_S) * (p_L - p_L_prev) / dt +
1479 Kpu * (u - u_prev) / dt;
1480
1481
1483 .noalias() += Kup * p_L;
1484}
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.12.0