## Abstract

In this paper, a thermal–mechanical augmented finite-element method (TM-AFEM) has been proposed, implemented and validated for steady-state and transient, coupled thermal–mechanical analyses of complex materials with explicit consideration of arbitrary evolving cracks. The method permits the derivation of explicit, fully condensed thermal–mechanical equilibrium equations which are of mathematical exactness in the piece-wise linear sense. The method has been implemented with a 4-node quadrilateral two-dimensional (2D) element and a 4-node tetrahedron three-dimensional (3D) element. It has been demonstrated, through several numerical examples that the new TM-AFEM can provide significantly improved numerical accuracy and efficiency when dealing with crack propagation problems in 2D and 3D solids under coupled thermal–mechanical loading conditions.

This article is part of the themed issue ‘Multiscale modelling of the structural integrity of composite materials’.

## 1. Introduction

Owing to their ultra-high temperature-resisting capabilities, ceramic matrix composites (CMCs) have been widely used in applications such as nuclear fusion reactors, gas turbines and aircraft engines, where the high heat flux environment is hostile and damaging for typical engineering materials. For example, integral textile CMCs have been fabricated into components, such as rocket nozzles, turbine engine combustor liners, thermal protection systems and hypersonic flow path components. However, CMCs are prone to cracking due to large thermal stresses [1,2]. Thermally induced cracking, especially those cracks generated in-service are extremely dangerous, because they provide pathways for further damaging processes such as oxidation and vapour-assisted corrosion. Failure of such materials in any of these applications can be catastrophic and highly undesirable from an engineering point of view.

High-fidelity thermal–mechanical analyses to high-temperature CMCs with consideration of arbitrary cracking are very challenging, because the highly heterogeneous nature and processing defects make it nearly impossible to know the crack locations *a priori*. Further, the discontinuous nature of localized cracks cannot be directly modelled by standard finite-element method (FEM). A major breakthrough enabling explicit crack modelling within an FEM framework is due to the seminal theory of partition of unity (PoU), which led to the development of the generalized FEM [3,4]. During the same time period, Belytschko and co-workers [5,6] independently developed the extended FEM (X-FEM), which permits arbitrary discontinuities to arise in a simulation without the need for remeshing. The PoU theory allows for the addition of *a priori* solution of a boundary-value problem into the approximation spaces through numerical enrichment schemes. Thus, an FE mesh does not have to conform to structural boundaries or crack surfaces [5–12].

In composite failure analysis, another class of methods, the phantom node method (PNM) [13–18] and the augmented finite-element method (A-FEM) [19–21], are more widely used. These methods are derived from the seminal work of Hansbo & Hansbo [22] and further developed in [13,14]. Other emerging methods for composite failure analyses include the regularized FEM (RxFEM) of Iarve *et al.* [23,24] and the continuum–decohesion FEM by Waas and co-workers [25].

In composite materials, correct and efficient treatment of crack coalescence and bifurcation is critical for simulating the complex, coupled multi-crack damage evolution. The original version of the A-FEM and the RxFEM have been successfully used to study multiple arbitrary cracking problems in laminated composites [21,26–28]. In these studies, the coupling between intra-ply cracks and delaminations were done through proper designation of different types of elements for different failure processes, i.e. A-FE or RxFE for intra-ply cracking and augmented cohesive zone elements for interface delamination. For more complex interactive multi-crack systems typically observed in CMCs (crack merging, bifurcation crossing, etc.), these methods are much less capable because the need of using multiple copies of nodes or d.f. makes them rather inflexible in handling the nonlinear interactions. More recently, a new version of the A-FEM, which allows for explicit inclusion of descriptions of cohesive cracks, without the need to change mesh adaptively or to add extra d.f. dynamically as the cracks propagate, has been developed and demonstrated with much improved numerical performance for a variety of fracture problems [2,29–31]. It has been shown repeatedly that the new A-FEM is able to account for multiple, arbitrary cracks and crack interactions in solids with much improved numerical efficiency [29,30] (2–3 orders of magnitude in many cases).

Motivated by the success in the new two-dimensional (2D) A-FEM for arbitrary cracking in solids under mechanical loading, this paper is aimed to further empower the A-FEM formulation with transient and steady-state temperature analysis capabilities, so that it can deal with multiple cohesive crack formation and their direct interaction under coupled thermal–mechanical loading environment. Such capabilities are essential for high-fidelity structural safety analyses for complex CMCs under extreme thermal–mechanical loadings.

