# Modelling thermal shock in functionally graded plates with finite element method.

1. IntroductionThe discovery of functionally gradient materials (FGMs) in the 1980s designed originally for an improvement of temperature resistance of constructive elements in nuclear reactors and chemical plants has stimulated a broad research activity for studying their behavior. An overview on mechanics of FGMs and their exploitation in modern engineering applications can be found in [1]. In general, FGMs are a new type of composite materials, which combine two or more constituent material phases. Unlike layered composite materials, FGMs are microscopically inhomogeneous with smoothly varying mechanical properties along one or more defined directions. Their structure is achieved by a gradual changing of the volume fraction of constituent materials [2].

The concept of FGMs allows one to establish the superior material properties compared with the constituent materials themselves. For instance, a typical metal/ceramic FGM possesses the mechanical strength of metal, but thermal resistance properties resemble ceramics [3]. These materials have been commonly used as thermal barrier coatings in engineering applications, undergoing a high temperature in-service environment. Moreover, they quite often experience an exposure to a very high temperature change in a very short period of time that is known as a thermal shock Such sever loading conditions result in high thermally induced stresses, which may very likely initiate a crack appearance in coatings. Therefore the primary interest of an analyst is to predict accurately a temperature field and is related to its temperature stresses at a design stage, since the latter have a critical relevance to trigger fracture mechanisms in FGM structures [4]. Hence, the goal of this research is to develop an efficient and reliable tool for studying thermal and mechanical response of a functionally graded metal/ceramic plate subjected to thermal shock.

A closed-form solution of the problem under consideration is very complicated and, in essence, it is possible in few one-dimensional cases with the simplest geometry, boundary conditions, and exponential gradation profiles for volume fractions, for example, [5-7]. While numerical techniques permit a look at more complex tasks under various boundary and loading conditions, material variations, including bidirectional FGMs, allow us to perform the crack growth analysis in FGM structures [8]. In this respect, the finite element method (FEM) is the most power tool [9]. So, we are going to use the FEM within the commercial engineering package ABAQUS [10] for carrying out thermomechanical and crack propagation analyses in a FGM plate under thermal shock.

The main issue encountered when applying the FEM to analyze graded materials is concerned with modelling their spatially dependent properties. The simplest way to do it involves the use of conventional homogeneous elements in successive layers of the mesh, where each layer contains its own material characteristics. This approach leads to a stepwise change in material properties in the direction of the material gradient. Those models have been already employed by a number of researchers and have enabled acquiring reasonable results, for example, [11-13]. However, this approach requires a fine mesh to achieve the accuracy; in turn it leads to an excessive computational cost. To overcome this disadvantage of homogeneous finite elements, gradient finite elements allowing the implementation of a gradient of material properties into a model at the element level so that the accuracy can be retained at a coarse mesh have been proposed in the literature. In [14], the authors have developed a two-dimensional (2D) graded finite element with material properties evaluated directly at the Gauss points of the element. An alternative graded element using a fully isoparametric element formulation with the shape functions the same as for displacement interpolations has been elaborated in [15].

In the present work, we use a two-dimensional plane strain graded finite element to model thermomechanical behavior and crack analysis of a metal/ceramic FGM plate. The element has been developed by programming appropriate user-defined subroutines within the ABAQUS software environment. An actual material gradation has been accomplished by sampling material quantities directly at the Gauss points of the element [15]. Moreover, for the sake of completeness of the study the finite element formulation underlying the coupled thermomechanical problem in 2D FGM structures is stated first. The virtual crack closure technique, used for modelling a crack growth under thermal loading, is incorporated into the finite element statement. A possible contact between crack lips during the crack advance is taken into account as well. Finally, the presented finite element framework is applied to provide an insight into the thermomechanical response and failure of a metal/ceramic FGM plate.

2. Governing Equations

Naturally the study of a structure behavior under temperature loading has to be conducted using a theory of thermomechanical problems. Here we assume that a plate is a 2D body made of an isotropic graded metal/ceramic material with the volume fraction varying along one selected direction as shown in Figure 1. Note that we neglect the crack existence inside the plate at the beginning of consideration.

