*Article* **A Phase Field Approach to Two-Dimensional Quasicrystals with Mixed Mode Cracks**

**Tong Li 1, Zhenting Yang 1, Chenghui Xu 2, Xinsheng Xu <sup>1</sup> and Zhenhuan Zhou 1,\***


**Abstract:** Quasicrystals (QCs) are representatives of a novel kind of material exhibiting a large number of remarkable specific properties. However, QCs are usually brittle, and crack propagation inevitably occurs in such materials. Therefore, it is of great significance to study the crack growth behaviors in QCs. In this work, the crack propagation of two-dimensional (2D) decagonal QCs is investigated by a fracture phase field method. In this method, a phase field variable is introduced to evaluate the damage of QCs near the crack. Thus, the crack topology is described by the phase field variable and its gradient. In this manner, it is unnecessary to track the crack tip, and therefore remeshing is avoided during the crack propagation. In the numerical examples, the crack propagation paths of 2D QCs are simulated by the proposed method, and the effects of the phason field on the crack growth behaviors of QCs are studied in detail. Furthermore, the interaction of the double cracks in QCs is also discussed.

**Keywords:** phase field model; decagonal quasicrystal; crack propagation; brittle fracture; mixed mode crack; finite element method

#### **1. Introduction**

Quasicrystals (QCs) are a new kind of material with perfect long-range order but no periodicity. Due to their unique atomic structures, QCs exhibit many excellent physical properties, such as high hardness, a low friction coefficient, and high resistance. Therefore, QCs are very promising materials for potential applications in corrosion-resistant coatings, hydrogen storage, photovoltaic solar cells, etc. [1]. However, QCs are usually brittle, and consequently, cracks, holes, and other defects will inevitably occur during daily use. If the crack propagates, the QC will fail, which may lead to a catastrophic accident. Therefore, the fracture analysis of QCs is of great importance, and it is very necessary to investigate crack propagation behaviors in cracked QCs.

Plenty of research work has been carried out on the fracture problems of QCs. Li et al. [2] extended the classical linear elasticity fracture mechanics to investigate a decagonal QC with a Griffith crack. The results indicated the phonon and phason stresses at the crack tip exhibit the well-known square root singularity, and the strain energy release rate was obtained. Zhou and Fan [3] developed the plane elasticity theory of two-dimensional (2D) octagonal QCs and obtained the exact analytic solution of a Mode I Griffith crack. Guo and Fan [4] studied the fracture problem of a Mode II Griffith crack in decagonal quasicrystals and obtained the corresponding stress intensity factors and strain energy release rate. Shen and Fan [5] calculated the stress intensity factors for an infinitely long strip of finite height containing two straight, semi-infinite collinear cracks. Li and Fan [6] obtained exact solutions for two semi-infinite collinear cracks in a strip of 1D hexagonal QCs. After that, they [7] further obtained an analytic solution for the elliptic notch problem

**Citation:** Li, T.; Yang, Z.; Xu, C.; Xu, X.; Zhou, Z. A Phase Field Approach to Two-Dimensional Quasicrystals with Mixed Mode Cracks. *Materials* **2023**, *16*, 3628. https://doi.org/ 10.3390/ma16103628

Academic Editors: Grzegorz Lesiuk and Dariusz Rozumek

Received: 31 March 2023 Revised: 5 May 2023 Accepted: 8 May 2023 Published: 9 May 2023

**Copyright:** © 2023 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/).

of the material in icosahedral QCs by using the complex variable method. The solution can be reduced to that of a Griffith crack problem. Fan et al. [8] presented linear, nonlinear, and dynamic fracture problems for different QCs. Li and Liu [9] obtained closed-form expressions for the elastic displacement and stress fields induced by a dislocation in icosahedral QCs. Sladek et al. [10] developed a meshless method based on the local Petrov-Galerkin approach for fracture analysis of decagonal QCs; both static and transient dynamic boundary value problems were considered. After that, they [11] present path-independent integrals for accurate evaluation of energy release and stress intensity factors in decagonal QCs. Wang et al. [12] obtained the phonon-phason coupling field in the half-space, which can be expressed in terms of elementary functions. These solutions could have applications in 3D contact and crack problems in QCs. Li et al. [13,14] derived solutions for elliptical crack and planar crack problems in 2D hexagonal QCs. Li and Shi [15] employed the method of potential function theory to solve plane defect problems originating from twodimensional decagonal QCs. Zhao et al. [16] derived the fundamental solutions for interface cracks in 1D hexagonal QC coatings under in-plane loads.

For crack propagation in QCs, some investigations have been conducted based on the atomistic model [17–29]. However, it is inconvenient to apply atomistic simulation to engineering. Therefore, Wang and Ricoeur [30] adopted the finite element method (FEM) to simulate the crack growth in 1D QCs and predicted the crack pattern for different boundary conditions. Fan et al. [31–39] investigated crack propagation based on the elastodynamic/hydrodynamic model.