The paper is organized as follows: in §2, we introduce the statement of the coupled thermal mechanical problem in strong and weak form, and the thermal and mechanical cohesive laws used in this paper to describe the kinematics of cohesive cracks; in §3, the detailed augmentation processes for a 2D 4-node quadrilateral element and a three-dimensional (3D) 4-node tetrahedron element are described with all necessary equations explicitly derived; §4 provides several numerical examples to validate and to demonstrate the numerical capabilities of the 2D and 3D TM-AFEs and finally, the paper is concluded in §5 with highlights.

## 2. Problem statement

### (a) Strong and weak statements for a thermo-elastic domain with a discontinuity

Suppose a domain *Ω* as shown in figure 1 is divided by a discontinuity, . The resulting subdomains can also be either the same or of different materials such that *Ω*=*Ω*^{+}∪*Ω*^{−}. In general, this domain is mechanically subjected to prescribed displacement and surface traction on the boundaries and respectively. Likewise, thermally, temperature and heat flux are prescribed on boundaries and For such a problem, the strong forms of the mechanical and thermal equilibrium equations under both mechanical and thermal boundary conditions are
2.1where *σ*^{+/−} is the Cauchy stress tensor, **q**^{+/−} is the heat flux, *ρ*^{+/−} is the density, *c*^{+/−} is the specific heat and *t* is the physical time. The superscripts ‘+/−’ denotes quantities in the subdomain *Ω*^{+} or *Ω*^{−}, respectively. The last two expressions in equation (2.1) come from the condition of stress and heat flux continuity across the discontinuity. In general, the traction and heat flux along the discontinuity, and , are functions of crack displacement and temperature jump across the discontinuity, i.e.
2.2Assuming small strain and displacement in the domain surrounding the discontinuity, and that the heat flow within the subdomains follows Fourier’s Law, the constitutive, kinematic and heat flow equations for the two subdomains are
2.3where **C**^{+/−} is the material stiffness matrix tensor, **k**^{+/−} is the thermal conductivity tensor, *ε*^{+/−}_{T} is the total strain and the thermal strain. ∇ is the gradient operator. Δ*θ*^{+/−}=*θ*^{+/−}−*θ*_{0} is the temperature change from a stress-free reference temperature *θ*_{0} to the current temperature *θ*^{+/−}. *α*^{+/−} is the thermal expansion coefficient.

Applying the principle of virtual work, the above strong form can be expressed in a weak form as follows: 2.4and 2.5Note that this set of thermo-mechanical equations is mechanically quasi-static and thermally transient but reduces to steady-state heat transfer when the time derivative term becomes zero. The last terms in equations (2.4) and (2.5) are due to the embedded discontinuity, which accounts for the work due to cohesive stresses and the heat flow at an interface or a cohesive crack surface.

### (b) A thermo-mechanical cohesive zone model

The thermo-mechanical cohesive zone model (TM-CZM) used in this paper to describe the cohesive stresses and heat flux across the crack surfaces as functions of crack displacements is composed of two parts: (i) a mixed-mode mechanical CZM with piece-wise linear traction–separation relationships for normal and shear and fracture modes as shown in figure 2*a*,*b* and (ii) a conductance–normal separation law as shown in figure 2*c*.

The mechanical part of the TM-CZM as shown in figure 2*a*,*b* has piece-wise linear traction–separation laws for normal and shear separation. Here *δ*_{nc} and *δ*_{sc} are the critical normal and tangential separations beyond which cohesive stresses become zero under the respective pure modes. The cohesive stresses and separations are defined in the local coordinate system, (*s*,*n*), where *s* and *n* are directions tangential and perpendicular to the crack plane as shown in figure 1. Each linear segment of the traction–separation laws is indexed with free indices *i*,*j*∈[1,4] as labelled. Each segment is defined in terms of characteristic stresses ( or ), characteristic crack displacements ( or ), and constant-valued cohesive stiffness ( or ), such that for the cohesive model
2.6where and . Each of the constant cohesive stiffness segments corresponds to a crack displacement range consistent with the indexed segments as follows:
2.7From these notations, the piece-wise linear cohesive law may be expressed as
2.8where sgn(⋅) is the sign function. We further note that the segment index-based cohesive stress calculation of equation (2.8) and the associated displacement range definition of equation (2.7) permit small but negative opening (i.e. *δ*_{n}<0). In such a case, the slope of segment 1 in the mode I traction–separation law, which is typically set to be as large as possible numerically but not to cause numerical instability, serves a purpose similar to the penalty method to minimize the interpenetration of cohesive crack surfaces under severe compression.

