diff --git a/inputFiles/thermoPoromechanics/ThermoDruckerPrager_1DCooling_fim_smoke.xml b/inputFiles/thermoPoromechanics/ThermoDruckerPrager_1DCooling_fim_smoke.xml new file mode 100644 index 00000000000..fa6d45549b0 --- /dev/null +++ b/inputFiles/thermoPoromechanics/ThermoDruckerPrager_1DCooling_fim_smoke.xml @@ -0,0 +1,159 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/thermoPoromechanics/ThermoElastic_1DCooling_fim_smoke.xml b/inputFiles/thermoPoromechanics/ThermoElastic_1DCooling_fim_smoke.xml new file mode 100644 index 00000000000..d03c8003a39 --- /dev/null +++ b/inputFiles/thermoPoromechanics/ThermoElastic_1DCooling_fim_smoke.xml @@ -0,0 +1,155 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/thermoPoromechanics/ThermoElastic_1DCooling_sq_smoke.xml b/inputFiles/thermoPoromechanics/ThermoElastic_1DCooling_sq_smoke.xml new file mode 100644 index 00000000000..6c56d1f51cf --- /dev/null +++ b/inputFiles/thermoPoromechanics/ThermoElastic_1DCooling_sq_smoke.xml @@ -0,0 +1,163 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/thermoPoromechanics/ThermoMech_1DCooling_base.xml b/inputFiles/thermoPoromechanics/ThermoMech_1DCooling_base.xml new file mode 100644 index 00000000000..858d025e01b --- /dev/null +++ b/inputFiles/thermoPoromechanics/ThermoMech_1DCooling_base.xml @@ -0,0 +1,116 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/thermoPoromechanics/thermoPoromechanics.ats b/inputFiles/thermoPoromechanics/thermoPoromechanics.ats index 443edf0817f..d59ce6bc506 100644 --- a/inputFiles/thermoPoromechanics/thermoPoromechanics.ats +++ b/inputFiles/thermoPoromechanics/thermoPoromechanics.ats @@ -64,6 +64,30 @@ decks = [ partitions=((1, 1, 1), (1, 2, 1)), restart_step=0, check_step=1, + restartcheck_params=RestartcheckParameters(**restartcheck_params)), + TestDeck( + name="ThermoDruckerPrager_1DCooling_fim_smoke", + description= + '1D thermomechanics fim 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_fim_smoke", + description= + '1D thermomechanics fim 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)), + TestDeck( + name="ThermoElastic_1DCooling_sq_smoke", + description= + '1D thermomechanics sequential 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)) ] diff --git a/src/coreComponents/constitutive/solid/Damage.hpp b/src/coreComponents/constitutive/solid/Damage.hpp index b54bb30f372..84ed4e28024 100644 --- a/src/coreComponents/constitutive/solid/Damage.hpp +++ b/src/coreComponents/constitutive/solid/Damage.hpp @@ -293,7 +293,6 @@ class DamageUpdates : public UPDATE_BASE stiffness.scaleParams( factor ); } - // 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..689d2289d92 100644 --- a/src/coreComponents/constitutive/solid/DamageSpectral.hpp +++ b/src/coreComponents/constitutive/solid/DamageSpectral.hpp @@ -268,7 +268,6 @@ class DamageSpectralUpdates : public DamageUpdates< UPDATE_BASE > smallStrainUpdate( k, q, timeIncrement, strainIncrement, stress, stiffness.m_c ); } - GEOS_HOST_DEVICE virtual real64 getStrainEnergyDensity( localIndex const k, localIndex const q ) const override final diff --git a/src/coreComponents/constitutive/solid/DamageVolDev.hpp b/src/coreComponents/constitutive/solid/DamageVolDev.hpp index 951eded451f..cad75993cbd 100644 --- a/src/coreComponents/constitutive/solid/DamageVolDev.hpp +++ b/src/coreComponents/constitutive/solid/DamageVolDev.hpp @@ -153,7 +153,6 @@ class DamageVolDevUpdates : public DamageUpdates< UPDATE_BASE > } } - GEOS_HOST_DEVICE virtual real64 getStrainEnergyDensity( localIndex const k, localIndex const q ) const override final diff --git a/src/coreComponents/constitutive/solid/DelftEgg.hpp b/src/coreComponents/constitutive/solid/DelftEgg.hpp index 167301468f0..5d9afab89f4 100644 --- a/src/coreComponents/constitutive/solid/DelftEgg.hpp +++ b/src/coreComponents/constitutive/solid/DelftEgg.hpp @@ -449,8 +449,6 @@ void DelftEggUpdates::smallStrainUpdate( localIndex const k, smallStrainUpdate( k, q, timeIncrement, strainIncrement, stress, stiffness.m_c ); } - - /** * @class DelftEgg * diff --git a/src/coreComponents/constitutive/solid/DruckerPrager.hpp b/src/coreComponents/constitutive/solid/DruckerPrager.hpp index e59494d7b06..02e2a44fcf7 100644 --- a/src/coreComponents/constitutive/solid/DruckerPrager.hpp +++ b/src/coreComponents/constitutive/solid/DruckerPrager.hpp @@ -339,8 +339,6 @@ void DruckerPragerUpdates::smallStrainUpdate( localIndex const k, smallStrainUpdate( k, q, timeIncrement, strainIncrement, stress, stiffness.m_c ); } - - /** * @class DruckerPrager * diff --git a/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp b/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp index dc79394677f..ebe6e8d14b9 100644 --- a/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp +++ b/src/coreComponents/constitutive/solid/DruckerPragerExtended.hpp @@ -116,7 +116,6 @@ class DruckerPragerExtendedUpdates : public ElasticIsotropicUpdates real64 ( &stress )[6], real64 ( &stiffness )[6][6] ) const override; - GEOS_HOST_DEVICE inline virtual void saveConvergedState( localIndex const k, @@ -365,8 +364,6 @@ void DruckerPragerExtendedUpdates::smallStrainUpdate( localIndex const k, smallStrainUpdate( k, q, timeIncrement, strainIncrement, stress, stiffness.m_c ); } - - /** * @class DruckerPragerExtended * diff --git a/src/coreComponents/constitutive/solid/PorousDamageSolid.hpp b/src/coreComponents/constitutive/solid/PorousDamageSolid.hpp index d32923f83de..99e85bc4223 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,45 @@ 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 + strainIncrementNoThermalStrain, + totalStress, stiffness ); + stiffness.computeTemperatureDerivative( thermalExpansionCoefficient, 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; @@ -300,13 +313,6 @@ class PorousDamageSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPor dTotalStress_dPressure[3] = 0; 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; } }; diff --git a/src/coreComponents/constitutive/solid/PorousSolid.hpp b/src/coreComponents/constitutive/solid/PorousSolid.hpp index d7f27433787..b586181b420 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,41 @@ 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 + 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( k, q, timeIncrement, - strainIncrement, - totalStress, // first effective stress increment accumulated + strainIncrementNoThermalStrain, + totalStress, stiffness ); + stiffness.computeTemperatureDerivative( thermalExpansionCoefficient, 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,13 +322,6 @@ 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; - } }; diff --git a/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsFullyAnisotropic.hpp b/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsFullyAnisotropic.hpp index 19a7db1bea3..101dd175048 100644 --- a/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsFullyAnisotropic.hpp +++ b/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsFullyAnisotropic.hpp @@ -65,6 +65,15 @@ struct SolidModelDiscretizationOpsFullyAnisotropic : public SolidModelDiscretiza LvArray::tensorOps::scale< 6, 6 >( m_c, scale ); } + GEOS_HOST_DEVICE + inline + void computeTemperatureDerivative( real64 const thermalExpansionCoefficient, real64 ( & dStress_dT )[6] ) const + { + for( int i = 0; i < 6; ++i ) + { + dStress_dT[i] = -thermalExpansionCoefficient * ( m_c[i][0] + m_c[i][1] + m_c[i][2] ); + } + } real64 m_c[6][6]; }; diff --git a/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsIsotropic.hpp b/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsIsotropic.hpp index 78172ee75f9..68a4922e1d7 100644 --- a/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsIsotropic.hpp @@ -103,6 +103,19 @@ struct SolidModelDiscretizationOpsIsotropic : public SolidModelDiscretizationOps m_shearModulus *= scale; } + GEOS_HOST_DEVICE + inline + void computeTemperatureDerivative( real64 const thermalExpansionCoefficient, real64 ( & dStress_dT )[6] ) const + { + real64 const val = -3.0 * thermalExpansionCoefficient * m_bulkModulus; + dStress_dT[0] = val; + dStress_dT[1] = val; + dStress_dT[2] = val; + dStress_dT[3] = 0; + dStress_dT[4] = 0; + dStress_dT[5] = 0; + } + real64 m_bulkModulus; ///< Bulk modulus real64 m_shearModulus; ///< Shear modulus }; diff --git a/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsOrthotropic.hpp b/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsOrthotropic.hpp index a16fdb8fb23..859fc8fc3e9 100644 --- a/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsOrthotropic.hpp +++ b/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsOrthotropic.hpp @@ -72,6 +72,18 @@ struct SolidModelDiscretizationOpsOrthotropic : public SolidModelDiscretizationO m_c66 *= scale; } + GEOS_HOST_DEVICE + inline + void computeTemperatureDerivative( real64 const thermalExpansionCoefficient, real64 ( & dStress_dT )[6] ) const + { + dStress_dT[0] = -thermalExpansionCoefficient * ( m_c11 + m_c12 + m_c13 ); + dStress_dT[1] = -thermalExpansionCoefficient * ( m_c12 + m_c22 + m_c23 ); + dStress_dT[2] = -thermalExpansionCoefficient * ( m_c13 + m_c23 + m_c33 ); + dStress_dT[3] = 0; + dStress_dT[4] = 0; + dStress_dT[5] = 0; + } + real64 m_c11; real64 m_c12; real64 m_c13; diff --git a/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsTransverseIsotropic.hpp b/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsTransverseIsotropic.hpp index 9cf78b133de..bfcb13b2202 100644 --- a/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsTransverseIsotropic.hpp +++ b/src/coreComponents/constitutive/solid/SolidModelDiscretizationOpsTransverseIsotropic.hpp @@ -74,6 +74,20 @@ struct SolidModelDiscretizationOpsTransverseIsotropic : public SolidModelDiscret m_c66 *= scale; } + // dStress_dT[i] = -alpha * sum_{j=0}^{2} C[i][j], where C12 = C11 - 2*C66 + // Shear rows are zero because TI symmetry decouples normal and shear components. + GEOS_HOST_DEVICE + inline + void computeTemperatureDerivative( real64 const thermalExpansionCoefficient, real64 ( & dStress_dT )[6] ) const + { + dStress_dT[0] = -thermalExpansionCoefficient * ( 2 * m_c11 - 2 * m_c66 + m_c13 ); + dStress_dT[1] = dStress_dT[0]; + dStress_dT[2] = -thermalExpansionCoefficient * ( 2 * m_c13 + m_c33 ); + dStress_dT[3] = 0; + dStress_dT[4] = 0; + dStress_dT[5] = 0; + } + real64 m_c11; ///< (1,1) element of TI stiffness matrix real64 m_c13; ///< (1,3) element of TI stiffness matrix real64 m_c33; ///< (3,3) element of TI stiffness matrix diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp index 9c919cffa96..4432bd4356a 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp @@ -49,6 +49,7 @@ #include "physicsSolvers/solidMechanics/kernels/SolidMechanicsKernelsDispatchTypeList.hpp" #include "physicsSolvers/solidMechanics/kernels/SolidMechanicsFixedStressThermoPoromechanicsKernelsDispatchTypeList.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBase.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" namespace geos { @@ -1000,6 +1001,16 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS solidMechanics::arrayView2dLayoutStrain avgPlasticStrain = subRegion.getField< solidMechanics::averagePlasticStrain >(); solidMechanics::arrayView2dLayoutAvgStress avgStress = subRegion.getField< solidMechanics::averageStress >(); + arrayView1d< real64 const > const temperature = + subRegion.hasField< fields::flow::temperature >() + ? subRegion.getField< fields::flow::temperature >().toViewConst() + : arrayView1d< real64 const >{}; + + arrayView1d< real64 const > const temperature_n = + subRegion.hasField< fields::flow::temperature_n >() + ? subRegion.getField< fields::flow::temperature_n >().toViewConst() + : arrayView1d< real64 const >{}; + constitutive::ConstitutivePassThru< SolidBase >::execute( constitutiveRelation, [&] ( auto & solidModel ) { using SOLID_TYPE = TYPEOFREF( solidModel ); @@ -1019,7 +1030,9 @@ void SolidMechanicsLagrangianFEM::implicitStepComplete( real64 const & GEOS_UNUS avgStrain, avgPlasticStrain, stress, - avgStress ); + avgStress, + temperature, + temperature_n ); } ); diff --git a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StressStrainAverageKernels.hpp b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StressStrainAverageKernels.hpp index 9765e7a262d..216ee15ee7f 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/kernels/StressStrainAverageKernels.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/kernels/StressStrainAverageKernels.hpp @@ -28,6 +28,7 @@ #include "mesh/CellElementSubRegion.hpp" #include "mesh/utilities/AverageOverQuadraturePointsKernel.hpp" #include "physicsSolvers/solidMechanics/SolidMechanicsFields.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" namespace geos { @@ -66,6 +67,8 @@ class AverageStressStrainOverQuadraturePoints : * @param avgStrain the strain averaged over quadrature points * @param stress the stress solution field * @param avgStress the stress averaged over quadrature points + * @param temperature the temperature field (may be empty for non-thermal simulations) + * @param temperature_n the temperature at the previous time step (may be empty for non-thermal simulations) */ AverageStressStrainOverQuadraturePoints( NodeManager & nodeManager, EdgeManager const & edgeManager, @@ -78,7 +81,9 @@ class AverageStressStrainOverQuadraturePoints : fields::solidMechanics::arrayView2dLayoutStrain const avgStrain, fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain, arrayView3d< real64 const, solid::STRESS_USD > const stress, - fields::solidMechanics::arrayView2dLayoutAvgStress const avgStress ): + fields::solidMechanics::arrayView2dLayoutAvgStress const avgStress, + arrayView1d< real64 const > const temperature, + arrayView1d< real64 const > const temperature_n ): Base( nodeManager, edgeManager, faceManager, @@ -90,7 +95,9 @@ class AverageStressStrainOverQuadraturePoints : m_avgStrain( avgStrain ), m_avgPlasticStrain( avgPlasticStrain ), m_stress( stress ), - m_avgStress( avgStress ) + m_avgStress( avgStress ), + m_temperature( temperature ), + m_temperature_n( temperature_n ) {} /** @@ -151,6 +158,11 @@ class AverageStressStrainOverQuadraturePoints : real64 elasticStrainInc[6] = {0.0}; m_solidUpdate.getElasticStrainInc( k, q, elasticStrainInc ); + real64 const thermalExpansionCoefficient = m_solidUpdate.getThermalExpansionCoefficient( k ); + real64 const deltaTemperature = ( m_temperature.size() > 0 ) + ? ( m_temperature[k] - m_temperature_n[k] ) + : 0.0; + real64 conversionFactor[6] = {1.0, 1.0, 1.0, 0.5, 0.5, 0.5}; // used for converting from engineering shear to tensor shear for( int icomp = 0; icomp < 6; ++icomp ) @@ -158,12 +170,17 @@ class AverageStressStrainOverQuadraturePoints : m_avgStrain[k][icomp] += conversionFactor[icomp]*detJxW*strain[icomp]/m_elementVolume[k]; m_avgStress[k][icomp] += detJxW*m_stress[k][q][icomp]/m_elementVolume[k]; + // Thermal strain is purely volumetric: subtract thermal strain from normal components so that + // only the mechanical (plastic) part is accumulated. + real64 const thermalStrainInc = ( icomp < 3 ) ? thermalExpansionCoefficient * deltaTemperature : 0.0; + real64 const mechanicalStrainInc = strainInc[icomp] - thermalStrainInc; + // This is a hack to handle boundary conditions such as those seen in plane-strain wellbore problems // Essentially, if bcs are constraining the strain (and thus total displacement), we do not accumulate any plastic strain (regardless // of stresses in material law) - if( std::abs( strainInc[icomp] ) > 1.0e-8 ) + if( std::abs( mechanicalStrainInc ) > 1.0e-8 ) { - m_avgPlasticStrain[k][icomp] += conversionFactor[icomp]*detJxW*(strainInc[icomp] - elasticStrainInc[icomp])/m_elementVolume[k]; + m_avgPlasticStrain[k][icomp] += conversionFactor[icomp]*detJxW*(mechanicalStrainInc - elasticStrainInc[icomp])/m_elementVolume[k]; } } } @@ -217,6 +234,12 @@ class AverageStressStrainOverQuadraturePoints : /// The average stress fields::solidMechanics::arrayView2dLayoutAvgStress const m_avgStress; + /// The temperature at the current time step (empty for non-thermal simulations) + arrayView1d< real64 const > const m_temperature; + + /// The temperature at the previous time step (empty for non-thermal simulations) + arrayView1d< real64 const > const m_temperature_n; + }; @@ -258,11 +281,14 @@ class AverageStressStrainOverQuadraturePointsKernelFactory fields::solidMechanics::arrayView2dLayoutStrain const avgStrain, fields::solidMechanics::arrayView2dLayoutStrain const avgPlasticStrain, arrayView3d< real64 const, solid::STRESS_USD > const stress, - fields::solidMechanics::arrayView2dLayoutAvgStress const avgStress ) + fields::solidMechanics::arrayView2dLayoutAvgStress const avgStress, + arrayView1d< real64 const > const temperature = {}, + arrayView1d< real64 const > const temperature_n = {} ) { AverageStressStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE > kernel( nodeManager, edgeManager, faceManager, elementSubRegion, finiteElementSpace, - solidModel, displacement, displacementInc, avgStrain, avgPlasticStrain, stress, avgStress ); + solidModel, displacement, displacementInc, avgStrain, avgPlasticStrain, stress, avgStress, + temperature, temperature_n ); AverageStressStrainOverQuadraturePoints< FE_TYPE, SOLID_TYPE >::template kernelLaunch< POLICY >( elementSubRegion.size(), kernel );