From 0d75919cda04990ccf534aa3a151f610b8f03c17 Mon Sep 17 00:00:00 2001 From: Fan Fei Date: Fri, 3 Oct 2025 12:35:49 -0700 Subject: [PATCH 1/4] initial commit to enable anelastic strain --- .../constitutive/solid/CeramicDamage.hpp | 17 ++++-- .../constitutive/solid/DelftEgg.hpp | 16 ++++-- .../constitutive/solid/DruckerPrager.hpp | 17 ++++-- .../solid/DruckerPragerExtended.hpp | 16 ++++-- .../constitutive/solid/ElasticIsotropic.hpp | 48 ++++++++++++++-- .../ElasticIsotropicPressureDependent.hpp | 20 +++++-- .../constitutive/solid/ElasticOrthotropic.hpp | 16 ++++-- .../solid/ElasticTransverseIsotropic.hpp | 16 ++++-- .../constitutive/solid/ModifiedCamClay.hpp | 17 ++++-- .../constitutive/solid/PerfectlyPlastic.hpp | 16 ++++-- .../constitutive/solid/SolidBase.cpp | 19 +++++++ .../constitutive/solid/SolidBase.hpp | 57 ++++++++++++++++++- .../constitutive/solid/SolidFields.hpp | 8 +++ .../ImplicitSmallStrainQuasiStatic.hpp | 8 ++- .../ImplicitSmallStrainQuasiStatic_impl.hpp | 4 +- 15 files changed, 248 insertions(+), 47 deletions(-) diff --git a/src/coreComponents/constitutive/solid/CeramicDamage.hpp b/src/coreComponents/constitutive/solid/CeramicDamage.hpp index ae1a9775be6..42a6d976379 100644 --- a/src/coreComponents/constitutive/solid/CeramicDamage.hpp +++ b/src/coreComponents/constitutive/solid/CeramicDamage.hpp @@ -67,8 +67,11 @@ class CeramicDamageUpdates : public ElasticIsotropicUpdates * @param[in] bulkModulus The ArrayView holding the bulk modulus data for each element. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. + * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. * @param[in] newStress The ArrayView holding the new stress data for each quadrature point. * @param[in] oldStress The ArrayView holding the old stress data for each quadrature point. + * @param[in] disableInelasticity Flag to disable plasticity for inelastic models + * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain */ CeramicDamageUpdates( arrayView2d< real64 > const & damage, arrayView2d< real64 > const & jacobian, @@ -80,10 +83,12 @@ class CeramicDamageUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, + arrayView1d< real64 const > const & anelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, - bool const & disableInelasticity ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, newStress, oldStress, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainMagnitude, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), m_damage( damage ), m_jacobian( jacobian ), m_lengthScale( lengthScale ), @@ -472,9 +477,11 @@ class CeramicDamage : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -498,9 +505,11 @@ class CeramicDamage : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } diff --git a/src/coreComponents/constitutive/solid/DelftEgg.hpp b/src/coreComponents/constitutive/solid/DelftEgg.hpp index 167301468f0..7f1688f4502 100644 --- a/src/coreComponents/constitutive/solid/DelftEgg.hpp +++ b/src/coreComponents/constitutive/solid/DelftEgg.hpp @@ -54,9 +54,11 @@ class DelftEggUpdates : public ElasticIsotropicUpdates * @param[in] bulkModulus The ArrayView holding the bulk modulus data for each element. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. + * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. * @param[in] newStress The ArrayView holding the new stress data for each quadrature point. * @param[in] oldStress The ArrayView holding the old stress data from the previous converged state for each point * @param[in] disableInelasticity Flag to disable plastic response + * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain */ DelftEggUpdates( arrayView1d< real64 const > const & recompressionIndex, arrayView1d< real64 const > const & virginCompressionIndex, @@ -67,10 +69,12 @@ class DelftEggUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, + arrayView1d< real64 const > const & anelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, - const bool & disableInelasticity ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, newStress, oldStress, disableInelasticity ), + const bool & disableInelasticity, + const integer & enableAnelasticStrain ): + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainMagnitude, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), m_recompressionIndex( recompressionIndex ), m_virginCompressionIndex( virginCompressionIndex ), m_cslSlope( cslSlope ), @@ -524,9 +528,11 @@ class DelftEgg : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -549,9 +555,11 @@ class DelftEgg : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } diff --git a/src/coreComponents/constitutive/solid/DruckerPrager.hpp b/src/coreComponents/constitutive/solid/DruckerPrager.hpp index e59494d7b06..f9627730cb5 100644 --- a/src/coreComponents/constitutive/solid/DruckerPrager.hpp +++ b/src/coreComponents/constitutive/solid/DruckerPrager.hpp @@ -53,8 +53,11 @@ class DruckerPragerUpdates : public ElasticIsotropicUpdates * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. + * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. * @param[in] newStress The ArrayView holding the new stress data for each quadrature point. * @param[in] oldStress The ArrayView holding the old stress data for each quadrature point. + * @param[in] disableInelasticity Flag to disable plasticity for inelastic models + * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain */ DruckerPragerUpdates( arrayView1d< real64 const > const & friction, arrayView1d< real64 const > const & dilation, @@ -64,10 +67,12 @@ class DruckerPragerUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, + arrayView1d< real64 const > const & anelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, - bool const & disableInelasticity ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, newStress, oldStress, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainMagnitude, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), m_friction( friction ), m_dilation( dilation ), m_hardening( hardening ), @@ -410,9 +415,11 @@ class DruckerPrager : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -434,9 +441,11 @@ class DruckerPrager : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } diff --git a/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp b/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp index dc79394677f..27fca639b00 100644 --- a/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp +++ b/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp @@ -46,7 +46,9 @@ class DruckerPragerExtendedUpdates : public ElasticIsotropicUpdates * @param[in] bulkModulus The ArrayView holding the bulk modulus data for each element. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. + * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. * @param[in] stress The ArrayView holding the stress data for each quadrature point. + * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain */ DruckerPragerExtendedUpdates( arrayView1d< real64 const > const & initialFriction, arrayView1d< real64 const > const & residualFriction, @@ -58,10 +60,12 @@ class DruckerPragerExtendedUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, + arrayView1d< real64 const > const & anelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, - bool const & disableInelasticity ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, newStress, oldStress, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainMagnitude, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), m_initialFriction( initialFriction ), m_residualFriction( residualFriction ), m_dilationRatio( dilationRatio ), @@ -441,9 +445,11 @@ class DruckerPragerExtended : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -467,9 +473,11 @@ class DruckerPragerExtended : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp index a6663232fe8..b2f327c7014 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp @@ -46,17 +46,21 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates * @param[in] bulkModulus The ArrayView holding the bulk modulus data for each element. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. + * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. * @param[in] newStress The ArrayView holding the new stress data for each quadrature point. * @param[in] oldStress The ArrayView holding the old stress data for each quadrature point. * @param[in] disableInelasticity Flag to disable plasticity for inelastic models + * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain */ ElasticIsotropicUpdates( arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, + arrayView1d< real64 const > const & anelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, - const bool & disableInelasticity ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity ), + const bool & disableInelasticity, + const integer & enableAnelasticStrain ): + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainMagnitude, enableAnelasticStrain ), m_bulkModulus( bulkModulus ), m_shearModulus( shearModulus ) {} @@ -157,6 +161,11 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates localIndex const q, real64 beta ) const override; + GEOS_HOST_DEVICE + virtual void stressModificationByAnelasticStain( localIndex const k, + localIndex const q, + real64 ( &stressModifier )[6] ) const override; + // TODO: confirm hyper stress/strain measures before activatiing /* @@ -374,6 +383,29 @@ void ElasticIsotropicUpdates::viscousStateUpdate( localIndex const k, GEOS_UNUSED_VAR( beta ); } +GEOS_HOST_DEVICE +GEOS_FORCE_INLINE +void ElasticIsotropicUpdates::stressModificationByAnelasticStain( localIndex const k, + localIndex const q, + real64 ( & stressModifier )[6] ) const +{ + if( m_enableAnelasticStrain == 0 ) + { + return; + } + + real64 const anelasticStrainDirection[6] = { 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0 }; // To make it an input + + real64 anelasticStrain[6]; + for( integer i = 0; i < 6; ++i ) + { + anelasticStrain[i] = getAnelasticStrainMagnitude( k ) * anelasticStrainDirection[i]; + } + + smallStrainNoStateUpdate_StressOnly( k, q, anelasticStrain, stressModifier ); +} + // TODO: need to confirm stress / strain measures before activating hyper inferface /* @@ -533,18 +565,22 @@ class ElasticIsotropic : public SolidBase return ElasticIsotropicUpdates( m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } else // for "no state" updates, pass empty views to avoid transfer of stress data to device { return ElasticIsotropicUpdates( m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD >(), arrayView3d< real64, solid::STRESS_USD >(), - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } } @@ -563,9 +599,11 @@ class ElasticIsotropic : public SolidBase m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } protected: diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp index 8db58f1ed5b..bc2698928d0 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp @@ -49,19 +49,23 @@ class ElasticIsotropicPressureDependentUpdates : public SolidBaseUpdates * @param[in] recompressionIndex The ArrayView holding the recompression index data for each element. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. + * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. * @param[in] newStress The ArrayView holding the new stress data for each quadrature point. * @param[in] oldStress The ArrayView holding the old stress data from the previous converged step for each quadrature point. * @param[in] disableInelasticity Flag to disable plastic response for inelastic models + * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain */ ElasticIsotropicPressureDependentUpdates( real64 const & refPressure, real64 const & refStrainVol, arrayView1d< real64 const > const & recompressionIndex, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, + arrayView1d< real64 const > const & anelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, - bool const & disableInelasticity ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainMagnitude, enableAnelasticStrain ), m_refPressure( refPressure ), m_refStrainVol( refStrainVol ), m_recompressionIndex( recompressionIndex ), @@ -583,9 +587,11 @@ class ElasticIsotropicPressureDependent : public SolidBase m_recompressionIndex, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } else // for "no state" updates, pass empty views to avoid transfer of stress data to device { @@ -594,9 +600,11 @@ class ElasticIsotropicPressureDependent : public SolidBase m_recompressionIndex, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD >(), arrayView3d< real64, solid::STRESS_USD >(), - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } } @@ -617,9 +625,11 @@ class ElasticIsotropicPressureDependent : public SolidBase m_recompressionIndex, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } diff --git a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp index e3f9075a73b..daa28c22435 100644 --- a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp @@ -52,9 +52,11 @@ class ElasticOrthotropicUpdates : public SolidBaseUpdates * @param[in] c55 The 55 component of the Voigt stiffness tensor. * @param[in] c66 The 66 component of the Voigt stiffness tensor. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. + * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. * @param[in] newStress The ArrayView holding the new stress data for each point. * @param[in] oldStress The ArrayView holding the old stress data for each point. * @param[in] disableInelasticity Flag to disable plastic response for inelastic models. + * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain */ ElasticOrthotropicUpdates( arrayView1d< real64 const > const & c11, arrayView1d< real64 const > const & c12, @@ -66,10 +68,12 @@ class ElasticOrthotropicUpdates : public SolidBaseUpdates arrayView1d< real64 const > const & c55, arrayView1d< real64 const > const & c66, arrayView1d< real64 const > const & thermalExpansionCoefficient, + arrayView1d< real64 const > const & anelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, - bool const & disableInelasticity ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainMagnitude, enableAnelasticStrain ), m_c11( c11 ), m_c12( c12 ), m_c13( c13 ), @@ -737,9 +741,11 @@ class ElasticOrthotropic : public SolidBase m_c55, m_c66, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -764,9 +770,11 @@ class ElasticOrthotropic : public SolidBase m_c55, m_c66, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } protected: diff --git a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp index b402a3fe816..9ea651a2ea2 100644 --- a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp @@ -50,9 +50,11 @@ class ElasticTransverseIsotropicUpdates : public SolidBaseUpdates * @param[in] c44 The 44 component of the Voigt stiffness tensor. * @param[in] c66 The 66 component of the Voigt stiffness tensor. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. + * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. * @param[in] newStress The ArrayView holding the new stress data for each point. * @param[in] oldStress The ArrayView holding the old stress data for each point. * @param[in] disableInelasticity Flag to disable plastic response for inelastic models. + * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain */ ElasticTransverseIsotropicUpdates( arrayView1d< real64 const > const & c11, arrayView1d< real64 const > const & c13, @@ -60,10 +62,12 @@ class ElasticTransverseIsotropicUpdates : public SolidBaseUpdates arrayView1d< real64 const > const & c44, arrayView1d< real64 const > const & c66, arrayView1d< real64 const > const & thermalExpansionCoefficient, + arrayView1d< real64 const > const & anelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, - bool const & disableInelasticity ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainMagnitude, enableAnelasticStrain ), m_c11( c11 ), m_c13( c13 ), m_c33( c33 ), @@ -593,9 +597,11 @@ class ElasticTransverseIsotropic : public SolidBase m_c44, m_c66, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -616,9 +622,11 @@ class ElasticTransverseIsotropic : public SolidBase m_c44, m_c66, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } protected: diff --git a/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp b/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp index 93b0d7da7e7..e65788616eb 100644 --- a/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp +++ b/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp @@ -54,10 +54,12 @@ class ModifiedCamClayUpdates : public ElasticIsotropicPressureDependentUpdates * for each quadrature point. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. + * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. * @param[in] newstress The ArrayView holding the new stress data for each quadrature point. * @param[in] oldstress The ArrayView holding the old stress data from the previous converged state for each quadrature * point. * @param[in] disableInelasticity Flag to disable plastic response/ + * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain */ ModifiedCamClayUpdates( real64 const & refPressure, real64 const & refStrainVol, @@ -68,10 +70,13 @@ class ModifiedCamClayUpdates : public ElasticIsotropicPressureDependentUpdates arrayView2d< real64 > const & oldPreConsolidationPressure, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, + arrayView1d< real64 const > const & anelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, - bool const & disableInelasticity ): - ElasticIsotropicPressureDependentUpdates( refPressure, refStrainVol, recompressionIndex, shearModulus, thermalExpansionCoefficient, newStress, oldStress, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + ElasticIsotropicPressureDependentUpdates( refPressure, refStrainVol, recompressionIndex, shearModulus, thermalExpansionCoefficient, anelasticStrainMagnitude, newStress, oldStress, + disableInelasticity, enableAnelasticStrain ), m_virginCompressionIndex( virginCompressionIndex ), m_cslSlope( cslSlope ), m_newPreConsolidationPressure( newPreConsolidationPressure ), @@ -535,9 +540,11 @@ class ModifiedCamClay : public ElasticIsotropicPressureDependent m_oldPreConsolidationPressure, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -560,9 +567,11 @@ class ModifiedCamClay : public ElasticIsotropicPressureDependent m_oldPreConsolidationPressure, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } protected: diff --git a/src/coreComponents/constitutive/solid/PerfectlyPlastic.hpp b/src/coreComponents/constitutive/solid/PerfectlyPlastic.hpp index 019ffcb0b9b..77412caf82e 100644 --- a/src/coreComponents/constitutive/solid/PerfectlyPlastic.hpp +++ b/src/coreComponents/constitutive/solid/PerfectlyPlastic.hpp @@ -48,17 +48,21 @@ class PerfectlyPlasticUpdates : public ElasticIsotropicUpdates * @param[in] bulkModulus The ArrayView holding the bulk modulus data for each element. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. + * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. * @param[in] newStress The ArrayView holding the new stress data for each quadrature point. * @param[in] oldStress The ArrayView holding the old stress data for each quadrature point. + * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain */ PerfectlyPlasticUpdates( arrayView1d< real64 const > const & yieldStress, arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, + arrayView1d< real64 const > const & anelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, - bool const & disableInelasticity ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, newStress, oldStress, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainMagnitude, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), m_yieldStress( yieldStress ) {} @@ -281,9 +285,11 @@ class PerfectlyPlastic : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -301,9 +307,11 @@ class PerfectlyPlastic : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainMagnitude, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } diff --git a/src/coreComponents/constitutive/solid/SolidBase.cpp b/src/coreComponents/constitutive/solid/SolidBase.cpp index 1e1a4dce868..75f3f1683de 100644 --- a/src/coreComponents/constitutive/solid/SolidBase.cpp +++ b/src/coreComponents/constitutive/solid/SolidBase.cpp @@ -40,6 +40,16 @@ SolidBase::SolidBase( string const & name, Group * const parent ): setInputFlag( InputFlags::OPTIONAL ). setDescription( "Default Linear Thermal Expansion Coefficient of the Solid Rock Frame" ); + registerWrapper( viewKeyStruct::defaultAnelasticStrainMagnitudeString(), &m_defaultAnelasticStrainMagnitude ). + setApplyDefaultValue( 0.0 ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Default anelastic strain magnitude" ); + + registerWrapper( viewKeyStruct::enableAnelasticStrainString(), &m_enableAnelasticStrain ). + setApplyDefaultValue( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Whether to enable stress modification due to anelastic strain. Can be 0 or 1" ); + // register fields string const voightLabels[6] = { "XX", "YY", "ZZ", "YZ", "XZ", "XY" }; @@ -53,6 +63,8 @@ SolidBase::SolidBase( string const & name, Group * const parent ): registerField< fields::solid::density >( &m_density ); registerField< fields::solid::thermalExpansionCoefficient >( &m_thermalExpansionCoefficient ); + + registerField< fields::solid::anelasticStrainMagnitude >( &m_anelasticStrainMagnitude ); } @@ -63,6 +75,13 @@ void SolidBase::postInputInitialization() getField< fields::solid::thermalExpansionCoefficient >(). setApplyDefaultValue( m_defaultThermalExpansionCoefficient ); + + getField< fields::solid::anelasticStrainMagnitude >(). + setApplyDefaultValue( m_defaultAnelasticStrainMagnitude ); + + GEOS_ERROR_IF( m_enableAnelasticStrain == 0 && m_defaultAnelasticStrainMagnitude > 0.0, + getDataContext() << ": enableAnelasticStrain flag must be 1 if a nonzero" + " AnelasticStrainMagnitude is used" ); } diff --git a/src/coreComponents/constitutive/solid/SolidBase.hpp b/src/coreComponents/constitutive/solid/SolidBase.hpp index f21092d7a61..45f3d7e2d91 100644 --- a/src/coreComponents/constitutive/solid/SolidBase.hpp +++ b/src/coreComponents/constitutive/solid/SolidBase.hpp @@ -55,15 +55,21 @@ class SolidBaseUpdates * @param[in] oldStress The old stress data from the constitutive model class. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. * @param[in] disableInelasticity Flag to disable inelastic response + * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. + * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain */ SolidBaseUpdates( arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, arrayView1d< real64 const > const & thermalExpansionCoefficient, - const bool & disableInelasticity ): + const bool & disableInelasticity, + arrayView1d< real64 const > const & anelasticStrainMagnitude, + const integer enableAnelasticStrain ): m_newStress( newStress ), m_oldStress( oldStress ), m_thermalExpansionCoefficient( thermalExpansionCoefficient ), - m_disableInelasticity ( disableInelasticity ) + m_disableInelasticity ( disableInelasticity ), + m_anelasticStrainMagnitude( anelasticStrainMagnitude ), + m_enableAnelasticStrain( enableAnelasticStrain ) {} /// Deleted default constructor @@ -101,6 +107,12 @@ class SolidBaseUpdates /// Flag to disable inelasticity const bool m_disableInelasticity; + /// The anelastic strain magnitude (i.e. chemistry, electrochemistry, etc.) + arrayView1d< real64 const > const m_anelasticStrainMagnitude; + + /// Flag to enable stress modification due to anelastic strain + const integer m_enableAnelasticStrain; + /** * @brief Get bulkModulus * @param[in] k Element index. @@ -126,6 +138,17 @@ class SolidBaseUpdates return m_thermalExpansionCoefficient[k]; } + /** + * @brief Get anelasticStrainMagnitude + * @param[in] k Element index. + * @return the nelasticStrainMagnitude of element k + */ + GEOS_HOST_DEVICE + real64 getAnelasticStrainMagnitude( localIndex const k ) const + { + return m_anelasticStrainMagnitude[k]; + } + /** * @brief Get shear modulus * @param[in] k Element index. @@ -371,6 +394,22 @@ class SolidBaseUpdates GEOS_ERROR( "viscousStateUpdate() not implemented for this model" ); } + /** + * @brief Calculate the stress modifier due to anelastic strain + * + * @param[out] stress New stress value (Cauchy stress) + */ + GEOS_HOST_DEVICE + virtual void stressModificationByAnelasticStain( localIndex const k, + localIndex const q, + real64 ( & stressModifier )[6] ) const + { + GEOS_UNUSED_VAR( k ); + GEOS_UNUSED_VAR( q ); + GEOS_UNUSED_VAR( stressModifier ); + GEOS_ERROR( "stressModificationByAnelasticStain() not implemented for this model" ); + } + /** * @brief Return the strain energy density at a given material point * @@ -571,6 +610,10 @@ class SolidBase : public constitutive::ConstitutiveBase static constexpr char const * defaultDensityString() { return "defaultDensity"; } // Default drained linear thermal expansion coefficient key static constexpr char const * defaultThermalExpansionCoefficientString() { return "defaultDrainedLinearTEC"; } + // Default anelastic strain magnitude key + static constexpr char const * defaultAnelasticStrainMagnitudeString() { return "defaultAnelasticStrainMagnitude"; } + // Enable stress modification due to anelastic strain key + static constexpr char const * enableAnelasticStrainString() { return "enableAnelasticStrain"; } }; /// Save state data in preparation for next timestep @@ -698,6 +741,16 @@ class SolidBase : public constitutive::ConstitutiveBase /// Flag to disable inelasticity (plasticity, damage, etc.) bool m_disableInelasticity = false; + + /// The anelastic strain magnitude (i.e. chemistry, electrochemistry, etc.) + array1d< real64 > m_anelasticStrainMagnitude; + + /// The default value of the anelastic strain magnitude + real64 m_defaultAnelasticStrainMagnitude; + + /// Flag to enable stress modification due to anelastic strain + integer m_enableAnelasticStrain; + }; } // namespace constitutive diff --git a/src/coreComponents/constitutive/solid/SolidFields.hpp b/src/coreComponents/constitutive/solid/SolidFields.hpp index a9d37e3db9b..86f28b8e138 100644 --- a/src/coreComponents/constitutive/solid/SolidFields.hpp +++ b/src/coreComponents/constitutive/solid/SolidFields.hpp @@ -420,6 +420,14 @@ DECLARE_FIELD( dInternalEnergy_dTemperature, WRITE_AND_READ, "Derivative of the solid internal energy w.r.t. temperature [J/(m^3.K)]" ); +DECLARE_FIELD( anelasticStrainMagnitude, + "anelasticStrainMagnitude", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Anelastic strain magnitude (i.e. chemistry, electrochemistry, etc.)" ); + } } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic.hpp index 0e9bcf87b18..c998fe55e0f 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic.hpp @@ -184,11 +184,15 @@ class ImplicitSmallStrainQuasiStatic : * @brief operator() no-op used for modifying the stress tensor prior to * integrating the divergence to produce nodal forces. * @param stress The stress array. + * @param stressModification The stressModification array. */ GEOS_HOST_DEVICE inline constexpr - void operator() ( real64 (& stress)[6] ) + void operator() ( real64 (& stress)[6], real64 const (&stressModification)[6] ) { - GEOS_UNUSED_VAR( stress ); + for( localIndex i = 0; i < 6; ++i ) + { + stress[i] -= stressModification[i]; + } } }; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic_impl.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic_impl.hpp index 12753ab393f..b29a1d9d47e 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic_impl.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic_impl.hpp @@ -124,14 +124,16 @@ void ImplicitSmallStrainQuasiStatic< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE real64 strainInc[6] = {0}; real64 stress[6] = {0}; + real64 stressModifierAnelasticStrain[6] = {0}; typename CONSTITUTIVE_TYPE::KernelWrapper::DiscretizationOps stiffness; FE_TYPE::symmetricGradient( dNdX, stack.uhat_local, strainInc ); m_constitutiveUpdate.smallStrainUpdate( k, q, m_dt, strainInc, stress, stiffness ); + m_constitutiveUpdate.stressModificationByAnelasticStain( k, q, stressModifierAnelasticStrain ); - stressModifier( stress ); + stressModifier( stress, stressModifierAnelasticStrain ); // #pragma unroll for( localIndex i=0; i<6; ++i ) { From 4d36d9a815fe9c640e9de96864442130b653aa79 Mon Sep 17 00:00:00 2001 From: Fan Fei Date: Tue, 7 Oct 2025 08:37:14 -0700 Subject: [PATCH 2/4] make anelastic strain increasing with time --- .../constitutive/solid/CeramicDamage.hpp | 15 +++++--- .../constitutive/solid/DelftEgg.hpp | 15 +++++--- .../constitutive/solid/DruckerPrager.hpp | 15 +++++--- .../solid/DruckerPragerExtended.hpp | 15 +++++--- .../constitutive/solid/ElasticIsotropic.hpp | 23 +++++++++---- .../ElasticIsotropicPressureDependent.hpp | 19 ++++++++--- .../constitutive/solid/ElasticOrthotropic.hpp | 15 +++++--- .../solid/ElasticTransverseIsotropic.hpp | 15 +++++--- .../constitutive/solid/ModifiedCamClay.hpp | 15 +++++--- .../constitutive/solid/PerfectlyPlastic.hpp | 15 +++++--- .../constitutive/solid/SolidBase.cpp | 17 +++++++--- .../constitutive/solid/SolidBase.hpp | 34 ++++++++++++------- .../constitutive/solid/SolidFields.hpp | 22 ++++++++++-- 13 files changed, 172 insertions(+), 63 deletions(-) diff --git a/src/coreComponents/constitutive/solid/CeramicDamage.hpp b/src/coreComponents/constitutive/solid/CeramicDamage.hpp index 42a6d976379..14765786d62 100644 --- a/src/coreComponents/constitutive/solid/CeramicDamage.hpp +++ b/src/coreComponents/constitutive/solid/CeramicDamage.hpp @@ -83,12 +83,15 @@ class CeramicDamageUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainMagnitude, + arrayView1d< real64 const > const & anelasticStrainIncrement, + arrayView1d< real64 > const & newAnelasticStrainMagnitude, + arrayView1d< real64 > const & oldAnelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainMagnitude, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, newStress, oldStress, + disableInelasticity, enableAnelasticStrain ), m_damage( damage ), m_jacobian( jacobian ), m_lengthScale( lengthScale ), @@ -477,7 +480,9 @@ class CeramicDamage : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, @@ -505,7 +510,9 @@ class CeramicDamage : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/DelftEgg.hpp b/src/coreComponents/constitutive/solid/DelftEgg.hpp index 7f1688f4502..03a37055b5b 100644 --- a/src/coreComponents/constitutive/solid/DelftEgg.hpp +++ b/src/coreComponents/constitutive/solid/DelftEgg.hpp @@ -69,12 +69,15 @@ class DelftEggUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainMagnitude, + arrayView1d< real64 const > const & anelasticStrainIncrement, + arrayView1d< real64 > const & newAnelasticStrainMagnitude, + arrayView1d< real64 > const & oldAnelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, const bool & disableInelasticity, const integer & enableAnelasticStrain ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainMagnitude, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, newStress, oldStress, + disableInelasticity, enableAnelasticStrain ), m_recompressionIndex( recompressionIndex ), m_virginCompressionIndex( virginCompressionIndex ), m_cslSlope( cslSlope ), @@ -528,7 +531,9 @@ class DelftEgg : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, @@ -555,7 +560,9 @@ class DelftEgg : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/DruckerPrager.hpp b/src/coreComponents/constitutive/solid/DruckerPrager.hpp index f9627730cb5..645c7cd1255 100644 --- a/src/coreComponents/constitutive/solid/DruckerPrager.hpp +++ b/src/coreComponents/constitutive/solid/DruckerPrager.hpp @@ -67,12 +67,15 @@ class DruckerPragerUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainMagnitude, + arrayView1d< real64 const > const & anelasticStrainIncrement, + arrayView1d< real64 > const & newAnelasticStrainMagnitude, + arrayView1d< real64 > const & oldAnelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainMagnitude, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, newStress, oldStress, + disableInelasticity, enableAnelasticStrain ), m_friction( friction ), m_dilation( dilation ), m_hardening( hardening ), @@ -415,7 +418,9 @@ class DruckerPrager : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, @@ -441,7 +446,9 @@ class DruckerPrager : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp b/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp index 27fca639b00..6829dc18f42 100644 --- a/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp +++ b/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp @@ -60,12 +60,15 @@ class DruckerPragerExtendedUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainMagnitude, + arrayView1d< real64 const > const & anelasticStrainIncrement, + arrayView1d< real64 > const & newAnelasticStrainMagnitude, + arrayView1d< real64 > const & oldAnelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainMagnitude, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, newStress, oldStress, + disableInelasticity, enableAnelasticStrain ), m_initialFriction( initialFriction ), m_residualFriction( residualFriction ), m_dilationRatio( dilationRatio ), @@ -445,7 +448,9 @@ class DruckerPragerExtended : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, @@ -473,7 +478,9 @@ class DruckerPragerExtended : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp index b2f327c7014..c0c68aea4aa 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp @@ -55,12 +55,15 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates ElasticIsotropicUpdates( arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainMagnitude, + arrayView1d< real64 const > const & anelasticStrainIncrement, + arrayView1d< real64 > const & newAnelasticStrainMagnitude, + arrayView1d< real64 > const & oldAnelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, const bool & disableInelasticity, const integer & enableAnelasticStrain ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainMagnitude, enableAnelasticStrain ), + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, + enableAnelasticStrain ), m_bulkModulus( bulkModulus ), m_shearModulus( shearModulus ) {} @@ -397,10 +400,12 @@ void ElasticIsotropicUpdates::stressModificationByAnelasticStain( localIndex con real64 const anelasticStrainDirection[6] = { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0 }; // To make it an input + m_newAnelasticStrainMagnitude[k] = m_oldAnelasticStrainMagnitude[k] + getAnelasticStrainIncrement( k ); + real64 anelasticStrain[6]; for( integer i = 0; i < 6; ++i ) { - anelasticStrain[i] = getAnelasticStrainMagnitude( k ) * anelasticStrainDirection[i]; + anelasticStrain[i] = m_newAnelasticStrainMagnitude[k] * anelasticStrainDirection[i]; } smallStrainNoStateUpdate_StressOnly( k, q, anelasticStrain, stressModifier ); @@ -565,7 +570,9 @@ class ElasticIsotropic : public SolidBase return ElasticIsotropicUpdates( m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, @@ -576,7 +583,9 @@ class ElasticIsotropic : public SolidBase return ElasticIsotropicUpdates( m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD >(), arrayView3d< real64, solid::STRESS_USD >(), m_disableInelasticity, @@ -599,7 +608,9 @@ class ElasticIsotropic : public SolidBase m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp index bc2698928d0..0f228494f44 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp @@ -60,12 +60,15 @@ class ElasticIsotropicPressureDependentUpdates : public SolidBaseUpdates arrayView1d< real64 const > const & recompressionIndex, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainMagnitude, + arrayView1d< real64 const > const & anelasticStrainIncrement, + arrayView1d< real64 > const & newAnelasticStrainMagnitude, + arrayView1d< real64 > const & oldAnelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainMagnitude, enableAnelasticStrain ), + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, + enableAnelasticStrain ), m_refPressure( refPressure ), m_refStrainVol( refStrainVol ), m_recompressionIndex( recompressionIndex ), @@ -587,7 +590,9 @@ class ElasticIsotropicPressureDependent : public SolidBase m_recompressionIndex, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, @@ -600,7 +605,9 @@ class ElasticIsotropicPressureDependent : public SolidBase m_recompressionIndex, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD >(), arrayView3d< real64, solid::STRESS_USD >(), m_disableInelasticity, @@ -625,7 +632,9 @@ class ElasticIsotropicPressureDependent : public SolidBase m_recompressionIndex, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp index daa28c22435..9a613a80193 100644 --- a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp @@ -68,12 +68,15 @@ class ElasticOrthotropicUpdates : public SolidBaseUpdates arrayView1d< real64 const > const & c55, arrayView1d< real64 const > const & c66, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainMagnitude, + arrayView1d< real64 const > const & anelasticStrainIncrement, + arrayView1d< real64 > const & newAnelasticStrainMagnitude, + arrayView1d< real64 > const & oldAnelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainMagnitude, enableAnelasticStrain ), + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, + enableAnelasticStrain ), m_c11( c11 ), m_c12( c12 ), m_c13( c13 ), @@ -741,7 +744,9 @@ class ElasticOrthotropic : public SolidBase m_c55, m_c66, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, @@ -770,7 +775,9 @@ class ElasticOrthotropic : public SolidBase m_c55, m_c66, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp index 9ea651a2ea2..7ec996679d5 100644 --- a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp @@ -62,12 +62,15 @@ class ElasticTransverseIsotropicUpdates : public SolidBaseUpdates arrayView1d< real64 const > const & c44, arrayView1d< real64 const > const & c66, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainMagnitude, + arrayView1d< real64 const > const & anelasticStrainIncrement, + arrayView1d< real64 > const & newAnelasticStrainMagnitude, + arrayView1d< real64 > const & oldAnelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainMagnitude, enableAnelasticStrain ), + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, + enableAnelasticStrain ), m_c11( c11 ), m_c13( c13 ), m_c33( c33 ), @@ -597,7 +600,9 @@ class ElasticTransverseIsotropic : public SolidBase m_c44, m_c66, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, @@ -622,7 +627,9 @@ class ElasticTransverseIsotropic : public SolidBase m_c44, m_c66, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp b/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp index e65788616eb..f771d81d836 100644 --- a/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp +++ b/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp @@ -70,12 +70,15 @@ class ModifiedCamClayUpdates : public ElasticIsotropicPressureDependentUpdates arrayView2d< real64 > const & oldPreConsolidationPressure, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainMagnitude, + arrayView1d< real64 const > const & anelasticStrainIncrement, + arrayView1d< real64 > const & newAnelasticStrainMagnitude, + arrayView1d< real64 > const & oldAnelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - ElasticIsotropicPressureDependentUpdates( refPressure, refStrainVol, recompressionIndex, shearModulus, thermalExpansionCoefficient, anelasticStrainMagnitude, newStress, oldStress, + ElasticIsotropicPressureDependentUpdates( refPressure, refStrainVol, recompressionIndex, shearModulus, thermalExpansionCoefficient, anelasticStrainIncrement, newAnelasticStrainMagnitude, + oldAnelasticStrainMagnitude, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), m_virginCompressionIndex( virginCompressionIndex ), m_cslSlope( cslSlope ), @@ -540,7 +543,9 @@ class ModifiedCamClay : public ElasticIsotropicPressureDependent m_oldPreConsolidationPressure, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, @@ -567,7 +572,9 @@ class ModifiedCamClay : public ElasticIsotropicPressureDependent m_oldPreConsolidationPressure, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/PerfectlyPlastic.hpp b/src/coreComponents/constitutive/solid/PerfectlyPlastic.hpp index 77412caf82e..22afbc73b69 100644 --- a/src/coreComponents/constitutive/solid/PerfectlyPlastic.hpp +++ b/src/coreComponents/constitutive/solid/PerfectlyPlastic.hpp @@ -57,12 +57,15 @@ class PerfectlyPlasticUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainMagnitude, + arrayView1d< real64 const > const & anelasticStrainIncrement, + arrayView1d< real64 > const & newAnelasticStrainMagnitude, + arrayView1d< real64 > const & oldAnelasticStrainMagnitude, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainMagnitude, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, newStress, oldStress, + disableInelasticity, enableAnelasticStrain ), m_yieldStress( yieldStress ) {} @@ -285,7 +288,9 @@ class PerfectlyPlastic : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, @@ -307,7 +312,9 @@ class PerfectlyPlastic : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainMagnitude, + m_anelasticStrainIncrement, + m_newAnelasticStrainMagnitude, + m_oldAnelasticStrainMagnitude, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/SolidBase.cpp b/src/coreComponents/constitutive/solid/SolidBase.cpp index 75f3f1683de..6192742fe3a 100644 --- a/src/coreComponents/constitutive/solid/SolidBase.cpp +++ b/src/coreComponents/constitutive/solid/SolidBase.cpp @@ -40,7 +40,7 @@ SolidBase::SolidBase( string const & name, Group * const parent ): setInputFlag( InputFlags::OPTIONAL ). setDescription( "Default Linear Thermal Expansion Coefficient of the Solid Rock Frame" ); - registerWrapper( viewKeyStruct::defaultAnelasticStrainMagnitudeString(), &m_defaultAnelasticStrainMagnitude ). + registerWrapper( viewKeyStruct::defaultAnelasticStrainIncrementString(), &m_defaultAnelasticStrainIncrement ). setApplyDefaultValue( 0.0 ). setInputFlag( InputFlags::OPTIONAL ). setDescription( "Default anelastic strain magnitude" ); @@ -64,7 +64,9 @@ SolidBase::SolidBase( string const & name, Group * const parent ): registerField< fields::solid::thermalExpansionCoefficient >( &m_thermalExpansionCoefficient ); - registerField< fields::solid::anelasticStrainMagnitude >( &m_anelasticStrainMagnitude ); + registerField< fields::solid::anelasticStrainIncrement >( &m_anelasticStrainIncrement ); + registerField< fields::solid::newAnelasticStrainMagnitude >( &m_newAnelasticStrainMagnitude ); + registerField< fields::solid::oldAnelasticStrainMagnitude >( &m_oldAnelasticStrainMagnitude ); } @@ -76,10 +78,10 @@ void SolidBase::postInputInitialization() getField< fields::solid::thermalExpansionCoefficient >(). setApplyDefaultValue( m_defaultThermalExpansionCoefficient ); - getField< fields::solid::anelasticStrainMagnitude >(). - setApplyDefaultValue( m_defaultAnelasticStrainMagnitude ); + getField< fields::solid::anelasticStrainIncrement >(). + setApplyDefaultValue( m_defaultAnelasticStrainIncrement ); - GEOS_ERROR_IF( m_enableAnelasticStrain == 0 && m_defaultAnelasticStrainMagnitude > 0.0, + GEOS_ERROR_IF( m_enableAnelasticStrain == 0 && m_defaultAnelasticStrainIncrement > 0.0, getDataContext() << ": enableAnelasticStrain flag must be 1 if a nonzero" " AnelasticStrainMagnitude is used" ); } @@ -104,8 +106,13 @@ void SolidBase::saveConvergedState() const arrayView3d< real64 const, solid::STRESS_USD > newStress = m_newStress; arrayView3d< real64, solid::STRESS_USD > oldStress = m_oldStress; + arrayView1d< real64 const > newAnelasticStrainMagnitude = m_newAnelasticStrainMagnitude; + arrayView1d< real64 > oldAnelasticStrainMagnitude = m_oldAnelasticStrainMagnitude; + forAll< parallelDevicePolicy<> >( numE, [=] GEOS_HOST_DEVICE ( localIndex const k ) { + oldAnelasticStrainMagnitude[k] = newAnelasticStrainMagnitude[k]; + for( localIndex q = 0; q < numQ; ++q ) { LvArray::tensorOps::copy< 6 >( oldStress[k][q], newStress[k][q] ); diff --git a/src/coreComponents/constitutive/solid/SolidBase.hpp b/src/coreComponents/constitutive/solid/SolidBase.hpp index 45f3d7e2d91..9ba9bf3e4a5 100644 --- a/src/coreComponents/constitutive/solid/SolidBase.hpp +++ b/src/coreComponents/constitutive/solid/SolidBase.hpp @@ -62,13 +62,17 @@ class SolidBaseUpdates arrayView3d< real64, solid::STRESS_USD > const & oldStress, arrayView1d< real64 const > const & thermalExpansionCoefficient, const bool & disableInelasticity, - arrayView1d< real64 const > const & anelasticStrainMagnitude, + arrayView1d< real64 const > const & anelasticStrainIncrement, + arrayView1d< real64 > const & newAnelasticStrainMagnitude, + arrayView1d< real64 > const & oldAnelasticStrainMagnitude, const integer enableAnelasticStrain ): m_newStress( newStress ), m_oldStress( oldStress ), m_thermalExpansionCoefficient( thermalExpansionCoefficient ), m_disableInelasticity ( disableInelasticity ), - m_anelasticStrainMagnitude( anelasticStrainMagnitude ), + m_anelasticStrainIncrement( anelasticStrainIncrement ), + m_newAnelasticStrainMagnitude( newAnelasticStrainMagnitude ), + m_oldAnelasticStrainMagnitude( oldAnelasticStrainMagnitude ), m_enableAnelasticStrain( enableAnelasticStrain ) {} @@ -108,7 +112,10 @@ class SolidBaseUpdates const bool m_disableInelasticity; /// The anelastic strain magnitude (i.e. chemistry, electrochemistry, etc.) - arrayView1d< real64 const > const m_anelasticStrainMagnitude; + arrayView1d< real64 const > const m_anelasticStrainIncrement; + + arrayView1d< real64 > const m_newAnelasticStrainMagnitude; + arrayView1d< real64 > const m_oldAnelasticStrainMagnitude; /// Flag to enable stress modification due to anelastic strain const integer m_enableAnelasticStrain; @@ -139,14 +146,14 @@ class SolidBaseUpdates } /** - * @brief Get anelasticStrainMagnitude + * @brief Get anelasticStrainIncrement * @param[in] k Element index. - * @return the nelasticStrainMagnitude of element k + * @return the anelasticStrainIncrement of element k */ GEOS_HOST_DEVICE - real64 getAnelasticStrainMagnitude( localIndex const k ) const + real64 getAnelasticStrainIncrement( localIndex const k ) const { - return m_anelasticStrainMagnitude[k]; + return m_anelasticStrainIncrement[k]; } /** @@ -611,7 +618,7 @@ class SolidBase : public constitutive::ConstitutiveBase // Default drained linear thermal expansion coefficient key static constexpr char const * defaultThermalExpansionCoefficientString() { return "defaultDrainedLinearTEC"; } // Default anelastic strain magnitude key - static constexpr char const * defaultAnelasticStrainMagnitudeString() { return "defaultAnelasticStrainMagnitude"; } + static constexpr char const * defaultAnelasticStrainIncrementString() { return "defaultAnelasticStrainIncrement"; } // Enable stress modification due to anelastic strain key static constexpr char const * enableAnelasticStrainString() { return "enableAnelasticStrain"; } }; @@ -742,11 +749,14 @@ class SolidBase : public constitutive::ConstitutiveBase /// Flag to disable inelasticity (plasticity, damage, etc.) bool m_disableInelasticity = false; - /// The anelastic strain magnitude (i.e. chemistry, electrochemistry, etc.) - array1d< real64 > m_anelasticStrainMagnitude; + /// The anelastic strain rate magnitude (i.e. chemistry, electrochemistry, etc.) + array1d< real64 > m_anelasticStrainIncrement; + + array1d< real64 > m_newAnelasticStrainMagnitude; + array1d< real64 > m_oldAnelasticStrainMagnitude; - /// The default value of the anelastic strain magnitude - real64 m_defaultAnelasticStrainMagnitude; + /// The default value of the anelastic strain rate magnitude + real64 m_defaultAnelasticStrainIncrement; /// Flag to enable stress modification due to anelastic strain integer m_enableAnelasticStrain; diff --git a/src/coreComponents/constitutive/solid/SolidFields.hpp b/src/coreComponents/constitutive/solid/SolidFields.hpp index 86f28b8e138..8064e969e1b 100644 --- a/src/coreComponents/constitutive/solid/SolidFields.hpp +++ b/src/coreComponents/constitutive/solid/SolidFields.hpp @@ -420,13 +420,29 @@ DECLARE_FIELD( dInternalEnergy_dTemperature, WRITE_AND_READ, "Derivative of the solid internal energy w.r.t. temperature [J/(m^3.K)]" ); -DECLARE_FIELD( anelasticStrainMagnitude, - "anelasticStrainMagnitude", +DECLARE_FIELD( anelasticStrainIncrement, + "anelasticStrainIncrement", array1d< real64 >, 0, LEVEL_0, WRITE_AND_READ, - "Anelastic strain magnitude (i.e. chemistry, electrochemistry, etc.)" ); + "Anelastic strain increment (i.e. chemistry, electrochemistry, etc.)" ); + +DECLARE_FIELD( newAnelasticStrainMagnitude, + "newAnelasticStrainMagnitude", + array1d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "New anelastic strain magnitude" ); + +DECLARE_FIELD( oldAnelasticStrainMagnitude, + "oldAnelasticStrainMagnitude", + array1d< real64 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Old anelastic strain magnitude" ); } From dcfae80e2da539c69090592f7ca1eedc01257ab4 Mon Sep 17 00:00:00 2001 From: Fan Fei Date: Mon, 3 Aug 2026 20:08:50 -0700 Subject: [PATCH 3/4] changed the input anelastic strain to a 3D vector, changed increment to rate, renamed magnitude --- .../constitutive/solid/CeramicDamage.hpp | 22 ++++---- .../constitutive/solid/DelftEgg.hpp | 22 ++++---- .../constitutive/solid/DruckerPrager.hpp | 22 ++++---- .../solid/DruckerPragerExtended.hpp | 22 ++++---- .../constitutive/solid/ElasticIsotropic.hpp | 43 +++++++-------- .../ElasticIsotropicPressureDependent.hpp | 28 +++++----- .../constitutive/solid/ElasticOrthotropic.hpp | 22 ++++---- .../solid/ElasticTransverseIsotropic.hpp | 22 ++++---- .../constitutive/solid/ModifiedCamClay.hpp | 24 ++++---- .../constitutive/solid/PerfectlyPlastic.hpp | 22 ++++---- .../constitutive/solid/SolidBase.cpp | 42 +++++++++----- .../constitutive/solid/SolidBase.hpp | 55 +++++++++++-------- .../constitutive/solid/SolidFields.hpp | 24 ++++---- .../ImplicitSmallStrainQuasiStatic_impl.hpp | 2 +- 14 files changed, 198 insertions(+), 174 deletions(-) diff --git a/src/coreComponents/constitutive/solid/CeramicDamage.hpp b/src/coreComponents/constitutive/solid/CeramicDamage.hpp index 14765786d62..f2fd9012e60 100644 --- a/src/coreComponents/constitutive/solid/CeramicDamage.hpp +++ b/src/coreComponents/constitutive/solid/CeramicDamage.hpp @@ -67,7 +67,7 @@ class CeramicDamageUpdates : public ElasticIsotropicUpdates * @param[in] bulkModulus The ArrayView holding the bulk modulus data for each element. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. - * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. + * @param[in] anelasticStrainRate The ArrayView holding the anelastic strain rate data for each element. * @param[in] newStress The ArrayView holding the new stress data for each quadrature point. * @param[in] oldStress The ArrayView holding the old stress data for each quadrature point. * @param[in] disableInelasticity Flag to disable plasticity for inelastic models @@ -83,14 +83,14 @@ class CeramicDamageUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainIncrement, - arrayView1d< real64 > const & newAnelasticStrainMagnitude, - arrayView1d< real64 > const & oldAnelasticStrainMagnitude, + arrayView2d< real64 const > const & anelasticStrainRate, + arrayView2d< real64 > const & newAnelasticStrain, + arrayView2d< real64 > const & oldAnelasticStrain, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, newStress, oldStress, + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), m_damage( damage ), m_jacobian( jacobian ), @@ -480,9 +480,9 @@ class CeramicDamage : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, @@ -510,9 +510,9 @@ class CeramicDamage : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/DelftEgg.hpp b/src/coreComponents/constitutive/solid/DelftEgg.hpp index 03a37055b5b..25b51d10da1 100644 --- a/src/coreComponents/constitutive/solid/DelftEgg.hpp +++ b/src/coreComponents/constitutive/solid/DelftEgg.hpp @@ -54,7 +54,7 @@ class DelftEggUpdates : public ElasticIsotropicUpdates * @param[in] bulkModulus The ArrayView holding the bulk modulus data for each element. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. - * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. + * @param[in] anelasticStrainRate The ArrayView holding the anelastic strain rate data for each element. * @param[in] newStress The ArrayView holding the new stress data for each quadrature point. * @param[in] oldStress The ArrayView holding the old stress data from the previous converged state for each point * @param[in] disableInelasticity Flag to disable plastic response @@ -69,14 +69,14 @@ class DelftEggUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainIncrement, - arrayView1d< real64 > const & newAnelasticStrainMagnitude, - arrayView1d< real64 > const & oldAnelasticStrainMagnitude, + arrayView2d< real64 const > const & anelasticStrainRate, + arrayView2d< real64 > const & newAnelasticStrain, + arrayView2d< real64 > const & oldAnelasticStrain, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, const bool & disableInelasticity, const integer & enableAnelasticStrain ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, newStress, oldStress, + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), m_recompressionIndex( recompressionIndex ), m_virginCompressionIndex( virginCompressionIndex ), @@ -531,9 +531,9 @@ class DelftEgg : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, @@ -560,9 +560,9 @@ class DelftEgg : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/DruckerPrager.hpp b/src/coreComponents/constitutive/solid/DruckerPrager.hpp index 645c7cd1255..caca60db597 100644 --- a/src/coreComponents/constitutive/solid/DruckerPrager.hpp +++ b/src/coreComponents/constitutive/solid/DruckerPrager.hpp @@ -53,7 +53,7 @@ class DruckerPragerUpdates : public ElasticIsotropicUpdates * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. - * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. + * @param[in] anelasticStrainRate The ArrayView holding the anelastic strain rate data for each element. * @param[in] newStress The ArrayView holding the new stress data for each quadrature point. * @param[in] oldStress The ArrayView holding the old stress data for each quadrature point. * @param[in] disableInelasticity Flag to disable plasticity for inelastic models @@ -67,14 +67,14 @@ class DruckerPragerUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainIncrement, - arrayView1d< real64 > const & newAnelasticStrainMagnitude, - arrayView1d< real64 > const & oldAnelasticStrainMagnitude, + arrayView2d< real64 const > const & anelasticStrainRate, + arrayView2d< real64 > const & newAnelasticStrain, + arrayView2d< real64 > const & oldAnelasticStrain, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, newStress, oldStress, + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), m_friction( friction ), m_dilation( dilation ), @@ -418,9 +418,9 @@ class DruckerPrager : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, @@ -446,9 +446,9 @@ class DruckerPrager : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp b/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp index 6829dc18f42..5aa94492fc6 100644 --- a/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp +++ b/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp @@ -46,7 +46,7 @@ class DruckerPragerExtendedUpdates : public ElasticIsotropicUpdates * @param[in] bulkModulus The ArrayView holding the bulk modulus data for each element. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. - * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. + * @param[in] anelasticStrainRate The ArrayView holding the anelastic strain rate data for each element. * @param[in] stress The ArrayView holding the stress data for each quadrature point. * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain */ @@ -60,14 +60,14 @@ class DruckerPragerExtendedUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainIncrement, - arrayView1d< real64 > const & newAnelasticStrainMagnitude, - arrayView1d< real64 > const & oldAnelasticStrainMagnitude, + arrayView2d< real64 const > const & anelasticStrainRate, + arrayView2d< real64 > const & newAnelasticStrain, + arrayView2d< real64 > const & oldAnelasticStrain, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, newStress, oldStress, + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), m_initialFriction( initialFriction ), m_residualFriction( residualFriction ), @@ -448,9 +448,9 @@ class DruckerPragerExtended : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, @@ -478,9 +478,9 @@ class DruckerPragerExtended : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp index ed38e788c81..58c966354d0 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp @@ -46,7 +46,7 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates * @param[in] bulkModulus The ArrayView holding the bulk modulus data for each element. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. - * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. + * @param[in] anelasticStrainRate The ArrayView holding the anelastic strain rate data for each element. * @param[in] newStress The ArrayView holding the new stress data for each quadrature point. * @param[in] oldStress The ArrayView holding the old stress data for each quadrature point. * @param[in] disableInelasticity Flag to disable plasticity for inelastic models @@ -55,14 +55,14 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates ElasticIsotropicUpdates( arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainIncrement, - arrayView1d< real64 > const & newAnelasticStrainMagnitude, - arrayView1d< real64 > const & oldAnelasticStrainMagnitude, + arrayView2d< real64 const > const & anelasticStrainRate, + arrayView2d< real64 > const & newAnelasticStrain, + arrayView2d< real64 > const & oldAnelasticStrain, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, const bool & disableInelasticity, const integer & enableAnelasticStrain ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, enableAnelasticStrain ), m_bulkModulus( bulkModulus ), m_shearModulus( shearModulus ) @@ -167,6 +167,7 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates GEOS_HOST_DEVICE virtual void stressModificationByAnelasticStain( localIndex const k, localIndex const q, + real64 const & timeIncrement, real64 ( &stressModifier )[6] ) const override; // TODO: confirm hyper stress/strain measures before activatiing @@ -390,6 +391,7 @@ GEOS_HOST_DEVICE GEOS_FORCE_INLINE void ElasticIsotropicUpdates::stressModificationByAnelasticStain( localIndex const k, localIndex const q, + real64 const & timeIncrement, real64 ( & stressModifier )[6] ) const { if( m_enableAnelasticStrain == 0 ) @@ -397,15 +399,12 @@ void ElasticIsotropicUpdates::stressModificationByAnelasticStain( localIndex con return; } - real64 const anelasticStrainDirection[6] = { 0.0, 1.0, 0.0, - 0.0, 0.0, 0.0 }; // To make it an input - - m_newAnelasticStrainMagnitude[k] = m_oldAnelasticStrainMagnitude[k] + getAnelasticStrainIncrement( k ); - - real64 anelasticStrain[6]; - for( integer i = 0; i < 6; ++i ) + arraySlice1d< real64 const > const anelasticStrainRate = getAnelasticStrainRate( k ); + real64 anelasticStrain[6] = {}; + for( integer i = 0; i < 3; ++i ) { - anelasticStrain[i] = m_newAnelasticStrainMagnitude[k] * anelasticStrainDirection[i]; + m_newAnelasticStrain[k][i] = m_oldAnelasticStrain[k][i] + anelasticStrainRate[i] * timeIncrement; + anelasticStrain[i] = m_newAnelasticStrain[k][i]; } smallStrainNoStateUpdate_StressOnly( k, q, anelasticStrain, stressModifier ); @@ -570,9 +569,9 @@ class ElasticIsotropic : public SolidBase return ElasticIsotropicUpdates( m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, @@ -583,9 +582,9 @@ class ElasticIsotropic : public SolidBase return ElasticIsotropicUpdates( m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, arrayView3d< real64, solid::STRESS_USD >(), arrayView3d< real64, solid::STRESS_USD >(), m_disableInelasticity, @@ -608,9 +607,9 @@ class ElasticIsotropic : public SolidBase m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp index 0f228494f44..3176f714ee6 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp @@ -49,7 +49,7 @@ class ElasticIsotropicPressureDependentUpdates : public SolidBaseUpdates * @param[in] recompressionIndex The ArrayView holding the recompression index data for each element. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. - * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. + * @param[in] anelasticStrainRate The ArrayView holding the anelastic strain rate data for each element. * @param[in] newStress The ArrayView holding the new stress data for each quadrature point. * @param[in] oldStress The ArrayView holding the old stress data from the previous converged step for each quadrature point. * @param[in] disableInelasticity Flag to disable plastic response for inelastic models @@ -60,14 +60,14 @@ class ElasticIsotropicPressureDependentUpdates : public SolidBaseUpdates arrayView1d< real64 const > const & recompressionIndex, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainIncrement, - arrayView1d< real64 > const & newAnelasticStrainMagnitude, - arrayView1d< real64 > const & oldAnelasticStrainMagnitude, + arrayView2d< real64 const > const & anelasticStrainRate, + arrayView2d< real64 > const & newAnelasticStrain, + arrayView2d< real64 > const & oldAnelasticStrain, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, enableAnelasticStrain ), m_refPressure( refPressure ), m_refStrainVol( refStrainVol ), @@ -590,9 +590,9 @@ class ElasticIsotropicPressureDependent : public SolidBase m_recompressionIndex, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, @@ -605,9 +605,9 @@ class ElasticIsotropicPressureDependent : public SolidBase m_recompressionIndex, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, arrayView3d< real64, solid::STRESS_USD >(), arrayView3d< real64, solid::STRESS_USD >(), m_disableInelasticity, @@ -632,9 +632,9 @@ class ElasticIsotropicPressureDependent : public SolidBase m_recompressionIndex, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp index 9a613a80193..0498f18b108 100644 --- a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp @@ -52,7 +52,7 @@ class ElasticOrthotropicUpdates : public SolidBaseUpdates * @param[in] c55 The 55 component of the Voigt stiffness tensor. * @param[in] c66 The 66 component of the Voigt stiffness tensor. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. - * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. + * @param[in] anelasticStrainRate The ArrayView holding the anelastic strain rate data for each element. * @param[in] newStress The ArrayView holding the new stress data for each point. * @param[in] oldStress The ArrayView holding the old stress data for each point. * @param[in] disableInelasticity Flag to disable plastic response for inelastic models. @@ -68,14 +68,14 @@ class ElasticOrthotropicUpdates : public SolidBaseUpdates arrayView1d< real64 const > const & c55, arrayView1d< real64 const > const & c66, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainIncrement, - arrayView1d< real64 > const & newAnelasticStrainMagnitude, - arrayView1d< real64 > const & oldAnelasticStrainMagnitude, + arrayView2d< real64 const > const & anelasticStrainRate, + arrayView2d< real64 > const & newAnelasticStrain, + arrayView2d< real64 > const & oldAnelasticStrain, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, enableAnelasticStrain ), m_c11( c11 ), m_c12( c12 ), @@ -744,9 +744,9 @@ class ElasticOrthotropic : public SolidBase m_c55, m_c66, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, @@ -775,9 +775,9 @@ class ElasticOrthotropic : public SolidBase m_c55, m_c66, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp index 7ec996679d5..de20bb92b38 100644 --- a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp @@ -50,7 +50,7 @@ class ElasticTransverseIsotropicUpdates : public SolidBaseUpdates * @param[in] c44 The 44 component of the Voigt stiffness tensor. * @param[in] c66 The 66 component of the Voigt stiffness tensor. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. - * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. + * @param[in] anelasticStrainRate The ArrayView holding the anelastic strain rate data for each element. * @param[in] newStress The ArrayView holding the new stress data for each point. * @param[in] oldStress The ArrayView holding the old stress data for each point. * @param[in] disableInelasticity Flag to disable plastic response for inelastic models. @@ -62,14 +62,14 @@ class ElasticTransverseIsotropicUpdates : public SolidBaseUpdates arrayView1d< real64 const > const & c44, arrayView1d< real64 const > const & c66, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainIncrement, - arrayView1d< real64 > const & newAnelasticStrainMagnitude, - arrayView1d< real64 > const & oldAnelasticStrainMagnitude, + arrayView2d< real64 const > const & anelasticStrainRate, + arrayView2d< real64 > const & newAnelasticStrain, + arrayView2d< real64 > const & oldAnelasticStrain, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, enableAnelasticStrain ), m_c11( c11 ), m_c13( c13 ), @@ -600,9 +600,9 @@ class ElasticTransverseIsotropic : public SolidBase m_c44, m_c66, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, @@ -627,9 +627,9 @@ class ElasticTransverseIsotropic : public SolidBase m_c44, m_c66, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp b/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp index f771d81d836..ce1e181db8e 100644 --- a/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp +++ b/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp @@ -54,7 +54,7 @@ class ModifiedCamClayUpdates : public ElasticIsotropicPressureDependentUpdates * for each quadrature point. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. - * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. + * @param[in] anelasticStrainRate The ArrayView holding the anelastic strain rate data for each element. * @param[in] newstress The ArrayView holding the new stress data for each quadrature point. * @param[in] oldstress The ArrayView holding the old stress data from the previous converged state for each quadrature * point. @@ -70,15 +70,15 @@ class ModifiedCamClayUpdates : public ElasticIsotropicPressureDependentUpdates arrayView2d< real64 > const & oldPreConsolidationPressure, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainIncrement, - arrayView1d< real64 > const & newAnelasticStrainMagnitude, - arrayView1d< real64 > const & oldAnelasticStrainMagnitude, + arrayView2d< real64 const > const & anelasticStrainRate, + arrayView2d< real64 > const & newAnelasticStrain, + arrayView2d< real64 > const & oldAnelasticStrain, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - ElasticIsotropicPressureDependentUpdates( refPressure, refStrainVol, recompressionIndex, shearModulus, thermalExpansionCoefficient, anelasticStrainIncrement, newAnelasticStrainMagnitude, - oldAnelasticStrainMagnitude, newStress, oldStress, + ElasticIsotropicPressureDependentUpdates( refPressure, refStrainVol, recompressionIndex, shearModulus, thermalExpansionCoefficient, anelasticStrainRate, newAnelasticStrain, + oldAnelasticStrain, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), m_virginCompressionIndex( virginCompressionIndex ), m_cslSlope( cslSlope ), @@ -543,9 +543,9 @@ class ModifiedCamClay : public ElasticIsotropicPressureDependent m_oldPreConsolidationPressure, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, @@ -572,9 +572,9 @@ class ModifiedCamClay : public ElasticIsotropicPressureDependent m_oldPreConsolidationPressure, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/PerfectlyPlastic.hpp b/src/coreComponents/constitutive/solid/PerfectlyPlastic.hpp index 22afbc73b69..5e64aa3590b 100644 --- a/src/coreComponents/constitutive/solid/PerfectlyPlastic.hpp +++ b/src/coreComponents/constitutive/solid/PerfectlyPlastic.hpp @@ -48,7 +48,7 @@ class PerfectlyPlasticUpdates : public ElasticIsotropicUpdates * @param[in] bulkModulus The ArrayView holding the bulk modulus data for each element. * @param[in] shearModulus The ArrayView holding the shear modulus data for each element. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. - * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. + * @param[in] anelasticStrainRate The ArrayView holding the anelastic strain rate data for each element. * @param[in] newStress The ArrayView holding the new stress data for each quadrature point. * @param[in] oldStress The ArrayView holding the old stress data for each quadrature point. * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain @@ -57,14 +57,14 @@ class PerfectlyPlasticUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, - arrayView1d< real64 const > const & anelasticStrainIncrement, - arrayView1d< real64 > const & newAnelasticStrainMagnitude, - arrayView1d< real64 > const & oldAnelasticStrainMagnitude, + arrayView2d< real64 const > const & anelasticStrainRate, + arrayView2d< real64 > const & newAnelasticStrain, + arrayView2d< real64 > const & oldAnelasticStrain, arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, bool const & disableInelasticity, const integer & enableAnelasticStrain ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainIncrement, newAnelasticStrainMagnitude, oldAnelasticStrainMagnitude, newStress, oldStress, + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, newStress, oldStress, disableInelasticity, enableAnelasticStrain ), m_yieldStress( yieldStress ) {} @@ -288,9 +288,9 @@ class PerfectlyPlastic : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, @@ -312,9 +312,9 @@ class PerfectlyPlastic : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, - m_anelasticStrainIncrement, - m_newAnelasticStrainMagnitude, - m_oldAnelasticStrainMagnitude, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, m_disableInelasticity, diff --git a/src/coreComponents/constitutive/solid/SolidBase.cpp b/src/coreComponents/constitutive/solid/SolidBase.cpp index 6192742fe3a..26db6d1fb3f 100644 --- a/src/coreComponents/constitutive/solid/SolidBase.cpp +++ b/src/coreComponents/constitutive/solid/SolidBase.cpp @@ -40,10 +40,10 @@ SolidBase::SolidBase( string const & name, Group * const parent ): setInputFlag( InputFlags::OPTIONAL ). setDescription( "Default Linear Thermal Expansion Coefficient of the Solid Rock Frame" ); - registerWrapper( viewKeyStruct::defaultAnelasticStrainIncrementString(), &m_defaultAnelasticStrainIncrement ). - setApplyDefaultValue( 0.0 ). + registerWrapper( viewKeyStruct::defaultAnelasticStrainRateString(), &m_defaultAnelasticStrainRate ). + setApplyDefaultValue( R1Tensor{} ). setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Default anelastic strain magnitude" ); + setDescription( "Default anelastic strain rate components" ); registerWrapper( viewKeyStruct::enableAnelasticStrainString(), &m_enableAnelasticStrain ). setApplyDefaultValue( 0 ). @@ -53,6 +53,7 @@ SolidBase::SolidBase( string const & name, Group * const parent ): // register fields string const voightLabels[6] = { "XX", "YY", "ZZ", "YZ", "XZ", "XY" }; + string const normalStrainLabels[3] = { "XX", "YY", "ZZ" }; registerField< fields::solid::stress >( &m_newStress ). setDimLabels( 2, voightLabels ); @@ -64,9 +65,12 @@ SolidBase::SolidBase( string const & name, Group * const parent ): registerField< fields::solid::thermalExpansionCoefficient >( &m_thermalExpansionCoefficient ); - registerField< fields::solid::anelasticStrainIncrement >( &m_anelasticStrainIncrement ); - registerField< fields::solid::newAnelasticStrainMagnitude >( &m_newAnelasticStrainMagnitude ); - registerField< fields::solid::oldAnelasticStrainMagnitude >( &m_oldAnelasticStrainMagnitude ); + registerField< fields::solid::anelasticStrainRate >( &m_anelasticStrainRate ). + setDimLabels( 1, normalStrainLabels ); + registerField< fields::solid::newAnelasticStrain >( &m_newAnelasticStrain ). + setDimLabels( 1, normalStrainLabels ); + registerField< fields::solid::oldAnelasticStrain >( &m_oldAnelasticStrain ). + setDimLabels( 1, normalStrainLabels ); } @@ -78,12 +82,13 @@ void SolidBase::postInputInitialization() getField< fields::solid::thermalExpansionCoefficient >(). setApplyDefaultValue( m_defaultThermalExpansionCoefficient ); - getField< fields::solid::anelasticStrainIncrement >(). - setApplyDefaultValue( m_defaultAnelasticStrainIncrement ); + getField< fields::solid::anelasticStrainRate >(). + setApplyDefaultValue( 0.0 ); - GEOS_ERROR_IF( m_enableAnelasticStrain == 0 && m_defaultAnelasticStrainIncrement > 0.0, + GEOS_ERROR_IF( m_enableAnelasticStrain == 0 && + LvArray::tensorOps::l2NormSquared< 3 >( m_defaultAnelasticStrainRate ) > 0.0, getDataContext() << ": enableAnelasticStrain flag must be 1 if a nonzero" - " AnelasticStrainMagnitude is used" ); + " AnelasticStrainRate is used" ); } @@ -93,8 +98,19 @@ void SolidBase::allocateConstitutiveData( Group & parent, localIndex const numPt m_density.resize( 0, numPts ); m_newStress.resize( 0, numPts, 6 ); m_oldStress.resize( 0, numPts, 6 ); + m_anelasticStrainRate.resize( 0, 3 ); + m_newAnelasticStrain.resize( 0, 3 ); + m_oldAnelasticStrain.resize( 0, 3 ); ConstitutiveBase::allocateConstitutiveData( parent, numPts ); + + for( localIndex ei = 0; ei < parent.size(); ++ei ) + { + for( integer dim = 0; dim < 3; ++dim ) + { + m_anelasticStrainRate[ei][dim] = m_defaultAnelasticStrainRate[dim]; + } + } } @@ -106,12 +122,12 @@ void SolidBase::saveConvergedState() const arrayView3d< real64 const, solid::STRESS_USD > newStress = m_newStress; arrayView3d< real64, solid::STRESS_USD > oldStress = m_oldStress; - arrayView1d< real64 const > newAnelasticStrainMagnitude = m_newAnelasticStrainMagnitude; - arrayView1d< real64 > oldAnelasticStrainMagnitude = m_oldAnelasticStrainMagnitude; + arrayView2d< real64 const > newAnelasticStrain = m_newAnelasticStrain; + arrayView2d< real64 > oldAnelasticStrain = m_oldAnelasticStrain; forAll< parallelDevicePolicy<> >( numE, [=] GEOS_HOST_DEVICE ( localIndex const k ) { - oldAnelasticStrainMagnitude[k] = newAnelasticStrainMagnitude[k]; + LvArray::tensorOps::copy< 3 >( oldAnelasticStrain[k], newAnelasticStrain[k] ); for( localIndex q = 0; q < numQ; ++q ) { diff --git a/src/coreComponents/constitutive/solid/SolidBase.hpp b/src/coreComponents/constitutive/solid/SolidBase.hpp index 9ba9bf3e4a5..a0e81e4e3c6 100644 --- a/src/coreComponents/constitutive/solid/SolidBase.hpp +++ b/src/coreComponents/constitutive/solid/SolidBase.hpp @@ -55,24 +55,24 @@ class SolidBaseUpdates * @param[in] oldStress The old stress data from the constitutive model class. * @param[in] thermalExpansionCoefficient The ArrayView holding the thermal expansion coefficient data for each element. * @param[in] disableInelasticity Flag to disable inelastic response - * @param[in] anelasticStrainMagnitude The ArrayView holding the anelastic strain magnitude data for each element. + * @param[in] anelasticStrainRate The ArrayView holding the anelastic strain rate data for each element. * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain */ SolidBaseUpdates( arrayView3d< real64, solid::STRESS_USD > const & newStress, arrayView3d< real64, solid::STRESS_USD > const & oldStress, arrayView1d< real64 const > const & thermalExpansionCoefficient, const bool & disableInelasticity, - arrayView1d< real64 const > const & anelasticStrainIncrement, - arrayView1d< real64 > const & newAnelasticStrainMagnitude, - arrayView1d< real64 > const & oldAnelasticStrainMagnitude, + arrayView2d< real64 const > const & anelasticStrainRate, + arrayView2d< real64 > const & newAnelasticStrain, + arrayView2d< real64 > const & oldAnelasticStrain, const integer enableAnelasticStrain ): m_newStress( newStress ), m_oldStress( oldStress ), m_thermalExpansionCoefficient( thermalExpansionCoefficient ), m_disableInelasticity ( disableInelasticity ), - m_anelasticStrainIncrement( anelasticStrainIncrement ), - m_newAnelasticStrainMagnitude( newAnelasticStrainMagnitude ), - m_oldAnelasticStrainMagnitude( oldAnelasticStrainMagnitude ), + m_anelasticStrainRate( anelasticStrainRate ), + m_newAnelasticStrain( newAnelasticStrain ), + m_oldAnelasticStrain( oldAnelasticStrain ), m_enableAnelasticStrain( enableAnelasticStrain ) {} @@ -111,11 +111,14 @@ class SolidBaseUpdates /// Flag to disable inelasticity const bool m_disableInelasticity; - /// The anelastic strain magnitude (i.e. chemistry, electrochemistry, etc.) - arrayView1d< real64 const > const m_anelasticStrainIncrement; + /// The anelastic strain rate (i.e. chemistry, electrochemistry, etc.) + arrayView2d< real64 const > const m_anelasticStrainRate; - arrayView1d< real64 > const m_newAnelasticStrainMagnitude; - arrayView1d< real64 > const m_oldAnelasticStrainMagnitude; + /// The current accumulated anelastic strain components. + arrayView2d< real64 > const m_newAnelasticStrain; + + /// The previous accumulated anelastic strain components. + arrayView2d< real64 > const m_oldAnelasticStrain; /// Flag to enable stress modification due to anelastic strain const integer m_enableAnelasticStrain; @@ -146,14 +149,14 @@ class SolidBaseUpdates } /** - * @brief Get anelasticStrainIncrement + * @brief Get anelasticStrainRate * @param[in] k Element index. - * @return the anelasticStrainIncrement of element k + * @return the anelasticStrainRate of element k */ GEOS_HOST_DEVICE - real64 getAnelasticStrainIncrement( localIndex const k ) const + arraySlice1d< real64 const > getAnelasticStrainRate( localIndex const k ) const { - return m_anelasticStrainIncrement[k]; + return m_anelasticStrainRate[k]; } /** @@ -404,15 +407,18 @@ class SolidBaseUpdates /** * @brief Calculate the stress modifier due to anelastic strain * + * @param[in] timeIncrement time increment used to integrate the anelastic strain rate * @param[out] stress New stress value (Cauchy stress) */ GEOS_HOST_DEVICE virtual void stressModificationByAnelasticStain( localIndex const k, localIndex const q, + real64 const & timeIncrement, real64 ( & stressModifier )[6] ) const { GEOS_UNUSED_VAR( k ); GEOS_UNUSED_VAR( q ); + GEOS_UNUSED_VAR( timeIncrement ); GEOS_UNUSED_VAR( stressModifier ); GEOS_ERROR( "stressModificationByAnelasticStain() not implemented for this model" ); } @@ -617,8 +623,8 @@ class SolidBase : public constitutive::ConstitutiveBase static constexpr char const * defaultDensityString() { return "defaultDensity"; } // Default drained linear thermal expansion coefficient key static constexpr char const * defaultThermalExpansionCoefficientString() { return "defaultDrainedLinearTEC"; } - // Default anelastic strain magnitude key - static constexpr char const * defaultAnelasticStrainIncrementString() { return "defaultAnelasticStrainIncrement"; } + // Default anelastic strain rate key + static constexpr char const * defaultAnelasticStrainRateString() { return "defaultAnelasticStrainRate"; } // Enable stress modification due to anelastic strain key static constexpr char const * enableAnelasticStrainString() { return "enableAnelasticStrain"; } }; @@ -749,14 +755,17 @@ class SolidBase : public constitutive::ConstitutiveBase /// Flag to disable inelasticity (plasticity, damage, etc.) bool m_disableInelasticity = false; - /// The anelastic strain rate magnitude (i.e. chemistry, electrochemistry, etc.) - array1d< real64 > m_anelasticStrainIncrement; + /// The anelastic strain rate (i.e. chemistry, electrochemistry, etc.) + array2d< real64 > m_anelasticStrainRate; + + /// The current accumulated anelastic strain components. + array2d< real64 > m_newAnelasticStrain; - array1d< real64 > m_newAnelasticStrainMagnitude; - array1d< real64 > m_oldAnelasticStrainMagnitude; + /// The previous accumulated anelastic strain components. + array2d< real64 > m_oldAnelasticStrain; - /// The default value of the anelastic strain rate magnitude - real64 m_defaultAnelasticStrainIncrement; + /// The default value of the anelastic strain rate + R1Tensor m_defaultAnelasticStrainRate = { 0.0, 0.0, 0.0 }; /// Flag to enable stress modification due to anelastic strain integer m_enableAnelasticStrain; diff --git a/src/coreComponents/constitutive/solid/SolidFields.hpp b/src/coreComponents/constitutive/solid/SolidFields.hpp index fd6670079e7..f680263e115 100644 --- a/src/coreComponents/constitutive/solid/SolidFields.hpp +++ b/src/coreComponents/constitutive/solid/SolidFields.hpp @@ -436,29 +436,29 @@ DECLARE_FIELD( dInternalEnergy_dTemperature, WRITE_AND_READ, "Derivative of the solid internal energy w.r.t. temperature [J/(m^3.K)]" ); -DECLARE_FIELD( anelasticStrainIncrement, - "anelasticStrainIncrement", - array1d< real64 >, +DECLARE_FIELD( anelasticStrainRate, + "anelasticStrainRate", + array2d< real64 >, 0, LEVEL_0, WRITE_AND_READ, - "Anelastic strain increment (i.e. chemistry, electrochemistry, etc.)" ); + "Anelastic strain rate (i.e. chemistry, electrochemistry, etc.)" ); -DECLARE_FIELD( newAnelasticStrainMagnitude, - "newAnelasticStrainMagnitude", - array1d< real64 >, +DECLARE_FIELD( newAnelasticStrain, + "newAnelasticStrain", + array2d< real64 >, 0, LEVEL_0, WRITE_AND_READ, - "New anelastic strain magnitude" ); + "New anelastic strain components" ); -DECLARE_FIELD( oldAnelasticStrainMagnitude, - "oldAnelasticStrainMagnitude", - array1d< real64 >, +DECLARE_FIELD( oldAnelasticStrain, + "oldAnelasticStrain", + array2d< real64 >, 0, NOPLOT, WRITE_AND_READ, - "Old anelastic strain magnitude" ); + "Old anelastic strain components" ); } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic_impl.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic_impl.hpp index f6119c20fb2..82b5ecd2bcd 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic_impl.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic_impl.hpp @@ -129,7 +129,7 @@ void ImplicitSmallStrainQuasiStatic< SUBREGION_TYPE, CONSTITUTIVE_TYPE, FE_TYPE finiteElement::feOps::symmetricGradient( dNdX, stack.uhat_local, strainInc ); m_constitutiveUpdate.smallStrainUpdate( k, q, m_dt, strainInc, stress, stiffness ); - m_constitutiveUpdate.stressModificationByAnelasticStain( k, q, stressModifierAnelasticStrain ); + m_constitutiveUpdate.stressModificationByAnelasticStain( k, q, m_dt, stressModifierAnelasticStrain ); stressModifier( stress, stressModifierAnelasticStrain ); // #pragma unroll From 9eea5517c5668dbf5553c169ff6ea88114efbd03 Mon Sep 17 00:00:00 2001 From: Fan Fei Date: Mon, 3 Aug 2026 20:09:05 -0700 Subject: [PATCH 4/4] added a xml --- .../PhaseFieldFracture_AnelasticStrain.xml | 214 ++++++++++++++++++ 1 file changed, 214 insertions(+) create mode 100644 inputFiles/phaseField/PhaseFieldFracture_AnelasticStrain.xml diff --git a/inputFiles/phaseField/PhaseFieldFracture_AnelasticStrain.xml b/inputFiles/phaseField/PhaseFieldFracture_AnelasticStrain.xml new file mode 100644 index 00000000000..10578450614 --- /dev/null +++ b/inputFiles/phaseField/PhaseFieldFracture_AnelasticStrain.xml @@ -0,0 +1,214 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +