1073{
1075
1076 auto const [p_L, u] =
localDOF(local_x);
1077 auto const [p_L_prev, u_prev] =
localDOF(local_x_prev);
1078
1080 typename ShapeMatricesTypeDisplacement::template MatrixType<
1085
1087 typename ShapeMatricesTypeDisplacement::template VectorType<
1090
1093 DisplacementDim)>::identity2;
1094
1096 ShapeMatricesTypePressure::NodalMatrixType::Zero(
pressure_size,
1098
1100 ShapeMatricesTypePressure::NodalMatrixType::Zero(
pressure_size,
1102
1104 ShapeMatricesTypePressure::NodalMatrixType::Zero(
pressure_size,
1106
1108 ShapeMatricesTypePressure::NodalMatrixType::Zero(
pressure_size,
1110
1111 typename ShapeMatricesTypeDisplacement::template MatrixType<
1113 Kup = ShapeMatricesTypeDisplacement::template MatrixType<
1116
1117 typename ShapeMatricesTypeDisplacement::template MatrixType<
1119 Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
1122
1123 auto const& medium =
1125 auto const& liquid_phase = medium->phase("AqueousLiquid");
1126 auto const& solid_phase = medium->phase("Solid");
1129
1130 unsigned const n_integration_points =
1132 for (unsigned ip = 0; ip < n_integration_points; ip++)
1133 {
1139
1140 auto const& w =
ip_data_[ip].integration_weight;
1141
1142 auto const& N_u =
ip_data_[ip].N_u;
1143 auto const& dNdx_u =
ip_data_[ip].dNdx_u;
1144
1145 auto const& N_p =
ip_data_[ip].N_p;
1146 auto const& dNdx_p =
ip_data_[ip].dNdx_p;
1147
1155
1156 auto const B =
1158 ShapeFunctionDisplacement::NPOINTS,
1161
1162 double p_cap_ip;
1164
1165 double p_cap_prev_ip;
1167
1170
1171
1173
1175 medium->property(MPL::PropertyType::reference_temperature)
1176 .template value<double>(variables, x_position, t, dt);
1178
1179 std::get<StrainData<DisplacementDim>>(SD).eps.noalias() = B * u;
1180
1182 t, dt, x_position,
ip_data_[ip], variables, variables_prev, medium,
1184 CapillaryPressureData<DisplacementDim>{
1185 p_cap_ip, p_cap_prev_ip,
1186 Eigen::Vector<double, DisplacementDim>::Zero()},
1187 CD, SD, SD_prev, this->
process_data_.micro_porosity_parameters,
1189
1190 {
1191 auto const& C = *std::get<StiffnessTensor<DisplacementDim>>(CD);
1192 local_Jac
1193 .template block<displacement_size, displacement_size>(
1195 .noalias() += B.transpose() * C * B * w;
1196 }
1197
1199
1200 {
1201 auto const& sigma_eff =
1204 .sigma_eff;
1205 double const rho = *std::get<Density>(CD);
1207 .noalias() -= (B.transpose() * sigma_eff -
1208 N_u_op(N_u).transpose() * rho * b) *
1209 w;
1210 }
1211
1212
1213
1214
1215
1216 double const alpha =
1217 *std::get<ProcessLib::ThermoRichardsMechanics::BiotData>(CD);
1218 double const dS_L_dp_cap =
1219 std::get<ProcessLib::ThermoRichardsMechanics::SaturationDataDeriv>(
1220 CD)
1221 .dS_L_dp_cap;
1222
1223 {
1224 double const chi_S_L =
1225 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
1226 .chi_S_L;
1227 Kup.noalias() +=
1228 B.transpose() *
alpha * chi_S_L * identity2 * N_p * w;
1229 double const dchi_dS_L =
1230 std::get<ProcessLib::ThermoRichardsMechanics::BishopsData>(CD)
1231 .dchi_dS_L;
1232
1233 local_Jac
1234 .template block<displacement_size, pressure_size>(
1236 .noalias() -= B.transpose() *
alpha *
1237 (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
1238 identity2 * N_p * w;
1239 }
1240
1241 double const phi =
1242 std::get<ProcessLib::ThermoRichardsMechanics::PorosityData>(CD).phi;
1243 double const rho_LR = *std::get<LiquidDensity>(CD);
1244 local_Jac
1245 .template block<displacement_size, pressure_size>(
1247 .noalias() +=
1248 N_u_op(N_u).transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259 if (!medium->hasProperty(MPL::PropertyType::saturation_micro) &&
1260 solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
1261 {
1262 using DimMatrix = Eigen::Matrix<double, 3, 3>;
1263 auto const dsigma_sw_dS_L =
1265 solid_phase
1266 .property(MPL::PropertyType::swelling_stress_rate)
1267 .template dValue<DimMatrix>(
1268 variables, variables_prev,
1269 MPL::Variable::liquid_saturation, x_position, t,
1270 dt));
1271 local_Jac
1272 .template block<displacement_size, pressure_size>(
1274 .noalias() +=
1275 B.transpose() * dsigma_sw_dS_L * dS_L_dp_cap * N_p * w;
1276 }
1277
1278
1279
1280 double const S_L =
1281 std::get<ProcessLib::ThermoRichardsMechanics::SaturationData>(
1283 .S_L;
1284 if (this->
process_data_.explicit_hm_coupling_in_unsaturated_zone)
1285 {
1286 double const chi_S_L_prev = std::get<PrevState<
1288 ->chi_S_L;
1289 Kpu.noalias() += N_p.transpose() * chi_S_L_prev * rho_LR *
alpha *
1290 identity2.transpose() * B * w;
1291 }
1292 else
1293 {
1294 Kpu.noalias() += N_p.transpose() * S_L * rho_LR *
alpha *
1295 identity2.transpose() * B * w;
1296 }
1297
1298
1299
1300
1301
1302 double const k_rel =
1304 DisplacementDim>>(CD)
1305 .k_rel;
1306 auto const& K_intrinsic =
1308 DisplacementDim>>(CD)
1309 .Ki;
1310 double const mu =
1311 *std::get<ProcessLib::ThermoRichardsMechanics::LiquidViscosityData>(
1312 CD);
1313
1315
1316 laplace_p.noalias() +=
1317 dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
1318
1319 auto const beta_LR =
1320 1 / rho_LR *
1321 liquid_phase.property(MPL::PropertyType::density)
1322 .template dValue<double>(variables,
1323 MPL::Variable::liquid_phase_pressure,
1324 x_position, t, dt);
1325
1326 double const beta_SR =
1327 std::get<
1329 CD)
1330 .beta_SR;
1331 double const a0 = (
alpha - phi) * beta_SR;
1332 double const specific_storage_a_p = S_L * (phi * beta_LR + S_L * a0);
1333 double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
1334
1335 double const dspecific_storage_a_p_dp_cap =
1336 dS_L_dp_cap * (phi * beta_LR + 2 * S_L * a0);
1337 double const dspecific_storage_a_S_dp_cap =
1338 -a0 * (S_L + p_cap_ip * dS_L_dp_cap);
1339
1340 storage_p_a_p.noalias() +=
1341 N_p.transpose() * rho_LR * specific_storage_a_p * N_p * w;
1342
1343 double const DeltaS_L_Deltap_cap =
1344 std::get<SaturationSecantDerivative>(CD).DeltaS_L_Deltap_cap;
1345 storage_p_a_S.noalias() -= N_p.transpose() * rho_LR *
1346 specific_storage_a_S * DeltaS_L_Deltap_cap *
1347 N_p * w;
1348
1349 local_Jac
1352 .noalias() += N_p.transpose() * (p_cap_ip - p_cap_prev_ip) / dt *
1353 rho_LR * dspecific_storage_a_p_dp_cap * N_p * w;
1354
1355 double const S_L_prev =
1356 std::get<
1357 PrevState<ProcessLib::ThermoRichardsMechanics::SaturationData>>(
1359 ->S_L;
1360 storage_p_a_S_Jpp.noalias() -=
1361 N_p.transpose() * rho_LR *
1362 ((S_L - S_L_prev) * dspecific_storage_a_S_dp_cap +
1363 specific_storage_a_S * dS_L_dp_cap) /
1364 dt * N_p * w;
1365
1366 if (!this->
process_data_.explicit_hm_coupling_in_unsaturated_zone)
1367 {
1368 local_Jac
1371 .noalias() -= N_p.transpose() * rho_LR * dS_L_dp_cap *
alpha *
1372 identity2.transpose() * B * (u - u_prev) / dt *
1373 N_p * w;
1374 }
1375
1376 double const dk_rel_dS_l =
1377 medium->property(MPL::PropertyType::relative_permeability)
1378 .template dValue<double>(variables,
1379 MPL::Variable::liquid_saturation,
1380 x_position, t, dt);
1382 grad_p_cap = -dNdx_p * p_L;
1383 local_Jac
1386 .noalias() += dNdx_p.transpose() * rho_Ki_over_mu * grad_p_cap *
1387 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1388
1389 local_Jac
1392 .noalias() += dNdx_p.transpose() * rho_LR * rho_Ki_over_mu * b *
1393 dk_rel_dS_l * dS_L_dp_cap * N_p * w;
1394
1395 local_rhs.template segment<pressure_size>(
pressure_index).noalias() +=
1396 dNdx_p.transpose() * rho_LR * k_rel * rho_Ki_over_mu * b * w;
1397
1398 if (medium->hasProperty(MPL::PropertyType::saturation_micro))
1399 {
1400 double const alpha_bar =
1402 ->mass_exchange_coefficient;
1403 auto const p_L_m =
1406 .noalias() -=
1407 N_p.transpose() * alpha_bar / mu * (-p_cap_ip - p_L_m) * w;
1408
1409 local_Jac
1412 .noalias() += N_p.transpose() * alpha_bar / mu * N_p * w;
1413 if (p_cap_ip != p_cap_prev_ip)
1414 {
1415 auto const p_L_m_prev = **std::get<PrevState<MicroPressure>>(
1417 local_Jac
1418 .template block<pressure_size, pressure_size>(
1420 .noalias() += N_p.transpose() * alpha_bar / mu *
1421 (p_L_m - p_L_m_prev) /
1422 (p_cap_ip - p_cap_prev_ip) * N_p * w;
1423 }
1424 }
1425 }
1426
1428 {
1429 storage_p_a_p = storage_p_a_p.colwise().sum().eval().asDiagonal();
1430 storage_p_a_S = storage_p_a_S.colwise().sum().eval().asDiagonal();
1431 storage_p_a_S_Jpp =
1432 storage_p_a_S_Jpp.colwise().sum().eval().asDiagonal();
1433 }
1434
1435
1436 local_Jac
1439 .noalias() += laplace_p + storage_p_a_p / dt + storage_p_a_S_Jpp;
1440
1441
1442 local_Jac
1443 .template block<pressure_size, displacement_size>(
pressure_index,
1445 .noalias() = Kpu / dt;
1446
1447
1448 local_rhs.template segment<pressure_size>(
pressure_index).noalias() -=
1449 laplace_p * p_L +
1450 (storage_p_a_p + storage_p_a_S) * (p_L - p_L_prev) / dt +
1451 Kpu * (u - u_prev) / dt;
1452
1453
1455 .noalias() += Kup * p_L;
1456}
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