From e49166da8f8d8e7932e1f00eeacaa71c03c419ff Mon Sep 17 00:00:00 2001 From: Anne van de Graaf Date: Tue, 28 May 2024 15:03:31 +0200 Subject: [PATCH] Towards element permeability matrices (#12408) Precalculate relative permeability values as preparation for calculating element permeability matrices. When needed, modify these values by accounting for the permeability update factors (i.e. to account for volume changes in geometrically linear elements). Other changes include: - Members `CalculateRetentionResponse` no longer compute the relative permeability. It is now pre-calculated. - Extracted members that calculate the relative permeability values for all integration points of an element. - Moved the calculation of permeability update factors to the transport equation utilities. - Extracted a utility function that calculates fluid pressures at all integration points of an element. - Inlined the wrapper member functions that calculate and add the permeability matrix, to avoid needless levels of indirection. - Removed member `PermeabilityUpdateFactor` from two `ElementVariables` data structures. - Removed the unused `ProcessInfo` data member from the retention laws. - Removed several comments that had no additional value. - Miscellaneous cleanup: * Reduced the scope of some variables. * Don't use `noalias` if no performance benefit is to be expected. * Prefer to use `this->` to refer to base members rather than a `using` statement. * Removed some unnecessary `KRATOS_TRY` and `KRATOS_CATCH` statements. * Removed several empty statements. * Unnamed some unused function parameters. --- .../thermal_dispersion_law.cpp | 5 +- .../thermal_dispersion_law.h | 2 +- .../thermal_filter_law.cpp | 2 +- .../custom_constitutive/thermal_filter_law.h | 2 +- .../custom_constitutive/thermal_law.h | 3 +- .../U_Pw_small_strain_FIC_element.cpp | 6 +- .../U_Pw_small_strain_element.cpp | 129 +++++++----------- .../U_Pw_small_strain_element.hpp | 8 +- .../U_Pw_small_strain_interface_element.cpp | 26 ++-- ...Pw_small_strain_link_interface_element.cpp | 3 +- .../U_Pw_updated_lagrangian_FIC_element.cpp | 6 +- .../U_Pw_updated_lagrangian_element.cpp | 6 +- .../small_strain_U_Pw_diff_order_element.cpp | 120 +++++++--------- .../small_strain_U_Pw_diff_order_element.hpp | 6 +- .../steady_state_Pw_element.cpp | 18 +-- .../steady_state_Pw_interface_element.cpp | 5 +- .../steady_state_Pw_interface_element.hpp | 1 - .../steady_state_Pw_piping_element.cpp | 5 +- .../steady_state_Pw_piping_element.hpp | 1 - .../custom_elements/transient_Pw_element.cpp | 55 +++----- .../custom_elements/transient_Pw_element.hpp | 5 +- .../transient_Pw_interface_element.cpp | 20 ++- .../transient_Pw_interface_element.hpp | 5 +- .../transient_Pw_line_element.h | 39 +++--- .../transient_thermal_element.h | 14 +- ...ted_lagrangian_U_Pw_diff_order_element.cpp | 6 +- .../custom_retention/retention_law.h | 9 +- .../transport_equation_utilities.hpp | 49 +++++-- .../cpp_tests/test_permeability_matrix.cpp | 8 +- .../cpp_tests/test_thermal_dispersion_law.cpp | 8 +- .../cpp_tests/test_thermal_filter_law.cpp | 3 +- .../test_transport_equation_utilities.cpp | 70 ++++++++++ 32 files changed, 326 insertions(+), 319 deletions(-) diff --git a/applications/GeoMechanicsApplication/custom_constitutive/thermal_dispersion_law.cpp b/applications/GeoMechanicsApplication/custom_constitutive/thermal_dispersion_law.cpp index 0420263c4271..0c1e3ce53b4d 100644 --- a/applications/GeoMechanicsApplication/custom_constitutive/thermal_dispersion_law.cpp +++ b/applications/GeoMechanicsApplication/custom_constitutive/thermal_dispersion_law.cpp @@ -26,12 +26,11 @@ GeoThermalDispersionLaw::GeoThermalDispersionLaw(std::size_t NumberOfDimensions) << "Got invalid number of dimensions: " << mNumberOfDimensions << std::endl; } -Matrix GeoThermalDispersionLaw::CalculateThermalDispersionMatrix(const Properties& rProp, - const ProcessInfo& rProcessInfo) const +Matrix GeoThermalDispersionLaw::CalculateThermalDispersionMatrix(const Properties& rProp) const { Matrix result = ZeroMatrix(mNumberOfDimensions, mNumberOfDimensions); - RetentionLaw::Parameters parameters(rProp, rProcessInfo); + RetentionLaw::Parameters parameters(rProp); auto retention_law = RetentionLawFactory::Clone(rProp); const double saturation = retention_law->CalculateSaturation(parameters); const double water_fraction = rProp[POROSITY] * saturation; diff --git a/applications/GeoMechanicsApplication/custom_constitutive/thermal_dispersion_law.h b/applications/GeoMechanicsApplication/custom_constitutive/thermal_dispersion_law.h index e3e8cdbd47e7..e61545557e6d 100644 --- a/applications/GeoMechanicsApplication/custom_constitutive/thermal_dispersion_law.h +++ b/applications/GeoMechanicsApplication/custom_constitutive/thermal_dispersion_law.h @@ -30,7 +30,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoThermalDispersionLaw : public Geo explicit GeoThermalDispersionLaw(SizeType NumberOfDimensions); - Matrix CalculateThermalDispersionMatrix(const Properties& rProp, const ProcessInfo& rProcessInfo) const override; + [[nodiscard]] Matrix CalculateThermalDispersionMatrix(const Properties& rProp) const override; private: std::size_t mNumberOfDimensions = 2; diff --git a/applications/GeoMechanicsApplication/custom_constitutive/thermal_filter_law.cpp b/applications/GeoMechanicsApplication/custom_constitutive/thermal_filter_law.cpp index 3ade36cd4aad..9f6bc3919a86 100644 --- a/applications/GeoMechanicsApplication/custom_constitutive/thermal_filter_law.cpp +++ b/applications/GeoMechanicsApplication/custom_constitutive/thermal_filter_law.cpp @@ -17,7 +17,7 @@ namespace Kratos { -Matrix GeoThermalFilterLaw::CalculateThermalDispersionMatrix(const Properties& rProp, const ProcessInfo&) const +Matrix GeoThermalFilterLaw::CalculateThermalDispersionMatrix(const Properties& rProp) const { return ScalarMatrix(1, 1, rProp[THERMAL_CONDUCTIVITY_WATER]); } diff --git a/applications/GeoMechanicsApplication/custom_constitutive/thermal_filter_law.h b/applications/GeoMechanicsApplication/custom_constitutive/thermal_filter_law.h index 795c6c02ca90..813ecb1716d1 100644 --- a/applications/GeoMechanicsApplication/custom_constitutive/thermal_filter_law.h +++ b/applications/GeoMechanicsApplication/custom_constitutive/thermal_filter_law.h @@ -23,7 +23,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoThermalFilterLaw : public GeoTher public: KRATOS_CLASS_POINTER_DEFINITION(GeoThermalFilterLaw); - Matrix CalculateThermalDispersionMatrix(const Properties& rProp, const ProcessInfo&) const override; + [[nodiscard]] Matrix CalculateThermalDispersionMatrix(const Properties& rProp) const override; private: friend class Serializer; diff --git a/applications/GeoMechanicsApplication/custom_constitutive/thermal_law.h b/applications/GeoMechanicsApplication/custom_constitutive/thermal_law.h index 2ed800ff56e5..cacdaa7a5e9e 100644 --- a/applications/GeoMechanicsApplication/custom_constitutive/thermal_law.h +++ b/applications/GeoMechanicsApplication/custom_constitutive/thermal_law.h @@ -27,8 +27,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoThermalLaw virtual ~GeoThermalLaw() = default; - virtual Matrix CalculateThermalDispersionMatrix(const Properties& rProp, - const ProcessInfo& rProcessInfo) const = 0; + [[nodiscard]] virtual Matrix CalculateThermalDispersionMatrix(const Properties& rProp) const = 0; private: friend class Serializer; diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp index 13ce812b0e00..615f096acb87 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp @@ -446,8 +446,7 @@ void UPwSmallStrainFICElement::CalculateAll(MatrixType& rLeftHa FICElementVariables FICVariables; this->InitializeFICElementVariables(FICVariables, Variables.DN_DXContainer, Geom, Prop, CurrentProcessInfo); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), CurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); const auto b_matrices = this->CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto integration_coefficients = @@ -462,6 +461,8 @@ void UPwSmallStrainFICElement::CalculateAll(MatrixType& rLeftHa strain_vectors, mStressVector, constitutive_matrices); const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients(constitutive_matrices, Prop); + const auto relative_permeability_values = this->CalculateRelativePermeabilityValues( + GeoTransportEquationUtilities::CalculateFluidPressures(Variables.NContainer, Variables.PressureVector)); for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { this->CalculateKinematics(Variables, GPoint); @@ -476,6 +477,7 @@ void UPwSmallStrainFICElement::CalculateAll(MatrixType& rLeftHa this->CalculateShapeFunctionsSecondOrderGradients(FICVariables, Variables); this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + Variables.RelativePermeability = relative_permeability_values[GPoint]; FICVariables.ShearModulus = CalculateShearModulus(Variables.ConstitutiveMatrix); diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp index 86dcda6cf6e1..bedd3bd6d209 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -149,8 +149,7 @@ void UPwSmallStrainElement::InitializeSolutionStep(const Proces ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - // Create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto deformation_gradients = CalculateDeformationGradients(); @@ -304,8 +303,7 @@ void UPwSmallStrainElement::FinalizeSolutionStep(const ProcessI ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - // Create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto deformation_gradients = CalculateDeformationGradients(); @@ -517,8 +515,7 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { // Compute Np, GradNpT, B and StrainVector @@ -641,12 +638,10 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints( deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, this->GetStressStatePolicy().GetVoigtSize()); - std::vector permeability_update_factors; - for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { - permeability_update_factors.push_back(this->CalculatePermeabilityUpdateFactor(strain_vectors[GPoint])); - } - - const auto fluid_fluxes = CalculateFluidFluxes(permeability_update_factors, rCurrentProcessInfo); + const auto fluid_fluxes = + CalculateFluidFluxes(GeoTransportEquationUtilities::CalculatePermeabilityUpdateFactors( + strain_vectors, this->GetProperties()), + rCurrentProcessInfo); for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { GeoElementUtilities::FillArray1dOutput(rOutput[GPoint], fluid_fluxes[GPoint]); } @@ -694,8 +689,7 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - // Create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); const Vector& VoigtVector = this->GetStressStatePolicy().GetVoigtVector(); @@ -920,8 +914,7 @@ void UPwSmallStrainElement::CalculateMassMatrix(MatrixType& rMa const auto N_container = r_geom.ShapeFunctionsValues(integration_method); const auto degrees_saturation = GeoTransportEquationUtilities::CalculateDegreesSaturation( - this->GetPressureSolutionVector(), N_container, mRetentionLawVector, this->GetProperties(), - rCurrentProcessInfo); + this->GetPressureSolutionVector(), N_container, mRetentionLawVector, this->GetProperties()); const auto solid_densities = GeoTransportEquationUtilities::CalculateSoilDensities(degrees_saturation, this->GetProperties()); @@ -968,8 +961,7 @@ void UPwSmallStrainElement::CalculateAll(MatrixType& rLe ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - // Create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto integration_coefficients = @@ -984,6 +976,13 @@ void UPwSmallStrainElement::CalculateAll(MatrixType& rLe strain_vectors, mStressVector, constitutive_matrices); const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients( constitutive_matrices, this->GetProperties()); + auto relative_permeability_values = this->CalculateRelativePermeabilityValues( + GeoTransportEquationUtilities::CalculateFluidPressures(Variables.NContainer, Variables.PressureVector)); + const auto permeability_update_factors = GeoTransportEquationUtilities::CalculatePermeabilityUpdateFactors( + strain_vectors, this->GetProperties()); + std::transform(relative_permeability_values.cbegin(), relative_permeability_values.cend(), + permeability_update_factors.cbegin(), relative_permeability_values.begin(), + std::multiplies{}); for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { this->CalculateKinematics(Variables, GPoint); @@ -997,13 +996,13 @@ void UPwSmallStrainElement::CalculateAll(MatrixType& rLe GeoElementUtilities::InterpolateVariableWithComponents( Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); - CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + Variables.RelativePermeability = relative_permeability_values[GPoint]; Variables.BiotCoefficient = biot_coefficients[GPoint]; Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse( Variables.BiotCoefficient, Variables.DegreeOfSaturation, Variables.DerivativeOfSaturation, this->GetProperties()); - Variables.PermeabilityUpdateFactor = this->CalculatePermeabilityUpdateFactor(Variables.StrainVector); Variables.IntegrationCoefficient = integration_coefficients[GPoint]; @@ -1034,54 +1033,31 @@ std::vector> UPwSmallStrainElement::Calc const PropertiesType& rProp = this->GetProperties(); - array_1d GradPressureTerm; - - // Create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + auto relative_permeability_values = this->CalculateRelativePermeabilityValues( + GeoTransportEquationUtilities::CalculateFluidPressures(Variables.NContainer, Variables.PressureVector)); + std::transform(relative_permeability_values.cbegin(), relative_permeability_values.cend(), + rPermeabilityUpdateFactors.cbegin(), relative_permeability_values.begin(), + std::multiplies<>{}); for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { this->CalculateKinematics(Variables, GPoint); - Variables.PermeabilityUpdateFactor = rPermeabilityUpdateFactors[GPoint]; GeoElementUtilities::InterpolateVariableWithComponents( Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); - RetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure( - Variables.Np, Variables.PressureVector)); - - Variables.RelativePermeability = - mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); + Variables.RelativePermeability = relative_permeability_values[GPoint]; - noalias(GradPressureTerm) = prod(trans(Variables.GradNpT), Variables.PressureVector); - noalias(GradPressureTerm) += PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration; + array_1d GradPressureTerm = prod(trans(Variables.GradNpT), Variables.PressureVector); + GradPressureTerm += PORE_PRESSURE_SIGN_FACTOR * rProp[DENSITY_WATER] * Variables.BodyAcceleration; FluidFluxes.push_back(PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse * - Variables.RelativePermeability * Variables.PermeabilityUpdateFactor * + Variables.RelativePermeability * prod(Variables.PermeabilityMatrix, GradPressureTerm)); } return FluidFluxes; } -template -double UPwSmallStrainElement::CalculatePermeabilityUpdateFactor(const Vector& rStrainVector) const -{ - KRATOS_TRY - - if (const auto& r_prop = this->GetProperties(); r_prop[PERMEABILITY_CHANGE_INVERSE_FACTOR] > 0.0) { - const double InverseCK = r_prop[PERMEABILITY_CHANGE_INVERSE_FACTOR]; - const double epsV = StressStrainUtilities::CalculateTrace(rStrainVector); - const double ePrevious = r_prop[POROSITY] / (1.0 - r_prop[POROSITY]); - const double eCurrent = (1.0 + ePrevious) * std::exp(epsV) - 1.0; - const double permLog10 = (eCurrent - ePrevious) * InverseCK; - return std::pow(10.0, permLog10); - } - - return 1.0; - - KRATOS_CATCH("") -} - template void UPwSmallStrainElement::InitializeElementVariables(ElementVariables& rVariables, const ProcessInfo& rCurrentProcessInfo) @@ -1172,7 +1148,12 @@ void UPwSmallStrainElement::CalculateAndAddLHS(MatrixType& rLef if (!rVariables.IgnoreUndrained) { this->CalculateAndAddCouplingMatrix(rLeftHandSideMatrix, rVariables); - this->CalculateAndAddPermeabilityMatrix(rLeftHandSideMatrix, rVariables); + + const auto permeability_matrix = + GeoTransportEquationUtilities::CalculatePermeabilityMatrix( + rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix, + rVariables.RelativePermeability, rVariables.IntegrationCoefficient); + GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, permeability_matrix); } KRATOS_CATCH("") @@ -1232,22 +1213,6 @@ void UPwSmallStrainElement::CalculateAndAddCompressibilityMatri KRATOS_CATCH("") } -template -void UPwSmallStrainElement::CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix, - const ElementVariables& rVariables) -{ - KRATOS_TRY - - const auto permeability_matrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix( - rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix, - rVariables.RelativePermeability, rVariables.PermeabilityUpdateFactor, rVariables.IntegrationCoefficient); - - // Distribute permeability block matrix into the elemental matrix - GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, permeability_matrix); - - KRATOS_CATCH("") -} - template void UPwSmallStrainElement::CalculateAndAddRHS(VectorType& rRightHandSideVector, ElementVariables& rVariables, @@ -1385,13 +1350,30 @@ void UPwSmallStrainElement::CalculatePermeabilityFlow( rPermeabilityMatrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix( rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix, - rVariables.RelativePermeability, rVariables.PermeabilityUpdateFactor, rVariables.IntegrationCoefficient); + rVariables.RelativePermeability, rVariables.IntegrationCoefficient); noalias(rPVector) = -prod(rPermeabilityMatrix, rVariables.PressureVector); KRATOS_CATCH("") } +template +std::vector UPwSmallStrainElement::CalculateRelativePermeabilityValues( + const std::vector& rFluidPressures) const +{ + KRATOS_ERROR_IF_NOT(rFluidPressures.size() == mRetentionLawVector.size()); + + auto retention_law_params = RetentionLaw::Parameters{this->GetProperties()}; + + auto result = std::vector{}; + std::transform(mRetentionLawVector.begin(), mRetentionLawVector.end(), rFluidPressures.begin(), + std::back_inserter(result), [&retention_law_params](auto pRetentionLaw, auto FluidPressure) { + retention_law_params.SetFluidPressure(FluidPressure); + return pRetentionLaw->CalculateRelativePermeability(retention_law_params); + }); + return result; +} + template void UPwSmallStrainElement::CalculateAndAddPermeabilityFlow(VectorType& rRightHandSideVector, ElementVariables& rVariables) @@ -1416,8 +1398,7 @@ void UPwSmallStrainElement::CalculateFluidBodyFlow(BoundedMatri noalias(rPDimMatrix) = prod(rVariables.GradNpT, rVariables.PermeabilityMatrix) * rVariables.IntegrationCoefficient; noalias(rPVector) = rVariables.DynamicViscosityInverse * rVariables.FluidDensity * - rVariables.RelativePermeability * rVariables.PermeabilityUpdateFactor * - prod(rPDimMatrix, rVariables.BodyAcceleration); + rVariables.RelativePermeability * prod(rPDimMatrix, rVariables.BodyAcceleration); KRATOS_CATCH("") } @@ -1543,8 +1524,6 @@ void UPwSmallStrainElement::InitializeProperties(ElementVariabl rVariables.Porosity = rProp[POROSITY]; GeoElementUtilities::FillPermeabilityMatrix(rVariables.PermeabilityMatrix, rProp); - rVariables.PermeabilityUpdateFactor = 1.0; - KRATOS_CATCH("") } @@ -1579,8 +1558,6 @@ void UPwSmallStrainElement::CalculateRetentionResponse(ElementV rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters); rVariables.DerivativeOfSaturation = mRetentionLawVector[GPoint]->CalculateDerivativeOfSaturation(rRetentionParameters); - rVariables.RelativePermeability = - mRetentionLawVector[GPoint]->CalculateRelativePermeability(rRetentionParameters); rVariables.BishopCoefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(rRetentionParameters); KRATOS_CATCH("") diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp index 9ffc5d2b333a..03a350c1e8f5 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp @@ -134,7 +134,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa double SolidDensity; double Density; double Porosity; - double PermeabilityUpdateFactor; double BiotCoefficient; double BiotModulusInverse; @@ -216,8 +215,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa virtual void CalculateKinematics(ElementVariables& rVariables, unsigned int PointNumber); - double CalculatePermeabilityUpdateFactor(const Vector& rStrainVector) const; - Matrix CalculateBMatrix(const Matrix& rDN_DX, const Vector& rN) const; std::vector CalculateBMatrices(const GeometryType::ShapeFunctionsGradientsType& rDN_DXContainer, const Matrix& rNContainer) const; @@ -230,9 +227,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa virtual void CalculateAndAddCompressibilityMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables); - virtual void CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix, - const ElementVariables& rVariables); - virtual void CalculateAndAddRHS(VectorType& rRightHandSideVector, ElementVariables& rVariables, unsigned int GPoint); void CalculateAndAddStiffnessForce(VectorType& rRightHandSideVector, @@ -249,8 +243,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa array_1d& rPVector, const ElementVariables& rVariables) const; + [[nodiscard]] std::vector CalculateRelativePermeabilityValues(const std::vector& rFluidPressures) const; virtual void CalculateAndAddPermeabilityFlow(VectorType& rRightHandSideVector, ElementVariables& rVariables); - virtual void CalculatePermeabilityFlow(BoundedMatrix& rPMatrix, array_1d& rPVector, const ElementVariables& rVariables) const; diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp index c0229a56a474..feda96d5a064 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp @@ -173,8 +173,7 @@ void UPwSmallStrainInterfaceElement::CalculateMassMatrix(Matrix InterfaceElementVariables Variables; this->InitializeElementVariables(Variables, Geom, Prop, rCurrentProcessInfo); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); // Defining necessary variables BoundedMatrix AuxDensityMatrix = ZeroMatrix(TDim, TNumNodes * TDim); @@ -253,8 +252,7 @@ void UPwSmallStrainInterfaceElement::InitializeSolutionStep(con unsigned int NumGPoints = mConstitutiveLawVector.size(); std::vector JointWidthContainer(NumGPoints); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { @@ -315,8 +313,7 @@ void UPwSmallStrainInterfaceElement::FinalizeSolutionStep(const unsigned int NumGPoints = mConstitutiveLawVector.size(); std::vector JointWidthContainer(NumGPoints); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { @@ -597,8 +594,8 @@ void UPwSmallStrainInterfaceElement::CalculateOnIntegrationPoin InterfaceElementVariables Variables; this->InitializeElementVariables(Variables, rGeom, this->GetProperties(), rCurrentProcessInfo); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); + // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { Variables.FluidPressure = GeoTransportEquationUtilities::CalculateFluidPressure( @@ -823,8 +820,7 @@ void UPwSmallStrainInterfaceElement::CalculateOnLobattoIntegrat InterfaceElementVariables Variables; this->InitializeElementVariables(Variables, Geom, Prop, rCurrentProcessInfo); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { @@ -990,8 +986,7 @@ void UPwSmallStrainInterfaceElement::CalculateOnLobattoIntegrat InterfaceElementVariables Variables; this->InitializeElementVariables(Variables, rGeom, rProp, rCurrentProcessInfo); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); Vector VoigtVector(mStressVector[0].size()); noalias(VoigtVector) = ZeroVector(VoigtVector.size()); @@ -1343,8 +1338,7 @@ void UPwSmallStrainInterfaceElement::CalculateAll(MatrixType& r array_1d RelDispVector; SFGradAuxVariables SFGradAuxVars; - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), CurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); const bool hasBiotCoefficient = Prop.Has(BIOT_COEFFICIENT); @@ -1382,7 +1376,7 @@ void UPwSmallStrainInterfaceElement::CalculateAll(MatrixType& r ModifyInactiveElementStress(Variables.JointWidth, mStressVector[GPoint]); - CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); this->InitializeBiotCoefficients(Variables, hasBiotCoefficient); @@ -1918,7 +1912,7 @@ void UPwSmallStrainInterfaceElement::CalculateAndAddPermeabilit rVariables.PPMatrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix( rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.LocalPermeabilityMatrix, - rVariables.RelativePermeability, rVariables.JointWidth, rVariables.IntegrationCoefficient); + rVariables.RelativePermeability * rVariables.JointWidth, rVariables.IntegrationCoefficient); // Distribute permeability block matrix into the elemental matrix GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, rVariables.PPMatrix); diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_link_interface_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_link_interface_element.cpp index 566ed00af83d..3d9871c71beb 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_link_interface_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_link_interface_element.cpp @@ -332,8 +332,7 @@ void UPwSmallStrainLinkInterfaceElement::CalculateAll(MatrixTyp array_1d RelDispVector; SFGradAuxVariables SFGradAuxVars; - // create general parametes of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), CurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); const bool hasBiotCoefficient = Prop.Has(BIOT_COEFFICIENT); diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.cpp index dd0383f72365..6bf60d74208f 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.cpp @@ -69,8 +69,7 @@ void UPwUpdatedLagrangianFICElement::CalculateAll(MatrixType& r FICElementVariables FICVariables; this->InitializeFICElementVariables(FICVariables, Variables.DN_DXContainer, Geom, Prop, rCurrentProcessInfo); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); const auto b_matrices = this->CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto integration_coefficients = @@ -85,6 +84,8 @@ void UPwUpdatedLagrangianFICElement::CalculateAll(MatrixType& r strain_vectors, mStressVector, constitutive_matrices); const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients(constitutive_matrices, Prop); + const auto relative_permeability_values = this->CalculateRelativePermeabilityValues( + GeoTransportEquationUtilities::CalculateFluidPressures(Variables.NContainer, Variables.PressureVector)); // Computing in all integrations points for (IndexType GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) { @@ -108,6 +109,7 @@ void UPwUpdatedLagrangianFICElement::CalculateAll(MatrixType& r this->CalculateShapeFunctionsSecondOrderGradients(FICVariables, Variables); this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + Variables.RelativePermeability = relative_permeability_values[GPoint]; // set shear modulus from stiffness matrix FICVariables.ShearModulus = CalculateShearModulus(Variables.ConstitutiveMatrix); diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp index fead3540251a..6e708af001bd 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp @@ -66,8 +66,7 @@ void UPwUpdatedLagrangianElement::CalculateAll(MatrixType& rLef ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); const auto b_matrices = this->CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto integration_coefficients = @@ -82,6 +81,8 @@ void UPwUpdatedLagrangianElement::CalculateAll(MatrixType& rLef strain_vectors, mStressVector, constitutive_matrices); const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients( constitutive_matrices, this->GetProperties()); + const auto relative_permeability_values = this->CalculateRelativePermeabilityValues( + GeoTransportEquationUtilities::CalculateFluidPressures(Variables.NContainer, Variables.PressureVector)); for (IndexType GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) { this->CalculateKinematics(Variables, GPoint); @@ -99,6 +100,7 @@ void UPwUpdatedLagrangianElement::CalculateAll(MatrixType& rLef Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + Variables.RelativePermeability = relative_permeability_values[GPoint]; Variables.BiotCoefficient = biot_coefficients[GPoint]; Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse( diff --git a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp index 75146b6718b5..fa2525fe2b2d 100644 --- a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp @@ -313,8 +313,7 @@ void SmallStrainUPwDiffOrderElement::InitializeSolutionStep(const ProcessInfo& r ConstitutiveParameters.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN); ConstitutiveParameters.Set(ConstitutiveLaw::INITIALIZE_MATERIAL_RESPONSE); // Note: this is for nonlocal damage - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(GetProperties()); const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto deformation_gradients = CalculateDeformationGradients(); @@ -455,7 +454,7 @@ void SmallStrainUPwDiffOrderElement::CalculateMassMatrix(MatrixType& rMassMatrix const PropertiesType& r_prop = this->GetProperties(); const auto degrees_saturation = GeoTransportEquationUtilities::CalculateDegreesSaturation( - this->GetPressureSolutionVector(), Np_container, mRetentionLawVector, r_prop, rCurrentProcessInfo); + this->GetPressureSolutionVector(), Np_container, mRetentionLawVector, r_prop); const auto solid_densities = GeoTransportEquationUtilities::CalculateSoilDensities(degrees_saturation, r_prop); @@ -531,8 +530,8 @@ void SmallStrainUPwDiffOrderElement::FinalizeSolutionStep(const ProcessInfo& rCu ConstitutiveLaw::Parameters ConstitutiveParameters(GetGeometry(), GetProperties(), rCurrentProcessInfo); ConstitutiveParameters.GetOptions().Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(GetProperties()); + const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto deformation_gradients = CalculateDeformationGradients(); const auto determinants_of_deformation_gradients = @@ -914,8 +913,7 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(GetProperties()); // Loop over integration points for (unsigned int GPoint = 0; GPoint < mRetentionLawVector.size(); ++GPoint) { @@ -1003,13 +1001,19 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(GetProperties(), rCurrentProcessInfo); const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto deformation_gradients = CalculateDeformationGradients(); const auto strain_vectors = StressStrainUtilities::CalculateStrains( deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain, mpStressStatePolicy->GetVoigtSize()); + auto relative_permeability_values = + CalculateRelativePermeabilityValues(GeoTransportEquationUtilities::CalculateFluidPressures( + Variables.NpContainer, Variables.PressureVector)); + const auto permeability_update_factors = + GeoTransportEquationUtilities::CalculatePermeabilityUpdateFactors(strain_vectors, GetProperties()); + std::transform(relative_permeability_values.cbegin(), relative_permeability_values.cend(), + permeability_update_factors.cbegin(), relative_permeability_values.begin(), + std::multiplies<>{}); // Loop over integration points for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { @@ -1028,16 +1032,11 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable BodyAcceleration[idim] += Variables.Nu[i] * Variables.BodyAcceleration[Index++]; } - RetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure( - Variables.Np, Variables.PressureVector)); - - const double RelativePermeability = - mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters); + const auto relative_permeability = relative_permeability_values[GPoint]; // Compute strain, need to update porosity Variables.F = deformation_gradients[GPoint]; Variables.StrainVector = strain_vectors[GPoint]; - Variables.PermeabilityUpdateFactor = this->CalculatePermeabilityUpdateFactor(Variables.StrainVector); Vector GradPressureTerm(Dim); noalias(GradPressureTerm) = prod(trans(Variables.DNp_DX), Variables.PressureVector); @@ -1046,8 +1045,7 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable Vector AuxFluidFlux = ZeroVector(Dim); AuxFluidFlux = PORE_PRESSURE_SIGN_FACTOR * Variables.DynamicViscosityInverse * - RelativePermeability * Variables.PermeabilityUpdateFactor * - prod(Variables.IntrinsicPermeability, GradPressureTerm); + relative_permeability * prod(Variables.IntrinsicPermeability, GradPressureTerm); Vector FluidFlux = ZeroVector(3); for (unsigned int idim = 0; idim < Dim; ++idim) @@ -1129,8 +1127,7 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable const Vector& VoigtVector = mpStressStatePolicy->GetVoigtVector(); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(GetProperties()); const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto deformation_gradients = CalculateDeformationGradients(); @@ -1255,8 +1252,7 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi if (CalculateResidualVectorFlag) ConstitutiveParameters.GetOptions().Set(ConstitutiveLaw::COMPUTE_STRESS); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(rProp, rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(rProp); // Loop over integration points const GeometryType::IntegrationPointsArrayType& IntegrationPoints = @@ -1275,6 +1271,14 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi strain_vectors, mStressVector, constitutive_matrices); const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients( constitutive_matrices, this->GetProperties()); + auto relative_permeability_values = + CalculateRelativePermeabilityValues(GeoTransportEquationUtilities::CalculateFluidPressures( + Variables.NpContainer, Variables.PressureVector)); + const auto permeability_update_factors = + GeoTransportEquationUtilities::CalculatePermeabilityUpdateFactors(strain_vectors, GetProperties()); + std::transform(relative_permeability_values.cbegin(), relative_permeability_values.cend(), + permeability_update_factors.cbegin(), relative_permeability_values.begin(), + std::multiplies<>{}); for (unsigned int GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) { this->CalculateKinematics(Variables, GPoint); @@ -1284,13 +1288,12 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi Variables.ConstitutiveMatrix = constitutive_matrices[GPoint]; CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + Variables.RelativePermeability = relative_permeability_values[GPoint]; Variables.BiotCoefficient = biot_coefficients[GPoint]; Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse( Variables.BiotCoefficient, Variables.DegreeOfSaturation, Variables.DerivativeOfSaturation, this->GetProperties()); - Variables.PermeabilityUpdateFactor = this->CalculatePermeabilityUpdateFactor(Variables.StrainVector); - Variables.IntegrationCoefficient = integration_coefficients[GPoint]; Variables.IntegrationCoefficientInitialConfiguration = this->CalculateIntegrationCoefficient( @@ -1443,9 +1446,6 @@ void SmallStrainUPwDiffOrderElement::InitializeElementVariables(ElementVariables rVariables.RelativePermeability = 1.0; rVariables.BishopCoefficient = 1.0; - // permeability change - rVariables.PermeabilityUpdateFactor = 1.0; - KRATOS_CATCH("") } @@ -1495,24 +1495,6 @@ void SmallStrainUPwDiffOrderElement::InitializeNodalVariables(ElementVariables& KRATOS_CATCH("") } -double SmallStrainUPwDiffOrderElement::CalculatePermeabilityUpdateFactor(const Vector& rStrainVector) const -{ - KRATOS_TRY - - if (const auto& r_prop = this->GetProperties(); r_prop[PERMEABILITY_CHANGE_INVERSE_FACTOR] > 0.0) { - const double InverseCK = r_prop[PERMEABILITY_CHANGE_INVERSE_FACTOR]; - const double epsV = StressStrainUtilities::CalculateTrace(rStrainVector); - const double ePrevious = r_prop[POROSITY] / (1.0 - r_prop[POROSITY]); - const double eCurrent = (1.0 + ePrevious) * std::exp(epsV) - 1.0; - const double permLog10 = (eCurrent - ePrevious) * InverseCK; - return std::pow(10.0, permLog10); - } - - return 1.0; - - KRATOS_CATCH("") -} - void SmallStrainUPwDiffOrderElement::InitializeProperties(ElementVariables& rVariables) { KRATOS_TRY @@ -1628,7 +1610,11 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddLHS(MatrixType& rLeftHandSid if (!rVariables.IgnoreUndrained) { this->CalculateAndAddCouplingMatrix(rLeftHandSideMatrix, rVariables); - this->CalculateAndAddPermeabilityMatrix(rLeftHandSideMatrix, rVariables); + + const auto permeability_matrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix( + rVariables.DNp_DX, rVariables.DynamicViscosityInverse, rVariables.IntrinsicPermeability, + rVariables.RelativePermeability, rVariables.IntegrationCoefficient); + GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, permeability_matrix); } KRATOS_CATCH("") @@ -1688,21 +1674,6 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddCompressibilityMatrix(Matrix KRATOS_CATCH("") } -void SmallStrainUPwDiffOrderElement::CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix, - const ElementVariables& rVariables) const -{ - KRATOS_TRY - - Matrix PermeabilityMatrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix( - rVariables.DNp_DX, rVariables.DynamicViscosityInverse, rVariables.IntrinsicPermeability, - rVariables.RelativePermeability, rVariables.PermeabilityUpdateFactor, rVariables.IntegrationCoefficient); - - // Distribute permeability block matrix into the elemental matrix - GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, PermeabilityMatrix); - - KRATOS_CATCH("") -} - void SmallStrainUPwDiffOrderElement::CalculateAndAddRHS(VectorType& rRightHandSideVector, ElementVariables& rVariables, unsigned int GPoint) @@ -1816,14 +1787,28 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddCompressibilityFlow(VectorTy KRATOS_CATCH("") } +std::vector SmallStrainUPwDiffOrderElement::CalculateRelativePermeabilityValues(const std::vector& rFluidPressures) const +{ + KRATOS_ERROR_IF_NOT(rFluidPressures.size() == mRetentionLawVector.size()); + + auto retention_law_params = RetentionLaw::Parameters{this->GetProperties()}; + + auto result = std::vector{}; + std::transform(mRetentionLawVector.begin(), mRetentionLawVector.end(), rFluidPressures.begin(), + std::back_inserter(result), [&retention_law_params](auto pRetentionLaw, auto FluidPressure) { + retention_law_params.SetFluidPressure(FluidPressure); + return pRetentionLaw->CalculateRelativePermeability(retention_law_params); + }); + return result; +} + void SmallStrainUPwDiffOrderElement::CalculateAndAddPermeabilityFlow(VectorType& rRightHandSideVector, ElementVariables& rVariables) const { KRATOS_TRY Matrix PermeabilityMatrix = - -PORE_PRESSURE_SIGN_FACTOR * rVariables.DynamicViscosityInverse * - rVariables.RelativePermeability * rVariables.PermeabilityUpdateFactor * + -PORE_PRESSURE_SIGN_FACTOR * rVariables.DynamicViscosityInverse * rVariables.RelativePermeability * prod(rVariables.DNp_DX, Matrix(prod(rVariables.IntrinsicPermeability, trans(rVariables.DNp_DX)))) * rVariables.IntegrationCoefficient; @@ -1840,10 +1825,9 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddFluidBodyFlow(VectorType& rR { KRATOS_TRY - Matrix GradNpTPerm = rVariables.DynamicViscosityInverse * GetProperties()[DENSITY_WATER] * - rVariables.RelativePermeability * rVariables.PermeabilityUpdateFactor * - prod(rVariables.DNp_DX, rVariables.IntrinsicPermeability) * - rVariables.IntegrationCoefficient; + Matrix GradNpTPerm = + rVariables.DynamicViscosityInverse * GetProperties()[DENSITY_WATER] * rVariables.RelativePermeability * + prod(rVariables.DNp_DX, rVariables.IntrinsicPermeability) * rVariables.IntegrationCoefficient; const GeometryType& rGeom = GetGeometry(); const SizeType Dim = rGeom.WorkingSpaceDimension(); @@ -1965,8 +1949,6 @@ void SmallStrainUPwDiffOrderElement::CalculateRetentionResponse(ElementVariables rVariables.DegreeOfSaturation = mRetentionLawVector[GPoint]->CalculateSaturation(rRetentionParameters); rVariables.DerivativeOfSaturation = mRetentionLawVector[GPoint]->CalculateDerivativeOfSaturation(rRetentionParameters); - rVariables.RelativePermeability = - mRetentionLawVector[GPoint]->CalculateRelativePermeability(rRetentionParameters); rVariables.BishopCoefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(rRetentionParameters); KRATOS_CATCH("") @@ -2031,4 +2013,4 @@ void SmallStrainUPwDiffOrderElement::CalculateAnyOfMaterialResponse( } } -} // Namespace Kratos \ No newline at end of file +} // Namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp index da783d4c8c3e..fd0c3e314539 100644 --- a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp @@ -204,7 +204,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub bool ConsiderGeometricStiffness; // stress/flow variables - double PermeabilityUpdateFactor; double BiotCoefficient; double BiotModulusInverse; double DynamicViscosityInverse; @@ -238,8 +237,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub void InitializeProperties(ElementVariables& rVariables); - double CalculatePermeabilityUpdateFactor(const Vector& rStrainVector) const; - virtual void CalculateKinematics(ElementVariables& rVariables, unsigned int GPoint); void CalculateDerivativesOnInitialConfiguration( @@ -258,8 +255,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub void CalculateAndAddCompressibilityMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables) const; - void CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix, const ElementVariables& rVariables) const; - void CalculateAndAddRHS(VectorType& rRightHandSideVector, ElementVariables& rVariables, unsigned int GPoint); void CalculateAndAddStiffnessForce(VectorType& rRightHandSideVector, @@ -273,6 +268,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub void CalculateAndAddCompressibilityFlow(VectorType& rRightHandSideVector, const ElementVariables& rVariables) const; + [[nodiscard]] std::vector CalculateRelativePermeabilityValues(const std::vector& rFluidPressures) const; void CalculateAndAddPermeabilityFlow(VectorType& rRightHandSideVector, ElementVariables& rVariables) const; void CalculateAndAddFluidBodyFlow(VectorType& rRightHandSideVector, ElementVariables& rVariables); diff --git a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.cpp b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.cpp index ef3f96f25813..13b4c251ab35 100644 --- a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.cpp @@ -12,6 +12,7 @@ // Application includes #include "custom_elements/steady_state_Pw_element.hpp" +#include "custom_utilities/transport_equation_utilities.hpp" namespace Kratos { @@ -152,12 +153,13 @@ void SteadyStatePwElement::CalculateAll(MatrixType& rLef ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); + const auto relative_permeability_values = this->CalculateRelativePermeabilityValues( + GeoTransportEquationUtilities::CalculateFluidPressures(Variables.NContainer, Variables.PressureVector)); const auto integration_coefficients = this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer); - + // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; GPoint++) { // Compute GradNpT, B and StrainVector @@ -169,6 +171,7 @@ void SteadyStatePwElement::CalculateAll(MatrixType& rLef Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + Variables.RelativePermeability = relative_permeability_values[GPoint]; Variables.IntegrationCoefficient = integration_coefficients[GPoint]; @@ -188,11 +191,10 @@ template void SteadyStatePwElement::CalculateAndAddLHS(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables) { - KRATOS_TRY; - - this->CalculateAndAddPermeabilityMatrix(rLeftHandSideMatrix, rVariables); - - KRATOS_CATCH(""); + const auto permeability_matrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix( + rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix, + rVariables.RelativePermeability, rVariables.IntegrationCoefficient); + rLeftHandSideMatrix += permeability_matrix; } //---------------------------------------------------------------------------------------- diff --git a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_interface_element.cpp b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_interface_element.cpp index 90d0e179e1fb..69e0201b16a9 100644 --- a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_interface_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_interface_element.cpp @@ -149,8 +149,7 @@ void SteadyStatePwInterfaceElement::CalculateAll(MatrixType& rL array_1d RelDispVector; SFGradAuxVariables SFGradAuxVars; - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), CurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); const auto integration_coefficients = this->CalculateIntegrationCoefficients(IntegrationPoints, detJContainer); @@ -171,7 +170,7 @@ void SteadyStatePwInterfaceElement::CalculateAll(MatrixType& rL InterfaceElementUtilities::FillPermeabilityMatrix( Variables.LocalPermeabilityMatrix, Variables.JointWidth, Prop[TRANSVERSAL_PERMEABILITY]); - CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); Variables.IntegrationCoefficient = integration_coefficients[GPoint]; diff --git a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_interface_element.hpp b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_interface_element.hpp index 03493024894c..c4f2d73490b5 100644 --- a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_interface_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_interface_element.hpp @@ -47,7 +47,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SteadyStatePwInterfaceElement /// The definition of the sizetype using SizeType = std::size_t; - using BaseType::CalculateRetentionResponse; using BaseType::mRetentionLawVector; using BaseType::mThisIntegrationMethod; diff --git a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.cpp b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.cpp index bf382761cf4f..196ede25f9c2 100644 --- a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.cpp @@ -211,8 +211,7 @@ void SteadyStatePwPipingElement::CalculateAll(MatrixType& rLeft array_1d RelDispVector; SFGradAuxVariables SFGradAuxVars; - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), CurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); const auto integration_coefficients = this->CalculateIntegrationCoefficients(IntegrationPoints, detJContainer); @@ -233,7 +232,7 @@ void SteadyStatePwPipingElement::CalculateAll(MatrixType& rLeft InterfaceElementUtilities::FillPermeabilityMatrix( Variables.LocalPermeabilityMatrix, Variables.JointWidth, Prop[TRANSVERSAL_PERMEABILITY]); - CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); Variables.IntegrationCoefficient = integration_coefficients[GPoint]; diff --git a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.hpp b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.hpp index 5248f3fd1caa..79ae0d90f880 100644 --- a/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_piping_element.hpp @@ -47,7 +47,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SteadyStatePwPipingElement /// The definition of the sizetype using SizeType = std::size_t; - using BaseType::CalculateRetentionResponse; using BaseType::mRetentionLawVector; using BaseType::mThisIntegrationMethod; diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp index 971d1bca4c21..e35855be6d4d 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp @@ -283,7 +283,7 @@ int TransientPwElement::Check(const ProcessInfo& rCurrentProces template void TransientPwElement::InitializeSolutionStep(const ProcessInfo& rCurrentProcessInfo) { - KRATOS_TRY; + KRATOS_TRY if (!mIsInitialised) this->Initialize(rCurrentProcessInfo); @@ -291,8 +291,7 @@ void TransientPwElement::InitializeSolutionStep(const ProcessIn const GeometryType& Geom = this->GetGeometry(); const unsigned int NumGPoints = Geom.IntegrationPointsNumber(this->GetIntegrationMethod()); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { @@ -303,30 +302,22 @@ void TransientPwElement::InitializeSolutionStep(const ProcessIn // reset hydraulic discharge this->ResetHydraulicDischarge(); - KRATOS_CATCH(""); + KRATOS_CATCH("") } //---------------------------------------------------------------------------------------------------- template -void TransientPwElement::InitializeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) +void TransientPwElement::InitializeNonLinearIteration(const ProcessInfo&) { - KRATOS_TRY; - // nothing - - KRATOS_CATCH(""); } //---------------------------------------------------------------------------------------------------- template -void TransientPwElement::FinalizeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) +void TransientPwElement::FinalizeNonLinearIteration(const ProcessInfo&) { - KRATOS_TRY; - // nothing - - KRATOS_CATCH(""); } //---------------------------------------------------------------------------------------------------- @@ -341,8 +332,7 @@ void TransientPwElement::FinalizeSolutionStep(const ProcessInfo const GeometryType& Geom = this->GetGeometry(); const unsigned int NumGPoints = Geom.IntegrationPointsNumber(this->GetIntegrationMethod()); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { @@ -450,9 +440,10 @@ void TransientPwElement::CalculateAll(MatrixType& rLeftH ElementVariables Variables; this->InitializeElementVariables(Variables, rCurrentProcessInfo); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); + const auto relative_permeability_values = this->CalculateRelativePermeabilityValues( + GeoTransportEquationUtilities::CalculateFluidPressures(Variables.NContainer, Variables.PressureVector)); const auto integration_coefficients = this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer); std::vector biot_coefficients(NumGPoints, Prop[BIOT_COEFFICIENT]); @@ -468,6 +459,7 @@ void TransientPwElement::CalculateAll(MatrixType& rLeftH Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + Variables.RelativePermeability = relative_permeability_values[GPoint]; Variables.BiotCoefficient = biot_coefficients[GPoint]; Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse( @@ -537,12 +529,16 @@ template void TransientPwElement::CalculateAndAddLHS(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables) { - KRATOS_TRY; + KRATOS_TRY this->CalculateAndAddCompressibilityMatrix(rLeftHandSideMatrix, rVariables); - this->CalculateAndAddPermeabilityMatrix(rLeftHandSideMatrix, rVariables); - KRATOS_CATCH(""); + const auto permeability_matrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix( + rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix, + rVariables.RelativePermeability, rVariables.IntegrationCoefficient); + rLeftHandSideMatrix += permeability_matrix; + + KRATOS_CATCH("") } //---------------------------------------------------------------------------------------- @@ -561,23 +557,6 @@ void TransientPwElement::CalculateAndAddCompressibilityMatrix(M KRATOS_CATCH(""); } -//---------------------------------------------------------------------------------------- -template -void TransientPwElement::CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix, - const ElementVariables& rVariables) -{ - KRATOS_TRY; - - const auto permeability_matrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix( - rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix, - rVariables.RelativePermeability, rVariables.PermeabilityUpdateFactor, rVariables.IntegrationCoefficient); - - // Distribute permeability block matrix into the elemental matrix - rLeftHandSideMatrix += permeability_matrix; - - KRATOS_CATCH(""); -} - //---------------------------------------------------------------------------------------- template void TransientPwElement::CalculateAndAddRHS(VectorType& rRightHandSideVector, diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.hpp b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.hpp index bd2c46485e79..233a40333cf3 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.hpp @@ -93,9 +93,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwElement : public UPwSmall void InitializeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override; - void InitializeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) override; + void InitializeNonLinearIteration(const ProcessInfo&) override; - void FinalizeNonLinearIteration(const ProcessInfo& rCurrentProcessInfo) override; + void FinalizeNonLinearIteration(const ProcessInfo&) override; void FinalizeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override; @@ -167,7 +167,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwElement : public UPwSmall void CalculateAndAddCompressibilityMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables) override; void CalculateAndAddPermeabilityFlow(VectorType& rRightHandSideVector, ElementVariables& rVariables) override; void CalculateAndAddFluidBodyFlow(VectorType& rRightHandSideVector, ElementVariables& rVariables) override; - void CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix, const ElementVariables& rVariables) override; void CalculateAndAddCompressibilityFlow(VectorType& rRightHandSideVector, ElementVariables& rVariables) override; unsigned int GetNumberOfDOF() const override; diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_interface_element.cpp b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_interface_element.cpp index 2fe776903e24..40aa31582bba 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_interface_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_interface_element.cpp @@ -143,12 +143,11 @@ void TransientPwInterfaceElement::CalculateMassMatrix(MatrixTyp } template -void TransientPwInterfaceElement::InitializeSolutionStep(const ProcessInfo& rCurrentProcessInfo) +void TransientPwInterfaceElement::InitializeSolutionStep(const ProcessInfo&) { KRATOS_TRY - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); // Loop over integration points for (unsigned int GPoint = 0; GPoint < mRetentionLawVector.size(); ++GPoint) { @@ -160,12 +159,11 @@ void TransientPwInterfaceElement::InitializeSolutionStep(const } template -void TransientPwInterfaceElement::FinalizeSolutionStep(const ProcessInfo& rCurrentProcessInfo) +void TransientPwInterfaceElement::FinalizeSolutionStep(const ProcessInfo&) { KRATOS_TRY - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); // Loop over integration points for (unsigned int GPoint = 0; GPoint < mRetentionLawVector.size(); ++GPoint) { @@ -295,8 +293,7 @@ void TransientPwInterfaceElement::CalculateOnLobattoIntegration // VG: Perhaps a new parameter to get join width and not minimum joint width const double& JointWidth = Prop[MINIMUM_JOINT_WIDTH]; - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { @@ -474,8 +471,7 @@ void TransientPwInterfaceElement::CalculateAll(MatrixType& rLef array_1d RelDispVector; SFGradAuxVariables SFGradAuxVars; - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), CurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); const bool hasBiotCoefficient = Prop.Has(BIOT_COEFFICIENT); @@ -498,7 +494,7 @@ void TransientPwInterfaceElement::CalculateAll(MatrixType& rLef InterfaceElementUtilities::FillPermeabilityMatrix( Variables.LocalPermeabilityMatrix, Variables.JointWidth, Prop[TRANSVERSAL_PERMEABILITY]); - CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); this->InitializeBiotCoefficients(Variables, hasBiotCoefficient); @@ -599,7 +595,7 @@ void TransientPwInterfaceElement::CalculateAndAddPermeabilityMa rVariables.PPMatrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix( rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.LocalPermeabilityMatrix, - rVariables.RelativePermeability, rVariables.JointWidth, rVariables.IntegrationCoefficient); + rVariables.RelativePermeability * rVariables.JointWidth, rVariables.IntegrationCoefficient); // Distribute permeability block matrix into the elemental matrix rLeftHandSideMatrix += rVariables.PPMatrix; diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_interface_element.hpp b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_interface_element.hpp index cb9060cb4bb7..2c7ca8298229 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_interface_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_interface_element.hpp @@ -44,7 +44,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwInterfaceElement /// The definition of the sizetype using SizeType = std::size_t; - using BaseType::CalculateRetentionResponse; using BaseType::mRetentionLawVector; using BaseType::mThisIntegrationMethod; @@ -109,8 +108,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwInterfaceElement void CalculateMassMatrix(MatrixType& rMassMatrix, const ProcessInfo& rCurrentProcessInfo) override; - void InitializeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override; - void FinalizeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override; + void InitializeSolutionStep(const ProcessInfo&) override; + void FinalizeSolutionStep(const ProcessInfo&) override; void CalculateOnIntegrationPoints(const Variable& rVariable, std::vector& rValues, diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.h b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.h index 429a232b6db1..156a67a0d0e3 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.h +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.h @@ -75,10 +75,10 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem const Matrix& r_N_container = GetGeometry().ShapeFunctionsValues(GetIntegrationMethod()); const auto integration_coefficients = CalculateIntegrationCoefficients(det_J_container); - const auto permeability_matrix = CalculatePermeabilityMatrix(dN_dX_container, integration_coefficients, rCurrentProcessInfo); - const auto compressibility_matrix = CalculateCompressibilityMatrix(r_N_container, integration_coefficients, rCurrentProcessInfo); + const auto permeability_matrix = CalculatePermeabilityMatrix(dN_dX_container, integration_coefficients); + const auto compressibility_matrix = CalculateCompressibilityMatrix(r_N_container, integration_coefficients); - const auto fluid_body_vector = CalculateFluidBodyVector(r_N_container, dN_dX_container, rCurrentProcessInfo, integration_coefficients); + const auto fluid_body_vector = CalculateFluidBodyVector(r_N_container, dN_dX_container, integration_coefficients); AddContributionsToLhsMatrix(rLeftHandSideMatrix, permeability_matrix, compressibility_matrix, rCurrentProcessInfo[DT_PRESSURE_COEFFICIENT]); AddContributionsToRhsVector(rRightHandSideVector, permeability_matrix, compressibility_matrix, fluid_body_vector); @@ -220,10 +220,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem BoundedMatrix CalculatePermeabilityMatrix( const GeometryType::ShapeFunctionsGradientsType& rShapeFunctionGradients, - const Vector& rIntegrationCoefficients, - const ProcessInfo& rCurrentProcessInfo) const + const Vector& rIntegrationCoefficients) const { - RetentionLaw::Parameters RetentionParameters(GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(GetProperties()); BoundedMatrix constitutive_matrix; const auto& r_properties = GetProperties(); GeoElementUtilities::FillPermeabilityMatrix(constitutive_matrix, r_properties); @@ -236,18 +235,16 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem double dynamic_viscosity_inverse = 1.0 / r_properties[DYNAMIC_VISCOSITY]; result += GeoTransportEquationUtilities::CalculatePermeabilityMatrix( rShapeFunctionGradients[integration_point_index], dynamic_viscosity_inverse, constitutive_matrix, - RelativePermeability, 1.0, rIntegrationCoefficients[integration_point_index]); + RelativePermeability, rIntegrationCoefficients[integration_point_index]); } return result; } - BoundedMatrix CalculateCompressibilityMatrix( - const Matrix& rNContainer, - const Vector& rIntegrationCoefficients, - const ProcessInfo& rCurrentProcessInfo) const + BoundedMatrix CalculateCompressibilityMatrix(const Matrix& rNContainer, + const Vector& rIntegrationCoefficients) const { const auto& r_properties = GetProperties(); - RetentionLaw::Parameters parameters(r_properties, rCurrentProcessInfo); + RetentionLaw::Parameters parameters(r_properties); auto retention_law = RetentionLawFactory::Clone(r_properties); auto result = BoundedMatrix{ZeroMatrix{TNumNodes, TNumNodes}}; @@ -255,7 +252,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem integration_point_index < GetGeometry().IntegrationPointsNumber(GetIntegrationMethod()); ++integration_point_index) { const auto N = Vector{row(rNContainer, integration_point_index)}; - const double BiotModulusInverse = CalculateBiotModulusInverse(rCurrentProcessInfo, integration_point_index); + const double BiotModulusInverse = CalculateBiotModulusInverse(integration_point_index); result += GeoTransportEquationUtilities::CalculateCompressibilityMatrix( N, BiotModulusInverse, rIntegrationCoefficients[integration_point_index]); } @@ -287,8 +284,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem } - double CalculateBiotModulusInverse(const ProcessInfo& rCurrentProcessInfo, - const unsigned int integrationPointIndex) const + double CalculateBiotModulusInverse(const unsigned int integrationPointIndex) const { const auto& r_properties = GetProperties(); const double biot_coefficient = r_properties[BIOT_COEFFICIENT]; @@ -300,7 +296,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem double result = (biot_coefficient - r_properties[POROSITY]) / r_properties[BULK_MODULUS_SOLID] + r_properties[POROSITY] / bulk_fluid; - RetentionLaw::Parameters RetentionParameters(GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(GetProperties()); const double degree_of_saturation = mRetentionLawVector[integrationPointIndex]->CalculateSaturation(RetentionParameters); const double derivative_of_saturation = mRetentionLawVector[integrationPointIndex]->CalculateDerivativeOfSaturation(RetentionParameters); @@ -310,18 +306,18 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem } - void InitializeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override + void InitializeSolutionStep(const ProcessInfo&) override { - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); for (auto retention_law : mRetentionLawVector) { retention_law->InitializeSolutionStep(RetentionParameters); } } - void FinalizeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override + void FinalizeSolutionStep(const ProcessInfo&) override { - RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(this->GetProperties()); for (auto retention_law : mRetentionLawVector) { retention_law->FinalizeSolutionStep(RetentionParameters); } @@ -331,7 +327,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem array_1d CalculateFluidBodyVector( const Matrix& rNContainer, const GeometryType::ShapeFunctionsGradientsType& rShapeFunctionGradients, - const ProcessInfo& rCurrentProcessInfo, const Vector& rIntegrationCoefficients) const { const std::size_t number_integration_points = GetGeometry().IntegrationPointsNumber(GetIntegrationMethod()); @@ -346,7 +341,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem BoundedMatrix constitutive_matrix; GeoElementUtilities::FillPermeabilityMatrix(constitutive_matrix, r_properties); - RetentionLaw::Parameters RetentionParameters(GetProperties(), rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(GetProperties()); array_1d volume_acceleration; GeoElementUtilities::GetNodalVariableVector(volume_acceleration, GetGeometry(), VOLUME_ACCELERATION); diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_thermal_element.h b/applications/GeoMechanicsApplication/custom_elements/transient_thermal_element.h index b65080a566bb..451dec406238 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_thermal_element.h +++ b/applications/GeoMechanicsApplication/custom_elements/transient_thermal_element.h @@ -88,8 +88,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientThermalElement : public Ele const auto integration_coefficients = CalculateIntegrationCoefficients(det_J_container); const auto conductivity_matrix = - CalculateConductivityMatrix(dN_dX_container, integration_coefficients, rCurrentProcessInfo); - const auto capacity_matrix = CalculateCapacityMatrix(integration_coefficients, rCurrentProcessInfo); + CalculateConductivityMatrix(dN_dX_container, integration_coefficients); + const auto capacity_matrix = CalculateCapacityMatrix(integration_coefficients); AddContributionsToLhsMatrix(rLeftHandSideMatrix, conductivity_matrix, capacity_matrix, rCurrentProcessInfo[DT_TEMPERATURE_COEFFICIENT]); @@ -266,11 +266,10 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientThermalElement : public Ele BoundedMatrix CalculateConductivityMatrix( const GeometryType::ShapeFunctionsGradientsType& rShapeFunctionGradients, - const Vector& rIntegrationCoefficients, - const ProcessInfo& rCurrentProcessInfo) const + const Vector& rIntegrationCoefficients) const { const auto law = CreateThermalLaw(); - const auto constitutive_matrix = law->CalculateThermalDispersionMatrix(GetProperties(), rCurrentProcessInfo); + const auto constitutive_matrix = law->CalculateThermalDispersionMatrix(GetProperties()); auto result = BoundedMatrix{ZeroMatrix{TNumNodes, TNumNodes}}; for (unsigned int integration_point_index = 0; @@ -285,11 +284,10 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientThermalElement : public Ele return result; } - BoundedMatrix CalculateCapacityMatrix(const Vector& rIntegrationCoefficients, - const ProcessInfo& rCurrentProcessInfo) const + BoundedMatrix CalculateCapacityMatrix(const Vector& rIntegrationCoefficients) const { const auto& r_properties = GetProperties(); - RetentionLaw::Parameters parameters(r_properties, rCurrentProcessInfo); + RetentionLaw::Parameters parameters(r_properties); auto retention_law = RetentionLawFactory::Clone(r_properties); const double saturation = retention_law->CalculateSaturation(parameters); const auto c_water = r_properties[POROSITY] * saturation * r_properties[DENSITY_WATER] * diff --git a/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.cpp b/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.cpp index a5b41c7fc351..fc7324687682 100644 --- a/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.cpp @@ -64,8 +64,7 @@ void UpdatedLagrangianUPwDiffOrderElement::CalculateAll(MatrixType& rLeft if (CalculateResidualVectorFlag) ConstitutiveParameters.GetOptions().Set(ConstitutiveLaw::COMPUTE_STRESS); - // create general parameters of retention law - RetentionLaw::Parameters RetentionParameters(rProp, rCurrentProcessInfo); + RetentionLaw::Parameters RetentionParameters(rProp); // Loop over integration points const GeometryType::IntegrationPointsArrayType& IntegrationPoints = @@ -84,6 +83,8 @@ void UpdatedLagrangianUPwDiffOrderElement::CalculateAll(MatrixType& rLeft strain_vectors, mStressVector, constitutive_matrices); const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients( constitutive_matrices, this->GetProperties()); + const auto relative_permeability_values = CalculateRelativePermeabilityValues( + GeoTransportEquationUtilities::CalculateFluidPressures(Variables.NpContainer, Variables.PressureVector)); for (IndexType GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) { this->CalculateKinematics(Variables, GPoint); @@ -94,6 +95,7 @@ void UpdatedLagrangianUPwDiffOrderElement::CalculateAll(MatrixType& rLeft Variables.ConstitutiveMatrix = constitutive_matrices[GPoint]; CalculateRetentionResponse(Variables, RetentionParameters, GPoint); + Variables.RelativePermeability = relative_permeability_values[GPoint]; Variables.BiotCoefficient = biot_coefficients[GPoint]; Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse( diff --git a/applications/GeoMechanicsApplication/custom_retention/retention_law.h b/applications/GeoMechanicsApplication/custom_retention/retention_law.h index d239aadcaa66..e4e87b992a16 100644 --- a/applications/GeoMechanicsApplication/custom_retention/retention_law.h +++ b/applications/GeoMechanicsApplication/custom_retention/retention_law.h @@ -67,8 +67,10 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw */ public: - Parameters(const Properties& rMaterialProperties, const ProcessInfo& rCurrentProcessInfo) - : mrCurrentProcessInfo(rCurrentProcessInfo), mrMaterialProperties(rMaterialProperties){}; + explicit Parameters(const Properties& rMaterialProperties) + : mrMaterialProperties(rMaterialProperties) + { + } ~Parameters() = default; @@ -82,13 +84,10 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLaw return mFluidPressure.value(); } - const ProcessInfo& GetProcessInfo() const { return mrCurrentProcessInfo; } - const Properties& GetMaterialProperties() const { return mrMaterialProperties; } private: std::optional mFluidPressure; - const ProcessInfo& mrCurrentProcessInfo; const Properties& mrMaterialProperties; }; // class Parameters end diff --git a/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp b/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp index 5b55ef8a9f1a..0b355334f9d0 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp +++ b/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp @@ -16,6 +16,7 @@ // Application includes #include "custom_retention/retention_law.h" +#include "custom_utilities/stress_strain_utilities.h" #include "geo_mechanics_application_variables.h" namespace Kratos @@ -30,24 +31,21 @@ class GeoTransportEquationUtilities double DynamicViscosityInverse, const BoundedMatrix& rMaterialPermeabilityMatrix, double RelativePermeability, - double PermeabilityUpdateFactor, double IntegrationCoefficient) { - return CalculatePermeabilityMatrix(rGradNpT, DynamicViscosityInverse, - rMaterialPermeabilityMatrix, RelativePermeability, - PermeabilityUpdateFactor, IntegrationCoefficient); + return CalculatePermeabilityMatrix(rGradNpT, DynamicViscosityInverse, rMaterialPermeabilityMatrix, + RelativePermeability, IntegrationCoefficient); } static inline Matrix CalculatePermeabilityMatrix(const Matrix& rGradNpT, double DynamicViscosityInverse, const Matrix& rMaterialPermeabilityMatrix, double RelativePermeability, - double PermeabilityUpdateFactor, double IntegrationCoefficient) { return -PORE_PRESSURE_SIGN_FACTOR * DynamicViscosityInverse * prod(rGradNpT, Matrix(prod(rMaterialPermeabilityMatrix, trans(rGradNpT)))) * - RelativePermeability * PermeabilityUpdateFactor * IntegrationCoefficient; + RelativePermeability * IntegrationCoefficient; } template @@ -100,10 +98,9 @@ class GeoTransportEquationUtilities static Vector CalculateDegreesSaturation(const Vector& rPressureSolution, const Matrix& rNContainer, const std::vector& rRetentionLawVector, - const Properties& rProp, - const ProcessInfo& rCurrentProcessInfo) + const Properties& rProp) { - RetentionLaw::Parameters retention_parameters(rProp, rCurrentProcessInfo); + RetentionLaw::Parameters retention_parameters(rProp); Vector result(rRetentionLawVector.size()); for (unsigned int g_point = 0; g_point < rRetentionLawVector.size(); ++g_point) { const double fluid_pressure = CalculateFluidPressure(row(rNContainer, g_point), rPressureSolution); @@ -118,6 +115,15 @@ class GeoTransportEquationUtilities return inner_prod(rN, rPressureVector); } + [[nodiscard]] static std::vector CalculateFluidPressures(const Matrix& rNContainer, const Vector& rPressureVector) + { + auto result = std::vector{}; + for (auto i = std::size_t{0}; i < rNContainer.size1(); ++i) { + result.emplace_back(CalculateFluidPressure(row(rNContainer, i), rPressureVector)); + } + return result; + } + [[nodiscard]] static double CalculateBiotModulusInverse(double BiotCoefficient, double DegreeOfSaturation, double DerivativeOfSaturation, @@ -153,6 +159,17 @@ class GeoTransportEquationUtilities return result; } + [[nodiscard]] static std::vector CalculatePermeabilityUpdateFactors(const std::vector& rStrainVectors, + const Properties& rProperties) + { + auto result = std::vector{}; + std::transform(rStrainVectors.cbegin(), rStrainVectors.cend(), std::back_inserter(result), + [&rProperties](const auto& rStrainVector) { + return CalculatePermeabilityUpdateFactor(rStrainVector, rProperties); + }); + return result; + } + private: [[nodiscard]] static double CalculateBiotCoefficient(const Matrix& rConstitutiveMatrix, const Properties& rProperties) @@ -162,5 +179,19 @@ class GeoTransportEquationUtilities : 1.0 - CalculateBulkModulus(rConstitutiveMatrix) / rProperties[BULK_MODULUS_SOLID]; } + [[nodiscard]] static double CalculatePermeabilityUpdateFactor(const Vector& rStrainVector, + const Properties& rProperties) + { + if (rProperties[PERMEABILITY_CHANGE_INVERSE_FACTOR] > 0.0) { + const double InverseCK = rProperties[PERMEABILITY_CHANGE_INVERSE_FACTOR]; + const double epsV = StressStrainUtilities::CalculateTrace(rStrainVector); + const double ePrevious = rProperties[POROSITY] / (1.0 - rProperties[POROSITY]); + const double eCurrent = (1.0 + ePrevious) * std::exp(epsV) - 1.0; + const double permLog10 = (eCurrent - ePrevious) * InverseCK; + return std::pow(10.0, permLog10); + } + + return 1.0; + } }; /* Class GeoTransportEquationUtilities*/ } /* namespace Kratos.*/ diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_permeability_matrix.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_permeability_matrix.cpp index 10e3860afe44..1463faa36023 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/test_permeability_matrix.cpp +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_permeability_matrix.cpp @@ -41,8 +41,8 @@ KRATOS_TEST_CASE_IN_SUITE(CalculatePermeabilityMatrix2D3NGivesCorrectResults, Kr const double RelativePermeability = 0.02; const double PermeabilityUpdateFactor = 1.5; PermeabilityMatrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix<2, 3>( - GradNpT, DynamicViscosityInverse, MaterialPermeabilityMatrix, RelativePermeability, - PermeabilityUpdateFactor, IntegrationCoefficient); + GradNpT, DynamicViscosityInverse, MaterialPermeabilityMatrix, + RelativePermeability * PermeabilityUpdateFactor, IntegrationCoefficient); BoundedMatrix PMatrix; // clang-format off @@ -76,8 +76,8 @@ KRATOS_TEST_CASE_IN_SUITE(CalculatePermeabilityMatrix3D4NGivesCorrectResults, Kr const double RelativePermeability = 0.1; const double PermeabilityUpdateFactor = 2.0; PermeabilityMatrix = GeoTransportEquationUtilities::CalculatePermeabilityMatrix<3, 4>( - GradNpT, DynamicViscosityInverse, MaterialPermeabilityMatrix, RelativePermeability, - PermeabilityUpdateFactor, IntegrationCoefficient); + GradNpT, DynamicViscosityInverse, MaterialPermeabilityMatrix, + RelativePermeability * PermeabilityUpdateFactor, IntegrationCoefficient); BoundedMatrix PMatrix; // clang-format off diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_thermal_dispersion_law.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_thermal_dispersion_law.cpp index 8067261ed9b4..fc6516a1f782 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/test_thermal_dispersion_law.cpp +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_thermal_dispersion_law.cpp @@ -33,11 +33,9 @@ KRATOS_TEST_CASE_IN_SUITE(CalculateThermalDispersionMatrix2D, KratosGeoMechanics const SizeType dimension = 2; GeoThermalDispersionLaw geo_thermal_dispersion_2D_law(dimension); - ProcessInfo info; const Matrix thermal_dispersion_matrix = - geo_thermal_dispersion_2D_law.CalculateThermalDispersionMatrix( - *cond_prop, info); + geo_thermal_dispersion_2D_law.CalculateThermalDispersionMatrix(*cond_prop); Matrix expected_solution = ZeroMatrix(2, 2); expected_solution(0, 0) = 1125.0; @@ -74,11 +72,9 @@ KRATOS_TEST_CASE_IN_SUITE(CalculateThermalDispersionMatrix3D, KratosGeoMechanics const SizeType dimension = 3; GeoThermalDispersionLaw geo_thermal_dispersion_3D_law(dimension); - ProcessInfo info; const Matrix thermal_dispersion_matrix = - geo_thermal_dispersion_3D_law.CalculateThermalDispersionMatrix( - *cond_prop, info); + geo_thermal_dispersion_3D_law.CalculateThermalDispersionMatrix(*cond_prop); Matrix expected_solution = ZeroMatrix(3, 3); expected_solution(0, 0) = 800.0; diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_thermal_filter_law.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_thermal_filter_law.cpp index 105c4c224c3e..1d745b3804e8 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/test_thermal_filter_law.cpp +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_thermal_filter_law.cpp @@ -27,9 +27,8 @@ KRATOS_TEST_CASE_IN_SUITE(CalculateThermalFilterLawMatrix, KratosGeoMechanicsFas p_cond_prop->SetValue(THERMAL_CONDUCTIVITY_WATER, 1000.0); GeoThermalFilterLaw geo_thermal_filter_law; - ProcessInfo info; - const Matrix thermal_filter_matrix = geo_thermal_filter_law.CalculateThermalDispersionMatrix(*p_cond_prop, info); + const Matrix thermal_filter_matrix = geo_thermal_filter_law.CalculateThermalDispersionMatrix(*p_cond_prop); Matrix expected_solution = ScalarMatrix(1, 1, 1000.0); diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_transport_equation_utilities.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_transport_equation_utilities.cpp index 0944c29bde02..e5e0fc372e9e 100644 --- a/applications/GeoMechanicsApplication/tests/cpp_tests/test_transport_equation_utilities.cpp +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_transport_equation_utilities.cpp @@ -14,6 +14,13 @@ #include "testing/testing.h" #include +namespace +{ + +constexpr auto tolerance = 1.0e-12; + +} + namespace Kratos::Testing { @@ -134,4 +141,67 @@ KRATOS_TEST_CASE_IN_SUITE(CalculateBiotCoefficients_GivesInfResults_ForZeroBulkM [](const double value) { return std::isinf(value); })) } +KRATOS_TEST_CASE_IN_SUITE(EachFluidPressureIsTheInnerProductOfShapeFunctionsAndPressure, KratosGeoMechanicsFastSuite) +{ + auto shape_function_values = Matrix{3, 3, 0.0}; + // clang-format off + shape_function_values <<= 0.8, 0.2, 0.3, + 0.1, 0.7, 0.4, + 0.2, 0.5, 0.6; + // clang-format on + + auto pore_water_pressures = Vector(3); + pore_water_pressures <<= 2.0, 3.0, 4.0; + + const auto expected_fluid_pressures = std::vector{3.4, 3.9, 4.3}; + KRATOS_EXPECT_VECTOR_NEAR(GeoTransportEquationUtilities::CalculateFluidPressures( + shape_function_values, pore_water_pressures), + expected_fluid_pressures, tolerance) +} + +KRATOS_TEST_CASE_IN_SUITE(PermeabilityUpdateFactorEqualsOneWhenChangeInverseFactorIsNotGiven, KratosGeoMechanicsFastSuite) +{ + const auto unused_strain_vectors = std::vector(3, Vector{}); + auto properties = Properties{}; + + const auto expected_factors = std::vector(unused_strain_vectors.size(), 1.0); + KRATOS_EXPECT_VECTOR_NEAR(GeoTransportEquationUtilities::CalculatePermeabilityUpdateFactors( + unused_strain_vectors, properties), + expected_factors, tolerance) +} + +KRATOS_TEST_CASE_IN_SUITE(PermeabilityUpdateFactorEqualsOneWhenChangeInverseFactorIsNonPositive, KratosGeoMechanicsFastSuite) +{ + const auto unused_strain_vectors = std::vector(3, Vector{}); + auto properties = Properties{}; + properties[PERMEABILITY_CHANGE_INVERSE_FACTOR] = -1.0; + + const auto expected_factors = std::vector(unused_strain_vectors.size(), 1.0); + KRATOS_EXPECT_VECTOR_NEAR(GeoTransportEquationUtilities::CalculatePermeabilityUpdateFactors( + unused_strain_vectors, properties), + expected_factors, tolerance) + + properties[PERMEABILITY_CHANGE_INVERSE_FACTOR] = 0.0; + + KRATOS_EXPECT_VECTOR_NEAR(GeoTransportEquationUtilities::CalculatePermeabilityUpdateFactors( + unused_strain_vectors, properties), + expected_factors, tolerance) +} + +KRATOS_TEST_CASE_IN_SUITE(PermeabilityUpdateFactorIsComputedFromStrainsAndPropertiesWhenChangeInverseFactorIsPositive, KratosGeoMechanicsFastSuite) +{ + auto test_strains = Vector{3}; + test_strains <<= 0.001, 0.002, 0.0; + auto strain_vectors = std::vector{test_strains, 2.0 * test_strains, 4.0 * test_strains}; + + auto properties = Properties{}; + properties[PERMEABILITY_CHANGE_INVERSE_FACTOR] = 0.5; + properties[POROSITY] = 0.2; + + const auto expected_factors = std::vector{1.00433, 1.0087, 1.01753}; + KRATOS_EXPECT_VECTOR_NEAR(GeoTransportEquationUtilities::CalculatePermeabilityUpdateFactors( + strain_vectors, properties), + expected_factors, 1e-5) +} + } // namespace Kratos::Testing \ No newline at end of file