Let the plate undergo small displacements and deformations, and it is under a plane strain state. Also we assume that the material behavior of the plate satisfies a linear elastic law. Thus, using a Lagrangian description, a thermoelasticity problem of an inhomogeneous 2D body A [member of] [R.sup.2] with boundary T [member of] [R.sup.1] at every point x = x(x, y) at a time instant t [member of] [0, [t.sub.*] can be stated as the following system of equations [16]:

[[sigma].sub.ji,j] + [b.sub.i] = [[rho][??].sub.i], in A x [0, [t.sub.*]], (1)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (2)

Here, we use notations usual in continuum mechanics and, in what follows, the subscript "," represents the partial differentiation, and the superscript "'" denotes the time derivative. The symbol [[sigma].sub.ij] stands for the components of Cauchy stress tensor, [u.sub.i] represents the displacements, [theta](x, t) = T - [T.sub.0] is a temperature change, where T(x, t) is an absolute temperature, while [T.sub.0] is a reference temperature at which the stress free state of the body takes place, [b.sub.i] are body forces, [q.sub.i] are the components of a surface heat flux, and R are internal heat generation sources. The material parameters such as the mass density [rho](x), the specific heat coefficient [c.sub.[upsilon]](x), and the stress-temperature modulus [beta](x) are functions of a spatial position in the Cartesian coordinate system.

In the equations of motion (1), the stress tensor satisfies the constitutive equations of thermoelasticity known as Duhamel-Neumann relations:

[[sigma].sub.ij] = [[lambda][epsilon].sub.kk][[delta].sub.ij] + 2[mu][[epsilon].sub.ij] - [beta][[delta].sub.ij][theta], (3)

where [lambda](x) and [mu](x) are the position dependent Lame constants and the linear relations between the stains [[epsilon].sub.ij] and the displacements [u.sub.i] are as follows:

[[epsilon].sub.ij] = 1/2([u.sub.i,j] + [u.sub.j,i]). (4)

Solving (3) for strains, we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (5)

where [alpha](x) is the coefficient of thermal expansion such that [beta](x) = (3[lambda] + 2[mu])[alpha] and [??] = -[lambda]/2[mu](3[lambda] + 2[mu]) and [??] = 1/4[mu]. The last term in (5) can be referred to as the thermal deformation which is the linear function of the temperature change [theta].

In the heat transfer equation (2), the heat flux is governed by the Fourier conduction law:

[q.sub.i] = -[kappa][[theta],.sub.i], (6)

where [kappa](x) is the coefficient of thermal conduction depending on the coordinates x.

To complete the initial-boundary value problem of thermoelasticity stated above, appropriate boundary and initial conditions have to be imposed. Given prescribed displacements [bar.[u.sub.i]] and surface loading [bar.[t.sub.i]], define the mechanical boundary conditions on the boundary T as follows:

[u.sub.i] = [bar.[u.sub.i]] on [T.sub.u] x [0, [t.sub.*], [[sigma].sub.ij][n.sub.j] = [bar.[t.sub.i]] on [T.sub.[sigma]] x [0, [t.sub.*]], (7)

where the corresponding boundary surfaces are such that [T.sub.u] [union] [T.sub.[sigma]] = T and [T.sub.u] [intersection] [T.sub.[sigma]] = [empty set], and [n.sub.j] is a unit outward normal to the surface T.

The thermal boundary conditions can be specified as a prescribed temperature [bar.T] and a prescribed heat flux involving a given surface heat flux [bar.q], volumetric heat flux R, and surface convection flux on appropriate nonoverlapping parts of the boundary T as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (8)

where h(x) is the heat film transfer coefficient and [T.sub.[infinity]] is an ambient temperature.

The initial conditions at t = 0 are assumed such that displacements [u.sub.i], velocities [[??].sub.i], and a temperature change [theta] are known functions:

[u.sub.i] (x, 0) = [bar.[u.sub.i.sup.0]] in A [union] T, [[??].sub.i] (x, 0) = [[??].sub.i.sup.0]] in A [union] T, [theta] (x, 0) = [[theta].sup.0]] in A [union] T. (9)

The system of (1)-(6) with boundary and initial conditions (7)-(16) are coupled because the equations of motion (1) contain a term defined by the temperature field through the thermal deformation (5), while the energy equation (2) is supplemented by a term corresponding to a rate of the dilatation of strains. Thereby, to find the solution, the system formulated above has to be solved for the displacements [u.sub.i] and the temperature field [theta] simultaneously.

Note that the dilatation term in the energy equation is important only for some type materials like rubber materials and it vanishes for most other materials when the external forces and heating are stationary. Besides, one can see that all material parameters involved in the constitutive laws depend on the coordinates. The latter feature considerably distinguishes the analysis of FGMs from the case of homogeneous materials.

3. Crack Modelling

If nonlinearities other than a crack propagation condition can be neglected, methods based on Linear Fracture Mechanics have been proven to be effective for crack modelling. In this respect, the virtual crack closure technique (VCCT) in conjunction with the FEM is one of the most commonly used approaches for determining the components of strain energy release rate (SERR) along the crack front [17]. This approach has been already successfully used for analyzing a crack propagation in layered composites, for example, in [18, 19], and has been considered as a suitable method for predicting a crack growth in FGMs, for example, [20, 21], because no assumptions of isotropy or homogeneity around the crack are necessary in the method.

The VCCT requires the introduction into a model of a predefined geometrical discontinuity. Let the crack of length a exist in the FGM plate as shown in Figure 1. In accordance with the Griffiths crack growth criterion, we assume that crack starts its advance, when the SERR G at the crack tip exceeds the critical strain energy release rate known as fracture toughness Gc at any point of the material. Using a 2D finite element model of the plate discretized with four-node elements, the work needed to extend (or that is equivalent to close) a crack along one element length can be defined with the VCCT as follows:

[DELTA]W = 1/2 ([Y.sub.1.sup.i][DELTA][u.sub.1.sup.j] + [Y.sub.2.sup.i][DELTA][u.sub.2.sup.j]). (10)

Here as shown in Figure 2(a), [Y.sub.1.sup.i] and [Y.sub.2.sup.i] are the shear and normal forces at the ith node, and [DELTA][u.sub.1.sup.j] and [DELTA][u.sub.2.sup.j] refer to the shear and opening displacements at the jth node. This work done by the reactive forces and, consequently, the SERR G = [DELTA]W/[DELTA]A, where [DELTA]A = [DELTA]a x b is an area of a new surface created by a crack extension [DELTA]a, can be calculated in one-step finite element analysis. Then the separation of individual components of the mixed-mode SERR in (10) can be found via the following equations [19]:

