rsta.royalsocietypublishing.org

Research Cite this article: Jung J, Do BC, Yang QD. 2016 Augmented finite-element method for arbitrary cracking and crack interaction in solids under thermo-mechanical loadings. Phil. Trans. R. Soc. A 374: 20150282. http://dx.doi.org/10.1098/rsta.2015.0282 Accepted: 7 March 2016 One contribution of 22 to a Theo Murphy meeting issue ‘Multiscale modelling of the structural integrity of composite materials’. Subject Areas: materials science, computer modelling and simulation Keywords: augmented finite-element method, fracture, cohesive zone models, composites Author for correspondence: Q. D. Yang e-mail: [email protected]

Augmented finite-element method for arbitrary cracking and crack interaction in solids under thermo-mechanical loadings J. Jung, B. C. Do and Q. D. Yang Department of Mechanical and Aerospace Engineering, University of Miami, Coral Gables, FL 33124, USA QDY, 0000-0003-0574-6461 In this paper, a thermal–mechanical augmented finiteelement 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

2016 The Author(s) Published by the Royal Society. All rights reserved.

2 .........................................................

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

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

3

mechanical

thermal

∇ · σ +/− = 0

in Ω +/−

+/− σ +/− · n+/− = ¯t

on Γt

u+/− = u¯ +/−

on Γu

+/−

tc

ρ +/− c+/−

in Ω +/−

+/−

q+/− · n+/− = q¯ +/−

on Γq

+/−

θ +/− = θ¯ +/−

on Γθ

+/−

= σ +/− · n+/−

∂θ +/− + ∇ · q+/− = 0 ∂t

+/−

on Γc

qc

+/−

(2.1)

+/− +/−

= q+/− · n+/−

on Γc

where σ +/− 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, tc = (t− c = −tc ) and qc (= qc = −qc ), are functions of crack displacement and temperature jump across the discontinuity, i.e. tc = tc (u) = tc (u+ − u− ) qc = qc (θc ) = qc (θ + − θ − )

on Γc .

(2.2)

Assuming 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 +/−

σ +/− = C+/− · (εT

+/−

− εθ

);

+/−

εT

+/−

= [∇u+/− + (∇u+/− )T ]/2;

εθ

= α +/− θ +/−

q+/− = −k+/− · ∇θ +/− ,

and (2.3) +/−

where 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:    Ω +/−

and

 Ω +/−

(ρc)

σ : ∇δu dΩ =

∂θ δθ dΩ + ∂t

+/−

Γt

¯t · δu dΓ +



+/−

Γc

(σ · nc ) · δu dΓ

 Ω +/−

q · ∇δθ dΩ =

(2.4)

 +/−

Γq

δθ(q · n) dΓ +

+/−

Γc

δθc (q · nc ) dΓ .

(2.5)

Note 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.

.........................................................

Suppose a domain Ω as shown in figure 1 is divided by a discontinuity, Γc = Γc+ ∪ Γc− . 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 u¯ +/− and surface +/− +/− +/− on the boundaries Γu and Γt , respectively. Likewise, thermally, temperature traction ¯t +/− +/− θ¯ +/− and heat flux q¯ +/− are prescribed on boundaries Γθ and Γq . For such a problem, the strong forms of the mechanical and thermal equilibrium equations under both mechanical and thermal boundary conditions are

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

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

n

4

s

G +t G +q

–– t

–q –

G –t

G –q

W–

G –q

+

q–c

G +u

G +q

G –u

– q–

G –c W

Gc

t+c

n+

G +c

n–

q+c

–u+

– q+

–u–

Figure 1. A thermo-elastic body with a discontinuity under general thermo-mechanical boundary conditions. (Online version in colour.) (a)

(b)

sˆ1 sˆ* sˆ2

tˆ1 tˆ* tˆ2

2 1

–dsc –ds2

4

dn1

dn*

dn2

dnc

1

hˆ0

2 3

GII ds1 ds* ds2 –tˆ2 –tˆ* –tˆ1

4

dn

hc

–ds1

3

GI

(c)

t

s

dsc

ds

1

hˆ1 hˆ

2

2

3 4

dqn1

dqn2

dqnc dn

Figure 2. Mode I (a) and Mode II (b) traction–separation laws for cohesive stresses between crack surfaces and (c) the conductance–normal separation law for heat transfer across a cohesive crack. The numbering labels index the respective linear segments in the piece-wise linear laws. For 3D CZM, the mode III traction–separation law is assumed to be identical to that of the mode II. (Online version in colour.)

(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 2a,b and (ii) a conductance–normal separation law as shown in figure 2c. The mechanical part of the TM-CZM as shown in figure 2a,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 (τˆ (i) or σˆ (j) ), characteristic (j)

(i)

(i)

(j)

crack displacements (δs or δn ), and constant-valued cohesive stiffness (αs or αn ), such that for the cohesive model σˆ (1) = σˆ 1 , (1)

δn = δn1 , (i) αs

=

σˆ (2) = σˆ 2 , (2)

δn = δn2 ,

τˆ (i) − τˆ (i−1) (i)

σˆ (3) = 0,

(i−1)

δs − δs

and

τˆ (1) = τˆ1 ,

τˆ (2) = τˆ2 ,

(1)

(2)

(3)

δn = δnc , (j) αn

=

δs = δs1 , δs = δs2 ,

σˆ (j) − σˆ (j−1) (j)

(j−1)

δn − δn

,

⎫ ⎪ ⎪ ⎪ ⎪ ⎪ (3) δs = δsc ,⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭

τˆ (3) = 0,

(2.6)

.........................................................

t–c

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

–+ t –q +

(0)

(0)

n

n

n

n

From these notations, the piece-wise linear cohesive law may be expressed as ⎫    ↔(i) (i) (i−1) ⎪ |δs | ∈ δ s && |δs | > δs∗ ,⎪ τ (δs ) = sgn (δs ) τˆ (i−1) + αs |δs | − δs , ⎬ (j)

σ (δn ) = σˆ (j−1) + αn



(j−1)

δn − δn



↔(j)

δn ∈ δ n && δn > δn∗ ,

,

⎪ ⎪ ⎭

(2.8)

where 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 (1) (αn > 105 ∼ 106 ) 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 δn∗ and δs∗ 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 ⎫ ↔(4) (4) ⎬ τ (δs ) = αs δs for |δs | ∈ δ s && |δs | ≤ δs∗ ⎪ (2.9) ⎪ ↔(4) ⎭ (4) σ (δn ) = αn δn for δn ∈ δ n && δn ≤ δn∗ , where the cohesive slopes and separation ranges for segment 4 are (4)

αs = ↔(4) δs

τˆ ∗ , δs∗

(4)

αn =



= 0, δs∗ ,

σˆ ∗ , δn∗

↔(4) δn

 ∗



(2.10)

= −∞, δn .

This 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 G is the sum of the mode I (opening) and mode II (shear) contributions, i.e. G(δn , δs ) = GI (δn ) + GII (δs ), where GI and GII are calculated from the traction–separation law as  δnf  δsf σ (δ) dδ; GII = τ (δ) dδ. GI = 0

(2.11)

(2.12)

0

Complete fracture is to occur when the following critical condition is satisfied: GI (δnf ) ΓIC

+

GII (δsf ) ΓIIC

= 1,

(2.13)

where Γ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.

.........................................................

n

5

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

where σˆ (0) = τˆ (0) = 0 and δn = δs = 0. Each of the constant cohesive stiffness segments corresponds to a crack displacement range consistent with the indexed segments as follows: ⎫   ↔(i) (i−1) (i) ⎪ (i = 1, 2, 3) for mode II⎪ , δs δ s = δs ⎬ (2.7)     ⎪ ↔(1) ⎪ (j−1) (j) (1) ↔(j) ⎭ ; (j = 2, 3) for mode I. = −∞, δ = δ , δ δ δ

(2.14)

where hc 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 2c that captures key aspects of heat transfer through an interface based on physical considerations is adopted.

3. Two- and three-dimensional thermal–mechanical augmented finite-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

 Ω +/−

  BT DB dΩ u =

+/−

Γt

NT ¯t dΓ +



 Γc

NT tc dΓ +

Ω +/−

BT DαN(θ − θ 0 ) dΩ

(3.1a)

and

 Ω +/−



 NT ρcNdΩ θ˙ +

Ω +/−

  BT kB dΩ θ =

 Γc

NT qc dΓ +

+/−

Γq

NT q¯ dΓ ,

(3.1b)

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. θ˙ = θ n − θ n−1 /tn with tn = tn − tn−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. (un ), i.e. solving equations (3.1a,b) simultaneously. However, as will be shown shortly, using a constant cohesive conductance hc 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 3a, the element is cut into two quadrilateral subdomains, and in figure 3b 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 le 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

.........................................................

qc = hc θc ,

6

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

For the thermal response of the TM-CZM, following [35–37], it is assumed that the heat flux qc across the cohesive crack can be computed as the product of a cohesive conductance coefficient hc and the temperature jump across the cohesive crack θ , i.e.

(a)

P4,q4 F4x, u4

y n

s

6¢

x le

6

4

5¢

6¢

W– 1

(b)

4 W+ 6¢

5¢ le

F4y, v4

6¢ Dv

y

6

5

6

W–

1

1

F1x, u1 P1, q1

2

2

F2x, u2 P2, q2

F1y, v1

dn

ds

Du

x W–

5¢

4 + W

3

ds

F2y, v2

F1y, v1

5

5 W–

P4, q4 F4x, u4

s f

n

6 1

F1x, u1 P1, q1

2

dn

Du

Dv

5

f

3 5¢

W+

3

2 F2y, v2

F3y, v3 P3, q3 F3x, u3

F2x, u2 P2, q2

Figure 3. (a) The quadrilateral–quadrilateral cut configuration and (b) the triangular–pentagonal cut configuration for a 4-node quadrilateral element with associated d.f., loads and heat flows. (Online version in colour.)

˜ as follows: (x, y) by a rotation matrix Q  ˜ = Q

cos φ − sin φ

 sin φ . cos φ

(3.2)

As 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 (u65 )n = {u6 , v6 , u5 , v5 }Tn , (u6 5 )n = {u6 , v6 , u5 , v5 }Tn , (θ 65 )n = {θ6 , θ5 }Tn , (θ 6 5 )n = {θ6 , θ5 }Tn , etc. The arrays for the external nodal d.f. are dependent on whether the subdomains are quadrilateral–quadrilateral (figure 3a) or triangular–pentagonal (figure 3b). − T For quadrilateral–quadrilateral cut, (u+ ext )n = (u34 )n = {u3 , v3 , u4 , v4 }n ; (uext )n = (u12 )n = {u1 , v1 , + − T T T u2 , v2 }n ; (θ ext )n = {θ3 , θ4 }n ; (θ ext )n = {θ1 , θ2 }n , etc. For triangular–pentagonal cut, (u+ ext )n = (u4 )n = − T ; (θ + ) = {θ }T ; (θ − ) = {θ , θ , θ }T , etc. Other {u4 , v4 }T } {u ; (u ) = (u ) = , v , u , v , u , v n n n n 123 1 1 2 2 3 3 4 1 2 3 n n n n ext ext ext 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: 

L+ 11 L+ 21

L+ 12 L+ 22



u+ ext u6 5

 n



F+ ext = F6 5

 n



a+ + ext a6 5

 , n

(3.3a)

7 .........................................................

W+

3

P3,q3 F3x, u3

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

4

F3y,v3

F4y,v4



R+ 12

 and

R− 11 R− 21

 where

L− 21

+/−

+/− 

R11

R12

R21

R22

+/−

+/−

 L+/− =

+/−

L11

+/− L21

 =

+/−

+/− 

K11

K12

K21

K22

+/−

+/− 

L12

+/− L22

=

+/−

IP1

 +

1 tn

+/−

C12

C21

C22

T +/− B k k=1 wk Jk Bk D



C+ 22

θ+ ext θ 6 5

 ,

(3.3b)

(3.3c) C− 12



C− 22

θ− ext θ 65

+/−

8

n−1

 (3.3d) n−1

+/− 

C11

+/−

C+ 12

, and

is the stiffness matrix for subdomain + or

− (IP1 is the number of integration points and wk and Jk are the weight and Jacobian for integration point k);  +/− +/−  K11 K12  T +/− B is the conduction matrix for K+/− = = IP1 k k=1 wk Jk Bk k +/− +/− K21 K22 subdomain + or −;   +/− +/− C C12  11 T +/− N = IP1 C+/− = k is the specific heat matrix for k=1 wk Jk Nk (ρc) +/− +/− C21 C22 subdomain + or −;  +/− T ¯+/− ) is the equivalent nodal force array at external nodes (Fext )n = IP2 n j=1 wj Jj Nj (tj integrated from surface traction t¯ (IP2 is the number of integration points used for surface traction integration and wj and Jj are the corresponding weight and Jacobian for integration point j);  T (F65 )n = (−F6 5 )n = IP2 j=1 wj Jj Nj (tcj )n is the equivalent nodal force array at internal nodes integrated from cohesive stresses (tc ) between the cohesive crack surfaces;  +/− T q ) is the equivalent nodal heat flow array at external nodes (pext )n = IP2 j n j=1 wj Jj Nj (¯ integrated from q;  T (p65 )n = (−p6 5 )n = IP2 j=1 wj Jj Nj (qcj )n is the equivalent nodal heat flow array at internal nodes integrated from the heat flux across the cohesive    crack;     + −   a aext = IP2 wj Jj BTj D+ α + Nj θ + − θ 0 and ext = IP2 wj Jj BTj D− α − Nj θ − − θ0 j=1 j=1 j j a6 5

a65 n n n n are the equivalent nodal force array in subdomain + or − due to thermally induced stresses.     +/− +/− and θ ext , are Note that in a traditional displacement-based FEM, the external d.f., uext n

n

typically given as trial quantities for element equilibrium approximation. However, the internal d.f., {u65 }n , {u6 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 ˆ . Using the secondconductance coefficient across the cohesive crack to be a constant, i.e. qc = hθ 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

.........................................................

R+ 21  − L11



rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

    + 1 C11 p+ θ+ ext ext = + θ 6 5

p6 5

tn C+ R+ 22 21 n n        − − − L− F a uext 12 = ext + ext u65 F65 a65 L− 22 n n n       − − − − R12 1 C11 p θ ext = ext + − − θ p t R22 65 n C21 65 n n

R+ 11

the internal nodes can be obtained as

9

ˆ h]; ˆ θ n = {θ 6 5 − θ 65 }n . H0 = Th0 N; h0 = Diag[h;

(3.4b)

where

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 − − + −1 + + −1 (θ 65 )n = −(κ B )−1 n R21 (θ ext )n − (κ B )n H0 (κ 22 )n R21 (θ ext )n   − − + −1 + + −1 (κ B )−1 1 n C21 (θ ext )n−1 + (κ B ) H0 (κ 22 )n C21 (θ ext )n−1 + − −1 + −1 + tn +(κ B )−1 n C22 (θ 65 )n−1 + (κ B )n H0 (κ 22 )n C22 (θ 6 5 )n−1

(3.5a)

and −1 − − −1 + + (θ 6 5 )n = −(κ A )n −1 H0 (κ − 22 )n R21 (θ ext )n − (κ A )n R21 (θ ext )n   − −1 − − −1 + + (κ A )−1 1 n H0 (κ 22 )n C21 (θ ext )n−1 + (κ A )n C21 (θ ext )n−1 + , − −1 − −1 + tn +(κ A )−1 n H0 (κ )n C (θ 65 )n−1 + (κ A )n C (θ 6 5 )n−1 22

22

(3.5b)

22

where − C+ C− 22 22 + H0 ; κ 22 n = K− + + H0 ; 22 tn tn − −1 + −1   (κ A )n = κ + (κ B )n = κ − 22 n − H0 κ 22 n H0 ; 22 n − H0 κ 22 n H0 .



κ+ 22



= K+ 22 + n

(3.5c)

With the temperature fields known in both subdomains, the mechanical equation of equations (3.3a,c) can be solved following exactly the same condensation detailed  procedure    a+ a− ext ext in [29,30]. The only difference is due to the thermal stress terms, i.e. and a6 5

a65 n n in equations (3.3a,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 §2b and the integration procedure detailed in [29,30], the equivalent nodal forces at the internal nodes can be obtained as (F65 )n = −(F6 5 )n = S0 + A0 (u6 5 − u65 )n

(3.6a)

where T

T

T

T

˜ ;Q ˜ ]Tσ 0 ; A0 = Diag[Q ˜ ;Q ˜ ]Tα 0 NDiag[Q; ˜ Q] ˜ S0 = Diag[Q   (i,j,k,l) (j−1) (j−1) (i) (i−1) (j−1) (k) (k−1) (l−1) (l) (l−1) T = τˆ (i−1) − αs δs , σˆ − αn δn , τˆ (k−1) − αs δs , σˆ − αn δn σ0 (i,j,k,l)

α0

(i)

(j)

(k)

(3.6b)

(l)

= Diag[αs ; αn ; αs ; αn ].

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

.........................................................

(3.4a)

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

(p65 )n = −(p6 5 )n = H0 θ n ,

as follows:

⎡ −1 − −A−1 A0 Ψ − L21 22 ⎣ + − −1 −B L21 0 22 21 ⎡ ⎤  − −1 −1 −1 A A0 (Ψ 22 ) A ⎦ a65 +⎣ , −1 a6 5

B−1 B−1 A0 (Ψ + ) n 22

n

(3.7)

where + Ψ+ 22 = L22 + A0 ;

− Ψ− 22 = L22 + A0 ,

− −1 A=Ψ+ 22 − A0 (Ψ 22 ) A0 ;

+ −1 B=Ψ− 22 − A0 (Ψ 22 ) A0 .

− Equation (3.7) remains highly nonlinear because the matrices S0 , A0 , A, B, Ψ + 22 and Ψ 22 are all nonlinear functions of crack separation u (= u6 5 − u65 ) 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 S0 , A0 , A, B, Ψ + 22 and − Ψ 22 are all established, the thermal–mechanical equilibrium equation can be readily derived by substituting equations (3.5) and (3.7) into equation (3.3) ⎤  ⎡  − + −1 + −1 − −1 R− −R− θ− 11 − R12 (κ B ) R21 12 (κ B ) H0 (κ 22 ) R21 ext ⎦ ⎣ − −1 − + −1 −1 ¯ + θ+ −R+ R+ ext n 12 (κ A ) H0 (κ 22 ) R21 11 − R12 (κ A ) R21 n ⎡ ⎤   −   − − + −1 + −1 − −1 −R− pext θ− 1 ⎣ C11 − R12 (κ B ) C21 12 (κ B ) H0 (κ 22 ) C21 ext ⎦ = + + tn −R+ (κ A )−1 H0 (κ − )−1 C− pext θ+ C+ − R+ (κ A )−1 C+ ext n

⎡

+ ⎡ ⎣

12

22

− − −1 − 1 ⎣ C12 − R12 (κ B ) C22 tn −R+ (κ A )−1 H0 (κ − )−1 C− 12 22 22

− −1 − L− 11 − L12 B L21

21

11

12

21

n

⎤   + −1 + −1 −R− θ 65 12 (κ B ) H0 (κ 22 ) C22 ⎦ θ 6 5

C+ − R+ (κ )−1 C+ 12

+ −1 + −1 −L− 12 B A0 (Ψ 22 ) L21

12

⎤ ⎦

A

u− ext

22



n

n−1

(3.8a)

n−1

+ −1 + u+ L+ ext n 11 − L12 A L21 ⎡ ⎤   −   −  −1 I − A (Ψ + )−1 L− 0 fext aext 0 12 B 22 ⎦ S0 = + + + −⎣ − −1  −1 S0 fext n aext n I − A A (Ψ ) 0 −L+ 0 12 22 ⎡ ⎤  + −1 −1 −1 L− L− a65 12 B 12 B A0 (Ψ 22 ) ⎣ ⎦ (3.8b) − − −1 −1 −1 a6 5

L+ L+ n 12 A A0 (Ψ 22 ) 12 A

− −1 − −1 −L+ 12 A A0 (Ψ 22 ) L21

(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 4a). 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 4b and (ii) the element is cut into two wedge subdomains as shown in figure 4c.

.........................................................

n

10

⎤  0 S0  −1 ⎦ S0 B−1 I − A0 Ψ + 22 ⎤  − −A−1 L+ 21 ⎦ uext u+ ext −B−1 A (Ψ + )−1 L+

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

⎡    −1 −A−1 I − A0 Ψ − u6 5

22 =⎣ u65 0

(a)

(b)

(c)

FZ

4

FY

s

1

1

2

s

FX

nZ t

FY

FX 2 Y

3 ts FX

7 FZ FX

FZ FY

FZ

6¢

ts

6

FZ

5¢

FZ

3 F Y

nz

FY

sn 8¢ 5 tt tt 7¢ 6 sn ts FY FX

n

b

8 1

g

t

nx

FZ 2

a

7 FY

ny

X

Y s

X

Figure 4. Element illustration for (a) a 4-node tetrahedron element with two possible configurations cut by a crack plane, (b) the A-FE with tetrahedron–wedge cut configuration, (c) the A-FE with wedge–wedge cut configuration and (d) definition of local coordinates for the crack plane. (Online version in colour.)

The tetrahedron–wedge configuration results in a triangular crack plane (figure 4b) 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 4c. 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 4d with the surface out-normal (N) as the local z-axis, the rotation matrix can be obtained as ⎤ ⎡ l1 m1 n1 ⎥ ˜ =⎢ (3.9a) Q ⎣l2 m2 n2 ⎦ , l3 m3 n3 where li , mi and ni (i = 1, 2, 3) are the direction cosines of the crack plane with N = l3 i + m3 j + n3 k being the out-normal direction, T = l2 i + m2 j + n2 k and S = l1 i + m1 j + n1 k being the two orthogonal in-plane shear directions. li = cos αi , mi = cos βi , ni = cos γi

(i = 1, 2, 3),

(3.9b)

and αi , βi and γi are the angles between the local and global axes as illustrated in figure 4d. 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 C0 -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 C0 -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 userdefined elements. A more in-depth discussion of the implementation can be found in [43]. In

.........................................................

Z

tt 6¢ s sn n tt t 5

FX 7¢

5¢

3

11

FY

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

FX

(d)

FZ 4

(a)

(b)

traction A

B

traction A

12

B

H

C

H

C

LC #3, Ref. [43] LC #1, 2D A-FEM LC #2, 2D A-FEM

L a G

LC #3, 2D A-FEM

D

G

D

LC #1, 3D A-FEM LC #2, 3D A-FEM LC #3, 3D A-FEM

F

E

F

E

Figure 5. (a) The geometry and thermal–mechanical boundary conditions of a cruciform plate with an initial corner crack; (b) Comparison of the 2D (dotted lines) and 3D (solid lines) TM-AFEM predicted crack paths against the BEM solution of [44] (dashed lines) for three different loading cases.

Table 1. Material properties of the cruciform plate. material property

value

Young’s modulus, E

218.4 kPa

Poisson ratio, ν

0.3

fracture toughness, KIC

1.5 × 10−3 Pa m−1/2

cohesive strength (σˆ = τˆ )

0.01 N m−2

coefficient of thermal expansion, α

1.67 × 10−5 ◦ C−1

specific heat capacity, c

1050 J kg−1 ◦ C

thermal conductivity, k

0.2 W m−1 ◦ C

heat transfer coefficient, h

1.0 W m−2 ◦ C

.......................................................................................................................................................................................................... .......................................................................................................................................................................................................... .......................................................................................................................................................................................................... .......................................................................................................................................................................................................... .......................................................................................................................................................................................................... .......................................................................................................................................................................................................... .......................................................................................................................................................................................................... ..........................................................................................................................................................................................................

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 5a 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.2L (L = 100 mm) and oriented 135◦ with respect to the horizontal direction (figure 5a). 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 5a. 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 lcoh ∼ EΓIC /σˆ 2 = 20 mm, about five times the numerical mesh size of le = 4 mm. We note that in

.........................................................

LC #2, Ref. [43]

L

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

LC #1, Ref. [43]

(a)

(b)

bond coating alloy

0.16 mm

asc

Hc

asc Hc

bond coat substrate

metallic substrate 12 mm

1.5 mm horizontal cracks (Hc)

Figure 6. (a) Illustration of the three-layer thermal barrier coating (TBC) system and (b) illustration of experimentally observed crack types reported in [45]). (Online version in colour.)

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 5b 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 5b).

(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 6a illustrates the three layers of the TBC: a top surface ceramic coating (yttria-stabilized zirconia, or YSZ) of thickness lc 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

.........................................................

lc

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

asc TBC ceramic coating

13

surface cracks (SC)

Table 2. Thermal–mechanical material properties used in this study. fracture toughness (J m−2 ) 20

cohesive strength (MPa) 20

..........................................................................................................................................................................................................

interface

586

5257

10.67

0.25

8.7

0.89

40

40

53.50

0.30

10.7

25.0

—

—

0.30

15.0

51.9

—

—

..........................................................................................................................................................................................................

bond coat

565

6291

metallic substrate

434

7850

..........................................................................................................................................................................................................

200

..........................................................................................................................................................................................................

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 6b. 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 6a. 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 lc = 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 lcoh ∼ 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 bondcoating 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 bondcoating 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].

.........................................................

specific elastic expansion heat density modulus Poisson’s coefficient conductivity ratio (×10−6 ◦ C−1 ) (W m−1 · ◦ C) (J kg−1 · ◦ C) (kg m−3 ) (GPa) 586 5257 10.67 0.25 8.7 0.89

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

material region ceramic coating

14

first surface crack

NT11

15

NT11 0 –4.792 × 10 –9.583 × 10 –1.438 × 102 –1.917 × 102 –2.396 × 102 –2.875 × 102 –3.354 × 102 –3.833 × 102 –4.313 × 102 –4.792 × 102 –5.271 × 102 –5.750 × 102

Dq = –575°C

NT11 0 –6.131 × 10 –1.226 × 102 –1.839 × 102 –2.452 × 102 –3.066 × 102 –3.679 × 102 –4.292 × 102 –4.905 × 102 –5.518 × 102 –6.131 × 102 –6.744 × 102 –7.357 × 102

Dq = –736°C

NT11 0 –6.545 × 10 –1.309 × 102 –1.964 × 102 –2.618 × 102 –3.273 × 102 –3.927 × 102 –4.582 × 102 –5.236 × 102 –5.891 × 102 –6.545 × 102 –7.200 × 102 –7.854 × 102

delamination connecting two surface cracks – spall off

surface cracks turn into delamination

Dq = –785°C

Figure 7. TM-AFEM predicted cracking patterns as a result of applied temperature drop for the TBCs with lc = 0.25 mm.

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 lc = 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 lc = 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

.........................................................

Dq = –250°C

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

0 –2.083 × 10 –4.167 × 10 –6.250 × 10 –8.333 × 10 –1.042 × 102 –1.250 × 102 –1.458 × 102 –1.667 × 102 –1.875 × 102 –2.083 × 102 –2.292 × 102 –2.500 × 102

(a)

(ii) 0.1% global strain (iii)

0.1%

Figure 8. (a) 3D TM-AFEM model generated from the µCT data; (b) the same mode as (a) with matrix element removed to show the mesh for the reinforcement fibre tows. (c)(i) textile preform coated with thin layer of matrix (C/SiC angle interlock); (ii) observations of micro-cracks by digital image correlation in composite formed by full infiltration of SiC matrix into preform in (i); (iii) 3D TM-AFEM predicted crack pattern.

thinner coating of lc = 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 8a,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 8c), 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

