[TRM] Implemented water vaporization terms

ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h

@@ -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.transpose() * (dNdx * T).transpose() *
                               Ki_over_mu * dNdx * 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.transpose() * storage_coefficient_by_water_vapor * N * w;
+
+            double const vapor_expansion =
+                phi * (1 - S_L) * dwaterVaporDensitydT(T_ip, p_L_ip, rho_LR);
+            M_pT.noalias() += N.transpose() * vapor_expansion * N * 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.transpose() * f_Tv_D_Tv * dNdx * w;
+
+            local_rhs.template segment<pressure_size>(pressure_index)
+                .noalias() -= f_Tv_D_Tv * dNdx.transpose() * (dNdx * T) * w;
+
+            laplace_p.noalias() +=
+                dNdx.transpose() *
+                Dpv(T_ip, p_L_ip, rho_LR, S_L, phi, tortuosity) * dNdx * w;
+        }
     }
 
     if (process_data_.apply_mass_lumping)