An a priori irreversible phase-field formulation for ductile fracture at finite strains based on the Allen–Cahn theory: a variational approach and FE-implementation

Abstract

In this paper, a new crack surface energy for the simulation of ductile fracture is proposed, which is based on the Allen–Cahn theory of diffuse interfaces. In contrast to existing fracture approaches, here, the crack surface energy density is a double-well potential based on a new interpretation of the crack surface. That is, the energy associated with the whole diffuse region between the fully cracked and intact regions is interpreted as crack surface energy. This kind of formulation, on the one hand, results in the balance of micromechanical forces and on the other hand, is a priori thermodynamically consistent. Furthermore, the proposed formulation is based on a gamma-convergent interface energy and it is in agreement with the classical solution of Irwin (Appl Mech Trans ASME E24:351–369, 1957). It is shown that in contrast to existing models, crack irreversibility is automatically fulfilled and no further constraints related to neither local nor global irreversibility are needed. To also account for potential plastic shear band localization, the approach is extended by a micromorphic plasticity model. By analyzing two different classical numerical benchmark problems, the proposed formulation is shown to enable mesh-independent results which are in agreement with the results of competing approaches.

Appendices

Appendix

A Continuum mechanical framework for ductile materials

Here, the continuum mechanical basis for the phase field model proposed in this paper as well as the considered elasto-plastic material model are described. Generally, a fictitiously undamaged strain energy density functional $\psi _0$, defined per unit reference volume, can be written as a function of a deformation gradient $\varvec{F}=\text {Grad}_{\varvec{X}}(\varvec{\varphi })=\nabla _{\varvec{X}}\varvec{\varphi }$. The deformation mapping $\varvec{\varphi }$ maps the point $\varvec{X}\in \Omega$ in the reference (undeformed) configuration to the point $\varvec{x}\in \varvec{\varphi }(\Omega )$ in the current (deformed) configuration. In addition, the plastic deformation is an inelastic stress-free process. Therefore, an intermediate stress-free configuration is considered, which results in a multiplicative decomposition of the deformation gradient into a plastic and elastic one, i.e.,

$$\begin{aligned} \varvec{F} =\varvec{F}^{\textrm{e}}\cdot \varvec{F}^{\textrm{p}}, \end{aligned}$$

(55)

cf. [40]. Based on this decomposition, it is clear that the plastic deformation gradient has no contribution to the stored elastic energy. Furthermore, in some materials such as metals, the plastic deformation is volume preserving, i.e., det$(\varvec{F^{\textrm{p}}})= 1$. Therefore, the plastic incompressibility condition has to be satisfied. The strain energy density can thus be decomposed additively into an elastic and plastic term

$$\begin{aligned} \psi _0(\varvec{F},\varvec{F}^{\textrm{p}},\alpha ) =\psi _0^{\textrm{e}}(\varvec{F},\varvec{F}^{\textrm{p}})+\psi _0^{\textrm{p}}(\alpha ). \end{aligned}$$

(56)

Herein, $\alpha$ is a strain-like internal variable. To fulfill objectivity, the elastic part is rewritten as a function of the elastic right Cauchy-Green tensor $\varvec{C}^{\textrm{e}}=\varvec{F}^{\mathrm {e^T}}\cdot \varvec{F}^{\textrm{e}} =\varvec{F}^{\mathrm {p^{-T}}}\cdot \varvec{C}\cdot \varvec{F}^{\mathrm {p^{-1}}}$ with the right Cauchy-Green tensor $\varvec{C}:=\varvec{F}^{\textrm{T}}\cdot \varvec{F}$, such that $\psi _0^{\textrm{e}}=\psi _0^{\textrm{e}}(\varvec{F}^{\mathrm {e^{T}}}\cdot \varvec{F}^{\textrm{e}}) = \psi _0(\varvec{C}^{\textrm{e}})$. By means of the postulate of minimum total potential energy and assuming conservative external forces, the deformation mapping is computed by

$$\begin{aligned} \varvec{\varphi }=\text {arg}\,\Biggl \{\text {inf}\left( \int _\Omega \psi _0\left( \varvec{C}^{\textrm{e}},\alpha \right) \,{\text {d}}V-\int _{\Omega }\rho _0\,\varvec{b}\cdot \varvec{\varphi } {\text {d}}V-\int _{\partial _{\varvec{N}}\Omega }\varvec{T}\cdot \varvec{\varphi }\,{\text {d}}A\right) \Biggr \}, \end{aligned}$$

