Commit 52374cbf authored by wenqing's avatar wenqing

[TRM] Implemented water vaporization terms

parent e31313ca
......@@ -17,11 +17,13 @@
#include "MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.h"
#include "MaterialLib/MPL/Utils/FormEigenTensor.h"
#include "MaterialLib/MPL/Utils/GetLiquidThermalExpansivity.h"
#include "MaterialLib/PhysicalConstant.h"
#include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
#include "MathLib/KelvinVector.h"
#include "NumLib/Function/Interpolation.h"
#include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
#include "ThermoRichardsMechanicsFEM.h"
#include "WaterVaporProperty.h"
namespace ProcessLib
{
......@@ -756,6 +758,48 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
N_p.transpose() * (dNdx_p * T).transpose() *
Ki_over_mu * dNdx_p * w;
}
if (process_data_.has_water_vaporization)
{
double const p_L_ip = -p_cap_ip;
double const rho_wv = waterVaporDensity(T_ip, p_L_ip, rho_LR);
double const storage_coefficient_by_water_vapor =
phi * rho_wv *
(dS_L_dp_cap +
(1 - S_L) / (rho_LR * T_ip *
MaterialLib::PhysicalConstant::
SpecificGasConstant::WaterVapour));
storage_p_a_p.noalias() +=
N_p.transpose() * storage_coefficient_by_water_vapor * N_p * w;
double const vapor_expansion =
phi * (1 - S_L) * dwaterVaporDensitydT(T_ip, p_L_ip, rho_LR);
M_pT.noalias() += N_p.transpose() * vapor_expansion * N_p * w;
auto const f_Tv =
liquid_phase
.property(MaterialPropertyLib::PropertyType::
thermal_diffusion_enhancement_factor)
.template value<double>(variables, x_position, t, dt);
auto const tortuosity =
medium->property(MaterialPropertyLib::PropertyType::tortuosity)
.template value<double>(variables, x_position, t, dt);
double const f_Tv_D_Tv =
f_Tv * DTv(T_ip, p_L_ip, rho_LR, S_L, phi, tortuosity);
local_Jac
.template block<pressure_size, temperature_size>(
pressure_index, temperature_index)
.noalias() += dNdx_p.transpose() * f_Tv_D_Tv * dNdx_p * w;
local_rhs.template segment<pressure_size>(pressure_index)
.noalias() -= f_Tv_D_Tv * dNdx_p.transpose() * (dNdx_p * T) * w;
laplace_p.noalias() +=
dNdx_p.transpose() *
Dpv(T_ip, p_L_ip, rho_LR, S_L, phi, tortuosity) * dNdx_p * w;
}
}
if (process_data_.apply_mass_lumping)
......
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