# Inelastic Fracture

We have developed and implemented the {ref}`viscoplasticity-rheological-models` in the figure.
These have been obtained via the parallel assembly of the rheological branches described in previous sections.
We have further introduced a rheological "fracture element" that is assembled in series and described via the phase-field modeling of fracture.

::: Note
As shown below, the Prandtl element is limited to perfect plasticity so that only the values of $E_{prandtl}$ and $\sigma_{0}$ are required.
Linear strain hardening can be assigned to the Hooke branch and non-linear effects captured using viscoelastic Maxwell branch(es).
:::

```{figure} /doc/img/viscoplasticity-rheological-models.png
---
align: center
name: viscoplasticity-rheological-models
---
Rheological models of inelastic fracture
```

## Implementation

The stress update for the Maxwell and Prandtl branches is carried out using degraded trial quantities and the admissibility check for perfect plasticity evaluated with yield surface $f = q^{\mathrm{tr}}(\phi) - \sigma_y = 0$.
Expressions in Eq. {eq}`viscoelasticity-deviatoric-strain-update-hencky` and Eq. {eq}`plasticity-deviatoric-strain-update` become

$$
\begin{aligned}
  \boldsymbol{\varepsilon}_\mathrm{dev}^{e}
   &=
  \begin{cases}
    \left(1 + \mu_{degr} \frac{\Delta t}{\eta_\mathrm{dev}}\right)^{-1}
    \boldsymbol{\varepsilon}_\mathrm{dev}^{e\,\mathrm{tr}},
    \qquad \mathrm{MAXWELL}, \\[6pt]
    \left(1 - 3 \mu_{degr}\frac{\Delta\gamma(\phi)}{q^{\mathrm{tr}}(\phi)}\right)
    \boldsymbol{\varepsilon}_\mathrm{dev}^{e\,\mathrm{tr}}
    =
    \left(\frac{\sigma_0}{ g(\phi) q^{\mathrm{tr}}}\right)
    \boldsymbol{\varepsilon}_\mathrm{dev}^{e\,\mathrm{tr}},
    \qquad \mathrm{PRANDTL - Perfect\, plastic \,regime}.
  \end{cases}
\end{aligned}
$$

where $\mu_{degr} =g(\phi)\mu$.
Notably, for a Prandtl branch in the perfect plastic regime, the deviatoric stresses are independent of $\phi$,

$$
\begin{aligned}
  \boldsymbol{\tau}_{\mathrm{dev} \, \mathrm{degr}}
  =
g(\phi)\left(
  \frac{\sigma_0}{g(\phi)\,q^{\mathrm{tr}}}\boldsymbol{\tau}_{\mathrm{dev}}^{\mathrm{tr}}
  \right)
  =
  \left(\frac{\sigma_0}{q^{\mathrm{tr}}}\right)\boldsymbol{\tau}_{\mathrm{dev}}^{\mathrm{tr}}  \qquad \mathrm{PRANDTL - Perfect\, plastic \,regime},
\end{aligned}
$$

consistent with the serial assembly of the rheological fracture element in the proposed rheological framework.

## Consistent tangent for the monolithic scheme

From Eq. {eq}`linear-degraded-stress`, the linearization of the
degraded Kirchhoff stress,

$$
\begin{aligned}
{\diff}\boldsymbol{\tau}_{\mathrm{degr}}
=
{\diff}_{\boldsymbol{u}}\boldsymbol{\tau}_{\mathrm{degr}}
+
{\diff}_{\phi}\boldsymbol{\tau}_{\mathrm{degr}},
\end{aligned}
$$

reads

$$
\begin{aligned}
  {\diff}\boldsymbol{\tau}_{\mathrm{degr}}
  &=
  \begin{cases}
    g(\phi)\,{\diff}_{\boldsymbol{u}}\boldsymbol{\tau}
    + g'(\phi)\,{\diff}\phi\,\boldsymbol{\tau}
    + g(\phi)\,{\diff}_{\phi}\boldsymbol{\tau}_{\mathrm{dev}},
    \qquad \mathrm{if } \operatorname{trace}\boldsymbol{\varepsilon} \ge 0, \\[6pt]
    {\diff}\boldsymbol{\tau}_{\mathrm{vol}}
    + g(\phi)\,{\diff}_{\boldsymbol{u}}\boldsymbol{\tau}_{\mathrm{dev}}
    + g'(\phi)\,{\diff}\phi\,\boldsymbol{\tau}_{\mathrm{dev}}
    + g(\phi)\,{\diff}_{\phi}\boldsymbol{\tau}_{\mathrm{dev}},
    \qquad \mathrm{if } \operatorname{trace}\boldsymbol{\varepsilon} < 0 .
  \end{cases}