From the above literature review, it is found that increasing attention has been paid to the fracture of QCs. However, most of them were concentrated on the derivation of exact solutions and determination of fracture parameters of cracked QCs, and the investigations on crack propagation were mainly based on the atomistic model [17–29]. The work on crack growth in QCs based on continuum mechanics is still insufficient. Additionally, in the traditional FEM, the crack topology is modeled by geometry. The remesh process is necessary when simulating crack propagation, which has a huge computational cost. In this paper, a phase field method is introduced to predict the crack propagation path in QCs. Unlike conventional discrete crack models (such as the FEM), the fracture phased field method employs diffusive cracks to avoid an explicit representation of kinematic discontinuities, and therefore the propagating cracks are tracked automatically without additional ad-hoc criteria in the classical Griffith fracture theory [40,41]. In this method, the fracture energy and degraded strain energy of QCs are formulated using the phase field variable. Subsequently, the total potential energy of QCs under the phase field framework is obtained, and the governing equations for the phase field model are derived by means of the Francfort-Marigo variational principle. The phase field evolution equation for QCs is constructed. Finally, the FEM is adopted to solve the phase field governing equations. The phase field variable and phonon/phason displacement of the entire model can be obtained, as well as the reaction force and the crack pattern.

This paper is organized as follows: The basic equations of 2D decagonal QCs are presented in Section 2.1. The phase field model for 2D QCs is formulated in Section 2.2. The finite element implementation for the phase field model is presented in Section 2.3. Several numerical examples are presented in Section 3. Conclusions are summarized in Section 4.

#### **2. Phase Field Method for 2D Decagonal QCs**

#### *2.1. The Basic Equations*

According to the elasticity of 2D decagonal QCs, the basic equations for the plane problem of QCs in the absence of body force and phason self-action are [42]:

$$\begin{cases} \frac{\partial \tau\_x}{\partial x} + \frac{\partial \tau\_{xy}}{\partial y} = 0\\ \frac{\partial \tau\_{yx}}{\partial x} + \frac{\partial \tau\_y}{\partial y} = 0' \begin{cases} \frac{\partial H\_x}{\partial x} + \frac{\partial H\_{xy}}{\partial y} = 0\\ \frac{\partial H\_{yx}}{\partial x} + \frac{\partial H\_y}{\partial y} = 0 \end{cases} \end{cases} \tag{1}$$

$$\begin{Bmatrix} \varepsilon\_x\\ \varepsilon\_y\\ \gamma\_{xy} \end{Bmatrix} = \begin{bmatrix} \frac{\partial}{\partial x} & & & \\ & \frac{\partial}{\partial y} & \\ & & \frac{\partial}{\partial x} \end{bmatrix} \begin{Bmatrix} u\_x\\ u\_y \end{Bmatrix} , \begin{Bmatrix} \omega\_x\\ \omega\_y\\ \omega\_{xy}\\ \omega\_{yx} \end{Bmatrix} = \begin{bmatrix} \frac{\partial}{\partial x} & & \\ & \frac{\partial}{\partial y} \\\ \frac{\partial}{\partial y} & & \\ & & \frac{\partial}{\partial x} \end{bmatrix} \begin{Bmatrix} w\_x\\ w\_y \end{Bmatrix} \tag{2}$$

$$
\begin{Bmatrix} \sigma\_x\\ \sigma\_y\\ \tau\_{xy} \end{Bmatrix} = \begin{bmatrix} \mathbf{C}\_{11} & \mathbf{C}\_{12} & & \\ \mathbf{C}\_{12} & \mathbf{C}\_{22} & & \\ & & \mathbf{C}\_{66} \end{bmatrix} \begin{Bmatrix} \varepsilon\_x\\ \varepsilon\_y\\ \gamma\_{xy} \end{Bmatrix} + \begin{bmatrix} R\_1 & R\_1 & R\_2 & -R\_2\\ -R\_1 & -R\_1 & -R\_2 & R\_2\\ R\_2 & R\_2 & -R\_1 & R\_1 \end{bmatrix} \begin{Bmatrix} \omega\_x\\ \omega\_y\\ \omega\_{xy}\\ \omega\_{yx} \end{Bmatrix} \tag{3}
$$

$$\begin{Bmatrix} H\_{\mathcal{X}} \\ H\_{\mathcal{Y}} \\ H\_{\mathcal{X}\mathbf{y}} \\ H\_{\mathcal{Y}\mathbf{x}} \end{Bmatrix} = \begin{bmatrix} R\_1 & -R\_1 & R\_2 \\ R\_1 & -R\_1 & R\_2 \\ R\_2 & -R\_2 & -R\_1 \\ -R\_2 & R\_2 & R\_1 \end{bmatrix} \begin{Bmatrix} \varepsilon\_{\mathcal{X}} \\ \varepsilon\_{\mathcal{Y}} \\ \gamma\_{\mathcal{X}\mathbf{y}} \end{Bmatrix} + \begin{bmatrix} K\_1 & K\_2 \\ K\_2 & K\_1 \\ & K\_1 & -K\_2 \\ & & -K\_2 & K\_1 \end{bmatrix} \begin{Bmatrix} \omega\_{\mathcal{X}} \\ \omega\_{\mathcal{Y}} \\ \omega\_{\mathcal{Y}\mathbf{x}} \\ \omega\_{\mathcal{Y}\mathbf{x}} \end{Bmatrix} \tag{4}$$

