Article

# Phase-Field Model for the Simulation of Brittle-Anisotropic and Ductile Crack Propagation in Composite Materials

Christoph Herrmann 1,*, Daniel Schneider 1,2 ID, Ephraim Schoof 1 ID, Felix Schwab 3 ID and Britta Nestler 1,2

<figure><img src="image_1.png" /></figure>

check for
updates

**Citation:** Herrmann, C.; Schneider, D.; Schoof, E.; Schwab, F.; Nestler, B. Phase-Field Model for the Simulation of Brittle-Anisotropic and Ductile Crack Propagation in Composite Materials. *Materials* **2021**, *14*, 4956. https://doi.org/10.3390/ma14174956

Academic Editor: Luisa Silva

Received: 26 July 2021
Accepted: 27 August 2021
Published: 30 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_2.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** In this work, a small-strain phase-field model is presented, which is able to predict crack propagation in systems with anisotropic brittle and ductile constituents. To model the anisotropic brittle crack propagation, an anisotropic critical energy release rate is used. The brittle constituents behave linear-elastically in a transversely isotropic manner. Ductile crack growth is realised by a special crack degradation function, depending on the accumulated plastic strain, which is calculated by following the J₂-plasticity theory. The mechanical jump conditions are applied in solid-solid phase transition regions. The influence of the relevant model parameters on a crack propagating through a planar brittle-ductile interface, and furthermore a crack developing in a domain with a single anisotropic brittle ellipsoid, embedded in a ductile matrix, is investigated. We demonstrate that important properties concerning the mechanical behaviour of grey cast iron, such as the favoured growth of cracks along the graphite lamellae and the tension–compression load asymmetry of the stress–strain response, are covered by the model. The behaviour is analysed on the basis of a simulation domain consisting of three differently oriented elliptical inclusions, embedded in a ductile matrix, which is subjected to tensile and compressive load. The material parameters used correspond to graphite lamellae and pearlite.

**Keywords:** phase-field; multiphase-field; grey cast iron; brittle fracture; ductile fracture; anisotropic fracture

## 1. Introduction

For the numerical prediction or analysis of fractures and the associated failure of components, the investigation of the formation and propagation of cracks is necessary. In the study of crack development, some physical models and numerical methods have gained acceptance over time. Important approaches of the “old guard”, for example, are discrete formulations, with cohesion zone models [1] as the most prominent representative, and continuum damage models [2,3]. With these approaches, it is, however, difficult to capture phenomena such as crack branching and merging as well as crack nucleation processes in arbitrarily complex systems. These phenomena can be intrinsically captured when using a phase-field approach to crack modelling, which is why such methods are enjoying growing popularity [4–6]. The evolution of cracks is described within the phase-field framework by a regularised representation of the Griffith criterion [7], where the degree of damage is given by a scalar continuous phase field. Subsequently, a phase-field model is then presented that allows the simulation of crack development processes in composite materials, consisting of brittle-anisotropic and ductile components. This is done by analogy with lamellar graphite cast iron, in which the matrix of the cast iron behaves

<footnote n="1">¹ Institute of Applied Materials (IAM-CMS), Karlsruhe Institute of Technology (KIT), Kaiserstrasse 12, D-76131 Karlsruhe, Germany; daniel.schneider@kit.edu (D.S.); ephraim.schoof@kit.edu (E.S.); britta.nestler@kit.edu (B.N.)</footnote>

<footnote n="2">² Institute of Digital Materials Science (IDM), Karlsruhe University of Applied Sciences, Moltkestrasse 30, D-76133 Karlsruhe, Germany</footnote>

<footnote n="3">³ Groupe de Physique des Matériaux, University of Rouen, CNRS, Avenue de l'Université, F-76800 Saint Etienne du Rouvray, France; felix.schwab@univ-rouen.fr</footnote>

<footnote>* Correspondence: christoph.herrmann2@kit.edu</footnote>
in a ductile manner and the inclusions in a brittle-anisotropic manner. Subsequently, the constituents of cast iron with lamellar graphite, the pearlitic matrix, and the graphite particles are represented in an idealised form. The ductile phase is considered to be purely pearlitic and is modelled as a homogeneous ductile isotropic phase, according to the $J_2$-plasticity theory with linear hardening. The effects of the microscopic structure of pearlite, consisting of ferrite and cementite, are only considered in a homogenised manner. Brittle phases are related to lamellar graphite. As inclusions, they are similar to graphite lamellae in terms of shape and mechanical properties. However, they are larger. On the atomistic level, lamellar graphite is made up of interconnected graphene layers. Within the graphene layers, covalent bonds exist between the carbon atoms. The different graphene layers are connected by weak van der Waals forces [8–10]. Due to its atomistic structure, lamellar graphite possesses a transverse isotropy, with regard to its elastic stiffness and critical energy release rate, and so do the anisotropic brittle phases. For multiphase materials, the proposed model is the first phase-field crack model we are aware of, which combines brittle-anisotropic and ductile crack propagation. For this purpose, two existing phase-field models for crack development are combined and introduced into the multiphase context: one for brittle-anisotropic crack development and one for ductile crack development.

The formulation of the previous publication by Prajapati et al. [11] serves as a model for the description of brittle-anisotropic crack development. Here, the critical energy release rate is directionally dependent. Alternative approaches for anisotropic crack growth are based on multiplying the divergence term of the evolution equation of the crack phase field with an anisotropy tensor, resulting in a directional dependency (see, e.g., [12–14]). In the context of gradient energy, however, the physical representation of such an anisotropy tensor is still unclear [14].

When choosing a possible ductile crack model, various approaches are available, which are briefly summarised below. The model of Duda et al. [15] was one of the first models to integrate plasticity into a phase-field crack model, so as to simulate brittle fractures in plastic materials. Some important phenomenological characteristics of ductile fracture, reported in experimental literature, could be reproduced in Ambati et al. [16] by coupling the degradation function with the plastic strain state. In Ambati et al. [17], the model was extended to finite strains. In the context of thermoplasticity, a thermodynamically consistent phase-field crack model for brittle to ductile fractures is introduced in Miehe et al. [18] at large deformations. As an extension to this, they also worked on porous-isotropic plasticity in Miehe et al. [19]. In Kuhn et al. [20], an elastoplastic phase-field fracture model, where a monolithic solution is possible, is proposed. In Miehe et al. [21], gradient plasticity is used at finite strains in order to model ductile fracture in a variational-based phase-field framework. In a recently published book article Alessi et al. [6], prominent phase-field models for ductile crack growth are compared, and a study of their predictive capabilities is conducted.

Based on the investigations of the anisotropic brittle crack model of our previously published work, Prajapati et al. [11], the focus and novelty of this work is on coupling the anisotropic brittle crack model with the ductile crack model of Ambati et al. [16]. As modifications of Prajapati et al. [11], plastic deformations are additionally considered and the degradation function of Ambati et al. [16] is introduced to map ductile fracture characteristics. The calculation of the elastic and plastic fields is based on Herrmann et al. [22], and thus the mechanical jump conditions are applied to obtain locally homogenised stresses and plastic strains.

In Section 2, the model extension of the work presented by Prajapati et al. [11] is introduced in order to allow the propagation of cracks in anisotropic brittle and ductile materials. Section 3 discusses all relevant numerical aspects, while Section 4 includes numerical studies and applications. In Section 4.1, a detailed investigation of a crack is given, which passes a planar diffuse interface between an anisotropic brittle and a ductile phase. In Section 4.2, crack nucleation and propagation in a domain with a single anisotropic ellipsoid, embedded in a ductile matrix, is analysed. To show that important
properties concerning the mechanical behaviour of grey cast iron can be predicted by the model, a simulation area, consisting of three elliptical, anisotropic inclusions, which are oriented differently in a ductile matrix, is subjected to tensile and compressive load in Section 4.3. Finally, Section 5 summarises the presented results and gives an outlook on upcoming applications.

## 2. Model Formulation

The model presented in Prajapati et al. [11] provides the basis for the present work. In order to allow simulations of crack nucleation and propagation in materials with anisotropic brittle and ductile constituents, the model is extended by an additional degradation function, which is sensitive to the accumulated plastic strain, as proposed by Ambati et al. [16]. The reasons for choosing this model from the numerous published ductile crack models are the following: In Ambati et al. [16], it is shown that important characteristics concerning ductile crack nucleation and propagation can be mapped by the model. This is made possible by a modified degradation function, which includes the accumulated plastic strain. The modification leads to crack initiation processes, which are consistent with experimental observations, despite the assumption of small deformations. Most of the features under consideration, e.g., the strain localisation in the middle of an I-shaped specimen, under tensile load, cannot be mapped by the majority of the other published ductile phase-field crack models. Compared to brittle phase-field crack models, the main difference in most other models that account for plastic deformations is that a plastic energy contribution, due to hardening, is introduced in the free energy functional, and that the stress calculation takes the plastic strains into account by splitting the total strain into an elastic and a plastic part [6]. These changes lead to stress–strain curves corresponding to a ductile material, but the plastic contribution in the energy functional, especially in the case of small deformations, is substantially smaller than the elastic strain energy and thus has only little influence on the crack nucleation and propagation. In Ambati et al. [16], the formulation of the ductile degradation function furthermore ensures that the experimentally determined yield points are represented accurately. If interested in alternative approaches of modelling ductile crack propagation using the phase-field method, the reader is referred to Alessi et al. [6]. In this work, popular phase-field models of fracture coupled with plasticity are compared by means of the resulting predictive capabilities for several well-defined problems.

In the classical understanding, phase-field models in which more than two phase fields can coexist are called multiphase-field models, which usually implies that each phase field possesses its own evolution equation. The existence of several evolution equations, however, only applies to the present model to a limited extent, since no solid-solid phase transformation is integrated. For this reason, the present model is not a full “multiphase model” but nevertheless builds on a multiphase-field representation of physically distinguishable regions and their interfaces. However, since a solid-solid phase evolution is to be integrated in further work and can be implemented straightforwardly, starting from the present model formulation, by applying a staggered scheme for the crack phase field and the solid-solid evolution, the solid-solid phase transitions were nevertheless discretised diffusely.

Within phase-field models, the transition between physically distinguishable solid phases is regularised, in order to allow an easy tracking of grain boundaries and to enable phase transformation processes. For this purpose, solid phase fields, $\phi^\alpha(x, t)$, $\alpha = 1...N$, are introduced, with $N$ as the number of solid phases. Within phase transition regions, the phase fields change continuously in the range of $0 < \phi^\alpha < 1$, while in material points, where only one phase $\alpha$ exists, $\phi^\alpha = 1$ applies. To simplify the mathematical representation in the following, all solid phase fields are collected in the $N$-tuple $\hat{\phi}(x, t)$, and their gradients are joined in $\nabla \hat{\phi}(x, t)$.

Within phase-field crack models in a domain $\Omega \subset \mathbb{R}^d, d \in \{1, 2, 3\}$, the critical energy release rate is modelled in a regularised manner by a diffuse transition between solid
phases and cracks, whereby on top of the phase fields, representing the solid phases, an additional crack phase field, $\phi_c(x, t)$ : $\Omega \times \mathbb{R}_{\ge 0} \to [0, 1]$, is introduced in the spatial point $x$, at the time $t$. The crack phase field possesses a continuous transition between intact, $\phi_c = 0$, and fractured, $\phi_c = 1$, material points.

Each phase field represents the corresponding volume fraction of the phase, so that the constraint

