Damage Model

Damage represents the progressive degradation of the material’s load-carrying capability due to microcracks, voids, or accumulated plastic deformation. As damage grows, the effective load-bearing area decreases, which increases local stresses and reduces the effective stiffness of the material. This degradation is incorporated into the constitutive model by scaling the elastic stiffness tensor and the critical resolved shear stress using a scalar damage variable (D in [0,1]), where (D=0) corresponds to an undamaged state and (D=1) represents complete loss of load-bearing capacity.

Key Classes and Concepts

• Damage — reads and stores damage properties; selects between the Creep and PEEQ damage models via $DamageModel.
• Damage Creep — stress-based creep damage (extension of the Cailletaud model; Kachanov's law for creep). See Damage Creep below.
• Damage PEEQ — plastic-strain-based damage; thresholds minPEEQ_$\alpha$ and minPEEQ_$\alpha$ per phase. See Damage PEEQ below.

Effective local Damage update the Critical Resolved Shear Stress ${\tau }_c$ and Elastic Stiffness $C$ as follows:

$$\begin{array}{}\text{(1)}& {\tau }_c^{eff}=(1-D){\tau }_c\end{array}$$

$$\begin{array}{}\text{(2)}& C^{eff}=(1-D)C\end{array}$$

General Damage Input Parameters

Parameter	Description
DamageModel	Damage Model Selection: Creep or PEEQ
DamageSwitch ($\alpha$)	Damage model activation switch for each phase: True or False
SmoothFlag	Allowed smoothening for localized damage: True or False. Default value: False
maxDamageAllowed	Max Value of allowed damage between 0-1. Default value: 0.95

Damage Creep

In order to capture damage evolution, an extension of the Cailletaud model [33] proposed in Reference 34 was used. The resolved shear stress was assumed to act as a driving force for damage accumulation, which is in accordance with Kachanov’s damage law for creep as follows:

$$\begin{array}{}\text{(3)}& {\dot{D}}_s=A({\tau }_s{)}^n\end{array}$$

where ${\dot{D}}_s$ is a measure for damage evolution in slip system $s$. $A$ and $s$ are additional material parameters. The resolved shear stresses acting in each slip system have to be corrected. This was achieved by double contracting the global damage stress with the Schmid tensor $M$:

$$\begin{array}{}\text{(4)}& {\tau }_s=\frac{1}{1-D_s}\sigma :\mathbf{M}\end{array}$$

In the preceding equation, $D_s$ is the scalar damage parameter for slip system $s$. In the current case, $D_s=0$ corresponds to an undamaged material state for slip system $s$, and $D_s=1$ represents end of life for that slip system.

Parameter	Description
$A_{\alpha }$	Material constant for each phase $\alpha$
$n_{\alpha }$	Stress exponent for each phase $\alpha$

Damage PEEQ

A damage model is used to represent the microscopic pores and cracks on a mesoscopic scale. Microcracks and pores effectively reduce the load bearing area of a hypothetical sample, thereby increasing the stresses acting on the remaining area. This effect is translated into the following model.

$$\begin{array}{}\text{(5)}& \overline{p}=\sqrt{\frac{2}{3}{\epsilon }_{pl}:{\epsilon }_{pl}}\end{array}$$

The locally accumulated plastic strain ${\epsilon }_{pl}$ is used to calculate the damage parameter $D$ (Eq. 1.30). $D$ is 0 for a plastic strain measure $\overline{p}$ below the minimum plastic strain $p_{min}$, between $p_{min}$ and $p_{max}$ it scales linearly with the plastic strain measure, and above $p_{max}$ it is 1.

$$\begin{array}{}\text{(6)}& D=\{\begin{array}{ll}0& \text{for }\overline{p}<p_{min}\\ \frac{\overline{p}-p_{min}}{p_{max}-p_{min}}& \text{for }{p}_{min}\le \overline{p}\le p_{max}\\ 1& \text{for }\overline{p}>p_{max}\end{array}\end{array}$$

The damage parameter $D$ is used to adjust the stiffness tensor $C$ and stress $\tau$ for the change in load bearing area. The reduction of the stiffness tensor with increasing damage results in effective softening of the material. If a plasticity model is applied in conjunction with the damage model, the critical resolved shear stress ${\tau }_c^{eff}$ has to be adjusted accordingly. This is the case because a softer material with the same critical resolved shear stress would still yield at the same stress, resulting in elastic softening but an unchanged plastic behavior. Thus, the critical resolved shear stress is also modified.

Parameter	Description
$p_{min,\alpha }$	Minimum plastic strain for damage initiation for each phase $\alpha$
$p_{max,\alpha }$	Maximum plastic strain for complete damage for each phase $\alpha$

---

Usage

Input

The full @Damage and @CrystalPlasticity blocks are listed below in Input Parameters used in the above example. The minimum set of keys is:

Key	Meaning
$DamageModel	Creep or PEEQ
$DamageSwitch_<α>	enable damage on phase $\alpha$ (Yes / No)
$SmoothFlag	smooth localised damage (default No)
$MaxDamage	maximum admissible damage value (default 0.95)

For Creep: $A_<α>, $n_<α>. For PEEQ: $minPEEQ_<α>, $maxPEEQ_<α>.

Output

• Damage evolution Solve — Damage::Solve(Phase, EP, CP, BC, dt) recomputes the damage field each call:

  void Damage::Solve(PhaseField& Phase, ElasticProperties& EP,
                     CrystalPlasticity& CP, BoundaryConditions& BC,
                     double dt);

• VTK output with optional distorted-grid rendering:

  void WriteVTK(const Settings& locSettings, const PhaseField& Phase,
                const ElasticProperties& EP, const int tStep,
                const int precision = 6, bool Distorted = false);

• Raw / restart I/O:

  bool Read (const Settings& locSettings, const BoundaryConditions& BC, const int tStep) override;
  bool Write(const Settings& locSettings, const int tStep) const override;

Example

This example simulates polycrystalline material under mechanical loading with crystal plasticity and damage. A Voronoi tessellation seeds the polycrystal; the time loop applies strain-rate loads, solves crystal plasticity, then updates damage. Stress, strain, and damage metrics are written as VTK and CSV.

This example demonstrates the simulation of polycrystalline material behaviour under mechanical loading, incorporating crystal plasticity and damage evolution. The simulation initializes a polycrystalline microstructure using Voronoi tessellation, stabilizes the interfaces, and then enters a time loop where mechanical loads as strain rates are applied, damage is computed, and crystal plasticity is solved. The results, including phase fields, stresses, strains, and damage metrics, are outputted in VTK format for visualization and analysis. Tabulated data for stress-strain and plastic strain evolution are outputted to text files for plotting purposes.

#include "Mechanics.h"
#include "RunTimeControl.h"
#include "DoubleObstacle.h"
#include "PhaseField.h"
#include "Initializations.h"
#include "BoundaryConditions.h"
#include "InterfaceProperties.h"
#include "Tools/TimeInfo.h"
#include "TextOutput.h"
#include "MechanicalLoads.h"

using namespace std;
using namespace openphase;

int main(int argc, char *argv[])
{
#ifdef MPI_PARALLEL
    int provided = 0;
    MPI_Init_thread(&argc, &argv, MPI_THREAD_FUNNELED, &provided);
    MPI_Comm_rank(MPI_COMM_WORLD, &MPI_RANK);
    MPI_Comm_size(MPI_COMM_WORLD, &MPI_SIZE);
#endif

    Settings                    OPSettings;
    OPSettings.ReadInput();

    RunTimeControl              RTC  (OPSettings);
    PhaseField                  Phi  (OPSettings);
    DoubleObstacle              DO   (OPSettings);
    InterfaceProperties         IP   (OPSettings);
    BoundaryConditions          BC   (OPSettings);
    ElasticProperties           EP   (OPSettings);
    ElasticitySolverSpectral    ES   (OPSettings);
    CrystalPlasticity           CP   (OPSettings);
    MechanicalLoads             ML   (OPSettings);
    Damage                      DP   (OPSettings);
    TimeInfo                    Timer(OPSettings, "Execution Time Statistics");

    size_t GrainsPhase = 0;
    int number_of_grains = 10;
    Initializations::VoronoiTessellation(Phi, BC, number_of_grains, GrainsPhase);

    for(size_t LocalStep = 0; LocalStep <= 100; LocalStep++)
    {
        IP .Set                          (Phi, BC);
        DO .CalculatePhaseFieldIncrements(Phi, IP, dG);
        Phi.NormalizeIncrements          (BC, RTC.dt);
        Phi.MergeIncrements              (BC, RTC.dt);

    }

    string StressStrainFile  = DefaultTextDir + "StressStrainFile.txt";
    string PlasticStrainFile = DefaultTextDir + "PlasticStrainFile.txt";
    remove(StressStrainFile. c_str());
    remove(PlasticStrainFile.c_str());

    EP.SetEffectiveProperties(Phi);
    CP.SetInitialHardening(Phi, BC);
    CP.SetGrainsProperties(Phi, EP);

//-------------------------------------------------------------------------------------------//
    ConsoleOutput::WriteLineInsert("Entering the Time Loop!!!");
//-------------------------------------------------------------------------------------------//
    for(RTC.tStep = RTC.tStart; RTC.tStep <= RTC.nSteps; RTC.IncrementTimeStep())
    {
        EP.SetEffectiveProperties(Phi, DP);

        ML.Activate(Phi, EP, RTC);
        ML.Apply(Phi, EP, RTC);

        DP.Solve(Phi, EP, CP, BC, RTC.dt);
        CP.Solve(Phi, BC, EP, ES, DP, RTC.dt);

        if (RTC.WriteVTK())
        {
            Phi.WriteVTK             (OPSettings, RTC.tStep);
            DP.WriteVTK              (OPSettings, Phi, EP,  RTC.tStep);
        }
        if (RTC.WriteToScreen())
        {
            double I_En = DO.AverageEnergyDensity(Phi, IP);
            std::string message  = ConsoleOutput::GetStandard("Interface energy density", I_En);
            ConsoleOutput::WriteTimeStep(RTC, message);
            Timer.PrintWallClockSummary();

            ConsoleOutput::WriteStandard("Average Strain", EP.AverageStrain);
            ConsoleOutput::WriteStandard("Plastic Strain", EP.AveragePlasticStrain());
            ConsoleOutput::WriteLine();
            ConsoleOutput::WriteStandard("Max Damage", DP.MaxDamage);
            ConsoleOutput::WriteStandard("Average Damage", DP.AverageDamage(Phi));
            ConsoleOutput::WriteLine("#");
            EP.WriteStressStrainData(RTC.tStep, RTC.SimulationTime, StressStrainFile);
            CP.WriteAveragePlasticStrainCSV(EP, PlasticStrainFile, EP.AverageStrain[0]);
        }
    }
#ifdef MPI_PARALLEL
    MPI_Finalize ();
#endif
    return 0;
}

Input Parameters used in the above example

Standard Open Phase Input File
!!!All values in MKS (or properly scaled) units please!!!

----------------------------------------------------------------------------
@Damage
----------------------------------------------------------------------------

Model Types : [ CREEP, PEEQ ]

$DamageModel      |  Damage model being used             |  : CREEP 
$MaxDamage        |  Max allowable damage                |  : 0.90 
$DamageSwitch_0   |  Enable/ disable damage model        |  : Yes 
$SmoothFlag       |  Enable / disable Damage smoothening |  : Yes

|--------------| CREEP model parameters| ----------------------------|

// --- in each phase, here only one phase is used which is 0

$A_0        Damage Rate value                      : 1e-18
$n_0        Exponent value                         : 2

|--------------| PEEQ model parameters| -----------------------------|

$minPEEQ_0  Minimum value of PEEQ                  : 0.02
$maxPEEQ_0  Minimum value of PEEQ                  : 0.2

----------------------------------------------------------------
--------------------| End of File |-- --------------------------
----------------------------------------------------------------

Dependencies

• Mechanics — continuum-mechanics and crystal-plasticity theory.
• Elastic Properties — the stiffness tensor that the damage scalar modifies.
• PhaseField — fields and increments the energy / mobility values are evaluated against. BoundaryConditions — applied throughout the simulation Settings — phase list and grid metadata.