Commit 03b07337 authored by wenqing's avatar wenqing

[RM] Improved the pressure calculation for the total stress for the RHS

parent 1571bbfc
......@@ -615,12 +615,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
ShapeMatricesTypePressure::NodalMatrixType::Zero(pressure_size,
pressure_size);
typename ShapeMatricesTypeDisplacement::template MatrixType<
displacement_size, pressure_size>
Kup = ShapeMatricesTypeDisplacement::template MatrixType<
displacement_size, pressure_size>::Zero(displacement_size,
pressure_size);
typename ShapeMatricesTypeDisplacement::template MatrixType<
pressure_size, displacement_size>
Kpu = ShapeMatricesTypeDisplacement::template MatrixType<
......@@ -844,15 +838,17 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
.template value<double>(variables, x_position, t, dt);
double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
typename BMatricesType::KelvinVectorType total_stress =
sigma_eff + sigma_sw - alpha * chi_S_L * (-p_cap_ip) * identity2;
local_rhs.template segment<displacement_size>(displacement_index)
.noalias() -= (B.transpose() * (sigma_eff + sigma_sw) -
N_u_op.transpose() * rho * b) *
w;
.noalias() -=
(B.transpose() * total_stress - N_u_op.transpose() * rho * b) * w;
//
// displacement equation, pressure part
//
Kup.noalias() += B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
auto const dchi_dS_L =
medium->property(MPL::PropertyType::bishops_effective_stress)
......@@ -984,10 +980,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
local_rhs.template segment<pressure_size>(pressure_index).noalias() -=
laplace_p * p_L + (storage_p_a_p + storage_p_a_S) * p_L_dot +
Kpu * u_dot;
// displacement equation
local_rhs.template segment<displacement_size>(displacement_index)
.noalias() += Kup * p_L;
}
template <int Components, typename StoreValuesFunction>
......
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