where *ux* and *uy* are the phonon displacements; *wx* and *wy* are the phason displacements; *εx*, *εy*, and *γxy* are the phonon strains; *ωx*, *ωy*, *ωxy*, and *ωyx* are the phason strains; *σx*, *σy*, and *σxy* are the phonon stresses; *Hx*, *Hy*, *Hxy*, and *Hyx* are the phason stresses; *C*11, *C*12, *C*22, and *C*<sup>66</sup> are the phonon moduli; *K*<sup>1</sup> and *K*<sup>2</sup> are the phason moduli; and *R*<sup>1</sup> and *R*<sup>2</sup> are the phonon-phason coupling coefficients.

The strain energy of QCs is:

$$\mathrm{II}\_{\mathfrak{s}} = \int\_{\varOmega} \psi(\mathfrak{s}, \,\,\mathfrak{w}) \,\mathrm{d}V = \int\_{\varOmega} \frac{1}{2} \mathfrak{s}^{\mathrm{T}} \mathfrak{o} + \frac{1}{2} \mathfrak{w}^{\mathrm{T}} \mathbf{H} \,\mathrm{d}V \tag{5}$$

where σ = - *<sup>σ</sup>x*, *<sup>σ</sup>y*, *<sup>σ</sup>xy*<sup>T</sup> and **<sup>H</sup>** = - *Hx*, *Hy*, *Hxy*, *Hyx*<sup>T</sup> are the stress vectors; ε = - *<sup>ε</sup>x*, *<sup>ε</sup>y*, *<sup>γ</sup>xy*<sup>T</sup> and <sup>ε</sup> = - *<sup>ω</sup>x*, *<sup>ω</sup>y*, *<sup>ω</sup>xy*, *<sup>ω</sup>yx*<sup>T</sup> are the strain vectors; and *<sup>ψ</sup>*(ε, <sup>ω</sup>) is the strain energy density.

#### *2.2. Phase Field Method*

Consider a 1D QC strip with a centered crack, as shown in Figure 1. The fracture energy is:

$$\mathrm{II}\_{\mathfrak{c}} = \int\_{\partial \Omega\_{\mathfrak{c}}} \mathrm{G}\_{\mathfrak{c}} \, \mathrm{d}S = \mathrm{G}\_{\mathfrak{c}} A\_{\mathfrak{c}} \tag{6}$$

where *Gc* is the critical energy release rate (CERR); *∂*Ω*<sup>c</sup>* is the crack surface; and *Ac* is the cross-sectional area of the strip.

**Figure 1.** A QC strip with a centered crack.

In a phase field model, the crack is supposed to be surrounded by a diffusive degraded zone, and a phase field variable *d* is introduced to describe the damage in the diffusive degraded zone:

$$d = e^{-\frac{|x|}{\tau\_c}}\tag{7}$$

where *lc* is an internal length scale that controls the width of the diffusive zone. Here, it is noted that, when *x* approaches infinity (±∞), *d* converges to zero, which indicates that the material is intact; when *x* is zero, *d* equals one, which indicates that the material

is totally broken and the crack surface is produced. Define a crack's surface density as functional [43]:

$$\gamma \begin{pmatrix} d, \ d' \end{pmatrix} = \frac{1}{2l\_c} \left( d^2 + l\_c^2 d'^2 \right) \tag{8}$$

The crack topology can then be described by the phase field variable.

Substituting Equation (8) into Equation (6) yields the fracture energy in the phase field model:

$$
\Pi\_{\mathfrak{c}} = \int\_{\Omega} \mathbb{G}\_{\mathfrak{c}} \gamma \begin{pmatrix} d\_{\mathfrak{c}} \ d' \end{pmatrix} \,\mathrm{d}V \tag{9}
$$

where Ω is the domain of the overall model. As observed in Equation (9), the integral over the crack path is transformed into the integral over the model. The crack topology is implicitly expressed in the framework of the phase field method. Therefore, it is unnecessary to track the crack's path while it is propagating.

Similarly, the 2D crack surface density functional can be defined as:

$$\gamma(d,\ \nabla d) = \frac{1}{2l\_c} \left( d^2 + l\_c^2 |\nabla d|^2 \right) \tag{10}$$

Therefore, the fracture energy of QCs in a plane problem is:

$$\Pi\_{\mathfrak{c}} = \int\_{\Omega} \mathcal{G}\_{\mathfrak{c}} \gamma(d\_{\succ} \,\, \nabla d) \,\, \mathrm{d}V \tag{11}$$

In the phase field model, the material at the crack is assumed to be softened. Therefore, a degradation function is introduced to evaluate the damage to the material in the diffusive degraded zone [43]:

