diff --git a/inputFiles/thermoPoromechanics/ThermoDruckerPrager_1DCooling_smoke.xml b/inputFiles/thermoPoromechanics/ThermoDruckerPrager_1DCooling_smoke.xml new file mode 100644 index 00000000000..51044274748 --- /dev/null +++ b/inputFiles/thermoPoromechanics/ThermoDruckerPrager_1DCooling_smoke.xml @@ -0,0 +1,119 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/thermoPoromechanics/ThermoElastic_1DCooling_smoke.xml b/inputFiles/thermoPoromechanics/ThermoElastic_1DCooling_smoke.xml new file mode 100644 index 00000000000..c16357117d3 --- /dev/null +++ b/inputFiles/thermoPoromechanics/ThermoElastic_1DCooling_smoke.xml @@ -0,0 +1,115 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/thermoPoromechanics/ThermoMech_1DCooling_base.xml b/inputFiles/thermoPoromechanics/ThermoMech_1DCooling_base.xml new file mode 100644 index 00000000000..0243f448fe4 --- /dev/null +++ b/inputFiles/thermoPoromechanics/ThermoMech_1DCooling_base.xml @@ -0,0 +1,155 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/thermoPoromechanics/thermoPoromechanics.ats b/inputFiles/thermoPoromechanics/thermoPoromechanics.ats index 443edf0817f..d11c5a236d2 100644 --- a/inputFiles/thermoPoromechanics/thermoPoromechanics.ats +++ b/inputFiles/thermoPoromechanics/thermoPoromechanics.ats @@ -65,6 +65,22 @@ decks = [ restart_step=0, check_step=1, restartcheck_params=RestartcheckParameters(**restartcheck_params)) + TestDeck( + name="ThermoDruckerPrager_1DCooling_smoke", + description= + '1D thermomechanics problem under cooling with Drucker Prager model', + partitions=((1, 1, 1), (1, 2, 1)), + restart_step=0, + check_step=1, + restartcheck_params=RestartcheckParameters(**restartcheck_params)) + TestDeck( + name="ThermoElastic_1DCooling_smoke", + description= + '1D thermomechanics problem under cooling with linear elastic isotropic model', + partitions=((1, 1, 1), (1, 2, 1)), + restart_step=0, + check_step=1, + restartcheck_params=RestartcheckParameters(**restartcheck_params)) ] generate_geos_tests(decks) diff --git a/src/coreComponents/constitutive/solid/Damage.hpp b/src/coreComponents/constitutive/solid/Damage.hpp index b54bb30f372..f89980206bd 100644 --- a/src/coreComponents/constitutive/solid/Damage.hpp +++ b/src/coreComponents/constitutive/solid/Damage.hpp @@ -108,6 +108,7 @@ class DamageUpdates : public UPDATE_BASE using UPDATE_BASE::smallStrainNoStateUpdate; using UPDATE_BASE::smallStrainUpdate; + using UPDATE_BASE::smallStrainUpdate_thermal; using UPDATE_BASE::smallStrainNoStateUpdate_StressOnly; using UPDATE_BASE::smallStrainUpdate_StressOnly; using UPDATE_BASE::saveConvergedState; @@ -293,6 +294,26 @@ class DamageUpdates : public UPDATE_BASE stiffness.scaleParams( factor ); } + GEOS_HOST_DEVICE + void smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( & stress )[6], + DiscretizationOps & stiffness, + real64 ( & dStress_dTemperature )[6] ) const override + { + this->smallStrainUpdate( k, q, timeIncrement, strainIncrementNoThermalStrain, stress, stiffness ); + + real64 const bulkMod = stiffness.m_bulkModulus; + + dStress_dTemperature[0] = 3 * bulkMod; + dStress_dTemperature[1] = 3 * bulkMod; + dStress_dTemperature[2] = 3 * bulkMod; + dStress_dTemperature[3] = 0; + dStress_dTemperature[4] = 0; + dStress_dTemperature[5] = 0; + } // TODO: The code below assumes the strain energy density will never be // evaluated in a non-converged / garbage configuration. diff --git a/src/coreComponents/constitutive/solid/DamageSpectral.hpp b/src/coreComponents/constitutive/solid/DamageSpectral.hpp index 03a5e3bff73..84f53c61e59 100644 --- a/src/coreComponents/constitutive/solid/DamageSpectral.hpp +++ b/src/coreComponents/constitutive/solid/DamageSpectral.hpp @@ -64,6 +64,7 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > using DiscretizationOps = SolidModelDiscretizationOpsFullyAnisotropic; using DamageUpdates< UPDATE_BASE >::smallStrainUpdate; + using DamageUpdates< UPDATE_BASE >::smallStrainUpdate_thermal; using DamageUpdates< UPDATE_BASE >::saveConvergedState; using DamageUpdates< UPDATE_BASE >::getDegradationValue; @@ -268,6 +269,28 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > smallStrainUpdate( k, q, timeIncrement, strainIncrement, stress, stiffness.m_c ); } + GEOS_HOST_DEVICE + virtual void smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( & stress )[6], + DiscretizationOps & stiffness, + real64 ( & dStress_dTemperature )[6] ) const final + { + smallStrainUpdate( k, q, timeIncrement, strainIncrementNoThermalStrain, stress, stiffness ); + + real64 const ( &stiffnessTensor )[6][6] = stiffness.m_c; + + for( int i=0; i<6; ++i ) + { + for( int j=0; j<3; ++j ) + { + dStress_dTemperature[i] += stiffnessTensor[i][j]; + } + } + } + GEOS_HOST_DEVICE virtual real64 getStrainEnergyDensity( localIndex const k, diff --git a/src/coreComponents/constitutive/solid/DamageVolDev.hpp b/src/coreComponents/constitutive/solid/DamageVolDev.hpp index 951eded451f..6a88556d2c9 100644 --- a/src/coreComponents/constitutive/solid/DamageVolDev.hpp +++ b/src/coreComponents/constitutive/solid/DamageVolDev.hpp @@ -60,6 +60,7 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > using DiscretizationOps = typename DamageUpdates< UPDATE_BASE >::DiscretizationOps; using DamageUpdates< UPDATE_BASE >::smallStrainUpdate; + using DamageUpdates< UPDATE_BASE >::smallStrainUpdate_thermal; using DamageUpdates< UPDATE_BASE >::saveConvergedState; using DamageUpdates< UPDATE_BASE >::getDegradationValue; @@ -153,6 +154,26 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > } } + GEOS_HOST_DEVICE + virtual void smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( & stress )[6], + DiscretizationOps & stiffness, + real64 ( & dStress_dTemperature )[6] ) const final + { + smallStrainUpdate( k, q, timeIncrement, strainIncrementNoThermalStrain, stress, stiffness ); + + real64 const bulkMod = stiffness.m_bulkModulus; + + dStress_dTemperature[0] = 3 * bulkMod; + dStress_dTemperature[1] = 3 * bulkMod; + dStress_dTemperature[2] = 3 * bulkMod; + dStress_dTemperature[3] = 0; + dStress_dTemperature[4] = 0; + dStress_dTemperature[5] = 0; + } GEOS_HOST_DEVICE virtual real64 getStrainEnergyDensity( localIndex const k, diff --git a/src/coreComponents/constitutive/solid/DelftEgg.hpp b/src/coreComponents/constitutive/solid/DelftEgg.hpp index 167301468f0..14746aa9e35 100644 --- a/src/coreComponents/constitutive/solid/DelftEgg.hpp +++ b/src/coreComponents/constitutive/solid/DelftEgg.hpp @@ -99,6 +99,7 @@ class DelftEggUpdates : public ElasticIsotropicUpdates // Bring in base implementations to prevent hiding warnings using ElasticIsotropicUpdates::smallStrainUpdate; + using ElasticIsotropicUpdates::smallStrainUpdate_thermal; GEOS_HOST_DEVICE void evaluateYield( real64 const p, @@ -135,6 +136,15 @@ class DelftEggUpdates : public ElasticIsotropicUpdates real64 ( &stress )[6], DiscretizationOps & stiffness ) const final; + GEOS_HOST_DEVICE + virtual void smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( &stress )[6], + DiscretizationOps & stiffness, + real64 ( &dStress_dTemperature )[6] ) const final; + GEOS_HOST_DEVICE inline virtual void saveConvergedState( localIndex const k, @@ -449,7 +459,28 @@ void DelftEggUpdates::smallStrainUpdate( localIndex const k, smallStrainUpdate( k, q, timeIncrement, strainIncrement, stress, stiffness.m_c ); } +GEOS_HOST_DEVICE +inline +void DelftEggUpdates::smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( & stress )[6], + DiscretizationOps & stiffness, + real64 ( & dStress_dTemperature )[6] ) const +{ + smallStrainUpdate( k, q, timeIncrement, strainIncrementNoThermalStrain, stress, stiffness ); + + real64 const ( &stiffnessTensor )[6][6] = stiffness.m_c; + for( int i=0; i<6; ++i ) + { + for( int j=0; j<3; ++j ) + { + dStress_dTemperature[i] += stiffnessTensor[i][j]; + } + } +} /** * @class DelftEgg diff --git a/src/coreComponents/constitutive/solid/DruckerPrager.hpp b/src/coreComponents/constitutive/solid/DruckerPrager.hpp index e59494d7b06..9241ac7f22f 100644 --- a/src/coreComponents/constitutive/solid/DruckerPrager.hpp +++ b/src/coreComponents/constitutive/solid/DruckerPrager.hpp @@ -95,6 +95,7 @@ class DruckerPragerUpdates : public ElasticIsotropicUpdates // Bring in base implementations to prevent hiding warnings using ElasticIsotropicUpdates::smallStrainUpdate; + using ElasticIsotropicUpdates::smallStrainUpdate_thermal; GEOS_HOST_DEVICE void smallStrainUpdate( localIndex const k, @@ -120,6 +121,15 @@ class DruckerPragerUpdates : public ElasticIsotropicUpdates real64 ( &stress )[6], real64 ( &stiffness )[6][6] ) const override; + GEOS_HOST_DEVICE + virtual void smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( &stress )[6], + DiscretizationOps & stiffness, + real64 ( &dStress_dTemperature )[6] ) const; + GEOS_HOST_DEVICE inline virtual void saveConvergedState( localIndex const k, @@ -339,7 +349,28 @@ void DruckerPragerUpdates::smallStrainUpdate( localIndex const k, smallStrainUpdate( k, q, timeIncrement, strainIncrement, stress, stiffness.m_c ); } +GEOS_HOST_DEVICE +inline +void DruckerPragerUpdates::smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( & stress )[6], + DiscretizationOps & stiffness, + real64 ( & dStress_dTemperature )[6] ) const +{ + smallStrainUpdate( k, q, timeIncrement, strainIncrementNoThermalStrain, stress, stiffness ); + + real64 const ( &stiffnessTensor )[6][6] = stiffness.m_c; + for( int i=0; i<6; ++i ) + { + for( int j=0; j<3; ++j ) + { + dStress_dTemperature[i] += stiffnessTensor[i][j]; + } + } +} /** * @class DruckerPrager diff --git a/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp b/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp index dc79394677f..cdc98fc8e13 100644 --- a/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp +++ b/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp @@ -91,6 +91,7 @@ class DruckerPragerExtendedUpdates : public ElasticIsotropicUpdates // Bring in base implementations to prevent hiding warnings using ElasticIsotropicUpdates::smallStrainUpdate; + using ElasticIsotropicUpdates::smallStrainUpdate_thermal; GEOS_HOST_DEVICE void smallStrainUpdate( localIndex const k, @@ -116,6 +117,14 @@ class DruckerPragerExtendedUpdates : public ElasticIsotropicUpdates real64 ( &stress )[6], real64 ( &stiffness )[6][6] ) const override; + GEOS_HOST_DEVICE + virtual void smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( &stress )[6], + DiscretizationOps & stiffness, + real64 ( &dStress_dTemperature )[6] ) const; GEOS_HOST_DEVICE inline @@ -365,7 +374,28 @@ void DruckerPragerExtendedUpdates::smallStrainUpdate( localIndex const k, smallStrainUpdate( k, q, timeIncrement, strainIncrement, stress, stiffness.m_c ); } +GEOS_HOST_DEVICE +inline +void DruckerPragerExtendedUpdates::smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( & stress )[6], + DiscretizationOps & stiffness, + real64 ( & dStress_dTemperature )[6] ) const +{ + smallStrainUpdate( k, q, timeIncrement, strainIncrementNoThermalStrain, stress, stiffness ); + + real64 const ( &stiffnessTensor )[6][6] = stiffness.m_c; + for( int i=0; i<6; ++i ) + { + for( int j=0; j<3; ++j ) + { + dStress_dTemperature[i] += stiffnessTensor[i][j]; + } + } +} /** * @class DruckerPragerExtended diff --git a/src/coreComponents/constitutive/solid/DuvautLionsSolid.hpp b/src/coreComponents/constitutive/solid/DuvautLionsSolid.hpp index a3a857b82ee..21f05884c48 100644 --- a/src/coreComponents/constitutive/solid/DuvautLionsSolid.hpp +++ b/src/coreComponents/constitutive/solid/DuvautLionsSolid.hpp @@ -60,6 +60,7 @@ class DuvautLionsSolidUpdates : public UPDATE_BASE using UPDATE_BASE::smallStrainUpdate; //using UPDATE_BASE::smallStrainUpdate_ElasticOnly; + using UPDATE_BASE::smallStrainUpdate_thermal; using UPDATE_BASE::saveConvergedState; using UPDATE_BASE::viscousStateUpdate; @@ -124,6 +125,28 @@ class DuvautLionsSolidUpdates : public UPDATE_BASE this->smallStrainUpdate( k, q, timeIncrement, strainIncrement, stress, stiffness.m_c ); } + GEOS_HOST_DEVICE + virtual void smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( & stress )[6], + DiscretizationOps & stiffness, + real64 ( & dStress_dTemperature )[6] ) const override final + { + this->smallStrainUpdate( k, q, timeIncrement, strainIncrementNoThermalStrain, stress, stiffness ); + + real64 const ( &stiffnessTensor )[6][6] = stiffness.m_c; + + for( int i=0; i<6; ++i ) + { + for( int j=0; j<3; ++j ) + { + dStress_dTemperature[i] += stiffnessTensor[i][j]; + } + } + } + real64 const m_relaxationTime; }; diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp index a6663232fe8..9fe62b64e91 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropic.hpp @@ -125,6 +125,15 @@ class ElasticIsotropicUpdates : public SolidBaseUpdates real64 ( &stress )[6], DiscretizationOps & stiffness ) const; + GEOS_HOST_DEVICE + virtual void smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( &stress )[6], + DiscretizationOps & stiffness, + real64 ( &dStress_dTemperature )[6] ) const; + GEOS_HOST_DEVICE virtual void getElasticStiffness( localIndex const k, localIndex const q, @@ -363,6 +372,26 @@ void ElasticIsotropicUpdates::smallStrainUpdate( localIndex const k, stiffness.m_shearModulus = m_shearModulus[k]; } +GEOS_HOST_DEVICE +inline +void ElasticIsotropicUpdates::smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( & stress )[6], + DiscretizationOps & stiffness, + real64 ( & dStress_dTemperature )[6] ) const +{ + smallStrainUpdate( k, q, timeIncrement, strainIncrementNoThermalStrain, stress, stiffness ); + + dStress_dTemperature[0] = 3 * m_bulkModulus[k]; + dStress_dTemperature[1] = 3 * m_bulkModulus[k]; + dStress_dTemperature[2] = 3 * m_bulkModulus[k]; + dStress_dTemperature[3] = 0; + dStress_dTemperature[4] = 0; + dStress_dTemperature[5] = 0; +} + GEOS_HOST_DEVICE GEOS_FORCE_INLINE void ElasticIsotropicUpdates::viscousStateUpdate( localIndex const k, diff --git a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp index 8db58f1ed5b..b14f960ae3d 100644 --- a/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp +++ b/src/coreComponents/constitutive/solid/ElasticIsotropicPressureDependent.hpp @@ -105,6 +105,15 @@ class ElasticIsotropicPressureDependentUpdates : public SolidBaseUpdates real64 ( &stress )[6], DiscretizationOps & stiffness ) const; + GEOS_HOST_DEVICE + virtual void smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( &stress )[6], + DiscretizationOps & stiffness, + real64 ( &dStress_dTemperature )[6] ) const; + GEOS_HOST_DEVICE virtual void getElasticStiffness( localIndex const k, localIndex const q, @@ -467,6 +476,30 @@ void ElasticIsotropicPressureDependentUpdates::smallStrainUpdate( localIndex con stiffness.m_shearModulus = m_shearModulus[k]; } +GEOS_HOST_DEVICE +inline +void ElasticIsotropicPressureDependentUpdates::smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( & stress )[6], + DiscretizationOps & stiffness, + real64 ( & dStress_dTemperature )[6] ) const +{ + smallStrainUpdate( k, q, timeIncrement, strainIncrementNoThermalStrain, stress, stiffness ); + + real64 stiffnessTensor[6][6] = {{}}; + getElasticStiffness( k, q, stiffnessTensor ); + + for( int i=0; i<6; ++i ) + { + for( int j=0; j<3; ++j ) + { + dStress_dTemperature[i] += stiffnessTensor[i][j]; + } + } +} + GEOS_HOST_DEVICE GEOS_FORCE_INLINE void ElasticIsotropicPressureDependentUpdates::viscousStateUpdate( localIndex const k, diff --git a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp index e3f9075a73b..e1948b60558 100644 --- a/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticOrthotropic.hpp @@ -147,6 +147,15 @@ class ElasticOrthotropicUpdates : public SolidBaseUpdates real64 ( &stress )[6], DiscretizationOps & stiffness ) const final; + GEOS_HOST_DEVICE + virtual void smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( &stress )[6], + DiscretizationOps & stiffness, + real64 ( &dStress_dTemperature )[6] ) const final; + GEOS_HOST_DEVICE virtual void getElasticStrain( localIndex const k, localIndex const q, @@ -396,6 +405,30 @@ void ElasticOrthotropicUpdates::smallStrainUpdate( localIndex const k, stiffness.m_c66 = m_c66[k]; } +GEOS_HOST_DEVICE +inline +void ElasticOrthotropicUpdates::smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( & stress )[6], + DiscretizationOps & stiffness, + real64 ( & dStress_dTemperature )[6] ) const +{ + smallStrainUpdate( k, q, timeIncrement, strainIncrementNoThermalStrain, stress, stiffness ); + + real64 stiffnessTensor[6][6] = {{}}; + getElasticStiffness( k, q, stiffnessTensor ); + + for( int i=0; i<6; ++i ) + { + for( int j=0; j<3; ++j ) + { + dStress_dTemperature[i] += stiffnessTensor[i][j]; + } + } +} + /** * @class ElasticOrthotropic * diff --git a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp index b402a3fe816..d5630e7d45f 100644 --- a/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/ElasticTransverseIsotropic.hpp @@ -139,6 +139,15 @@ class ElasticTransverseIsotropicUpdates : public SolidBaseUpdates real64 ( &stress )[6], DiscretizationOps & stiffness ) const final; + GEOS_HOST_DEVICE + virtual void smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( &stress )[6], + DiscretizationOps & stiffness, + real64 ( &dStress_dTemperature )[6] ) const final; + // miscellaneous getters GEOS_HOST_DEVICE @@ -360,6 +369,30 @@ void ElasticTransverseIsotropicUpdates::smallStrainUpdate( localIndex const k, stiffness.m_c66 = m_c66[k]; } +GEOS_HOST_DEVICE +inline +void ElasticTransverseIsotropicUpdates::smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( & stress )[6], + DiscretizationOps & stiffness, + real64 ( & dStress_dTemperature )[6] ) const +{ + smallStrainUpdate( k, q, timeIncrement, strainIncrementNoThermalStrain, stress, stiffness ); + + real64 stiffnessTensor[6][6] = {{}}; + getElasticStiffness( k, q, stiffnessTensor ); + + for( int i=0; i<6; ++i ) + { + for( int j=0; j<3; ++j ) + { + dStress_dTemperature[i] += stiffnessTensor[i][j]; + } + } +} + /** * @class ElasticTransverseIsotropic * diff --git a/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp b/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp index 93b0d7da7e7..d8aefdbff7e 100644 --- a/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp +++ b/src/coreComponents/constitutive/solid/ModifiedCamClay.hpp @@ -98,6 +98,7 @@ class ModifiedCamClayUpdates : public ElasticIsotropicPressureDependentUpdates // Bring in base implementations to prevent hiding warnings using ElasticIsotropicPressureDependentUpdates::smallStrainUpdate; + using ElasticIsotropicPressureDependentUpdates::smallStrainUpdate_thermal; GEOS_HOST_DEVICE void evaluateYield( real64 const p, @@ -141,6 +142,15 @@ class ModifiedCamClayUpdates : public ElasticIsotropicPressureDependentUpdates real64 ( &stress )[6], real64 ( &stiffness )[6][6] ) const override; + GEOS_HOST_DEVICE + virtual void smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( &stress )[6], + DiscretizationOps & stiffness, + real64 ( &dStress_dTemperature )[6] ) const; + GEOS_HOST_DEVICE virtual real64 getBulkModulus( localIndex const k ) const override final { @@ -468,6 +478,29 @@ void ModifiedCamClayUpdates::smallStrainUpdate( localIndex const k, smallStrainUpdate( k, q, timeIncrement, strainIncrement, stress, stiffness.m_c ); } +GEOS_HOST_DEVICE +inline +void ModifiedCamClayUpdates::smallStrainUpdate_thermal( localIndex const k, + localIndex const q, + real64 const & timeIncrement, + real64 const ( &strainIncrementNoThermalStrain )[6], + real64 ( & stress )[6], + DiscretizationOps & stiffness, + real64 ( & dStress_dTemperature )[6] ) const +{ + smallStrainUpdate( k, q, timeIncrement, strainIncrementNoThermalStrain, stress, stiffness ); + + real64 const ( &stiffnessTensor )[6][6] = stiffness.m_c; + + for( int i=0; i<6; ++i ) + { + for( int j=0; j<3; ++j ) + { + dStress_dTemperature[i] += stiffnessTensor[i][j]; + } + } +} + /** * @class ModifiedCamClay * diff --git a/src/coreComponents/constitutive/solid/PorousDamageSolid.hpp b/src/coreComponents/constitutive/solid/PorousDamageSolid.hpp index d32923f83de..077123084fd 100644 --- a/src/coreComponents/constitutive/solid/PorousDamageSolid.hpp +++ b/src/coreComponents/constitutive/solid/PorousDamageSolid.hpp @@ -59,7 +59,7 @@ class PorousDamageSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPor real64 const & timeIncrement, real64 const & pressure_n, real64 const & pressure, - real64 const & temperature, + real64 const & deltaTemperature, real64 const & deltaTemperatureFromLastStep, real64 const ( &strainIncrement )[6], real64 ( & totalStress )[6], @@ -80,7 +80,8 @@ class PorousDamageSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPor timeIncrement, pressure_n, pressure, - temperature, + deltaTemperature, + deltaTemperatureFromLastStep, strainIncrement, totalStress, dTotalStress_dPressure, @@ -266,33 +267,44 @@ class PorousDamageSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPor real64 const & timeIncrement, real64 const & pressure_n, real64 const & pressure, - real64 const & temperature, + real64 const & deltaTemperature, + real64 const & deltaTemperatureFromLastStep, real64 const ( &strainIncrement )[6], real64 ( & totalStress )[6], real64 ( & dTotalStress_dPressure )[6], real64 ( & dTotalStress_dTemperature )[6], DiscretizationOps & stiffness ) const { + GEOS_UNUSED_VAR( deltaTemperature ); + updateBiotCoefficientAndAssignModuli( k ); + real64 const thermalExpansionCoefficient = m_solidUpdate.getThermalExpansionCoefficient( k ); + + real64 strainIncrementNoThermalStrain[6]{}; + strainIncrementNoThermalStrain[0] = strainIncrement[0] - thermalExpansionCoefficient * deltaTemperatureFromLastStep; + strainIncrementNoThermalStrain[1] = strainIncrement[1] - thermalExpansionCoefficient * deltaTemperatureFromLastStep; + strainIncrementNoThermalStrain[2] = strainIncrement[2] - thermalExpansionCoefficient * deltaTemperatureFromLastStep; + strainIncrementNoThermalStrain[3] = strainIncrement[3]; + strainIncrementNoThermalStrain[4] = strainIncrement[4]; + strainIncrementNoThermalStrain[5] = strainIncrement[5]; + // Compute total stress increment and its derivative w.r.t. pressure - m_solidUpdate.smallStrainUpdate( k, - q, - timeIncrement, - strainIncrement, - totalStress, // first effective stress increment accumulated - stiffness ); + m_solidUpdate.smallStrainUpdate_thermal( k, + q, + timeIncrement, + strainIncrementNoThermalStrain, + totalStress, // first effective stress increment accumulated + stiffness, + dTotalStress_dTemperature ); // Add the contributions of pressure and temperature to the total stress real64 const biotCoefficient = m_porosityUpdate.getBiotCoefficient( k ); - real64 const thermalExpansionCoefficient = m_solidUpdate.getThermalExpansionCoefficient( k ); - real64 const bulkModulus = m_solidUpdate.getBulkModulus( k ); - real64 const thermalExpansionCoefficientTimesBulkModulus = thermalExpansionCoefficient * bulkModulus; real64 const pressureDamage = m_solidUpdate.pressureDamageFunction( k, q ); real64 const damagedBiotCoefficient = pressureDamage * biotCoefficient; - LvArray::tensorOps::symAddIdentity< 3 >( totalStress, pressureDamage * pressure_n - damagedBiotCoefficient * pressure - 3 * thermalExpansionCoefficientTimesBulkModulus * temperature ); + LvArray::tensorOps::symAddIdentity< 3 >( totalStress, pressureDamage * pressure_n - damagedBiotCoefficient * pressure ); dTotalStress_dPressure[0] = -damagedBiotCoefficient; dTotalStress_dPressure[1] = -damagedBiotCoefficient; @@ -301,12 +313,7 @@ class PorousDamageSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPor dTotalStress_dPressure[4] = 0; dTotalStress_dPressure[5] = 0; - dTotalStress_dTemperature[0] = -3 * thermalExpansionCoefficientTimesBulkModulus; - dTotalStress_dTemperature[1] = -3 * thermalExpansionCoefficientTimesBulkModulus; - dTotalStress_dTemperature[2] = -3 * thermalExpansionCoefficientTimesBulkModulus; - dTotalStress_dTemperature[3] = 0; - dTotalStress_dTemperature[4] = 0; - dTotalStress_dTemperature[5] = 0; + LvArray::tensorOps::scale< 6 >( dTotalStress_dTemperature, -thermalExpansionCoefficient ); } }; diff --git a/src/coreComponents/constitutive/solid/PorousSolid.hpp b/src/coreComponents/constitutive/solid/PorousSolid.hpp index 1edd995d78b..c956b1fa436 100644 --- a/src/coreComponents/constitutive/solid/PorousSolid.hpp +++ b/src/coreComponents/constitutive/solid/PorousSolid.hpp @@ -98,6 +98,7 @@ class PorousSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPorosity, timeIncrement, pressure, deltaTemperature, + deltaTemperatureFromLastStep, strainIncrement, totalStress, dTotalStress_dPressure, @@ -150,11 +151,13 @@ class PorousSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPorosity, // Compute total stress increment and its derivative real64 const deltaTemperature = temperature - referenceTemperature; + real64 const deltaTemperatureFromLastStep = temperature - temperature_n; computeTotalStress( k, q, timeIncrement, pressure, deltaTemperature, + deltaTemperatureFromLastStep, strainIncrement, totalStress, dTotalStress_dPressure, // To pass something here @@ -275,29 +278,40 @@ class PorousSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPorosity, real64 const & timeIncrement, real64 const & pressure, real64 const & deltaTemperature, + real64 const & deltaTemperatureFromLastStep, real64 const ( &strainIncrement )[6], real64 ( & totalStress )[6], real64 ( & dTotalStress_dPressure )[6], real64 ( & dTotalStress_dTemperature )[6], DiscretizationOps & stiffness ) const { + GEOS_UNUSED_VAR( deltaTemperature ); + updateBiotCoefficientAndAssignModuli( k ); - // Compute total stress increment and its derivative w.r.t. pressure - m_solidUpdate.smallStrainUpdate( k, - q, - timeIncrement, - strainIncrement, - totalStress, // first effective stress increment accumulated - stiffness ); + real64 const thermalExpansionCoefficient = m_solidUpdate.getThermalExpansionCoefficient( k ); + + real64 strainIncrementNoThermalStrain[6]{}; + strainIncrementNoThermalStrain[0] = strainIncrement[0] - thermalExpansionCoefficient * deltaTemperatureFromLastStep; + strainIncrementNoThermalStrain[1] = strainIncrement[1] - thermalExpansionCoefficient * deltaTemperatureFromLastStep; + strainIncrementNoThermalStrain[2] = strainIncrement[2] - thermalExpansionCoefficient * deltaTemperatureFromLastStep; + strainIncrementNoThermalStrain[3] = strainIncrement[3]; + strainIncrementNoThermalStrain[4] = strainIncrement[4]; + strainIncrementNoThermalStrain[5] = strainIncrement[5]; + + // Compute total stress increment and its derivative w.r.t. pressure and temperature + m_solidUpdate.smallStrainUpdate_thermal( k, + q, + timeIncrement, + strainIncrementNoThermalStrain, + totalStress, // first effective stress increment accumulated + stiffness, + dTotalStress_dTemperature ); // Add the contributions of pressure and temperature to the total stress real64 const biotCoefficient = m_porosityUpdate.getBiotCoefficient( k ); - real64 const thermalExpansionCoefficient = m_solidUpdate.getThermalExpansionCoefficient( k ); - real64 const bulkModulus = m_solidUpdate.getBulkModulus( k ); - real64 const thermalExpansionCoefficientTimesBulkModulus = thermalExpansionCoefficient * bulkModulus; - LvArray::tensorOps::symAddIdentity< 3 >( totalStress, -biotCoefficient * pressure - 3 * thermalExpansionCoefficientTimesBulkModulus * deltaTemperature ); + LvArray::tensorOps::symAddIdentity< 3 >( totalStress, -biotCoefficient * pressure ); // Compute derivatives of total stress dTotalStress_dPressure[0] = -biotCoefficient; @@ -307,12 +321,7 @@ class PorousSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPorosity, dTotalStress_dPressure[4] = 0; dTotalStress_dPressure[5] = 0; - dTotalStress_dTemperature[0] = -3 * thermalExpansionCoefficientTimesBulkModulus; - dTotalStress_dTemperature[1] = -3 * thermalExpansionCoefficientTimesBulkModulus; - dTotalStress_dTemperature[2] = -3 * thermalExpansionCoefficientTimesBulkModulus; - dTotalStress_dTemperature[3] = 0; - dTotalStress_dTemperature[4] = 0; - dTotalStress_dTemperature[5] = 0; + LvArray::tensorOps::scale< 6 >( dTotalStress_dTemperature, -thermalExpansionCoefficient ); }