Here and in equation (2.8) are the maximum normal and shear separations ever experienced by the mode I and mode II fracture process, respectively. These maximum crack separations are solution-dependent variables used to account for the irreversibility of the cohesive model associated with the unloading/reloading segment indicated by index 4. They have corresponding normal and shear cohesive stresses and that can be calculated from equation (2.8). Thus, in segment 4, the traction–separation law is given by 2.9where the cohesive slopes and separation ranges for segment 4 are 2.10This indexing procedure applies to any piece-wise linear function and thus other cohesive laws that can be linearized and expressed in a similar form.

Under mixed-mode conditions, it is assumed that the total energy release rate is the sum of the mode I (opening) and mode II (shear) contributions, i.e.
2.11where and are calculated from the traction–separation law as
2.12Complete fracture is to occur when the following critical condition is satisfied:
2.13where *Γ*_{IC} and *Γ*_{IIC} are the mode I and mode II fracture toughness in the linear elastic fracture mechanics (LEFM) context, which are equal to the total areas under the respective traction–separation curves. Here *δ*_{nf} and *δ*_{sf} are the normal and shear crack displacements when the failure condition equation (2.13) is met. This mixed-mode mechanical cohesive law is further explained in [32–34]. It guarantees correct mode mixedness when LEFM conditions are satisfied.

For the thermal response of the TM-CZM, following [35–37], it is assumed that the heat flux *q*_{c} across the cohesive crack can be computed as the product of a cohesive conductance coefficient *h*_{c} and the temperature jump across the cohesive crack Δ*θ*, i.e.
2.14where *h*_{c} can be calibrated from micromechanics consideration or proper experimental tests of fracture damage, and it is dependent on the physical process. Here, a simple, piece-wise linear conductance–separation law shown in figure 2*c* that captures key aspects of heat transfer through an interface based on physical considerations is adopted.

## 3. Two- and three-dimensional thermal–mechanical augmentedfinite-element method formulation

In this section, we will describe how to augment a standard FE to account for intra-elemental cohesive cracks under general thermo-mechanical loading. From equations (2.4) and (2.5), the discretized weak form of the thermo-mechanical equations is
3.1*a*and
3.1*b*where **u** and ** θ** are the nodal displacement (mechanical) and temperature (thermal) d.f., respectively.

**N**is the shape function matrix and

**B**is the shape function derivative matrix.

**D**and

**k**are the material stiffness and conduction matrix, respectively. For the temperature derivative, a backward finite difference scheme is used to evaluate the time integration for the transient heat transfer process, i.e. with Δ

*t*

_{n}=

*t*

_{n}−

*t*

_{n−1}. Here, the subscript

*n*denotes the current time step.

Note that the last term in the right-hand side of equation (3.1a) represents the thermal stress induced strain energy, and it is a result of the thermal–mechanical coupling. In general, the temperature at the current time (*θ*_{n}) has to be solved together with the mechanical d.f. (**u**_{n}), i.e. solving equations (3.1a,*b*) simultaneously. However, as will be shown shortly, using a constant cohesive conductance *h*_{c} will allow us to solve the heat transfer equation of equation (3.1b) first, and then to solve the mechanical equation of equation (3.1a). This is the reason we place this term in the right-hand side of equation (3.1a). This enables an efficient procedure to solve the mechanical and heat transfer equations separately and sequentially.

### (a) Augmentation scheme for a two-dimensional 4-node quadrilateral element