$$ \phi_c + \sum_{\alpha=1}^{N} \phi^{\alpha} = 1 \qquad (1) $$

has to be fulfilled, which applies to the solid phase fields as well as the crack phase field. The understanding that each phase field reflects the volume fraction of a particular phase in a material point leads to the fact that during the growth of the crack phase field, i.e., the increase of the volume fraction of the crack phase, the corresponding volume fraction of the solid phases dissipates, which is comparable to continuum damage models. In terms of mass conservation, this is only acceptable as long as a negligible amount of the total mass of a system dissipates. If solid phases are converted into the crack phase, the relative composition of the solid phases is maintained with respect to the absolute volume fraction of the solid phases. Therefore, the transformation of the solid phase fields is given by

$$ \dot{\phi}^{\alpha}(x, t) = -h_s^{\alpha} \dot{\phi}_c. \qquad (2) $$

Thus, it is sufficient to only solve the phase-field equation of the crack phase explicitly, into which the mechanical forces of a multiphase field problem are incorporated. To identify the relative volume fraction of a solid phase in relation to the absolute fraction of solid phases, it is required to formulate an interpolation function for solid phases:

$$ h_s^{\alpha}(\dot{\boldsymbol{\phi}}) = \frac{\phi^{\alpha}}{\sum_{\beta=1}^{N} \phi^{\beta}}, \qquad (3) $$

which satisfies the condition $\sum_{\alpha=1}^{N} h_s^{\alpha} = 1$.

According to Griffith [7], in the case of linear fracture mechanics, crack growth occurs when the energy release rate of a material point, which corresponds to the elastic strain energy, exceeds the necessary free energy of an emerging free surface, given by a critical energy release rate $G_c(x, t)$. For the body $\Omega$, this energetic view is captured in the phase-field context in a regularised way by using a free energy functional of the form

$$ \mathcal{F}(\dot{\boldsymbol{\phi}}, \phi_c, \nabla \phi_c, \hat{\boldsymbol{\varepsilon}}_e, \boldsymbol{\varepsilon}_{\text{acc}}) = \int_{\Omega} \underbrace{\frac{3}{8} \bar{G}_c \left( \frac{\phi_c}{\bar{l}} + \bar{l} |\nabla \phi_c|^2 \right)}_{=f(\dot{\boldsymbol{\phi}}, \phi_c, \nabla \phi_c, \hat{\boldsymbol{\varepsilon}}_e, \boldsymbol{\varepsilon}_{\text{acc}})} + f_e + f_p \, dV, \qquad (4) $$

which is comparable with [6,23], for example. The free energy functional includes the interpolated critical energy release rate, $\bar{G}_c$, the effective elastic strain energy,

$$ f_e(\hat{\boldsymbol{\varepsilon}}_e, \boldsymbol{\varepsilon}_{\text{acc}}, \dot{\boldsymbol{\phi}}) = \sum_{\alpha=1}^{N} h_s^{\alpha} f_e^{\alpha}(\boldsymbol{\varepsilon}_e^{\alpha}, \boldsymbol{\varepsilon}_{\text{acc}}^{\alpha}, \phi_c), \qquad (5) $$

and the effective plastic energy contribution, due to hardening,

$$ f_p(\boldsymbol{\varepsilon}_{\text{acc}}, \dot{\boldsymbol{\phi}}) = \sum_{\alpha=1}^{N} h_s^{\alpha} f_p^{\alpha}(\boldsymbol{\varepsilon}_{\text{acc}}^{\alpha}), \qquad (6) $$
by applying the so-called AT-1 representation, according to [24], with regard to the geometrical shape of the diffuse transition between the solid phases and the crack. As a typing aid, all phase-inherent elastic tensors, $\epsilon_e^\alpha$, are collected in $\epsilon_e$, and the phase-inherent accumulated plastic strains, $\epsilon_{acc}^\alpha$, are contained in $\epsilon_{acc}$. Phase-inherent quantities are quantities associated with each solid phase, i.e., each solid phase has its own mechanical fields. The width of the regularised transition between damaged and undamaged material points is determined by the interpolated regularisation parameter:

$$ \bar{l}(\hat{l}, \hat{\phi}) = \sum_{\alpha=1}^{N} h_s^{\alpha} l^{\alpha}. \qquad (7) $$

The regularisation parameter is considered as a material property and is chosen as proposed by Tanné et al. [25]. According to Tanné et al., the choice of the regularisation parameter guarantees that crack nucleation only takes place when the stresses in an isotropic brittle material reach the tensile strength. The calculation procedure for the determination of the regularisation parameter, shown in Tanné et al., is not readily applicable to the model used in this paper, which has an anisotropic and ductile material behaviour. Nevertheless, the correlations from Tanné et al. are used to calculate the regularisation parameters of each phase, $l^\alpha$, which are collected in $\hat{l}$, in order to use the best possible trade-off available. It is particularly important to mention that the regularisation parameters of transversely isotropic phases are related to the critical energy release rate and Young's modulus of the weakest mechanical material direction.

It is known that the tension-compression asymmetry [26] of cast iron with lamellar graphite is controlled by the graphite lamellae [9,27]. While under tensile load, the graphite lamellae have a negligible loading capacity and act as microcracks. They tend to close under compressive forces and transmit the applied load [9,28]. To map this property, an unphysical crack propagation under compressive load must be avoided. Therefore, only the positive part of the phase-inherent elastic strain energy should lead to nucleation and growth of the crack phase. The elastic strain energy is decomposed into its positive and negative energy part:

$$ f_e^\alpha(\epsilon_e^\alpha, \epsilon_{acc}^\alpha, \phi_c) = h_c^\alpha[f_e^\alpha]^+ + [f_e^\alpha]^-, \qquad (8) $$

using the spectral decomposition for isotropic phases, as proposed by Miehe et al. [4], and following Teichtmeister et al. [12], so as to allow a decomposition of transversely isotropic materials. Here, the crack degradation function $h_c^\alpha$ reflects the dissipation of the solid phases on the energetic level, caused by crack growth, and is applied exclusively to the positive energy contribution. To allow the transition of cracks between brittle and ductile phases, two different crack degradation functions are used:

$$ h_c^\alpha(\epsilon_{acc}^\alpha, \phi_c) = \begin{cases} (1 - \phi_c)^2, & \text{brittle} \\ (1 - \phi_c)^2 \frac{\epsilon_{acc}^\alpha}{\epsilon_{acc,crit}^\alpha}, & \text{ductile.} \end{cases} \qquad (9) $$

The brittle degradation function is very common for brittle phase-field crack models. The ductile degradation function was introduced by Ambati et al. [16]. By incorporating the accumulated plastic strain into the ductile crack degradation function, the prediction of crack nucleation is possible, despite the small-strain formulation, which is in line with experimental findings [16]. Furthermore, it guarantees that the phase-inherent stresses of a ductile phase are not degraded, as long as no plastic deformation has occurred in a material point, even though the crack phase does exist. This ensures that plastic deformations only arise when the phase-inherent stress fulfils the yield condition. To control the ductile crack propagation, the threshold value $\epsilon_{acc,crit} \in \mathbb{R}_{\ge 0}$ [16] is applied.

The local homogenisation approach of Herrmann et al. [22] is used as an elastoplastic multiphase-field model. Thus, the mechanical jump conditions of solid-solid phase
transition regions are applied, which makes an interpolation of mechanical material parameters obsolete, apart from the critical energy release rate. The clean separation between elastoplastic and purely elastic phases is only possible without any difficulties, when this homogenisation approach is used. The homogenisation approach enables the calculation of phase-inherent stresses, $\sigma^\alpha$, and plastic strains, $\epsilon_p^\alpha$, and therefore avoids the occurrence of plastic strains within brittle phases. Another advantage of the elastoplastic homogenisation approach is that an anisotropic elastic behaviour can be realised in one phase, even though an isotropic nonlinear constitutive model is applied in another phase. This holds especially in multiphase regions, which, in contrast, is not possible with models using interpolated material parameters. In multiphase regions, this is particularly required for the considered composite material, as the volumetrically interpolated stiffnesses of the isotropic pearlite and the transversely isotropic ellipsoids would not allow the $J_2$-plasticity model to be used in combination with a return mapping algorithm.

The following considerations are limited to small deformations; the total strain is expressed as $\epsilon(H) = \frac{1}{2}(H + H^T)$, using the displacement gradient $H = \nabla u(x, t)$; and an additive decomposition of the strains $\epsilon = \epsilon_e + \epsilon_p$ is valid. The total, elastic, and plastic strains are identified as $\epsilon(x, t)$, $\epsilon_e(x, t)$, and $\epsilon_p(x, t)$. Within the simulations, sharp crack tips, caused by the brittle ellipsoids, act as stress concentrators. Thus, the local strains of >2% even occur when the macroscopic strain is small. Under compressive load, grey cast iron specimens fail at a macroscopic strain of >8%. Under the assumption of small deformations, the geometrical nonlinearities, which actually occur under the described states, are neglected for the sake of simplicity.

In material points with multiple solid phases, the mechanical jump conditions, given by the sharp interface theory of bounded solid-solid transitions [29], are used. On the one hand, this is the static balance of the linear momentum on a singular surface:

$$ (\sigma^\alpha - \sigma^\beta) \cdot n^{\alpha\beta} = 0, \qquad (10) $$

stating that the jump of the stresses of two phases $\alpha$ and $\beta$ vanishes in the normal direction of a singular surface between the solid phases, $n^{\alpha\beta}(\nabla\hat{\phi})$. On the other hand, the Hadamard jump condition

$$ (H^\alpha - H^\beta) = a^{\alpha\beta} \otimes n^{\alpha\beta} \qquad (11) $$

represents a no-slip condition in the tangential directions of singular surfaces. The rank-1 tensor $a^{\alpha\beta}$ defines the amplitude of the jump of the displacement gradient in the direction of the surface normal vector, as prescribed by the Hadamard lemma. In regions where only two phases coexist, the unit normal vectors between the two phases can be written as $n^{\alpha\beta} = \nabla\phi^\alpha / |\nabla\phi^\alpha| = -\nabla\phi^\beta / |\nabla\phi^\beta|$, according to [30]. In multiphase regions, the normal vector between two solid phases reads as $n^{\alpha\beta} = (\nabla\phi^\alpha - \nabla\phi^\beta) / |\nabla\phi^\alpha - \nabla\phi^\beta|$, as described in [31]. Applying the jump conditions with respect to the phase with the largest volume fraction, e.g., phase 1, in regions with $N$ coexisting solid phases, $N-1$ unknown rank-1 tensors, $a^{12}, a^{13}, \dots, a^{1N}$, have to be determined to describe the jumps of the phase-inherent displacement gradients at the $\alpha-\beta$ transition in the $n^{\alpha\beta}$-direction. For the sake of simplicity, all rank-1 tensors are collected in the tuple

$$ \hat{a} = (a^{12}, a^{13}, \dots, a^{1N})^T. \qquad (12) $$

The effective displacement gradient of multiphase material points, $H^{\text{eff}}$, is defined by all phase-inherent displacement gradients, collected in $\hat{H}$, by means of a volumetric decomposition:

$$ H^{\text{eff}}(\hat{H}, \hat{\phi}) = \sum_{\alpha=1}^{N} h_s^\alpha H^\alpha, \qquad (13) $$
where $H^\alpha$ are the phase-dependent displacement gradients. By rearranging the effective displacement gradient in combination with Equation (11), the displacement gradient of phase 1 can be written as

$$H^1 = H^{\text{eff}} + \sum_{\alpha=2}^{N} h_s^\alpha a^{1\alpha} \otimes n^{1\alpha}. \quad (14)$$

For the remaining deformation gradients, this results in

$$H^\alpha = H^1 - a^{1\alpha} \otimes n^{1\alpha}. \quad (15)$$

All unknown rank-1 tensors, $a^{1\alpha}$, are calculated by solving the equation system

$$\hat{g}(\hat{a}) = \begin{pmatrix} (\sigma^1 - \sigma^2) \cdot n^{12} \\ (\sigma^1 - \sigma^3) \cdot n^{13} \\ \vdots \\ (\sigma^1 - \sigma^{\tilde{N}}) \cdot n^{1\tilde{N}} \end{pmatrix} = \begin{pmatrix} (h_c^1 \mathbb{C}^1 [\varepsilon_e^1(\hat{a})] - h_c^2 \mathbb{C}^2 [\varepsilon_e^2(\hat{a})]) \cdot n^{12} \\ (h_c^1 \mathbb{C}^1 [\varepsilon_e^1(\hat{a})] - h_c^3 \mathbb{C}^3 [\varepsilon_e^3(\hat{a})]) \cdot n^{13} \\ \vdots \\ (h_c^1 \mathbb{C}^1 [\varepsilon_e^1(\hat{a})] - h_c^{\tilde{N}} \mathbb{C}^{\tilde{N}} [\varepsilon_e^{\tilde{N}}(\hat{a})]) \cdot n^{1\tilde{N}} \end{pmatrix} = 0, \quad (16)$$

composed of the pairwise static momentum balances on singular interfaces, Equation (10). The phase-inherent elastic strains are given by $\varepsilon_e^\alpha = \varepsilon^\alpha - \varepsilon_p^\alpha$, with $\varepsilon^\alpha = \frac{1}{2}(H^\alpha + (H^\alpha)^T)$ taking Equations (14) and (15) into account. For the phase-inherent stresses, this results in

$$\sigma^\alpha = h_c^\alpha \mathbb{C}^\alpha [\varepsilon_e^\alpha], \quad (17)$$

where in the elastic stiffness tensor of phase $\alpha$ is denoted by $\mathbb{C}^\alpha$. It can be observed that these stresses cannot be derived from the elastic strain energy, Equation (8), in the sense of a hyperelastic material model, since, unlike the elastic strain energy, the entire stress tensor is degraded, and not only its positive part. A similar approach to formulate a phase-field crack model was introduced in Ambati et al. [5], where this is referred to as a hybrid model. A $J_2$-plasticity model [32] with linear isotropic hardening, resulting in the yield function

$$f_y^\alpha(s^\alpha, \varepsilon_{\text{acc}}^\alpha) = \sqrt{\frac{3}{2}} |s^\alpha| - (\sigma_{y,0}^\alpha + H^\alpha \varepsilon_{\text{acc}}^\alpha) \le 0, \quad (18)$$

is applied to compute the phase-inherent plastic strains $\varepsilon_p^\alpha$. Here, $s^\alpha = \sigma^\alpha - 1/3 \text{tr}(\sigma^\alpha)I$ represents the deviatoric stress tensor, and $|s^\alpha| = \sqrt{s_{ij}^\alpha s_{ij}^\alpha}$ its Frobenius norm, where the Einstein summation convention is valid. The first term of the yield function is the equivalent von Mises stress, which is compared with the yield stress, expressed by the yield strength $\sigma_{y,0}^\alpha$ and the strain hardening $H^\alpha \varepsilon_{\text{acc}}^\alpha$, which includes the phase-inherent linear hardening parameter $H^\alpha$. If the yield function is violated, plastic yielding occurs. Using the associative flow rule

$$\dot{\varepsilon}_p^\alpha = \dot{\gamma}^\alpha \frac{\partial f_y^\alpha(s^\alpha, \varepsilon_{\text{acc}}^\alpha)}{\partial s^\alpha}, \quad (19)$$

the evolution of the plastic strain is described. The phase-inherent accumulated plastic strain $\varepsilon_{\text{acc}}^\alpha = \int_t \sqrt{2/3} |\dot{\varepsilon}_p^\alpha| dt$ serves as an internal variable. As a result, the internal plastic energy contribution of phase $\alpha$ [32] reads as

$$f_p^\alpha(\varepsilon_{\text{acc}}^\alpha) = \frac{1}{2} H^\alpha \varepsilon_{\text{acc}}^\alpha. \quad (20)$$
To establish the consistency condition $f_y^\alpha \stackrel{!}{=} 0$ in the case of plastic yielding, a two-step return mapping algorithm is employed. Applying a volumetric interpolation, the effective accumulated plastic strain

$$ \epsilon_{\text{acc}}^{\text{eff}}(\epsilon_{\text{acc}}, \hat{\boldsymbol{\phi}}) = \sum_{\alpha=1}^{N} h_{\text{s}}^{\alpha} \epsilon_{\text{acc}}^{\alpha} \qquad (21) $$

is introduced.

To solve the static balance of linear momentum,

$$ \operatorname{div} \sigma^{\text{eff}} = 0, \qquad (22) $$

for the displacement field $u(x, t)$, the degraded volume-averaged effective stresses

$$ \sigma^{\text{eff}}(\hat{\sigma}, \hat{\boldsymbol{\phi}}) = \sum_{\alpha=1}^{N} h_{\text{s}}^{\alpha} \sigma^{\alpha} \qquad (23) $$

are used, which result from all phase-inherent stresses, collected in $\hat{\sigma}$.

An important part of a phase-field crack model with diffusely overlapping transition regions of solid phases is the volumetrically interpolated critical energy release rate

$$ \bar{G}_{\text{c}}(\hat{\boldsymbol{\phi}}, \nabla \phi_{\text{c}}) = \sum_{\alpha=1}^{N} h_{\text{s}}^{\alpha} G_{\text{c}}^{\alpha}, \qquad (24) $$

allowing the transition between phases with different crack characteristics and anisotropy orientations. A transverse anisotropy resulting for, e.g., lamellar graphite from its atomistic structure, as indicated in the introduction, can be modelled for phase $\alpha$, with an anisotropic critical energy release rate, such as

$$ G_{\text{c}}^{\alpha}(\nabla \phi_{\text{c}}) = G_{\text{c},0}^{\alpha} \left( (\tilde{n}_{x}^{\alpha})^2 + F_{\text{aniso}}^{\alpha} \left( (\tilde{n}_{y}^{\alpha})^2 + (\tilde{n}_{z}^{\alpha})^2 \right) \right). \qquad (25) $$

Through the phase-inherent anisotropy factor $F_{\text{aniso}}^\alpha \in \mathbb{R}_{>1}^\alpha$, the referential critical energy release rate $G_{\text{c},0}^\alpha$ is deformed into an ellipsoid. If $F_{\text{aniso}}^\alpha = 1$ is chosen, the critical energy release rate is isotropic, which is the case for ductile phases. The orientations of the phase-dependent anisotropies are described by the phase-inherent Euler angles and are taken into account by the transformation matrices $Q_{\text{aniso}}^\alpha(\varphi_{\text{aniso}}^\alpha, \theta_{\text{aniso}}^\alpha, \psi_{\text{aniso}}^\alpha)$, which are applied to transform the outward-pointing normal vector of the crack phase field from the reference $x, y, z$-coordinate system into the hexagonal lattice $\tilde{x}, \tilde{y}, \tilde{z}$-coordinate systems of graphite lamellae:

$$ \tilde{n}^{\alpha}(\nabla \phi_{\text{c}}) = \begin{pmatrix} \tilde{n}_{x}^{\alpha}(\nabla \phi_{\text{c}}) \\ \tilde{n}_{y}^{\alpha}(\nabla \phi_{\text{c}}) \\ \tilde{n}_{z}^{\alpha}(\nabla \phi_{\text{c}}) \end{pmatrix} = -Q_{\text{aniso}}^{\alpha} n_{\text{c}} = -Q_{\text{aniso}}^{\alpha} \frac{\nabla \phi_{\text{c}}}{|\nabla \phi_{\text{c}}|}. \qquad (26) $$

Graphite lamellae show a preferred growth direction along their graphene planes (a-direction of their hexagonal crystal system) [33]. Therefore, their spatial expansion within the basal planes is much larger than orthogonal to it [34]. Based on this fact, the orientation of the longitudinal axis of a graphite lamella can be used to define its transformation matrices $Q_{\text{aniso}}^\alpha$. From the formulation of Equation (25), it can be concluded that the crystalline a-direction of the hexagonal crystal system of a graphite lamella points in the same direction as the $\tilde{x}$-axis.
Following Allen-Cahn [35], the temporal and spatial evolution of the crack phase field is given as

$$ \dot{\phi}_c(x, t) = -M \frac{\delta F(\hat{\phi}, \phi_c, \nabla \phi_c, \epsilon_e, \epsilon_{acc})}{\delta \phi_c}, \quad (27) $$

including the variational derivative of the total free energy of the system, with respect to the crack phase field and the kinetic coefficient $M$. To correctly predict the kinetics of crack propagation, the kinetic coefficient has to be chosen according to experimental findings. In this work, the explicit Euler scheme is applied to calculate the time derivatives in order to determine an equilibrium state of the crack phase field for each mechanical load increment, $\dot{\phi}_c \approx 0$. Thus, there is only a fictive time dependence, and the kinetic coefficient can be chosen in such a way that the simulations are numerically stable. The performed variational derivation of Equation (4) results in

$$ \begin{aligned} \dot{\phi}_c &= -M \left( \frac{\partial f}{\partial \phi_c} - \nabla \cdot \frac{\partial f}{\partial \nabla \phi_c} \right) \\ &= -M \left[ \underbrace{\frac{3}{8} \left( \bar{G}_c \left( \frac{1}{\bar{l}} - 2\bar{l}\nabla \cdot \nabla \phi_c \right) - \nabla \cdot \frac{\partial \bar{G}_c}{\partial \nabla \phi_c} \left( \frac{\phi_c}{\bar{l}} + \bar{l}|\nabla \phi_c|^2 \right) \right)}_{g_c(\hat{\phi}, \phi_c, \nabla \phi_c, \epsilon_e, \epsilon_{acc}) \le 0} + \frac{\partial f_e}{\partial \phi_c} \right], \end{aligned} \quad (28) $$

where by

$$ \frac{\partial n_c(\nabla \phi_c)}{\partial \nabla \phi_c} = \frac{I - n_c \otimes n_c}{|\nabla \phi_c|} \quad (29) $$

applies for the derivation of the normal vector of the crack phase field [36], which is necessary for the calculation of $\partial \bar{G}_c / \partial \nabla \phi_c$. It should be considered that all vanishing terms of the variational derivative are already excluded in the equation above. To prevent crack healing, $\dot{\phi}_c \ge 0$ is enforced. Therefore, the presented model is classified as a damage model, as the evolution of the crack phase field is irreversible.

All governing equations of the crack model are summarised in Table 2. The meaning of the symbols is summarized in Table 2.

### Model Restrictions

The starting point of developing this model was to study the evolution of cracks in cast iron with lamellar graphite. However, this is only possible to a limited extent, as the new model cannot be applied to the actual length scale, the micrometer range, necessary to resolve the microstructure of cast iron with lamellar graphite. However, with the presented model, the step towards a model that is applicable on the microscopic length scale is no longer far off. As an adjustment, only one extension is necessary, so that the model formulation is independent of the choice of the parameter to define the width of the diffuse transition between the solid phases and the cracks. This could be done in a similar way as the procedure published by Wu and Nguyen [37], but with the additional difficulty that instead of a brittle-isotropic model formulation, a brittle-anisotropic model formulation, combined with a ductile formulation, would have to be considered.
<table><caption>**Table 1.** Governing equations of the phase-field crack model.</caption><tbody><tr><th colspan="2">Mechanical equilibrium</th></tr><tr><td>Static balance of momentum</td><td>divσ = 0</td></tr><tr><td>Static balance of momentum, at a singular surface</td><td>(σ<sup>α</sup> - σ<sup>β</sup>)n<sup>αβ</sup> = 0</td></tr><tr><td>Hadamard jump condition</td><td>(H<sup>α</sup> - H<sup>β</sup>) = a<sup>αβ</sup> ⊗ n<sup>αβ</sup></td></tr><tr><th colspan="2">Plasticity</th></tr><tr><td>Karush-Kuhn-Tucker system</td><td><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 2px;"><span style="font-size: 2.5em; line-height: 1;">{</span><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 2px;"><span style="font-size: 2em; line-height: 1;">f</span><span style="font-size: 2em; line-height: 1;">y</span><span style="font-size: 2em; line-height: 1;">&le;</span><span style="font-size: 2em; line-height: 1;">0</span></div></div><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 2px;"><span style="font-size: 2em; line-height: 1;">&gamma;</span><span style="font-size: 2em; line-height: 1;">&alpha;</span><span style="font-size: 2em; line-height: 1;">&ge;</span><span style="font-size: 2em; line-height: 1;">0</span></div></td></tr><tr><td>Flow rule</td><td>ε<sup>α</sup><sub>P</sub> = &gamma;<sup>α</sup> &part;f<sup>α</sup><sub>y</sub> / &part;s<sup>α</sup></td></tr><tr><td>Yield function</td><td>f<sup>α</sup><sub>y</sub> = &radic;(3/2 |s<sup>α</sup>| - (σ<sup>α</sup><sub>y,0</sub> + H<sup>α</sup>ε<sup>α</sup><sub>acc</sub>))</td></tr><tr><th colspan="2">Crack phase field evolution</th></tr><tr><td>Karush-Kuhn-Tucker system</td><td><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 2px;"><span style="font-size: 2.5em; line-height: 1;">{</span><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 2px;"><span style="font-size: 2em; line-height: 1;">g</span><span style="font-size: 2em; line-height: 1;">c</span><span style="font-size: 2em; line-height: 1;">&ge;</span><span style="font-size: 2em; line-height: 1;">0</span></div></div><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 2px;"><span style="font-size: 2em; line-height: 1;">&phi;</span><span style="font-size: 2em; line-height: 1;">c</span><span style="font-size: 2em; line-height: 1;">&ge;</span><span style="font-size: 2em; line-height: 1;">0</span></div></td></tr><tr><td>Flow rule</td><td>φ<sub>c</sub> = -Mg<sub>c</sub></td></tr><tr><td>Yield function</td><td>g<sub>c</sub> = (3/8) (G&#773;<sub>c</sub> (1/&#773; - 2I&#773;∇ &sdot; ∇φ<sub>c</sub>) - &nabla; &sdot; (∂G&#773;<sub>c</sub>/∂&nabla;φ<sub>c</sub>) (φ<sub>c</sub>/&#773; + &#773;|∇φ<sub>c</sub>|<sup>2</sup>)) + ∂f<sub>e</sub>/∂φ<sub>c</sub></td></tr><tr><th colspan="2">Solid phase evolution</th></tr><tr><td>Flow rule</td><td>φ<sup>α</sup> = -h<sup>α</sup><sub>S</sub>φ<sub>c</sub></td></tr><tr><th colspan="2">Volume fraction conditions of multiphase field approach</th></tr><tr><td>Sum condition of multiphase-field model</td><td>h<sup>α</sup><sub>S</sub> = φ<sup>α</sup> / (N &sum;<sub>β=1</sub><sup>N</sup> φ<sup>β</sup>)</td></tr><tr><td>Sum condition of solid phases</td><td>N &sum;<sub>α=1</sub><sup>N</sup> h<sup>α</sup><sub>S</sub> = 1</td></tr></tbody></table>

<table><caption>**Table 2.** Summary of the most important symbols, categorically grouped.</caption><thead><tr><th>Symbol</th><th>Meaning</th></tr></thead><tbody><tr><td>φ<sup>α</sup></td><td>Solid phase fields</td></tr><tr><td>φ&#775;</td><td>Phase field-N-tuple</td></tr><tr><td>h<sup>α</sup><sub>S</sub></td><td>Interpolation function for solid phases</td></tr><tr><td>F</td><td>Free energy functional</td></tr><tr><td>M</td><td>Mobility kinetic coefficient</td></tr><tr><td>φ<sub>c</sub></td><td>Crack phase field</td></tr><tr><td>h<sup>α</sup><sub>c</sub></td><td>Crack degradation function</td></tr><tr><td>l<sup>α</sup></td><td>Phase-inherent regularisation parameter</td></tr><tr><td>l&#773;</td><td>Interpolated regularisation parameter</td></tr><tr><td>l&#773;</td><td>Collection of phase-inherent regularisation parameters</td></tr><tr><td>G<sup>α</sup><sub>c</sub></td><td>Phase-inherent critical energy release rate</td></tr><tr><td>G&#773;<sub>c</sub></td><td>Interpolated critical energy release rate</td></tr><tr><td>G<sup>α</sup><sub>c,0</sub></td><td>Phase-inherent referential critical energy release rate</td></tr><tr><td>F<sup>α</sup><sub>aniso</sub></td><td>Phase-inherent anisotropy factor</td></tr><tr><td>ε<sub>acc,crit</sub></td><td>Plastic threshold value</td></tr></tbody></table>
<table><caption>**Table 2.** Cont.</caption><thead><tr><th>Symbol</th><th>Meaning</th></tr></thead><tbody><tr><td><b>H</b><sup>α</sup></td><td>Phase-inherent displacement gradients</td></tr><tr><td><b>Ĥ</b></td><td>Collection of phase-inherent displacement gradients</td></tr><tr><td><b>H</b><sup>eff</sup></td><td>Effective displacement gradient</td></tr><tr><td>ε</td><td>Total strain tensor</td></tr><tr><td>ε<sub>e</sub></td><td>Elastic strain tensor</td></tr><tr><td>ε<sub>e</sub><sup>α</sup></td><td>Phase-inherent elastic strain tensor</td></tr><tr><td>ε̂<sub>e</sub></td><td>Collection of phase-inherent elastic strains tensors</td></tr><tr><td>f<sub>e</sub></td><td>Effective elastic strain energy</td></tr><tr><td>σ</td><td>Stress tensor</td></tr><tr><td>σ<sup>α</sup></td><td>Phase-inherent stress tensor</td></tr><tr><td>σ̂</td><td>Collection of phase-inherent stress tensors</td></tr><tr><td>σ<sup>eff</sup></td><td>Effective stress tensor</td></tr><tr><td>s<sup>α</sup></td><td>Phase-inherent deviatoric stress tensor</td></tr><tr><td>ε<sub>p</sub></td><td>Plastic strain tensor</td></tr><tr><td>ε<sub>p</sub><sup>α</sup></td><td>Phase-inherent plastic strain tensors</td></tr><tr><td>ε<sub>acc</sub><sup>α</sup></td><td>Phase-inherent accumulated plastic strain</td></tr><tr><td>ε<sub>acc</sub><sup>eff</sup></td><td>Effective accumulated plastic strain</td></tr><tr><td>f<sub>y</sub><sup>α</sup></td><td>Phase-inherent yield function</td></tr><tr><td>σ<sup>α</sup><sub>y,0</sub></td><td>Phase-inherent yield strength</td></tr><tr><td><b>H</b><sup>α</sup></td><td>Phase-inherent linear hardening parameter</td></tr><tr><td>f<sub>p</sub><sup>α</sup></td><td>Phase-inherent plastic energy</td></tr><tr><td>f<sub>p</sub></td><td>Effective plastic energy</td></tr><tr><td><b>n</b><sup>αβ</sup></td><td>Unit normal vectors between the two phases</td></tr><tr><td><b>a</b></td><td>Rank-1 tensor</td></tr><tr><td><b>Q</b><sup>α</sup><sub>aniso</sub></td><td>Transformation matrices</td></tr></tbody></table>

### 3. Numerical Aspects

In this section, a brief overview of the numerical treatment is provided, which is applied in this publication.

All simulations were executed with the in-house software package PACE3D [38] (Parallel Algorithms for Crystal Evolution in 3D). The equations are implemented for general 3D problems, although simulations in quasi-2D domains, with one cell in z-direction of space, were performed exclusively, due to computational efforts. An equidistant orthogonal grid is used, since curvatures can be well approximated, using a diffuse transition between different phases. For the mapping of curvatures, a more complex discretisation grid is not necessary. Initially, a diffuse interface is established between all occurring solid phases. The calculation of the phase fields and the mechanical fields is processed in a staggered manner, which means that no monolithic solution is used, where the fields serve as an input to each other and are interpolated into the required positions of the discrete cells.

The phase-field Equation (27) is solved at central positions of the numerical cells by applying the explicit Euler scheme for the time derivative and the finite-difference method (FDM), using second-order accurate central differences for the spatial derivatives. If, due to the explicit Euler scheme, values of the crack phase field occur that are greater than one, all solid phase fields are set to zero, and $\phi_c = 1$ holds. In points where the gradient of the crack phase field is very smaller, isotropic material behaviour of the critical energy release rate, $F_{\text{aniso}} = 1$, is assumed. This is done to avoid numerical inaccuracies in the calculation of normal vectors.

On the discretisation grid, the positions of the stresses and strains correspond to an FE mesh with linear elements and a full integration. Using a predictor-corrector two-step return mapping scheme, the displacement fields are calculated locally with respect to the
mechanical jump conditions (see Section 2). The global adjustment of the displacement fields is done by a Newton-Raphson algorithm. First, initial predictor displacement fields are determined iteratively, fulfilling the mechanical equilibrium condition, given by Equation (22), including the plastic strains of the previous time step/initial values. After the calculation of the unknown plastic strains and the accumulated plastic strains, based on the elastic prediction, the mechanical equilibrium condition might be violated. From the violated mechanical equilibrium condition, a correction of the displacement fields is calculated and applied to the current displacement fields. This procedure describes one iteration step of the global Newton scheme and is repeated until the mechanical equilibrium condition holds for the new plastic strains.

For each load increment of the mechanical boundary conditions, the phase fields are brought into a steady state. Using fixed mechanical boundary conditions and constant mechanical fields, the phase-field evolution Equation (27) is solved explicitly for constant normal vectors of the crack phase field, until the criterion $\dot{\phi}_c \approx 0$ is fulfilled. After fulfilling the equilibrium condition, updated normal vectors of the crack phase field are calculated. As soon as the steady state of the phase fields holds for the updated normal vectors, the mechanical load is increased. Through the described procedure, the time increment $\Delta t$ and the kinetic coefficient $M$ of the phase-field evolution Equation (27) have only a negligible influence on the simulation results and are chosen in such a way that no numerical instabilities occur. Since the simulation results are time-independent, it is only possible to make conclusions about possible crack paths, but not about the kinetics of the crack growth.

## 4. Numerical Examples

In this section, the developed model is validated and verified. As the anisotropic brittle crack model has already been validated in our previous publication Prajapati et al. [11], and the ductile crack model for single solid phases in Ambati et al. [16], the focus of the following examinations is on the crack nucleation and propagation in domains with brittle, transversely isotropic and ductile solid phases. Section 4.1 discusses an initial crack, passing the planar diffuse interface of a brittle, transversely isotropic and ductile solid phase. In Section 4.2, the crack nucleation and propagation of a domain with a single brittle, transversely isotropic ellipsoid, which imitates a graphite lamella, embedded in a pearlitic matrix is analysed. Finally, crack nucleation and growth are simulated in a simulation area consisting of three differently oriented elliptical inclusions, embedded in a ductile matrix, which is subjected to tensile and compressive load. The simulation areas are assumed to have a physical size of 80 mm × 80 mm in order to ensure the correct use of the regularisation parameter between the solid phases and crack, where $\Delta x = \Delta y = \Delta z = 0.265$ mm is chosen as grid spacing.

The model allows the simulation of combined anisotropic and ductile properties of a single phase, which, however, is not necessary for the fictive materials under consideration, since the imitated graphite lamellae are anisotropic and brittle, while the pearlitic matrix is isotropic and ductile. In the present study, a combination of anisotropic and ductile properties is not considered.

Table 3 summarises the applied material parameters for the performed simulations. Young's modulus and the fracture toughness of the graphite lamellae are taken from Pickup et al. [39], while the Poisson's ratio is given by Fishlock et al. [40]. The pearlitic material parameters are standard values for steel. Through the introduced material parameters, the critical energy release rate $G_c$ is given as $G_c = K_{ic}^2(1 - v^2)/E$, according to a plane strain state. For the calculation of the phase-dependent regularisation parameters $l^\alpha$, according to Tanné et al. [25], tensile strengths from Zhang et al. [41] and Boyer [42] are used.
<table><caption>**Table 3.** Mechanical parameters for pearlite and graphite, used in the simulations.</caption><thead><tr><th>Properties</th><th>Pearlite</th><th>Graphite</th></tr></thead><tbody><tr><td>Young's modulus <i>E</i></td><td>210.0 GPa</td><td>12.40 GPa</td></tr><tr><td>Poisson's ratio <i>ν</i></td><td>0.3</td><td>0.1</td></tr><tr><td>Yield stress <i>σ</i><sub>y,0</sub></td><td>230.0 MPa</td><td></td></tr><tr><td>Hardening modulus <i>H</i></td><td>2.1 GPa</td><td></td></tr><tr><td>Fracture toughness <i>K</i><sub>Ic</sub></td><td>50.0 MPa√mm</td><td>1.34 MPa√mm</td></tr><tr><td>Tensile strength <i>σ</i><sub>m</sub></td><td>738.0 MPa</td><td>27.6 MPa</td></tr></tbody></table>

Transversely isotropic material stiffnesses are obtained with the help of Young's modulus, Poisson's ratio, and the anisotropy factor of Equation (25). Considering the atomic structure of graphite, made from many layered basal planes, a transversally isotropic stiffness tensor of the following form is obtained within the hexagonal lattice $\tilde{x}, \tilde{y}, \tilde{z}$-coordinate systems of graphite lamellae:

$$ \mathcal{C}_{\text{tran.}} = \begin{pmatrix} C_{22}(1 + F_{\text{aniso}}) & \lambda & C_{12}(1 + F_{\text{aniso}}) & 0 & 0 & 0 \\ \lambda & \lambda + 2\mu & C_{12} & 0 & 0 & 0 \\ & \text{sym.} & C_{22}(1 + F_{\text{aniso}}) & \frac{C_{22}-C_{12}}{2} & 0 & 0 \\ & & & \frac{C_{11}-C_{13}}{2} & 0 & 0 \\ & & & & C_{44} & \end{pmatrix}_{\tilde{x}, \tilde{y}, \tilde{z}} . \quad (30) $$

Herein, the Lamé constants are given by

$$ \lambda = \frac{\nu E}{(1 - 2\nu)(1 + \nu)}; \quad \mu = \frac{E}{2(1 + \nu)}. \qquad (31) $$

Again, it is assumed that the basal planes, pointing in the a-direction of the hexagonal crystal system, are parallel to the $\tilde{x}, \tilde{z}$-plane, as in the case of the anisotropic critical energy release rate (25). The procedure used here achieves a correlation between the anisotropy strength of the critical energy release rate and the material stiffness.

Each simulation was carried out with the same mechanical boundary conditions, which are shown in the following Section 4.1. The boundaries in x-direction are defined as stress-free. On the lower and upper boundaries of the domain, an orthogonal displacement $u$ is applied in an incremental manner. The orthogonal displacement boundary condition allows a free contraction in directions tangential to normal vectors of free surfaces. To create a condition of the quasi-2D domain corresponding to a plane strain state, the orthogonal displacement boundary condition is also applied in the z-direction, where the orthogonal displacements are equal to zero. After the condition $\dot{\phi}_c < 10^{-4}$ holds for the updated normal vectors of the crack phase field, the mechanical load is increased.

The visualisation of the crack phase field $\phi_c$ is done by overlaying the solid phases and a continuous transition from opaque (black), $\phi_c = 1$, to transparent, $\phi_c = 0$.

### 4.1. Planar Interface between Brittle-Anisotropic and Ductile Solid Phase

The investigated problem consists of the transition of an anisotropic brittle phase into a ductile phase, by means of a planar interface, whereby the left half of the domain corresponds to the anisotropic brittle phase and the right half to the ductile phase. On the left side of the simulation domain, a pre-existing 20 mm long initial crack is located. In Figure 1, the simulation setup and the mechanical boundary conditions are depicted. Only uniaxial tensile tests were carried out, whose mechanical boundary conditions are specified at the beginning of this Section 4.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Simulation domain for a crack passing a vertical planar diffuse interface, accompanied by the mechanical boundary conditions. The dotted lines indicate the diffuse transition between the two solid phases.</figcaption></figure>

In this example, both phases use the material properties of pearlite (Table 3) to ensure transcrystalline crack growth. However, the brittle phase still makes use of the transversely isotropic stiffness tensor (Equation (30)) with no plastic deformations occurring.

In the following studies, only the influence of a single parameter is examined at once. Therefore, the following set of standard model parameters is defined, which is kept constant, while one parameter is varied:

$$ \epsilon = 3\Delta x; \quad \epsilon_{\text{acc,crit}}^{\text{d}} = 0.08; \quad F_{\text{aniso}}^{\text{b}} = 2; \quad \psi_{\text{aniso}}^{\text{b}} = 30^{\circ}. \qquad (32) $$

The influences of the following parameters on the crack propagation and the homogenised macroscopic stress–strain response are investigated: the parameter $\epsilon$, which defines the width of the diffuse transition between the solid phases; the rotation angle $\psi_{\text{aniso}}^b$; and the anisotropy factor $F_{\text{aniso}}^b$ of the transverse isotropy around the z-axis of the reference x, y, z-coordinate system, as well as the threshold value $\epsilon_{\text{acc,crit}}^d$ of the ductile phase.

First of all, the influence of the length parameter for the definition of the width of the diffuse solid-solid phase transitions is investigated for the values $\epsilon = \{3, 5, 7\}\Delta x$, in terms of the homogenised macroscopic stress–strain response and the crack morphology at fracture (see Figure 2a). The corresponding contour plots of the equivalent von Mises stresses of Equation (23), $\sigma_{\text{v.M.}}^{\text{eff}}$, and the effective accumulated plastic strain, $\epsilon_{\text{acc}}^{\text{eff}}$, for the highlighted states in Figure 2a, as well as the set of standard parameters, Equation (32), are shown in Figure 3. It can be observed that a linear elastic range is present in the macroscopic stress–strain diagrams ($t_0 - t_2$) before the stress value drops abruptly after its peak ($t_2$), in the same manner as brittle materials, without showing a macroscopic plastic deformation. The steep descent ($t_2 - t_3$) is caused by an unstable crack growth within the brittle anisotropic phase. As can be seen in the polar plots, overlaying the final crack paths in Figure 2a, the angle between the crack in the anisotropic phase and the x-axis of the reference coordinate system is $\psi = 19.2^{\circ}$. The discrepancy between the inclination of the crack and the orientation of the anisotropy occurs because the elastic strain energy tends to target a perfectly horizontal crack ($\psi = 0^{\circ}$) when exposed to a load of Mode-I, whereas the formulation of the critical energy release rate forces a slope, corresponding to $\psi_{\text{aniso}}^b$. These two counterparts cause an actual inclination of the crack between $\psi = 0^{\circ}$ and $\psi_{\text{aniso}}^b$ [11]. As soon as the crack reaches the diffuse interface between the solid phases, the characteristics of the macroscopic stress–strain plots and the crack morphology change. The further course
of the macroscopic stress–strain response is strongly nonlinear ($t_4$) and is associated with plastic deformations. The crack path loses its inclination and spreads horizontally in a pure Mode-I manner. After the sample is broken ($t_5$), the macroscopic stress drops to zero, while on the level of the contour plots, residual stresses remain locally, caused by the plastic deformations. In Figure 2a, it can be seen that the width of the diffuse transition between the solid phases has almost no influence on the macroscopic stress–strain curve and the final crack morphology in the considered range for $\epsilon$.

Next, the anisotropy orientation for the values $\psi_{\text{aniso}}^{\text{b}} = \{0, 15, 30, 45, -30\}^{\circ}$ is varied. The corresponding stress–strain curves and the final crack paths are depicted in Figure 2b. Due to the anisotropic stiffness tensor of the brittle phase, the peak stress value of the diagram increases with an increasing anisotropy angle $\psi_{\text{aniso}}^{\text{b}}$. In analogy to graphite, only the van der Waals forces are stressed at a rotation angle of $\psi_{\text{aniso}}^{\text{b}} = 0^{\circ}$, whereas only the covalent bonds are stressed at $\psi_{\text{aniso}}^{\text{b}} = 90^{\circ}$, resulting in higher macroscopic stress values. As the rotation angle of the anisotropy increases, the deflection of the resulting crack also increases as expected, whereby it should again be noted that the anisotropic formulation and the strain energy act as opponents. When comparing $\psi_{\text{aniso}}^{\text{b}} = 30^{\circ}$ and $\psi_{\text{aniso}}^{\text{b}} = -30^{\circ}$ in Figure 2b, it can be observed that the resulting stress–strain diagrams for a change of sign of the rotation angle are almost equivalent. The small differences that occur are due to numerical inaccuracies that arise because the mechanical load is only applied at the upper edge of the simulation area. For varying orientations of anisotropy, the ductile part of the stress–strain curves does not differ significantly; all curves are nearly identical, except for $\psi_{\text{aniso}}^{\text{b}} = 45^{\circ}$. Due to geometric reasons, the cracks become wider with increasing deflection within the brittle phase, as they pass through the solid phase transition region. This phenomenon has an effect on the stress–strain curve of the simulation with an anisotropy orientation of $\psi_{\text{aniso}}^{\text{b}} = 45^{\circ}$.

In Figure 4a, the influence of the anisotropy factor on the mechanical behaviour is discussed for the values $F_{\text{aniso}}^{\text{b}} = \{1, 1.5, 2, 2.5, 3\}$. The stress–strain curve shows higher stress peaks as the anisotropy factor increases and thus also the effective stiffness of the simulation domain. The reason for this is that the anisotropy factor increases the stiffness of the brittle phase. In analogy to lamellar graphite, the increase in the elastic modulus occurs within the basal planes along the a-direction of the hexagonal crystal system. If the anisotropy factor is chosen as $F_{\text{aniso}}^{\text{b}} = 1$, the brittle phase behaves isotropically, and the crack profile has no slope, $\psi = 0^{\circ}$. The deflection of the crack increases with increasing anisotropy factor, due to the anisotropic critical energy release rate.

Finally, the influence of the threshold value $\epsilon_{\text{acc,crit}}^{\text{d}}$ of the ductile phase is investigated for the values $\epsilon_{\text{acc,crit}}^{\text{d}} = \{0.8, 0.6, 0.4, 0.2\}$. The corresponding results are given in Figure 4b. As long as the crack propagates within the brittle phase, a change of the threshold value has no influence on the material behaviour, neither with regard to the stress–strain behaviour nor with regard to the crack morphology. In the ductile phase, the crack morphology is also not significantly altered by the threshold value. However, as expected, there is a strong influence on the course of the stress–strain diagram. As the threshold increases, the elongation until complete failure rises, because the ductile degradation function depends on $\epsilon_{\text{acc,crit}}^{\text{d}}$, which means that higher plastic deformations are required to cause crack propagation in the ductile phase as $\epsilon_{\text{acc,crit}}^{\text{d}}$ increases.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Homogenised macroscopic stress–strain responses in the case of uniaxial tension for different values of $\epsilon = \{3, 5, 7\}\Delta x$, for the simulation setup with a planar diffuse interface between two solid phases. The dots in the stress–strain diagram correspond to the strain states of the contour plots of Figure 3. Below the diagram, contour plots of the crack phase $\phi_c$ at fracture are given, overlaying the solid phases for the various values of $\epsilon$. (b) Homogenised macroscopic stress–strain responses in the case of uniaxial tension for different values of $\psi_{\text{aniso}}^b = \{0, 15, 30, 45, -30\}^\circ$, for the simulation setup with a planar diffuse interface between two solid phases. The contour plots of the crack phase $\phi_c$ at fracture, overlaying the solid phases, for the different values of $\psi_{\text{aniso}}^b$ are given below and beneath the stress–strain diagram.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Contour plots of the evolution of the crack phase field $\phi_c$ (first column) and the resulting von Mises stresses $\sigma_{\text{v.M.}}^{\text{eff}}$ (second column), together with the effective accumulated plastic strain $\varepsilon_{\text{acc}}^{\text{eff}}$ (third column) for the simulation setup with a planar diffuse interface between two solid phases for the set of standard parameters, Equation (32). The corresponding macroscopic uniaxial strain states are given in Figure 2a.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Homogenised macroscopic stress–strain responses in the case of uniaxial tension for different values of $F_{\text{aniso}}^{\text{b}} = \{1, 1.5, 2, 2.5, 3\}$, for the simulation setup with a planar diffuse interface between two solid phases. The contour plots of the crack phase $\phi_c$ at fracture, overlaying the solid phases, for the different values of $F_{\text{aniso}}^{\text{b}}$ are given below and beneath the stress–strain diagram. (b) Homogenised macroscopic stress–strain responses in the case of uniaxial tension for different values of $\epsilon_{\text{acc,crit}}^{\text{d}} = \{0.08, 0.06, 0.04, 0.02\}$, for the simulation setup with a planar diffuse interface between two solid phases. The contour plots of the crack phase $\phi_c$ at fracture, overlaying the solid phases, for the different values of $\epsilon_{\text{acc,crit}}^{\text{d}}$ are given below and beneath the stress–strain diagram.</figcaption></figure>

## 4.2. Single Elliptic, Transversely Isotropic Brittle Inclusion in Ductile Matrix

In this section, a **single elliptic, transversely isotropic brittle inclusion**, embedded in a **ductile matrix**, is used as a simulation setup. The longitudinal axis of the ellipsoid, which is equivalent to the $\tilde{x}$-direction of the anisotropic critical energy release rate, is rotated
with respect to the z-axis of the reference coordinate system. Its length is 64 mm and its width 4 mm. In Figure 5, the simulation setup and the mechanical boundary conditions are depicted. For the uniaxial tensile tests, the mechanical boundary conditions, introduced at the beginning of this Section 4, are applied.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Setup of the crack development investigation in a domain consisting of an elliptic, transversely isotropic brittle inclusion, embedded in a ductile matrix, accompanied by the mechanical boundary conditions. The ductile matrix is represented in **blue**, while the brittle inclusion phase is represented in **orange**.</figcaption></figure>

As material parameters for the anisotropic brittle phase, the values of graphite given in Table 3 are used and applied to Equation (30). The ductile matrix possesses the mechanical properties of pearlite as given in Table 3. In the following investigations, only the influence of one parameter is examined at a time, as was done previously. Therefore, the set of standard model parameters, Equation (32), is used, which is kept constant in case another parameter is varied.

First of all, the length parameter to define the width of the diffuse transitions between solid phases is varied for the values $\epsilon = \{3, 5, 7\}\Delta x$. The influence of the length parameter is investigated in terms of the stress–strain behaviour and the crack morphology (see Figure 6a). In Figure 7, the corresponding contour plots of the equivalent von Mises stress of Equation (23), $\sigma_{\text{v.M.}}^{\text{eff}}$, and the effective accumulated plastic strain, $\epsilon_{\text{acc}}^{\text{eff}}$, for the highlighted states in Figure 6a, as well as the set of standard parameters, are shown. A linear elastic region is present in the stress–strain diagrams, up to a strain of $\epsilon \approx 0.075\%$ ($t_0-t_1$), before crack nucleation takes place. Based on the evolution of the crack phase, as well as on the stress–strain plots, it can be observed that there are three stages of crack growth. First, a crack spreads within the brittle phase. This crack growth is accompanied by a drop of the stress–strain plot ($t_1-t_2$), which proceeds in a brittle manner. After the stress reaches a critical value ($t_3$), the crack tips start to propagate through the diffuse interface. During this stage, the stress–strain plot has an almost constant stress value ($t_3-t_4$). Afterwards, the crack propagates from the upper and lower ends of the ellipsoid and grows horizontally until a complete failure occurs, where the stress–strain plots drop to zero ($t_4-t_6$). After the rupture, residual stresses exist, due to the presence of plastic strains. During crack propagation in the ductile phase, a pronounced nonlinear behaviour can be observed.

The name grey cast iron originates from the fractured surfaces of broken grey cast iron components, which usually have a grey shimmer. Such greyish fracture surfaces are the result of a preferred crack growth along graphite lamellae, due to the significantly
lower critical energy release rate of graphite lamellae, compared to the metallic matrix [26]. Under tensile load, graphite lamellae have a negligible loading capacity and therefore act as microcracks even at very small strains [9,28]. The edges of the ruptured graphite lamellae act as stress raisers, which induce local plastic zones [9]. Only when the graphite lamellae have lost their load-bearing capacity does the ductile matrix begin to break. The simulation results show that this behaviour can be mapped by the introduced model. Crack nucleation exclusively takes place within the brittle phase, associated with graphite lamellae, whereby the direction of the crack growth along the brittle phase is favoured by the anisotropic formulation of the critical energy release rate, before the crack propagation takes place within the ductile matrix. The crack propagation within the ductile phase is accompanied by local plastic strains occurring in front of the crack tips (see Figure 7) ($t_3$). Within the matrix, a Mode-I-type crack growth predominantly occurs in the horizontal direction.

When defining the diffuse width of the solid phase transitions, a slight dependence of the stress–strain curves on the length parameter can be seen, whereby a convergence of the material behaviour can be guessed, when the length parameter increases. The crack morphology at fracture is independent of the chosen solid-solid phase transition width, as can be seen by the final crack paths in Figure 6a—bottom. In areas where the simulation domain has multiple coexisting phases, the calculation of the mechanical fields is significantly more expensive, compared to areas with only one phase. Since the length parameter has no strong influence on the simulation results, it is sufficient to choose a length parameter value that is as small as possible but still represents a good compromise between accuracy and calculation time.

Subsequently, the orientation of the brittle ellipsoid is varied. The rotation of the longitudinal axis of the ellipsoid around the z-axis of the reference coordinate system for the values $\psi_{\text{aniso}}^{\text{b}} = \{0, 15, 30, 45, -30\}^{\circ}$ is correlated to the rotation of the $\tilde{x}, \tilde{y}, \tilde{z}$-coordinate system of the anisotropy formulation. In Figure 6b, the corresponding stress–strain curves and crack paths at the fracture are shown. As one can see, the stress–strain curves are heavily dependent on the orientation of the inclusion. The simulation with a rotation angle of $\psi_{\text{aniso}}^{\text{b}} = 0^{\circ}$ shows the smallest maximum stress value. In analogy to lamellar graphite, the loading takes place exactly orthogonally to the basal planes, so that the highest value of the elastic strain energy coincides with the direction of the lowest crack resistance. The gap between the maximum elastic strain energy and the minimum crack resistance grows as the rotation angle increases, which results in the fact that greater stresses are required to cause cracking. For this reason, the simulation with a rotation angle of $\psi_{\text{aniso}}^{\text{b}} = 45^{\circ}$ has the highest maximum stress. Regardless of the rotation angle, the crack always develops along the longitudinal axes of the brittle particles. When comparing $\psi_{\text{aniso}}^{\text{b}} = 30^{\circ}$ and $\psi_{\text{aniso}}^{\text{b}} = -30^{\circ}$ in Figure 6b, it can be observed that the resulting stress–strain diagrams for a change of sign of the rotation angle are matching.

The influence of the anisotropy factor on the mechanical behaviour is shown in Figure 8a. The results correspond to the values $F_{\text{aniso}}^{\text{b}} = \{1, 1.5, 2, 2.5, 3\}$. The stress–strain curves are almost identical. As in the previous example, the increase in the stiffness components of the brittle phase, which is due to the anisotropy factor, only has a small influence on the macroscopic material behaviour, as a microcrack already develops inside the ellipsoid at very small strains. The only difference is observed between the maximum values of the stresses. In the nonlinear area of the diagram, which is characterised by plastic deformations, the graphs are identical. On the scale of the crack phase field, it is noticeable that the crack becomes thinner as the anisotropy factor increases, since the strength of the directional dependence increases.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Homogenised macroscopic stress–strain responses in the case of uniaxial tension for different values of $\epsilon = \{3, 5, 7\}\Delta x$, for the simulation setup with a single elliptic, transversely isotropic, brittle inclusion in a ductile matrix. The dots in the stress–strain diagram correspond to the strain states of the contour plots of Figure 7. Below the diagram, contour plots of the crack phase $\phi_c$ at fracture are given, overlaying the solid phases for the various values of $\epsilon$. (b) Homogenised macroscopic stress–strain responses in the case of uniaxial tension for different values of $\psi_{\text{aniso}}^b = \{0, 15, 30, 45, -30\}^\circ$, for the simulation setup with a single elliptic, transversely isotropic, brittle inclusion in a ductile matrix. The contour plots of the crack phase $\phi_c$ at fracture, overlaying the solid phases, for the different values of $\psi_{\text{aniso}}^b$ are given below and beneath the stress–strain diagram.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Contour plots of the evolution of the crack phase field $\phi_c$ (first column) and the resulting von Mises stresses $\sigma_{\text{v.M.}}^{\text{eff}}$ (second column), together with the effective accumulated plastic strain $\varepsilon_{\text{acc}}^{\text{eff}}$ (third column) for the simulation setup with a single elliptic, transversely isotropic, brittle inclusion, embedded in a ductile matrix, for the set of standard parameters, Equation (32). The corresponding macroscopic uniaxial strain states are given in Figure 6a.</figcaption></figure>

Finally, the influence of the threshold value $\varepsilon_{\text{acc,crit}}^{\text{d}}$ of the ductile phase is investigated for the values $\varepsilon_{\text{acc,crit}}^{\text{d}} = \{0.8, 0.6, 0.4, 0.2\}$. The corresponding results are given in Figure 8b. Again, a change of the threshold value has no influence on the material behaviour, as long as the crack propagates within the brittle phase, neither with regard to the stress–strain behaviour nor with regard to the crack morphology. In the ductile phase, the crack morphology is also not changed significantly by the threshold value. With an increasing plastic threshold, the elongation at the break rises, as the ductile degradation function depends on $\varepsilon_{\text{acc,crit}}^{\text{d}}$, which means that higher plastic deformations are required to cause crack propagation in the ductile phase when $\varepsilon_{\text{acc,crit}}^{\text{d}}$ does increase.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** (a) Homogenised macroscopic stress–strain responses in the case of uniaxial tension for different values of $F_{\text{aniso}}^{\text{b}} = \{1, 1.5, 2, 2.5, 3\}$, for the simulation setup with a single elliptic, transversely isotropic, brittle inclusion in a ductile matrix. The contour plots of the crack phase $\phi_c$ at fracture, overlaying the solid phases, for the different values of $F_{\text{aniso}}^{\text{b}}$ are given below and beneath the stress–strain diagram. (b) Homogenised macroscopic stress–strain responses in the case of uniaxial tension for different values of $\epsilon_{\text{acc,crit}}^{\text{d}} = \{0.08, 0.06, 0.04, 0.02\}$, for the simulation setup with a single elliptic, transversely isotropic, brittle inclusion in a ductile matrix. The contour plots of the crack phase $\phi_c$ at fracture, overlaying the solid phases, for the different values of $\epsilon_{\text{acc,crit}}^{\text{d}}$ are given below and beneath the stress–strain diagram.</figcaption></figure>

### 4.3. Multiphase Simulations in an Idealised Grey Cast Iron Microstructure

As a more realistic application example, the developed model is applied to an idealised grey cast iron microstructure. As simplified representations of graphite lamellae,
the simulation domain contains three differently oriented elliptic, transversely isotropic, brittle inclusions in a ductile matrix. The length of 26.67 mm and the width of 3.33 mm are the same for all inclusions. Figure 9 shows the simulation setup and the applied mechanical boundary conditions (see Section 4), which have already been applied in the previous validation cases. Between the solid phases, a diffuse interface, corresponding to a width of $\epsilon = 3\Delta x$, is generated, while no solid-solid phase transformation takes place during the simulations. Since the simulation domain is not a representative volume element in terms of a realistic cast iron with the lamellar graphite microstructure, the brittle particles extend almost through the entire domain.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Simulation domain with three differently oriented elliptic, transversely isotropic, brittle inclusions in a ductile matrix, accompanied by the applied mechanical boundary conditions. The anisotropic brittle phases are represented in **orange**, while the ductile matrix is represented in **blue**.</figcaption></figure>

The material parameters for the ductile matrix and the brittle inclusions are given in Table 3. The anisotropy factor $F_{aniso}^b = 1.5$ and the plastic threshold value $\epsilon_{acc,crit}^d = 0.08$ were selected as model parameters. With the anisotropy factor, the unrotated transversal isotropic stiffness tensor of the anisotropic ellipsoids is obtained according to Equation (30).

Figures 10 and 11 show the evolution of the crack phase field under uniaxial tensile and compressive load. In addition, contour plots of the von Mises stresses $\sigma_{v.M.}^{eff}$ and the effective accumulated plastic strain, $\epsilon_{acc}^{eff}$, are given. All corresponding states of Figures 10 and 11 are highlighted in the stress-strain plot in Figure 12. In both loading cases, the crack nucleation takes place in the ellipsoid, whose longitudinal axis is oriented orthogonally to the loading direction. The weakest critical energy release rate lies exactly along that longitudinal axis. In analogy to physical reality, in this case, the weak van der Waals bonds of this graphite lamella are mainly stressed, which results in the observed growth behaviour. The crack growth within the horizontal ellipsoid is followed by local stress peaks at the ellipsoid's edges and associated with plastic deformations of the pearlitic matrix. The described behaviour is consistent with experimental observations [9] and has already been observed in Section 4.2. After the first ellipsoid is fractured, the ellipsoid which shows a rotation of the longitudinal axis by $\psi = -45^\circ$ with respect to the reference coordinate system starts to break. Under compressive load, a crack also develops within the vertical ellipsoid. When pressure is applied, significantly higher plastic deformations occur, as expected. After the specimens are broken, the zones under plastic strain induce a complex distribution of residual stresses, which is typical for grey cast iron [9]. Under tensile load, the crack morphology within the ductile matrix is narrower than under compressive load. Overall, the behaviour observed in the simulations—that microcracks first
form within the graphite particles and that these microcracks subsequently coalesce into macrocracks—is consistent with experimental results (see, e.g., [43–45]).

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Evolution of the crack phase field $\phi_c$ (first column) and the resulting von Mises stresses $\sigma_{\text{v.M.}}^{\text{eff}}$ (second column), together with the effective accumulated plastic strain $\varepsilon_{\text{acc}}^{\text{eff}}$ (third column) of an idealised grey cast iron microstructure. Uniaxial tensile load was applied. The corresponding macroscopic uniaxial strain states are given in Figure 12.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Evolution of the crack phase field $\phi_c$ (first column) and the resulting von Mises stresses $\sigma_{\text{v.M.}}^{\text{eff}}$ (second column), together with the effective accumulated plastic strain $\epsilon_{\text{acc}}^{\text{eff}}$ (third column) of an idealised grey cast iron microstructure. Uniaxial compressive load was applied. The corresponding macroscopic uniaxial strain states are given in Figure 12.</figcaption></figure>

With regard to the type of load, important characteristics of the complicated constitutive behaviour of grey cast iron can be found in the simulation results. This includes the nonexistence of a linear elastic regime under tensile load and the tension–compression–stress asymmetry [46]. Figure 12 clearly depicts the asymmetric behaviour under tensile and compressive load, which is typical for cast iron with lamellar graphite on the macroscopic scale [47,48]. Under tensile load, the tested material possesses a brittle behaviour: at an elongation of $\epsilon \approx 0.2\%$ and a stress peak of $\sigma \approx 95$ MPa, crack nucleation takes place. Crack growth occurs without the occurrence of large macroscopic plastic deformations and leads to a fracture at an elongation of $\epsilon \approx 1.63\%$. Keeping the size of the simulation area in mind, these results are in relatively good agreement with the experimental findings documented in the literature. According to Noguchi and Shimizu [49], grey cast iron does not show a pronounced yield strength under tensile load and breaks at a strain between $\epsilon = 0.5–1.0\%$, without visible necking. In grey cast irons, tensile fracture stresses are
in the range of $\sigma = 100–500$ MPa [50]. Under compressive load, a drop in the stress–strain diagram occurs after a strain of $\varepsilon \approx 0.2\%$ and a stress peak of $\sigma \approx 134.5$ MPa. Before this point, the course of the curve is linear. Afterwards, a large segment occurs, which shows nonlinear material behaviour, caused by the plastic deformation of the ductile matrix. The elongation at fracture under compressive load is $\varepsilon > 5\%$ and can therefore no longer be properly represented when small deformations are assumed. This different behaviour under tensile and compressive load is mainly driven by the tension–compression split of the elastic strain energy. Under compressive load, the elastic strain energy is considerably smaller compared to the tensile example; therefore, failure of the ductile matrix occurs at significantly higher strains.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Homogenised macroscopic stress–strain response of a simulation domain, inspired by an idealised cast iron microstructure with lamellar graphite. The diagrams show the stress and strain in the case of uniaxial tension and compression. The dots in the stress–strain diagram correspond to the strain states of the contour plots of Figures 10 and 11. Even though the simulation domain is not a representative volume element, a clear tension–compression dependency of the stress–strain plot can be observed when looking at the peak stresses and the elongations at failure.</figcaption></figure>

**5. Summary and Outlook**

In this work, a small-strain phase-field model capable of predicting crack propagation in systems with anisotropic brittle and ductile constituents was presented. To distinguish between the brittle and ductile crack characteristics, two different formulations of the degradation function are used. As an extension to our previous model formulation [11], the ductile degradation function by Ambati et al. [16] was also introduced. To simulate nonlinear mechanical behaviour, a model, based on the $J_2$-plasticity theory is applied, which fulfils the mechanical jump conditions for the phase with the largest volume fraction. Idealised simulation setups were used to investigate the influences of the relevant model parameters on cracks passing through a planar interface between a brittle and a ductile solid phase and for crack development in a domain of a single brittle ellipsoid, embedded in a ductile matrix. Afterwards, uniaxial tension and compression tests were performed in a domain, inspired by a grey cast iron microstructure. It was shown that important mechanical grey cast iron properties can be mapped by the model, which makes it suitable for the simulation of that material group. The model captures the initiation of microcracks within graphite lamellae and the subsequent coalescence of these microcracks by crack growth processes within the ductile matrix to form a fracture. On the macroscopic length
scale, the tension–compression load asymmetry of the uniaxial stress–strain response
is mapped.

An important and necessary improvement of the model is to make it independent of
the regularisation parameter, which defines the width of the diffuse transition between the
solid phases and the crack. With such a modification of the model, the crack development
in cast iron microstructures could be simulated on its true physical length scale in the
millimetre range. This could be done in a similar way as published by Wu and Nguyen [37]
but with the additional difficulty that a brittle-anisotropic model formulation, combined
with a ductile model formulation, would have to be considered, instead of a brittle-isotropic
model formulation.

In future work, we intend to investigate cracking processes in grey cast iron materials under realistic thermomechanical load, following braking operations of truck brake discs. Furthermore, we want to investigate the crack formation in grey cast iron that has undergone the martensitic transformation process. In addition, the use of a finite deformation model would be of interest. Furthermore, the model can be applied or extended to other types of material composites, composed of multiple phases with different mechanical properties such as fibre-reinforced polymer structures, fuel cells, or battery materials.

**Author Contributions:** Conceptualization, C.H., D.S., E.S. and F.S.; Investigation, C.H.; Methodology, C.H.; Software, C.H., D.S., E.S. and F.S.; Supervision, B.N.; Validation, C.H.; Visualization, C.H.; Writing—original draft, C.H.; Writing—review and editing, D.S., E.S., F.S. and B.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** The authors thank the DFG ("Deutsche Forschungsgemeinschaft") for funding the research work through the transfer project NE 822/16-1 of IRTG 2078, the Gottfried-Wilhelm Leibniz prize NE822/31-1, and POLiS Cluster of Excellence, Project ID 390874152.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the
corresponding author.

**Acknowledgments:** Modelling concepts have also been developed in collaboration with the DFG Research Training Group 2078. To explore the directions of further applications, a discussion was furthermore established within the Helmholtz programme Materials System Engineering (MSE) (43.31.01). The work contributes to the research performed at CELEST (Centre for Electrochemical Energy Storage Ulm & Karlsruhe). The authors gratefully acknowledge the editorial support by Leon Geisen.

**Conflicts of Interest:** The authors declare no conflict of interest.

**References**

1. Elices, M.G.G.V.; Guinea, G.V.; Gomez, J.; Planas, J. The cohesive zone model: Advantages, limitations and challenges. *Eng. Fract. Mech.* **2002**, *69*, 137–163. [CrossRef]

2. Brocks, W.; Cornec, A.; Scheider, I. *Computational Aspects of Nonlinear Fracture Mechanics*; GKSS Forschungszentrum Geesthacht GmbH: Geesthacht, Germany, 2003.

3. Besson, J. Continuum models of ductile fracture: A review. *Int. J. Damage Mech.* **2010**, *19*, 3–52. [CrossRef]

4. Miehe, C.; Welschinger, F.; Hofacker, M. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. *Int. J. Numer. Methods Eng.* **2010**, *83*, 1273–1311. [CrossRef]

5. Ambati, M.; Gerasimov, T.; De Lorenzis, L. A review on phase-field models of brittle fracture and a new fast hybrid formulation. *Comput. Mech.* **2015**, *55*, 383–405. [CrossRef]

6. Alessi, R.; Ambati, M.; Gerasimov, T.; Vidoli, S.; De Lorenzis, L., Comparison of phase-field models of fracture coupled with plasticity. In *Advances in Computational Plasticity: A Book in Honour of D. Roger J. Owen*; Springer International Publishing: Cham, Switzerland, 2018; pp. 1–21.

7. Griffith, A.A. VI. The phenomena of rupture and flow in solids. *Philos. Trans. R. Soc. London. Ser. A Contain. Pap. A Math. Phys. Character* **1921**, *221*, 163–198.

8. Kelly, B.T. *Physics of Graphite*; Applied Science Publishers Ltd.: London, UK, 1981.

9. Haenny, L.; Zambelli, G. Strain mechanisms in grey cast iron. *Eng. Fract. Mech.* **1983**, *18*, 377–387. [CrossRef]
10. Pina, J.C.; Kouznetsova, V.G.; Geers, M.G.D. Thermo-mechanical analyses of heterogeneous materials with a strongly anisotropic phase: The case of cast iron. *Int. J. Solids Struct.* **2015**, *63*, 153–166. [CrossRef]

11. Prajapati, N.; Herrmann, C.; Späth, M.; Schneider, D.; Selzer, M.; Nestler, B. Brittle anisotropic fracture propagation in quartz sandstone: Insights from phase-field simulations. *Comput. Geosci.* **2020**, *24*, 1–16. [CrossRef]

12. Teichtmeister, S.; Kienle, D.; Aldakheel, F.; Keip, M.A. Phase field modeling of fracture in anisotropic brittle solids. *Int. J. Non-Linear Mech.* **2017**, *97*, 1–21. [CrossRef]

13. Nguyen, T.T.; Réthoré, J.; Yvonnet, J.; Baietto, M.C. Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials. *Comput. Mech.* **2017**, *60*, 289–314. [CrossRef]

14. Bleyer, J.; Alessi, R. Phase-field modeling of anisotropic brittle fracture including several damage mechanisms. *Comput. Methods Appl. Mech. Eng.* **2018**, *336*, 213–236. [CrossRef]

15. Duda, F.P.; Ciarbonetti, A.; Sánchez, P.J.; Huespe, A.E. A phase-field/gradient damage model for brittle fracture in elastic-plastic solids. *Int. J. Plast.* **2015**, *65*, 269–296. [CrossRef]

16. Ambati, M.; Gerasimov, T.; De Lorenzis, L. Phase-field modeling of ductile fracture. *Comput. Mech.* **2015**, *55*, 1017–1040. [CrossRef]

17. Ambati, M.; Kruse, R.; De Lorenzis, L. A phase-field model for ductile fracture at finite strains and its experimental verification. *Comput. Mech.* **2016**, *57*, 149–167. [CrossRef]

18. Miehe, C.; Hofacker, M.; Schänzel, L.M.; Aldakheel, F. Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids. *Comput. Methods Appl. Mech. Eng.* **2015**, *294*, 486–522. [CrossRef]

19. Miehe, C.; Kienle, D.; Aldakheel, F.; Teichtmeister, S. Phase field modeling of fracture in porous plasticity: A variational gradient-extended eulerian framework for the macroscopic analysis of ductile failure. *Comput. Methods Appl. Mech. Eng.* **2016**, *312*, 3–50. [CrossRef]

20. Kuhn, C.; Noll, T.; Müller, R. On phase field modeling of ductile fracture. *GAMM-Mitteilungen* **2016**, *39*, 35–54. [CrossRef]

21. Miehe, C.; Aldakheel, F.; Raina, A. Phase field modeling of ductile fracture at finite strains: A variational gradient-extended plasticity-damage theory. *Int. J. Plast.* **2016**, *84*, 1–32. [CrossRef]

22. Herrmann, C.; Schoof, E.; Schneider, D.; Schwab, F.; Reiter, A.; Selzer, M.; Nestler, B. Multiphase-field model of small strain elasto-plasticity according to the mechanical jump conditions. *Comput. Mech.* **2018**, *62*, 1399–1412. [CrossRef]

23. Schneider, D.; Schoof, E.; Huang, Y.; Selzer, M.; Nestler, B. Phase-field modeling of crack propagation in multiphase systems. *Comput. Methods Appl. Mech. Eng.* **2016**, *312*, 186–195. [CrossRef]

24. Ambrosio, L.; Tortorelli, V.M. Approximation of functional depending on jumps by elliptic functional via Γ-convergence. *Commun. Pure Appl. Math.* **1990**, *43*, 999–1036. [CrossRef]

25. Tanné, E.; Li, T.; Bourdin, B.; Marigo, J.J.; Maurini, C. Crack nucleation in variational phase-field models of brittle fracture. *J. Mech. Phys. Solids* **2018**, *110*, 80–99. [CrossRef]

26. Macherauch, E.; Zoch, H.W. Gusseisen unter Zug- und Druckbeanspruchung. In *Praktikum in Werkstoffkunde: 91 ausführliche Versuche aus wichtigen Gebieten der Werkstofftechnik*; Vieweg+Teubner: Wiesbaden, Germany, 2011.

27. Wiese, J.W.; Dantzig, J.A. Modeling stress development during the solidification of gray iron castings. *Metall. Trans. A* **1990**, *21*, 489–497. [CrossRef]

28. Collini, L.; Nicoletto, G.; Konečná, R. Microstructure and mechanical properties of pearlitic gray cast iron. *Mater. Sci. Eng. A* **2008**, *488*, 529–539. [CrossRef]

29. Silhavy, M. *The Mechanics and Thermodynamics of Continuous Media*; Springer: Berlin/Heidelberg, Germany, 1997.

30. Beckermann, C.; Diepers, H.J.; Steinbach, I.; Karma, A.; Tong, X. Modeling melt convection in phase-field simulations of solidification. *J. Comput. Phys.* **1999**, *154*, 468–496. [CrossRef]

31. Moelans, N.; Blanpain, B.; Wollants, P. Quantitative analysis of grain boundary properties in a generalized phase field model for grain growth in anisotropic systems. *Phys. Rev. B* **2008**, *78*, 024113. [CrossRef]

32. Simo, J.C.; Hughes, T.J.R. *Computational Inelasticity*; Springer: Berlin/Heidelberg, Germany, 2013.

33. Holmgren, D.; Källbom, R.; Svensson, I.L. Influences of the graphite growth direction on the thermal conductivity of cast iron. *Metall. Mater. Trans. A* **2007**, *38*, 268–275. [CrossRef]

34. Shebatinov, M.P. A study of the fine structure of graphite inclusions in gray cast irons by means of the scanning electron microscope. *Met. Sci. Heat Treat.* **1974**, *16*, 288–293. [CrossRef]

35. Cahn, J.W.; Allen, S.M. A microscopic theory for domain wall motion and its experimental verification in Fe-Al alloy domain growth kinetics. *J. Phys. Colloq.* **1977**, *38*, C7-51–C7-54. [CrossRef]

36. Schneider, D.; Schwab, F.; Schoof, E.; Reiter, A.; Herrmann, C.; Selzer, M.; Böhlke, T.; Nestler, B. On the stress calculation within phase-field approaches: A model for finite deformations. *Comput. Mech.* **2017**, *60*, 203–217. [CrossRef]

37. Wu, J.Y.; Nguyen, V.P. A length scale insensitive phase-field damage model for brittle fracture. *J. Mech. Phys. Solids* **2018**, *119*, 20–42. [CrossRef]

38. Hötzer, J.; Reiter, A.; Hierl, H.; Steinmetz, P.; Selzer, M.; Nestler, B. The parallel multi-physics phase-field framework Pace3D. *J. Comput. Sci.* **2018**, *26*, 1–12. [CrossRef]

39. Pickup, I.M.; McEnaney, B.; Cooke, R.G. Fracture processes in graphite and the effects of oxidation. *Carbon* **1986**, *24*, 535–543. [CrossRef]
40. Fishlock, S.J.; Grech, D.; McBride, J.W.; Chong, H.M.H.; Pu, S.H. Mechanical characterisation of nanocrystalline graphite using micromechanical structures. *Microelectron. Eng.* **2016**, *159*, 184–189. [CrossRef]

41. Zhang, X.; Yi, Y.; Zhu, H.; Liu, G.; Sun, L.; Shi, L.; Jiang, H.; Ma, S. Measurement of tensile strength of nuclear graphite based on ring compression test. *J. Nucl. Mater.* **2018**, *511*, 134–140. [CrossRef]

42. Boyer, H.E. *Atlas of Stress-Strain Curves*; ASM International: Novelty, OH, USA, 1987.

43. Pirgazi, H.; Ghodrat, S.; Kestens, L. Three-dimensional EBSD characterization of thermo-mechanical fatigue crack morphology in compacted graphite iron. *Mater. Charact.* **2014**, *90*, 13–20. [CrossRef]

44. Norman, V.; Skoglund, P.; Moverare, J. Damage evolution in compacted graphite iron during thermomechanical fatigue testing. *Int. J. Cast Met. Res.* **2016**, *29*, 26–33. [CrossRef]

45. Hosdez, J.; Limodin, N.; Najjar, D.; Witz, J.F.; Charkaluk, E.; Osmond, P.; Forré, A.; Szmytka, F. Fatigue crack growth in compacted and spheroidal graphite cast irons. *Int. J. Fatigue* **2020**, *131*, 105319. [CrossRef]

46. Norman, V.; Calmunger, M. On the micro-and macroscopic elastoplastic deformation behaviour of cast iron when subjected to cyclic loading. *Int. J. Plast.* **2019**, *115*, 200–215. [CrossRef]

47. Hjelm, H.E. Yield surface for grey cast iron under biaxial stress. *J. Eng. Mater. Technol.* **1994**, *116*, 148–154. [CrossRef]

48. Seifert, T.; Riedel, H. Mechanism-based thermomechanical fatigue life prediction of cast iron. Part I: Models. *Int. J. Fatigue* **2010**, *32*, 1358–1367. [CrossRef]

49. Noguchi, T.; Shimizu, K. Accurate evaluation of the mechanical properties of grey cast iron. *Cast Met.* **1993**, *6*, 146–152. [CrossRef]

50. Rundman, K.B.; Iacoviello, F. Cast irons. In *Encyclopedia of Materials: Science and Technology*, 2nd ed.; Elsevier: Oxford, UK, 2001; pp. 1003–1010.