.........................................................

(i)

16

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

(c)

(b)

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. W911NF13-1-0211). B.C.D. gratefully acknowledges a doctoral fellowship support from the Florida Space Grant Consortium (FSGC).

References 1. Marshall DB, Cox BN. 2008 Integral textile ceramic structures. Annu. Rev. Mater. Res. 38, 425–443. (doi:10.1146/annurev.matsci.38.060407.130214) 2. Cox BN et al. 2014 Stochastic virtual tests for high-temperature ceramic matrix composites. Annu. Rev. Mater. Res. 44, 479–529. (doi:10.1146/annurev-matsci-122013-025024) 3. Melenk JM, Babuška I. 1996 The partition of unity finite element method: basic theory and applications. Comp. Methods Appl. Mech. Eng. 139, 289–314. (doi:10.1016/S0045-7825 (96)01087-0) 4. Babuška I, Melenk JM. 1997 The partition of unity method. Int. J. Numer. Methods Eng. 40, 727–758. (doi:10.1002/(SICI)1097-0207(19970228)40:43.3.CO;2-E) 5. Moës N, Belytschko T. 2002 Extended finite element method for cohesive crack growth. Eng. Fract. Mech. 69, 813–833. (doi:10.1016/S0013-7944(01)00128-X)

.........................................................

