diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Advection.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Advection.cpp new file mode 100644 index 00000000000..a236f888782 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/Advection.cpp @@ -0,0 +1,75 @@ +/** + * \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 "Advection.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +template +void AdvectionModel::eval( + ConstituentDensityData const& constituent_density_data, + PermeabilityData const& permeability_data, + PureLiquidDensityData const& rho_W_LR, + ViscosityData const& viscosity_data, + AdvectionData& advection_data) const +{ + GlobalDimMatrix const k_over_mu_G = + permeability_data.Ki * permeability_data.k_rel_G / viscosity_data.mu_GR; + GlobalDimMatrix const k_over_mu_L = + permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR; + + advection_data.advection_C_G = + constituent_density_data.rho_C_GR * k_over_mu_G; + advection_data.advection_C_L = + constituent_density_data.rho_C_LR * k_over_mu_L; + advection_data.advection_W_G = + constituent_density_data.rho_W_GR * k_over_mu_G; + advection_data.advection_W_L = rho_W_LR() * k_over_mu_L; +} + +template +void AdvectionModel::dEval( + ConstituentDensityData const& constituent_density_data, + PermeabilityData const& permeability_data, + ViscosityData const& viscosity_data, + SaturationDataDeriv const& dS_L_dp_cap, + PhaseTransitionData const& phase_transition_data, + AdvectionDerivativeData& advection_d_data) const +{ + GlobalDimMatrix const k_over_mu_G = + permeability_data.Ki * permeability_data.k_rel_G / viscosity_data.mu_GR; + GlobalDimMatrix const k_over_mu_L = + permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR; + + GlobalDimMatrix const dk_over_mu_G_dp_cap = + permeability_data.Ki * permeability_data.dk_rel_G_dS_L * dS_L_dp_cap() / + viscosity_data.mu_GR; + GlobalDimMatrix const dk_over_mu_L_dp_cap = + permeability_data.Ki * permeability_data.dk_rel_L_dS_L * dS_L_dp_cap() / + viscosity_data.mu_LR; + + advection_d_data.dadvection_C_dp_GR = + phase_transition_data.drho_C_GR_dp_GR * k_over_mu_G + // + rhoCGR * (dk_over_mu_G_dp_GR = 0) + // + rhoCLR * (dk_over_mu_L_dp_GR = 0) + + phase_transition_data.drho_C_LR_dp_GR * k_over_mu_L; + + advection_d_data.dadvection_C_dp_cap = + //(drho_C_GR_dp_cap = 0) * k_over_mu_G + constituent_density_data.rho_C_GR * dk_over_mu_G_dp_cap + + (-phase_transition_data.drho_C_LR_dp_LR) * k_over_mu_L + + constituent_density_data.rho_C_LR * dk_over_mu_L_dp_cap; +} + +template struct AdvectionModel<2>; +template struct AdvectionModel<3>; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Advection.h b/ProcessLib/TH2M/ConstitutiveRelations/Advection.h new file mode 100644 index 00000000000..316f52b82ab --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/Advection.h @@ -0,0 +1,61 @@ +/** + * \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 + */ + +#pragma once + +#include "Base.h" +#include "ConstitutiveDensity.h" +#include "PermeabilityData.h" +#include "PhaseTransitionData.h" +#include "PureLiquidDensity.h" +#include "Saturation.h" +#include "Viscosity.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +template +struct AdvectionData +{ + GlobalDimMatrix advection_C_G; + GlobalDimMatrix advection_C_L; + GlobalDimMatrix advection_W_G; + GlobalDimMatrix advection_W_L; +}; + +template +struct AdvectionDerivativeData +{ + GlobalDimMatrix dadvection_C_dp_GR; + GlobalDimMatrix dadvection_C_dp_cap; +}; + +template +struct AdvectionModel +{ + void eval(ConstituentDensityData const& constituent_density_data, + PermeabilityData const& permeability_data, + PureLiquidDensityData const& rho_W_LR, + ViscosityData const& viscosity_data, + AdvectionData& advection_data) const; + + void dEval( + ConstituentDensityData const& constituent_density_data, + PermeabilityData const& permeability_data, + ViscosityData const& viscosity_data, + SaturationDataDeriv const& dS_L_dp_cap, + PhaseTransitionData const& phase_transition_data, + AdvectionDerivativeData& advection_d_data) const; +}; + +extern template struct AdvectionModel<2>; +extern template struct AdvectionModel<3>; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Base.h b/ProcessLib/TH2M/ConstitutiveRelations/Base.h index 0667e3730d9..70c0467e6cf 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/Base.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/Base.h @@ -59,6 +59,23 @@ using ReferenceTemperatureData = using GasPressureData = BaseLib::StrongType; using CapillaryPressureData = BaseLib::StrongType; +template +using GasPressureGradientData = + BaseLib::StrongType, + struct GasPressureGradientTag>; +template +using CapillaryPressureGradientData = + BaseLib::StrongType, + struct CapillaryPressureGradientTag>; +template +using TemperatureGradientData = + BaseLib::StrongType, + struct TemperatureGradientTag>; + +template +using SpecificBodyForceData = + BaseLib::StrongType, + struct SpecificBodyForceTag>; } // namespace ConstitutiveRelations } // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/CEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/CEquation.cpp new file mode 100644 index 00000000000..6242c36efd4 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/CEquation.cpp @@ -0,0 +1,467 @@ +/** + * \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 "CEquation.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +template +void FC1Model::eval( + AdvectionData const& advection_data, + FluidDensityData const& fluid_density_data, + FC1Data& fC_1) const +{ + fC_1.A = advection_data.advection_C_G * fluid_density_data.rho_GR + + advection_data.advection_C_L * fluid_density_data.rho_LR; +} + +template struct FC1Model<2>; +template struct FC1Model<3>; + +void FC2aModel::eval(BiotData const biot_data, + CapillaryPressureData const pCap, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + SolidCompressibilityData const beta_p_SR, + FC2aData& fC_2a) const +{ + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + double const rho_C_FR = S_G * constituent_density_data.rho_C_GR + + S_L * constituent_density_data.rho_C_LR; + fC_2a.a = + porosity_data.phi * (constituent_density_data.rho_C_LR - + constituent_density_data.rho_C_GR) - + rho_C_FR * pCap() * (biot_data() - porosity_data.phi) * beta_p_SR(); +} + +void FC2aModel::dEval(BiotData const& biot_data, + CapillaryPressureData const pCap, + ConstituentDensityData const& constituent_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + PorosityDerivativeData const& porosity_d_data, + SaturationData const& S_L_data, + SaturationDataDeriv const& dS_L_dp_cap, + SolidCompressibilityData const& beta_p_SR, + FC2aDerivativeData& dfC_2a) const +{ + double const S_L = S_L_data.S_L; + double const S_G = 1. - S_L; + + double const drho_C_FR_dp_GR = + /*(dS_G_dp_GR = 0) * constituent_density_data.rho_C_GR +*/ + S_G * phase_transition_data.drho_C_GR_dp_GR + + /*(dS_L_dp_GR = 0) * constituent_density_data.rho_C_LR +*/ + S_L * phase_transition_data.drho_C_LR_dp_GR; + + dfC_2a.dp_GR = -porosity_data.phi * phase_transition_data.drho_C_GR_dp_GR - + drho_C_FR_dp_GR * pCap() * + (biot_data() - porosity_data.phi) * beta_p_SR(); + + double const dS_G_dp_cap = -dS_L_dp_cap(); + double const rho_C_FR = S_G * constituent_density_data.rho_C_GR + + S_L * constituent_density_data.rho_C_LR; + + // TODO (naumov) Extend for partially saturated media. + constexpr double drho_C_GR_dp_cap = 0; + + double const drho_C_FR_dp_cap = + dS_G_dp_cap * constituent_density_data.rho_C_GR + + S_G * drho_C_GR_dp_cap + + dS_L_dp_cap() * constituent_density_data.rho_C_LR - + S_L * phase_transition_data.drho_C_LR_dp_LR; + + dfC_2a.dp_cap = + porosity_data.phi * + (-phase_transition_data.drho_C_LR_dp_LR - drho_C_GR_dp_cap) - + drho_C_FR_dp_cap * pCap() * (biot_data() - porosity_data.phi) * + beta_p_SR() + + rho_C_FR * (biot_data() - porosity_data.phi) * beta_p_SR(); + + double const drho_C_FR_dT = S_G * phase_transition_data.drho_C_GR_dT + + S_L * phase_transition_data.drho_C_LR_dT; + dfC_2a.dT = porosity_d_data.dphi_dT * (constituent_density_data.rho_C_LR - + constituent_density_data.rho_C_GR) + + porosity_data.phi * (phase_transition_data.drho_C_LR_dT - + phase_transition_data.drho_C_GR_dT) - + drho_C_FR_dT * pCap() * (biot_data() - porosity_data.phi) * + beta_p_SR() + + rho_C_FR * pCap() * porosity_d_data.dphi_dT * beta_p_SR(); +} + +void FC3aModel::eval( + double const dt, + ConstituentDensityData const& constituent_density_data, + PrevState const& constituent_density_data_prev, + SaturationData const& S_L_data, + FC3aData& fC_3a) const +{ + if (dt == 0.) + { + fC_3a.a = 0; + return; + } + + double const rho_C_GR_dot = (constituent_density_data.rho_C_GR - + constituent_density_data_prev->rho_C_GR) / + dt; + double const rho_C_LR_dot = (constituent_density_data.rho_C_LR - + constituent_density_data_prev->rho_C_LR) / + dt; + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + fC_3a.a = S_G * rho_C_GR_dot + S_L * rho_C_LR_dot; +} + +void FC3aModel::dEval( + double const dt, + ConstituentDensityData const& constituent_density_data, + PrevState const& constituent_density_data_prev, + PhaseTransitionData const& phase_transition_data, + SaturationData const& S_L_data, + SaturationDataDeriv const& dS_L_dp_cap, + FC3aDerivativeData& dfC_3a) const +{ + if (dt == 0.) + { + dfC_3a.dp_GR = 0.; + dfC_3a.dp_cap = 0.; + dfC_3a.dT = 0.; + return; + } + double const rho_C_GR_dot = (constituent_density_data.rho_C_GR - + constituent_density_data_prev->rho_C_GR) / + dt; + double const rho_C_LR_dot = (constituent_density_data.rho_C_LR - + constituent_density_data_prev->rho_C_LR) / + dt; + + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + dfC_3a.dp_GR = + /*(dS_G_dp_GR = 0) * rho_C_GR_dot +*/ + S_G * phase_transition_data.drho_C_GR_dp_GR / dt + + /*(dS_L_dp_GR = 0) * rho_C_LR_dot +*/ + S_L * phase_transition_data.drho_C_LR_dp_GR / dt; + + double const dS_G_dp_cap = -dS_L_dp_cap(); + // TODO (naumov) Extend for partially saturated media. + constexpr double drho_C_GR_dp_cap = 0; + + dfC_3a.dp_cap = dS_G_dp_cap * rho_C_GR_dot + S_G * drho_C_GR_dp_cap / dt + + dS_L_dp_cap() * rho_C_LR_dot - + S_L * phase_transition_data.drho_C_LR_dp_LR / dt; + + dfC_3a.dT = S_G * phase_transition_data.drho_C_GR_dT / dt + + S_L * phase_transition_data.drho_C_LR_dT / dt; +} + +template +void FC4LCpGModel::eval( + AdvectionData const& advection_data, + FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + FC4LCpGData& fC_4_LCpG) const +{ + GlobalDimMatrix const advection_C = + advection_data.advection_C_G + advection_data.advection_C_L; + + double const sD_G = phase_transition_data.diffusion_coefficient_vapour; + double const sD_L = phase_transition_data.diffusion_coefficient_solute; + + double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi; + double const phi_L = S_L_data.S_L * porosity_data.phi; + + double const diffusion_CGpGR = -phi_G * fluid_density_data.rho_GR * sD_G * + phase_transition_data.dxmWG_dpGR; + double const diffusion_CLpGR = -phi_L * fluid_density_data.rho_LR * sD_L * + phase_transition_data.dxmWL_dpGR; + + double const diffusion_C_pGR = diffusion_CGpGR + diffusion_CLpGR; + + auto const I = + Eigen::Matrix::Identity(); + fC_4_LCpG.L.noalias() = diffusion_C_pGR * I + advection_C; +} + +template +void FC4LCpGModel::dEval( + PermeabilityData const& permeability_data, + ViscosityData const& viscosity_data, + PhaseTransitionData const& phase_transition_data, + AdvectionDerivativeData const& advection_d_data, + FC4LCpGDerivativeData& dfC_4_LCpG) const +{ + dfC_4_LCpG.dp_GR = advection_d_data.dadvection_C_dp_GR + // + ddiffusion_C_p_dp_GR TODO (naumov) + ; + + dfC_4_LCpG.dp_cap = advection_d_data.dadvection_C_dp_cap + // + ddiffusion_C_p_dp_cap TODO (naumov) + ; + + GlobalDimMatrix const k_over_mu_G = + permeability_data.Ki * permeability_data.k_rel_G / viscosity_data.mu_GR; + GlobalDimMatrix const k_over_mu_L = + permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR; + + dfC_4_LCpG.dT = phase_transition_data.drho_C_GR_dT * k_over_mu_G + + phase_transition_data.drho_C_LR_dT * k_over_mu_L + // + ddiffusion_C_p_dT TODO (naumov) + ; +} + +template struct FC4LCpGModel<2>; +template struct FC4LCpGModel<3>; + +template +void FC4LCpCModel::eval( + AdvectionData const& advection_data, + FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + FC4LCpCData& fC_4_LCpC) const +{ + double const sD_G = phase_transition_data.diffusion_coefficient_vapour; + double const sD_L = phase_transition_data.diffusion_coefficient_solute; + + double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi; + double const phi_L = S_L_data.S_L * porosity_data.phi; + + double const diffusion_CGpCap = -phi_G * fluid_density_data.rho_GR * sD_G * + phase_transition_data.dxmWG_dpCap; + double const diffusion_CLpCap = -phi_L * fluid_density_data.rho_LR * sD_L * + phase_transition_data.dxmWL_dpCap; + + double const diffusion_C_pCap = diffusion_CGpCap + diffusion_CLpCap; + + auto const I = + Eigen::Matrix::Identity(); + + fC_4_LCpC.L.noalias() = diffusion_C_pCap * I - advection_data.advection_C_L; +} + +template +void FC4LCpCModel::dEval( + ConstituentDensityData const& constituent_density_data, + PermeabilityData const& permeability_data, + PhaseTransitionData const& phase_transition_data, + SaturationDataDeriv const& dS_L_dp_cap, + ViscosityData const& viscosity_data, + FC4LCpCDerivativeData& dfC_4_LCpC) const +{ + ////// Diffusion Part ///// + // TODO (naumov) d(diffusion_C_pCap)/dX for dxmW*/d* != 0 + + ////// Advection part ///// + GlobalDimMatrix const k_over_mu_L = + permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR; + + dfC_4_LCpC.dp_GR = phase_transition_data.drho_C_LR_dp_GR * k_over_mu_L + //+ rhoCLR * (dk_over_mu_L_dp_GR = 0) + ; + + auto const dk_over_mu_L_dp_cap = permeability_data.Ki * + permeability_data.dk_rel_L_dS_L * + dS_L_dp_cap() / viscosity_data.mu_LR; + + dfC_4_LCpC.dp_cap = -phase_transition_data.drho_C_LR_dp_LR * k_over_mu_L + + constituent_density_data.rho_C_LR * dk_over_mu_L_dp_cap; + + dfC_4_LCpC.dT = phase_transition_data.drho_W_LR_dT * k_over_mu_L + //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T)) + ; +} + +template struct FC4LCpCModel<2>; +template struct FC4LCpCModel<3>; + +template +void FC4LCTModel::eval( + FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + FC4LCTData& fC_4_LCT) const +{ + double const sD_G = phase_transition_data.diffusion_coefficient_vapour; + double const sD_L = phase_transition_data.diffusion_coefficient_solute; + + double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi; + double const phi_L = S_L_data.S_L * porosity_data.phi; + + double const diffusion_C_G_T = -phi_G * fluid_density_data.rho_GR * sD_G * + phase_transition_data.dxmWG_dT; + double const diffusion_C_L_T = -phi_L * fluid_density_data.rho_LR * sD_L * + phase_transition_data.dxmWL_dT; + + double const diffusion_C_T = diffusion_C_G_T + diffusion_C_L_T; + + auto const I = + Eigen::Matrix::Identity(); + + fC_4_LCT.L.noalias() = diffusion_C_T * I; +} + +template struct FC4LCTModel<2>; +template struct FC4LCTModel<3>; + +void FC4MCpGModel::eval(BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + SolidCompressibilityData const& beta_p_SR, + FC4MCpGData& fC_4_MCpG) const +{ + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + double const rho_C_FR = S_G * constituent_density_data.rho_C_GR + + S_L * constituent_density_data.rho_C_LR; + + fC_4_MCpG.m = rho_C_FR * (biot_data() - porosity_data.phi) * beta_p_SR(); +} + +void FC4MCpGModel::dEval(BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + PorosityDerivativeData const& porosity_d_data, + SaturationData const& S_L_data, + SolidCompressibilityData const& beta_p_SR, + FC4MCpGDerivativeData& dfC_4_MCpG) const +{ + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + double const rho_C_FR = S_G * constituent_density_data.rho_C_GR + + S_L * constituent_density_data.rho_C_LR; + + double const drho_C_FR_dp_GR = + /*(dS_G_dp_GR = 0) * constituent_density_data.rho_C_GR +*/ + S_G * phase_transition_data.drho_C_GR_dp_GR + + /*(dS_L_dp_GR = 0) * constituent_density_data.rho_C_LR +*/ + S_L * phase_transition_data.drho_C_LR_dp_GR; + + dfC_4_MCpG.dp_GR = + drho_C_FR_dp_GR * (biot_data() - porosity_data.phi) * beta_p_SR(); + + double const drho_C_FR_dT = S_G * phase_transition_data.drho_C_GR_dT + + S_L * phase_transition_data.drho_C_LR_dT; + dfC_4_MCpG.dT = + drho_C_FR_dT * (biot_data() - porosity_data.phi) * beta_p_SR() - + rho_C_FR * porosity_d_data.dphi_dT * beta_p_SR(); +} + +void FC4MCpCModel::eval(BiotData const& biot_data, + CapillaryPressureData const pCap, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + PrevState const& S_L_data_prev, + SaturationData const& S_L_data, + SolidCompressibilityData const& beta_p_SR, + FC4MCpCData& fC_4_MCpC) const +{ + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + double const rho_C_FR = S_G * constituent_density_data.rho_C_GR + + S_L * constituent_density_data.rho_C_LR; + + fC_4_MCpC.m = + -rho_C_FR * (biot_data() - porosity_data.phi) * beta_p_SR() * S_L; + + fC_4_MCpC.ml = + (porosity_data.phi * (constituent_density_data.rho_C_LR - + constituent_density_data.rho_C_GR) - + rho_C_FR * pCap() * (biot_data() - porosity_data.phi) * beta_p_SR()) * + (S_L - S_L_data_prev->S_L); +} + +template +void FC4MCTModel::eval( + BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + SolidThermalExpansionData const& s_therm_exp_data, + FC4MCTData& fC_4_MCT) const +{ + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + double const rho_C_FR = S_G * constituent_density_data.rho_C_GR + + S_L * constituent_density_data.rho_C_LR; + + fC_4_MCT.m = -rho_C_FR * (biot_data() - porosity_data.phi) * + s_therm_exp_data.beta_T_SR; +} + +template +void FC4MCTModel::dEval( + BiotData const& biot_data, + [[maybe_unused]] ConstituentDensityData const& constituent_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + [[maybe_unused]] PorosityDerivativeData const& porosity_d_data, + SaturationData const& S_L_data, + SolidThermalExpansionData const& s_therm_exp_data, + FC4MCTDerivativeData& dfC_4_MCT) const +{ + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; +#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION + double const rho_C_FR = S_G * constituent_density_data.rho_C_GR + + S_L * constituent_density_data.rho_C_LR; +#endif + + double const drho_C_FR_dT = S_G * phase_transition_data.drho_C_GR_dT + + S_L * phase_transition_data.drho_C_LR_dT; + + dfC_4_MCT.dT = drho_C_FR_dT * (biot_data() - porosity_data.phi) * + s_therm_exp_data.beta_T_SR +#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION + + rho_C_FR * (biot_data() - porosity_d_data.dphi_dT) * + s_therm_exp_data.beta_T_SR +#endif + ; +} + +template struct FC4MCTModel<2>; +template struct FC4MCTModel<3>; + +void FC4MCuModel::eval(BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + SaturationData const& S_L_data, + FC4MCuData& fC_4_MCu) const +{ + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + double const rho_C_FR = S_G * constituent_density_data.rho_C_GR + + S_L * constituent_density_data.rho_C_LR; + + fC_4_MCu.m = rho_C_FR * biot_data(); +} + +void FC4MCuModel::dEval(BiotData const& biot_data, + PhaseTransitionData const& phase_transition_data, + SaturationData const& S_L_data, + FC4MCuDerivativeData& dfC_4_MCu) const +{ + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + double const drho_C_FR_dT = S_G * phase_transition_data.drho_C_GR_dT + + S_L * phase_transition_data.drho_C_LR_dT; + dfC_4_MCu.dT = drho_C_FR_dT * biot_data(); +} +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/CEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/CEquation.h new file mode 100644 index 00000000000..c4de9c36a64 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/CEquation.h @@ -0,0 +1,301 @@ +/** + * \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 + */ + +#pragma once + +#include "Advection.h" +#include "Base.h" +#include "Biot.h" +#include "ConstitutiveDensity.h" +#include "FluidDensity.h" +#include "PermeabilityData.h" +#include "PhaseTransitionData.h" +#include "Porosity.h" +#include "Saturation.h" +#include "SolidCompressibility.h" +#include "Viscosity.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +template +struct FC1Data +{ + GlobalDimMatrix A; +}; + +template +struct FC1Model +{ + void eval(AdvectionData const& advection_data, + FluidDensityData const& fluid_density_data, + FC1Data& fC_1) const; +}; + +extern template struct FC1Model<2>; +extern template struct FC1Model<3>; + +struct FC2aData +{ + double a = nan; +}; + +struct FC2aDerivativeData +{ + double dp_GR = nan; + double dp_cap = nan; + double dT = nan; +}; + +struct FC2aModel +{ + void eval(BiotData const biot_data, + CapillaryPressureData const pCap, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + SolidCompressibilityData const beta_p_SR, + FC2aData& fC_2a) const; + + void dEval(BiotData const& biot_data, + CapillaryPressureData const pCap, + ConstituentDensityData const& constituent_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + PorosityDerivativeData const& porosity_d_data, + SaturationData const& S_L_data, + SaturationDataDeriv const& dS_L_dp_cap, + SolidCompressibilityData const& beta_p_SR, + FC2aDerivativeData& dfC_2a) const; +}; + +struct FC3aData +{ + double a = nan; +}; + +struct FC3aDerivativeData +{ + double dp_GR = nan; + double dp_cap = nan; + double dT = nan; +}; + +struct FC3aModel +{ + void eval( + double const dt, + ConstituentDensityData const& constituent_density_data, + PrevState const& constituent_density_data_prev, + SaturationData const& S_L_data, + FC3aData& fC_3a) const; + + void dEval( + double const dt, + ConstituentDensityData const& constituent_density_data, + PrevState const& constituent_density_data_prev, + PhaseTransitionData const& phase_transition_data, + SaturationData const& S_L_data, + SaturationDataDeriv const& dS_L_dp_cap, + FC3aDerivativeData& dfC_3a) const; +}; + +template +struct FC4LCpGData +{ + GlobalDimMatrix L; +}; + +template +struct FC4LCpGDerivativeData +{ + GlobalDimMatrix dp_GR; + GlobalDimMatrix dp_cap; + GlobalDimMatrix dT; +}; + +template +struct FC4LCpGModel +{ + void eval(AdvectionData const& advection_data, + FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + FC4LCpGData& fC_4_LCpG) const; + + void dEval(PermeabilityData const& permeability_data, + ViscosityData const& viscosity_data, + PhaseTransitionData const& phase_transition_data, + AdvectionDerivativeData const& advection_d_data, + FC4LCpGDerivativeData& dfC_4_LCpG) const; +}; + +extern template struct FC4LCpGModel<2>; +extern template struct FC4LCpGModel<3>; + +template +struct FC4LCpCData +{ + GlobalDimMatrix L; +}; + +template +struct FC4LCpCDerivativeData +{ + GlobalDimMatrix dp_GR; + GlobalDimMatrix dp_cap; + GlobalDimMatrix dT; +}; + +template +struct FC4LCpCModel +{ + void eval(AdvectionData const& advection_data, + FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + FC4LCpCData& fC_4_LCpC) const; + + void dEval(ConstituentDensityData const& constituent_density_data, + PermeabilityData const& permeability_data, + PhaseTransitionData const& phase_transition_data, + SaturationDataDeriv const& dS_L_dp_cap, + ViscosityData const& viscosity_data, + FC4LCpCDerivativeData& dfC_4_LCpC) const; +}; + +extern template struct FC4LCpCModel<2>; +extern template struct FC4LCpCModel<3>; + +template +struct FC4LCTData +{ + GlobalDimMatrix L; +}; + +template +struct FC4LCTModel +{ + void eval(FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + FC4LCTData& fC_4_LCT) const; +}; + +struct FC4MCpGData +{ + double m = nan; +}; + +struct FC4MCpGDerivativeData +{ + double dp_GR = nan; + double dT = nan; +}; + +struct FC4MCpGModel +{ + void eval(BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + SolidCompressibilityData const& beta_p_SR, + FC4MCpGData& fC_4_MCpG) const; + + void dEval(BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + PorosityDerivativeData const& porosity_d_data, + SaturationData const& S_L_data, + SolidCompressibilityData const& beta_p_SR, + FC4MCpGDerivativeData& dfC_4_MCpG) const; +}; + +struct FC4MCpCData +{ + double m = nan; + double ml = nan; +}; + +struct FC4MCpCModel +{ + void eval(BiotData const& biot_data, + CapillaryPressureData const pCap, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + PrevState const& S_L_data_prev, + SaturationData const& S_L_data, + SolidCompressibilityData const& beta_p_SR, + FC4MCpCData& fC_4_MCpC) const; +}; + +struct FC4MCTData +{ + double m = nan; +}; + +struct FC4MCTDerivativeData +{ + double dT = nan; +}; + +template +struct FC4MCTModel +{ + void eval( + BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + SolidThermalExpansionData const& s_therm_exp_data, + FC4MCTData& fC_4_MCT) const; + + void dEval( + BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + PorosityDerivativeData const& porosity_d_data, + SaturationData const& S_L_data, + SolidThermalExpansionData const& s_therm_exp_data, + FC4MCTDerivativeData& dfC_4_MCT) const; +}; + +struct FC4MCuData +{ + double m = nan; +}; + +struct FC4MCuDerivativeData +{ + double dT = nan; +}; + +struct FC4MCuModel +{ + void eval(BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + SaturationData const& S_L_data, + FC4MCuData& fC_4_MCu) const; + + void dEval(BiotData const& biot_data, + PhaseTransitionData const& phase_transition_data, + SaturationData const& S_L_data, + FC4MCuDerivativeData& dfC_4_MCu) const; +}; + +extern template struct FC4MCTModel<2>; +extern template struct FC4MCTModel<3>; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h index 7c4611cea70..38f09a60085 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveData.h @@ -9,8 +9,10 @@ #pragma once +#include "Advection.h" #include "Biot.h" #include "Bishops.h" +#include "CEquation.h" #include "ConstitutiveDensity.h" #include "DarcyVelocity.h" #include "DiffusionVelocity.h" @@ -18,6 +20,7 @@ #include "Enthalpy.h" #include "EquivalentPlasticStrainData.h" #include "FluidDensity.h" +#include "Gravity.h" #include "InternalEnergy.h" #include "MassMoleFractions.h" #include "MechanicalStrain.h" @@ -35,10 +38,13 @@ #include "SolidMechanics.h" #include "SolidThermalExpansion.h" #include "Swelling.h" +#include "TEquation.h" #include "ThermalConductivity.h" #include "TotalStress.h" +#include "UEquation.h" #include "VapourPartialPressure.h" #include "Viscosity.h" +#include "WEquation.h" namespace ProcessLib::TH2M { @@ -99,7 +105,8 @@ struct OutputData { ProcessLib::ConstitutiveRelations::StrainData eps_data; PermeabilityData permeability_data; - EnthalpyData enthalpy_data; + FluidEnthalpyData fluid_enthalpy_data; + SolidEnthalpyData solid_enthalpy_data; MassMoleFractionsData mass_mole_fractions_data; FluidDensityData fluid_density_data; VapourPartialPressureData vapour_pressure_data; @@ -114,7 +121,8 @@ struct OutputData return Reflection::reflectWithoutName(&Self::eps_data, &Self::permeability_data, - &Self::enthalpy_data, + &Self::fluid_enthalpy_data, + &Self::solid_enthalpy_data, &Self::mass_mole_fractions_data, &Self::fluid_density_data, &Self::vapour_pressure_data, @@ -141,73 +149,76 @@ struct ConstitutiveTempData ElasticTangentStiffnessData C_el_data; BiotData biot_data; SolidCompressibilityData beta_p_SR; - SaturationDataDeriv dS_L_dp_cap; BishopsData chi_S_L; SolidThermalExpansionData s_therm_exp_data; TotalStressData total_stress_data; EquivalentPlasticStrainData equivalent_plastic_strain_data; ViscosityData viscosity_data; PhaseTransitionData phase_transition_data; - PorosityDerivativeData porosity_d_data; - SolidDensityDerivativeData solid_density_d_data; SolidHeatCapacityData solid_heat_capacity_data; ThermalConductivityData thermal_conductivity_data; EffectiveVolumetricEnthalpy effective_volumetric_enthalpy_data; - EffectiveVolumetricEnthalpyDerivatives effective_volumetric_enthalpy_d_data; + AdvectionData advection_data; + VolumetricBodyForce volumetric_body_force; + FC1Data fC_1; + FC2aData fC_2a; + FC3aData fC_3a; + FC4LCpGData fC_4_LCpG; + FC4LCpCData fC_4_LCpC; + FC4LCTData fC_4_LCT; + FC4MCpGData fC_4_MCpG; + FC4MCpCData fC_4_MCpC; + FC4MCTData fC_4_MCT; + FC4MCuData fC_4_MCu; + + FW1Data fW_1; + FW2Data fW_2; + FW3aData fW_3a; + FW4LWpGData fW_4_LWpG; + FW4LWpCData fW_4_LWpC; + FW4LWTData fW_4_LWT; + FW4MWpGData fW_4_MWpG; + FW4MWpCData fW_4_MWpC; + FW4MWTData fW_4_MWT; + FW4MWuData fW_4_MWu; + + FT1Data fT_1; + FT2Data fT_2; + FT3Data fT_3; + + FU2KUpCData fu_2_KupC; +}; + +/// Data that stores intermediate values (derivatives), which are not needed +/// outside the constitutive setting. +template +struct DerivativesData +{ + SaturationDataDeriv dS_L_dp_cap; + AdvectionDerivativeData advection_d_data; + PorosityDerivativeData porosity_d_data; + ThermalConductivityDerivativeData + thermal_conductivity_d_data; + SolidDensityDerivativeData solid_density_d_data; EffectiveVolumetricInternalEnergyDerivatives effective_volumetric_internal_energy_d_data; - - using DisplacementDimVector = Eigen::Matrix; - using DisplacementDimMatrix = - Eigen::Matrix; - - DisplacementDimVector drho_GR_h_w_eff_dp_GR_Npart; - DisplacementDimMatrix drho_GR_h_w_eff_dp_GR_gradNpart; - DisplacementDimVector drho_LR_h_w_eff_dp_cap_Npart; - DisplacementDimMatrix drho_LR_h_w_eff_dp_cap_gradNpart; - DisplacementDimVector drho_GR_h_w_eff_dT; - DisplacementDimMatrix dfW_4_LWpG_a_dp_GR; - DisplacementDimMatrix dfW_4_LWpG_a_dp_cap; - DisplacementDimMatrix dfW_4_LWpG_a_dT; - DisplacementDimMatrix dfW_4_LWpG_d_dp_GR; - DisplacementDimMatrix dfW_4_LWpG_d_dp_cap; - DisplacementDimMatrix dfW_4_LWpG_d_dT; - DisplacementDimMatrix dfW_4_LWpC_a_dp_GR; - DisplacementDimMatrix dfW_4_LWpC_a_dp_cap; - DisplacementDimMatrix dfW_4_LWpC_a_dT; - DisplacementDimMatrix dfW_4_LWpC_d_dp_GR; - DisplacementDimMatrix dfW_4_LWpC_d_dp_cap; - DisplacementDimMatrix dfW_4_LWpC_d_dT; - DisplacementDimMatrix dfC_4_LCpG_dT; - DisplacementDimMatrix dfC_4_LCpC_a_dp_GR; - DisplacementDimMatrix dfC_4_LCpC_a_dp_cap; - DisplacementDimMatrix dfC_4_LCpC_a_dT; - DisplacementDimMatrix dfC_4_LCpC_d_dp_GR; - DisplacementDimMatrix dfC_4_LCpC_d_dp_cap; - DisplacementDimMatrix dfC_4_LCpC_d_dT; - DisplacementDimMatrix dadvection_C_dp_GR; - DisplacementDimMatrix dadvection_C_dp_cap; - DisplacementDimMatrix dk_over_mu_G_dp_cap; - DisplacementDimMatrix dk_over_mu_L_dp_cap; - double dfC_4_MCpG_dp_GR = std::numeric_limits::quiet_NaN(); - double dfC_4_MCpG_dT = std::numeric_limits::quiet_NaN(); - double dfC_4_MCT_dT = std::numeric_limits::quiet_NaN(); - double dfC_4_MCu_dT = std::numeric_limits::quiet_NaN(); - double dfC_3a_dp_GR = std::numeric_limits::quiet_NaN(); - double dfC_3a_dp_cap = std::numeric_limits::quiet_NaN(); - double dfC_3a_dT = std::numeric_limits::quiet_NaN(); - double dfC_2a_dp_GR = std::numeric_limits::quiet_NaN(); - double dfC_2a_dp_cap = std::numeric_limits::quiet_NaN(); - double dfC_2a_dT = std::numeric_limits::quiet_NaN(); - double dfW_2a_dp_GR = std::numeric_limits::quiet_NaN(); - double dfW_2b_dp_GR = std::numeric_limits::quiet_NaN(); - double dfW_2a_dp_cap = std::numeric_limits::quiet_NaN(); - double dfW_2b_dp_cap = std::numeric_limits::quiet_NaN(); - double dfW_2a_dT = std::numeric_limits::quiet_NaN(); - double dfW_2b_dT = std::numeric_limits::quiet_NaN(); - double dfW_3a_dp_GR = std::numeric_limits::quiet_NaN(); - double dfW_3a_dp_cap = std::numeric_limits::quiet_NaN(); - double dfW_3a_dT = std::numeric_limits::quiet_NaN(); + EffectiveVolumetricEnthalpyDerivatives effective_volumetric_enthalpy_d_data; + FC2aDerivativeData dfC_2a; + FC3aDerivativeData dfC_3a; + FC4LCpGDerivativeData dfC_4_LCpG; + FC4LCpCDerivativeData dfC_4_LCpC; + FC4MCpGDerivativeData dfC_4_MCpG; + FC4MCTDerivativeData dfC_4_MCT; + FC4MCuDerivativeData dfC_4_MCu; + FW2DerivativeData dfW_2; + FW3aDerivativeData dfW_3a; + FW4LWpGDerivativeData dfW_4_LWpG; + FW4LWpCDerivativeData dfW_4_LWpC; + FT1DerivativeData dfT_1; + FT2DerivativeData dfT_2; + FU1KUTDerivativeData dfu_1_KuT; + FU2KUpCDerivativeData dfu_2_KupC; }; + } // namespace ConstitutiveRelations } // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h index e37d7c906d2..277093d023f 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ConstitutiveModels.h @@ -9,9 +9,15 @@ #pragma once +#include "Advection.h" #include "Biot.h" #include "Bishops.h" +#include "CEquation.h" +#include "DiffusionVelocity.h" #include "ElasticTangentStiffnessModel.h" +#include "Enthalpy.h" +#include "Gravity.h" +#include "InternalEnergy.h" #include "MechanicalStrain.h" #include "PermeabilityModel.h" #include "PhaseTransitionModel.h" @@ -24,9 +30,12 @@ #include "SolidMechanics.h" #include "SolidThermalExpansion.h" #include "Swelling.h" +#include "TEquation.h" #include "ThermalConductivity.h" #include "TotalStress.h" +#include "UEquation.h" #include "Viscosity.h" +#include "WEquation.h" namespace ProcessLib::TH2M { @@ -74,6 +83,41 @@ struct ConstitutiveModels #endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION SolidHeatCapacityModel solid_heat_capacity_model; ThermalConductivityModel thermal_conductivity_model; + AdvectionModel advection_model; + InternalEnergyModel internal_energy_model; + EffectiveVolumetricEnthalpyModel effective_volumetric_enthalpy_model; + GravityModel gravity_model; + DiffusionVelocityModel diffusion_velocity_model; + SolidEnthalpyModel solid_enthalpy_model; + DarcyVelocityModel darcy_velocity_model; + FC1Model fC_1_model; + FC2aModel fC_2a_model; + FC3aModel fC_3a_model; + FC4LCpGModel fC_4_LCpG_model; + FC4LCpCModel fC_4_LCpC_model; + FC4LCTModel fC_4_LCT_model; + FC4MCpGModel fC_4_MCpG_model; + FC4MCpCModel fC_4_MCpC_model; + FC4MCTModel fC_4_MCT_model; + FC4MCuModel fC_4_MCu_model; + + FW1Model fW_1_model; + FW2Model fW_2_model; + FW3aModel fW_3a_model; + FW4LWpGModel fW_4_LWpG_model; + FW4LWpCModel fW_4_LWpC_model; + FW4LWTModel fW_4_LWT_model; + FW4MWpGModel fW_4_MWpG_model; + FW4MWpCModel fW_4_MWpC_model; + FW4MWTModel fW_4_MWT_model; + FW4MWuModel fW_4_MWu_model; + + FT1Model fT_1_model; + FT2Model fT_2_model; + FT3Model fT_3_model; + + FU1KUTModel fu_1_KuT_model; + FU2KUpCModel fu_2_KupC_model; }; } // namespace ConstitutiveRelations } // namespace ProcessLib::TH2M 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 +void DarcyVelocityModel::eval( + CapillaryPressureGradientData const& grad_p_cap, + FluidDensityData const& fluid_density_data, + GasPressureGradientData const& grad_p_GR, + PermeabilityData const& permeability_data, + SpecificBodyForceData const& specific_body_force, + ViscosityData const& viscosity_data, + DarcyVelocityData& 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 +struct DarcyVelocityModel +{ + void eval(CapillaryPressureGradientData const& grad_p_cap, + FluidDensityData const& fluid_density_data, + GasPressureGradientData const& grad_p_GR, + PermeabilityData const& permeability_data, + SpecificBodyForceData const& specific_body_force, + ViscosityData const& viscosity_data, + DarcyVelocityData& darcy_velocity_data) const; +}; + +extern template struct DarcyVelocityModel<2>; +extern template struct DarcyVelocityModel<3>; } // namespace ConstitutiveRelations } // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.cpp new file mode 100644 index 00000000000..20bd5917026 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.cpp @@ -0,0 +1,90 @@ +/** + * \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 "DiffusionVelocity.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +template +void DiffusionVelocityModel::eval( + CapillaryPressureGradientData const& grad_p_cap, + GasPressureGradientData const& grad_p_GR, + MassMoleFractionsData const& mass_mole_fractions_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + TemperatureGradientData const& grad_T, + DiffusionVelocityData& diffusion_velocity_data) const +{ + auto const gradxmWG = phase_transition_data.dxmWG_dpGR * grad_p_GR() + + phase_transition_data.dxmWG_dpCap * grad_p_cap() + + phase_transition_data.dxmWG_dT * grad_T(); + auto const gradxmCG = -gradxmWG; + + auto const gradxmWL = phase_transition_data.dxmWL_dpGR * grad_p_GR() + + phase_transition_data.dxmWL_dpCap * grad_p_cap() + + phase_transition_data.dxmWL_dT * grad_T(); + auto const gradxmCL = -gradxmWL; + + double const phi_L = S_L_data.S_L * porosity_data.phi; + double const phi_G = (1. - S_L_data.S_L) * porosity_data.phi; + + double const phi_G_D_vapour = + phi_G * phase_transition_data.diffusion_coefficient_vapour; + double const phi_L_D_solute = + phi_L * phase_transition_data.diffusion_coefficient_solute; + + if (mass_mole_fractions_data.xmCG == 0) + { + diffusion_velocity_data.d_CG = GlobalDimVector::Zero(); + } + else + { + diffusion_velocity_data.d_CG = + -phi_G_D_vapour / mass_mole_fractions_data.xmCG * gradxmCG; + } + + if (mass_mole_fractions_data.xmCG == 1) + { + diffusion_velocity_data.d_WG = GlobalDimVector::Zero(); + } + else + { + diffusion_velocity_data.d_WG = + -phi_G_D_vapour / (1 - mass_mole_fractions_data.xmCG) * gradxmWG; + } + + if (mass_mole_fractions_data.xmWL == 1) + + { + diffusion_velocity_data.d_CL = GlobalDimVector::Zero(); + } + else + { + diffusion_velocity_data.d_CL = + -phi_L_D_solute / (1. - mass_mole_fractions_data.xmWL) * gradxmCL; + } + + if (mass_mole_fractions_data.xmWL == 0) + { + diffusion_velocity_data.d_WL = GlobalDimVector::Zero(); + } + else + { + diffusion_velocity_data.d_WL = + -phi_L_D_solute / mass_mole_fractions_data.xmWL * gradxmWL; + } +} + +template struct DiffusionVelocityModel<2>; +template struct DiffusionVelocityModel<3>; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h index 25a203ca2c0..7ec45405ed8 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/DiffusionVelocity.h @@ -10,7 +10,11 @@ #pragma once #include "Base.h" +#include "MassMoleFractions.h" +#include "PhaseTransitionData.h" +#include "Porosity.h" #include "ProcessLib/Reflection/ReflectionData.h" +#include "Saturation.h" namespace ProcessLib::TH2M { @@ -38,5 +42,22 @@ struct DiffusionVelocityData &Self::d_WL)}; } }; + +template +struct DiffusionVelocityModel +{ + void eval( + CapillaryPressureGradientData const& grad_p_cap, + GasPressureGradientData const& grad_p_GR, + MassMoleFractionsData const& mass_mole_fractions_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + TemperatureGradientData const& grad_T, + DiffusionVelocityData& diffusion_velocity_data) const; +}; + +extern template struct DiffusionVelocityModel<2>; +extern template struct DiffusionVelocityModel<3>; } // namespace ConstitutiveRelations } // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.cpp new file mode 100644 index 00000000000..046bb5eb2fb --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.cpp @@ -0,0 +1,101 @@ +/** + * \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 "InternalEnergy.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +void EffectiveVolumetricEnthalpyModel::eval( + FluidDensityData const& fluid_density_data, + FluidEnthalpyData const& fluid_enthalpy_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + SolidDensityData const& solid_density_data, + SolidEnthalpyData const& solid_enthalpy_data, + EffectiveVolumetricEnthalpy& effective_volumetric_enthalpy_data) const +{ + auto const phi_L = S_L_data.S_L * porosity_data.phi; + auto const phi_G = (1. - S_L_data.S_L) * porosity_data.phi; + double const phi_S = 1. - porosity_data.phi; + + effective_volumetric_enthalpy_data.rho_h_eff = + phi_G * fluid_density_data.rho_GR * fluid_enthalpy_data.h_G + + phi_L * fluid_density_data.rho_LR * fluid_enthalpy_data.h_L + + phi_S * solid_density_data.rho_SR * solid_enthalpy_data.h_S; +} + +void EffectiveVolumetricEnthalpyModel::dEval( + FluidDensityData const& fluid_density_data, + FluidEnthalpyData const& fluid_enthalpy_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + PorosityDerivativeData const& porosity_d_data, + SaturationData const& S_L_data, + SolidDensityData const& solid_density_data, + SolidDensityDerivativeData const& solid_density_d_data, + SolidEnthalpyData const& solid_enthalpy_data, + SolidHeatCapacityData const& solid_heat_capacity_data, + EffectiveVolumetricEnthalpyDerivatives& + effective_volumetric_enthalpy_d_data) const +{ + auto const phi_L = S_L_data.S_L * porosity_data.phi; + auto const phi_G = (1. - S_L_data.S_L) * porosity_data.phi; + double const phi_S = 1. - porosity_data.phi; + + // From p_LR = p_GR - p_cap it follows for + // drho_LR/dp_GR = drho_LR/dp_LR * dp_LR/dp_GR + // = drho_LR/dp_LR * (dp_GR/dp_GR - dp_cap/dp_GR) + // = drho_LR/dp_LR * (1 - 0) + double const drho_LR_dp_GR = phase_transition_data.drho_LR_dp_LR; + double const drho_LR_dp_cap = -phase_transition_data.drho_LR_dp_LR; + // drho_GR_dp_cap = 0; + + effective_volumetric_enthalpy_d_data.drho_h_eff_dp_GR = + /*(dphi_G_dp_GR = 0) * fluid_density_data.rho_GR * + fluid_enthalpy_data.h_G +*/ + phi_G * phase_transition_data.drho_GR_dp_GR * fluid_enthalpy_data.h_G + + /*(dphi_L_dp_GR = 0) * fluid_density_data.rho_LR * + fluid_enthalpy_data.h_L +*/ + phi_L * drho_LR_dp_GR * fluid_enthalpy_data.h_L; + effective_volumetric_enthalpy_d_data.drho_h_eff_dp_cap = + porosity_d_data.dphi_L_dp_cap * fluid_density_data.rho_GR * + fluid_enthalpy_data.h_G + + /*phi_G * (drho_GR_dp_cap = 0) * fluid_enthalpy_data.h_G +*/ + porosity_d_data.dphi_L_dp_cap * fluid_density_data.rho_LR * + fluid_enthalpy_data.h_L + + phi_L * drho_LR_dp_cap * fluid_enthalpy_data.h_L; + + // TODO (naumov) Extend for temperature dependent porosities. + constexpr double dphi_G_dT = 0; + constexpr double dphi_L_dT = 0; + effective_volumetric_enthalpy_d_data.drho_h_eff_dT = + dphi_G_dT * fluid_density_data.rho_GR * fluid_enthalpy_data.h_G + + phi_G * phase_transition_data.drho_GR_dT * fluid_enthalpy_data.h_G + + phi_G * fluid_density_data.rho_GR * phase_transition_data.dh_G_dT + + dphi_L_dT * fluid_density_data.rho_LR * fluid_enthalpy_data.h_L + + phi_L * phase_transition_data.drho_LR_dT * fluid_enthalpy_data.h_L + + phi_L * fluid_density_data.rho_LR * phase_transition_data.dh_L_dT - + porosity_d_data.dphi_dT * solid_density_data.rho_SR * + solid_enthalpy_data.h_S + + phi_S * solid_density_d_data.drho_SR_dT * solid_enthalpy_data.h_S + + phi_S * solid_density_data.rho_SR * solid_heat_capacity_data(); +} + +void SolidEnthalpyModel::eval( + SolidHeatCapacityData const& solid_heat_capacity_data, + TemperatureData const& T_data, + SolidEnthalpyData& solid_enthalpy_data) const +{ + solid_enthalpy_data.h_S = solid_heat_capacity_data() * T_data.T; +} + +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h index b7a1686618b..4bff673ce9a 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/Enthalpy.h @@ -10,26 +10,44 @@ #pragma once #include "Base.h" +#include "Enthalpy.h" +#include "FluidDensity.h" +#include "PhaseTransitionData.h" +#include "Porosity.h" #include "ProcessLib/Reflection/ReflectionData.h" +#include "Saturation.h" +#include "SolidDensity.h" +#include "SolidHeatCapacity.h" namespace ProcessLib::TH2M { namespace ConstitutiveRelations { -struct EnthalpyData +struct FluidEnthalpyData { double h_G = nan; double h_L = nan; - double h_S = nan; static auto reflect() { - using Self = EnthalpyData; + using Self = FluidEnthalpyData; namespace R = ProcessLib::Reflection; return std::tuple{R::makeReflectionData("enthalpy_gas", &Self::h_G), - R::makeReflectionData("enthalpy_liquid", &Self::h_L), - R::makeReflectionData("enthalpy_solid", &Self::h_S)}; + R::makeReflectionData("enthalpy_liquid", &Self::h_L)}; + } +}; + +struct SolidEnthalpyData +{ + double h_S = nan; + + static auto reflect() + { + using Self = SolidEnthalpyData; + namespace R = ProcessLib::Reflection; + + return std::tuple{R::makeReflectionData("enthalpy_solid", &Self::h_S)}; } }; @@ -45,5 +63,36 @@ struct EffectiveVolumetricEnthalpyDerivatives double drho_h_eff_dp_cap = nan; }; +struct EffectiveVolumetricEnthalpyModel +{ + void eval( + FluidDensityData const& fluid_density_data, + FluidEnthalpyData const& fluid_enthalpy_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + SolidDensityData const& solid_density_data, + SolidEnthalpyData const& solid_enthalpy_data, + EffectiveVolumetricEnthalpy& effective_volumetric_enthalpy_data) const; + + void dEval(FluidDensityData const& fluid_density_data, + FluidEnthalpyData const& fluid_enthalpy_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + PorosityDerivativeData const& porosity_d_data, + SaturationData const& S_L_data, + SolidDensityData const& solid_density_data, + SolidDensityDerivativeData const& solid_density_d_data, + SolidEnthalpyData const& solid_enthalpy_data, + SolidHeatCapacityData const& solid_heat_capacity_data, + EffectiveVolumetricEnthalpyDerivatives& + effective_volumetric_enthalpy_d_data) const; +}; + +struct SolidEnthalpyModel +{ + void eval(SolidHeatCapacityData const& solid_heat_capacity_data, + TemperatureData const& T_data, + SolidEnthalpyData& solid_enthalpy_data) const; +}; } // namespace ConstitutiveRelations } // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Gravity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Gravity.cpp new file mode 100644 index 00000000000..8a4e994a876 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/Gravity.cpp @@ -0,0 +1,39 @@ +/** + * \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 "Gravity.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +template +void GravityModel::eval( + FluidDensityData const& fluid_density_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + SolidDensityData const& solid_density_data, + SpecificBodyForceData const& specific_body_force, + VolumetricBodyForce& volumetric_body_force) const +{ + auto const phi_G = (1. - S_L_data.S_L) * porosity_data.phi; + auto const phi_L = S_L_data.S_L * porosity_data.phi; + auto const phi_S = 1. - porosity_data.phi; + + auto const rhoGR = fluid_density_data.rho_GR; + auto const rhoLR = fluid_density_data.rho_LR; + auto const rho = + phi_G * rhoGR + phi_L * rhoLR + phi_S * solid_density_data.rho_SR; + *volumetric_body_force = rho * specific_body_force(); +} + +template struct GravityModel<2>; +template struct GravityModel<3>; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Gravity.h b/ProcessLib/TH2M/ConstitutiveRelations/Gravity.h new file mode 100644 index 00000000000..dc1aac245e6 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/Gravity.h @@ -0,0 +1,42 @@ +/** + * \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 + */ + +#pragma once + +#include "Base.h" +#include "FluidDensity.h" +#include "Porosity.h" +#include "Saturation.h" +#include "SolidDensity.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +template +using VolumetricBodyForce = + BaseLib::StrongType, + struct VolumetricBodyForceTag>; + +template +struct GravityModel +{ + void eval( + FluidDensityData const& fluid_density_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + SolidDensityData const& solid_density_data, + SpecificBodyForceData const& specific_body_force, + VolumetricBodyForce& volumetric_body_force) const; +}; + +extern template struct GravityModel<2>; +extern template struct GravityModel<3>; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.cpp b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.cpp new file mode 100644 index 00000000000..2aa2b1d7552 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.cpp @@ -0,0 +1,91 @@ +/** + * \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 "InternalEnergy.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +void InternalEnergyModel::eval(FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + SolidDensityData const& solid_density_data, + SolidEnthalpyData const& solid_enthalpy_data, + InternalEnergyData& internal_energy_data) const +{ + auto const phi_L = S_L_data.S_L * porosity_data.phi; + auto const phi_G = (1. - S_L_data.S_L) * porosity_data.phi; + double const phi_S = 1. - porosity_data.phi; + + auto const u_S = solid_enthalpy_data.h_S; + + internal_energy_data() = + phi_G * fluid_density_data.rho_GR * phase_transition_data.uG + + phi_L * fluid_density_data.rho_LR * phase_transition_data.uL + + phi_S * solid_density_data.rho_SR * u_S; +} + +void InternalEnergyModel::dEval( + FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + PorosityDerivativeData const& porosity_d_data, + SaturationData const& S_L_data, + SolidDensityData const& solid_density_data, + SolidDensityDerivativeData const& solid_density_d_data, + SolidEnthalpyData const& solid_enthalpy_data, + SolidHeatCapacityData const& solid_heat_capacity_data, + EffectiveVolumetricInternalEnergyDerivatives& + effective_volumetric_internal_energy_d_data) const +{ + auto const phi_L = S_L_data.S_L * porosity_data.phi; + auto const phi_G = (1. - S_L_data.S_L) * porosity_data.phi; + double const phi_S = 1. - porosity_data.phi; + + auto const u_S = solid_enthalpy_data.h_S; + effective_volumetric_internal_energy_d_data.drho_u_eff_dT = + phi_G * phase_transition_data.drho_GR_dT * phase_transition_data.uG + + phi_G * fluid_density_data.rho_GR * phase_transition_data.du_G_dT + + phi_L * phase_transition_data.drho_LR_dT * phase_transition_data.uL + + phi_L * fluid_density_data.rho_LR * phase_transition_data.du_L_dT + + phi_S * solid_density_d_data.drho_SR_dT * u_S + + phi_S * solid_density_data.rho_SR * solid_heat_capacity_data() - + porosity_d_data.dphi_dT * solid_density_data.rho_SR * u_S; + + // From p_LR = p_GR - p_cap it follows for + // drho_LR/dp_GR = drho_LR/dp_LR * dp_LR/dp_GR + // = drho_LR/dp_LR * (dp_GR/dp_GR - dp_cap/dp_GR) + // = drho_LR/dp_LR * (1 - 0) + double const drho_LR_dp_GR = phase_transition_data.drho_LR_dp_LR; + double const drho_LR_dp_cap = -phase_transition_data.drho_LR_dp_LR; + // drho_GR_dp_cap = 0; + + effective_volumetric_internal_energy_d_data.drho_u_eff_dp_GR = + /*(dphi_G_dp_GR = 0) * fluid_density_data.rho_GR * + phase_transition_data.uG +*/ + phi_G * phase_transition_data.drho_GR_dp_GR * phase_transition_data.uG + + phi_G * fluid_density_data.rho_GR * phase_transition_data.du_G_dp_GR + + /*(dphi_L_dp_GR = 0) * fluid_density_data.rho_LR * + phase_transition_data.uL +*/ + phi_L * drho_LR_dp_GR * phase_transition_data.uL + + phi_L * fluid_density_data.rho_LR * phase_transition_data.du_L_dp_GR; + + effective_volumetric_internal_energy_d_data.drho_u_eff_dp_cap = + -porosity_d_data.dphi_L_dp_cap * fluid_density_data.rho_GR * + phase_transition_data.uG + + /*phi_G * (drho_GR_dp_cap = 0) * phase_transition_data.uG +*/ + porosity_d_data.dphi_L_dp_cap * fluid_density_data.rho_LR * + phase_transition_data.uL + + phi_L * drho_LR_dp_cap * phase_transition_data.uL + + phi_L * fluid_density_data.rho_LR * phase_transition_data.du_L_dp_cap; +} +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h index 73644146b34..d74dfd71677 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/InternalEnergy.h @@ -11,6 +11,13 @@ #include "Base.h" #include "BaseLib/StrongType.h" +#include "Enthalpy.h" +#include "FluidDensity.h" +#include "PhaseTransitionData.h" +#include "Porosity.h" +#include "Saturation.h" +#include "SolidDensity.h" +#include "SolidHeatCapacity.h" namespace ProcessLib::TH2M { @@ -25,5 +32,28 @@ struct EffectiveVolumetricInternalEnergyDerivatives double drho_u_eff_dp_GR = nan; double drho_u_eff_dp_cap = nan; }; + +struct InternalEnergyModel +{ + void eval(FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + SolidDensityData const& solid_density_data, + SolidEnthalpyData const& solid_enthalpy_data, + InternalEnergyData& internal_energy_data) const; + + void dEval(FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + PorosityDerivativeData const& porosity_d_data, + SaturationData const& S_L_data, + SolidDensityData const& solid_density_data, + SolidDensityDerivativeData const& solid_density_d_data, + SolidEnthalpyData const& solid_enthalpy_data, + SolidHeatCapacityData const& solid_heat_capacity_data, + EffectiveVolumetricInternalEnergyDerivatives& + effective_volumetric_internal_energy_d_data) const; +}; } // namespace ConstitutiveRelations } // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp index 782b8b60ecb..72e475c0434 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.cpp @@ -43,7 +43,7 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, CapillaryPressureData const& p_cap, TemperatureData const& T_data, PureLiquidDensityData const& rho_W_LR, - EnthalpyData& enthalpy_data, + FluidEnthalpyData& fluid_enthalpy_data, MassMoleFractionsData& mass_mole_fractions_data, FluidDensityData& fluid_density_data, VapourPartialPressureData& vapour_pressure_data, @@ -108,14 +108,14 @@ void NoPhaseTransition::eval(SpaceTimeData const& x_t, .template value(variables, x_t.x, x_t.t, x_t.dt); // specific phase enthalpies - enthalpy_data.h_G = cpG * T; - enthalpy_data.h_L = cpL * T; + fluid_enthalpy_data.h_G = cpG * T; + fluid_enthalpy_data.h_L = cpL * T; cv.dh_G_dT = cpG; cv.dh_L_dT = cpL; // specific inner energies - cv.uG = enthalpy_data.h_G - pGR / fluid_density_data.rho_GR; - cv.uL = enthalpy_data.h_L; + cv.uG = fluid_enthalpy_data.h_G - pGR / fluid_density_data.rho_GR; + cv.uL = fluid_enthalpy_data.h_L; cv.drho_GR_dT = gas_phase[MaterialPropertyLib::PropertyType::density] diff --git a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.h b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.h index 918723e9060..102443321b9 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/NoPhaseTransition.h @@ -28,7 +28,7 @@ struct NoPhaseTransition : PhaseTransitionModel GasPressureData const& p_GR, CapillaryPressureData const& p_cap, TemperatureData const& T_data, PureLiquidDensityData const& rho_W_LR, - EnthalpyData& enthalpy_data, + FluidEnthalpyData& fluid_enthalpy_data, MassMoleFractionsData& mass_mole_fractions_data, FluidDensityData& fluid_density_data, VapourPartialPressureData& vapour_pressure_data, diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp index 7f71fcb2492..597bdca0440 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.cpp @@ -139,7 +139,7 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, CapillaryPressureData const& p_cap, TemperatureData const& T_data, PureLiquidDensityData const& rho_W_LR, - EnthalpyData& enthalpy_data, + FluidEnthalpyData& fluid_enthalpy_data, MassMoleFractionsData& mass_mole_fractions_data, FluidDensityData& fluid_density_data, VapourPartialPressureData& vapour_pressure_data, @@ -333,11 +333,12 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, cv.hWG = cpWG * T + dh_evap; // specific enthalpy of gas phase - enthalpy_data.h_G = mass_mole_fractions_data.xmCG * cv.hCG + xmWG * cv.hWG; + fluid_enthalpy_data.h_G = + mass_mole_fractions_data.xmCG * cv.hCG + xmWG * cv.hWG; cv.dh_G_dT = 0; // specific inner energies of gas phase - cv.uG = enthalpy_data.h_G - pGR / fluid_density_data.rho_GR; + cv.uG = fluid_enthalpy_data.h_G - pGR / fluid_density_data.rho_GR; cv.du_G_dT = 0; cv.du_G_dp_GR = 0; @@ -505,11 +506,11 @@ void PhaseTransition::eval(SpaceTimeData const& x_t, // specific enthalpy of liquid phase and its components double const hCL = cpCL * T + dh_sol; double const hWL = cpWL * T; - enthalpy_data.h_L = xmCL * hCL + mass_mole_fractions_data.xmWL * hWL; + fluid_enthalpy_data.h_L = xmCL * hCL + mass_mole_fractions_data.xmWL * hWL; cv.dh_L_dT = 0; // specific inner energies of liquid phase - cv.uL = enthalpy_data.h_L; + cv.uL = fluid_enthalpy_data.h_L; cv.du_L_dT = 0; // diffusion diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.h index 69d621b3dad..58766650f53 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransition.h @@ -28,7 +28,7 @@ struct PhaseTransition : PhaseTransitionModel GasPressureData const& p_GR, CapillaryPressureData const& p_cap, TemperatureData const& T_data, PureLiquidDensityData const& rho_W_LR, - EnthalpyData& enthalpy_data, + FluidEnthalpyData& fluid_enthalpy_data, MassMoleFractionsData& mass_mole_fractions_data, FluidDensityData& fluid_density_data, VapourPartialPressureData& vapour_pressure_data, diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h index a0042c6afbb..b6e27f07574 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionData.h @@ -51,10 +51,6 @@ struct PhaseTransitionData double dxmWL_dpGR = nan; double dxmWL_dpCap = nan; double dxmWL_dT = nan; - // TODO (naumov) This is zero in both models but used in the assembly. - // Remove it and simplify assembly or correct the expressions in the phase - // transition models? - static constexpr double dxmWL_dpLR = 0.; double diffusion_coefficient_vapour = nan; double diffusion_coefficient_solute = nan; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h index 8cf92e94f7e..0642ce42ca8 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/PhaseTransitionModel.h @@ -52,7 +52,7 @@ struct PhaseTransitionModel CapillaryPressureData const& p_cap, TemperatureData const& T_data, PureLiquidDensityData const& rho_W_LR, - EnthalpyData& enthalpy_data, + FluidEnthalpyData& fluid_enthalpy_data, MassMoleFractionsData& mass_mole_fractions_data, FluidDensityData& fluid_density_data, VapourPartialPressureData& vapour_pressure_data, diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp index 508dc689d3a..2688dd02698 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.cpp @@ -16,8 +16,7 @@ namespace ConstitutiveRelations void PorosityModel::eval(SpaceTimeData const& x_t, MediaData const& media_data, - PorosityData& porosity_data, - PorosityDerivativeData& porosity_d_data) const + PorosityData& porosity_data) const { MaterialPropertyLib::VariableArray variables; @@ -26,10 +25,23 @@ void PorosityModel::eval(SpaceTimeData const& x_t, porosity_data.phi = mpl_porosity.template value(variables, x_t.x, x_t.t, x_t.dt); +} + +void PorosityModel::dEval(SpaceTimeData const& x_t, MediaData const& media_data, + PorosityData const& porosity_data, + SaturationDataDeriv const& dS_L_dp_cap, + PorosityDerivativeData& porosity_d_data) const +{ + MaterialPropertyLib::VariableArray variables; + + auto const& mpl_porosity = + media_data.medium[MaterialPropertyLib::PropertyType::porosity]; porosity_d_data.dphi_dT = mpl_porosity.template dValue( variables, MaterialPropertyLib::Variable::temperature, x_t.x, x_t.t, x_t.dt); + + porosity_d_data.dphi_L_dp_cap = dS_L_dp_cap() * porosity_data.phi; } template @@ -39,8 +51,7 @@ void PorosityModelNonConstantSolidPhaseVolumeFraction::eval( BiotData const& biot, StrainData const& strain_data, SolidThermalExpansionData const& s_therm_exp_data, - PorosityData& porosity_data, - PorosityDerivativeData& porosity_d_data) const + PorosityData& porosity_data) const { MaterialPropertyLib::VariableArray variables; @@ -60,6 +71,31 @@ void PorosityModelNonConstantSolidPhaseVolumeFraction::eval( (1. + s_therm_exp_data.thermal_volume_strain - biot() * div_u); porosity_data.phi = 1. - phi_S; +} + +template +void PorosityModelNonConstantSolidPhaseVolumeFraction::dEval( + SpaceTimeData const& x_t, + MediaData const& media_data, + PorosityData const& porosity_data, + SaturationDataDeriv const& dS_L_dp_cap, + BiotData const& biot, + SolidThermalExpansionData const& s_therm_exp_data, + StrainData const& strain_data, + PorosityDerivativeData& porosity_d_data) const +{ + MaterialPropertyLib::VariableArray variables; + + auto const& mpl_porosity = + media_data.medium[MaterialPropertyLib::PropertyType::porosity]; + + double const phi_0 = + mpl_porosity.template value(variables, x_t.x, x_t.t, x_t.dt); + + static int const KelvinVectorSize = + MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim); + using Invariants = MathLib::KelvinVector::Invariants; + double const div_u = Invariants::trace(strain_data.eps); auto const dphi_0_dT = mpl_porosity.template dValue( variables, MaterialPropertyLib::Variable::temperature, x_t.x, x_t.t, @@ -69,6 +105,8 @@ void PorosityModelNonConstantSolidPhaseVolumeFraction::eval( dphi_0_dT * (1. + s_therm_exp_data.thermal_volume_strain - biot() * div_u) - (1. - phi_0) * s_therm_exp_data.beta_T_SR; + + porosity_d_data.dphi_L_dp_cap = dS_L_dp_cap() * porosity_data.phi; } template struct PorosityModelNonConstantSolidPhaseVolumeFraction<2>; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h index f7e0d5d8f09..49884ed45ed 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/Porosity.h @@ -13,6 +13,7 @@ #include "Biot.h" #include "ProcessLib/ConstitutiveRelations/StrainData.h" #include "ProcessLib/Reflection/ReflectionData.h" +#include "Saturation.h" #include "SolidThermalExpansion.h" namespace ProcessLib::TH2M @@ -22,6 +23,10 @@ namespace ConstitutiveRelations struct PorosityDerivativeData { double dphi_dT = nan; + // dphi_G_dp_GR = -ds_L_dp_GR * phi = 0; + // dphi_L_dp_GR = ds_L_dp_GR * phi = 0; + double dphi_L_dp_cap = nan; + // dphi_G_dp_cap = -dphi_L_dp_cap }; struct PorosityData @@ -41,8 +46,13 @@ struct PorosityModel { void eval(SpaceTimeData const& x_t, MediaData const& media_data, - PorosityData& porosity_data, - PorosityDerivativeData& porosity_d_data) const; + PorosityData& porosity_data) const; + + void dEval(SpaceTimeData const& x_t, + MediaData const& media_data, + PorosityData const& porosity_data, + SaturationDataDeriv const& dS_L_dp_cap, + PorosityDerivativeData& porosity_d_data) const; }; template @@ -54,7 +64,16 @@ struct PorosityModelNonConstantSolidPhaseVolumeFraction BiotData const& biot, StrainData const& strain_data, SolidThermalExpansionData const& s_therm_exp_data, - PorosityData& porosity_data, + PorosityData& porosity_data) const; + + void dEval( + SpaceTimeData const& x_t, + MediaData const& media_data, + PorosityData const& porosity_data, + SaturationDataDeriv const& dS_L_dp_cap, + BiotData const& biot, + SolidThermalExpansionData const& s_therm_exp_data, + StrainData const& strain_data, PorosityDerivativeData& porosity_d_data) const; }; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Saturation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/Saturation.cpp index 5f4693174a6..8ca73adaa1b 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/Saturation.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/Saturation.cpp @@ -16,8 +16,7 @@ namespace ConstitutiveRelations void SaturationModel::eval(SpaceTimeData const& x_t, MediaData const& media_data, CapillaryPressureData const& p_cap, - SaturationData& S_L_data, - SaturationDataDeriv& dS_L_data) const + SaturationData& S_L_data) const { namespace MPL = MaterialPropertyLib; MPL::VariableArray variables; @@ -27,6 +26,18 @@ void SaturationModel::eval(SpaceTimeData const& x_t, S_L_data.S_L = medium.property(MPL::PropertyType::saturation) .template value(variables, x_t.x, x_t.t, x_t.dt); +} + +void SaturationModel::dEval(SpaceTimeData const& x_t, + MediaData const& media_data, + CapillaryPressureData const& p_cap, + SaturationDataDeriv& dS_L_data) const +{ + namespace MPL = MaterialPropertyLib; + MPL::VariableArray variables; + variables.capillary_pressure = p_cap(); + + auto const& medium = media_data.medium; dS_L_data() = medium.property(MPL::PropertyType::saturation) .template dValue( diff --git a/ProcessLib/TH2M/ConstitutiveRelations/Saturation.h b/ProcessLib/TH2M/ConstitutiveRelations/Saturation.h index f8ab3071673..66dbbea5dcc 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/Saturation.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/Saturation.h @@ -34,8 +34,12 @@ using SaturationDataDeriv = struct SaturationModel { void eval(SpaceTimeData const& x_t, MediaData const& media_data, - CapillaryPressureData const& p_cap, SaturationData& S_L_data, - SaturationDataDeriv& dS_L_data) const; + CapillaryPressureData const& p_cap, + SaturationData& S_L_data) const; + + void dEval(SpaceTimeData const& x_t, MediaData const& media_data, + CapillaryPressureData const& p_cap, + SaturationDataDeriv& dS_L_data) const; }; } // namespace ConstitutiveRelations } // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.cpp index 9978c151690..5a764bc2015 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.cpp @@ -14,11 +14,25 @@ namespace ProcessLib::TH2M namespace ConstitutiveRelations { -void SolidDensityModel::eval( +void SolidDensityModel::eval(SpaceTimeData const& x_t, + MediaData const& media_data, + TemperatureData const& T_data, + SolidDensityData& solid_density_data) const +{ + MaterialPropertyLib::VariableArray variables; + variables.temperature = T_data.T; + + auto const& mpl_solid_density = + media_data.solid[MaterialPropertyLib::PropertyType::density]; + + solid_density_data.rho_SR = mpl_solid_density.template value( + variables, x_t.x, x_t.t, x_t.dt); +} + +void SolidDensityModel::dEval( SpaceTimeData const& x_t, MediaData const& media_data, TemperatureData const& T_data, - SolidDensityData& solid_density_data, SolidDensityDerivativeData& solid_density_d_data) const { MaterialPropertyLib::VariableArray variables; @@ -27,9 +41,6 @@ void SolidDensityModel::eval( auto const& mpl_solid_density = media_data.solid[MaterialPropertyLib::PropertyType::density]; - solid_density_data.rho_SR = mpl_solid_density.template value( - variables, x_t.x, x_t.t, x_t.dt); - solid_density_d_data.drho_SR_dT = mpl_solid_density.template dValue( variables, MaterialPropertyLib::Variable::temperature, x_t.x, x_t.t, x_t.dt); @@ -43,8 +54,7 @@ void SolidDensityModelNonConstantSolidPhaseVolumeFraction:: BiotData const& biot, StrainData const& strain_data, SolidThermalExpansionData const& s_therm_exp_data, - SolidDensityData& solid_density_data, - SolidDensityDerivativeData& solid_density_d_data) const + SolidDensityData& solid_density_data) const { MaterialPropertyLib::VariableArray variables; variables.temperature = T_data.T; @@ -63,6 +73,31 @@ void SolidDensityModelNonConstantSolidPhaseVolumeFraction:: solid_density_data.rho_SR = rho_ref_SR * (1. - s_therm_exp_data.thermal_volume_strain + (biot() - 1.) * div_u); +} + +template +void SolidDensityModelNonConstantSolidPhaseVolumeFraction:: + dEval(SpaceTimeData const& x_t, + MediaData const& media_data, + TemperatureData const& T_data, + BiotData const& biot, + StrainData const& strain_data, + SolidThermalExpansionData const& s_therm_exp_data, + SolidDensityDerivativeData& solid_density_d_data) const +{ + MaterialPropertyLib::VariableArray variables; + variables.temperature = T_data.T; + + static int const KelvinVectorSize = + MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim); + using Invariants = MathLib::KelvinVector::Invariants; + double const div_u = Invariants::trace(strain_data.eps); + + auto const& mpl_solid_density = + media_data.solid[MaterialPropertyLib::PropertyType::density]; + + auto const rho_ref_SR = mpl_solid_density.template value( + variables, x_t.x, x_t.t, x_t.dt); solid_density_d_data.drho_SR_dT = mpl_solid_density.template dValue( diff --git a/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h index ad14ab20ff8..d9f2d4890c4 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/SolidDensity.h @@ -43,8 +43,12 @@ struct SolidDensityModel void eval(SpaceTimeData const& x_t, MediaData const& media_data, TemperatureData const& T_data, - SolidDensityData& solid_density_data, - SolidDensityDerivativeData& solid_density_d_data) const; + SolidDensityData& solid_density_data) const; + + void dEval(SpaceTimeData const& x_t, + MediaData const& media_data, + TemperatureData const& T_data, + SolidDensityDerivativeData& solid_density_d_data) const; }; template @@ -57,7 +61,15 @@ struct SolidDensityModelNonConstantSolidPhaseVolumeFraction BiotData const& biot, StrainData const& strain_data, SolidThermalExpansionData const& s_therm_exp_data, - SolidDensityData& solid_density_data, + SolidDensityData& solid_density_data) const; + + void dEval( + SpaceTimeData const& x_t, + MediaData const& media_data, + TemperatureData const& T_data, + BiotData const& biot, + StrainData const& strain_data, + SolidThermalExpansionData const& s_therm_exp_data, SolidDensityDerivativeData& solid_density_d_data) const; }; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/TEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/TEquation.cpp new file mode 100644 index 00000000000..0dd5b7b814b --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/TEquation.cpp @@ -0,0 +1,147 @@ +/** + * \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 "TEquation.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +void FT1Model::eval( + double const dt, + InternalEnergyData const& internal_energy_data, + PrevState const& internal_energy_data_prev, + FT1Data& fT_1) const +{ + if (dt == 0) + { + fT_1.m = 0; + return; + } + + auto const rho_u_eff_dot = + (internal_energy_data() - **internal_energy_data_prev) / dt; + fT_1.m = rho_u_eff_dot; +} + +void FT1Model::dEval(double const dt, + EffectiveVolumetricInternalEnergyDerivatives const& + effective_volumetric_internal_energy_d_data, + FT1DerivativeData& dfT_1) const +{ + if (dt == 0) + { + dfT_1.dp_GR = 0; + dfT_1.dp_cap = 0; + dfT_1.dT = 0; + return; + } + + dfT_1.dp_GR = + effective_volumetric_internal_energy_d_data.drho_u_eff_dp_GR / dt; + + dfT_1.dp_cap = + effective_volumetric_internal_energy_d_data.drho_u_eff_dp_cap / dt; + + dfT_1.dT = effective_volumetric_internal_energy_d_data.drho_u_eff_dT / dt; +} + +template +void FT2Model::eval( + DarcyVelocityData const& darcy_velocity_data, + FluidDensityData const& fluid_density_data, + FluidEnthalpyData const& fluid_enthalpy_data, + FT2Data& fT_2) const +{ + fT_2.A.noalias() = fluid_density_data.rho_GR * fluid_enthalpy_data.h_G * + darcy_velocity_data.w_GS + + fluid_density_data.rho_LR * fluid_enthalpy_data.h_L * + darcy_velocity_data.w_LS; +} + +template +void FT2Model::dEval( + DarcyVelocityData const& darcy_velocity_data, + FluidDensityData const& fluid_density_data, + FluidEnthalpyData const& fluid_enthalpy_data, + PermeabilityData const& permeability_data, + PhaseTransitionData const& phase_transition_data, + SpecificBodyForceData const& specific_body_force, + ViscosityData const& viscosity_data, + FT2DerivativeData& dfT_2) 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; + + dfT_2.dp_GR_Npart = phase_transition_data.drho_GR_dp_GR * + fluid_enthalpy_data.h_G * darcy_velocity_data.w_GS + + fluid_density_data.rho_GR * fluid_enthalpy_data.h_G * + k_over_mu_G * phase_transition_data.drho_GR_dp_GR * + specific_body_force(); + dfT_2.dp_GR_gradNpart = + fluid_density_data.rho_GR * fluid_enthalpy_data.h_G * k_over_mu_G - + fluid_density_data.rho_LR * fluid_enthalpy_data.h_L * k_over_mu_L; + + // From p_LR = p_GR - p_cap it follows for + // drho_LR/dp_GR = drho_LR/dp_LR * dp_LR/dp_GR + // = drho_LR/dp_LR * (dp_GR/dp_GR - dp_cap/dp_GR) + // = drho_LR/dp_LR * (1 - 0) + double const drho_LR_dp_cap = -phase_transition_data.drho_LR_dp_LR; + + dfT_2.dp_cap_Npart = + -drho_LR_dp_cap * fluid_enthalpy_data.h_L * darcy_velocity_data.w_LS - + fluid_density_data.rho_LR * fluid_enthalpy_data.h_L * k_over_mu_L * + drho_LR_dp_cap * specific_body_force(); + dfT_2.dp_cap_gradNpart = + fluid_density_data.rho_LR * fluid_enthalpy_data.h_L * k_over_mu_L; + + dfT_2.dT = phase_transition_data.drho_GR_dT * fluid_enthalpy_data.h_G * + darcy_velocity_data.w_GS + + fluid_density_data.rho_GR * phase_transition_data.dh_G_dT * + darcy_velocity_data.w_GS + + phase_transition_data.drho_LR_dT * fluid_enthalpy_data.h_L * + darcy_velocity_data.w_LS + + fluid_density_data.rho_LR * phase_transition_data.dh_L_dT * + darcy_velocity_data.w_LS; + // TODO (naumov) + k_over_mu_G * drho_GR_dT * b + k_over_mu_L * + // drho_LR_dT * b +} + +template struct FT2Model<2>; +template struct FT2Model<3>; + +template +void FT3Model::eval( + ConstituentDensityData const& constituent_density_data, + DarcyVelocityData const& darcy_velocity_data, + DiffusionVelocityData const& diffusion_velocity_data, + FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + SpecificBodyForceData const& specific_body_force, + FT3Data& fT_3) const +{ + fT_3.N = + (fluid_density_data.rho_GR * darcy_velocity_data.w_GS.transpose() + + fluid_density_data.rho_LR * darcy_velocity_data.w_LS.transpose()) * + specific_body_force(); + + fT_3.gradN.noalias() = + constituent_density_data.rho_C_GR * phase_transition_data.hCG * + diffusion_velocity_data.d_CG + + constituent_density_data.rho_W_GR * phase_transition_data.hWG * + diffusion_velocity_data.d_WG; +} + +template struct FT3Model<2>; +template struct FT3Model<3>; + +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/TEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/TEquation.h new file mode 100644 index 00000000000..94d80aaf05f --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/TEquation.h @@ -0,0 +1,112 @@ +/** + * \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 + */ + +#pragma once + +#include "Base.h" +#include "ConstitutiveDensity.h" +#include "DarcyVelocity.h" +#include "DiffusionVelocity.h" +#include "Enthalpy.h" +#include "FluidDensity.h" +#include "InternalEnergy.h" +#include "PermeabilityData.h" +#include "Viscosity.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +struct FT1Data +{ + double m = nan; +}; + +struct FT1DerivativeData +{ + double dp_GR = nan; + double dp_cap = nan; + double dT = nan; +}; + +struct FT1Model +{ + void eval(double const dt, + InternalEnergyData const& internal_energy_data, + PrevState const& internal_energy_data_prev, + FT1Data& fT_1) const; + + void dEval(double const dt, + EffectiveVolumetricInternalEnergyDerivatives const& + effective_volumetric_internal_energy_d_data, + FT1DerivativeData& dfT_1) const; +}; + +template +struct FT2Data +{ + GlobalDimVector A; +}; + +template +struct FT2DerivativeData +{ + GlobalDimVector dp_GR_Npart; + GlobalDimMatrix dp_GR_gradNpart; + GlobalDimVector dp_cap_Npart; + GlobalDimMatrix dp_cap_gradNpart; + GlobalDimVector dT; +}; + +template +struct FT2Model +{ + void eval(DarcyVelocityData const& darcy_velocity_data, + FluidDensityData const& fluid_density_data, + FluidEnthalpyData const& fluid_enthalpy_data, + FT2Data& fT_2) const; + + void dEval( + DarcyVelocityData const& darcy_velocity_data, + FluidDensityData const& fluid_density_data, + FluidEnthalpyData const& fluid_enthalpy_data, + PermeabilityData const& permeability_data, + PhaseTransitionData const& phase_transition_data, + SpecificBodyForceData const& specific_body_force, + ViscosityData const& viscosity_data, + FT2DerivativeData& dfT_2) const; +}; + +extern template struct FT2Model<2>; +extern template struct FT2Model<3>; + +template +struct FT3Data +{ + double N = nan; + GlobalDimVector gradN; +}; + +template +struct FT3Model +{ + void eval( + ConstituentDensityData const& constituent_density_data, + DarcyVelocityData const& darcy_velocity_data, + DiffusionVelocityData const& diffusion_velocity_data, + FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + SpecificBodyForceData const& specific_body_force, + FT3Data& fT_3) const; +}; + +extern template struct FT3Model<2>; +extern template struct FT3Model<3>; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp index abd9e0cc9f3..358ee797cc7 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp +++ b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.cpp @@ -19,8 +19,7 @@ template void ThermalConductivityModel::eval( SpaceTimeData const& x_t, MediaData const& media_data, TemperatureData const& T_data, PorosityData const& porosity_data, - PorosityDerivativeData const& porosity_d_data, - SaturationData const& S_L_data, SaturationDataDeriv const& dS_L_dp_cap, + SaturationData const& S_L_data, ThermalConductivityData& thermal_conductivity_data) const { namespace MPL = MaterialPropertyLib; @@ -34,6 +33,22 @@ void ThermalConductivityModel::eval( thermal_conductivity_data.lambda = MPL::formEigenTensor( mpl_thermal_conductivity.value(variables, x_t.x, x_t.t, x_t.dt)); +} + +template +void ThermalConductivityModel::dEval( + SpaceTimeData const& x_t, MediaData const& media_data, + TemperatureData const& T_data, PorosityData const& porosity_data, + PorosityDerivativeData const& porosity_d_data, + SaturationData const& S_L_data, + ThermalConductivityDerivativeData& + thermal_conductivity_d_data) const +{ + namespace MPL = MaterialPropertyLib; + MPL::VariableArray variables; + variables.temperature = T_data.T; + variables.porosity = porosity_data.phi; + variables.liquid_saturation = S_L_data.S_L; // Derivatives computed here and not in the MPL property because various // derivatives are not available in the VariableArray. @@ -83,13 +98,9 @@ void ThermalConductivityModel::eval( x_t.t, x_t.dt)) : MPL::formEigenTensor(0.); - // dphi_G_dp_GR = -ds_L_dp_GR * phi = 0; - double const dphi_G_dp_cap = -dS_L_dp_cap() * porosity_data.phi; - // dphi_L_dp_GR = ds_L_dp_GR * phi = 0; - double const dphi_L_dp_cap = dS_L_dp_cap() * porosity_data.phi; - - thermal_conductivity_data.dlambda_dp_cap = - dphi_G_dp_cap * lambdaGR + dphi_L_dp_cap * lambdaLR; + thermal_conductivity_d_data.dlambda_dp_cap = + -porosity_d_data.dphi_L_dp_cap * lambdaGR + + porosity_d_data.dphi_L_dp_cap * lambdaLR; double const phi_L = S_L_data.S_L * porosity_data.phi; double const phi_G = (1. - S_L_data.S_L) * porosity_data.phi; @@ -99,13 +110,14 @@ void ThermalConductivityModel::eval( // dphi_G_dT = -dS_L/dT * phi + (1 - S_L) * dphi_dT = (1 - S_L) * dphi_dT // dphi_L_dT = dS_L/dT * phi + S_L * dphi_dT S_L * dphi_dT // dphi_S_dT = -dphi_dT -dphi_dT - thermal_conductivity_data.dlambda_dT = + thermal_conductivity_d_data.dlambda_dT = (1 - S_L_data.S_L) * porosity_d_data.dphi_dT * lambdaGR + phi_G * dlambda_GR_dT + S_L_data.S_L * porosity_d_data.dphi_dT * lambdaLR + +phi_L * dlambda_LR_dT - porosity_d_data.dphi_dT * lambdaSR + phi_S * dlambda_SR_dT; } + template struct ThermalConductivityModel<2>; template struct ThermalConductivityModel<3>; } // namespace ConstitutiveRelations diff --git a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h index fbf34021e55..35720adac3f 100644 --- a/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h +++ b/ProcessLib/TH2M/ConstitutiveRelations/ThermalConductivity.h @@ -22,6 +22,11 @@ template struct ThermalConductivityData { GlobalDimMatrix lambda; +}; + +template +struct ThermalConductivityDerivativeData +{ // Currently unused, but there is a comment in TH2MFEM-impl.h referring to // this matrix // GlobalDimMatrix dlambda_dp_GR; @@ -34,11 +39,16 @@ struct ThermalConductivityModel { void eval(SpaceTimeData const& x_t, MediaData const& media_data, TemperatureData const& T_data, PorosityData const& porosity_data, - PorosityDerivativeData const& porosity_d_data, SaturationData const& S_L_data, - SaturationDataDeriv const& dS_L_dp_cap, ThermalConductivityData& thermal_conductivity_data) const; + + void dEval(SpaceTimeData const& x_t, MediaData const& media_data, + TemperatureData const& T_data, PorosityData const& porosity_data, + PorosityDerivativeData const& porosity_d_data, + SaturationData const& S_L_data, + ThermalConductivityDerivativeData& + thermal_conductivity_d_data) const; }; extern template struct ThermalConductivityModel<2>; diff --git a/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp new file mode 100644 index 00000000000..dcea844e0a4 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.cpp @@ -0,0 +1,47 @@ +/** + * \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 "UEquation.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +template +void FU1KUTModel::dEval( + SolidMechanicsDataStateless const& s_mech_data, + SolidThermalExpansionData const& s_therm_exp_data, + FU1KUTDerivativeData& dfu_1_KuT) const +{ + dfu_1_KuT.dT = s_mech_data.stiffness_tensor * + s_therm_exp_data.solid_linear_thermal_expansivity_vector; +} + +template struct FU1KUTModel<2>; +template struct FU1KUTModel<3>; + +void FU2KUpCModel::eval(BiotData const& biot_data, + BishopsData const& chi_S_L, + FU2KUpCData& fu_2_KupC) const +{ + fu_2_KupC.m = biot_data() * chi_S_L.chi_S_L; +} + +void FU2KUpCModel::dEval(BiotData const& biot_data, + BishopsData const& chi_S_L, + CapillaryPressureData const& p_cap, + SaturationDataDeriv const& dS_L_dp_cap, + FU2KUpCDerivativeData& dfu_2_KupC) const +{ + dfu_2_KupC.dp_cap = + biot_data() * chi_S_L.dchi_dS_L * dS_L_dp_cap() * p_cap(); +} + +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h new file mode 100644 index 00000000000..70388bf50d2 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/UEquation.h @@ -0,0 +1,64 @@ +/** + * \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 + */ + +#pragma once + +#include "Base.h" +#include "Biot.h" +#include "Bishops.h" +#include "Saturation.h" +#include "SolidMechanics.h" +#include "SolidThermalExpansion.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +template +struct FU1KUTDerivativeData +{ + KelvinVector dT; +}; + +template +struct FU1KUTModel +{ + void dEval( + SolidMechanicsDataStateless const& s_mech_data, + SolidThermalExpansionData const& s_therm_exp_data, + FU1KUTDerivativeData& dfu_1_KuT) const; +}; + +extern template struct FU1KUTModel<2>; +extern template struct FU1KUTModel<3>; + +struct FU2KUpCData +{ + double m = nan; +}; + +struct FU2KUpCDerivativeData +{ + double dp_cap = nan; +}; + +struct FU2KUpCModel +{ + void eval(BiotData const& biot_data, + BishopsData const& chi_S_L, + FU2KUpCData& fu_2_KupC) const; + + void dEval(BiotData const& biot_data, + BishopsData const& chi_S_L, + CapillaryPressureData const& p_cap, + SaturationDataDeriv const& dS_L_dp_cap, + FU2KUpCDerivativeData& dfu_2_KupC) const; +}; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp new file mode 100644 index 00000000000..370e22cd5d0 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.cpp @@ -0,0 +1,449 @@ +/** + * \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 "WEquation.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +template +void FW1Model::eval( + AdvectionData const& advection_data, + FluidDensityData const& fluid_density_data, + FW1Data& fW_1) const +{ + fW_1.A = advection_data.advection_W_G * fluid_density_data.rho_GR + + advection_data.advection_W_L * fluid_density_data.rho_LR; +} + +template struct FW1Model<2>; +template struct FW1Model<3>; + +void FW2Model::eval(BiotData const biot_data, + CapillaryPressureData const pCap, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + SolidCompressibilityData const beta_p_SR, + FW2Data& fW_2) const +{ + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + double const rho_W_FR = + S_G * constituent_density_data.rho_W_GR + S_L * rho_W_LR(); + + fW_2.a = + porosity_data.phi * (rho_W_LR() - constituent_density_data.rho_W_GR) - + rho_W_FR * pCap() * (biot_data() - porosity_data.phi) * beta_p_SR(); +} + +void FW2Model::dEval(BiotData const& biot_data, + CapillaryPressureData const pCap, + ConstituentDensityData const& constituent_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + PorosityDerivativeData const& porosity_d_data, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + SaturationDataDeriv const& dS_L_dp_cap, + SolidCompressibilityData const& beta_p_SR, + FW2DerivativeData& dfW_2) const +{ + double const S_L = S_L_data.S_L; + double const S_G = 1. - S_L; + + double const drho_C_FR_dp_GR = + /*(dS_G_dp_GR = 0) * constituent_density_data.rho_C_GR +*/ + S_G * phase_transition_data.drho_C_GR_dp_GR + + /*(dS_L_dp_GR = 0) * constituent_density_data.rho_C_LR +*/ + S_L * phase_transition_data.drho_C_LR_dp_GR; + + dfW_2.dp_GR = -porosity_data.phi * phase_transition_data.drho_C_GR_dp_GR - + drho_C_FR_dp_GR * pCap() * (biot_data() - porosity_data.phi) * + beta_p_SR(); + + double const dfW_2a_dp_cap = + porosity_data.phi * (-phase_transition_data.drho_W_LR_dp_LR - + phase_transition_data.drho_W_GR_dp_cap); + double const rho_W_FR = + S_G * constituent_density_data.rho_W_GR + S_L * rho_W_LR(); + + double const drho_W_FR_dp_cap = + -dS_L_dp_cap() * constituent_density_data.rho_W_GR + + S_G * phase_transition_data.drho_W_GR_dp_cap + + dS_L_dp_cap() * rho_W_LR() - + S_L * phase_transition_data.drho_W_LR_dp_LR; + + double const dfW_2b_dp_cap = + drho_W_FR_dp_cap * pCap() * (biot_data() - porosity_data.phi) * + beta_p_SR() + + rho_W_FR * (biot_data() - porosity_data.phi) * beta_p_SR(); + + dfW_2.dp_cap = dfW_2a_dp_cap - dfW_2b_dp_cap; + + double const drho_W_FR_dT = S_G * phase_transition_data.drho_W_GR_dT + + S_L * phase_transition_data.drho_W_LR_dT; + + double const dfW_2a_dT = + porosity_d_data.dphi_dT * + (rho_W_LR() - constituent_density_data.rho_W_GR) + + porosity_data.phi * (phase_transition_data.drho_W_LR_dT - + phase_transition_data.drho_W_GR_dT); + double const dfW_2b_dT = + drho_W_FR_dT * pCap() * (biot_data() - porosity_data.phi) * + beta_p_SR() - + rho_W_FR * pCap() * porosity_d_data.dphi_dT * beta_p_SR(); + + dfW_2.dT = dfW_2a_dT - dfW_2b_dT; +} + +void FW3aModel::eval( + double const dt, + ConstituentDensityData const& constituent_density_data, + PrevState const& constituent_density_data_prev, + PrevState const& rho_W_LR_prev, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + FW3aData& fW_3a) const +{ + if (dt == 0.) + { + fW_3a.a = 0; + return; + } + + double const rho_W_GR_dot = (constituent_density_data.rho_W_GR - + constituent_density_data_prev->rho_W_GR) / + dt; + double const rho_W_LR_dot = (rho_W_LR() - **rho_W_LR_prev) / dt; + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + fW_3a.a = S_G * rho_W_GR_dot + S_L * rho_W_LR_dot; +} + +void FW3aModel::dEval( + double const dt, + ConstituentDensityData const& constituent_density_data, + PhaseTransitionData const& phase_transition_data, + PrevState const& constituent_density_data_prev, + PrevState const& rho_W_LR_prev, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + SaturationDataDeriv const& dS_L_dp_cap, + FW3aDerivativeData& dfW_3a) const +{ + if (dt == 0.) + { + dfW_3a.dp_GR = 0.; + dfW_3a.dp_cap = 0.; + dfW_3a.dT = 0.; + return; + } + + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + + double const rho_W_GR_dot = (constituent_density_data.rho_W_GR - + constituent_density_data_prev->rho_W_GR) / + dt; + double const rho_W_LR_dot = (rho_W_LR() - **rho_W_LR_prev) / dt; + + dfW_3a.dp_GR = /*(ds_G_dp_GR = 0) * rho_W_GR_dot +*/ + S_G * phase_transition_data.drho_W_GR_dp_GR / dt + + /*(ds_L_dp_GR = 0) * rho_W_LR_dot +*/ + S_L * phase_transition_data.drho_W_LR_dp_GR / dt; + + dfW_3a.dp_cap = -dS_L_dp_cap() * rho_W_GR_dot + + S_G * phase_transition_data.drho_W_GR_dp_cap / dt + + dS_L_dp_cap() * rho_W_LR_dot - + S_L * phase_transition_data.drho_W_LR_dp_LR / dt; + dfW_3a.dT = S_G * phase_transition_data.drho_W_GR_dT / dt + + S_L * phase_transition_data.drho_W_LR_dT / dt; +} + +template +void FW4LWpGModel::eval( + AdvectionData const& advection_data, + FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + FW4LWpGData& fW_4_LWpG) const +{ + GlobalDimMatrix const advection_W = + advection_data.advection_W_G + advection_data.advection_W_L; + + double const sD_G = phase_transition_data.diffusion_coefficient_vapour; + double const sD_L = phase_transition_data.diffusion_coefficient_solute; + + double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi; + double const phi_L = S_L_data.S_L * porosity_data.phi; + + double const diffusion_WGpGR = phi_G * fluid_density_data.rho_GR * sD_G * + phase_transition_data.dxmWG_dpGR; + double const diffusion_WLpGR = phi_L * fluid_density_data.rho_LR * sD_L * + phase_transition_data.dxmWL_dpGR; + double const diffusion_W_pGR = diffusion_WGpGR + diffusion_WLpGR; + + auto const I = + Eigen::Matrix::Identity(); + fW_4_LWpG.L.noalias() = diffusion_W_pGR * I + advection_W; +} + +template +void FW4LWpGModel::dEval( + ConstituentDensityData const& constituent_density_data, + PermeabilityData const& permeability_data, + PhaseTransitionData const& phase_transition_data, + PureLiquidDensityData const& rho_W_LR, + SaturationDataDeriv const& dS_L_dp_cap, + ViscosityData const& viscosity_data, + FW4LWpGDerivativeData& dfW_4_LWpG) const +{ + ////// Diffusion Part ///// + // TODO (naumov) d(diffusion_W_p)/dX for dxmW*/d* != 0 + + 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; + + // dk_over_mu_G_dp_GR = ip_out.permeability_data.Ki * + // ip_out.permeability_data.dk_rel_G_dS_L * + // (ds_L_dp_GR = 0) / + // ip_cv.viscosity_data.mu_GR = 0; + // dk_over_mu_L_dp_GR = ip_out.permeability_data.Ki * + // ip_out.permeability_data.dk_rel_L_dS_L * + // (ds_L_dp_GR = 0) / + // ip_cv.viscosity_data.mu_LR = 0; + auto const dk_over_mu_G_dp_cap = permeability_data.Ki * + permeability_data.dk_rel_G_dS_L * + dS_L_dp_cap() / viscosity_data.mu_GR; + + auto const dk_over_mu_L_dp_cap = permeability_data.Ki * + permeability_data.dk_rel_L_dS_L * + dS_L_dp_cap() / viscosity_data.mu_LR; + + dfW_4_LWpG.dp_GR = phase_transition_data.drho_W_GR_dp_GR * k_over_mu_G + // + rhoWGR * (dk_over_mu_G_dp_GR = 0) + + phase_transition_data.drho_W_LR_dp_GR * k_over_mu_L + // + rhoWLR * (dk_over_mu_L_dp_GR = 0) + ; + + dfW_4_LWpG.dp_cap = + phase_transition_data.drho_W_GR_dp_cap * k_over_mu_G + + constituent_density_data.rho_W_GR * dk_over_mu_G_dp_cap + + -phase_transition_data.drho_W_LR_dp_LR * k_over_mu_L + + rho_W_LR() * dk_over_mu_L_dp_cap; + + dfW_4_LWpG.dT = phase_transition_data.drho_W_GR_dT * k_over_mu_G + //+ rhoWGR * (dk_over_mu_G_dT != 0 TODO for mu_G(T)) + + phase_transition_data.drho_W_LR_dT * k_over_mu_L + //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_G(T)) + ; +} + +template struct FW4LWpGModel<2>; +template struct FW4LWpGModel<3>; + +template +void FW4LWpCModel::eval( + AdvectionData const& advection_data, + FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + FW4LWpCData& fW_4_LWpC) const +{ + double const sD_G = phase_transition_data.diffusion_coefficient_vapour; + double const sD_L = phase_transition_data.diffusion_coefficient_solute; + + double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi; + double const phi_L = S_L_data.S_L * porosity_data.phi; + + double const diffusion_WGpCap = phi_G * fluid_density_data.rho_GR * sD_G * + phase_transition_data.dxmWG_dpCap; + double const diffusion_WLpCap = phi_L * fluid_density_data.rho_LR * sD_L * + phase_transition_data.dxmWL_dpCap; + + double const diffusion_W_pCap = diffusion_WGpCap + diffusion_WLpCap; + + auto const I = + Eigen::Matrix::Identity(); + + fW_4_LWpC.L.noalias() = diffusion_W_pCap * I - advection_data.advection_W_L; +} + +template +void FW4LWpCModel::dEval( + AdvectionData const& advection_data, + FluidDensityData const& fluid_density_data, + PermeabilityData const& permeability_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + SaturationDataDeriv const& dS_L_dp_cap, + ViscosityData const& viscosity_data, + FW4LWpCDerivativeData& dfW_4_LWpC) const +{ + ////// Diffusion Part ///// + // TODO (naumov) d(diffusion_W_pCap)/dX for dxmW*/d* != 0 + + ////// Advection part ///// + GlobalDimMatrix const k_over_mu_L = + permeability_data.Ki * permeability_data.k_rel_L / viscosity_data.mu_LR; + + dfW_4_LWpC.dp_GR = phase_transition_data.drho_W_LR_dp_GR * k_over_mu_L + //+ rhoWLR * (dk_over_mu_L_dp_GR = 0) + ; + + double const sD_G = phase_transition_data.diffusion_coefficient_vapour; + double const sD_L = phase_transition_data.diffusion_coefficient_solute; + + double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi; + double const phi_L = S_L_data.S_L * porosity_data.phi; + + double const diffusion_WGpCap = phi_G * fluid_density_data.rho_GR * sD_G * + phase_transition_data.dxmWG_dpCap; + double const diffusion_WLpCap = phi_L * fluid_density_data.rho_LR * sD_L * + phase_transition_data.dxmWL_dpCap; + + double const diffusion_W_pCap = diffusion_WGpCap + diffusion_WLpCap; + + auto const I = + Eigen::Matrix::Identity(); + + dfW_4_LWpC.dp_cap = diffusion_W_pCap * I - advection_data.advection_W_L; + + auto const dk_over_mu_L_dp_cap = permeability_data.Ki * + permeability_data.dk_rel_L_dS_L * + dS_L_dp_cap() / viscosity_data.mu_LR; + dfW_4_LWpC.dp_cap = -phase_transition_data.drho_W_LR_dp_LR * k_over_mu_L + + rho_W_LR() * dk_over_mu_L_dp_cap; + + dfW_4_LWpC.dT = phase_transition_data.drho_W_LR_dT * k_over_mu_L + //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T)) + ; +} + +template struct FW4LWpCModel<2>; +template struct FW4LWpCModel<3>; + +template +void FW4LWTModel::eval( + FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + FW4LWTData& fW_4_LWT) const +{ + double const sD_G = phase_transition_data.diffusion_coefficient_vapour; + double const sD_L = phase_transition_data.diffusion_coefficient_solute; + + double const phi_G = (1 - S_L_data.S_L) * porosity_data.phi; + double const phi_L = S_L_data.S_L * porosity_data.phi; + + double const diffusion_W_G_T = phi_G * fluid_density_data.rho_GR * sD_G * + phase_transition_data.dxmWG_dT; + double const diffusion_W_L_T = phi_L * fluid_density_data.rho_LR * sD_L * + phase_transition_data.dxmWL_dT; + + double const diffusion_W_T = diffusion_W_G_T + diffusion_W_L_T; + + auto const I = + Eigen::Matrix::Identity(); + + fW_4_LWT.L.noalias() = diffusion_W_T * I; +} + +template struct FW4LWTModel<2>; +template struct FW4LWTModel<3>; + +void FW4MWpGModel::eval(BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + SolidCompressibilityData const& beta_p_SR, + FW4MWpGData& fW_4_MWpG) const +{ + double const S_L = S_L_data.S_L; + double const S_G = 1. - S_L; + + double const rho_W_FR = + S_G * constituent_density_data.rho_W_GR + S_L * rho_W_LR(); + + fW_4_MWpG.m = rho_W_FR * (biot_data() - porosity_data.phi) * beta_p_SR(); +} + +void FW4MWpCModel::eval(BiotData const& biot_data, + CapillaryPressureData const pCap, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + PrevState const& S_L_data_prev, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + SolidCompressibilityData const& beta_p_SR, + FW4MWpCData& fW_4_MWpC) const +{ + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + double const rho_W_FR = + S_G * constituent_density_data.rho_W_GR + S_L * rho_W_LR(); + + fW_4_MWpC.m = + -rho_W_FR * (biot_data() - porosity_data.phi) * beta_p_SR() * S_L; + + fW_4_MWpC.ml = + (porosity_data.phi * (rho_W_LR() - constituent_density_data.rho_W_GR) - + rho_W_FR * pCap() * (biot_data() - porosity_data.phi) * beta_p_SR()) * + (S_L - S_L_data_prev->S_L); +} + +template +void FW4MWTModel::eval( + BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + SolidThermalExpansionData const& s_therm_exp_data, + FW4MWTData& fW_4_MWT) const +{ + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + double const rho_W_FR = + S_G * constituent_density_data.rho_W_GR + S_L * rho_W_LR(); + + fW_4_MWT.m = -rho_W_FR * (biot_data() - porosity_data.phi) * + s_therm_exp_data.beta_T_SR; +} + +template struct FW4MWTModel<2>; +template struct FW4MWTModel<3>; + +void FW4MWuModel::eval(BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + FW4MWuData& fW_4_MWu) const +{ + auto const S_L = S_L_data.S_L; + auto const S_G = 1. - S_L; + double const rho_W_FR = + S_G * constituent_density_data.rho_W_GR + S_L * rho_W_LR(); + + fW_4_MWu.m = rho_W_FR * biot_data(); +} + +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h new file mode 100644 index 00000000000..09f1b09fc44 --- /dev/null +++ b/ProcessLib/TH2M/ConstitutiveRelations/WEquation.h @@ -0,0 +1,277 @@ +/** + * \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 + */ + +#pragma once + +#include "Advection.h" +#include "Base.h" +#include "Biot.h" +#include "ConstitutiveDensity.h" +#include "FluidDensity.h" +#include "PermeabilityData.h" +#include "PhaseTransitionData.h" +#include "Porosity.h" +#include "Saturation.h" +#include "SolidCompressibility.h" +#include "Viscosity.h" + +namespace ProcessLib::TH2M +{ +namespace ConstitutiveRelations +{ +template +struct FW1Data +{ + GlobalDimMatrix A; +}; + +template +struct FW1Model +{ + void eval(AdvectionData const& advection_data, + FluidDensityData const& fluid_density_data, + FW1Data& fW_1) const; +}; + +extern template struct FW1Model<2>; +extern template struct FW1Model<3>; + +struct FW2Data +{ + double a = nan; +}; + +struct FW2DerivativeData +{ + double dp_GR = nan; + double dp_cap = nan; + double dT = nan; +}; + +struct FW2Model +{ + void eval(BiotData const biot_data, + CapillaryPressureData const pCap, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + SolidCompressibilityData const beta_p_SR, + FW2Data& fW_2) const; + + void dEval(BiotData const& biot_data, + CapillaryPressureData const pCap, + ConstituentDensityData const& constituent_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + PorosityDerivativeData const& porosity_d_data, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + SaturationDataDeriv const& dS_L_dp_cap, + SolidCompressibilityData const& beta_p_SR, + FW2DerivativeData& dfW_2) const; +}; + +struct FW3aData +{ + double a = nan; +}; + +struct FW3aDerivativeData +{ + double dp_GR = nan; + double dp_cap = nan; + double dT = nan; +}; + +struct FW3aModel +{ + void eval( + double const dt, + ConstituentDensityData const& constituent_density_data, + PrevState const& constituent_density_data_prev, + PrevState const& rho_W_LR_prev, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + FW3aData& fW_3a) const; + + void dEval( + double const dt, + ConstituentDensityData const& constituent_density_data, + PhaseTransitionData const& phase_transition_data, + PrevState const& constituent_density_data_prev, + PrevState const& rho_W_LR_prev, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + SaturationDataDeriv const& dS_L_dp_cap, + FW3aDerivativeData& dfW_3a) const; +}; + +template +struct FW4LWpGData +{ + GlobalDimMatrix L; +}; + +template +struct FW4LWpGDerivativeData +{ + GlobalDimMatrix dp_GR; + GlobalDimMatrix dp_cap; + GlobalDimMatrix dT; +}; + +template +struct FW4LWpGModel +{ + void eval(AdvectionData const& advection_data, + FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + FW4LWpGData& fW_4_LWpG) const; + + void dEval(ConstituentDensityData const& constituent_density_data, + PermeabilityData const& permeability_data, + PhaseTransitionData const& phase_transition_data, + PureLiquidDensityData const& rho_W_LR, + SaturationDataDeriv const& dS_L_dp_cap, + ViscosityData const& viscosity_data, + FW4LWpGDerivativeData& dfW_4_LWpG) const; +}; + +extern template struct FW4LWpGModel<2>; +extern template struct FW4LWpGModel<3>; + +template +struct FW4LWpCData +{ + GlobalDimMatrix L; +}; + +template +struct FW4LWpCDerivativeData +{ + GlobalDimMatrix dp_GR; + GlobalDimMatrix dp_cap; + GlobalDimMatrix dT; +}; + +template +struct FW4LWpCModel +{ + void eval(AdvectionData const& advection_data, + FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + FW4LWpCData& fW_4_LWpC) const; + + void dEval(AdvectionData const& advection_data, + FluidDensityData const& fluid_density_data, + PermeabilityData const& permeability_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + SaturationDataDeriv const& dS_L_dp_cap, + ViscosityData const& viscosity_data, + FW4LWpCDerivativeData& dfW_4_LWpC) const; +}; + +extern template struct FW4LWpCModel<2>; +extern template struct FW4LWpCModel<3>; + +template +struct FW4LWTData +{ + GlobalDimMatrix L; +}; + +template +struct FW4LWTModel +{ + void eval(FluidDensityData const& fluid_density_data, + PhaseTransitionData const& phase_transition_data, + PorosityData const& porosity_data, + SaturationData const& S_L_data, + FW4LWTData& fW_4_LWT) const; +}; + +struct FW4MWpGData +{ + double m = nan; +}; + +struct FW4MWpGModel +{ + void eval(BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + SolidCompressibilityData const& beta_p_SR, + FW4MWpGData& fW_4_MWpG) const; +}; + +struct FW4MWpCData +{ + double m = nan; + double ml = nan; +}; + +struct FW4MWpCModel +{ + void eval(BiotData const& biot_data, + CapillaryPressureData const pCap, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + PrevState const& S_L_data_prev, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + SolidCompressibilityData const& beta_p_SR, + FW4MWpCData& fW_4_MWpC) const; +}; + +struct FW4MWTData +{ + double m = nan; +}; + +template +struct FW4MWTModel +{ + void eval( + BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + PorosityData const& porosity_data, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + SolidThermalExpansionData const& s_therm_exp_data, + FW4MWTData& fW_4_MWT) const; +}; + +extern template struct FW4MWTModel<2>; +extern template struct FW4MWTModel<3>; + +struct FW4MWuData +{ + double m = nan; +}; + +struct FW4MWuModel +{ + void eval(BiotData const& biot_data, + ConstituentDensityData const& constituent_density_data, + PureLiquidDensityData const& rho_W_LR, + SaturationData const& S_L_data, + FW4MWuData& fW_4_MWu) const; +}; +} // namespace ConstitutiveRelations +} // namespace ProcessLib::TH2M diff --git a/ProcessLib/TH2M/TH2MFEM-impl.h b/ProcessLib/TH2M/TH2MFEM-impl.h index 2f14daaaeeb..154248a716b 100644 --- a/ProcessLib/TH2M/TH2MFEM-impl.h +++ b/ProcessLib/TH2M/TH2MFEM-impl.h @@ -152,16 +152,19 @@ TH2MLocalAssemblerprocess_data_.reference_temperature(t, pos)[0]}; - GlobalDimVectorType const gradpGR = gradNp * gas_pressure; - GlobalDimVectorType const gradpCap = gradNp * capillary_pressure; - GlobalDimVectorType const gradT = gradNp * temperature; + ConstitutiveRelations::GasPressureGradientData const + grad_p_GR{gradNp * gas_pressure}; + ConstitutiveRelations::CapillaryPressureGradientData< + DisplacementDim> const grad_p_cap{gradNp * capillary_pressure}; + ConstitutiveRelations::TemperatureGradientData const + grad_T{gradNp * temperature}; // medium properties models.elastic_tangent_stiffness_model.eval({pos, t, dt}, T_data, @@ -175,10 +178,9 @@ TH2MLocalAssembler( gradNu, Nu, x_coord, this->is_axially_symmetric_); - auto& eps = ip_out.eps_data.eps; - eps.noalias() = Bu * displacement; + ip_out.eps_data.eps.noalias() = Bu * displacement; models.S_L_model.eval({pos, t, dt}, media_data, pCap_data, - current_state.S_L_data, ip_cv.dS_L_dp_cap); + current_state.S_L_data); models.chi_S_L_model.eval({pos, t, dt}, media_data, current_state.S_L_data, ip_cv.chi_S_L); @@ -223,7 +225,7 @@ TH2MLocalAssembler{ + this->process_data_.specific_body_force}, + ip_cv.volumetric_body_force); + + models.diffusion_velocity_model.eval(grad_p_cap, + grad_p_GR, + ip_out.mass_mole_fractions_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + current_state.S_L_data, + grad_T, + ip_out.diffusion_velocity_data); + + models.solid_enthalpy_model.eval(ip_cv.solid_heat_capacity_data, T_data, + ip_out.solid_enthalpy_data); + + models.internal_energy_model.eval(ip_out.fluid_density_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + current_state.S_L_data, + ip_out.solid_density_data, + ip_out.solid_enthalpy_data, + current_state.internal_energy_data); + + models.effective_volumetric_enthalpy_model.eval( + ip_out.fluid_density_data, + ip_out.fluid_enthalpy_data, + ip_out.porosity_data, + current_state.S_L_data, + ip_out.solid_density_data, + ip_out.solid_enthalpy_data, + ip_cv.effective_volumetric_enthalpy_data); + + models.fC_1_model.eval(ip_cv.advection_data, ip_out.fluid_density_data, + ip_cv.fC_1); - 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; - - // dk_over_mu_G_dp_GR = ip_out.permeability_data.Ki * - // ip_out.permeability_data.dk_rel_G_dS_L * - // (ds_L_dp_GR = 0) / - // ip_cv.viscosity_data.mu_GR = 0; - // dk_over_mu_L_dp_GR = ip_out.permeability_data.Ki * - // ip_out.permeability_data.dk_rel_L_dS_L * - // (ds_L_dp_GR = 0) / - // ip_cv.viscosity_data.mu_LR = 0; - ip_cv.dk_over_mu_G_dp_cap = ip_out.permeability_data.Ki * - ip_out.permeability_data.dk_rel_G_dS_L * - ip_cv.dS_L_dp_cap() / - ip_cv.viscosity_data.mu_GR; - ip_cv.dk_over_mu_L_dp_cap = ip_out.permeability_data.Ki * - ip_out.permeability_data.dk_rel_L_dS_L * - ip_cv.dS_L_dp_cap() / - 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; - - ip_cv.drho_GR_h_w_eff_dp_GR_Npart = - c.drho_GR_dp_GR * ip_out.enthalpy_data.h_G * - ip_out.darcy_velocity_data.w_GS + - ip_out.fluid_density_data.rho_GR * ip_out.enthalpy_data.h_G * - k_over_mu_G * c.drho_GR_dp_GR * b; - ip_cv.drho_GR_h_w_eff_dp_GR_gradNpart = - ip_out.fluid_density_data.rho_GR * ip_out.enthalpy_data.h_G * - k_over_mu_G - - ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L * - k_over_mu_L; - - ip_cv.drho_LR_h_w_eff_dp_cap_Npart = - -drho_LR_dp_cap * ip_out.enthalpy_data.h_L * - ip_out.darcy_velocity_data.w_LS - - ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L * - k_over_mu_L * drho_LR_dp_cap * b; - ip_cv.drho_LR_h_w_eff_dp_cap_gradNpart = - // TODO (naumov) why the minus sign?????? - ip_out.fluid_density_data.rho_LR * ip_out.enthalpy_data.h_L * - k_over_mu_L; - - ip_cv.drho_GR_h_w_eff_dT = - c.drho_GR_dT * ip_out.enthalpy_data.h_G * - ip_out.darcy_velocity_data.w_GS + - ip_out.fluid_density_data.rho_GR * c.dh_G_dT * - ip_out.darcy_velocity_data.w_GS + - c.drho_LR_dT * ip_out.enthalpy_data.h_L * - ip_out.darcy_velocity_data.w_LS + - ip_out.fluid_density_data.rho_LR * c.dh_L_dT * - ip_out.darcy_velocity_data.w_LS; - // TODO (naumov) + k_over_mu_G * drho_GR_dT * b + k_over_mu_L * - // drho_LR_dT * b - - // Derivatives of s_G * rho_C_GR_dot + s_L * rho_C_LR_dot abbreviated - // here with S_rho_C_eff. - double const s_L = current_state.S_L_data.S_L; - double const s_G = 1. - s_L; - double const rho_C_FR = - s_G * current_state.constituent_density_data.rho_C_GR + - s_L * current_state.constituent_density_data.rho_C_LR; - double const rho_W_FR = - s_G * current_state.constituent_density_data.rho_W_GR + - s_L * current_state.rho_W_LR(); - // TODO (naumov) Extend for partially saturated media. - constexpr double drho_C_GR_dp_cap = 0; - if (dt == 0.) + if (!this->process_data_.apply_mass_lumping) { - ip_cv.dfC_3a_dp_GR = 0.; - ip_cv.dfC_3a_dp_cap = 0.; - ip_cv.dfC_3a_dT = 0.; + models.fC_2a_model.eval(ip_cv.biot_data, + pCap_data, + current_state.constituent_density_data, + ip_out.porosity_data, + current_state.S_L_data, + ip_cv.beta_p_SR, + ip_cv.fC_2a); } - else + models.fC_3a_model.eval(dt, + current_state.constituent_density_data, + prev_state.constituent_density_data, + current_state.S_L_data, + ip_cv.fC_3a); + + models.fC_4_LCpG_model.eval(ip_cv.advection_data, + ip_out.fluid_density_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + current_state.S_L_data, + ip_cv.fC_4_LCpG); + + models.fC_4_LCpC_model.eval(ip_cv.advection_data, + ip_out.fluid_density_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + current_state.S_L_data, + ip_cv.fC_4_LCpC); + + models.fC_4_LCT_model.eval(ip_out.fluid_density_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + current_state.S_L_data, + ip_cv.fC_4_LCT); + + models.fC_4_MCpG_model.eval(ip_cv.biot_data, + current_state.constituent_density_data, + ip_out.porosity_data, + current_state.S_L_data, + ip_cv.beta_p_SR, + ip_cv.fC_4_MCpG); + + models.fC_4_MCpC_model.eval(ip_cv.biot_data, + pCap_data, + current_state.constituent_density_data, + ip_out.porosity_data, + prev_state.S_L_data, + current_state.S_L_data, + ip_cv.beta_p_SR, + ip_cv.fC_4_MCpC); + + models.fC_4_MCT_model.eval(ip_cv.biot_data, + current_state.constituent_density_data, + ip_out.porosity_data, + current_state.S_L_data, + ip_cv.s_therm_exp_data, + ip_cv.fC_4_MCT); + + models.fC_4_MCu_model.eval(ip_cv.biot_data, + current_state.constituent_density_data, + current_state.S_L_data, + ip_cv.fC_4_MCu); + + models.fW_1_model.eval(ip_cv.advection_data, ip_out.fluid_density_data, + ip_cv.fW_1); + + if (!this->process_data_.apply_mass_lumping) { - double const rho_C_GR_dot = - (current_state.constituent_density_data.rho_C_GR - - prev_state.constituent_density_data->rho_C_GR) / - dt; - double const rho_C_LR_dot = - (current_state.constituent_density_data.rho_C_LR - - prev_state.constituent_density_data->rho_C_LR) / - dt; - ip_cv.dfC_3a_dp_GR = - /*(ds_G_dp_GR = 0) * rho_C_GR_dot +*/ s_G * c.drho_C_GR_dp_GR / - dt + - /*(ds_L_dp_GR = 0) * rho_C_LR_dot +*/ s_L * c.drho_C_LR_dp_GR / - dt; - ip_cv.dfC_3a_dp_cap = ds_G_dp_cap * rho_C_GR_dot + - s_G * drho_C_GR_dp_cap / dt + - ip_cv.dS_L_dp_cap() * rho_C_LR_dot - - s_L * c.drho_C_LR_dp_LR / dt; - ip_cv.dfC_3a_dT = - s_G * c.drho_C_GR_dT / dt + s_L * c.drho_C_LR_dT / dt; + models.fW_2_model.eval(ip_cv.biot_data, + pCap_data, + current_state.constituent_density_data, + ip_out.porosity_data, + current_state.rho_W_LR, + current_state.S_L_data, + ip_cv.beta_p_SR, + ip_cv.fW_2); } + models.fW_3a_model.eval(dt, + current_state.constituent_density_data, + prev_state.constituent_density_data, + prev_state.rho_W_LR, + current_state.rho_W_LR, + current_state.S_L_data, + ip_cv.fW_3a); + + models.fW_4_LWpG_model.eval(ip_cv.advection_data, + ip_out.fluid_density_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + current_state.S_L_data, + ip_cv.fW_4_LWpG); + + models.fW_4_LWpC_model.eval(ip_cv.advection_data, + ip_out.fluid_density_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + current_state.S_L_data, + ip_cv.fW_4_LWpC); + + models.fW_4_LWT_model.eval(ip_out.fluid_density_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + current_state.S_L_data, + ip_cv.fW_4_LWT); + + models.fW_4_MWpG_model.eval(ip_cv.biot_data, + current_state.constituent_density_data, + ip_out.porosity_data, + current_state.rho_W_LR, + current_state.S_L_data, + ip_cv.beta_p_SR, + ip_cv.fW_4_MWpG); + + models.fW_4_MWpC_model.eval(ip_cv.biot_data, + pCap_data, + current_state.constituent_density_data, + ip_out.porosity_data, + prev_state.S_L_data, + current_state.rho_W_LR, + current_state.S_L_data, + ip_cv.beta_p_SR, + ip_cv.fW_4_MWpC); + + models.fW_4_MWT_model.eval(ip_cv.biot_data, + current_state.constituent_density_data, + ip_out.porosity_data, + current_state.rho_W_LR, + current_state.S_L_data, + ip_cv.s_therm_exp_data, + ip_cv.fW_4_MWT); + + models.fW_4_MWu_model.eval(ip_cv.biot_data, + current_state.constituent_density_data, + current_state.rho_W_LR, + current_state.S_L_data, + ip_cv.fW_4_MWu); + + models.fT_1_model.eval(dt, + current_state.internal_energy_data, + prev_state.internal_energy_data, + ip_cv.fT_1); + + // --------------------------------------------------------------------- + // Derivatives for Jacobian + // --------------------------------------------------------------------- + + models.darcy_velocity_model.eval( + grad_p_cap, + ip_out.fluid_density_data, + grad_p_GR, + ip_out.permeability_data, + ConstitutiveRelations::SpecificBodyForceData{ + 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.fluid_density_data, + ip_out.fluid_enthalpy_data, + ip_cv.fT_2); + + models.fT_3_model.eval( + current_state.constituent_density_data, + ip_out.darcy_velocity_data, + ip_out.diffusion_velocity_data, + ip_out.fluid_density_data, + ip_cv.phase_transition_data, + ConstitutiveRelations::SpecificBodyForceData{ + this->process_data_.specific_body_force}, + ip_cv.fT_3); + + models.fu_2_KupC_model.eval(ip_cv.biot_data, ip_cv.chi_S_L, + ip_cv.fu_2_KupC); + } - double const drho_C_FR_dp_GR = - /*(ds_G_dp_GR = 0) * current_state.constituent_density_data.rho_C_GR - +*/ - s_G * c.drho_C_GR_dp_GR + - /*(ds_L_dp_GR = 0) * current_state.constituent_density_data.rho_C_LR - +*/ - s_L * c.drho_C_LR_dp_GR; - ip_cv.dfC_4_MCpG_dp_GR = - drho_C_FR_dp_GR * (ip_cv.biot_data() - ip_out.porosity_data.phi) * - ip_cv.beta_p_SR(); - - double const drho_C_FR_dT = s_G * c.drho_C_GR_dT + s_L * c.drho_C_LR_dT; - ip_cv.dfC_4_MCpG_dT = - drho_C_FR_dT * (ip_cv.biot_data() - ip_out.porosity_data.phi) * - ip_cv.beta_p_SR() - - rho_C_FR * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR(); - - ip_cv.dfC_4_MCT_dT = - drho_C_FR_dT * (ip_cv.biot_data() - ip_out.porosity_data.phi) * - ip_cv.s_therm_exp_data.beta_T_SR + return {ip_constitutive_data, ip_constitutive_variables}; +} + +template +std::vector> +TH2MLocalAssembler:: + updateConstitutiveVariablesDerivatives( + Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev, + double const t, double const dt, + std::vector< + ConstitutiveRelations::ConstitutiveData> const& + ip_constitutive_data, + std::vector< + ConstitutiveRelations::ConstitutiveTempData> const& + ip_constitutive_variables, + ConstitutiveRelations::ConstitutiveModels const& + models) +{ + [[maybe_unused]] auto const matrix_size = + gas_pressure_size + capillary_pressure_size + temperature_size + + displacement_size; + + assert(local_x.size() == matrix_size); + + auto const temperature = + local_x.template segment(temperature_index); + auto const temperature_prev = + local_x_prev.template segment(temperature_index); + + auto const capillary_pressure = + local_x.template segment( + capillary_pressure_index); + + auto const& medium = + *this->process_data_.media_map.getMedium(this->element_.getID()); + ConstitutiveRelations::MediaData media_data{medium}; + + unsigned const n_integration_points = + this->integration_method_.getNumberOfPoints(); + + std::vector> + ip_d_data(n_integration_points); + + for (unsigned ip = 0; ip < n_integration_points; ip++) + { + auto const& ip_data = _ip_data[ip]; + auto& ip_dd = ip_d_data[ip]; + auto const& ip_cd = ip_constitutive_data[ip]; + auto const& ip_cv = ip_constitutive_variables[ip]; + auto const& ip_out = this->output_data_[ip]; + auto const& current_state = this->current_states_[ip]; + auto const& prev_state = this->prev_states_[ip]; + + auto const& Nu = ip_data.N_u; + auto const& Np = ip_data.N_p; + auto const& NT = Np; + + ParameterLib::SpatialPosition const pos{ + std::nullopt, this->element_.getID(), ip, + MathLib::Point3d( + NumLib::interpolateCoordinates( + this->element_, Nu))}; + + double const T = NT.dot(temperature); + double const T_prev = NT.dot(temperature_prev); + ConstitutiveRelations::TemperatureData const T_data{T, T_prev}; + ConstitutiveRelations::CapillaryPressureData const pCap_data{ + Np.dot(capillary_pressure)}; + + models.S_L_model.dEval({pos, t, dt}, media_data, pCap_data, + ip_dd.dS_L_dp_cap); + + models.advection_model.dEval(current_state.constituent_density_data, + ip_out.permeability_data, + ip_cv.viscosity_data, + ip_dd.dS_L_dp_cap, + ip_cv.phase_transition_data, + ip_dd.advection_d_data); + + models.porosity_model.dEval( + {pos, t, dt}, media_data, ip_out.porosity_data, ip_dd.dS_L_dp_cap, #ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION - + rho_C_FR * (ip_cv.biot_data() - ip_cv.porosity_d_data.dphi_dT) * - ip_cv.s_therm_exp_data.beta_T_SR -#endif - ; - - ip_cv.dfC_4_MCu_dT = drho_C_FR_dT * ip_cv.biot_data(); - - ip_cv.dfC_2a_dp_GR = - -ip_out.porosity_data.phi * c.drho_C_GR_dp_GR - - drho_C_FR_dp_GR * pCap * - (ip_cv.biot_data() - ip_out.porosity_data.phi) * - ip_cv.beta_p_SR(); - - double const drho_C_FR_dp_cap = - ds_G_dp_cap * current_state.constituent_density_data.rho_C_GR + - s_G * drho_C_GR_dp_cap + - ip_cv.dS_L_dp_cap() * - current_state.constituent_density_data.rho_C_LR - - s_L * c.drho_C_LR_dp_LR; - - ip_cv.dfC_2a_dp_cap = - ip_out.porosity_data.phi * (-c.drho_C_LR_dp_LR - drho_C_GR_dp_cap) - - drho_C_FR_dp_cap * pCap * - (ip_cv.biot_data() - ip_out.porosity_data.phi) * - ip_cv.beta_p_SR() + - rho_C_FR * (ip_cv.biot_data() - ip_out.porosity_data.phi) * - ip_cv.beta_p_SR(); - - ip_cv.dfC_2a_dT = - ip_cv.porosity_d_data.dphi_dT * - (current_state.constituent_density_data.rho_C_LR - - current_state.constituent_density_data.rho_C_GR) + - ip_out.porosity_data.phi * (c.drho_C_LR_dT - c.drho_C_GR_dT) - - drho_C_FR_dT * pCap * - (ip_cv.biot_data() - ip_out.porosity_data.phi) * - ip_cv.beta_p_SR() + - rho_C_FR * pCap * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR(); - - ip_cv.dadvection_C_dp_GR = c.drho_C_GR_dp_GR * k_over_mu_G - // + rhoCGR * (dk_over_mu_G_dp_GR = 0) - // + rhoCLR * (dk_over_mu_L_dp_GR = 0) - + c.drho_C_LR_dp_GR * k_over_mu_L; - - ip_cv.dadvection_C_dp_cap = - //(drho_C_GR_dp_cap = 0) * k_over_mu_G - current_state.constituent_density_data.rho_C_GR * - ip_cv.dk_over_mu_G_dp_cap + - (-c.drho_C_LR_dp_LR) * k_over_mu_L + - current_state.constituent_density_data.rho_C_LR * - ip_cv.dk_over_mu_L_dp_cap; - - ip_cv.dfC_4_LCpG_dT = - c.drho_C_GR_dT * k_over_mu_G + c.drho_C_LR_dT * k_over_mu_L - // + ip_cv.ddiffusion_C_p_dT TODO (naumov) - ; - - double const drho_W_FR_dp_GR = - /*(ds_G_dp_GR = 0) * current_state.constituent_density_data.rho_W_GR - +*/ - s_G * c.drho_W_GR_dp_GR + - /*(ds_L_dp_GR = 0) * current_state.rho_W_LR() +*/ - s_L * c.drho_W_LR_dp_GR; - double const drho_W_FR_dp_cap = - ds_G_dp_cap * current_state.constituent_density_data.rho_W_GR + - s_G * c.drho_W_GR_dp_cap + - ip_cv.dS_L_dp_cap() * current_state.rho_W_LR() - - s_L * c.drho_W_LR_dp_LR; - double const drho_W_FR_dT = s_G * c.drho_W_GR_dT + s_L * c.drho_W_LR_dT; - - ip_cv.dfW_2a_dp_GR = - ip_out.porosity_data.phi * (c.drho_W_LR_dp_GR - c.drho_W_GR_dp_GR); - ip_cv.dfW_2b_dp_GR = drho_W_FR_dp_GR * pCap * - (ip_cv.biot_data() - ip_out.porosity_data.phi) * - ip_cv.beta_p_SR(); - ip_cv.dfW_2a_dp_cap = ip_out.porosity_data.phi * - (-c.drho_W_LR_dp_LR - c.drho_W_GR_dp_cap); - ip_cv.dfW_2b_dp_cap = - drho_W_FR_dp_cap * pCap * - (ip_cv.biot_data() - ip_out.porosity_data.phi) * - ip_cv.beta_p_SR() + - rho_W_FR * (ip_cv.biot_data() - ip_out.porosity_data.phi) * - ip_cv.beta_p_SR(); - - ip_cv.dfW_2a_dT = - ip_cv.porosity_d_data.dphi_dT * - (current_state.rho_W_LR() - - current_state.constituent_density_data.rho_W_GR) + - ip_out.porosity_data.phi * (c.drho_W_LR_dT - c.drho_W_GR_dT); - ip_cv.dfW_2b_dT = - drho_W_FR_dT * pCap * - (ip_cv.biot_data() - ip_out.porosity_data.phi) * - ip_cv.beta_p_SR() - - rho_W_FR * pCap * ip_cv.porosity_d_data.dphi_dT * ip_cv.beta_p_SR(); - - if (dt == 0.) + ip_cv.biot_data, ip_out.eps_data, ip_cv.s_therm_exp_data, +#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION + ip_dd.porosity_d_data); + + models.thermal_conductivity_model.dEval( + {pos, t, dt}, media_data, T_data, ip_out.porosity_data, + ip_dd.porosity_d_data, current_state.S_L_data, + ip_dd.thermal_conductivity_d_data); + + models.solid_density_model.dEval({pos, t, dt}, media_data, T_data, +#ifdef NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION + ip_cv.biot_data, ip_out.eps_data, + ip_cv.s_therm_exp_data, +#endif // NON_CONSTANT_SOLID_PHASE_VOLUME_FRACTION + ip_dd.solid_density_d_data); + + models.internal_energy_model.dEval( + ip_out.fluid_density_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + ip_dd.porosity_d_data, + current_state.S_L_data, + ip_out.solid_density_data, + ip_dd.solid_density_d_data, + ip_out.solid_enthalpy_data, + ip_cv.solid_heat_capacity_data, + ip_dd.effective_volumetric_internal_energy_d_data); + + models.effective_volumetric_enthalpy_model.dEval( + ip_out.fluid_density_data, + ip_out.fluid_enthalpy_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + ip_dd.porosity_d_data, + current_state.S_L_data, + ip_out.solid_density_data, + ip_dd.solid_density_d_data, + ip_out.solid_enthalpy_data, + ip_cv.solid_heat_capacity_data, + ip_dd.effective_volumetric_enthalpy_d_data); + if (!this->process_data_.apply_mass_lumping) { - ip_cv.dfW_3a_dp_GR = 0.; - ip_cv.dfW_3a_dp_cap = 0.; - ip_cv.dfW_3a_dT = 0.; + models.fC_2a_model.dEval(ip_cv.biot_data, + pCap_data, + current_state.constituent_density_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + ip_dd.porosity_d_data, + current_state.S_L_data, + ip_dd.dS_L_dp_cap, + ip_cv.beta_p_SR, + ip_dd.dfC_2a); } - else + models.fC_3a_model.dEval(dt, + current_state.constituent_density_data, + prev_state.constituent_density_data, + ip_cv.phase_transition_data, + current_state.S_L_data, + ip_dd.dS_L_dp_cap, + ip_dd.dfC_3a); + + models.fC_4_LCpG_model.dEval(ip_out.permeability_data, + ip_cv.viscosity_data, + ip_cv.phase_transition_data, + ip_dd.advection_d_data, + ip_dd.dfC_4_LCpG); + + models.fC_4_LCpC_model.dEval(current_state.constituent_density_data, + ip_out.permeability_data, + ip_cv.phase_transition_data, + ip_dd.dS_L_dp_cap, + ip_cv.viscosity_data, + ip_dd.dfC_4_LCpC); + + models.fC_4_MCpG_model.dEval(ip_cv.biot_data, + current_state.constituent_density_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + ip_dd.porosity_d_data, + current_state.S_L_data, + ip_cv.beta_p_SR, + ip_dd.dfC_4_MCpG); + + models.fC_4_MCT_model.dEval(ip_cv.biot_data, + current_state.constituent_density_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + ip_dd.porosity_d_data, + current_state.S_L_data, + ip_cv.s_therm_exp_data, + ip_dd.dfC_4_MCT); + + models.fC_4_MCu_model.dEval(ip_cv.biot_data, + ip_cv.phase_transition_data, + current_state.S_L_data, + ip_dd.dfC_4_MCu); + + if (!this->process_data_.apply_mass_lumping) { - double const rho_W_GR_dot = - (current_state.constituent_density_data.rho_W_GR - - prev_state.constituent_density_data->rho_W_GR) / - dt; - double const rho_W_LR_dot = - (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt; - - ip_cv.dfW_3a_dp_GR = - /*(ds_G_dp_GR = 0) * rho_W_GR_dot +*/ s_G * c.drho_W_GR_dp_GR / - dt + - /*(ds_L_dp_GR = 0) * rho_W_LR_dot +*/ s_L * c.drho_W_LR_dp_GR / - dt; - ip_cv.dfW_3a_dp_cap = ds_G_dp_cap * rho_W_GR_dot + - s_G * c.drho_W_GR_dp_cap / dt + - ip_cv.dS_L_dp_cap() * rho_W_LR_dot - - s_L * c.drho_W_LR_dp_LR / dt; - ip_cv.dfW_3a_dT = - s_G * c.drho_W_GR_dT / dt + s_L * c.drho_W_LR_dT / dt; + models.fW_2_model.dEval(ip_cv.biot_data, + pCap_data, + current_state.constituent_density_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + ip_dd.porosity_d_data, + current_state.rho_W_LR, + current_state.S_L_data, + ip_dd.dS_L_dp_cap, + ip_cv.beta_p_SR, + ip_dd.dfW_2); } - ip_cv.dfW_4_LWpG_a_dp_GR = c.drho_W_GR_dp_GR * k_over_mu_G - // + rhoWGR * (dk_over_mu_G_dp_GR = 0) - + c.drho_W_LR_dp_GR * k_over_mu_L - // + rhoWLR * (dk_over_mu_L_dp_GR = 0) - ; - ip_cv.dfW_4_LWpG_a_dp_cap = - c.drho_W_GR_dp_cap * k_over_mu_G + - current_state.constituent_density_data.rho_W_GR * - ip_cv.dk_over_mu_G_dp_cap + - -c.drho_W_LR_dp_LR * k_over_mu_L + - current_state.rho_W_LR() * ip_cv.dk_over_mu_L_dp_cap; - - ip_cv.dfW_4_LWpG_a_dT = - c.drho_W_GR_dT * k_over_mu_G - //+ rhoWGR * (dk_over_mu_G_dT != 0 TODO for mu_G(T)) - + c.drho_W_LR_dT * k_over_mu_L - //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_G(T)) - ; - - // TODO (naumov) for dxmW*/d* != 0 - ip_cv.dfW_4_LWpG_d_dp_GR = - Eigen::Matrix::Zero(); - ip_cv.dfW_4_LWpG_d_dp_cap = - Eigen::Matrix::Zero(); - ip_cv.dfW_4_LWpG_d_dT = - Eigen::Matrix::Zero(); - - ip_cv.dfW_4_LWpC_a_dp_GR = c.drho_W_LR_dp_GR * k_over_mu_L - //+ rhoWLR * (dk_over_mu_L_dp_GR = 0) - ; - ip_cv.dfW_4_LWpC_a_dp_cap = - -c.drho_W_LR_dp_LR * k_over_mu_L + - current_state.rho_W_LR() * ip_cv.dk_over_mu_L_dp_cap; - ip_cv.dfW_4_LWpC_a_dT = c.drho_W_LR_dT * k_over_mu_L - //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T)) - ; - - // TODO (naumov) for dxmW*/d* != 0 - ip_cv.dfW_4_LWpC_d_dp_GR = - Eigen::Matrix::Zero(); - ip_cv.dfW_4_LWpC_d_dp_cap = - Eigen::Matrix::Zero(); - ip_cv.dfW_4_LWpC_d_dT = - Eigen::Matrix::Zero(); - - ip_cv.dfC_4_LCpC_a_dp_GR = c.drho_C_LR_dp_GR * k_over_mu_L - //+ rhoCLR * (dk_over_mu_L_dp_GR = 0) - ; - ip_cv.dfC_4_LCpC_a_dp_cap = - -c.drho_C_LR_dp_LR * k_over_mu_L + - current_state.constituent_density_data.rho_C_LR * - ip_cv.dk_over_mu_L_dp_cap; - ip_cv.dfC_4_LCpC_a_dT = c.drho_W_LR_dT * k_over_mu_L - //+ rhoWLR * (dk_over_mu_L_dT != 0 TODO for mu_L(T)) - ; - - // TODO (naumov) for dxmW*/d* != 0 - ip_cv.dfC_4_LCpC_d_dp_GR = - Eigen::Matrix::Zero(); - ip_cv.dfC_4_LCpC_d_dp_cap = - Eigen::Matrix::Zero(); - ip_cv.dfC_4_LCpC_d_dT = - Eigen::Matrix::Zero(); + models.fW_3a_model.dEval(dt, + current_state.constituent_density_data, + ip_cv.phase_transition_data, + prev_state.constituent_density_data, + prev_state.rho_W_LR, + current_state.rho_W_LR, + current_state.S_L_data, + ip_dd.dS_L_dp_cap, + ip_dd.dfW_3a); + + models.fW_4_LWpG_model.dEval(current_state.constituent_density_data, + ip_out.permeability_data, + ip_cv.phase_transition_data, + current_state.rho_W_LR, + ip_dd.dS_L_dp_cap, + ip_cv.viscosity_data, + ip_dd.dfW_4_LWpG); + + models.fW_4_LWpC_model.dEval(ip_cv.advection_data, + ip_out.fluid_density_data, + ip_out.permeability_data, + ip_cv.phase_transition_data, + ip_out.porosity_data, + current_state.rho_W_LR, + current_state.S_L_data, + ip_dd.dS_L_dp_cap, + ip_cv.viscosity_data, + ip_dd.dfW_4_LWpC); + + models.fT_1_model.dEval( + dt, ip_dd.effective_volumetric_internal_energy_d_data, ip_dd.dfT_1); + + models.fT_2_model.dEval( + ip_out.darcy_velocity_data, + ip_out.fluid_density_data, + ip_out.fluid_enthalpy_data, + ip_out.permeability_data, + ip_cv.phase_transition_data, + ConstitutiveRelations::SpecificBodyForceData{ + this->process_data_.specific_body_force}, + ip_cv.viscosity_data, + ip_dd.dfT_2); + + models.fu_1_KuT_model.dEval(ip_cd.s_mech_data, ip_cv.s_therm_exp_data, + ip_dd.dfu_1_KuT); + + models.fu_2_KupC_model.dEval(ip_cv.biot_data, + ip_cv.chi_S_L, + pCap_data, + ip_dd.dS_L_dp_cap, + ip_dd.dfu_2_KupC); } - return {ip_constitutive_data, ip_constitutive_variables}; + return ip_d_data; } template ( @@ -1073,299 +1079,117 @@ void TH2MLocalAssembler< typename BMatricesType::BMatrixType>( gradNu, Nu, x_coord, this->is_axially_symmetric_); - auto const BuT = Bu.transpose().eval(); + auto const NTN = (Np.transpose() * Np).eval(); + auto const BTI2N = (Bu.transpose() * Invariants::identity2 * Np).eval(); double const pCap = Np.dot(capillary_pressure); double const pCap_prev = Np.dot(capillary_pressure_prev); - double const beta_T_SR = ip_cv.s_therm_exp_data.beta_T_SR; - - auto const I = - Eigen::Matrix::Identity(); - - double const sD_G = - ip_cv.phase_transition_data.diffusion_coefficient_vapour; - double const sD_L = - ip_cv.phase_transition_data.diffusion_coefficient_solute; - - GlobalDimMatrixType const D_C_G = sD_G * I; - GlobalDimMatrixType const D_W_G = sD_G * I; - GlobalDimMatrixType const D_C_L = sD_L * I; - GlobalDimMatrixType const D_W_L = sD_L * I; - auto const s_L = current_state.S_L_data.S_L; - auto const s_G = 1. - s_L; auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt; - auto& alpha_B = ip_cv.biot_data(); - auto& beta_p_SR = ip_cv.beta_p_SR(); - auto const& b = this->process_data_.specific_body_force; - // volume fraction - auto const phi_G = s_G * ip_out.porosity_data.phi; - auto const phi_L = s_L * ip_out.porosity_data.phi; - auto const phi_S = 1. - ip_out.porosity_data.phi; - - auto const rhoGR = ip_out.fluid_density_data.rho_GR; - auto const rhoLR = ip_out.fluid_density_data.rho_LR; - - // effective density - auto const rho = phi_G * rhoGR + phi_L * rhoLR + - phi_S * ip_out.solid_density_data.rho_SR; - - // abbreviations - const double rho_C_FR = - s_G * current_state.constituent_density_data.rho_C_GR + - s_L * current_state.constituent_density_data.rho_C_LR; - const double rho_W_FR = - s_G * current_state.constituent_density_data.rho_W_GR + - s_L * current_state.rho_W_LR(); - - auto const rho_C_GR_dot = - (current_state.constituent_density_data.rho_C_GR - - prev_state.constituent_density_data->rho_C_GR) / - dt; - auto const rho_C_LR_dot = - (current_state.constituent_density_data.rho_C_LR - - prev_state.constituent_density_data->rho_C_LR) / - dt; - auto const rho_W_GR_dot = - (current_state.constituent_density_data.rho_W_GR - - prev_state.constituent_density_data->rho_W_GR) / - dt; - auto const rho_W_LR_dot = - (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt; - - 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; - // --------------------------------------------------------------------- // C-component equation // --------------------------------------------------------------------- - MCpG.noalias() += NpT * rho_C_FR * - (alpha_B - ip_out.porosity_data.phi) * beta_p_SR * - Np * w; - MCpC.noalias() -= NpT * rho_C_FR * - (alpha_B - ip_out.porosity_data.phi) * beta_p_SR * - s_L * Np * w; + MCpG.noalias() += NTN * (ip_cv.fC_4_MCpG.m * w); + MCpC.noalias() += NTN * (ip_cv.fC_4_MCpC.m * w); if (this->process_data_.apply_mass_lumping) { if (pCap - pCap_prev != 0.) // avoid division by Zero { MCpC.noalias() += - NpT * - (ip_out.porosity_data.phi * - (current_state.constituent_density_data.rho_C_LR - - current_state.constituent_density_data.rho_C_GR) - - rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) * - beta_p_SR) * - s_L_dot * dt / (pCap - pCap_prev) * Np * w; + NTN * (ip_cv.fC_4_MCpC.ml / (pCap - pCap_prev) * w); } } - MCT.noalias() -= NpT * rho_C_FR * (alpha_B - ip_out.porosity_data.phi) * - beta_T_SR * Np * w; - MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w; - - using DisplacementDimMatrix = - Eigen::Matrix; - - DisplacementDimMatrix const advection_C_G = - current_state.constituent_density_data.rho_C_GR * k_over_mu_G; - DisplacementDimMatrix const advection_C_L = - current_state.constituent_density_data.rho_C_LR * k_over_mu_L; - - DisplacementDimMatrix const diffusion_CGpGR = - -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpGR; - DisplacementDimMatrix const diffusion_CLpGR = - -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpGR; - - DisplacementDimMatrix const diffusion_CGpCap = - -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpCap; - DisplacementDimMatrix const diffusion_CLpCap = - -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpCap; - - DisplacementDimMatrix const diffusion_CGT = - -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dT; - DisplacementDimMatrix const diffusion_CLT = - -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dT; + MCT.noalias() += NTN * (ip_cv.fC_4_MCT.m * w); + MCu.noalias() += BTI2N.transpose() * (ip_cv.fC_4_MCu.m * w); - DisplacementDimMatrix const advection_C = advection_C_G + advection_C_L; - DisplacementDimMatrix const diffusion_C_pGR = - diffusion_CGpGR + diffusion_CLpGR; - DisplacementDimMatrix const diffusion_C_pCap = - diffusion_CGpCap + diffusion_CLpCap; + LCpG.noalias() += gradNpT * ip_cv.fC_4_LCpG.L * gradNp * w; - DisplacementDimMatrix const diffusion_C_T = - diffusion_CGT + diffusion_CLT; + LCpC.noalias() += gradNpT * ip_cv.fC_4_LCpC.L * gradNp * w; - LCpG.noalias() += - gradNpT * (advection_C + diffusion_C_pGR) * gradNp * w; + LCT.noalias() += gradNpT * ip_cv.fC_4_LCT.L * gradNp * w; - LCpC.noalias() += - gradNpT * (diffusion_C_pCap - advection_C_L) * gradNp * w; - - LCT.noalias() += gradNpT * (diffusion_C_T)*gradNp * w; - - fC.noalias() += - gradNpT * (advection_C_G * rhoGR + advection_C_L * rhoLR) * b * w; + fC.noalias() += gradNpT * ip_cv.fC_1.A * b * w; if (!this->process_data_.apply_mass_lumping) { - fC.noalias() -= - NpT * - (ip_out.porosity_data.phi * - (current_state.constituent_density_data.rho_C_LR - - current_state.constituent_density_data.rho_C_GR) - - rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) * - beta_p_SR) * - s_L_dot * w; + fC.noalias() -= NpT * (ip_cv.fC_2a.a * s_L_dot * w); } // fC_III - fC.noalias() -= NpT * ip_out.porosity_data.phi * - (s_G * rho_C_GR_dot + s_L * rho_C_LR_dot) * w; + fC.noalias() -= NpT * (ip_out.porosity_data.phi * ip_cv.fC_3a.a * w); // --------------------------------------------------------------------- // W-component equation // --------------------------------------------------------------------- - MWpG.noalias() += NpT * rho_W_FR * - (alpha_B - ip_out.porosity_data.phi) * beta_p_SR * - Np * w; - MWpC.noalias() -= NpT * rho_W_FR * - (alpha_B - ip_out.porosity_data.phi) * beta_p_SR * - s_L * Np * w; + MWpG.noalias() += NTN * (ip_cv.fW_4_MWpG.m * w); + MWpC.noalias() += NTN * (ip_cv.fW_4_MWpC.m * w); if (this->process_data_.apply_mass_lumping) { if (pCap - pCap_prev != 0.) // avoid division by Zero { MWpC.noalias() += - NpT * - (ip_out.porosity_data.phi * - (current_state.rho_W_LR() - - current_state.constituent_density_data.rho_W_GR) - - rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) * - beta_p_SR) * - s_L_dot * dt / (pCap - pCap_prev) * Np * w; + NTN * (ip_cv.fW_4_MWpC.ml / (pCap - pCap_prev) * w); } } - MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) * - beta_T_SR * Np * w; - - MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w; - - DisplacementDimMatrix const advection_W_G = - current_state.constituent_density_data.rho_W_GR * k_over_mu_G; - DisplacementDimMatrix const advection_W_L = - current_state.rho_W_LR() * k_over_mu_L; - - DisplacementDimMatrix const diffusion_WGpGR = - phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpGR; - DisplacementDimMatrix const diffusion_WLpGR = - phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpGR; - - DisplacementDimMatrix const diffusion_WGpCap = - phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpCap; - DisplacementDimMatrix const diffusion_WLpCap = - phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpCap; - - DisplacementDimMatrix const diffusion_WGT = - phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dT; - DisplacementDimMatrix const diffusion_WLT = - phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dT; - - DisplacementDimMatrix const advection_W = advection_W_G + advection_W_L; - DisplacementDimMatrix const diffusion_W_pGR = - diffusion_WGpGR + diffusion_WLpGR; - DisplacementDimMatrix const diffusion_W_pCap = - diffusion_WGpCap + diffusion_WLpCap; + MWT.noalias() += NTN * (ip_cv.fW_4_MWT.m * w); - DisplacementDimMatrix const diffusion_W_T = - diffusion_WGT + diffusion_WLT; + MWu.noalias() += BTI2N.transpose() * (ip_cv.fW_4_MWu.m * w); - LWpG.noalias() += - gradNpT * (advection_W + diffusion_W_pGR) * gradNp * w; + LWpG.noalias() += gradNpT * ip_cv.fW_4_LWpG.L * gradNp * w; - LWpC.noalias() += - gradNpT * (diffusion_W_pCap - advection_W_L) * gradNp * w; + LWpC.noalias() += gradNpT * ip_cv.fW_4_LWpC.L * gradNp * w; - LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w; + LWT.noalias() += gradNpT * ip_cv.fW_4_LWT.L * gradNp * w; - fW.noalias() += - gradNpT * (advection_W_G * rhoGR + advection_W_L * rhoLR) * b * w; + fW.noalias() += gradNpT * ip_cv.fW_1.A * b * w; if (!this->process_data_.apply_mass_lumping) { - fW.noalias() -= - NpT * - (ip_out.porosity_data.phi * - (current_state.rho_W_LR() - - current_state.constituent_density_data.rho_W_GR) - - rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) * - beta_p_SR) * - s_L_dot * w; + fW.noalias() -= NpT * (ip_cv.fW_2.a * s_L_dot * w); } - fW.noalias() -= NpT * ip_out.porosity_data.phi * - (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w; + fW.noalias() -= NpT * (ip_out.porosity_data.phi * ip_cv.fW_3a.a * w); // --------------------------------------------------------------------- // - temperature equation // --------------------------------------------------------------------- - MTu.noalias() += NTT * - ip_cv.effective_volumetric_enthalpy_data.rho_h_eff * - mT * Bu * w; + MTu.noalias() += + BTI2N.transpose() * + (ip_cv.effective_volumetric_enthalpy_data.rho_h_eff * w); KTT.noalias() += gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w; - auto const rho_u_eff_dot = (current_state.internal_energy_data() - - **prev_state.internal_energy_data) / - dt; - fT.noalias() -= NTT * rho_u_eff_dot * w; - - fT.noalias() += gradNTT * - (rhoGR * ip_out.enthalpy_data.h_G * - ip_out.darcy_velocity_data.w_GS + - rhoLR * ip_out.enthalpy_data.h_L * - ip_out.darcy_velocity_data.w_LS) * - w; - - fT.noalias() += gradNTT * - (current_state.constituent_density_data.rho_C_GR * - ip_cv.phase_transition_data.hCG * - ip_out.diffusion_velocity_data.d_CG + - current_state.constituent_density_data.rho_W_GR * - ip_cv.phase_transition_data.hWG * - ip_out.diffusion_velocity_data.d_WG) * - w; - - fT.noalias() += NTT * - (rhoGR * ip_out.darcy_velocity_data.w_GS.transpose() + - rhoLR * ip_out.darcy_velocity_data.w_LS.transpose()) * - b * w; + fT.noalias() -= NTT * (ip_cv.fT_1.m * w); + + fT.noalias() += gradNTT * ip_cv.fT_2.A * w; + + fT.noalias() += gradNTT * ip_cv.fT_3.gradN * w; + + fT.noalias() += NTT * (ip_cv.fT_3.N * w); // --------------------------------------------------------------------- // - displacement equation // --------------------------------------------------------------------- - KUpG.noalias() -= (BuT * alpha_B * m * Np) * w; + KUpG.noalias() -= BTI2N * (ip_cv.biot_data() * w); - KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w; + KUpC.noalias() += BTI2N * (ip_cv.fu_2_KupC.m * w); - fU.noalias() -= (BuT * current_state.eff_stress_data.sigma - - N_u_op(Nu).transpose() * rho * b) * - w; + fU.noalias() -= + (Bu.transpose() * current_state.eff_stress_data.sigma - + N_u_op(Nu).transpose() * ip_cv.volumetric_body_force()) * + w; if (this->process_data_.apply_mass_lumping) { @@ -1515,10 +1339,17 @@ void TH2MLocalAssembler(local_x.data(), local_x.size()), + Eigen::Map(local_x_prev.data(), + local_x_prev.size()), + t, dt, ip_constitutive_data, ip_constitutive_variables, models); + for (unsigned int_point = 0; int_point < n_integration_points; int_point++) { auto& ip = _ip_data[int_point]; auto& ip_cd = ip_constitutive_data[int_point]; + auto& ip_dd = ip_d_data[int_point]; auto& ip_cv = ip_constitutive_variables[int_point]; auto& ip_out = this->output_data_[int_point]; auto& current_state = this->current_states_[int_point]; @@ -1546,9 +1377,6 @@ void TH2MLocalAssembler( @@ -1560,7 +1388,8 @@ void TH2MLocalAssembler( gradNu, Nu, x_coord, this->is_axially_symmetric_); - auto const BuT = Bu.transpose().eval(); + auto const NTN = (Np.transpose() * Np).eval(); + auto const BTI2N = (Bu.transpose() * Invariants::identity2 * Np).eval(); double const div_u_dot = Invariants::trace(Bu * (displacement - displacement_prev) / dt); @@ -1576,402 +1405,232 @@ void TH2MLocalAssembler::Identity(); - - double const sD_G = - ip_cv.phase_transition_data.diffusion_coefficient_vapour; - double const sD_L = - ip_cv.phase_transition_data.diffusion_coefficient_solute; - - GlobalDimMatrixType const D_C_G = sD_G * I; - GlobalDimMatrixType const D_W_G = sD_G * I; - GlobalDimMatrixType const D_C_L = sD_L * I; - GlobalDimMatrixType const D_W_L = sD_L * I; - - auto& s_L = current_state.S_L_data.S_L; - auto const s_G = 1. - s_L; + auto const& s_L = current_state.S_L_data.S_L; auto const s_L_dot = (s_L - prev_state.S_L_data->S_L) / dt; - auto const alpha_B = ip_cv.biot_data(); - auto const beta_p_SR = ip_cv.beta_p_SR(); - auto const& b = this->process_data_.specific_body_force; - // volume fraction - auto const phi_G = s_G * ip_out.porosity_data.phi; - auto const phi_L = s_L * ip_out.porosity_data.phi; - auto const phi_S = 1. - ip_out.porosity_data.phi; - - auto const rhoGR = ip_out.fluid_density_data.rho_GR; - auto const rhoLR = ip_out.fluid_density_data.rho_LR; - - // effective density - auto const rho = phi_G * rhoGR + phi_L * rhoLR + - phi_S * ip_out.solid_density_data.rho_SR; - - // abbreviations - const double rho_C_FR = - s_G * current_state.constituent_density_data.rho_C_GR + - s_L * current_state.constituent_density_data.rho_C_LR; - const double rho_W_FR = - s_G * current_state.constituent_density_data.rho_W_GR + - s_L * current_state.rho_W_LR(); - - auto const rho_C_GR_dot = - (current_state.constituent_density_data.rho_C_GR - - prev_state.constituent_density_data->rho_C_GR) / - dt; - auto const rho_C_LR_dot = - (current_state.constituent_density_data.rho_C_LR - - prev_state.constituent_density_data->rho_C_LR) / - dt; - auto const rho_W_GR_dot = - (current_state.constituent_density_data.rho_W_GR - - prev_state.constituent_density_data->rho_W_GR) / - dt; - auto const rho_W_LR_dot = - (current_state.rho_W_LR() - **prev_state.rho_W_LR) / dt; - - 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; - // --------------------------------------------------------------------- // C-component equation // --------------------------------------------------------------------- - MCpG.noalias() += NpT * rho_C_FR * - (alpha_B - ip_out.porosity_data.phi) * beta_p_SR * - Np * w; - MCpC.noalias() -= NpT * rho_C_FR * - (alpha_B - ip_out.porosity_data.phi) * beta_p_SR * - s_L * Np * w; + MCpG.noalias() += NTN * (ip_cv.fC_4_MCpG.m * w); + MCpC.noalias() += NTN * (ip_cv.fC_4_MCpC.m * w); if (this->process_data_.apply_mass_lumping) { if (pCap - pCap_prev != 0.) // avoid division by Zero { MCpC.noalias() += - NpT * - (ip_out.porosity_data.phi * - (current_state.constituent_density_data.rho_C_LR - - current_state.constituent_density_data.rho_C_GR) - - rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) * - beta_p_SR) * - s_L_dot * dt / (pCap - pCap_prev) * Np * w; + NTN * (ip_cv.fC_4_MCpC.ml / (pCap - pCap_prev) * w); } } - MCT.noalias() -= NpT * rho_C_FR * (alpha_B - ip_out.porosity_data.phi) * - beta_T_SR * Np * w; + MCT.noalias() += NTN * (ip_cv.fC_4_MCT.m * w); // d (fC_4_MCT * T_dot)/d T local_Jac .template block(C_index, temperature_index) - .noalias() += NpT * ip_cv.dfC_4_MCT_dT * (T - T_prev) / dt * NT * w; + .noalias() += NTN * (ip_dd.dfC_4_MCT.dT * (T - T_prev) / dt * w); - MCu.noalias() += NpT * rho_C_FR * alpha_B * mT * Bu * w; + MCu.noalias() += BTI2N.transpose() * (ip_cv.fC_4_MCu.m * w); // d (fC_4_MCu * u_dot)/d T local_Jac .template block(C_index, temperature_index) - .noalias() += NpT * ip_cv.dfC_4_MCu_dT * div_u_dot * NT * w; - - GlobalDimMatrixType const advection_C_G = - current_state.constituent_density_data.rho_C_GR * k_over_mu_G; - GlobalDimMatrixType const advection_C_L = - current_state.constituent_density_data.rho_C_LR * k_over_mu_L; - GlobalDimMatrixType const diffusion_C_G_p = - -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dpGR; - GlobalDimMatrixType const diffusion_C_L_p = - -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dpLR; - GlobalDimMatrixType const diffusion_C_G_T = - -phi_G * rhoGR * D_C_G * ip_cv.phase_transition_data.dxmWG_dT; - GlobalDimMatrixType const diffusion_C_L_T = - -phi_L * rhoLR * D_C_L * ip_cv.phase_transition_data.dxmWL_dT; - - GlobalDimMatrixType const advection_C = advection_C_G + advection_C_L; - GlobalDimMatrixType const diffusion_C_p = - diffusion_C_G_p + diffusion_C_L_p; - GlobalDimMatrixType const diffusion_C_T = - diffusion_C_G_T + diffusion_C_L_T; - - LCpG.noalias() += gradNpT * (advection_C + diffusion_C_p) * gradNp * w; + .noalias() += NTN * (ip_dd.dfC_4_MCu.dT * div_u_dot * w); + + LCpG.noalias() += gradNpT * ip_cv.fC_4_LCpG.L * gradNp * w; // d (fC_4_LCpG * grad p_GR)/d p_GR local_Jac.template block(C_index, C_index).noalias() += - gradNpT * - (ip_cv.dadvection_C_dp_GR - // + ip_cv.ddiffusion_C_p_dp_GR TODO (naumov) - ) * - gradpGR * Np * w; + gradNpT * ip_dd.dfC_4_LCpG.dp_GR * gradpGR * Np * w; // d (fC_4_LCpG * grad p_GR)/d p_cap local_Jac.template block(C_index, W_index).noalias() += - gradNpT * - (ip_cv.dadvection_C_dp_cap - // + ip_cv.ddiffusion_C_p_dp_GR TODO (naumov) - ) * - gradpGR * Np * w; + gradNpT * ip_dd.dfC_4_LCpG.dp_cap * gradpGR * Np * w; // d (fC_4_LCpG * grad p_GR)/d T local_Jac .template block(C_index, temperature_index) - .noalias() += gradNpT * ip_cv.dfC_4_LCpG_dT * gradpGR * NT * w; + .noalias() += gradNpT * ip_dd.dfC_4_LCpG.dT * gradpGR * NT * w; // d (fC_4_MCpG * p_GR_dot)/d p_GR local_Jac.template block(C_index, C_index).noalias() += - NpT * ip_cv.dfC_4_MCpG_dp_GR * (pGR - pGR_prev) / dt * Np * w; + NTN * (ip_dd.dfC_4_MCpG.dp_GR * (pGR - pGR_prev) / dt * w); // d (fC_4_MCpG * p_GR_dot)/d T local_Jac .template block(C_index, temperature_index) .noalias() += - NpT * ip_cv.dfC_4_MCpG_dT * (pGR - pGR_prev) / dt * NT * w; + NTN * (ip_dd.dfC_4_MCpG.dT * (pGR - pGR_prev) / dt * w); - LCpC.noalias() -= - gradNpT * (advection_C_L + diffusion_C_L_p) * gradNp * w; + LCpC.noalias() -= gradNpT * ip_cv.fC_4_LCpC.L * gradNp * w; /* TODO (naumov) This part is not tested by any of the current ctests. // d (fC_4_LCpC * grad p_cap)/d p_GR local_Jac.template block(C_index, C_index).noalias() += - gradNpT * - (ip_cv.dfC_4_LCpC_a_dp_GR - // + ip_cv.dfC_4_LCpC_d_dp_GR TODO (naumov) - ) * - gradpCap * Np * w; + gradNpT * ip_dd.dfC_4_LCpC.dp_GR * gradpCap * Np * w; // d (fC_4_LCpC * grad p_cap)/d p_cap local_Jac.template block(C_index, W_index).noalias() += - gradNpT * - (ip_cv.dfC_4_LCpC_a_dp_cap - // + ip_cv.dfC_4_LCpC_d_dp_cap TODO (naumov) - ) * - gradpCap * Np * w; + gradNpT * ip_dd.dfC_4_LCpC.dp_cap * gradpCap * Np * w; local_Jac .template block(C_index, temperature_index) - .noalias() += gradNpT * - (ip_cv.dfC_4_LCpC_a_dT - // + ip_cv.dfC_4_LCpC_d_dT TODO (naumov) - ) * - gradpCap * Np * w; + .noalias() += gradNpT * ip_dd.dfC_4_LCpC.dT * gradpCap * Np * w; */ - LCT.noalias() += gradNpT * diffusion_C_T * gradNp * w; + LCT.noalias() += gradNpT * ip_cv.fC_4_LCT.L * gradNp * w; // fC_1 - fC.noalias() += - gradNpT * (advection_C_G * rhoGR + advection_C_L * rhoLR) * b * w; + fC.noalias() += gradNpT * ip_cv.fC_1.A * b * w; if (!this->process_data_.apply_mass_lumping) { // fC_2 = \int a * s_L_dot - auto const a = - ip_out.porosity_data.phi * - (current_state.constituent_density_data.rho_C_LR - - current_state.constituent_density_data.rho_C_GR) - - rho_C_FR * pCap * (alpha_B - ip_out.porosity_data.phi) * - beta_p_SR; - fC.noalias() -= NpT * a * s_L_dot * w; + fC.noalias() -= NpT * (ip_cv.fC_2a.a * s_L_dot * w); local_Jac.template block(C_index, C_index) .noalias() += - NpT * - (ip_cv.dfC_2a_dp_GR * s_L_dot /*- a * (ds_L_dp_GR = 0) / dt*/) * - Np * w; + NTN * ((ip_dd.dfC_2a.dp_GR * s_L_dot + /*- ip_cv.fC_2a.a * (ds_L_dp_GR = 0) / dt*/) * + w); local_Jac.template block(C_index, W_index) .noalias() += - NpT * - (ip_cv.dfC_2a_dp_cap * s_L_dot + a * ip_cv.dS_L_dp_cap() / dt) * - Np * w; + NTN * ((ip_dd.dfC_2a.dp_cap * s_L_dot + + ip_cv.fC_2a.a * ip_dd.dS_L_dp_cap() / dt) * + w); local_Jac .template block(C_index, temperature_index) - .noalias() += NpT * ip_cv.dfC_2a_dT * s_L_dot * NT * w; + .noalias() += NTN * (ip_dd.dfC_2a.dT * s_L_dot * w); } { // fC_3 = \int phi * a - double const a = s_G * rho_C_GR_dot + s_L * rho_C_LR_dot; - fC.noalias() -= NpT * ip_out.porosity_data.phi * a * w; + fC.noalias() -= + NpT * (ip_out.porosity_data.phi * ip_cv.fC_3a.a * w); local_Jac.template block(C_index, C_index) .noalias() += - NpT * ip_out.porosity_data.phi * ip_cv.dfC_3a_dp_GR * Np * w; + NTN * (ip_out.porosity_data.phi * ip_dd.dfC_3a.dp_GR * w); local_Jac.template block(C_index, W_index) .noalias() += - NpT * ip_out.porosity_data.phi * ip_cv.dfC_3a_dp_cap * Np * w; + NTN * (ip_out.porosity_data.phi * ip_dd.dfC_3a.dp_cap * w); local_Jac .template block(C_index, temperature_index) - .noalias() += NpT * - (ip_cv.porosity_d_data.dphi_dT * a + - ip_out.porosity_data.phi * ip_cv.dfC_3a_dT) * - NT * w; + .noalias() += + NTN * ((ip_dd.porosity_d_data.dphi_dT * ip_cv.fC_3a.a + + ip_out.porosity_data.phi * ip_dd.dfC_3a.dT) * + w); } // --------------------------------------------------------------------- // W-component equation // --------------------------------------------------------------------- - MWpG.noalias() += NpT * rho_W_FR * - (alpha_B - ip_out.porosity_data.phi) * beta_p_SR * - Np * w; - MWpC.noalias() -= NpT * rho_W_FR * - (alpha_B - ip_out.porosity_data.phi) * beta_p_SR * - s_L * Np * w; + MWpG.noalias() += NTN * (ip_cv.fW_4_MWpG.m * w); + MWpC.noalias() += NTN * (ip_cv.fW_4_MWpC.m * w); if (this->process_data_.apply_mass_lumping) { if (pCap - pCap_prev != 0.) // avoid division by Zero { MWpC.noalias() += - NpT * - (ip_out.porosity_data.phi * - (current_state.rho_W_LR() - - current_state.constituent_density_data.rho_W_GR) - - rho_W_FR * pCap * (alpha_B - ip_out.porosity_data.phi) * - beta_p_SR) * - s_L_dot * dt / (pCap - pCap_prev) * Np * w; + NTN * (ip_cv.fW_4_MWpC.ml / (pCap - pCap_prev) * w); } } - MWT.noalias() -= NpT * rho_W_FR * (alpha_B - ip_out.porosity_data.phi) * - beta_T_SR * Np * w; - - MWu.noalias() += NpT * rho_W_FR * alpha_B * mT * Bu * w; + MWT.noalias() += NTN * (ip_cv.fW_4_MWT.m * w); - GlobalDimMatrixType const advection_W_G = - current_state.constituent_density_data.rho_W_GR * k_over_mu_G; - GlobalDimMatrixType const advection_W_L = - current_state.rho_W_LR() * k_over_mu_L; - GlobalDimMatrixType const diffusion_W_G_p = - phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dpGR; - GlobalDimMatrixType const diffusion_W_L_p = - phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dpLR; - GlobalDimMatrixType const diffusion_W_G_T = - phi_G * rhoGR * D_W_G * ip_cv.phase_transition_data.dxmWG_dT; - GlobalDimMatrixType const diffusion_W_L_T = - phi_L * rhoLR * D_W_L * ip_cv.phase_transition_data.dxmWL_dT; + MWu.noalias() += BTI2N.transpose() * (ip_cv.fW_4_MWu.m * w); - GlobalDimMatrixType const advection_W = advection_W_G + advection_W_L; - GlobalDimMatrixType const diffusion_W_p = - diffusion_W_G_p + diffusion_W_L_p; - GlobalDimMatrixType const diffusion_W_T = - diffusion_W_G_T + diffusion_W_L_T; - - LWpG.noalias() += gradNpT * (advection_W + diffusion_W_p) * gradNp * w; + LWpG.noalias() += gradNpT * ip_cv.fW_4_LWpG.L * gradNp * w; // fW_4 LWpG' parts; LWpG = \int grad (a + d) grad local_Jac.template block(W_index, C_index).noalias() += - gradNpT * (ip_cv.dfW_4_LWpG_a_dp_GR + ip_cv.dfW_4_LWpG_d_dp_GR) * - gradpGR * Np * w; + gradNpT * ip_dd.dfW_4_LWpG.dp_GR * gradpGR * Np * w; local_Jac.template block(W_index, W_index).noalias() += - gradNpT * (ip_cv.dfW_4_LWpG_a_dp_cap + ip_cv.dfW_4_LWpG_d_dp_cap) * - gradpGR * Np * w; + gradNpT * ip_dd.dfW_4_LWpG.dp_cap * gradpGR * Np * w; local_Jac .template block(W_index, temperature_index) - .noalias() += gradNpT * - (ip_cv.dfW_4_LWpG_a_dT + ip_cv.dfW_4_LWpG_d_dT) * - gradpGR * NT * w; + .noalias() += gradNpT * ip_dd.dfW_4_LWpG.dT * gradpGR * NT * w; - LWpC.noalias() -= - gradNpT * (advection_W_L + diffusion_W_L_p) * gradNp * w; + LWpC.noalias() += gradNpT * ip_cv.fW_4_LWpC.L * gradNp * w; // fW_4 LWp_cap' parts; LWpC = \int grad (a + d) grad local_Jac.template block(W_index, C_index).noalias() -= - gradNpT * (ip_cv.dfW_4_LWpC_a_dp_GR + ip_cv.dfW_4_LWpC_d_dp_GR) * - gradpCap * Np * w; + gradNpT * ip_dd.dfW_4_LWpC.dp_GR * gradpCap * Np * w; local_Jac.template block(W_index, W_index).noalias() -= - gradNpT * (ip_cv.dfW_4_LWpC_a_dp_cap + ip_cv.dfW_4_LWpC_d_dp_cap) * - gradpCap * Np * w; + gradNpT * ip_dd.dfW_4_LWpC.dp_cap * gradpCap * Np * w; local_Jac .template block(W_index, temperature_index) - .noalias() -= gradNpT * - (ip_cv.dfW_4_LWpC_a_dT + ip_cv.dfW_4_LWpC_d_dT) * - gradpCap * NT * w; + .noalias() -= gradNpT * ip_dd.dfW_4_LWpC.dT * gradpCap * NT * w; - LWT.noalias() += gradNpT * (diffusion_W_T)*gradNp * w; + LWT.noalias() += gradNpT * ip_cv.fW_4_LWT.L * gradNp * w; // fW_1 - fW.noalias() += - gradNpT * (advection_W_G * rhoGR + advection_W_L * rhoLR) * b * w; + fW.noalias() += gradNpT * ip_cv.fW_1.A * b * w; - // fW_2 = \int (f - g) * s_L_dot + // fW_2 = \int a * s_L_dot if (!this->process_data_.apply_mass_lumping) { - double const f = ip_out.porosity_data.phi * - (current_state.rho_W_LR() - - current_state.constituent_density_data.rho_W_GR); - double const g = rho_W_FR * pCap * - (alpha_B - ip_out.porosity_data.phi) * beta_p_SR; - - fW.noalias() -= NpT * (f - g) * s_L_dot * w; + fW.noalias() -= NpT * (ip_cv.fW_2.a * s_L_dot * w); local_Jac.template block(W_index, C_index) - .noalias() += NpT * (ip_cv.dfW_2a_dp_GR - ip_cv.dfW_2b_dp_GR) * - s_L_dot * Np * w; + .noalias() += NTN * (ip_dd.dfW_2.dp_GR * s_L_dot * w); // sign negated because of dp_cap = -dp_LR // TODO (naumov) Had to change the sign to get equal Jacobian WW // blocks in A2 Test. Where is the error? local_Jac.template block(W_index, W_index) - .noalias() += - NpT * - ((ip_cv.dfW_2a_dp_cap - ip_cv.dfW_2b_dp_cap) * s_L_dot + - (f - g) * ip_cv.dS_L_dp_cap() / dt) * - Np * w; + .noalias() += NTN * ((ip_dd.dfW_2.dp_cap * s_L_dot + + ip_cv.fW_2.a * ip_dd.dS_L_dp_cap() / dt) * + w); local_Jac .template block(W_index, temperature_index) - .noalias() += - NpT * (ip_cv.dfW_2a_dT - ip_cv.dfW_2b_dT) * s_L_dot * Np * w; + .noalias() += NTN * (ip_dd.dfW_2.dT * s_L_dot * w); } // fW_3 = \int phi * a - fW.noalias() -= NpT * ip_out.porosity_data.phi * - (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) * w; + fW.noalias() -= NpT * (ip_out.porosity_data.phi * ip_cv.fW_3a.a * w); local_Jac.template block(W_index, C_index).noalias() += - NpT * ip_out.porosity_data.phi * ip_cv.dfW_3a_dp_GR * Np * w; + NTN * (ip_out.porosity_data.phi * ip_dd.dfW_3a.dp_GR * w); local_Jac.template block(W_index, W_index).noalias() += - NpT * ip_out.porosity_data.phi * ip_cv.dfW_3a_dp_cap * Np * w; + NTN * (ip_out.porosity_data.phi * ip_dd.dfW_3a.dp_cap * w); local_Jac .template block(W_index, temperature_index) - .noalias() += NpT * - (ip_cv.porosity_d_data.dphi_dT * - (s_G * rho_W_GR_dot + s_L * rho_W_LR_dot) + - ip_out.porosity_data.phi * ip_cv.dfW_3a_dT) * - NT * w; + .noalias() += + NTN * ((ip_dd.porosity_d_data.dphi_dT * ip_cv.fW_3a.a + + ip_out.porosity_data.phi * ip_dd.dfW_3a.dT) * + w); // --------------------------------------------------------------------- // - temperature equation // --------------------------------------------------------------------- - MTu.noalias() += NTT * - ip_cv.effective_volumetric_enthalpy_data.rho_h_eff * - mT * Bu * w; + MTu.noalias() += + BTI2N.transpose() * + (ip_cv.effective_volumetric_enthalpy_data.rho_h_eff * w); // dfT_4/dp_GR // d (MTu * u_dot)/dp_GR @@ -1979,8 +1638,8 @@ void TH2MLocalAssembler(temperature_index, C_index) .noalias() += - NTT * ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_GR * - div_u_dot * NT * w; + NTN * (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_GR * + div_u_dot * w); // dfT_4/dp_cap // d (MTu * u_dot)/dp_cap @@ -1988,8 +1647,9 @@ void TH2MLocalAssembler(temperature_index, W_index) .noalias() -= - NTT * ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_cap * - div_u_dot * NT * w; + NTN * + (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dp_cap * + div_u_dot * w); // dfT_4/dT // d (MTu * u_dot)/dT @@ -1997,8 +1657,8 @@ void TH2MLocalAssembler( temperature_index, temperature_index) .noalias() += - NTT * ip_cv.effective_volumetric_enthalpy_d_data.drho_h_eff_dT * - div_u_dot * NT * w; + NTN * (ip_dd.effective_volumetric_enthalpy_d_data.drho_h_eff_dT * + div_u_dot * w); KTT.noalias() += gradNTT * ip_cv.thermal_conductivity_data.lambda * gradNT * w; @@ -2022,57 +1682,41 @@ void TH2MLocalAssembler(temperature_index, W_index) .noalias() += gradNTT * - ip_cv.thermal_conductivity_data.dlambda_dp_cap * + ip_dd.thermal_conductivity_d_data.dlambda_dp_cap * gradT * Np * w; // d KTT/dT * T local_Jac .template block( temperature_index, temperature_index) - .noalias() += gradNTT * ip_cv.thermal_conductivity_data.dlambda_dT * - gradT * NT * w; + .noalias() += gradNTT * + ip_dd.thermal_conductivity_d_data.dlambda_dT * gradT * + NT * w; // fT_1 - auto const rho_u_eff_dot = (current_state.internal_energy_data() - - **prev_state.internal_energy_data) / - dt; - fT.noalias() -= NTT * rho_u_eff_dot * w; + fT.noalias() -= NTT * (ip_cv.fT_1.m * w); // dfT_1/dp_GR local_Jac .template block(temperature_index, C_index) - .noalias() += NpT * Np * - (ip_cv.effective_volumetric_internal_energy_d_data - .drho_u_eff_dp_GR / - dt * w); + .noalias() += NTN * (ip_dd.dfT_1.dp_GR * w); // dfT_1/dp_cap local_Jac .template block(temperature_index, W_index) - .noalias() += NpT * Np * - (ip_cv.effective_volumetric_internal_energy_d_data - .drho_u_eff_dp_cap / - dt * w); + .noalias() += NTN * (ip_dd.dfT_1.dp_cap * w); // dfT_1/dT // MTT local_Jac .template block( temperature_index, temperature_index) - .noalias() += - NTT * NT * - (ip_cv.effective_volumetric_internal_energy_d_data.drho_u_eff_dT / - dt * w); + .noalias() += NTN * (ip_dd.dfT_1.dT * w); // fT_2 - fT.noalias() += gradNTT * - (rhoGR * ip_out.enthalpy_data.h_G * - ip_out.darcy_velocity_data.w_GS + - rhoLR * ip_out.enthalpy_data.h_L * - ip_out.darcy_velocity_data.w_LS) * - w; + fT.noalias() += gradNTT * ip_cv.fT_2.A * w; // dfT_2/dp_GR local_Jac @@ -2080,9 +1724,9 @@ void TH2MLocalAssembler( temperature_index, temperature_index) - .noalias() -= gradNTT * ip_cv.drho_GR_h_w_eff_dT * NT * w; + .noalias() -= gradNTT * ip_dd.dfT_2.dT * NT * w; // fT_3 - fT.noalias() += NTT * - (rhoGR * ip_out.darcy_velocity_data.w_GS.transpose() + - rhoLR * ip_out.darcy_velocity_data.w_LS.transpose()) * - b * w; - - fT.noalias() += gradNTT * - (current_state.constituent_density_data.rho_C_GR * - ip_cv.phase_transition_data.hCG * - ip_out.diffusion_velocity_data.d_CG + - current_state.constituent_density_data.rho_W_GR * - ip_cv.phase_transition_data.hWG * - ip_out.diffusion_velocity_data.d_WG) * - w; + fT.noalias() += NTT * (ip_cv.fT_3.N * w); + + fT.noalias() += gradNTT * ip_cv.fT_3.gradN * w; // --------------------------------------------------------------------- // - displacement equation // --------------------------------------------------------------------- - KUpG.noalias() -= (BuT * alpha_B * m * Np) * w; + KUpG.noalias() -= BTI2N * (ip_cv.biot_data() * w); // dfU_2/dp_GR = dKUpG/dp_GR * p_GR + KUpG. The former is zero, the // latter is handled below. - KUpC.noalias() += (BuT * alpha_B * ip_cv.chi_S_L.chi_S_L * m * Np) * w; + KUpC.noalias() += BTI2N * (ip_cv.fu_2_KupC.m * w); // dfU_2/dp_cap = dKUpC/dp_cap * p_cap + KUpC. The former is handled // here, the latter below. local_Jac .template block(displacement_index, W_index) - .noalias() += BuT * alpha_B * ip_cv.chi_S_L.dchi_dS_L * - ip_cv.dS_L_dp_cap() * pCap * m * Np * w; + .noalias() += BTI2N * (ip_dd.dfu_2_KupC.dp_cap * w); local_Jac .template block( displacement_index, displacement_index) - .noalias() += BuT * ip_cd.s_mech_data.stiffness_tensor * Bu * w; + .noalias() += + Bu.transpose() * ip_cd.s_mech_data.stiffness_tensor * Bu * w; // fU_1 - fU.noalias() -= (BuT * current_state.eff_stress_data.sigma - - N_u_op(Nu).transpose() * rho * b) * - w; + fU.noalias() -= + (Bu.transpose() * current_state.eff_stress_data.sigma - + N_u_op(Nu).transpose() * ip_cv.volumetric_body_force()) * + w; // KuT local_Jac .template block( displacement_index, temperature_index) - .noalias() -= - BuT * - (ip_cd.s_mech_data.stiffness_tensor * - ip_cv.s_therm_exp_data.solid_linear_thermal_expansivity_vector) * - NT * w; + .noalias() -= Bu.transpose() * ip_dd.dfu_1_KuT.dT * NT * w; /* TODO (naumov) Test with gravity needed to check this Jacobian part. local_Jac diff --git a/ProcessLib/TH2M/TH2MFEM.h b/ProcessLib/TH2M/TH2MFEM.h index 0d88646ae1f..bad9ef32fd5 100644 --- a/ProcessLib/TH2M/TH2MFEM.h +++ b/ProcessLib/TH2M/TH2MFEM.h @@ -198,6 +198,19 @@ class TH2MLocalAssembler : public LocalAssemblerInterface ConstitutiveRelations::ConstitutiveModels const& models); + std::vector> + updateConstitutiveVariablesDerivatives( + Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev, + double const t, double const dt, + std::vector< + ConstitutiveRelations::ConstitutiveData> const& + ip_constitutive_data, + std::vector< + ConstitutiveRelations::ConstitutiveTempData> const& + ip_constitutive_variables, + ConstitutiveRelations::ConstitutiveModels const& + models); + private: using BMatricesType = BMatrixPolicyType; diff --git a/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp b/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp index bc74bda9cf0..7fd5c61d739 100644 --- a/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp +++ b/Tests/ProcessLib/TH2M/TestTH2MNoPhaseTransition.cpp @@ -108,7 +108,7 @@ TEST(ProcessLib, TH2MNoPhaseTransition) rhoWLR); ASSERT_NEAR(density_water, rhoWLR(), 1e-10); - CR::EnthalpyData enthalpy; + CR::FluidEnthalpyData enthalpy; CR::MassMoleFractionsData mass_mole_fractions; CR::FluidDensityData fluid_density; CR::VapourPartialPressureData vapour_pressure; diff --git a/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp b/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp index f636edf49f6..5dcc72fe30c 100644 --- a/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp +++ b/Tests/ProcessLib/TH2M/TestTH2MPhaseTransition.cpp @@ -218,7 +218,7 @@ TEST(ProcessLib, TH2MPhaseTransition) CR::PureLiquidDensityData rhoWLR; CR::PureLiquidDensityModel rhoWLR_model; - CR::EnthalpyData enthalpy; + CR::FluidEnthalpyData enthalpy; CR::MassMoleFractionsData mass_mole_fractions; CR::FluidDensityData fluid_density; CR::VapourPartialPressureData vapour_pressure; @@ -480,7 +480,7 @@ TEST(ProcessLib, TH2MPhaseTransitionConstRho) CR::PureLiquidDensityData rhoWLR; CR::PureLiquidDensityModel rhoWLR_model; - CR::EnthalpyData enthalpy; + CR::FluidEnthalpyData enthalpy; CR::MassMoleFractionsData mass_mole_fractions; CR::FluidDensityData fluid_density; CR::VapourPartialPressureData vapour_pressure;