Now consider a 4-node quadrilateral element as shown in figure 3. There are two possible cut configurations if the element is traversed by a cohesive crack. In figure 3*a*, the element is cut into two quadrilateral subdomains, and in figure 3*b* it is cut into a triangular and a pentagonal subdomain. In either case, it is augmented with internal nodes 5, 6, 5^{′} and 6^{′} to facilitate subdomain stiffness and crack displacement calculation. The augmentation procedure assumes that the crack tip always resides at an element boundary during growth and that the crack tip is non-singular due to the existence of the cohesive tractions. The intra-element crack has a length *l*_{e} and is oriented at an angle *ϕ* to the global Cartesian coordinate system. Thus, the displacements and nodal forces in the local coordinate system (*s*, *n*) are related to those in the global coordinates (*x*, *y*) by a rotation matrix as follows:
3.2As a convention the internal nodes will always be assigned as 6^{′} and 5^{′} for subdomain *Ω*^{+} and, 6 and 5 for subdomain *Ω*^{−}. The associated internal d.f., loads and heat flows at the current time step *n* can thus be grouped as , etc. The arrays for the external nodal d.f. are dependent on whether the subdomains are quadrilateral–quadrilateral (figure 3*a*) or triangular–pentagonal (figure 3*b*).

For quadrilateral–quadrilateral cut, ; ; ; , etc. For triangular–pentagonal cut, ; ; ; , etc. Other loads and heat flux arrays are also defined similar to the displacement d.f.

Following the discretization scheme detailed in [29,30], the following thermal–mechanical equilibrium equations for the two subdomains can be derived:
3.3*a*
3.3*b*
3.3*c*
3.3*d*where , and

is the stiffness matrix for subdomain + or − (IP1 is the number of integration points and

*w*_{k}and*J*_{k}are the weight and Jacobian for integration point*k*);is the conduction matrix for subdomain + or −;

is the specific heat matrix for subdomain + or −;

is the equivalent nodal force array

*at external nodes*integrated from surface traction (IP2 is the number of integration points used for surface traction integration and*w*_{j}and*J*_{j}are the corresponding weight and Jacobian for integration point*j*);is the equivalent nodal force array

*at internal nodes*integrated from cohesive stresses (**t**_{c}) between the cohesive crack surfaces;is the equivalent nodal heat flow array

*at external nodes*integrated from*q*;is the equivalent nodal heat flow array

*at internal nodes*integrated from the heat flux across the cohesive crack;and are the equivalent nodal force array in subdomain + or − due to thermally induced stresses.

Note that in a traditional displacement-based FEM, the external d.f., and , are typically given as trial quantities for element equilibrium approximation. However, the internal d.f., {**u**_{65}}_{n}, {**u**_{6′5′}}_{n}, {*θ*_{65}}_{n} and {*θ*_{6′5′}}_{n}, are not known because they are not part of the global solution. These internal d.f. have to be solved as part of the elemental equilibrium consideration (i.e. element condensation), which will be detailed in the next section.

### (b) Element condensation

Across the crack planes, the stresses and heat flux remain continuous, while the displacements and temperature are not. We first solve the heat transfer equation by assuming the heat conductance coefficient across the cohesive crack to be a constant, i.e. . Using the second-order Gaussian integration scheme (which can be conveniently done through an integration matrix **T** and interpolation matrix **N**—see, for example, [38] for more details) the heat flow at the internal nodes can be obtained as
3.4*a*where
3.4*b*Here Diag[⋅] defines a square, diagonal matrix where the diagonal elements are composed of the elements enclosed in the brackets. Here {*θ*_{6′5′}}_{n} and {*θ*_{65}}_{n} can be readily solved by substituting equation (3.4a) into equations (3.3b,*d*), which results in
3.5*a*and
3.5*b*where
3.5*c*With the temperature fields known in both subdomains, the mechanical equation of equations (3.3a,*c*) can be solved following exactly the same condensation procedure detailed in [29,30]. The only difference is due to the thermal stress terms, i.e. and in equations (*a*,*c*). However, these terms are now known because the temperature field has already been solved from equation (3.5). Based on the mechanical CZM described in §2*b* and the integration procedure detailed in [29,30], the equivalent nodal forces at the internal nodes can be obtained as
3.6*a*where
3.6*b*Here, the indices *i*, *j*, *k*, *l* (each can range from 1 to 4) designate the segment of the cohesive traction–separation laws to which the characteristic coefficients and slopes belong (four such free indices needed for two fracture modes at two integration points). Following [29,30], the internal displacements can be derived as explicit functions of the external displacements and temperatures as follows:
3.7where
Equation (3.7) remains highly nonlinear because the matrices **S**_{0}, **A**_{0}, **A**, **B**, and are all nonlinear functions of crack separation Δ**u** (=**u**_{6′5′}−**u**_{65}) through *α*_{0} and *σ*_{0} in equation (3.6). However, for the piece-wise linear cohesive laws in figure 2, a consistency-check-based solution algorithm has been established in [29,30]. It is directly applicable to solve equation (3.7) and we shall not repeat the solving procedure here.