5. Conclusion

17

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

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 8c). 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.

18 .........................................................

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

6. Moës N, Dolbow J, Belytschko T. 1999 A finite element method for crack growth without remeshing. Int. J. Numer. Methods Eng. 46, 131–150. (doi:10.1002/(SICI)1097-0207 (19990910)46:13.0.CO;2-J) 7. Daux C, Moës N, Dolbow J, Sukumar N, Belytschko T. 2000 Arbitrary branched and intersecting cracks with the extended finite element method. Int. J. Numer. Methods Eng. 48, 1741–1760. (doi:10.1002/1097-0207(20000830)48:123.0.CO;2-L) 8. Huynh DBP, Belytschko T. 2009 The extended finite element method for fracture in composite materials. Int. J. Numer. Methods Eng. 77, 214–239. (doi:10.1002/nme.2411) 9. Areias PMA, Belytschko T. 2005 Analysis of the 3D crack initiation and propagation using the extended finite element method. Int. J. Numer. Methods Eng. 63, 760–788. (doi:10.1002/nme.1305) 10. Belytschko T, Gracia R, Ventura G. 2009 A review of extended/generalized finite element methods for material modeling. Int. J. Numer. Methods Eng. 86, 637–666. (doi:10.1088/0965-0393/17/4/043001) 11. Wang YX, Waisman H. 2015 Progressive delamination analysis of composite materials using XFEM and a discrete damage zone model. Comput. Mech. 55, 1–26. (doi:10.1007/s00466-014-1079-0) 12. Zamani A, Gracie R, Eslami MR. 2012 Cohesive and noncohesive fracture by higher-order enrichment of XFEM. Int. J. Numer. Methods Eng. 90, 452–483. (doi:10.1002/nme.3329) 13. Song JH, Areias PMA, Belytschko T. 2006 A method for dynamic crack and shear band propagation with phantom nodes. Int. J. Numer. Methods Eng. 67, 868–893. (doi:10.1002/ nme.1652) 14. Van de Meer FP, Sluys LJ. 2009 A phantom node formulation with mixed mode cohesive law for splitting in laminates. Int. J. Fract. 158, 107–124. (doi:10.1007/s10704-009-9344-5) 15. Remmers JJC, de Borst R, Needleman A. 2003 A cohesive segments method for the simulation of crack growth. Comput. Mech. 31, 69–77. (doi:10.1007/s00466-002-0394-z) 16. Jager P, Steinmann P, Kuhl E. 2008 On local tracking algorithms for the simulation of threedimensional discontinuities. Comput. Mech. 42, 395–406. (doi:10.1007/s00466-008-0249-3) 17. Jager P, Steinmann P, Kuhl E. 2008 Modeling three-dimensional crack propagation—a comparison of crack path tracking strategies. Numer. Methods Eng. 76, 1328–1352. (doi:10.1002/nme.2353) 18. Chen BY, Pinho ST, De Carvalho NV, Baiz PM, Tay TE. 2014 A floating node method for the modelling of discontinuities in composites. Eng. Fract. Mech. 127, 104–134. (doi:10.1016/j.engfracmech.2014.05.018) 19. Ling DS, Bu LF, Tu FB, Yang QD, Chen YM. 2014 A finite element method with mesh-separation based approximation technique and its application in modeling crack propagation with adaptive mesh refinement. Int. J. Numer. Methods Eng. 99, 487–521. (doi:10.1002/nme.4689) 20. Ling DS, Fang XJ, Cox BN, Yang QD. 2011 Nonlinear fracture analysis of delamination crack jumps in laminated composites. J. Aerospace Eng. 24, 181–188. (doi:10.1061/ (ASCE)AS.1943-5525.0000008) 21. Ling DS, Yang QD, Cox BN. 2009 An augmented finite element method for modeling arbitrary discontinuities in composite materials. Int. J. Fract. 156, 53–73. (doi:10.1007/ s10704-009-9347-2) 22. Hansbo A, Hansbo P. 2004 A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Comput. Methods Appl. Mech. Eng. 193, 3523–3540. (doi:10.1016/j.cma.2003.12.041) 23. Iarve EV, Gurvich MR, Mollenhauer D, Rose CA, Davila CG. 2011 Mesh independent matrix cracking and delamination modelling in laminated composites. Int. J. Numer. Methods Eng. 88, 749–773. (doi:10.1002/nme.3195) 24. Iarve EV, Hoos KH, Mollenhauer D (eds). 2014 Damage initiation and propagation modeling in laminated composites under fatigue loading. In American Society for Composites 29th Technical Conf., 8–10 September 2014; La Jolla, CA, USA. 25. Prabhakar P, Waas A. 2013 A novel continuum-decohesive finite element for modeling in-plane fracture in fiber reinforced composites. Compos. Sci. Technol. 83, 1–10. (doi:10.1016/j.compscitech.2013.03.022) 26. Zhou ZQ, Fang XJ, Cox BN, Yang QD. 2010 The evolution of a transverse intra-ply crack coupled to delamination cracks. Int. J. Fract. 165, 77–92. (doi:10.1007/s10704-010-9506-5)

