Commit 79b1dcd0 authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov
Browse files

[MPL] Replace effective_pore_pressure_rate.

Using increment version.
parent 46f9c898
......@@ -60,13 +60,17 @@ PropertyDataType PorosityFromMassBalance::value(
double const e_dot = std::get<double>(
variable_array[static_cast<int>(Variable::volumetric_strain_rate)]);
double const p_eff_dot = std::get<double>(variable_array[static_cast<int>(
Variable::effective_pore_pressure_rate)]);
double const p_eff = std::get<double>(
variable_array[static_cast<int>(Variable::effective_pore_pressure)]);
double const p_eff_prev =
std::get<double>(variable_array_prev[static_cast<int>(
Variable::effective_pore_pressure)]);
double const delta_p_eff = p_eff - p_eff_prev;
double const phi_prev = std::get<double>(
variable_array_prev[static_cast<int>(Variable::porosity)]);
double const w = dt * (e_dot + p_eff_dot * beta_SR);
double const w = dt * e_dot + delta_p_eff * beta_SR;
return std::clamp((phi_prev + alpha_b * w) / (1 + w), phi_min_, phi_max_);
}
......
......@@ -50,8 +50,12 @@ PropertyDataType TransportPorosityFromMassBalance::value(
double const e_dot = std::get<double>(
variable_array[static_cast<int>(Variable::volumetric_strain_rate)]);
double const p_eff_dot = std::get<double>(variable_array[static_cast<int>(
Variable::effective_pore_pressure_rate)]);
double const p_eff = std::get<double>(
variable_array[static_cast<int>(Variable::effective_pore_pressure)]);
double const p_eff_prev =
std::get<double>(variable_array_prev[static_cast<int>(
Variable::effective_pore_pressure)]);
double const delta_p_eff = p_eff - p_eff_prev;
double const phi =
std::get<double>(variable_array[static_cast<int>(Variable::porosity)]);
......@@ -59,7 +63,7 @@ PropertyDataType TransportPorosityFromMassBalance::value(
double const phi_tr_prev = std::get<double>(
variable_array_prev[static_cast<int>(Variable::transport_porosity)]);
double const w = dt * (e_dot + p_eff_dot * beta_SR);
double const w = dt * e_dot + delta_p_eff * beta_SR;
return std::clamp(phi_tr_prev + (alpha_b - phi) * w, phi_min_, phi_max_);
}
......
......@@ -46,7 +46,7 @@ enum class Variable : int
concentration,
density,
displacement,
effective_pore_pressure_rate,
effective_pore_pressure,
grain_compressibility,
liquid_saturation,
phase_pressure,
......
......@@ -380,11 +380,11 @@ void RichardsMechanicsLocalAssembler<
double const chi_S_L = chi(S_L);
double const chi_S_L_prev = chi(S_L_prev);
variables[static_cast<int>(
MPL::Variable::effective_pore_pressure_rate)] =
(chi_S_L * (-p_cap_ip) -
chi_S_L_prev * (-p_cap_ip + p_cap_dot_ip * dt)) /
dt;
variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
-chi_S_L * p_cap_ip;
variables_prev[static_cast<int>(
MPL::Variable::effective_pore_pressure)] =
-chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
// Set volumetric strain rate for the general case without swelling.
double const div_u_dot = identity2.transpose() * B * u_dot;
......@@ -708,11 +708,11 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
double const chi_S_L = chi(S_L);
double const chi_S_L_prev = chi(S_L_prev);
variables[static_cast<int>(
MPL::Variable::effective_pore_pressure_rate)] =
(chi_S_L * (-p_cap_ip) -
chi_S_L_prev * (-p_cap_ip + p_cap_dot_ip * dt)) /
dt;
variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
-chi_S_L * p_cap_ip;
variables_prev[static_cast<int>(
MPL::Variable::effective_pore_pressure)] =
-chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
// Set volumetric strain rate for the general case without swelling.
double const div_u_dot = identity2.transpose() * B * u_dot;
......@@ -1434,11 +1434,11 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables[static_cast<int>(MPL::Variable::grain_compressibility)] =
beta_SR;
variables[static_cast<int>(
MPL::Variable::effective_pore_pressure_rate)] =
(chi_S_L * (-p_cap_ip) -
chi_S_L_prev * (-p_cap_ip + p_cap_dot_ip * dt)) /
dt;
variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
-chi_S_L * p_cap_ip;
variables_prev[static_cast<int>(
MPL::Variable::effective_pore_pressure)] =
-chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
// Set volumetric strain rate for the general case without swelling.
double const div_u_dot = identity2.transpose() * B * u_dot;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment