Skip to content

Commit abdbdec

Browse files
authored
[GeoMechanicsApplication] implement the "Undrained" behaviour, for the constitutive tensor (#14372)
* corrected error message * taken GetSkemptonB, GetUndrainedPoissonsRatio, GetUndrainedYoungsModulus with unit tests * added MakeContinuumConstitutiveTensor and unit test * added checks for denominators * changed variable names * renamed CalculateElasticMatrix as CalculateElasticConstitutiveTensor * added GEO_POISSON_UNDRAINED and GEO_SKEMPTON_B * renamed MakeInterfaceConstitutiveMatrix as MakeInterfaceConstitutiveTensor * used MakeContinuumConstitutiveTensor * added unit tests for MakeInterfaceConstitutiveTensor * added ConstitutiveLawUtitlities_CalculateK0NCFromFrictionAngleGivesK0NC * added unit tests that throw errors * minor for future * minor changes in variable names * more tests for GetUndrainedPoissonsRatio * added check and tests for GetSkemptonB * add warning in GetUndrainedPoissonsRatio * added comments and minor changes * fixed a code smell * renamed GetUndrainedYoungsModulus to CalculateUndrainedYoungsModulus * response to the review comments * fixed unused variables * used hard-coded expected values
1 parent a030ba3 commit abdbdec

24 files changed

Lines changed: 428 additions & 93 deletions

applications/GeoMechanicsApplication/custom_constitutive/constitutive_law_dimension.h

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -28,12 +28,12 @@ class ConstitutiveLawDimension
2828
public:
2929
virtual ~ConstitutiveLawDimension() = default;
3030

31-
[[nodiscard]] virtual Matrix CalculateElasticMatrix(const Properties& rProperties) const = 0;
32-
[[nodiscard]] virtual std::unique_ptr<ConstitutiveLawDimension> Clone() const = 0;
33-
[[nodiscard]] virtual std::size_t GetStrainSize() const = 0;
34-
[[nodiscard]] virtual std::size_t GetDimension() const = 0;
35-
[[nodiscard]] virtual std::size_t GetNumberOfNormalComponents() const = 0;
36-
[[nodiscard]] virtual Flags GetSpatialType() const = 0;
31+
[[nodiscard]] virtual Matrix CalculateElasticConstitutiveTensor(const Properties& rProperties) const = 0;
32+
[[nodiscard]] virtual std::unique_ptr<ConstitutiveLawDimension> Clone() const = 0;
33+
[[nodiscard]] virtual std::size_t GetStrainSize() const = 0;
34+
[[nodiscard]] virtual std::size_t GetDimension() const = 0;
35+
[[nodiscard]] virtual std::size_t GetNumberOfNormalComponents() const = 0;
36+
[[nodiscard]] virtual Flags GetSpatialType() const = 0;
3737

3838
private:
3939
friend class Serializer;

applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.cpp

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -59,7 +59,8 @@ Matrix& GeoIncrementalLinearElasticInterfaceLaw::CalculateValue(ConstitutiveLaw:
5959
Matrix& rValue)
6060
{
6161
if (rThisVariable == CONSTITUTIVE_MATRIX) {
62-
rValue = mpConstitutiveLawDimension->CalculateElasticMatrix(rParameterValues.GetMaterialProperties());
62+
rValue = mpConstitutiveLawDimension->CalculateElasticConstitutiveTensor(
63+
rParameterValues.GetMaterialProperties());
6364
} else {
6465
KRATOS_ERROR << "Can't calculate value of " << rThisVariable.Name() << ": unsupported variable\n";
6566
}
@@ -88,7 +89,7 @@ void GeoIncrementalLinearElasticInterfaceLaw::CalculateMaterialResponseCauchy(Co
8889
{
8990
rValues.GetStressVector() =
9091
mPreviousTraction +
91-
prod(mpConstitutiveLawDimension->CalculateElasticMatrix(rValues.GetMaterialProperties()),
92+
prod(mpConstitutiveLawDimension->CalculateElasticConstitutiveTensor(rValues.GetMaterialProperties()),
9293
rValues.GetStrainVector() - mPreviousRelativeDisplacement);
9394
}
9495

applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_law.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -125,7 +125,7 @@ void GeoIncrementalLinearElasticLaw::CalculateElasticMatrix(Matrix& C, Constitut
125125
{
126126
KRATOS_TRY
127127

128-
C = mpConstitutiveDimension->CalculateElasticMatrix(rValues.GetMaterialProperties());
128+
C = mpConstitutiveDimension->CalculateElasticConstitutiveTensor(rValues.GetMaterialProperties());
129129

130130
if (this->GetConsiderDiagonalEntriesOnlyAndNoShear()) {
131131
SetEntriesAboveDiagonalToZero(C);

applications/GeoMechanicsApplication/custom_constitutive/interface_coulomb_with_tension_cut_off.cpp

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -161,7 +161,7 @@ void InterfaceCoulombWithTensionCutOff::CalculateMaterialResponseCauchy(Paramete
161161

162162
if (!mCoulombWithTensionCutOffImpl.IsAdmissibleStressState(trial_sigma_tau)) {
163163
mapped_sigma_tau = mCoulombWithTensionCutOffImpl.DoReturnMapping(
164-
trial_sigma_tau, mpConstitutiveDimension->CalculateElasticMatrix(r_properties),
164+
trial_sigma_tau, mpConstitutiveDimension->CalculateElasticConstitutiveTensor(r_properties),
165165
Geo::PrincipalStresses::AveragingType::NO_AVERAGING);
166166
if (negative) mapped_sigma_tau.Tau() *= -1.0;
167167
}
@@ -176,7 +176,7 @@ Geo::SigmaTau InterfaceCoulombWithTensionCutOff::CalculateTrialTractionVector(co
176176
{
177177
constexpr auto number_of_normal_components = std::size_t{1};
178178
return Geo::SigmaTau{mTractionVectorFinalized +
179-
prod(ConstitutiveLawUtilities::MakeInterfaceConstitutiveMatrix(
179+
prod(ConstitutiveLawUtilities::MakeInterfaceElasticConstitutiveTensor(
180180
NormalStiffness, ShearStiffness, GetStrainSize(), number_of_normal_components),
181181
rRelativeDisplacementVector - mRelativeDisplacementVectorFinalized)};
182182
}
@@ -194,7 +194,7 @@ Matrix& InterfaceCoulombWithTensionCutOff::CalculateValue(Parameters& rConstitut
194194
if (rVariable == CONSTITUTIVE_MATRIX) {
195195
const auto& r_properties = rConstitutiveLawParameters.GetMaterialProperties();
196196
constexpr auto number_of_normal_components = std::size_t{1};
197-
rValue = ConstitutiveLawUtilities::MakeInterfaceConstitutiveMatrix(
197+
rValue = ConstitutiveLawUtilities::MakeInterfaceElasticConstitutiveTensor(
198198
r_properties[INTERFACE_NORMAL_STIFFNESS], r_properties[INTERFACE_SHEAR_STIFFNESS],
199199
GetStrainSize(), number_of_normal_components);
200200
} else {

applications/GeoMechanicsApplication/custom_constitutive/interface_plane_strain.cpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -20,9 +20,9 @@
2020
namespace Kratos
2121
{
2222

23-
Matrix InterfacePlaneStrain::CalculateElasticMatrix(const Properties& rProperties) const
23+
Matrix InterfacePlaneStrain::CalculateElasticConstitutiveTensor(const Properties& rProperties) const
2424
{
25-
return ConstitutiveLawUtilities::MakeInterfaceConstitutiveMatrix(
25+
return ConstitutiveLawUtilities::MakeInterfaceElasticConstitutiveTensor(
2626
rProperties[INTERFACE_NORMAL_STIFFNESS], rProperties[INTERFACE_SHEAR_STIFFNESS],
2727
GetStrainSize(), GetNumberOfNormalComponents());
2828
}

applications/GeoMechanicsApplication/custom_constitutive/interface_plane_strain.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@ namespace Kratos
2121
class KRATOS_API(GEO_MECHANICS_APPLICATION) InterfacePlaneStrain : public ConstitutiveLawDimension
2222
{
2323
public:
24-
[[nodiscard]] Matrix CalculateElasticMatrix(const Properties& rProperties) const override;
24+
[[nodiscard]] Matrix CalculateElasticConstitutiveTensor(const Properties& rProperties) const override;
2525
[[nodiscard]] std::unique_ptr<ConstitutiveLawDimension> Clone() const override;
2626
[[nodiscard]] std::size_t GetStrainSize() const override;
2727
[[nodiscard]] std::size_t GetDimension() const override;

applications/GeoMechanicsApplication/custom_constitutive/interface_three_dimensional_surface.cpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -20,9 +20,9 @@
2020
namespace Kratos
2121
{
2222

23-
Matrix InterfaceThreeDimensionalSurface::CalculateElasticMatrix(const Properties& rProperties) const
23+
Matrix InterfaceThreeDimensionalSurface::CalculateElasticConstitutiveTensor(const Properties& rProperties) const
2424
{
25-
auto result = ConstitutiveLawUtilities::MakeInterfaceConstitutiveMatrix(
25+
auto result = ConstitutiveLawUtilities::MakeInterfaceElasticConstitutiveTensor(
2626
rProperties[INTERFACE_NORMAL_STIFFNESS], rProperties[INTERFACE_SHEAR_STIFFNESS],
2727
GetStrainSize(), GetNumberOfNormalComponents());
2828

applications/GeoMechanicsApplication/custom_constitutive/interface_three_dimensional_surface.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@ namespace Kratos
2121
class KRATOS_API(GEO_MECHANICS_APPLICATION) InterfaceThreeDimensionalSurface : public ConstitutiveLawDimension
2222
{
2323
public:
24-
[[nodiscard]] Matrix CalculateElasticMatrix(const Properties& rProperties) const override;
24+
[[nodiscard]] Matrix CalculateElasticConstitutiveTensor(const Properties& rProperties) const override;
2525
[[nodiscard]] std::unique_ptr<ConstitutiveLawDimension> Clone() const override;
2626
[[nodiscard]] std::size_t GetStrainSize() const override;
2727
[[nodiscard]] std::size_t GetDimension() const override;

applications/GeoMechanicsApplication/custom_constitutive/linear_elastic_law.cpp

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -158,8 +158,7 @@ int GeoLinearElasticLaw::Check(const Properties& rMaterialProperties,
158158
check_properties.Check(YOUNG_MODULUS);
159159
constexpr auto min_value_poisson_ratio = -1.0;
160160
constexpr auto max_value_poisson_ratio = 0.5;
161-
check_properties.SingleUseBounds(CheckProperties::Bounds::InclusiveLowerAndExclusiveUpper)
162-
.Check(POISSON_RATIO, min_value_poisson_ratio, max_value_poisson_ratio);
161+
check_properties.Check(POISSON_RATIO, min_value_poisson_ratio, max_value_poisson_ratio);
163162

164163
return 0;
165164
}

applications/GeoMechanicsApplication/custom_constitutive/mohr_coulomb_with_tension_cutoff.cpp

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -200,7 +200,8 @@ void MohrCoulombWithTensionCutOff::CalculateMaterialResponseCauchy(ConstitutiveL
200200
const auto& r_properties = rParameters.GetMaterialProperties();
201201

202202
if (rParameters.GetOptions().Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) {
203-
rParameters.GetConstitutiveMatrix() = mpConstitutiveDimension->CalculateElasticMatrix(r_properties);
203+
rParameters.GetConstitutiveMatrix() =
204+
mpConstitutiveDimension->CalculateElasticConstitutiveTensor(r_properties);
204205
}
205206
if (!rParameters.GetOptions().Is(ConstitutiveLaw::COMPUTE_STRESS)) {
206207
return;
@@ -215,7 +216,7 @@ void MohrCoulombWithTensionCutOff::CalculateMaterialResponseCauchy(ConstitutiveL
215216
} else {
216217
mCoulombWithTensionCutOffImpl.SaveKappaOfCoulombYieldSurface();
217218
auto mapped_principal_stresses = mCoulombWithTensionCutOffImpl.DoReturnMapping(
218-
trial_principal_stresses, mpConstitutiveDimension->CalculateElasticMatrix(r_properties),
219+
trial_principal_stresses, mpConstitutiveDimension->CalculateElasticConstitutiveTensor(r_properties),
219220
Geo::PrincipalStresses::AveragingType::NO_AVERAGING);
220221

221222
// For interchanging principal stresses, retry mapping with averaged principal stresses.
@@ -226,7 +227,7 @@ void MohrCoulombWithTensionCutOff::CalculateMaterialResponseCauchy(ConstitutiveL
226227
mCoulombWithTensionCutOffImpl.RestoreKappaOfCoulombYieldSurface();
227228
mapped_principal_stresses = mCoulombWithTensionCutOffImpl.DoReturnMapping(
228229
averaged_principal_trial_stress_vector,
229-
mpConstitutiveDimension->CalculateElasticMatrix(r_properties), averaging_type);
230+
mpConstitutiveDimension->CalculateElasticConstitutiveTensor(r_properties), averaging_type);
230231
mapped_principal_stresses.Values()[1] =
231232
mapped_principal_stresses.Values()[AveragingTypeToArrayIndex(averaging_type)];
232233
}
@@ -240,7 +241,7 @@ void MohrCoulombWithTensionCutOff::CalculateMaterialResponseCauchy(ConstitutiveL
240241
Vector MohrCoulombWithTensionCutOff::CalculateTrialStressVector(const Vector& rStrainVector,
241242
const Properties& rProperties) const
242243
{
243-
return mStressVectorFinalized + prod(mpConstitutiveDimension->CalculateElasticMatrix(rProperties),
244+
return mStressVectorFinalized + prod(mpConstitutiveDimension->CalculateElasticConstitutiveTensor(rProperties),
244245
rStrainVector - mStrainVectorFinalized);
245246
}
246247

0 commit comments

Comments
 (0)