(57)

where $\rho _0$, $\rho _0\varvec{b}$, $\varvec{T}$ are the material density in reference configuration, body forces and the traction vector, respectively. Considering the second law of thermodynamics for isothermal conditions

$$\begin{aligned} \mathscr {D}=\varvec{P}_0:\dot{\varvec{F}}-\dot{\psi }_0\ge 0, \end{aligned}$$

(58)

with the first Piola–Kirchhoff stress tensor $\varvec{P}_0$. Inserting the time derivative of $\psi _0$ yields

$$\begin{aligned} \left( \varvec{P}_0 - \dfrac{\partial \psi _0}{\partial \varvec{F}}\right) :\dot{\varvec{F}} - \dfrac{\partial \psi _0}{\partial \varvec{F}^{\textrm{p}}}:\dot{\varvec{F}}^{\textrm{p}} - \dfrac{\partial \psi _0}{\partial \alpha }\dot{\alpha } \ge 0. \end{aligned}$$

(59)

Applying the Coleman–Noll procedure leads to the constitutive equation for the first Piola–Kirchhoff stress tensor

$$\begin{aligned} \varvec{P}_0=\dfrac{\partial \psi _0^{\textrm{e}}}{\partial {\varvec{F}}} = \varvec{P}_0^{\textrm{e}}\cdot \varvec{F}^{\textrm{p}^{\mathrm {-T}}} \end{aligned}$$

(60)

wherein $\varvec{P}^{\textrm{e}}_0:=\partial _{\varvec{F}^{\textrm{e}}}\psi _0^{\textrm{e}}$ and $\partial _{\varvec{F}}\varvec{F}^{\textrm{e}} = \varvec{F}^{\textrm{p}^{\mathrm {-T}}}$. Defining the Mandel stress tensor [44] in the intermediate configuration $\varvec{\Sigma }^{\textrm{e}}=\varvec{F}^{\mathrm {e^T}}\cdot \varvec{P}^{\textrm{e}}$ and the plastic spatial velocity gradient $\varvec{L}^{\textrm{p}}=\dot{\varvec{F}^{\textrm{p}}}\cdot \varvec{F}^{{\textrm{p}}^{-1}}$ yields the reduced dissipation inequality

$$\begin{aligned} \mathscr {D}=\varvec{\Sigma }^{\textrm{e}}:\varvec{L}^{\textrm{p}} -Q\dot{\alpha }\ge 0, \end{aligned}$$

(61)

where Q is the thermodynamic force of the internal variable $\alpha$, i.e.,

$$\begin{aligned} Q:= \dfrac{\partial \psi _0}{\partial \alpha }=\dfrac{\partial \psi _0^{\textrm{p}}}{\partial \alpha }. \end{aligned}$$

(62)

For the description of elasto-plasticity, an admissible domain is defined in terms of the yield function $\phi \le 0$. If the yield function is less than zero, i.e., $\phi <0$, no plastic deformation takes place [43]. Furthermore, $\phi =0$ represents the yield surface, where plastic deformations evolve. In case of metal plasticity with isotropic hardening, this state function is considered to have the form

$$\begin{aligned} \phi (\varvec{\Sigma ^{\textrm{e}}},Q):=||dev\varvec{\Sigma }^{\textrm{e}}||-Q. \end{aligned}$$

(63)

Application of the principle of maximum dissipation and considering an associative flow rule, the evolution equations are obtained as

$$\begin{aligned} \varvec{L}^{\textrm{p}}=\dot{\lambda }\varvec{N}, \,\,\,\,\, \dot{\lambda } = \dot{\alpha }, \end{aligned}$$

(64)

where $\varvec{N}:=\partial _{\varvec{\Sigma }^{\textrm{e}}}\phi$ is a second-order tensor imposing a constraint on the direction of plastic deformation and $\lambda$ is the plastic multiplier. Applying the principle of maximum dissipation yields the well-known Kuhn-Tucker conditions

$$\begin{aligned} \dot{\lambda }\ge 0, \,\,\,\,\, \dot{\lambda }\dot{\phi }\ge 0, \end{aligned}$$

(65)

