Skip to content

To account for user-defined inelastic strain in solid mechanics #3845

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
frankfeifan wants to merge 5 commits into
base: develop
Choose a base branch
from
Draft

To account for user-defined inelastic strain in solid mechanics #3845

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
0d75919
initial commit to enable anelastic strain
frankfeifan Oct 3, 2025
4d36d9a
make anelastic strain increasing with time
frankfeifan Oct 7, 2025
32c75d5
Merge remote-tracking branch 'origin/develop' into feature/frankfei/e…
frankfeifan Aug 3, 2026
dcfae80
changed the input anelastic strain to a 3D vector, changed increment …
frankfeifan Aug 4, 2026
9eea551
added a xml
frankfeifan Aug 4, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
214 changes: 214 additions & 0 deletions inputFiles/phaseField/PhaseFieldFracture_AnelasticStrain.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,214 @@
<?xml version="1.0" ?>

<Problem>
<Solvers
gravityVector="{ 0.0, 0.0, 0.0 }">
<PhaseFieldFracture
name="PhaseFieldSolve"
solidSolverName="SolidMechSolve"
damageSolverName="DamageSolve"
logLevel="1"
targetRegions="{ DamageRegion, ShrinkRegion }">
<NonlinearSolverParameters
couplingType="Sequential"
lineSearchAction="None"/>
</PhaseFieldFracture>

<SolidMechanicsLagrangianFEM
name="SolidMechSolve"
logLevel="1"
discretization="FE1"
targetRegions="{ DamageRegion, ShrinkRegion }"
>
<NonlinearSolverParameters
newtonMaxIter="20"
newtonTol="1.0e-6"/>
<LinearSolverParameters
solverType="gmres"
preconditionerType="amg"/>
</SolidMechanicsLagrangianFEM>

<PhaseFieldDamageFEM
name="DamageSolve"
discretization="FE1"
timeIntegrationOption="SteadyState"
fieldName="Damage"
localDissipation="Quadratic"
irreversibilityFlag="1"
damageUpperBound="0.98"
logLevel="1"
targetRegions="{ DamageRegion, ShrinkRegion }">
<NonlinearSolverParameters
newtonMaxIter="20"
newtonTol="1.0e-6"/>
<LinearSolverParameters
solverType="gmres"
preconditionerType="amg"/>
</PhaseFieldDamageFEM>

<SurfaceGenerator
name="SurfaceGen"
targetRegions="{ DamageRegion }"
initialRockToughness="0"
mpiCommOrder="1"/>
</Solvers>

<Events
maxTime="100.0">

<SoloEvent
name="preFracture"
target="/Solvers/SurfaceGen"/>

<PeriodicEvent
name="vtkOutput"
cycleFrequency="5"
target="/Outputs/VTKOutput"/>

<PeriodicEvent
name="solverApplications"
maxEventDt="0.1"
target="/Solvers/PhaseFieldSolve"/>

<PeriodicEvent
name="restarts"
cycleFrequency="50"
target="/Outputs/restartOutput"/>
</Events>

<Mesh>
<InternalMesh
name="mesh1"
elementTypes="{ C3D8 }"
xCoords="{ 0, 30 }"
yCoords="{ -5, -2.0, 0, 2.0, 5 }"
zCoords="{ 0, 0.1 }"
nx="{ 240 }"
ny="{ 20, 80, 80, 20 }"
nz="{ 1 }"
yBias="{0.75, 0, 0, -0.75}"
cellBlockNames="{ bottom, middle1, middle2, top }"/>
</Mesh>

<NumericalMethods>
<FiniteElements>
<FiniteElementSpace
name="FE1"
order="1"/>
</FiniteElements>
</NumericalMethods>

<ElementRegions>
<CellElementRegion
name="DamageRegion"
cellBlocks="{ middle1, middle2 }"
materialList="{ damageRock }"/>

<CellElementRegion
name="ShrinkRegion"
cellBlocks="{ bottom, top }"
materialList="{ rock }"/>

<SurfaceElementRegion
name="Fracture"
defaultAperture="1e-6"
materialList="{ damageRock }"/>
</ElementRegions>

<Constitutive>
<DamageElasticIsotropic
name="damageRock"
defaultDensity="2700"
defaultBulkModulus="1.73333333e4"
defaultShearModulus="8e3"
lengthScale="0.05"
defaultCriticalFractureEnergy="1.0"
criticalStrainEnergy="0"/>

<DamageElasticIsotropic
name="rock"
defaultDensity="2700"
defaultBulkModulus="1.73333333e4"
defaultShearModulus="8e3"
enableAnelasticStrain="1"
defaultAnelasticStrainRate="{ 0.0, -1e-4, 0.0 }"
lengthScale="0.05"
defaultCriticalFractureEnergy="1e20"
criticalStrainEnergy="0"/>
</Constitutive>

<FieldSpecifications>
<FieldSpecification
name="Gc"
initialCondition="1"
setNames="{ all }"
objectPath="ElementRegions/DamageRegion"
fieldName="damageRock_criticalFractureEnergy"
scale="1.0"/>

<FieldSpecification
name="GcShrinkage"
initialCondition="1"
setNames="{ all }"
objectPath="ElementRegions/ShrinkRegion"
fieldName="rock_criticalFractureEnergy"
scale="1e20"/>

<FieldSpecification
name="fixedY"
fieldName="totalDisplacement"
component="1"
objectPath="nodeManager"
scale="0.0"
setNames="{ yneg, ypos }"/>

<FieldSpecification
name="fixedX"
fieldName="totalDisplacement"
component="0"
objectPath="nodeManager"
scale="0.0"
setNames="{ xpos }"/>

<FieldSpecification
name="fixedZ"
fieldName="totalDisplacement"
component="2"
objectPath="nodeManager"
scale="0.0"
setNames="{ all }"/>

<FieldSpecification
name="separableFace"
initialCondition="1"
setNames="{ fracture }"
objectPath="faceManager"
fieldName="isFaceSeparable"
scale="1"/>

<FieldSpecification
name="frac"
initialCondition="1"
setNames="{ fracture }"
objectPath="faceManager"
fieldName="ruptureState"
scale="1"/>
</FieldSpecifications>

<Outputs>
<VTK
name="VTKOutput"
plotFileRoot="damageAnelasticStrain"/>

<Restart
name="restartOutput"/>
</Outputs>

<Geometry>
<Box
name="fracture"
xMin="{ -0.001, -0.001, -0.001 }"
xMax="{ 5.001, 0.001, 1.001 }"/>
</Geometry>

</Problem>
24 changes: 20 additions & 4 deletions src/coreComponents/constitutive/solid/CeramicDamage.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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 ),
Expand Down Expand Up @@ -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 );
}

/**
Expand All @@ -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 );
}


Expand Down
23 changes: 19 additions & 4 deletions src/coreComponents/constitutive/solid/DelftEgg.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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 ),
Expand Down Expand Up @@ -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 );
}

/**
Expand All @@ -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 );
}


Expand Down
24 changes: 20 additions & 4 deletions src/coreComponents/constitutive/solid/DruckerPrager.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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 ),
Expand Down Expand Up @@ -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 );
}

/**
Expand All @@ -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 );
}


Expand Down
Loading
Loading