\end{aligned}
$$

For the viscoelastic (Maxwell) branch,

$$
\begin{aligned}
  {\diff}_{\phi}\boldsymbol{\tau}_{\mathrm{dev}}
  = -2\mu \frac{\mu \Delta t}{\eta_\mathrm{dev}}\frac{g'(\phi)}{\left(1 + \frac{\mu \Delta t}{\eta_\mathrm{dev}}\, g(\phi)\right)^2} \boldsymbol{\varepsilon}^{e\,\mathrm{tr}}_\mathrm{dev}
  = -\frac{\mu \Delta t}{\eta_\mathrm{dev}} \frac{g'(\phi)}{1 + \dfrac{\mu \Delta t}{\eta_\mathrm{dev}} g(\phi)}\boldsymbol{\tau}_\mathrm{dev}
\end{aligned}
$$

For the elastoplastic (Prandtl) element in a perfect plasticity regime,

$$
\begin{aligned}
  {\diff}_{\phi}\boldsymbol{\tau}_{\mathrm{dev}}
  = -2\mu \frac{\sigma_0}{q^{\mathrm{tr}}}\,\frac{g'(\phi)}{g(\phi)^2}\,\boldsymbol{\varepsilon}^{e\,\mathrm{tr}}_\mathrm{dev}
  = -\frac{g'(\phi)}{g(\phi)}\boldsymbol{\tau}_\mathrm{dev}
\end{aligned}
$$

which gives

$$
\begin{aligned}
    {\diff}_{\phi} \boldsymbol{\tau}_{\mathrm{dev} \, \mathrm{degr}} = g'(\phi)\boldsymbol{\tau}_{\mathrm{dev}}
+ g(\phi)\diff_{\phi} \boldsymbol{\tau}_{\mathrm{dev}}
  =\boldsymbol{0}
 \qquad \mathrm{PRANDTL - Perfect\, plastic \,regime}
\end{aligned}
$$

consistent with the Prandtl contribution to $\boldsymbol{\tau}_{\mathrm{dev}}$ being independent of $\phi$ in such regime.

Unlike the (Hooke) branch case, for which ${\diff}_{\phi}\psi^{+} = 0$ and thus omitted in Eq. {eq}`linear-Psi-plus-linearization`, we obtain

$$
  \begin{aligned}
  {\diff}_{\phi}\psi^{+} = {\diff}_\phi \psi_{\mathrm{dev}}
  &= \boldsymbol{\tau}_\mathrm{dev}: \frac{\partial \boldsymbol{\varepsilon}_\mathrm{dev}^e}{\partial \phi} \\
  &=-\boldsymbol{\tau}_\mathrm{dev}: \frac{\mu \Delta t}{\eta_\mathrm{dev}}{g'(\phi)}{\left(1 + \frac{\mu \Delta t}{\eta_\mathrm{dev}}\, g(\phi)\right)^{-2}} \boldsymbol{\varepsilon}^{e \,\mathrm{tr}}_\mathrm{dev} \\
  &=-\boldsymbol{\tau}_\mathrm{dev}: \frac{\mu \Delta t}{\eta_\mathrm{dev}}{g'(\phi)}{\left(1 + \frac{\mu \Delta t}{\eta_\mathrm{dev}}\, g(\phi)\right)^{-1}} \boldsymbol{\varepsilon}^{e}_\mathrm{dev} \\
  &= - 2\frac{\mu \Delta t}{\eta_\mathrm{dev}}{g'(\phi)}{\left(1 + \frac{\mu \Delta t}{\eta_\mathrm{dev}}\, g(\phi)\right)^{-1}} \psi_{\mathrm{dev}}
  \end{aligned}
$$

and

$$
  \begin{aligned}
  {\diff}_{\phi}\psi^{+} = {\diff}_\phi \psi_{\mathrm{dev}}
  &= \boldsymbol{\tau}_\mathrm{dev}: \frac{\partial \boldsymbol{\varepsilon}_\mathrm{dev}^e}{\partial \phi} \\
  &=\boldsymbol{\tau}_\mathrm{dev}: \left[ - g'(\phi)\frac{\sigma_0}{g(\phi)^2 q^{\mathrm{tr}}}\right] \boldsymbol{\varepsilon}^{e \,\mathrm{tr}}_\mathrm{dev} \\
  &=\boldsymbol{\tau}_\mathrm{dev}:  - \frac{g'(\phi)}{g(\phi)} \boldsymbol{\varepsilon}^{e}_\mathrm{dev} \\
  &= -2 \frac{g'(\phi)}{g(\phi)} \psi_\mathrm{dev}
  \end{aligned}
$$

for the viscoelastic (Maxwell) branch and the elastoplastic (Prandtl) element in a perfect plasticity regime, respectively.

### Command-line interface

To enable the Hencky viscoelastic, elastoplastic, and visco-elastoplastic models with fracture, respectively use the model option (for current conf.): `-model viscolasticity-hencky-damage-principal-current`, `-model plasticity-hencky-damage-principal-current`, and, `-model viscoplasticity-hencky-damage-principal-current`.
Set the corresponding material parameters listed in {ref}`viscoelastic-hencky-damage-options`, {ref}`plastic-hencky-damage-options`, {ref}`viscoplastic-hencky-damage-options`.
Any parameter without a default option is required.

(viscoelastic-hencky-damage-options)=

:::{list-table} Hencky viscoelasticity model with damage options
:header-rows: 1
:widths: 30 55 15

* * Option
  * Description
  * Default Value

* * `-model viscoelasticity-hencky-damage-current | viscoelasticity-hencky-damage-current`
  * Required to enable the Hencky viscoelasticity with damage model, `|` denotes an option for the user.
  *

* * `-E_hooke [real]`
  * [Young's modulus in Hooke branch (i.e. a single spring)](https://en.wikipedia.org/wiki/Young%27s_modulus), $E_{hooke} > 0$
  *

* * `-nu_hooke [real]`
  * [Poisson's ratio in Hooke branch](https://en.wikipedia.org/wiki/Poisson%27s_ratio), $\nu_{hooke} \leq 0.5$.
  *

* * `-E_maxwell [real]`
  * [Young's modulus in Maxwell branch](https://en.wikipedia.org/wiki/Young%27s_modulus), $E_{maxwell} > 0$
  *

* * `-nu_maxwell [real]`
  * [Poisson's ratio in Maxwell branch](https://en.wikipedia.org/wiki/Poisson%27s_ratio), $\nu_{maxwell} \leq 0.5$.
  *

* * `-deviatoric_viscosity [real]`
  * [Deviatoric viscosity for Maxwell branch], $\eta_\mathrm{dev} >= 0$.
  *

* * `-fracture_toughness [real]`
  * [Fracture toughness expressed as critical energy release rate](https://en.wikipedia.org/wiki/Energy_release_rate_(fracture_mechanics)), $G_c > 0$
  *

* * `-characteristic_length [real]`
  * Length scale parameter of phase-field model, $l_0 > 0$
  *

* * `-residual_stiffness [real]`
  * Residual stiffness in full fracture state ($\phi = 1$), $\eta > 0$
  *

* * `-damage_viscosity [real]`
  * Viscosity for viscous regularization of damage problem, $\zeta >= 0$
  *

* * `-damage_scaling [real]`
  * Scaling parameter, $c_0 > 0$
  *

* * `-use_AT1 [bool]`
  * True: AT1 model, False: AT2
  * True

* * `-use_offdiagonal [bool]`
  * True: full Jacobian, False: Jacobian diagonal terms only (for quasi-Newton methods),
  * True

:::

An example using the Hencky viscoelasticity with damage model can be run via

```console
./bin/ratel-quasistatic -options_file examples/ymls/ex02/viscoelasticity-hencky-damage-initial-principal.yml
```

(plastic-hencky-damage-options)=

:::{list-table} Hencky plasticity model with damage options
:header-rows: 1
:widths: 30 55 15

* * Option
  * Description
  * Default Value

* * `-model plasticity-hencky-damage-current | plasticity-hencky-damage-current`
  * Required to enable the Hencky plasticity with damage model, `|` denotes an option for the user.
  *

* * `-E_hooke [real]`
  * [Young's modulus in Hooke branch (i.e. a single spring)](https://en.wikipedia.org/wiki/Young%27s_modulus), $E_{hooke} > 0$
  *

* * `-nu_hooke [real]`
  * [Poisson's ratio in Hooke branch](https://en.wikipedia.org/wiki/Poisson%27s_ratio), $\nu_{hooke} \leq 0.5$.
  *

* * `-E_prandtl [real]`
  * [Young's modulus in Prandtl branch](https://en.wikipedia.org/wiki/Young%27s_modulus), $E_{prandtl} > 0$
  *

* * `-nu_prandtl [real]`
  * [Poisson's ratio](https://en.wikipedia.org/wiki/Poisson%27s_ratio), $\nu_{prandtl} \leq 0.5$.
  *

* * `-yield_stress [real]`
  * [Initial yield stress threshold](https://en.wikipedia.org/wiki/Plasticity_(physics)), $\sigma_0 > 0$
  *

* * `-fracture_toughness [real]`
  * [Fracture toughness expressed as critical energy release rate](https://en.wikipedia.org/wiki/Energy_release_rate_(fracture_mechanics)), $G_c > 0$
  *

* * `-characteristic_length [real]`
  * Length scale parameter of phase-field model, $l_0 > 0$
  *

* * `-residual_stiffness [real]`
  * Residual stiffness in full fracture state ($\phi = 1$), $\eta > 0$
  *

* * `-damage_viscosity [real]`
  * Viscosity for viscous regularization of damage problem, $\zeta >= 0$
  *

* * `-damage_scaling [real]`
  * Scaling parameter, $c_0 > 0$
  *

* * `-use_AT1 [bool]`
  * True: AT1 model, False: AT2
  * True

* * `-use_offdiagonal [bool]`
  * True: full Jacobian, False: Jacobian diagonal terms only (for quasi-Newton methods),
  * True

:::

An example using the Hencky plasticity with damage model can be run via

```console
./bin/ratel-quasistatic -options_file examples/ymls/ex02/plasticity-hencky-damage-initial-principal.yml
```

(viscoplastic-hencky-damage-options)=

:::{list-table} Hencky viscoplasticity model with damage options
:header-rows: 1
:widths: 30 55 15

* * Option
  * Description
  * Default Value

* * `-model viscoplasticity-hencky-damage-current | viscoelasticity-hencky-damage-current`
  * Required to enable the Hencky viscoplasticity with damage model, `|` denotes an option for the user.
  *

* * `-E_hooke [real]`
  * [Young's modulus in Hooke branch (i.e. a single spring)](https://en.wikipedia.org/wiki/Young%27s_modulus), $E_{hooke} > 0$
  *

* * `-nu_hooke [real]`
  * [Poisson's ratio in Hooke branch](https://en.wikipedia.org/wiki/Poisson%27s_ratio), $\nu_{hooke} \leq 0.5$.
  *

* * `-E_maxwell [real]`
  * [Young's modulus in Maxwell branch](https://en.wikipedia.org/wiki/Young%27s_modulus), $E_{maxwell} > 0$
  *

* * `-nu_maxwell [real]`
  * [Poisson's ratio in Maxwell branch](https://en.wikipedia.org/wiki/Poisson%27s_ratio), $\nu_{maxwell} \leq 0.5$.
  *

* * `-E_prandtl [real]`
  * [Young's modulus in Prandtl branch](https://en.wikipedia.org/wiki/Young%27s_modulus), $E_{prandtl} > 0$
  *

* * `-nu_prandtl [real]`
  * [Poisson's ratio](https://en.wikipedia.org/wiki/Poisson%27s_ratio), $\nu_{prandtl} \leq 0.5$.
  *

* * `-deviatoric_viscosity [real]`
  * [Deviatoric viscosity for Maxwell branch], $\eta_\mathrm{d} >= 0$.
  *

* * `-yield_stress [real]`
  * [Initial yield stress threshold](https://en.wikipedia.org/wiki/Plasticity_(physics)), $\sigma_0 > 0$
  *

* * `-fracture_toughness [real]`
  * [Fracture toughness expressed as critical energy release rate](https://en.wikipedia.org/wiki/Energy_release_rate_(fracture_mechanics)), $G_c > 0$
  *

* * `-characteristic_length [real]`
  * Length scale parameter of phase-field model, $l_0 > 0$
  *

* * `-residual_stiffness [real]`
  * Residual stiffness in full fracture state ($\phi = 1$), $\eta > 0$
  *

* * `-damage_viscosity [real]`
  * Viscosity for viscous regularization of damage problem, $\zeta >= 0$
  *

* * `-damage_scaling [real]`
  * Scaling parameter, $c_0 > 0$
  *

* * `-use_AT1 [bool]`
  * True: AT1 model, False: AT2
  * True

* * `-use_offdiagonal [bool]`
  * True: full Jacobian, False: Jacobian diagonal terms only (for quasi-Newton methods),
  * True

:::

An example using the Hencky viscoplasticity with damage model can be run via

```console
./bin/ratel-quasistatic -options_file examples/ymls/ex02/viscoplasticity-hencky-damage-initial-principal.yml
```