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 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/constitutive/solid/CeramicDamage.hpp b/src/coreComponents/constitutive/solid/CeramicDamage.hpp index ae1a9775be6..f2fd9012e60 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] 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 + * @param[in] enableAnelasticStrain Flag to enable stress modification due to anelastic strain */ CeramicDamageUpdates( arrayView2d< real64 > const & damage, arrayView2d< real64 > const & jacobian, @@ -80,10 +83,15 @@ class CeramicDamageUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, + 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 ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, newStress, oldStress, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, newStress, oldStress, + disableInelasticity, enableAnelasticStrain ), m_damage( damage ), m_jacobian( jacobian ), m_lengthScale( lengthScale ), @@ -472,9 +480,13 @@ class CeramicDamage : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -498,9 +510,13 @@ class CeramicDamage : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, 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..25b51d10da1 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] 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 + * @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,15 @@ class DelftEggUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, + 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 ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, newStress, oldStress, disableInelasticity ), + const bool & disableInelasticity, + const integer & enableAnelasticStrain ): + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, newStress, oldStress, + disableInelasticity, enableAnelasticStrain ), m_recompressionIndex( recompressionIndex ), m_virginCompressionIndex( virginCompressionIndex ), m_cslSlope( cslSlope ), @@ -524,9 +531,13 @@ class DelftEgg : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -549,9 +560,13 @@ class DelftEgg : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, 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..caca60db597 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] 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 + * @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,15 @@ class DruckerPragerUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, + 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 ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, newStress, oldStress, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, newStress, oldStress, + disableInelasticity, enableAnelasticStrain ), m_friction( friction ), m_dilation( dilation ), m_hardening( hardening ), @@ -410,9 +418,13 @@ class DruckerPrager : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -434,9 +446,13 @@ class DruckerPrager : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, 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..5aa94492fc6 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] 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 */ DruckerPragerExtendedUpdates( arrayView1d< real64 const > const & initialFriction, arrayView1d< real64 const > const & residualFriction, @@ -58,10 +60,15 @@ class DruckerPragerExtendedUpdates : public ElasticIsotropicUpdates arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, + 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 ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, newStress, oldStress, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, newStress, oldStress, + disableInelasticity, enableAnelasticStrain ), m_initialFriction( initialFriction ), m_residualFriction( residualFriction ), m_dilationRatio( dilationRatio ), @@ -441,9 +448,13 @@ class DruckerPragerExtended : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -467,9 +478,13 @@ class DruckerPragerExtended : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, 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 79c186a5c06..58c966354d0 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp @@ -46,17 +46,24 @@ 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] 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 + * @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, + 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 ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity ), + const bool & disableInelasticity, + const integer & enableAnelasticStrain ): + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, + enableAnelasticStrain ), m_bulkModulus( bulkModulus ), m_shearModulus( shearModulus ) {} @@ -157,6 +164,12 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates localIndex const q, real64 beta ) const override; + 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 /* @@ -374,6 +387,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 const & timeIncrement, + real64 ( & stressModifier )[6] ) const +{ + if( m_enableAnelasticStrain == 0 ) + { + return; + } + + arraySlice1d< real64 const > const anelasticStrainRate = getAnelasticStrainRate( k ); + real64 anelasticStrain[6] = {}; + for( integer i = 0; i < 3; ++i ) + { + m_newAnelasticStrain[k][i] = m_oldAnelasticStrain[k][i] + anelasticStrainRate[i] * timeIncrement; + anelasticStrain[i] = m_newAnelasticStrain[k][i]; + } + + smallStrainNoStateUpdate_StressOnly( k, q, anelasticStrain, stressModifier ); +} + // TODO: need to confirm stress / strain measures before activating hyper inferface /* @@ -533,18 +569,26 @@ class ElasticIsotropic : public SolidBase return ElasticIsotropicUpdates( m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, 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_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, arrayView3d< real64, solid::STRESS_USD >(), arrayView3d< real64, solid::STRESS_USD >(), - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } } @@ -563,9 +607,13 @@ class ElasticIsotropic : public SolidBase m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, 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..3176f714ee6 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp @@ -49,19 +49,26 @@ 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] 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 + * @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, + 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 ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, + enableAnelasticStrain ), m_refPressure( refPressure ), m_refStrainVol( refStrainVol ), m_recompressionIndex( recompressionIndex ), @@ -583,9 +590,13 @@ class ElasticIsotropicPressureDependent : public SolidBase m_recompressionIndex, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, 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 +605,13 @@ class ElasticIsotropicPressureDependent : public SolidBase m_recompressionIndex, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, arrayView3d< real64, solid::STRESS_USD >(), arrayView3d< real64, solid::STRESS_USD >(), - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } } @@ -617,9 +632,13 @@ class ElasticIsotropicPressureDependent : public SolidBase m_recompressionIndex, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, 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..0498f18b108 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] 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. + * @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,15 @@ class ElasticOrthotropicUpdates : public SolidBaseUpdates arrayView1d< real64 const > const & c55, arrayView1d< real64 const > const & c66, arrayView1d< real64 const > const & thermalExpansionCoefficient, + 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 ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, + enableAnelasticStrain ), m_c11( c11 ), m_c12( c12 ), m_c13( c13 ), @@ -737,9 +744,13 @@ class ElasticOrthotropic : public SolidBase m_c55, m_c66, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -764,9 +775,13 @@ class ElasticOrthotropic : public SolidBase m_c55, m_c66, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, 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..de20bb92b38 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] 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. + * @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,15 @@ class ElasticTransverseIsotropicUpdates : public SolidBaseUpdates arrayView1d< real64 const > const & c44, arrayView1d< real64 const > const & c66, arrayView1d< real64 const > const & thermalExpansionCoefficient, + 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 ): - SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + SolidBaseUpdates( newStress, oldStress, thermalExpansionCoefficient, disableInelasticity, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, + enableAnelasticStrain ), m_c11( c11 ), m_c13( c13 ), m_c33( c33 ), @@ -593,9 +600,13 @@ class ElasticTransverseIsotropic : public SolidBase m_c44, m_c66, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -616,9 +627,13 @@ class ElasticTransverseIsotropic : public SolidBase m_c44, m_c66, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, 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..ce1e181db8e 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] 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. * @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,16 @@ class ModifiedCamClayUpdates : public ElasticIsotropicPressureDependentUpdates arrayView2d< real64 > const & oldPreConsolidationPressure, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, + 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 ): - ElasticIsotropicPressureDependentUpdates( refPressure, refStrainVol, recompressionIndex, shearModulus, thermalExpansionCoefficient, newStress, oldStress, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + ElasticIsotropicPressureDependentUpdates( refPressure, refStrainVol, recompressionIndex, shearModulus, thermalExpansionCoefficient, anelasticStrainRate, newAnelasticStrain, + oldAnelasticStrain, newStress, oldStress, + disableInelasticity, enableAnelasticStrain ), m_virginCompressionIndex( virginCompressionIndex ), m_cslSlope( cslSlope ), m_newPreConsolidationPressure( newPreConsolidationPressure ), @@ -535,9 +543,13 @@ class ModifiedCamClay : public ElasticIsotropicPressureDependent m_oldPreConsolidationPressure, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -560,9 +572,13 @@ class ModifiedCamClay : public ElasticIsotropicPressureDependent m_oldPreConsolidationPressure, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, 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..5e64aa3590b 100644 --- a/src/coreComponents/constitutive/solid/PerfectlyPlastic.hpp +++ b/src/coreComponents/constitutive/solid/PerfectlyPlastic.hpp @@ -48,17 +48,24 @@ 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] 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 */ PerfectlyPlasticUpdates( arrayView1d< real64 const > const & yieldStress, arrayView1d< real64 const > const & bulkModulus, arrayView1d< real64 const > const & shearModulus, arrayView1d< real64 const > const & thermalExpansionCoefficient, + 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 ): - ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, newStress, oldStress, disableInelasticity ), + bool const & disableInelasticity, + const integer & enableAnelasticStrain ): + ElasticIsotropicUpdates( bulkModulus, shearModulus, thermalExpansionCoefficient, anelasticStrainRate, newAnelasticStrain, oldAnelasticStrain, newStress, oldStress, + disableInelasticity, enableAnelasticStrain ), m_yieldStress( yieldStress ) {} @@ -281,9 +288,13 @@ class PerfectlyPlastic : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, m_newStress, m_oldStress, - m_disableInelasticity ); + m_disableInelasticity, + m_enableAnelasticStrain ); } /** @@ -301,9 +312,13 @@ class PerfectlyPlastic : public ElasticIsotropic m_bulkModulus, m_shearModulus, m_thermalExpansionCoefficient, + m_anelasticStrainRate, + m_newAnelasticStrain, + m_oldAnelasticStrain, 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..26db6d1fb3f 100644 --- a/src/coreComponents/constitutive/solid/SolidBase.cpp +++ b/src/coreComponents/constitutive/solid/SolidBase.cpp @@ -40,9 +40,20 @@ SolidBase::SolidBase( string const & name, Group * const parent ): setInputFlag( InputFlags::OPTIONAL ). setDescription( "Default Linear Thermal Expansion Coefficient of the Solid Rock Frame" ); + registerWrapper( viewKeyStruct::defaultAnelasticStrainRateString(), &m_defaultAnelasticStrainRate ). + setApplyDefaultValue( R1Tensor{} ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Default anelastic strain rate components" ); + + 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" }; + string const normalStrainLabels[3] = { "XX", "YY", "ZZ" }; registerField< fields::solid::stress >( &m_newStress ). setDimLabels( 2, voightLabels ); @@ -53,6 +64,13 @@ SolidBase::SolidBase( string const & name, Group * const parent ): registerField< fields::solid::density >( &m_density ); registerField< fields::solid::thermalExpansionCoefficient >( &m_thermalExpansionCoefficient ); + + 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 ); } @@ -63,6 +81,14 @@ void SolidBase::postInputInitialization() getField< fields::solid::thermalExpansionCoefficient >(). setApplyDefaultValue( m_defaultThermalExpansionCoefficient ); + + getField< fields::solid::anelasticStrainRate >(). + setApplyDefaultValue( 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" + " AnelasticStrainRate is used" ); } @@ -72,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]; + } + } } @@ -85,8 +122,13 @@ void SolidBase::saveConvergedState() const arrayView3d< real64 const, solid::STRESS_USD > newStress = m_newStress; arrayView3d< real64, solid::STRESS_USD > oldStress = m_oldStress; + arrayView2d< real64 const > newAnelasticStrain = m_newAnelasticStrain; + arrayView2d< real64 > oldAnelasticStrain = m_oldAnelasticStrain; + forAll< parallelDevicePolicy<> >( numE, [=] GEOS_HOST_DEVICE ( localIndex const k ) { + LvArray::tensorOps::copy< 3 >( oldAnelasticStrain[k], newAnelasticStrain[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 f21092d7a61..a0e81e4e3c6 100644 --- a/src/coreComponents/constitutive/solid/SolidBase.hpp +++ b/src/coreComponents/constitutive/solid/SolidBase.hpp @@ -55,15 +55,25 @@ 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] 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 ): + const bool & disableInelasticity, + 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_disableInelasticity ( disableInelasticity ), + m_anelasticStrainRate( anelasticStrainRate ), + m_newAnelasticStrain( newAnelasticStrain ), + m_oldAnelasticStrain( oldAnelasticStrain ), + m_enableAnelasticStrain( enableAnelasticStrain ) {} /// Deleted default constructor @@ -101,6 +111,18 @@ class SolidBaseUpdates /// Flag to disable inelasticity const bool m_disableInelasticity; + /// The anelastic strain rate (i.e. chemistry, electrochemistry, etc.) + arrayView2d< real64 const > const m_anelasticStrainRate; + + /// 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; + /** * @brief Get bulkModulus * @param[in] k Element index. @@ -126,6 +148,17 @@ class SolidBaseUpdates return m_thermalExpansionCoefficient[k]; } + /** + * @brief Get anelasticStrainRate + * @param[in] k Element index. + * @return the anelasticStrainRate of element k + */ + GEOS_HOST_DEVICE + arraySlice1d< real64 const > getAnelasticStrainRate( localIndex const k ) const + { + return m_anelasticStrainRate[k]; + } + /** * @brief Get shear modulus * @param[in] k Element index. @@ -371,6 +404,25 @@ class SolidBaseUpdates GEOS_ERROR( "viscousStateUpdate() not implemented for this model" ); } + /** + * @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" ); + } + /** * @brief Return the strain energy density at a given material point * @@ -571,6 +623,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 rate key + static constexpr char const * defaultAnelasticStrainRateString() { return "defaultAnelasticStrainRate"; } + // 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 +754,22 @@ class SolidBase : public constitutive::ConstitutiveBase /// Flag to disable inelasticity (plasticity, damage, etc.) bool m_disableInelasticity = false; + + /// The anelastic strain rate (i.e. chemistry, electrochemistry, etc.) + array2d< real64 > m_anelasticStrainRate; + + /// The current accumulated anelastic strain components. + array2d< real64 > m_newAnelasticStrain; + + /// The previous accumulated anelastic strain components. + array2d< real64 > m_oldAnelasticStrain; + + /// 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; + }; } // namespace constitutive diff --git a/src/coreComponents/constitutive/solid/SolidFields.hpp b/src/coreComponents/constitutive/solid/SolidFields.hpp index 8ed466eadc7..f680263e115 100644 --- a/src/coreComponents/constitutive/solid/SolidFields.hpp +++ b/src/coreComponents/constitutive/solid/SolidFields.hpp @@ -436,6 +436,30 @@ DECLARE_FIELD( dInternalEnergy_dTemperature, WRITE_AND_READ, "Derivative of the solid internal energy w.r.t. temperature [J/(m^3.K)]" ); +DECLARE_FIELD( anelasticStrainRate, + "anelasticStrainRate", + array2d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "Anelastic strain rate (i.e. chemistry, electrochemistry, etc.)" ); + +DECLARE_FIELD( newAnelasticStrain, + "newAnelasticStrain", + array2d< real64 >, + 0, + LEVEL_0, + WRITE_AND_READ, + "New anelastic strain components" ); + +DECLARE_FIELD( oldAnelasticStrain, + "oldAnelasticStrain", + array2d< real64 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Old anelastic strain components" ); + } } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic.hpp index b5f1d5a05d6..1fd1f03d95c 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic.hpp @@ -179,11 +179,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 1acbcf75f78..82b5ecd2bcd 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic_impl.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/ImplicitSmallStrainQuasiStatic_impl.hpp @@ -122,14 +122,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; finiteElement::feOps::symmetricGradient( dNdX, stack.uhat_local, strainInc ); m_constitutiveUpdate.smallStrainUpdate( k, q, m_dt, strainInc, stress, stiffness ); + m_constitutiveUpdate.stressModificationByAnelasticStain( k, q, m_dt, stressModifierAnelasticStrain ); - stressModifier( stress ); + stressModifier( stress, stressModifierAnelasticStrain ); // #pragma unroll for( localIndex i=0; i<6; ++i ) {