Once the true solution to equation (3.7) is found the nonlinear matrices **S**_{0}, **A**_{0}, **A**, **B**, and are all established, the thermal–mechanical equilibrium equation can be readily derived by substituting equations (3.5) and (3.7) into equation (3.3)
3.8*a*
3.8*b*

### (c) An augmented three-dimensional 4-node tetrahedron element

Formulation-wise extension of the 2D thermal–mechanical augmented finite-element method (TM-AFEM) to 3D is relatively straightforward. Figure 4 shows the augmentation process of a 4-node tetrahedron element (figure 4*a*). There are two possible configurations if the element is cut by a crack: (i) the element is cut into one tetrahedron subdomain and one wedge subdomain as shown in figure 4*b* and (ii) the element is cut into two wedge subdomains as shown in figure 4*c*.

The tetrahedron–wedge configuration results in a triangular crack plane (figure 4*b*) and six internal nodes 5, 6, 7, 5′, 6′ and 7′ are assigned to the top and bottom surfaces. The wedge–wedge configuration has a quadrilateral crack plane and eight internal nodes 5, 6, 7, 8, 5′, 6′, 7′ and 8′ are assigned as shown in figure 4*c*. Thus, the displacement jumps across the cohesive crack are simply the difference in displacements between node-pairs 5–5′, 6–6′, 7–7′ and 8–8′. Defining a local coordinate system as shown in figure 4*d* with the surface out-normal (** N**) as the local

*z*-axis, the rotation matrix can be obtained as 3.9

*a*where

*l*

_{i},

*m*

_{i}and

*n*

_{i}(

*i*=1,2,3) are the direction cosines of the crack plane with

**=**

*N**l*

_{3}

**+**

*i**m*

_{3}

**+**

*j**n*

_{3}

**being the out-normal direction,**

*k***=**

*T**l*

_{2}

**+**

*i**m*

_{2}

**+**

*j**n*

_{2}

**and**

*k***=**

*S**l*

_{1}

**+**

*i**m*

_{1}

**+**

*j**n*

_{1}

**being the two orthogonal in-plane shear directions. 3.9**

*k**b*and

*α*

_{i},

*β*

_{i}and

*γ*

_{i}are the angles between the local and global axes as illustrated in figure 4

*d*.

The rest of formulae are identical to those of the 2D TM-AFEM (equations (3.3)–((3.8))) except that dimensions of the respective matrices are different due to different numbers of d.f. at internal and external nodes. The consistency-check-based condensation algorithm is also directly applicable to the 3D TM-AFEM.

One of the biggest challenges in 3D fracture modelling is to maintain the conformity of an evolving crack surface. A crack surface tracking scheme is necessary and proper numerical treatments are needed to guarantee that the crack surface, which may cut multiple elements before they merge at some point, maintains at least C^{0}-continuity. Several crack tracking schemes have been developed in the literature [9,39–41]. Here, we adopt the local tracking algorithm proposed by Areias & Belytschko [9]. The detailed algorithm will not be given here but interested readers are referred to [9]. The necessary numerical treatment to guarantee C^{0}-continuity of a crack surface can be found in [42].

## 4. Numerical examples

The TM-AFEM described above can be integrated into any general purposed FEM programmes. In this study, it is integrated into a commercial software package ABAQUS (v. 6.11) as user-defined elements. A more in-depth discussion of the implementation can be found in [43]. In the remainder of this section, we shall use numerical examples to demonstrate the simulation capabilities of the proposed 2D and 3D TM-AFEM.

### (a) Thermal and mechanical crack propagation in a cruciform plate

In this section, the crack propagation in a cruciform plate with an initial corner crack as shown in figure 5*a* is investigated under three different loading conditions (LCs): a pure thermal loading, a pure mechanical loading and a combined thermo-mechanical loading. The plate has an initial traction-free corner crack of length *a*=0.2*L* (*L*=100 mm) and oriented 135° with respect to the horizontal direction (figure 5*a*). The boundary element method (BEM) solutions of crack paths based on LEFM were first reported in [44], and it will be used here for validating the TM-AFEM predictions. The LEFM solution indicates that the crack path is heavily dependent on the thermal–mechanical boundary conditions imposed on the four plate edges (AB, CD, EF and GH) illustrated in figure 5*a*.