$$
\omega(d) = \left(1 - d\right)^2 + k \tag{12}
$$

where *k* is a small positive number to ensure the non-singularity of the matrix. It can be observed that *ω*(*d*) satisfies *ω*(0) = 1 and *ω*(1) = 0.

The strain energy for QCs is modified as follows:

$$
\Pi\_5 = \int\_{\Omega} \omega(d)\psi(\mathfrak{e}, \,\,\mathfrak{w}) \,\,\mathrm{d}V \tag{13}
$$

It should be pointed out that the crack will not propagate if the crack face is under compression. Therefore, the strain energy density *ψ* should be decomposed into a tensile part and a compressed part, and the energy degradation only occurs on the tensile part [43]:

$$\Pi\_{\mathbb{S}} = \int\_{\Omega} \omega(d)\psi\_{+} + \psi\_{-} \,\mathrm{d}V \tag{14}$$

where *ψ*<sup>+</sup> = *λε*<sup>1</sup> + *ε*<sup>2</sup> + *ε*3 2 <sup>+</sup>/2 + *μ ε*1 2 <sup>+</sup> + *ε*2 2 <sup>+</sup> + *ε*3 2 + and *ψ*<sup>−</sup> = *ψ* − *ψ*<sup>+</sup> are the tensile and compressed strain energy densities, respectively; <sup>+</sup> is defined as *x*<sup>+</sup> = (|*x*| + *x*)/2; and *λ* and *μ* are the Lamé constants.

The potential energy of the overall model contains the strain energy, fracture energy, and potential energy of the external force:

$$\begin{array}{lcl} \Pi\_{\mathcal{V}} &= \Pi\_{\mathfrak{s}} + \Pi\_{\mathfrak{c}} - \Pi\_{\mathfrak{c}} \\ &= \int\_{\Omega} \omega(d) \psi\_{+} + \psi\_{-} \operatorname{d}V + \int\_{\Omega} \operatorname{G}\_{\mathfrak{c}} \gamma(d, \,\nabla d) \, \operatorname{d}V \\ &- \int\_{\partial\Omega} \mathbf{u}^{\mathsf{T}} \mathbf{t}^{\mathsf{u}} \operatorname{d}S - \int\_{\partial\Omega} \mathbf{w}^{\mathsf{T}} \mathbf{t}^{\mathsf{w}} \operatorname{d}S \end{array} \tag{15}$$

where **u** = - *ux*, *uy* T ; **w** = - *wx*, *wy* T ; and **t** *<sup>u</sup>* and **t** *<sup>w</sup>* are the external forces in the phonon and phason fields, respectively.

The Francfort-Marigo variational principle states that the real displacement **u** and the phase field variable *d* will minimize the potential energy.

*δ*II*<sup>p</sup>* = *<sup>∂</sup>*Ω*<sup>u</sup> <sup>ω</sup>*(*d*)*σ*<sup>+</sup> *ij niδuj* + *σ*<sup>−</sup> *ij niδuj* <sup>d</sup>*<sup>S</sup>* <sup>+</sup> *<sup>∂</sup>*Ω*<sup>w</sup> <sup>ω</sup>*(*d*)*H*<sup>+</sup> *ij niδwj* + *H*<sup>−</sup> *ij niδwj* d*S* − Ω *∂ ω*(*d*)*σ*<sup>+</sup> *ij <sup>∂</sup>xi <sup>δ</sup>uj* + *∂σ*− *ij <sup>∂</sup>xi <sup>δ</sup>uj* <sup>d</sup>*<sup>V</sup>* <sup>−</sup> Ω *∂ ω*(*d*)*H*<sup>+</sup> *ij <sup>∂</sup>xi <sup>δ</sup>wj* + *∂H*− *ij <sup>∂</sup>xi δwj* d*V* <sup>Ω</sup> *ω*- (*d*)*ψ*+*δd* d*V* + Ω *Gc lc <sup>d</sup>δ<sup>d</sup>* <sup>d</sup>*<sup>V</sup>* <sup>+</sup> *<sup>∂</sup>*<sup>Ω</sup> *Gclc <sup>∂</sup><sup>d</sup> ∂xi niδ<sup>d</sup>* <sup>d</sup>*<sup>S</sup>* <sup>−</sup> <sup>Ω</sup> *Gclc ∂*2*d ∂x*<sup>2</sup> *i δd* d*V* − *<sup>∂</sup>*<sup>Ω</sup> *δujt u <sup>j</sup>* <sup>d</sup>*<sup>S</sup>* <sup>−</sup> *<sup>∂</sup>*<sup>Ω</sup> *δwjt w <sup>j</sup>* d*S* = 0 (16)

where *σ*<sup>+</sup> *ij* and *<sup>H</sup>*<sup>+</sup> *ij* are the stresses induced by stretch, while *σ*<sup>−</sup> *ij* and *H*<sup>−</sup> *ij* are induced by compression. Equation (16) is valid for arbitrary *δui*, *δwi*, and *δd*. Hence, the governing equations for the phase field model of QCs are:

$$\nabla \left[ \omega(d) \sigma\_{ij}^{+} + \sigma\_{ij}^{-} \right] = 0 \text{, } \Omega \tag{17}$$

$$\nabla \left[ \omega(d) H\_{ij}^+ + \partial H\_{ij}^- \right] = 0 \text{, } \Omega \tag{18}$$

$$(2d - 2)\psi\_{+} + \frac{G\_{c}}{l\_{c}}d - G\_{c}l\_{c}\Delta d = 0, \Omega \tag{19}$$

$$\left[\omega(d)\sigma\_{ij}^{+} + \sigma\_{ij}^{-}\right]n\_{i} = t\_{j}^{\mu}\,\,\,\partial\Omega\tag{20}$$

$$\left[\omega(d)H\_{ij}^+ + H\_{ij}^-\right]n\_i = t\_j^w,\ \partial\Omega\tag{21}$$

$$d\_{,i}n\_i = 0,\ \partial\Omega\tag{22}$$

It should be mentioned that crack growth is an irreversible process. Therefore, Equation (19) should be modified by considering the history of the load [44]:

$$(2d - 2)H + \frac{G\_c}{l\_c}d - G\_c l\_c \Delta d = 0\tag{23}$$

where *H* = max [0, *t*] {*ψ*+} is a history variable that is the maximum strain energy during the crack propagation. This history variable ensures that the crack face does not close under compression. Equation (23) is the evolution law of the crack phase field for QCs. Cracks grow only if this equation is valid.

#### *2.3. Finite Element Implementation*

Due to the strong nonlinearity of the governing equations in the phase field model, the FEM is often adopted to solve the problem. In the FEM, the phonon/phason displacements and the phase field variable are approximated by the shape functions:

$$\left\{ \boldsymbol{\mu}\_{\mathcal{X}\_{\mathcal{I}}} \; \boldsymbol{\mu}\_{\mathcal{Y}} \right\}^{\mathsf{T}} = \sum\_{i=1}^{4} \mathbf{N}\_{i}^{\mathsf{u}} \mathbf{u}\_{i} \tag{24}$$

$$\left\{w\_{\mathbf{x}} \mid w\_{\mathbf{y}}\right\}^{\mathrm{T}} = \sum\_{i=1}^{4} \mathbf{N}\_{i}^{w} \mathbf{w}\_{i} \tag{25}$$

$$d = \sum\_{i=1}^{4} N\_i d\_i \tag{26}$$

where *Ni* is the shape function; **N***<sup>u</sup> <sup>i</sup>* = **<sup>N</sup>***<sup>w</sup> <sup>i</sup>* = diag(*Ni*, *Ni*); and **u***i*, **w***i*, and *d* are the nodal phonon and phason displacements and phase field variables, respectively.

Therefore, the phonon and phason strains and the gradient of the phase field variable are, respectively:

$$\left\{ \begin{array}{cccc} \varepsilon\_{X} & \varepsilon\_{Y} & \gamma\_{XY} \end{array} \right\}^{\mathrm{T}} = \sum\_{i=1}^{4} \mathbf{B}\_{i}^{\mu} \mathbf{u}\_{i} \tag{27}$$

$$\left\{\omega\_{\mathbf{x}} \mid \omega\_{\mathbf{y}} \mid \omega\_{\mathbf{x}\mathbf{y}} \mid \omega\_{\mathbf{y}\mathbf{x}}\right\}^{\mathrm{T}} = \sum\_{i=1}^{4} \mathbf{B}\_{i}^{\mathrm{w}} \mathbf{w}\_{i} \tag{28}$$

$$\nabla d = \sum\_{i=1}^{4} \mathbf{B}\_i^d \mathbf{d}\_i \tag{29}$$

where**B***<sup>u</sup> <sup>i</sup>* = *∂Ni*/*∂x ∂Ni*/*∂y ∂Ni*/*∂y ∂Ni*/*∂x* T , **B***<sup>w</sup> <sup>i</sup>* = *∂Ni*/*∂x ∂Ni*/*∂y ∂Ni*/*∂y ∂Ni*/*∂x* T , and **B***<sup>d</sup> <sup>i</sup>* <sup>=</sup> *∂Ni*/*∂x ∂Ni*/*∂y* T . Substituting Equations (24)–(29) into Equation (15) yields the residuals:

$$\mathbf{r}^{a} = \int\_{\partial\Omega} \mathbf{N}^{\mathrm{T}} \mathbf{t} \, \mathrm{d}S - \int\_{\Omega} (1 - d)^{2} \mathbf{B}^{\mathrm{T}} \mathbf{D}\_{\mathrm{QC}} \mathbf{B} \mathbf{a} \, \mathrm{d}V = 0 \tag{30}$$

$$\mathbf{r}^d = -\int\_{\Omega} 2(d-1) H^d \left(\mathbf{N}^d\right)^\mathrm{T} \mathrm{d}V - \int\_{\Omega} \frac{\mathrm{G}\_\mathbf{c}}{\mathcal{E}\_0} \left[\frac{2d}{l\_\mathbf{c}} \left(\mathbf{N}^d\right)^\mathrm{T} + 2l\_\mathbf{c} \left(\mathbf{B}^d\right)^\mathrm{T} \nabla d\right] \mathrm{d}V = 0 \tag{31}$$

