[MPL] Added a new function to get fluid thermal expansity

MaterialLib/MPL/Components/GetThermalExpansivity.cpp

@@ -0,0 +1,43 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * Created on August 16, 2019, 3:40 PM
+ */
+
+#include "GetThermalExpansivity.h"
+
+#include "MaterialLib/MPL/Phase.h"
+
+namespace MaterialPropertyLib
+{
+class Phase;
+
+double getThermalExpansivity(Phase const& phase,
+                             VariableArray const& vars,
+                             const double density)
+{
+    auto const thermal_expansivity_ptr =
+        &phase.property(MaterialPropertyLib::PropertyType::thermal_expansivity);
+
+    // The thermal expansivity is explicitly given in the project file.
+    if (thermal_expansivity_ptr)
+    {
+        return (*thermal_expansivity_ptr).template value<double>(vars);
+    }
+
+    // The thermal expansivity calculated by the density model directly.
+    return (density == 0.0)
+               ? 0.0
+               : -phase.property(MaterialPropertyLib::PropertyType::density)
+                         .template dValue<double>(
+                             vars, MaterialPropertyLib::Variable::temperature) /
+                     density;
+}
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Components/GetThermalExpansivity.h

@@ -0,0 +1,44 @@
+/**
+ * \file
+ *
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ *
+ * Created on August 16, 2019, 3:40 PM
+ */
+
+#pragma once
+
+#include "MaterialLib/MPL/VariableType.h"  // for VariableArray
+
+namespace MaterialPropertyLib
+{
+class Phase;
+
+/**
+ * It gets the thermal expansion coefficient.
+ *
+ * If the the thermal expansion coefficient is given in the project file via
+ * the media property of thermal_expansivity, e.g
+ * \verbatim
+ *     <property>
+ *       <name>thermal_expansivity</name>
+ *       <type>Constant</type>
+ *       <value>2.07e-4</value>
+ *     </property>
+ * \endverbatim
+ * it returns the value of the given property. Otherwise it returns the value
+ * computed from the density model by the following formula
+ * \f[
+ *      (\frac{\partial \rho}{\partial T})/\rho
+ * \f]
+ * where \f$\rho\f$ is the density, \f$T\f$ is the temperature.
+ */
+double getThermalExpansivity(Phase const& phase,
+                             VariableArray const& vars,
+                             const double density);
+}  // namespace MaterialPropertyLib

ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h

@@ -13,6 +13,8 @@
 #include "ThermoHydroMechanicsFEM.h"
 
 #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
+
+#include "MaterialLib/MPL/Components/GetThermalExpansivity.h"
 #include "MaterialLib/MPL/Medium.h"
 #include "MaterialLib/MPL/Property.h"
 #include "MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.h"
@@ -256,22 +258,9 @@ void ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
             liquid_phase.property(MaterialPropertyLib::PropertyType::density)
                 .template value<double>(vars);
 
-        auto const fluid_thermal_expansivity_ptr = &liquid_phase.property(
-            MaterialPropertyLib::PropertyType::thermal_expansivity);
-
         double const fluid_volumetric_thermal_expansion_coefficient =
-            fluid_thermal_expansivity_ptr
-                ? (*fluid_thermal_expansivity_ptr).template value<double>(vars)
-                : (fluid_density == 0.0)
-                      ? 0.0
-                      : -liquid_phase
-                                .property(
-                                    MaterialPropertyLib::PropertyType::density)
-                                .template dValue<double>(
-                                    vars,
-                                    MaterialPropertyLib::Variable::
-                                        temperature) /
-                            fluid_density;
+            MaterialPropertyLib::getThermalExpansivity(liquid_phase, vars,
+                                                       fluid_density);
 
         // Use the viscosity model to compute the viscosity
         auto const viscosity =
@@ -286,7 +275,7 @@ void ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
             MathLib::KelvinVector::KelvinVectorDimensions<
                 DisplacementDim>::value>::identity2;
 
-         //TODO: Change dT to time step wise increment
+        // TODO (Wenqing) : Change dT to time step wise increment
         double const delta_T(T_int_pt - T0);
         double const thermal_strain =
             solid_linear_thermal_expansion_coefficient * delta_T;