The cross-plate is discretized into 3380 quadrilateral plane stress elements. The elastic and thermal material properties of the cruciform plate as listed in table 1 are directly from [44]. The cohesive strength and fracture toughness are chosen so that the estimated cohesive zone size is , about five times the numerical mesh size of *l*_{e}=4 mm. We note that in the A-FEM simulations below, while the initial cracking at the corner is due to the pre-planted existing crack, further propagation of the crack is determined entirely by the A-FEM using the maximum principal stress criterion.

Three different sets of LCs are simulated using the newly developed TM-AFEM. For the pure thermal loading case (LC #1): the top (AB) and bottom (EF) edge are prescribed with a uniform temperature change of 10°C and −10°C, respectively, whereas the left (GH) and right (CD) edges are roller-supported in the horizontal direction; LC #2 is a pure mechanical loading with top (AB) subjected to a 10 N mm^{−2} normal traction and all other edges roller-supported; LC #3 is a fully coupled thermal mechanical loading: the top (AB) and bottom (EF) are prescribed with a temperature change of 10°C and −10°C, respectively, whereas the left (GH) and right (CD) edges are assigned a constant temperature change of −5°C. In addition, the top edge (AB) is subjected to a normal traction of 10 N mm^{−2}.

The LEFM/BEM solutions of [44] for the three LCs are reproduced in figure 5*b* by the dashed lines. Under the thermal loading case (LC #1), the crack initially deviates from 135° to about 90° and propagates upwards due to the horizontal tension caused by the thermal contraction of the bottom edge (EF) and the horizontal confinement at edges CD and GH. However, the propagation is increasingly difficult because the upper half of the cruciform is under horizontal compression, caused by the thermal expansion from the higher temperature edge AB and the horizontal confinement of CD and GH. As a result, the crack further deflects to roughly the horizontal direction.

Under the pure mechanical loading of LC #2, it is expected that the crack will orient itself to a direction that is perpendicular to the loading direction, which in this case is in line with the maximum principal stress direction. Under the combined thermal and mechanical loading, the crack path is in between the pure thermal and the pure mechanical crack path.

The TM-AFEM-simulated crack paths are all in excellent agreement with the LEFM/BEM solutions. Moreover, with the help of our inertia-based stabilizing algorithm, the 2D TM-AFEM is able to predict complete crack paths until the specimen is entirely severed.

The same problem was also modelled by the 3D TM-AFEM with 16 900 tetrahedron elements and similar results were obtained (figure 5*b*).

### (b) Prediction of crack patterns in a thermal barrier coating system

In this section, we use the newly developed TM-AFEM to analyse the multiple crack development in a thermal barrier coating (TBC) system reported in [45]. Figure 6*a* illustrates the three layers of the TBC: a top surface ceramic coating (yttria-stabilized zirconia, or YSZ) of thickness *l*_{c} ranging from 0.25 to 1.5 mm, a bond-coating alloy of thickness approximately 0.16 mm joining the top ceramic coating and the steel substrate of thickness 1.5 mm.

During the experiments, the top surface was raised to 1000–1600°C after being subjected to intensive laser heating. During the cooling period, a complex multiple crack system was observed in the top ceramic coating layer (but not in the bond coating layer nor in the substrate): the significant difference in material properties of metal substrate and the ceramic coating results in the initiation of surface cracks from the top surface propagating towards the interface, gradually deflecting along the interface, and sometimes the interface cracks joining neighbouring surface cracks to form local spalling of the protective coating, leading to eventual failure of the material system. Some of the typical cracking patterns observed experimentally are illustrated in figure 6*b*. Crack patterns that result have been found to depend on the coating thickness, elastic moduli mismatch among the layers and the fracture properties. Here, the surface cracks refer to vertical cracks that appear on the ceramic coating surface propagating towards the interfaces.

The TM-AFEM model and geometric dimensions for the three-layer TBC system is shown in figure 6*a*. The dimensions approximate the section of the larger TBC test specimen subjected to approximately uniform high heat flux laser heating in [46]. Two numerical models are made for ceramic coating thicknesses of *l*_{c}=0.25 and 0.75 mm with an element size of 0.05×0.1 mm for meshes of 2400 and 3000 TM-AFEM elements, respectively. The relatively small mesh size of 0.05 mm in the thickness direction is used to guarantee that a surface crack can propagate and deviate in the surface coating layer before joining with interface delamination cracks. The left and right edges of the model are constrained in the *x*-direction while the bottom surface is constrained in the *y*-direction, mimicking the confined local heating in the experiments.

The thermal–mechanical material properties listed in table 2 are summarized from [46,47]. Note that the interface toughness and strength were estimated in [47] through a detailed LEFM analysis of a numerical model, assuming idealized crack configurations based on experimental observation and zero-thickness interfaces. It was reported that interface toughness and strength values are approximately twice the respective values of the surface ceramic coating (YSZ). The fracture properties indicate that the cohesive zone size is *l*_{coh}∼0.25 mm, which can be well resolved by the numerical models described earlier.

Since according to the experimental observation, there is no crack development within the bond-coating alloy and the substrate, only the top ceramic YSZ coating was modelled with the TM-AFEM elements. The bottom layer of TM-AFEM elements, which is in junction with the bond-coating alloy, is designated to be the interface elements, which have twice the fracture toughness and cohesive strength when compared with those of the ceramic coating.

The top surface of the TBC system is prescribed with a temperature drop of Δ*θ*_{top_surface}=−1000 K, whereas the bottom substrate surface has the stress-free reference temperature of *θ*_{bottom_surface}=0 C.

The TM-AFEM predicted steady-state temperature field is shown in figure 7. Because the top YSZ coating has a very small thermal conductivity while the conductivity coefficients of the bond-coating alloy and the metal substrate are orders of magnitude larger, the temperature gradient is mostly in the top YSZ coating, causing significant thermal stresses mainly in this surface coating. The temperature in the bond-coating alloy and the substrate is not much different from the reference temperature. This is in good agreement with the measured temperature profiles reported in [47].

The TM-AFEM predicted a progressive crack formation process as temperature drop in the TBCs with coating thicknesses of 0.25 mm is shown in figure 7. For the thinner coating of *l*_{c}=0.25 mm, the A-FEM predicts that a first surface crack initiates in the middle of the top surface at a temperature drop of Δ*θ*_{top_surface}=−250°C. As the temperature further drops to Δ*θ*_{top_surface}=−575°C, more surface cracks are initiated and they all propagate towards the interface. However, the surface crack propagation becomes increasingly difficult because (i) the bond-coating alloy and the substrate is much stiffer than the YSZ coating and (ii) the interface toughness is much larger than the fracture toughness of the YSZ coating. (The reason that a hard substrate combined with a tougher interface deters crack propagation towards the interface has been well understood, owing to the work of Hutchinson and co-workers, e.g. [48,49].) Thus, some of the surface cracks start to deviate from the vertical direction towards the horizontal direction.

As the temperature drops further to Δ*θ*_{top_surface}=−736°C, one of the left corner surface cracks has transitioned into a delamination crack, and the delamination crack propagates along the interface quickly and links to a neighbouring surface crack, forming a spall-off region. Thus, this region is isolated thermally from the rest of the domain and a large discontinuity in temperature exists between the spall-off region and the rest of the domain. Further dropping the surface temperature to Δ*θ*_{top_surface}=−785°C, more delamination cracks form as shown in the bottom micrograph of figure 7.

The predicted final saturated crack spacing is 0.75 mm, which is also in good agreement with the experimental report of 0.70±0.24 mm [46]. Thus, the TM-AFEM predicted crack formation process is qualitatively in agreement with the experimental observation of Choules *et al.* [46], who reported that the horizontal interface cracks form at larger temperature drops and specifically underneath the surface cracks.

For the thicker coating thickness of *l*_{c}=0.75 mm, it is predicted that the first surface crack initiation occurs at Δ*θ*_{top_surface}=−292°C. As the surface temperature drops further to −495°C, more surface cracks develop. The predicted final saturated crack spacing is 1.2 mm, consistent with the experimental value of 1.42±0.37 mm [46]. The crack spacing is larger than that in the thinner coating of *l*_{c}=0.25 mm, because it is less effective to build up the horizontal tensile stress in a thicker layer between two surface cracks to form a new surface crack in between.

For both cases, we also examined the thickness effects of the bond-coating alloy by varying its thickness from 0.1 to 0.3 mm and assuming a linear distribution of elastic and thermal properties from the YSZ–bond coat alloy interface to the bond coat alloy–substrate interface. No significant changes in cracking patterns were observed numerically.

### (c) Complex damage evolution in a three-dimensional textile composite

In this example, the 3D TM-AFEM is used to simulate the multiple crack initiation in a C/SiC angle interlock weave shown in figure 8*a*,*b*, which is a 3D reconstruction from a high-resolution μCT measurement reported in [2,50,51]. The calibrated tow elasticity is numerically similar to the elasticity deduced for a similar C/SiC material, reported in [52]. Nonlinear cohesive laws were calibrated from the same tests by matching the predicted matrix cracking stress with the experiment. In this preliminary demonstration, we shall focus on the early stage of matrix cracking only, which occurred at a global strain level of approximately 0.1% according to the experimental observation in [2,52]. This is due to a severe difficulty in explicit meshing of the two/matrix interfaces. Therefore, the interfaces are assumed to be rigidly bonded (i.e. tow elements and matrix elements sharing the same nodes).

Key features of the crack pattern, including the approximate periodicity of the micro-cracks, which form above near-surface weft tows, and the alignment of micro-cracks in diagonal bands (one of which is outlined by the ellipse in figure 8*c*), are well replicated by the 3D TM-AFEM. The locations of the cracks relative to underlying weft tows are also well predicted by the 3D A-FEM simulation. The micro-cracks seen in experiments form above weft tows and propagate both into the composite and across the surface of the composite, tracking the weft tow. However, the propagation across the surface is limited: when the micro-crack approaches the region where the weft tow dips beneath a warp tow, propagation is arrested. Thus, the micro-cracks remain of finite length (figure 8*c*). This behaviour is replicated well by the 3D A-FEM simulations.

The phenomenon of crack arrest is correlated with lowering of stress in the superficial matrix by straightening of the warp tow. (Earlier evidence for tow straightening in this class of composites was reported in [53].) Both these stress heterogeneity and the critical stress for micro-cracking depend on geometrical details such as the thickness of the superficial matrix layer and the orientations of tows. They will therefore be affected by stochastic variance in the tow architecture.

## 5. Conclusion

In this paper, a TM-AFEM has been proposed, implemented and validated for steady-state or transient, coupled thermal–mechanical analyses of complex materials with explicit consideration of arbitrary evolving cracks. A detailed augmenting process for a 2D 4-node quadrilateral element and a 3D 4-node tetrahedron element has been derived. It has been shown that one of the distinct features of the current TM-AFEM is that it can account for arbitrary intra-element cracks and their interaction without the need for additional external d.f. as in X-FEM or extra nodes as in PNM–FEM. The new TM-AFEM introduces internal node-pairs with regular displacements as internal nodal d.f., which are eventually condensed at elemental level. Thus, the crack displacements and temperature jump across the crack plane are natural outcomes of the elemental equilibrium consideration. The TM-AFEM is further empowered by a novel elemental condensation algorithm that can provide solution to the element equilibrium without the numerically expensive iterative solving process.

Through the numerical examples in §4, it has been demonstrated that the new TM-AFEM achieved very significant improvements in numerical accuracy, efficiency and stability. The capability of the TM-AFEM in dealing with coupled multiple crack development under thermal loads has also been demonstrated with a multi-layered thermal barrier coating (TBC) system. All of the experimentally observed crack types, including surface cracks, interface delamination cracks and their interactions were reproduced and in good agreement with experiments. Finally, the potential of the A-FEM’s capabilities in high-fidelity simulations of interactive cohesive cracks in 3D complex heterogeneous materials has also been demonstrated through the multiple, 3D crack modelling in a 3D textile CMC specimen.

## Competing interests

We declare we have no competing interests.

## Funding

The authors are grateful to the support from the US Army Research Office (ARO grant no. W911NF-13-1-0211). B.C.D. gratefully acknowledges a doctoral fellowship support from the Florida Space Grant Consortium (FSGC).

## Footnotes

One contribution of 22 to a Theo Murphy meeting issue ‘Multiscale modelling of the structural integrity of composite materials’.

- Accepted March 7, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.