[PL/TRM/TRF] write vapor advection in terms of flux

ufz · Nov 28, 2023 · e1f8858 · e1f8858
1 parent 05a2995
commit e1f8858
Show file tree

Hide file tree

Showing 4 changed files with 23 additions and 27 deletions.
diff --git a/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h
@@ -541,8 +541,8 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
             double const D_pv = D_v * drho_wv_dp;
 
             GlobalDimVectorType const grad_T = dNdx * T;
-            GlobalDimVectorType const vapour_velocity =
-                -(f_Tv_D_Tv * grad_T - D_pv * grad_p_cap) / rho_LR;
+            GlobalDimVectorType const vapour_flux =
+                -(f_Tv_D_Tv * grad_T - D_pv * grad_p_cap);
             double const specific_heat_capacity_vapour =
                 gas_phase->property(MaterialPropertyLib::specific_heat_capacity)
                     .template value<double>(variables, x_position, t, dt);
@@ -551,9 +551,8 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
                 w * (rho_wv * specific_heat_capacity_vapour * (1 - S_L) * phi) *
                 N.transpose() * N;
 
-            K_TT.noalias() += N.transpose() * vapour_velocity.transpose() *
-                                  dNdx * rho_wv *
-                                  specific_heat_capacity_vapour * w;
+            K_TT.noalias() += N.transpose() * vapour_flux.transpose() * dNdx *
+                              specific_heat_capacity_vapour * w;
 
             double const storage_coefficient_by_water_vapor =
                 phi * (rho_wv * dS_L_dp_cap + (1 - S_L) * drho_wv_dp);
@@ -577,8 +576,7 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
             //
             // Latent heat term
             //
-            if (gas_phase->hasProperty(
-                    MPL::PropertyType::specific_latent_heat))
+            if (gas_phase->hasProperty(MPL::PropertyType::specific_latent_heat))
             {
                 double const factor = phi * (1 - S_L) / rho_LR;
                 // The volumetric latent heat of vaporization of liquid water
@@ -1000,11 +998,11 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::assemble(
 
             GlobalDimVectorType const grad_T = dNdx * T;
             GlobalDimVectorType const grad_p_cap = -dNdx * p_L;
-            GlobalDimVectorType const vapour_velocty =
-                    -(f_Tv_D_Tv * grad_T - D_pv * grad_p_cap) / rho_LR;
+            GlobalDimVectorType const vapour_flux =
+                -(f_Tv_D_Tv * grad_T - D_pv * grad_p_cap);
             double const specific_heat_capacity_vapour =
                 gas_phase->property(MaterialPropertyLib::specific_heat_capacity)
-                        .template value<double>(variables, x_position, t, dt);
+                    .template value<double>(variables, x_position, t, dt);
 
             local_M
                 .template block<temperature_size, temperature_size>(
@@ -1016,8 +1014,8 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::assemble(
             local_K
                 .template block<temperature_size, temperature_size>(
                     temperature_index, temperature_index)
-                .noalias() += N.transpose() * vapour_velocty.transpose() *
-                                dNdx * rho_wv * specific_heat_capacity_vapour * w;
+                .noalias() += N.transpose() * vapour_flux.transpose() * dNdx *
+                              specific_heat_capacity_vapour * w;
 
             double const storage_coefficient_by_water_vapor =
                 phi * (rho_wv * dS_L_dp_cap + (1 - S_L) * drho_wv_dp);

diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EqT.cpp b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/EqT.cpp
@@ -20,9 +20,8 @@ void EqTModel<DisplacementDim>::eval(
 {
     out.M_TT_X_NTN = heat_data.M_TT_X_NTN + vap_data.M_TT_X_NTN;
 
-    out.K_TT_NT_V_dN =
-        heat_data.advective_heat_flux_contribution_to_K_liquid +
-        vap_data.vapor_velocity * vap_data.volumetric_heat_capacity_vapor;
+    out.K_TT_NT_V_dN = heat_data.advective_heat_flux_contribution_to_K_liquid +
+                       vap_data.vapor_flux * vap_data.heat_capacity_vapor;
 }
 
 template struct EqTModel<2>;

diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMVaporDiffusion.cpp b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMVaporDiffusion.cpp
@@ -15,8 +15,8 @@ namespace ProcessLib::ThermoRichardsMechanics
 template <int DisplacementDim>
 void TRMVaporDiffusionData<DisplacementDim>::setZero()
 {
-    volumetric_heat_capacity_vapor = 0;
-    vapor_velocity = GlobalDimVector<DisplacementDim>::Zero(DisplacementDim);
+    heat_capacity_vapor = 0;
+    vapor_flux = GlobalDimVector<DisplacementDim>::Zero(DisplacementDim);
     storage_coefficient_by_water_vapor = 0;
 
     J_pT_X_dNTdN = 0;
@@ -81,17 +81,16 @@ void TRMVaporDiffusionModel<DisplacementDim>::eval(
         out.J_pT_X_dNTdN = f_Tv * D_v * drho_wv_dT;
         out.K_pp_X_dNTdN = D_v * drho_wv_dp;
 
-        out.vapor_velocity = -(out.J_pT_X_dNTdN * T_data.grad_T -
-                      out.K_pp_X_dNTdN * p_cap_data.grad_p_cap) /
-                    rho_L_data.rho_LR;
-        double const specific_heat_capacity_vapor =
-            gas_phase.property(MaterialPropertyLib::specific_heat_capacity)
+        out.vapor_flux = -(out.J_pT_X_dNTdN * T_data.grad_T -
+                           out.K_pp_X_dNTdN * p_cap_data.grad_p_cap);
+        out.heat_capacity_vapor =
+            gas_phase
+                ->property(
+                    MaterialPropertyLib::PropertyType::specific_heat_capacity)
                 .template value<double>(variables, x_t.x, x_t.t, x_t.dt);
 
-        out.volumetric_heat_capacity_vapor =
-            rho_wv * specific_heat_capacity_vapor;
         out.M_TT_X_NTN +=
-            out.volumetric_heat_capacity_vapor * (1 - S_L_data.S_L) * phi;
+            out.heat_capacity_vapor * rho_wv * (1 - S_L_data.S_L) * phi;
 
         out.storage_coefficient_by_water_vapor =
             phi *

diff --git a/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMVaporDiffusion.h b/ProcessLib/ThermoRichardsMechanics/ConstitutiveCommon/TRMVaporDiffusion.h
@@ -18,8 +18,8 @@ namespace ProcessLib::ThermoRichardsMechanics
 template <int DisplacementDim>
 struct TRMVaporDiffusionData
 {
-    double volumetric_heat_capacity_vapor;
-    GlobalDimVector<DisplacementDim> vapor_velocity;
+    double heat_capacity_vapor;
+    GlobalDimVector<DisplacementDim> vapor_flux;
     double storage_coefficient_by_water_vapor;
 
     double J_pT_X_dNTdN;