where **a** = - *ux*1, *uy*1, ..., *ux*4, *uy*4, *wx*1, *wy*1, ..., *wx*4, *wy*<sup>4</sup> T ; **N** = **<sup>N</sup>**<sup>0</sup> **<sup>0</sup>**2×<sup>8</sup> **0**2×<sup>8</sup> **N**<sup>0</sup> is the shape function matrix where **N**<sup>0</sup> = *N*<sup>1</sup> 0 *N*<sup>2</sup> 0 *N*<sup>3</sup> 0 *N*<sup>4</sup> 0 0 *N*<sup>1</sup> 0 *N*<sup>2</sup> 0 *N*<sup>3</sup> 0 *N*<sup>4</sup> ; **B** = **LN** is ⎡ *∂ <sup>∂</sup><sup>x</sup>* <sup>0</sup> *<sup>∂</sup> <sup>∂</sup><sup>y</sup>* 0000 ⎤ T

$$\text{the strain matrix where } \mathbf{L} = \begin{vmatrix} \mathbf{u} & \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 & 0 & 0 & 0\\ 0 & \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & \frac{\partial}{\partial x} & 0 & \frac{\partial}{\partial y} & 0\\ 0 & 0 & 0 & 0 & \frac{\partial}{\partial y} & 0 & \frac{\partial}{\partial x} \end{vmatrix}; \text{and } \mathbf{t} = \left\{ t\_{\mathbf{x}'}^{\mathbf{u}} \, t\_{\mathbf{y}'}^{\mathbf{u}} \, t\_{\mathbf{x}'}^{\mathbf{w}} \, t\_{\mathbf{y}}^{\mathbf{w}} \right\}^{\mathbf{T}}$$

is the load vector.

Equations (30) and (31) can be solved by this iteration method:

$$
\begin{Bmatrix} \mathbf{a}\_{n+1} \\ \mathbf{d}\_{n+1} \end{Bmatrix} = \begin{Bmatrix} \mathbf{a}\_{n} \\ \mathbf{d}\_{n} \end{Bmatrix} + \begin{bmatrix} \mathbf{K}\_{n}^{aa} & \mathbf{K}\_{n}^{ad} \\ \mathbf{K}\_{n}^{da} & \mathbf{K}\_{n}^{dd} \end{bmatrix}^{-1} \begin{Bmatrix} \mathbf{r}\_{n}^{a} \\ \mathbf{r}\_{n}^{d} \end{Bmatrix} \tag{32}$$

where **K***aa* = <sup>Ω</sup> *<sup>ω</sup>*(*d*)**B**T**DB** <sup>d</sup>*V*, **<sup>K</sup>***ad* = **K***da*<sup>T</sup> = <sup>Ω</sup> <sup>2</sup>(*<sup>d</sup>* <sup>−</sup> <sup>1</sup>)**B**Tσ**<sup>N</sup>** <sup>d</sup>*V*, and **K***dd* = Ω 2*H* + *Gc l* **N***dT***N***ddV* + <sup>Ω</sup> *lGc***B***dT***B***ddV*.

Due to the high nonlinearity of Equation (32), there is a convergence problem using the classic Newton iteration. A highly robust staggered algorithm is usually adopted to solve Equation (32) [41,44]. In this algorithm, one of the two unknowns (displacement and phase field variable) is assumed to be unchanged while solving the other unknown during one iteration, which yields:

$$\mathbf{a}\_{n+1} = \mathbf{a}\_n + \left(\mathbf{K}\_n^{aa}\right)^{-1} \mathbf{r}\_n^a \tag{33}$$

$$\mathbf{d}\_{n+1} = \mathbf{d}\_n + \left(\mathbf{K}\_n^{dd}\right)^{-1} \mathbf{r}\_n^d \tag{34}$$

Finally, the nodal phonon/phason displacements and the phase field variables can be calculated by Equations (33) and (34), respectively.

#### **3. Numerical Results**

In this section, a few numerical examples are presented to illustrate the application of the phase field method to the fracture of 2D decagonal QCs. The material parameters are tabulated in Table 1 [45–47].

**Table 1.** Material properties of 2D decagonal QCs.


According to Fan's work [42], the expression of CERR is:

$$G\_{\mathbb{C}} = \frac{\lambda \left(K\_1 + K\_2\right) + 2\left(R\_1^2 + R\_2^2\right)}{8\left(\lambda + M\right)c} K\_{IC}^2 \tag{35}$$

where *M* = (*C*<sup>11</sup> − *C*12)/2 and *c* = *M*(*K*<sup>1</sup> + *K*2) − 2 *R*2 <sup>1</sup> + *<sup>R</sup>*<sup>2</sup> 2 . The fracture toughness of Al-Ni-Co QCs is *KIC* <sup>=</sup> 1MPa√m. Therefore, the CERR in this paper is selected as *Gc* = 5.56 × <sup>10</sup>−<sup>4</sup> N/mm, according to the material constants in Table 1. Although no experimental method has been reported to apply a constant phason displacement on the surface of QCs, some investigations [48,49] reveal that some ways can cause the disorder of the phason field. Therefore, different phason displacement loads are considered in the following calculation to investigate the effect of the phason field on the fracture behavior of QCs.

#### *3.1. The Rectangular QCs with Edge Crack*

A rectangular QC model is considered, as shown in Figure 2. The geometry of the model is width *W*, height *L*, and crack length *a*. The upper edge is constrained to have the same displacement in phonon and phason fields, and the lower edge is only constrained in the vertical direction. A concentrated phonon force *P<sup>σ</sup>* is applied to the upper edge. The angle between *Pσ* and the horizon is *ϕ*.

**Figure 2.** A rectangular QC with an edge crack.

At the current stage, the crack growth of 2D QCs has not been reported in the open literature. To demonstrate the accuracy of the present method, an elastic material is selected by degenerating all the material constants in phason field, i.e., *K*<sup>1</sup> = *K*<sup>2</sup> = *R*<sup>1</sup> = *R*<sup>2</sup> = 0. The geometrical parameters are taken as *a* = 0.5 mm and *W* = *L* = 1 mm. The elastic material constants come from Ref. [44]. The variation of the vertical concentrated force *P<sup>σ</sup>* (*ϕ* = *π*/2) versus the displacement *u* is plotted in Figure 3. As observed, the present results are in excellent agreement with the reference data from Ref. [44]. Furthermore, it is noted that the peak force *P<sup>σ</sup>* shows an increasing trend as the length scale *lc* decreases, which indicates the crack grows slowly as the length scale declines. In the previous study, to ensure the high resolution of the crack topology, the minimum size of the element was required to meet the condition of *h* < *lc*/2 [43]. Therefore, in the following calculations, the length scale is selected as *lc* = 1 mm, and the element size is selected as *h* < *lc*/5.

**Figure 3.** Variation of the force *Pσ* versus the displacement *u* [44].

Subsequently, the effect of the phason field on the crack growth is studied in Figure 2 by considering an initial phason displacement *w*<sup>0</sup> on the upper edge of the model. The parameters are selected as *a* = 50 mm and *W* = *L* = 100 mm. The variation of the applied force versus the displacement for different *w*<sup>0</sup> values is illustrated in Figure 4. The results show that the peak force decreases as the initial phason displacement increases, which indicates that the initial stretch phason load will undermine the strength of the model. In addition, as the initial phason load gets larger, the crack grows slower. The crack pattern with *<sup>w</sup>*<sup>0</sup> = <sup>6</sup> × <sup>10</sup>−<sup>4</sup> mm is plotted in Figure 5.

**Figure 4.** Force-displacement relations for different initial phason loads *w*0.

**Figure 5.** The crack pattern with *<sup>w</sup>*<sup>0</sup> <sup>=</sup> <sup>6</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> mm.

Finally, the QC model subjected to a concentrated shear load at its upper edge is considered. The variation of the shear phonon force *P<sup>σ</sup>* (*ϕ* = *π*) versus the horizontal displacement at the upper edge is illustrated in Figure 6. Two peak forces were observed, and they decrease as the initial vertical phason displacement *w*<sup>0</sup> increases. The crack patterns with different initial phason displacements are shown in Figure 7. Clearly, the increasing initial phason displacement will lead to a significant crack deflection. Therefore, it is concluded that the phason load is a key influencing factor for the peak force and crack patterns under the applied shear loads.

**Figure 6.** The shear force-displacement relations for different initial tensile phason loads.

**Figure 7.** The crack pattern with shear load for different initial phason loads (edge crack): (**a**) *<sup>w</sup>*<sup>0</sup> <sup>=</sup> 0 mm; (**b**) *<sup>w</sup>*<sup>0</sup> <sup>=</sup> <sup>4</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> mm; and (**c**) *<sup>w</sup>*<sup>0</sup> <sup>=</sup> <sup>8</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> mm.

#### *3.2. The Rectangular QCs with an Internal Crack*

In this example, a rectangular QC with an internal crack is considered in Figure 8. The length of the crack is *a*, and the angle between the crack and the horizon is *θ*. The middle of the crack is centered in the model. The boundary conditions of the QC are the same as those in Section 3.1. The parameters are selected as *a* = 25 mm and *W* = *L* = 100 mm. The QC model is subjected to a vertical initial phason displacement *w*<sup>0</sup> at the upper edge.

**Figure 8.** A rectangular QC with an internal crack.

Firstly, the QC model is subjected to a tensile load in the phonon field (*ϕ* = *π*/2). The peaks of the phonon forces for different angles *θ* are shown in Figure 9. It can be observed that the peak force monotonously increases with the increase in *θ*, while it shows an opposite trend as the initial phason displacement increases. The force-displacement relation for different angles *θ* is shown in Figure 10. Clearly, the slope of the force-displacement curve increases as the angle *θ* increases. The crack patterns with different angles are plotted in Figure 11. As observed, the crack grows along a straight line to the edge of the model when subjected to a tensile load.

**Figure 9.** The peak phonon forces for different angles *θ* (edge crack).

**Figure 10.** The force-displacement relation for different angles *θ*.

**Figure 11.** The crack patterns with *<sup>w</sup>*<sup>0</sup> <sup>=</sup> <sup>6</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> mm for different angles *<sup>θ</sup>*: (**a**) *<sup>θ</sup>* <sup>=</sup> 0; (**b**) *<sup>θ</sup>* <sup>=</sup> *<sup>π</sup>*/6; and (**c**) *θ* = *π*/3.

Subsequently, the model is subjected to shear loads in the phonon field at the upper edge. The peak phonon forces for different angles *θ* are illustrated in Figure 12. As observed, the peak force first increases and then decreases as the angle *θ* increases. The crack patterns are plotted in Figure 13. Similar to the observation in Figure 7, the increasing initial phason displacement has a significant influence on the crack propagation path.

**Figure 12.** The peak forces for different angles *θ* (internal crack).

**Figure 13.** The crack pattern with shear load for different initial phason loads (internal crack): (**a**) *<sup>w</sup>*<sup>0</sup> <sup>=</sup> <sup>8</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> mm; (**b**) *<sup>w</sup>*<sup>0</sup> <sup>=</sup> 1.2 <sup>×</sup> <sup>10</sup>−<sup>3</sup> mm; and (**c**) *<sup>w</sup>*<sup>0</sup> <sup>=</sup> 1.6 <sup>×</sup> <sup>10</sup>−<sup>3</sup> mm.

#### *3.3. The Rectangular QCs with Double Cracks*

As the last example, the QC model with double cracks is investigated in Figure 14. The parameters are selected as *a*<sup>1</sup> = *a*<sup>2</sup> = 25 mm and *W* = *L* = 100 mm. The left crack is horizontal. The angle between the right crack and the horizon is *θ*. The lengths of the left and right cracks are *a*<sup>1</sup> and *a*2, respectively. Point *B* is located at the center of the right crack. The distance between the right edge of the model and point *B* is Δ<sup>1</sup> = 27.5 mm. The distance between the left crack and point *B* is Δ2. The model is subjected to a tensile phonon/phason displacements. The peak forces for different angles *θ* with Δ<sup>2</sup> = 32.5 mm are plotted in Figure 15. It can be found that, as *θ* increases, the peak force monotonously increases for *w*<sup>0</sup> = 0 mm and *w*<sup>0</sup> = 0.0004 mm, while it increases first and then decreases for *w*<sup>0</sup> = 0.0008 mm. The crack patterns for different angles *θ* and distances Δ<sup>2</sup> are illustrated in Figures 16 and 17. As observed, the interaction of two cracks has a big contribution to their crack propagation path. The two cracks eventually merge together as the applied load increases.

**Figure 14.** A rectangular QC with double internal cracks.

**Figure 15.** The peak forces for different angles *θ* (double cracks).

**Figure 16.** The crack pattern with Δ<sup>2</sup> = 32.5 mm for different angles (double cracks): (**a**) *θ* = *π*/12; (**b**) *θ* = *π*/4; (**c**) *θ* = 5*π*/12.

**Figure 17.** The crack pattern with *θ* = *π*/4 for different distances Δ2: (**a**) Δ<sup>2</sup> = 32.5 mm; (**b**) Δ<sup>2</sup> = 22.5 mm; and (**c**) Δ<sup>2</sup> = 12.5 mm.

#### **4. Conclusions**

In this paper, a phase field model is developed to predict crack propagation in 2D decagonal QCs. The contribution of the phason field to the potential energy is considered, and the evolution law of the crack phase field for QCs is established. Therefore, the crack topology in QCs is described by the phase field variable, which can be solved by the FEM. In this manner, the crack propagation in QCs can be accurately simulated without remeshing, and the evolution of the crack can be vividly illustrated. Numerical examples illustrate that the proposed model can predict accurate results for the crack propagation of QCs, and the phason field has a big contribution to both the force-displacement relation and the crack pattern.

**Author Contributions:** Conceptualization, Z.Z. and Z.Y.; methodology, Z.Y. and T.L.; data curation, T.L.; writing—preparation of original draft, T.L., Z.Y. and C.X.; writing—proofreading and editing, Z.Z. and X.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Fundamental Research Funds for the Central Universities, [grant number DUT21LK35, DUT22LK16]; State Key Laboratory of Structural Analysis, Optimization and CAE Software for Industrial Equipment, [grant number GZ22109].

**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 authors.

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

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

#### MDPI

St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Materials* Editorial Office E-mail: materials@mdpi.com www.mdpi.com/journal/materials

Academic Open Access Publishing

www.mdpi.com ISBN 978-3-0365-8125-5