From ceb378f60986453f51b0cb0aaf7295f13af36b6b Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Tue, 16 Apr 2024 19:31:46 +0200
Subject: [PATCH] [PL/TH2M] Extract darcy velocity model
---
ProcessLib/TH2M/ConstitutiveRelations/Base.h | 8 ++++
.../ConstitutiveModels.h | 1 +
.../ConstitutiveRelations/DarcyVelocity.cpp | 43 +++++++++++++++++++
.../ConstitutiveRelations/DarcyVelocity.h | 18 ++++++++
ProcessLib/TH2M/TH2MFEM-impl.h | 28 ++++++------
5 files changed, 83 insertions(+), 15 deletions(-)
create mode 100644 ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.cpp
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Base.h b/ProcessLib/TH2M/ConstitutiveRelations/Base.h
index 8739cfbab38..e5ec55f00b4 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/Base.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/Base.h
@@ -59,6 +59,14 @@ using ReferenceTemperatureData =
using GasPressureData = BaseLib::StrongType<double, struct GasPressureTag>;
using CapillaryPressureData =
BaseLib::StrongType<double, struct CapillaryPressureTag>;
+template <int DisplacementDim>
+using GasPressureGradientData =
+ BaseLib::StrongType<GlobalDimVector<DisplacementDim>,
+ struct GasPressureGradientTag>;
+template <int DisplacementDim>
+using CapillaryPressureGradientData =
+ BaseLib::StrongType<GlobalDimVector<DisplacementDim>,
+ struct CapillaryPressureGradientTag>;
template <int DisplacementDim>
using SpecificBodyForceData =
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
index df7bee1cd8a..3ac2a484076 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h
@@ -85,6 +85,7 @@ struct ConstitutiveModels
InternalEnergyModel internal_energy_model;
EffectiveVolumetricEnthalpyModel effective_volumetric_enthalpy_model;
GravityModel<DisplacementDim> gravity_model;
+ DarcyVelocityModel<DisplacementDim> darcy_velocity_model;
FC1Model<DisplacementDim> fC_1_model;
FC2aModel fC_2a_model;
FC3aModel fC_3a_model;
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.cpp
new file mode 100644
index 00000000000..5083d8c65c2
--- /dev/null
+++ b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.cpp
@@ -0,0 +1,43 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+
+#include "DarcyVelocity.h"
+
+namespace ProcessLib::TH2M
+{
+namespace ConstitutiveRelations
+{
+template <int DisplacementDim>
+void DarcyVelocityModel<DisplacementDim>::eval(
+ CapillaryPressureGradientData<DisplacementDim> const& grad_p_cap,
+ FluidDensityData const& fluid_density_data,
+ GasPressureGradientData<DisplacementDim> const& grad_p_GR,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ SpecificBodyForceData<DisplacementDim> const& specific_body_force,
+ ViscosityData const& viscosity_data,
+ DarcyVelocityData<DisplacementDim>& darcy_velocity_data) const
+{
+ auto const k_over_mu_G =
+ permeability_data.Ki * permeability_data.k_rel_G / viscosity_data.mu_GR;
+ auto const k_over_mu_L =
+ permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR;
+
+ darcy_velocity_data.w_GS =
+ k_over_mu_G * fluid_density_data.rho_GR * specific_body_force() -
+ k_over_mu_G * grad_p_GR();
+ darcy_velocity_data.w_LS =
+ k_over_mu_L * grad_p_cap() +
+ k_over_mu_L * fluid_density_data.rho_LR * specific_body_force() -
+ k_over_mu_L * grad_p_GR();
+}
+
+template struct DarcyVelocityModel<2>;
+template struct DarcyVelocityModel<3>;
+} // namespace ConstitutiveRelations
+} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h
index 84e66d31d02..91751b830f4 100644
--- a/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h
+++ b/ProcessLib/TH2M/ConstitutiveRelations/DarcyVelocity.h
@@ -10,7 +10,10 @@
#pragma once
#include "Base.h"
+#include "FluidDensity.h"
+#include "PermeabilityData.h"
#include "ProcessLib/Reflection/ReflectionData.h"
+#include "Viscosity.h"
namespace ProcessLib::TH2M
{
@@ -32,5 +35,20 @@ struct DarcyVelocityData
R::makeReflectionData("velocity_liquid", &Self::w_LS)};
}
};
+
+template <int DisplacementDim>
+struct DarcyVelocityModel
+{
+ void eval(CapillaryPressureGradientData<DisplacementDim> const& grad_p_cap,
+ FluidDensityData const& fluid_density_data,
+ GasPressureGradientData<DisplacementDim> const& grad_p_GR,
+ PermeabilityData<DisplacementDim> const& permeability_data,
+ SpecificBodyForceData<DisplacementDim> const& specific_body_force,
+ ViscosityData const& viscosity_data,
+ DarcyVelocityData<DisplacementDim>& darcy_velocity_data) const;
+};
+
+extern template struct DarcyVelocityModel<2>;
+extern template struct DarcyVelocityModel<3>;
} // namespace ConstitutiveRelations
} // namespace ProcessLib::TH2M
diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h
index cc118746435..c87165fd500 100644
--- a/ProcessLib/TH2M/TH2MFEM-impl.h
+++ b/ProcessLib/TH2M/TH2MFEM-impl.h
@@ -162,6 +162,10 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
GlobalDimVectorType const gradpGR = gradNp * gas_pressure;
GlobalDimVectorType const gradpCap = gradNp * capillary_pressure;
GlobalDimVectorType const gradT = gradNp * temperature;
+ ConstitutiveRelations::GasPressureGradientData<DisplacementDim> const
+ grad_p_GR{gradpGR};
+ ConstitutiveRelations::CapillaryPressureGradientData<
+ DisplacementDim> const grad_p_cap{gradpCap};
// medium properties
models.elastic_tangent_stiffness_model.eval({pos, t, dt}, T_data,
@@ -487,21 +491,15 @@ TH2MLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
// Derivatives for Jacobian
// ---------------------------------------------------------------------
- auto const& b = this->process_data_.specific_body_force;
- GlobalDimMatrixType const k_over_mu_G =
- ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_G /
- ip_cv.viscosity_data.mu_GR;
- GlobalDimMatrixType const k_over_mu_L =
- ip_out.permeability_data.Ki * ip_out.permeability_data.k_rel_L /
- ip_cv.viscosity_data.mu_LR;
-
- ip_out.darcy_velocity_data.w_GS =
- k_over_mu_G * ip_out.fluid_density_data.rho_GR * b -
- k_over_mu_G * gradpGR;
- ip_out.darcy_velocity_data.w_LS =
- k_over_mu_L * gradpCap +
- k_over_mu_L * ip_out.fluid_density_data.rho_LR * b -
- k_over_mu_L * gradpGR;
+ models.darcy_velocity_model.eval(
+ grad_p_cap,
+ ip_out.fluid_density_data,
+ grad_p_GR,
+ ip_out.permeability_data,
+ ConstitutiveRelations::SpecificBodyForceData<DisplacementDim>{
+ this->process_data_.specific_body_force},
+ ip_cv.viscosity_data,
+ ip_out.darcy_velocity_data);
models.fT_2_model.eval(ip_out.darcy_velocity_data,
ip_out.enthalpy_data,
--
GitLab