Commit b87ea9d8 authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov Committed by Dmitry Yu. Naumov
Browse files

[PL/RM] Use phi for local name of porosity.

parent 3c9c07b8
......@@ -391,15 +391,15 @@ void RichardsMechanicsLocalAssembler<
variables[static_cast<int>(MPL::Variable::volumetric_strain_rate)]
.emplace<double>(div_u_dot);
auto& porosity = _ip_data[ip].porosity;
auto& phi = _ip_data[ip].porosity;
{ // Porosity update
variables_prev[static_cast<int>(MPL::Variable::porosity)] =
_ip_data[ip].porosity_prev;
porosity = solid_phase.property(MPL::PropertyType::porosity)
.template value<double>(variables, variables_prev,
x_position, t, dt);
phi = solid_phase.property(MPL::PropertyType::porosity)
.template value<double>(variables, variables_prev,
x_position, t, dt);
}
variables[static_cast<int>(MPL::Variable::porosity)] = porosity;
variables[static_cast<int>(MPL::Variable::porosity)] = phi;
// Swelling and possibly volumetric strain rate update.
auto& sigma_sw = _ip_data[ip].sigma_sw;
......@@ -449,7 +449,7 @@ void RichardsMechanicsLocalAssembler<
else
{
variables[static_cast<int>(MPL::Variable::transport_porosity)] =
porosity;
phi;
}
}
......@@ -490,18 +490,17 @@ void RichardsMechanicsLocalAssembler<
displacement_index, displacement_index)
.noalias() += B.transpose() * C * B * w;
double const rho = rho_SR * (1 - porosity) + S_L * porosity * rho_LR;
double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
rhs.template segment<displacement_size>(displacement_index).noalias() +=
N_u_op.transpose() * rho * b * w - B.transpose() * sigma_sw * w;
//
// pressure equation, pressure part.
//
double const a0 = S_L * (alpha - porosity) * beta_SR;
double const a0 = S_L * (alpha - phi) * beta_SR;
// Volumetric average specific storage of the solid and fluid phases.
double const specific_storage =
dS_L_dp_cap * (p_cap_ip * a0 - porosity) +
S_L * (porosity / K_LR + a0);
dS_L_dp_cap * (p_cap_ip * a0 - phi) + S_L * (phi / K_LR + a0);
M.template block<pressure_size, pressure_size>(pressure_index,
pressure_index)
.noalias() += N_p.transpose() * rho_LR * specific_storage * N_p * w;
......@@ -722,15 +721,15 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables[static_cast<int>(MPL::Variable::volumetric_strain_rate)]
.emplace<double>(div_u_dot);
auto& porosity = _ip_data[ip].porosity;
auto& phi = _ip_data[ip].porosity;
{ // Porosity update
variables_prev[static_cast<int>(MPL::Variable::porosity)] =
_ip_data[ip].porosity_prev;
porosity = solid_phase.property(MPL::PropertyType::porosity)
.template value<double>(variables, variables_prev,
x_position, t, dt);
variables[static_cast<int>(MPL::Variable::porosity)] = porosity;
phi = solid_phase.property(MPL::PropertyType::porosity)
.template value<double>(variables, variables_prev,
x_position, t, dt);
variables[static_cast<int>(MPL::Variable::porosity)] = phi;
}
// Swelling and possibly volumetric strain rate update.
......@@ -778,7 +777,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
else
{
variables[static_cast<int>(MPL::Variable::transport_porosity)] =
porosity;
phi;
}
}
......@@ -815,7 +814,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
displacement_index, displacement_index)
.noalias() += B.transpose() * C * B * w;
auto const phi = _ip_data[ip].porosity;
double const p_FR = -chi_S_L * p_cap_ip;
// p_SR
variables[static_cast<int>(MPL::Variable::solid_grain_pressure)] =
......@@ -824,7 +822,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
solid_phase.property(MPL::PropertyType::density)
.template value<double>(variables, x_position, t, dt);
double const rho = rho_SR * (1 - porosity) + S_L * porosity * rho_LR;
double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
local_rhs.template segment<displacement_size>(displacement_index)
.noalias() -= (B.transpose() * (sigma_eff + sigma_sw) -
N_u_op.transpose() * rho * b) *
......@@ -851,7 +849,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
.template block<displacement_size, pressure_size>(
displacement_index, pressure_index)
.noalias() +=
N_u_op.transpose() * porosity * rho_LR * dS_L_dp_cap * b * N_p * w;
N_u_op.transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
{
......@@ -881,12 +879,12 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
laplace_p.noalias() +=
dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
double const a0 = (alpha - porosity) * beta_SR;
double const specific_storage_a_p = S_L * (porosity / K_LR + S_L * a0);
double const specific_storage_a_S = porosity - p_cap_ip * S_L * a0;
double const a0 = (alpha - phi) * beta_SR;
double const specific_storage_a_p = S_L * (phi / K_LR + S_L * a0);
double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
double const dspecific_storage_a_p_dp_cap =
dS_L_dp_cap * (porosity / K_LR + 2 * S_L * a0);
dS_L_dp_cap * (phi / K_LR + 2 * S_L * a0);
double const dspecific_storage_a_S_dp_cap =
-a0 * (S_L + p_cap_ip * dS_L_dp_cap);
......@@ -1451,14 +1449,14 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables[static_cast<int>(MPL::Variable::volumetric_strain_rate)]
.emplace<double>(div_u_dot);
auto& porosity = _ip_data[ip].porosity;
auto& phi = _ip_data[ip].porosity;
{ // Porosity update
variables_prev[static_cast<int>(MPL::Variable::porosity)] =
_ip_data[ip].porosity_prev;
porosity = solid_phase.property(MPL::PropertyType::porosity)
.template value<double>(variables, variables_prev,
x_position, t, dt);
variables[static_cast<int>(MPL::Variable::porosity)] = porosity;
phi = solid_phase.property(MPL::PropertyType::porosity)
.template value<double>(variables, variables_prev,
x_position, t, dt);
variables[static_cast<int>(MPL::Variable::porosity)] = phi;
}
// Swelling and possibly volumetric strain rate update.
......@@ -1505,7 +1503,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
else
{
variables[static_cast<int>(MPL::Variable::transport_porosity)] =
porosity;
phi;
}
}
auto const mu =
......@@ -1532,7 +1530,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
auto const phi = _ip_data[ip].porosity;
auto const& sigma_eff = _ip_data[ip].sigma_eff;
double const p_FR = -chi_S_L * p_cap_ip;
// p_SR
......@@ -1556,9 +1553,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
-K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b;
saturation_avg += S_L;
porosity_avg +=
_ip_data[ip].porosity; // Note, this is not updated, because needs
// xdot and dt to be passed.
porosity_avg += phi;
sigma_avg += sigma_eff;
}
saturation_avg /= n_integration_points;
......
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