which need to be fulfilled. In order to fulfill the plastic incompressibility condition observed in metals, exponential integration is often used, cf. de Souza Neto [14]. An alternative can be obtained by noting that the plastic incompressibility condition is automatically fulfilled if $\varvec{N}$ is traceless [27]. In our model this will be the case. As a specification of the model, we consider a Neo-Hookean type energy density function

$$\begin{aligned} \psi _0^{\textrm{e}} = \dfrac{\mu }{2}\left( {\text {tr}}\left( \varvec{C}^{\textrm{e}}\right) -3\right) +\lambda \dfrac{J^{{\textrm{e}}^{2}} -1}{4} -\left( \dfrac{\lambda }{2}+\mu \right) \ln (J^{\textrm{e}}), \,\,\,\,\,\,\, J^{\textrm{e}}={\text {det}}(\varvec{F}^{\textrm{e}})>0. \end{aligned}$$

(66)

The plastic energy density is assumed to have the form [68]

$$\begin{aligned} \psi _0^{\textrm{p}}=\dfrac{1}{2} h \alpha ^2 + (\sigma _{\infty }-\sigma _y)\left[ \alpha +\dfrac{\exp (-\zeta \alpha )-1}{\zeta }\right] + \sigma _y\alpha , \end{aligned}$$

(67)

considering a saturation parameter $\zeta$, a linear hardening parameter h, the initial yield stress $\sigma _y$, and a further parameter $\sigma _{\infty }$, which is associated with the modeling of the transition between an initially negative exponential hardening to a linear hardening. Note that the method proposed in this paper is not restricted to this choice of $\psi _0^{\textrm{e}}$ and $\psi _0^{\textrm{p}}$.

B Details of phase-field modeling for brittle fracture

For the case of a hyperelastic material, the stored energy can be obtained by multiplying the degradation function with the Helmholtz free energy density of the intact material $\psi _0$. This leads to the stored energy of the body given by

$$\begin{aligned} E (\varvec{F},d) = \int _{\Omega } g(d)\psi _0(\varvec{F}) \,{\text {d}}V, \end{aligned}$$

(68)

Moreover, using the chain rule and the time derivative of (68), the rate of degraded stored energy reads

$$\begin{aligned} \mathscr {E}(\varvec{F},d, \dot{\varvec{F}}, \dot{d}) =\int _{\Omega } \left( \varvec{P}:\dot{\varvec{F}} -f^f \,\dot{d}\right) \,{\text {d}}V, \end{aligned}$$

(69)

where the first Piola–Kirchhoff stress tensor is given by $\varvec{P} = g(d)\varvec{P}_0$. The fracture driving force $f^f$ being the work conjugate to the phase-field parameter [22, 47] reads

$$\begin{aligned} f^f:=-\dfrac{\partial \psi }{\partial d} =-g^{\prime }(d)\psi _0. \end{aligned}$$

(70)

In addition, considering the long and short-range forces acting on the body, the external power is obtained as

$$\begin{aligned} \mathscr {P} = \int _{\Omega } \rho _0\, \varvec{b} \cdot \dot{\varvec{\varphi }} \,{\text {d}}V + \int _{\partial _{\varvec{N}}\Omega }\varvec{T}\cdot \dot{\varvec{\varphi }} \,{\text {d}}A. \end{aligned}$$

(71)

Furthermore, the regularized crack surface energy is defined using the critical fracture energy constant $g_c$. This energy is required to convert a fully intact matter into a fully cracked one. Considering constant $g_c$ and the crack surface density functional, the crack energy is obtained

$$\begin{aligned} D (d) = \int _{\Omega } g_c\gamma (d, \nabla d)\, {\text {d}}V. \end{aligned}$$

(72)

Therefore, using the time derivative and the chain rule, the dissipation functional $\mathscr {D}$ for elastic materials reads

$$\begin{aligned} \mathscr {D}(\dot{d}) = \int _{\Omega } g_c\dot{\gamma }\, {\text {d}}V = \int _{\Omega } g_c \frac{\partial \gamma }{\partial d}\dot{d} \,{\text {d}}V, \end{aligned}$$

(73)

where according to the second law of thermodynamics, only the non-negative values of the dissipation functional are admissible, i.e., $\mathscr {D}(\dot{d})\ge 0$. Furthermore, Miehe et al. [47] enforced an irreversibility condition locally by defining a ramp-type energy function, that explodes for a negative evolution of the phase-field parameter. That means, the phase-field parameter is not allowed to reduce at any material point. Furthermore, the derivative of the crack surface density functional with respect to the phase-field parameter reads