[G.sub.I] = 1/2bd[Y.sub.2.sup.i]([u.sub.2.sup.j"] - [u.sub.2.sup.j']), (11)

[G.sub.II] = 1/2bd[Y.sub.1.sup.i]([u.sub.1.sup.j"] - [u.sub.1.sup.j']), (12)

where b is the plane strain/stress thickness and d is the length of the 2D element in the mesh at the crack front. This value coincides with the crack extension [DELTA]a. Then, the crack will propagate due to mode I loading conditions, if the force and displacement are achieved critical values as illustrated in Figure 2(b).

In the case of mixed-mode loading, the node at the crack tip is released, if the criterion of the equivalent SERR such as [G.sub.eq] [less than or equal to] [G.sub.eq.sup.cr] was met. The total SERR [G.sub.T] = [G.sub.I] + [G.sub.II] + [G.sub.III], where Gjjj = 0 in the 2D case, is often considered for this purpose. In the forthcoming calculations, we will use the power mixity law as follows [10]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (13)

Thereby, the finite element model should be provided by the test data such as fracture toughness modes [G.sub.IC], [G.sub.IIC], and [G.sub.IIIC], and the fitting coefficients [k.sub.1], [k.sub.2], and [k.sub.3].

If contact between crack lips occurs, appropriate contact conditions along the interface of their interactions [T.sub.c] [member of] A should be defined. Then the impenetrability and friction constraints can be imposed by establishing the KarushKuhn-Tucker conditions on [T.sub.c] = [T.sub.c.sup.-] [union] [T.sub.c.sup.+] in terms of displacements and tractions in the following form [22]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (14)

In the equations above we denote as in [22, 23] that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is a gap function defining minimal distance between a point of the slave surface x [member of] [T.sub.c.sup.-] and its orthogonal projection on the master surface [bar.x] [member of] [T.sub.c.sup.+] with a unit normal [bar.[n.sup.+]]. Similarly, [g.sub.T] = [[[u.sub.1]].sub.[tau]] is a slip function describing a relative movement of those points on the contacting lines in the tangential direction [tau] = [e.sub.3] x [bar.[n.sup.+]] in Figure 3. Also [t.sub.N] and [t.sub.T] are normal and shear scalar parameters of the contact traction vector [t.sub.c] = [t.sub.N][bar.[n.sup.+]] + [t.sub.T][tau], whereas [[tau].sub.cr] is a threshold of the tangential contact traction up to which the lines are retained together; otherwise the slip between the lines is governed by Coulomb's law of friction as [[tau].sub.cr] = [[gamma]t.sub.N] with a given coefficient of friction [gamma].

4. Weak Form of the IBVP

To obtain a numerical solution of the initial-boundary value problem (IBVP) formulated in the previous sections, its weak or variational form should be stated. The weak form of the IBVP can be obtained by applying Hamilton's principle to the body under consideration. Then, the coupled thermoelastic problem of the cracked FGM body in the current configuration corresponds to the following variational equalities:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (15)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (16)

where the test functions [delta][u.sub.i] and [delta][theta] reside in appropriate vector spaces of kinematically admissible displacements U and the admissible temperature field V such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (17)

Note that in (15) the last term comes from the contribution of the contact traction caused by the discontinuity in the displacement fields given on the boundary [T.sub.c]. If the contact conditions (14) are satisfied exactly, the contact term in (15) adds nothing to the total energy. However, it is possible only by a true solution, but it is not necessary by the arbitrary test functions. Because the impenetrably dictates fulfilling the inequality [delta][g.sub.N][t.sub.N] [greater than or equal to] 0, whence, the variational equality (15) can be reorganized into a variational inequality as in [22]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (18)

subject to constraints (14). The final form of the functional (15) depends on the method which is used to impose the contact constraints. For this purpose, several different techniques are available in the literature such as Lagrange multipliers and penalty method [22]. Also the contact integral provides additional conditions to which the stress-strain state of the body has to obey. In doing so, at each moment of time, the stress state (or reaction forces), recovered from the calculated displacement solutions, is checked on the fulfilling the fracture criterion (13). Here we assume that the existing crack does not cause a discontinuity in the temperature field. Thus, functional (16) is not affected by the contact conditions anyhow.

5. Finite Element Formulation

The variational problem (15)-(16) can be solved using the FEM. So, the plate domain A is discretized by a number of nonoverlapping elements [A.sup.(e)] such that the union of the elements comprises the total domain; that is, A = [[union].sub.e=1.sup.NE] [A.sup.(e)]. For each finite element, unknown displacements and temperature field can be approximated by functions residing in the finite dimensional counterparts of the infinite vector spaces U and V as follows:

[u.sub.i.sup.(e)] = [N.sup.I] (x)[U.sub.i.sup.I](t), [[theta](.sup.e)] = [[??].sup.P] (x)[[THETA].sup.P](t), (19)

where I, P = 1, 2, ..., [n.sup.(e)] are nodal points and [n.sup.(e)] is a number of nodes in the basic element, [N.sup.I](x) and [[??].sup.P](x) are shape functions associated with nodes I and P, and [U.sub.i.sup.I](t) and [[THETA].sup.P](t) are nodal unknown displacements and nodal unknown temperature field, respectively. The summation over the repeated index I and P is used.

Using the standard Galerkin approach, the finite element equation can be formulated by substituting interpolations (19) and their variations into the system of variational equations (15) and (16). Then, the system of semidiscrete finite element equations can be expressed as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (20)

where [M.sub.uu], [K.sub.uu], [K.sub.[theta][theta]] , and [C.sub.[theta][theta]] are global mass, stiffness, conductivity, and capacity matrices, [K.sub.u[theta]] and [C.sub.[theta]u] are coupled matrices, U and [THETA] are global vectors of unknown displacements and temperature and their appropriate time derivatives, and [F.sub.u], [F.sub.[theta]], and [F.sub.cont] are global vectors of external mechanical, thermal, and contact forces. Herein an explicit form of the vector of nodal contact forces [F.sub.cont] is constructed depending on the method chosen to implement the contact constraints into (15). It needs also to note that in (20) both the vector of unknowns and the contact forces should be found. Moreover, contact surfaces, on which contact constraints are enforced, are not known a priori as well. Thus (20) is a system of high nonlinear partial differential equations.

Discretizing a time domain t = [[union].sub.n=0.sup.N-1][[t.sup.n], [t.sup.n+1]], where it is assumed that a solution of (20) exists, and assuming an implicit time integration scheme to step forward the solution from a time instant t = [t.sup.n] to an instant t + [DELTA]t = [t.sup.n+1], we come to a nonlinear system of algebraic equations expressed as [22]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (21)

where M, C, and K are partitioned matrices conditionally referring mass, damping, and stiffness matrices, respectively, combining respective matrices of mechanical and thermal parts; d = [{[U.sub.u], [U.sub.[theta]]}.sup.T] is a combined vector of global nodal displacements and temperatures and the dots over this vector denote its appropriate time derivatives; F is a global thermal-mechanical force vector and [F.sub.cont] is a global vector of contact forces.

The Newton-Raphson method reduces (21) to a linearized form, which at each time step [DELTA]t has the following form:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (22)

In (22), [[??].sub.IJ] are submatrices of the coupled Jacobian matrix, [DELTA][U.sub.u] and [DELTA][U.sub.[theta]] are corrections referring to the solution vectors [U.sub.u] and [U.sub.[theta]], accordingly, and [R.sub.u] and [R.sub.[theta]] are components of the residual force vector R which is defined from (21) as follows [9]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (23)

where the superscript "n" represents a discrete time instant and the subscript "i" denotes an iteration.

It is important to notice that the system Jacobian matrix is nonsymmetric. Therefore, the system of linear algebraic equations (22) should be solved using a nonsymmetric matrix storage and a nonsymmetric solution scheme at each iteration within each time step.

6. Finite Element Modelling

The finite element formulation stated in previous sections is employed to develop the model of a cracked FGM plate, shown in Figure 1 within the engineering FEM package ABAQUS/Standard. As said earlier, one of basic issues referring to the use of the FEM in simulations of structures made of inhomogeneous materials is to model a gradient of material properties. In the present paper, we provide an effective way for implementing the continuous variation of thermal and mechanical properties into the model of FGM plate by developing a graded finite element. This is achieved by programming appropriate user-defined subroutines for finite elements available in the ABAQUS library. As a base element, a quadratic eight-node plane strain temperature-displacement element (CPE8T) is chosen. The element assumes a biquadratic displacement interpolation but a bilinear temperature interpolation. Either full or reduced integration over the element area is allowed.

The material subroutine UMAT is used to incorporate gradients into the element's Young's modulus, Poisson's ratio, and coefficient of thermal expansion, whereas the subroutine UMATHT is applied to assign a spatial variation to the other element's parameters such as conductivity and specific heat. The gradient of the element mass density is defined by the subroutine USDFLD. It should be noticed that all the subroutines specify the given material properties by sampling them directly at the Gauss points of the element.

A sequential temperature-stress analysis is carried out to model steady distributions of the temperature field and the related thermal stresses. In this analysis, the temperature field is firstly calculated by solving the heat equation and, then, the known temperature is inputted into the mechanical equation to compute thermally induced stresses. Thus, the mechanical and thermal equations in (20) are solved consequently and independently of each other. This procedure is numerically inexpensive. Once the steady thermomechanical solution is known, a transient heat transfer in conjunction with fracture analysis should be undertaken for examining the crack growth. A fully coupled temperature-stress analysis is carried out in this case. This means that both the heat equation and the motion equation in (20) are solved simultaneously with the same size of the time increment. Moreover, thermally induced stresses calculated in the transient analysis are used as driving forces for advancing the existing crack. Hence, the fracture criterion (13) is checked at each time step of the analysis. If the fracture criterion was met, the algorithm of the VCCT runs and crack is growing at this time instant. Herewith, a contact algorithm tracks a contact status and if contact occurs, the impenetrability of crack lips into each other is provided by the algorithm. The frictionless contact conditions are realized in ABAQUS by using the general contact scheme. The surface-to-surface contact behavior is governed by the "hard contact" model within the small displacement tracking algorithm [10]. The contact constraints are implemented via the penalty algorithm. Another contact model suitable for the crack propagation analysis can be found, for example, in [24, 25].

The VCCT does not require any assumptions of the form of the stresses and displacements around the crack tip [17]. Therefore, any singularity elements are not required to be embedded into the crack tip area. However, special crack tip elements have been proposed in the literature [18]. This allows one to minimize computer efforts for inducing accurately 1/ [square root of r]-singularity of the stress field at the crack tip because of mesh refinement only in a local small region. Thus, in the fracture analysis to calculate the static thermal stress intensity factor (TSIF), we employ a crack tip two-dimensional element with quarter-point nodes obtained by collapsing one side of the rectangular element [26]. A typical crack tip rosette of singular quarter-point finite elements used to induce the required singularity is illustrated in Figure 4.

The radius of ring of the singular elements around the crack tip is taken as 0.05a in the calculations. The TSIFs are computed by using the J-integral analysis employing the domain integral method [10]. For the sake of comparisons we applied the VCCT with and without the singular element to calculate the TSITs as well. In this case the node release option of the VCCT was suppressed in the calculations. Herewith the equations for the eight-node singular elements given in [27] have been used for evaluating mode I and mode II components of the SERR instead of formulae (11) and (12), as shown in Figure 4. From these results, the element size producing the required singularity in the mesh with regular elements only was found.

In fracture propagation analysis with the VCCT, we do not use the singular elements in the mesh to avoid a very time consuming procedure of remeshing for tracking a crack growth. However, to keep the required singularity properties which follow from the TSIF analysis with the singular elements, the functionally graded elements used in the analysis were refined up to a characteristic length of 0.01a in a small region containing the predefined crack path. The size of remaining elements in the mesh was gradually increasing toward the top and bottom edges of the plate, keeping the symmetry, as shown in Figure 5.

Although the crack propagation is, in essence, a dynamic phenomenon itself [28] and the dynamic solution of the problem is more appropriate, we assume a quasi-static character of the crack growth in all forthcoming predictions. This was done to simplify the solution procedure of the complex elastodynamic task stated as a coupled thermomechanical problem in the previous sections. Thereby, we will neglect the inertia effect on the crack evolution in the calculations. Besides, the strain dilatation term is not taken into account in the calculations as well. It means that the coupling between thermal and mechanical parts is assumed to be through the temperature field only. Moreover, for the sake of simplicity, we suppose that there is no variation of fracture toughness along the material gradation in thermal cracking. The value of fracture toughness is spatially independent and is defined by toughness of ceramics. It can be interpreted as a lower bound of the thermal cracking resistance of the FGM plate.

7. Numerical Results

In this section, the results of two simulations of the FGM plate shown in Figure 1 are presented. The first analysis concerns a steady state heat transfer problem in the FGM plate of the dimensionless width w of 1 and length 2/ of 2 under quasistatic thermal loading. It is assumed that the plate is made of a functional graded super alloy (Rene-41)/zirconia material. The thermal and mechanical properties of the constituents of the material are given in Table 1 as those in [6]. There, the mechanical parameters are the following functions of the x-coordinate:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (24)

where the subscripts "c" and "m" refer to ceramic and metal components, respectively.

It is supposed that the plate is under a stress-free state at the reference temperature [T.sub.0] and is subjected to the surface temperatures [T.sub.1] and [T.sub.2] at x = 0 and x = w, respectively. Because of the symmetry, a half plate is modelled only. The crack of length a/w = 0.2 is embedded into the plate as shown in Figure 6(a). The radius of ring of the singular elements around the crack tip is taken as 0.05a in the following calculations and is illustrated in Figure 6(b). Two conditions of thermal loading and groups of material parameters are considered in this analysis.

Figure 7(a) shows comparisons of the TSIFs normalized by [[sigma].sub.0] = [E.sub.c][[alpha].sub.c][T.sub.0]/(1 - v) between the analytical solutions in [29] and results calculated numerically with the graded elements in the cracked FGM Rene-41/zirconia plate. The material properties of the plate are taken as [E.sub.m]/[E.sub.c] = 5, [[alpha].sub.m]/[[alpha].sub.c] = 2, and [[kappa].sub.m]/[[kappa].sub.c] = 1. The plate undergoes a uniform dimensionless temperature [T.sub.1] = [T.sub.2], which is less than T0 = 10 and corresponds to a variety of values [T.sub.1]/[T.sub.0] such as 0.5, 0.2, 0.1, and 0.005.

Analogously, the normalized TSIFs found analytically in [29] and calculated numerically with the graded elements are compared in Figure 7(b) for the cracked FGM Rene-41/zirconia plate with the material properties as [E.sub.m]/[E.sub.c] = 10, [[alpha].sub.m]/[[alpha].sub.c] = 2, and [[kappa].sub.m]/[[kappa].sub.c] = 10, which is undergoing a uniform dimensionless temperature [T.sub.2] = 0.5[T.sub.0], where [T.sub.0] = 10 and the ratio [T.sub.1]/[T.sub.0] is equal to 0.5, 0.2, 0.1, and 0.005. In Figure 7, the lines refer to the analytic solutions, whereas the markers indicate the numerical results. One can see an excellent agreement between both solutions for two cases of thermal loading and material parameters.

In the other analysis, a plate of width w of 10 mm and length 2l of 40 mm made of metal/ceramic functionally graded material is modelled. It is assumed that the material of the plate consists of particles of titanium alloy Ti-6Al-4V and zirconium dioxide Zr[O.sub.2] ceramic. Mechanical properties of the constituent materials are listed in Table 2 and are those as given in [30].

Also the authors there supposed that the volume fraction of the metal phase in the direction of x-coordinate, as shown in Figure 1, varies in accordance with a power function:

[V.sub.m] = [(x/w).sup.P]. (25)

From (25), it follows that the FGM is rich in ceramic when the parameter p > 1 and rich in metal when the parameter p < 1. The effective properties of the metal/ceramic Ti-6Al-4V/Zr[O.sub.2] plate with negligible porosity have been estimated by means of the rule of mixtures of a two-phase material as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (26)

It is worth noticing that the rule of mixtures similar to (26) is only the simplest approximation of the effective properties of a nonhomogeneous material. The reasons of inaccuracy of this approach have been widely discussed earlier, for example, in [31, 32]. Instead, homogenization techniques based on microstructural analysis and empirical observations should be used for more accurate predictions in material properties of FGMs as shown in [2]. However, this discussion is out of the scope of the present paper. Only material parameters known in the literature are used in the calculations below.

As a thermal load, one cycle of thermal shock is applied to the metal/ceramic FGM plate, as shown in Figure 8. The plate is initially assumed to be stress-free at the temperature [T.sub.0]; then the plate is heated on the ceramic surface by a high temperature [T.sub.H], while the metallic surface remains at a low temperature [T.sub.L] equal to [T.sub.0]. After that, when the temperature field within the plate reaches a steady state, the ceramic surface undergoes a fast cooling due to an intensive forced convection. The temperature of the ceramic surface quickly decreases to [T.sub.L]. As known this cooling step promotes a propagation of an existing crack that we are aiming to investigate numerically herein.

In the analysis, we assume that [T.sub.0] = 300 K, [T.sub.H] = 1300 K, and [T.sub.L] = 300 K. The film convection coefficient is h = 2000 W/[mm.sup.2] K and the temperature of a surrounding medium is [T.sub.[infinity]] = 300 K. The mesh contains a total of 2047 two-dimensional graded finite elements with 14730 numbers of degrees of freedom in the model (Figure 5). A preexistence crack of the relative length a/w of 0.2 is modelled by an actual small gap between the finite elements located along the crack lips. The fracture criterion (13) with the input parameters taken as [G.sub.IC] = [G.sub.IIC] = [G.sub.IIIC] = 0.19 N/[mm.sup.2] and [k.sub.1] = [k.sub.2] = [k.sub.3] = 1 is exploited in the fracture analysis.

The loading conditions dictate the types of analyses that should be performed in this study. First a steady-state thermomechanical analysis will be used to simulate a heating phase. Second, a transient thermomechanical analysis, which is accompanied by estimations of crack onset and propagation, will be applied to describe cooling. The contact analysis is performed in both cases once contact is detected.

The contour plots of the distributions of steady state and transient temperatures in the cracked FGM plate with p = 0.5 are illustrated in Figures 9(a) and 9(b), respectively. One can see that due to heating the highest temperature is achieved on the ceramic surface of the plate. However, as a result of cooling, the temperature is redistributed with time such that its highest value is reached in the central region of the plate. Hence, an initial high level of thermally induced stresses in the ceramic surface arising during heating will be followed by their fast redistributions within the plate width under cooling.

For the sake of demonstration of a stress state induced by the corresponding temperature fields due to heating and cooling phases shown in Figure 9, the distributions of Mises stresses are displayed in Figure 10. One can see that the stress state under heating phase is more uniform than that in cooling. The fast temperature dropping gives rise to stress redistribution inside the plate; as a result a pronounced nonuniformity of the Mises stress occurs. This clearly proves a dangerous character of cooling from the standpoint of strength and the need to examine a fracture resistance of FGM plates under thermal shock conditions.

To provide a deeper insight into the thermal cracking behavior in the FGM plate, the stresses caused by the transient temperature distributions under cooling are further examined. The series of contour plots of the distributions of transient thermal stress [[sigma].sub.yy] at different nondimensional moment of time [tau] equal to a ratio of a current simulation time to a total analysis time are presented in Figure 11. From the calculations carried out with the proposed finite element model, one can conclude that crack opening due to cooling occurs because the regions of compressive stresses inside the plate appear. These compressive stresses act as a bending moment ahead of a crack tip resulting in tensile stresses in a small region near the crack lips. If these stresses are high enough to keep driving the crack, then the crack can become self-driven. Hence, as soon as the TSIF induced by the temperature field exceeds the critical stress intensity factor of the FGM, the crack grows. The crack stops because compressive stresses appear in the domain of the plate behind the crack tip.

8. Conclusions

A finite element formulation of the coupled thermomechanical problem in a functionally graded plate undergoing plane strain conditions is presented. A fully coupled thermal stress analysis in conjunction with the VCCT is used to examine the thermal crack growth in the FGM plate with a preexistence crack. The finite element analysis is carried out using a graded finite element developed within the ABAQUS code via the combination of user-defined subroutines. With the graded element, a material gradation of mechanical and thermal properties of the FGM is incorporated into the finite element model. Both the sequential thermomechanical analysis and the coupled thermomechanical analysis associated with crack extension procedure and contact conditions are carried out. The obtained results showed that the proposed finite element model based on the graded finite element is effective in thermomechanical simulations of FGM plates and enables performing accurate predictions of thermal and mechanical behavior of FGM plates under thermal shock.

http://dx.doi.org/10.1155/2016/7514638

Competing Interests

The author declares that he has no competing interests.

Acknowledgments

The author acknowledges the Erasmus Mundus exchange program ACTIVE (Grant no. AC/TG2/SOTON/PD/23/2015) for the support of his stay in the University of Southampton. Also he would like to thank Professor Atul Bhaskar for fruitful discussions.

References

[1] V. Birman and L. W. Byrd, "Modeling and analysis of functionally graded materials and structures," Applied Mechanics Reviews, vol. 60, no. 1-6, pp. 195-216, 2007.

[2] Y. Miyamoto, W. A. Kaysser, B. H. Rabin, A. Kawasaki, and R. G. Ford, Functionally Graded Materials: Design, Processing and Applications, Springer, New York, NY, USA, 1999.

[3] M. S. El-Wazery and A. R. El-Desouky, "A review on functionally graded ceramic-metal materials" Journal of Materials and Environmental Science, vol. 6, no. 5, pp. 1369-1376, 2015.

[4] T. Sadowski, "Non-symmetric thermal shock in ceramic matrix composite (CMC) materials" Solid Mechanics and Its Applications, vol. 154, pp. 99-148, 2009.

[5] Y.-D. Lee and F. Erdogan, "Residual/thermal stresses in FGM and laminated thermal barrier coatings," International Journal of Fracture, vol. 69, no. 2, pp. 145-165, 1994.

[6] F. Erdogan and B. H. Wu, "Crack problems in FGM layers under thermal stresses" Journal of Thermal Stresses, vol. 19, no. 3, pp. 237-265, 1996.

[7] Z.-H. Jin and G. H. Paulino, "Transient thermal stress analysis of an edge crack in a functionally graded material" International Journal of Fracture, vol. 107, no. 1, pp. 73-98, 2001.

[8] I. V. Ivanov, T. Sadowski, and D. Pietras, "Crack propagation in functionally graded strip under thermal shock" The European Physical Journal: Special Topics, vol. 222, no. 7, pp. 1587-1595, 2013.

[9] T. Belytschko, W. K. Liu, and B. Morgan, Nonlinear Finite Elements for Continua and Structures, John Wiley & Sons, Chester, UK, 2000.

[10] ABAQUS 6.12 User's Manuals, Dassault Systemes Simulia Corp., Providence, RI, USA, 2012.

[11] R. L. Williamson, B. H. Rabin, and J. T. Drake, "Finite element analysis of thermal residual stresses at graded ceramic-metal interfaces. Part I. Model description and geometrical effects" Journal of Applied Physics, vol. 74, no. 2, pp. 1310-1320, 1993.

[12] G. Anlas, M. H. Santare, and J. Lambros, "Numerical calculation of stress intensity factors in functionally graded materials," International Journal of Fracture, vol. 104, no. 2, pp. 131-143, 2000.

[13] J. C. Lee, J. H. Park, S. H. Ryu et al., "Reduction of functionally graded material layers for [Si.sub.3][N.sub.4]-[Al.sub.2][O.sub.3] system using three-dimensional finite element modeling," Materials Transactions, vol. 49, no. 4, pp. 829-834, 2008.

[14] M. H. Santare and J. Lambros, "Use of graded finite elements to model the behavior of nonhomogeneous materials," Journal of Applied Mechanics, vol. 67, no. 4, pp. 819-822, 2000.

[15] J.-H. Kim and G. H. Paulino, "Isoparametric graded finite elements for nonhomogeneous isotropic and orthotropic materials," Journal of Applied Mechanics, Transactions ASME, vol. 69, no. 4, pp. 502-514, 2002.

[16] R. B. Hetnarski and M. R. Eslami, Thermal Stresses--Advanced Theory and Applications, Springer Science & Business Media, 2009.

[17] E. F. Rybicki and M. F. Kanninen, "A finite element calculation of stress intensity factors by a modified crack closure integral," Engineering Fracture Mechanics, vol. 9, no. 4, pp. 931-938, 1977.

[18] R. Krueger, "Virtual crack closure technique: history, approach, and applications," Applied Mechanics Reviews, vol. 57, no. 2, pp. 109-143, 2004.

[19] V N. Burlayenko and T. Sadowski, "FE modeling of delamination growth in interlaminar fracture specimens," Budownictwo i Architektura, vol. 2, no. 1, pp. 95-109, 2008.

[20] H. L. Zhou, "Implementation of crack problem of functionally graded materials with ABAQUS[TM]," Advanced Materials Research, vol. 284-286, pp. 297-300, 2011.

[21] V N. Burlayenko, H. Altenbach, T. Sadowski, and S. D. Dimitrova, "Computational simulations of thermal shock cracking by the virtual crack closure technique in a functionally graded plate," Computational Materials Science, vol. 116, pp. 11-21, 2016.

[22] T. A. Laursen, Computational Contact and Impact Mechanics: Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis, Springer, Berlin, Germany, 2002.

[23] V. N. Burlayenko and T. Sadowski, "Finite element nonlinear dynamic analysis of sandwich plates with partially detached facesheet and core," Finite Elements in Analysis and Design, vol. 62, pp. 49-64, 2012.

[24] V N. Burlayenko and T. Sadowski, "Modeling the dynamic debonding growth of sandwich plate," in Proceedings of the 4th International Conference on Nonlinear Dynamics, Yu. V Mikhlin and M. V Perepelkin, Eds., pp. 225-230, Tochka, Sevastopol, Ukraine, June 2013.

[25] V N. Burlayenko and T. Sadowski, "Simulations of post-impact skin/core debond growth in sandwich plates under impulsive loading," Journal of Applied Nonlinear Dynamics, vol. 3, no. 4, pp. 369-379, 2014.

[26] R. S. Barsoum, "Triangular quarter-point elements as elastic and perfectly-plastic crack tip elements," International Journal for Numerical Methods in Engineering, vol. 11, no. 1, pp. 85-98, 1977

[27] I. S. Raju, "Calculation of strain-energy release rates with higher order and singular finite elements," Engineering Fracture Mechanics, vol. 28, no. 3, pp. 251-274, 1987

[28] L. B. Freund, Dynamic Fracture Mechanics, Cambridge Monographs on Mechanics and Applied Mathematics, Cambridge University Press, 1990.

[29] J.-H. Kim and G. H. Paulino, "Mixed-mode fracture of orthotropic functionally graded materials using finite elements and the modified crack closure method," Engineering Fracture Mechanics, vol. 69, no. 14-16, pp. 1557-1586, 2002.

[30] T. Fujimoto and N. Noda, "Influence of the compositional profile of functionally graded material on the crack path under thermal shock," Journal of the American Ceramic Society, vol. 84, no. 7, pp. 1480-1486, 2001.

[31] R. Hill, "Elastic properties of reinforced solids: some theoretical principles," Journal of the Mechanics and Physics of Solids, vol. 11, no. 5, pp. 357-372, 1963.

[32] L. J. Walpole, "On bounds for the overall elastic moduli of inhomogeneous systems--I," Journal of the Mechanics and Physics of Solids, vol. 14, no. 3, pp. 151-162, 1966.

Vyacheslav N. Burlayenko

Department of Applied Mathematics, National Technical University "KhPI", 21 Frunze Street, Kharkiv 61002, Ukraine

Correspondence should be addressed to Vyacheslav N. Burlayenko; burlayenko@yahoo.com

Received 15 April 2016; Revised 24 May 2016; Accepted 2 June 2016

Academic Editor: Sergei V. Panin

Caption: Figure 1: A functionally graded plate.

Caption: Figure 2: The VCCT: (a) for four-node elements and (b) with a crack growth criterion.

Caption: Figure 3: Two-dimensional contact constraints.

Caption: Figure 4: A crack tip rosette of singular quarter-point finite elements.

Caption: Figure 5: A finite element mesh used in the crack propagation analysis.

Caption: Figure 6: A model of the cracked plate.

Caption: Figure 7: Comparison of the normalized TSIF for the cracked FGM plate under thermal loading: (a) [T.sub.1] = [T.sub.2] = T and T < [T.sub.0]; and (b) [T.sub.2] = 0.5[T.sub.0].

Caption: Figure 8: A scheme of thermal shock.

Caption: Figure 9: Contour plots of the temperature distributions: (a) steady state under heating and (b) transient state due to cooling at the end of analysis.

Caption: Figure 10: Contour plots of the distributions of Mises stress: (a) under heating and ((b) and (c)) due to cooling at a moment of a small period of time after the beginning and at the end of transient analysis.

Caption: Figure 11: Contour plots of transient stress [[sigma].sub.yy] at different moments of time: (a) scale level, (b)[tau] = 0.1 (c)[tau] = 0.5, and (d) r = 1.0.

Table 1: Material properties of Rene-41 and Zr[O.sub.2] [6] Property Constituents Rene-41 Zr[O.sub.2] Young's modulus E, GPa 219.7 151 Mass density [rho], kg/[m.sup.3] 8250 5680 Poisson's ratio v 0.33 0.33 Coefficient of thermal expansion 16.7 10 [alpha], x [10.sup.-6] 1/K Thermal conductivity [kappa], W/mK 25.53 2.09 Table 2: Material properties of Ti-6Al-4V and Zr[O.sub.2] [28]. Property Constituents Ti-6Al-4V Zr[O.sub.2] Young's modulus E, GPa 66.2 117 Poisson's ratio v 0.32 0.333 Mass density p, kg/[m.sup.3] 4420 5600 Coefficient of thermal expansion 10.3 7.11 [alpha], x [10.sup.-6] 1/K Thermal conductivity [kappa], W/mK 18.1 2.036 Specific heat [c.sub.[upsilon], J/kgK 808.3 615.6 Fracture toughness [K.sub.IC] 60 6

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Research Article |
---|---|

Author: | Burlayenko, Vyacheslav N. |

Publication: | Advances in Materials Science and Engineering |

Article Type: | Report |

Date: | Jan 1, 2016 |

Words: | 7529 |

Previous Article: | Dynamic behaviours of a single soft rock-socketed shaft subjected to axial cyclic loading. |

Next Article: | Influence of pore size on the optical and electrical properties of screen printed Ti[O.sub.2] thin films. |

Topics: |