From 27604c2f60b475fc4d3d2dc0126cd155090e3610 Mon Sep 17 00:00:00 2001
From: Boyan Meng <meng.boyan@ufz.de>
Date: Wed, 30 Oct 2019 15:56:30 +0100
Subject: [PATCH] consider the liquid phase and change the calculation of the
 volumtric heat capacity accordingly

---
 .../HeatTransportBHELocalAssemblerSoil-impl.h | 34 ++++++++++++-------
 1 file changed, 22 insertions(+), 12 deletions(-)

diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil-impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil-impl.h
index 1ae416e3dae..1550f7fdde4 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil-impl.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil-impl.h
@@ -14,6 +14,7 @@
 #include <vector>
 
 #include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.h"
 #include "MaterialLib/MPL/Utils/FormEigenTensor.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
@@ -113,37 +114,45 @@ void HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>::
                     MaterialPropertyLib::PropertyType::thermal_conductivity)
                 .template value<double>(vars, pos, t);
 
+        auto const density_s =
+            solid_phase.property(MaterialPropertyLib::PropertyType::density)
+                .template value<double>(vars, pos, t);
+
         auto const heat_capacity_s =
             solid_phase
                 .property(
                     MaterialPropertyLib::PropertyType::specific_heat_capacity)
                 .template value<double>(vars, pos, t);
 
-        auto const heat_capacity_f =
+        auto const k_f =
             liquid_phase
                 .property(
-                    MaterialPropertyLib::PropertyType::specific_heat_capacity)
-                .template value<double>(vars, pos, t);
-
-        auto const density_s =
-            solid_phase.property(MaterialPropertyLib::PropertyType::density)
+                    MaterialPropertyLib::PropertyType::thermal_conductivity)
                 .template value<double>(vars, pos, t);
 
         auto const density_f =
             liquid_phase.property(MaterialPropertyLib::PropertyType::density)
                 .template value<double>(vars, pos, t);
-        // for now only using the solid phase parameters
 
-        /*auto const velocity =
+        auto const heat_capacity_f =
             liquid_phase
-                .property(MaterialPropertyLib::PropertyType::phase_velocity)
-                .template value<GlobalDimVectorType>(vars, pos, t);*/
+                .property(
+                    MaterialPropertyLib::PropertyType::specific_heat_capacity)
+                .template value<double>(vars, pos, t);
+
+        auto const porosity =
+            medium.property(MaterialPropertyLib::PropertyType::porosity)
+                .template value<double>(vars, pos, t);
+
+        // for now only using the solid and liquid phase parameters
+
         auto const velocity_arr =
             liquid_phase
                 .property(MaterialPropertyLib::PropertyType::phase_velocity)
                 .value(vars, pos, t);
         auto velocity = Eigen::Map<Eigen::Matrix<double, 1, 3> const>(
             std::get<MaterialPropertyLib::Vector>(velocity_arr).data(), 1, 3);
+
         // assemble Conductance matrix
         local_K.noalias() +=
             (dNdx.transpose() * dNdx * k_s +
@@ -151,8 +160,9 @@ void HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>::
             w;
 
         // assemble Mass matrix
-        local_M.noalias() +=
-            N.transpose() * N * w * density_s * heat_capacity_s;
+        local_M.noalias() += N.transpose() * N * w *
+                             (density_s * heat_capacity_s * (1 - porosity) +
+                              density_f * heat_capacity_f * porosity);
     }
 
     // debugging
-- 
GitLab