$$\begin{aligned} \frac{\partial \gamma }{\partial d} = \frac{1}{l_s} \left( d - l_f^2 \Delta d\right) , \end{aligned}$$

(74)

where $\Delta d$ is the material Laplacian of the phase-field parameter. At this stage, using the rate of stored energy functional (69), the external power (71) and the dissipation functional (73), the balance of mechanical power can be described as

$$\begin{aligned} \Pi (\dot{\varvec{\varphi }}, \dot{d}):= \mathscr {E}(\dot{\varvec{\varphi }}, \dot{d}) + \mathscr {D}(\dot{d}) - \mathscr {P}(\dot{\varvec{\varphi }}). \end{aligned}$$

(75)

Hence, the rates of deformation and damage parameter can be obtained from the variational principle

$$\begin{aligned} (\dot{\varvec{\varphi }}, \dot{d}) = \arg \{\inf _{\dot{\varvec{\varphi }}\in \mathscr {W}_{\varvec{\varphi }}}\,\inf _{\dot{d}\in \mathscr {W}_{d}}\,\Pi (\dot{\varvec{\varphi }}, \dot{d}) \}, \end{aligned}$$

(76)

where the following Dirichlet-type boundary condition is satisfied for the state variables

$$\begin{aligned} \varvec{\dot{\varphi }} \in \mathscr {W}_{\varvec{\varphi }}:=\Biggl \{\varvec{\dot{\varphi }}| \varvec{\dot{\varphi }}=\varvec{0}\, \,\textrm{on}\, \partial _{\varvec{N}} \Omega _{\varvec{\varphi }}\Bigg \} \,\,\,\,\,\textrm{and} \,\,\,\,\, \dot{d}\in \mathscr {W}_{d}:=\Biggl \{\dot{d}| \dot{d}=0\, \,\textrm{on}\, \partial _{\varvec{N}} \Omega _{d}\Bigg \}. \end{aligned}$$

(77)

The variation of the balance of mechanical power leads to two equations: the balance of linear momentum

$$\begin{aligned} \text {Div}\, \varvec{P} +\rho _0\, \varvec{b} = \varvec{0}. \end{aligned}$$

(78)

and the Kuhn–Tucker complementary conditions $\dot{d} \ge 0$, $f^f - g_c \delta _d \gamma \le 0$, and $(f^f - g_c \delta _d \gamma )\dot{d}=0$. That means the crack does not propagate as soon as $f^f - g_c \delta _d \gamma < 0$. On the other hand, as soon as the driving force reaches the critical value $f^f = g_c \delta _d \gamma$, the crack propagates. Remembering $f^f = -g^{\prime }(d)\psi _0$, i.e., (70), the case distinction can be formulated as

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{d} = 0 \,\,\,\,\,\,\, \text {if} \,\,\,\,\,\,\, -g^{\prime }(d)\psi _0 - g_c \delta _d \gamma < 0, \\ \dot{d} > 0 \,\,\,\,\,\,\, \text {if} \,\,\,\,\,\,\, -g^{\prime }(d)\psi _0 = g_c \delta _d \gamma . \end{array}\right. } \end{aligned}$$

(79)

Due to the local irreversibility, this type of formulation may lead to an unrealistic evolution of the phase-field parameter far away from the localization area in the case of ductile fracture. Moreover, May et al. [45] have shown numerically that discretized forms of such type of formulation are not necessarily $\Gamma$-convergent. This is also rooted in the local irreversibility condition, see [42]. This flaw was to some extent removed in Miehe et al. [50] by considering a fracture energy threshold. This threshold is imposed by a material constant $w_c$. That means, the phase-field parameter would evolve just in the material points, which possess energies beyond this threshold. Furthermore, the local irreversibility constraint is only enforced within this region which will vary throughout the crack formation in the general case. Taking this threshold into account, the thermodynamic force in (70) is modified to

$$\begin{aligned} f^f:=-g^{\prime }(d)(\psi _0 - w_c). \end{aligned}$$

(80)

This formulation is sensitive to the definition of the material parameter $w_c$. In other words, the phase-field parameter would only evolve in the region, primarily determined by this material parameter. Not only this, but also the damage parameter is not allowed to get more localized within this region. Additionally, the discrete form of this formulation is not proven to be $\Gamma$-convergent.