19 .........................................................

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

27. Fang XJ, Zhou ZQ, Cox BN, Yang QD. 2011 High-fidelity simulations of multiple fracture processes in a laminated composites in tension. J. Mech. Phys. Solids 59, 1355–1373. (doi:10.1016/j.jmps.2011.04.007) 28. Fang XJ, Yang QD, Cox BN, Zhou ZQ. 2011 An augmented cohesive zone element for arbitrary crack coalescence and bifurcation in heterogeneous materials. Int. J. Numer. Methods Eng. 88, 841–861. (doi:10.1002/nme.3200) 29. Liu W, Yang QD, Mohammadizadeh S, Su XY. 2014 An efficient augmented finite element method (A-FEM) for arbitrary cracking and crack interaction in solids. Int. J. Numer. Methods Eng. 99, 438–468. (doi:10.1002/nme.4697) 30. Liu W, Yang QD, Mohammadizadeh S, Su XY, Ling DS. 2013 An accurate and efficient augmented finite element method for arbitrary crack interactions. J. Appl. Mech. 80, 041033. (doi:10.1115/1.4007970) 31. Liu W, Schesser D, Yang QD, Ling DS. 2015 A consistency-check based algorithm for element condensation in augmented finite element methods for fracture analysis. Eng. Fract. Mech. 139, 78–97. (doi:10.1016/j.engfracmech.2015.03.038) 32. Yang QD, Thouless MD. 2001 Mixed mode fracture of plastically-deforming adhesive joints. Int. J. Fract. 110, 175–187. (doi:10.1023/A:1010869706996) 33. Yang QD, Cox BN. 2005 Cohesive models for damage evolution in laminated composites. Int. J. Fract. 133, 107–137. (doi:10.1007/s10704-005-4729-6) 34. Yang QD, Rosakis AJ, Cox BN. 2006 Dynamic fiber sliding along debonded frictional interface. Proc. R. Soc. A 462, 1081–1106. (doi:10.1098/rspa.2005.1602) 35. Hattiangadi A, Siegmund T. 2002 Bridging effects in cracked laminates under thermal gradients. Mech. Res. Commun. 29, 457–464. (doi:10.1016/S0093-6413(02)00300-2) 36. Hattiangadi A, Siegmund T. 2004 A thermomechanical cohesive zone model for bridged delamination cracks. J. Mech. Phys. Solids 52, 533–566. (doi:10.1016/S0022-5096(03)00122-4) 37. Benabou L, Sun Z, Dahoo PR. 2013 A thermo-mechanical cohesive zone model for solder joint lifetime prediction. Int. J. Fatigue 49, 18–30. (doi:10.1016/j.ijfatigue.2012.12.008) 38. Do BC, Liu W, Yang QD, Su XY. 2013 Improved cohesive stress integration schemes for cohesive zone elements. Eng. Fract. Mech. 107, 14–28. (doi:10.1016/j.engfracmech.2013.04.009) 39. Gasser TC, Holzapfel GA. 2006 3D Crack propagation in unreinforced concrete. A twostep algorithm for tracking 3D crack paths. Comput. Method Appl. Mech. Eng 195, 5198–5219. (doi:10.1016/j.cma.2005.10.023) 40. Oliver JHA. 2002 On strategies for tracking strong discontinuities in computational failure mechanics. In Proc. of the Fifth World Congress on Computational Mechanics (WCCM V) (eds HA Mang, FG Rammersdorfer, J Erberhardsteiner). Vienna, Austria. 41. Oliver J, Huespe AE, Samaniego E, Chaves EWV. 2004 Continuum approach to the numerical simulation of material failure in concrete. Int. J. Numer. Anal. Met. 28, 609–632. (doi:10.1002/nag.365) 42. Naderi M, Jung J, Yang QD. 2016 A three-dimensional augmented finite element for modeling arbitrary cracking in solids. Int. J. Fract. 192, 147–168. (doi:10.1007/s10704-016-0072-3) 43. Do BC. 2014 Extending the novel A-FEM to model arbitrary cracking in thermo-elastic solids. Coral Gables, FL: University of Miami. 44. Prasad NNV, Aliabadi MH, Rooke DP. 1994 Incremental crack growth in thermoelastic problems. Int. J. Fract. 66, R45–R50. (doi:10.1007/BF00042591) 45. Rangaraji S, Kokini K. 2004 Fracture in single-layer zirconia (YSZ)–bond coat alloy (NiCoCrAlY) composite coatings under thermal shock. Acta Mater. 52, 455–465. (doi:10.1016/j.actamat.2003.09.029) 46. Choules BD, Kokini K, Taylor TA. 2001 Thermal fracture of ceramic thermal barrier coatings under high heat flux with time-dependent behavior. Mater. Sci. Eng. A 299, 296–304. (doi:10.1016/S0921-5093(00)01393-9) 47. Gilbert A, Kokini K, Sankarasubramanian S. 2008 Thermal fracture of zirconia–mullite composite thermal barrier coatings under thermal shock: a numerical study. Surf. Coat. Technol. 203, 91–98. (doi:10.1016/j.surfcoat.2008.08.003) 48. Hutchinson JW, Suo Z. 1992 Mixed mode cracking in layered materials. Adv. Appl. Mech. 29, 63–191. (doi:10.1016/S0065-2156(08)70164-9) 49. He M-Y, Evans AG, Hutchinson JW. 1994 Crack deflection at an interface between dissimilar elastic materials: role of residual stresses. Int. J. Solids Struct. 31, 3443–3455. (doi:10.1016/0020-7683(94)90025-6)

