Abstract
Thermomechanical behavior and crack propagation in a functionally graded metal/ceramic plate undergoing thermal shock are analyzed by using the finite element method. A two-dimensional plane strain functionally graded finite element has been developed within the ABAQUS software environment for this purpose. An actual material gradation has been accomplished by sampling material quantities directly at the Gauss points of the element via programming appropriate user-defined subroutines. The virtual crack closure technique is used to model a crack growth under thermal loading. Contact possible between crack lips during the crack advance is taken into account in thermomechanical simulations as well. The paper shows that the presented finite element model can be applied to provide an insight into the thermomechanical respond and failure of the metal/ceramic plate.
1. Introduction
The 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 with boundary at every point at a time instant can be stated as the following system of equations [16]: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 stands for the components of Cauchy stress tensor, represents the displacements, is a temperature change, where is an absolute temperature, while is a reference temperature at which the stress free state of the body takes place, are body forces, are the components of a surface heat flux, and are internal heat generation sources. The material parameters such as the mass density , the specific heat coefficient , and the stress-temperature modulus 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:where and are the position dependent Lamé constants and the linear relations between the stains and the displacements are as follows: Solving (3) for strains, we getwhere is the coefficient of thermal expansion such that and and . The last term in (5) can be referred to as the thermal deformation which is the linear function of the temperature change .
In the heat transfer equation (2), the heat flux is governed by the Fourier conduction law:where 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 and surface loading , define the mechanical boundary conditions on the boundary as follows: where the corresponding boundary surfaces are such that and , and is a unit outward normal to the surface .
The thermal boundary conditions can be specified as a prescribed temperature and a prescribed heat flux involving a given surface heat flux , volumetric heat flux , and surface convection flux on appropriate nonoverlapping parts of the boundary as follows:where is the heat film transfer coefficient and is an ambient temperature.
The initial conditions at are assumed such that displacements , velocities , and a temperature change are known functions:
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 and the temperature field 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 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 at the crack tip exceeds the critical strain energy release rate known as fracture toughness 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:Here as shown in Figure 2(a), and are the shear and normal forces at the th node, and and refer to the shear and opening displacements at the th node. This work done by the reactive forces and, consequently, the SERR , where is an area of a new surface created by a crack extension , 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]:where is the plane strain/stress thickness and is the length of the 2D element in the mesh at the crack front. This value coincides with the crack extension . 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).
(a)
(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 was met. The total SERR , where in the 2D case, is often considered for this purpose. In the forthcoming calculations, we will use the power mixity law as follows [10]:Thereby, the finite element model should be provided by the test data such as fracture toughness modes , , and , and the fitting coefficients , , and .
If contact between crack lips occurs, appropriate contact conditions along the interface of their interactions should be defined. Then the impenetrability and friction constraints can be imposed by establishing the Karush-Kuhn-Tucker conditions on in terms of displacements and tractions in the following form [22]: In the equations above we denote as in [22, 23] that is a gap function defining minimal distance between a point of the slave surface and its orthogonal projection on the master surface with a unit normal . Similarly, is a slip function describing a relative movement of those points on the contacting lines in the tangential direction in Figure 3. Also and are normal and shear scalar parameters of the contact traction vector , whereas 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 with a given coefficient of friction .
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:where the test functions and reside in appropriate vector spaces of kinematically admissible displacements and the admissible temperature field such that
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 . 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 , whence, the variational equality (15) can be reorganized into a variational inequality as in [22]: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 is discretized by a number of nonoverlapping elements such that the union of the elements comprises the total domain; that is, . 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 and as follows:where are nodal points and is a number of nodes in the basic element, and are shape functions associated with nodes and , and and are nodal unknown displacements and nodal unknown temperature field, respectively. The summation over the repeated index and 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: where , , , and are global mass, stiffness, conductivity, and capacity matrices, and are coupled matrices, and are global vectors of unknown displacements and temperature and their appropriate time derivatives, and , , and are global vectors of external mechanical, thermal, and contact forces. Herein an explicit form of the vector of nodal contact forces 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 , 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 to an instant , we come to a nonlinear system of algebraic equations expressed as [22] where , , and are partitioned matrices conditionally referring mass, damping, and stiffness matrices, respectively, combining respective matrices of mechanical and thermal parts; is a combined vector of global nodal displacements and temperatures and the dots over this vector denote its appropriate time derivatives; is a global thermal-mechanical force vector and is a global vector of contact forces.
The Newton-Raphson method reduces (21) to a linearized form, which at each time step has the following form:In (22), are submatrices of the coupled Jacobian matrix, and are corrections referring to the solution vectors and , accordingly, and and are components of the residual force vector which is defined from (21) as follows [9]:where the superscript “” represents a discrete time instant and the subscript “” 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 -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.05 in the calculations. The TSIFs are computed by using the -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.01 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 of 1 and length of 2 under quasi-static 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 -coordinate: where the subscripts “” and “” refer to ceramic and metal components, respectively.
It is supposed that the plate is under a stress-free state at the reference temperature and is subjected to the surface temperatures and at and , respectively. Because of the symmetry, a half plate is modelled only. The crack of length 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.05 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.
(a)
(b)
Figure 7(a) shows comparisons of the TSIFs normalized by 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 , , and . The plate undergoes a uniform dimensionless temperature , which is less than and corresponds to a variety of values such as 0.5, 0.2, 0.1, and 0.005.
(a)
(b)
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 , , and , which is undergoing a uniform dimensionless temperature , where and the ratio 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 of 10 mm and length 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 ZrO2 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 -coordinate, as shown in Figure 1, varies in accordance with a power function:From (25), it follows that the FGM is rich in ceramic when the parameter and rich in metal when the parameter . The effective properties of the metal/ceramic Ti-6Al-4V/ZrO2 plate with negligible porosity have been estimated by means of the rule of mixtures of a two-phase material as follows:
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 ; then the plate is heated on the ceramic surface by a high temperature , while the metallic surface remains at a low temperature equal to . 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 . 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 K, K, and K. The film convection coefficient is W/mm2 K and the temperature of a surrounding medium is 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 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 N/mm2 and 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 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.
(a)
(b)
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.
(a)
(b)
(c)
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 at different nondimensional moment of time 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.
(a)
(b)
(c)
(d)
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.
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.