Merge branch 'fixdensities' into 'master'

Fix latent heat K_Tp terms in TRF/TRM See merge request ogs/ogs!4455
ufz · Nov 28, 2023 · c147d93 · c147d93
2 parents 7ddd91c + e1f8858
commit c147d93
Show file tree

Hide file tree

Showing 13 changed files with 105 additions and 198 deletions.
diff --git a/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h b/ProcessLib/ThermoRichardsFlow/ThermoRichardsFlowFEM-impl.h
@@ -541,9 +541,8 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
             double const D_pv = D_v * drho_wv_dp;
 
             GlobalDimVectorType const grad_T = dNdx * T;
-            // Vapour velocity
-            GlobalDimVectorType const q_v =
-                -(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);
@@ -552,7 +551,7 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
                 w * (rho_wv * specific_heat_capacity_vapour * (1 - S_L) * phi) *
                 N.transpose() * N;
 
-            K_TT.noalias() += N.transpose() * q_v.transpose() * dNdx * rho_wv *
+            K_TT.noalias() += N.transpose() * vapour_flux.transpose() * dNdx *
                               specific_heat_capacity_vapour * w;
 
             double const storage_coefficient_by_water_vapor =
@@ -577,12 +576,12 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
             //
             // Latent heat term
             //
-            if (gas_phase->hasProperty(MPL::PropertyType::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
                 double const L0 =
-                    gas_phase->property(MPL::PropertyType::latent_heat)
+                    gas_phase->property(MPL::PropertyType::specific_latent_heat)
                         .template value<double>(variables, x_position, t, dt) *
                     rho_LR;
 
@@ -608,7 +607,7 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
                 K_TT.noalias() +=
                     L0 * f_Tv_D_Tv * dNdx.transpose() * dNdx * w / rho_LR;
                 // temperature equation, pressure part
-                dK_TT_dp.noalias() +=
+                K_Tp.noalias() +=
                     L0 * D_pv * dNdx.transpose() * dNdx * w / rho_LR;
             }
         }