20 .........................................................

rsta.royalsocietypublishing.org Phil. Trans. R. Soc. A 374: 20150282

50. Bale H, Blacklock M, Begley MR, Marshall DB, Cox BN, Ritchie RO. 2011 Characterizing three-dimensional textile ceramic composites using synchrotron X-ray micro-computedtomography. J. Am. Ceram. Soc. 95, 392–402. (doi:10.1111/j.1551-2916.2011.04802.x) 51. Rinaldi R, Blacklock M, Bale H, Begley MR, Cox BN. 2012 Generating virtual textile composite specimens using statistical data from micro-computed tomography: 3D tow representations. J. Mech. Phys. Solids 60, 1561–1581. (doi:10.1016/j.jmps.2012.02.008) 52. Yang QD, Rugg KL, Cox BN, Marshall DB. 2005 Evaluation of macroscopic and local strains in a 3D woven C/SiC composite. J. Am. Ceram. Soc. 88, 719–725. (doi:10.1111/ j.1551-2916.2005.00156.x) 53. Berbon MZ, Rugg KL, Dadkhah MS, Marshall DB. 2002 Effect of weave architecture on tensile properties and local strain heterogeneity in thin-sheet C-SiC composites. J. Am. Ceram. Soc. 85, 2039–2048. (doi:10.1111/j.1151-2916.2002.tb00401.x)