@@ -664,7 +663,7 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::
         M_pT * (T - T_prev) / dt;
     if (gas_phase)
     {
-        if (gas_phase->hasProperty(MPL::PropertyType::latent_heat))
+        if (gas_phase->hasProperty(MPL::PropertyType::specific_latent_heat))
         {
             // Jacobian: temperature equation, pressure part
             local_Jac
@@ -999,13 +998,10 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::assemble(
 
             GlobalDimVectorType const grad_T = dNdx * T;
             GlobalDimVectorType const grad_p_cap = -dNdx * p_L;
-            // Vapour velocity
-            GlobalDimVectorType const q_v =
-                -(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::PropertyType::
-                                   specific_heat_capacity)
+                gas_phase->property(MaterialPropertyLib::specific_heat_capacity)
                     .template value<double>(variables, x_position, t, dt);
 
             local_M
@@ -1018,7 +1014,7 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::assemble(
             local_K
                 .template block<temperature_size, temperature_size>(
                     temperature_index, temperature_index)
-                .noalias() += N.transpose() * q_v.transpose() * dNdx * rho_wv *
+                .noalias() += N.transpose() * vapour_flux.transpose() * dNdx *
                               specific_heat_capacity_vapour * w;
 
             double const storage_coefficient_by_water_vapor =
@@ -1046,12 +1042,12 @@ void ThermoRichardsFlowLocalAssembler<ShapeFunction, GlobalDim>::assemble(
             //
             // Latent heat term
             //
-            if (gas_phase->hasProperty(MPL::PropertyType::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
                 double const L0 =
-                    gas_phase->property(MPL::PropertyType::latent_heat)
+                    gas_phase->property(MPL::PropertyType::specific_latent_heat)
                         .template value<double>(variables, x_position, t, dt) *
                     rho_LR;
 

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,20 +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;
 
-        // Vapour velocity
-        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 =
+        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 *
@@ -105,12 +101,12 @@ void TRMVaporDiffusionModel<DisplacementDim>::eval(
         //
         // Latent heat term
         //
-        if (gas_phase->hasProperty(MPL::PropertyType::latent_heat))
+        if (gas_phase->hasProperty(MPL::PropertyType::specific_latent_heat))
         {
             double const factor = phi * (1 - S_L_data.S_L) / rho_L_data.rho_LR;
             // The volumetric latent heat of vaporization of liquid water
             double const L0 =
-                gas_phase->property(MPL::PropertyType::latent_heat)
+                gas_phase->property(MPL::PropertyType::specific_latent_heat)
                     .template value<double>(variables, x_t.x, x_t.t, x_t.dt) *
                 rho_L_data.rho_LR;
 

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;

diff --git a/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h b/ProcessLib/ThermoRichardsMechanics/ThermoRichardsMechanicsFEM-impl.h
@@ -413,9 +413,10 @@ void ThermoRichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
     out.dK_TT_dp.noalias() =
         N.transpose() * (CD.heat_data.K_Tp_NT_V_dN.transpose() * dNdx) +
-        CD.heat_data.K_Tp_X_NTN * NTN + CD.vap_data.K_Tp_X_dNTdN * dNTdN;
+        CD.heat_data.K_Tp_X_NTN * NTN;
     out.K_Tp.noalias() =
-        dNdx.transpose() * CD.th_osmosis_data.K_Tp_Laplace * dNdx;
+        dNdx.transpose() * CD.th_osmosis_data.K_Tp_Laplace * dNdx +
+        CD.vap_data.K_Tp_X_dNTdN * dNTdN;
 
     out.K_pp.noalias() = dNdx.transpose() * CD.eq_p_data.K_pp_Laplace * dNdx +
                          CD.vap_data.K_pp_X_dNTdN * dNTdN;

diff --git a/Tests/Data/ThermoRichardsFlow/TaskCDECOVALEX2023/Decovalex-0-TRF.prj b/Tests/Data/ThermoRichardsFlow/TaskCDECOVALEX2023/Decovalex-0-TRF.prj
@@ -43,7 +43,7 @@
                             <tortuosity>0.8</tortuosity>
                         </property>
                         <property>
-                            <name>latent_heat</name>
+                            <name>specific_latent_heat</name>
                             <type>LinearWaterVapourLatentHeat</type>
                         </property>
                         <property>
@@ -192,7 +192,7 @@
                             <tortuosity>0.8</tortuosity>
                         </property>
                         <property>
-                            <name>latent_heat</name>
+                            <name>specific_latent_heat</name>
                             <type>LinearWaterVapourLatentHeat</type>
                         </property>
                         <property>
@@ -343,7 +343,7 @@
                             <tortuosity>0.8</tortuosity>
                         </property>
                         <property>
-                            <name>latent_heat</name>
+                            <name>specific_latent_heat</name>
                             <type>LinearWaterVapourLatentHeat</type>
                         </property>
                         <property>

diff --git a/...ermoRichardsFlow/TaskCDECOVALEX2023/WithPicardNonLinearSolverAndPETSc/Decovalex-0-TRF.prj b/...ermoRichardsFlow/TaskCDECOVALEX2023/WithPicardNonLinearSolverAndPETSc/Decovalex-0-TRF.prj
@@ -43,7 +43,7 @@
                             <tortuosity>0.8</tortuosity>
                         </property>
                         <property>
-                            <name>latent_heat</name>
+                            <name>specific_latent_heat</name>
                             <type>LinearWaterVapourLatentHeat</type>
                         </property>
                         <property>
@@ -192,7 +192,7 @@
                             <tortuosity>0.8</tortuosity>
                         </property>
                         <property>
-                            <name>latent_heat</name>
+                            <name>specific_latent_heat</name>
                             <type>LinearWaterVapourLatentHeat</type>
                         </property>
                         <property>
@@ -343,7 +343,7 @@
                             <tortuosity>0.8</tortuosity>
                         </property>
                         <property>
-                            <name>latent_heat</name>
+                            <name>specific_latent_heat</name>
                             <type>LinearWaterVapourLatentHeat</type>
                         </property>
                         <property>

diff --git a/...Data/ThermoRichardsFlow/TaskCDECOVALEX2023/expected_Decovalex-0_ts_10_t_864000.000000.vtu b/...Data/ThermoRichardsFlow/TaskCDECOVALEX2023/expected_Decovalex-0_ts_10_t_864000.000000.vtu
diff --git a/Tests/Data/ThermoRichardsMechanics/TaskCDECOVALEX2023/Decovalex-0.prj b/Tests/Data/ThermoRichardsMechanics/TaskCDECOVALEX2023/Decovalex-0.prj
@@ -68,7 +68,7 @@
                             <tortuosity>0.8</tortuosity>
                         </property>
                         <property>
-                            <name>latent_heat</name>
+                            <name>specific_latent_heat</name>
                             <type>LinearWaterVapourLatentHeat</type>
                         </property>
                         <property>
@@ -207,7 +207,7 @@
                             <tortuosity>0.8</tortuosity>
                         </property>
                         <property>
-                            <name>latent_heat</name>
+                            <name>specific_latent_heat</name>
                             <type>LinearWaterVapourLatentHeat</type>
                         </property>
                         <property>
@@ -348,7 +348,7 @@
                             <tortuosity>0.8</tortuosity>
                         </property>
                         <property>
-                            <name>latent_heat</name>
+                            <name>specific_latent_heat</name>
                             <type>LinearWaterVapourLatentHeat</type>
                         </property>
                         <property>