Thermomechanically coupled sensitivity analysis and design optimization of functionally graded materials

Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911 www.elsevier.com/locate/cma

Thermomechanically coupled sensitivity analysis and design optimization of functionally graded materials Biaosong Chen, Liyong Tong

*

School of Aerospace, Mechanical and Mechatronic Engineering, University of Sydney, Bldg. J11, NSW 2006, Australia Received 9 September 2003

Abstract This paper presents a systematic numerical technique for performing sensitivity analysis of coupled thermomechanical problem of functionally graded materials (FGMs). General formulations are presented based on finite element model by using the direct method and the adjoint method. In the modeling of spatial variances of material properties, the graded finite element method is employed to conduct the heat transfer analysis and structural analysis and their sensitivity analysis. The design variables are the volume fractions of FGMs constituents and structural shape parameters. The design optimization model is then constructed and solved by the sequential linear programming (SLP). Numerical examples are presented to demonstrate the accuracy and the applicability of the present method.
2004 Elsevier B.V. All rights reserved.

1. Introduction In recent years, functionally graded materials (FGMs) have been shown to possess superior advantages when employed in high temperature environment. Tanigawa [35], Koizumi [15], Suresh and Mortensen [31], Noda [27] reviewed the development history of FGMs from the aspects of heat conduction analysis, thermal stress related problems, research activities and manufacturing techniques. One of the most important characteristics of FGMs is that they are designable. One can design the structure of FGMs both to resist the severe thermal loading due to high temperature and to maintain the structural rigidity and stiffness in a desired way by tailoring design parameters.

*

Corresponding author. Tel.: +61 2 93516949; fax: +61 2 93514841. E-mail address: [email protected] (L. Tong).

0045-7825/$ - see front matter
2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.07.005

1892

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

In respect of sensitivity analysis and design optimization, the research for general themoelastic problems has placed a solid foundation for the investigation of FGMs where both the heat conduction and the structural analysis can be taken into account at the same time. The work of sensitivity analysis started from the middle of 1980s. Meric [23,24] deduced the sensitivity equations for static problem in which the design variables were the material properties and the parameters of external loads. Dems and Mroz [8], Meric [25] investigated the sensitivity analysis for static and quasi-static problems with respect to size and shape design variables. Tortorelli et al. [36] presented the sensitivity analysis with adjoint method for the coupled constitutive model. Hou et al. [12], discussed the shape sensitivity analysis with both the direct and the adjoint method for static problem in conjunction with the finite element method. The similar work for 3D problem was conducted by Yang [39]. Most of the contributions mentioned above are based on the continuum model where the sensitivity equations were deduced directly from the partial differential equations. Another approach is to derive the sensitivity equations based on discrete model, in which the differential equations of a problem are discretized firstly in space by numerical techniques, such as the finite element method and then the sensitivity analysis is performed based on the discretized equations. The formulations on the discrete model are concise and facilitate the program implementation. Based on the finite element model, Chen et al. [3] and Liu [20] presented the direct method and the adjoint method, respectively, for both static and quasi-static problems. On the other hand, different kinds of techniques have been applied to design optimization for general thermoelastic problem. The design optimization work based on sensitivity analysis were given by Meric [25], Hou et al. [12], Chen et al. [3] and Liu [20] for general thermoelastic problems, Michaleris et al. [26] for weldment problem and Autio [1] for laminated plates. By using boundary integral equation, Lee and Kwak [17,18] investigated the shape design optimization for two-dimensional and axisymmetric problems. Kok et al. [16] employed response surface method in their research. Li et al. [19] performed the topology design optimization for thermoelastic problem by using evolutionary structural optimization (ESO). Xu and Grandhi [38] studied the multi-points approximation techniques in the design optimization of thermostructural problems. Due to the benefits of the application of FGMs, attentions have been focused on the design optimization of FGMs (see Tanigawa [35], Markworth et al. [21], Markworth and Saunders [22] and Noda [27]). Tanaka et al. [32–34] fulfilled a series of work on design optimization of volume fractions of FGMs by using the incremental finite element analysis and the direct sensitivity method. Cho and Ha [6,7] performed the design optimization of 1D and 2D volume fraction of FGMs. They conducted the sensitivity analysis by using the global finite difference method (FDM) and then employed the interior penalty function method to solve the optimal model. Huang et al. [13] studied design optimization through employing the strength and the energy storage capability of the FGMs structure as two objectives. Trueltaub [37] discussed the optimal control and optimization of FGMs in which the objective function can reflect the control problem or the optimization problem or both by tailoring its parameters. In his study, the sensitivity analysis was performed based on the continuum model and optimal criterion method was employed. In addition to the above work, some researchers have presented different schemes on the subject. Markworth and Saunders [22] proposed a simple optimization model for one-dimensional problem. Based on the analytical solutions, Ootao et al. investigated the application of some bionic algorithms to the design optimization of FGMs, such as the neural network method for the hollow sphere [28] and the hollow circular cylinder [29], and the genetic algorithm for the laminated composite plate [30]. In this paper, we apply the techniques of structural optimization for thermoelasticity [3,20] to the design optimization for FGMs which treat the volume fraction and structural shape design variables. Firstly, based on finite element model, sensitivity analysis schemes, in conjunction with the direct method and the adjoint method, are presented for both static and quasi-static problems. And then the graded finite element method proposed by Kim and Paulino [14] is employed to conduct the heat transfer and structural analysis. The advantage of the graded finite element method is its compatibility with the prob-

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

1893

lem of spatial varying material properties of FGMs. Detailed formulations are developed for the sensitivity analysis. In the program implementation section, the semi-analytical sensitivity analysis method are discussed and the sequential linear programming (SLP) method is employed to solve the optimal model. The software platform used is JIFEX [10] which is a multi-purpose system of finite element analysis and design optimization. Numerical examples are presented to demonstrate the accuracy and the applicability of the present method.

2. Sensitivity analysis of static problem The static problems are referred to as when both the heat conduction and structural analysis are static problems. The finite element equation of heat conduction can be expressed as KT ¼ R;

ð1Þ

where K, T, R are the heat conductance matrix, the nodal temperature vector and the heat load vector respectively. The finite element equation of structural problem takes the following form: K mu ¼ f ;

ð2Þ

where Km, u, f are the structural stiffness matrix, the nodal displacement vector and the external load vector. Suppose that in the optimal model, we have a index function which can be viewed as the constraint function or the objective function and has following the expression g ¼ gðx; T; uÞ;

ð3Þ

where x is the design variable (here only one design variable is concerned, however the following formulations are valid for the problem with multi-design variables), then its sensitivity with respect to the design variable x is dg og og dT og du ¼ þ þ : dx ox oT dx ou dx

ð4Þ

There are two methods to compute the sensitivity dg/dx: the direct method and the adjoint method. (1) Direct method: By differentiating Eqs. (1) and (2), we obtain K

dT dR dK ¼  T; dx dx dx

Km

du of of dT dK m ¼ þ  u: dx ox oT dx dx

ð5Þ ð6Þ

By substituting the solution to the above equations to Eq. (4), we can obtain the sensitivity of the index function. (2) Adjoint method: Firstly, we introduce the adjoint vector of heat conduction Kh and the adjoint vector of structural problem Km. By multiplying ðKTh Þ to Eq. (5) and ðKTm Þ to Eq. (6), and then combining these two equations with Eq. (4), we obtain

 dg og og dT og du dT dR dK du of of dT dK m ¼ þ þ  KTh K  þ T  KTm K m   þ u : ð7Þ dx ox oT dx ou dx dx dx dx dx ox oT dx dx

1894

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

The above equation can be re-arranged as

 dg og dR dK of dK m ¼ þ KTh  T þ KTm  u dx ox dx dx ox dx

 og df dT og du þ  KTh K þ KTm þ  KTm K m : oT dT dx ou dx

ð8Þ

To avoid solving for dT/dx and du/dx, let their coefficient terms be zero, i.e., thus we obtain the adjoint equation of the structural problem
T og m K Km ¼ ð9Þ ou and the adjoint equation of the heat conduction
T
T og df KKh ¼ þ Km : oT dT

ð10Þ

The solution to the adjoint equation of structural problem should be obtained prior to that of Eq. (10). The solution sequence of adjoint equations here is different from that in system analysis where the solution to heat transfer is obtained prior to that of structural problem. The sensitivity of index function can be simplified as

 dg og dK dK m T dR T of ¼ þ Kh  T þ Km  u : ð11Þ dx ox dx dx ox dx The difference between the direct method and the adjoint method lies in the computational expense. For a problem with multiple design variables and multiple index functions, in order to get the sensitivity results, the direct method needs to solve the sensitivity equations of Eqs. (5) and (6) for each design variable, whereas the adjoint method requires to solve the adjoint equations of Eqs. (9) and (10) for each index function. Therefore if the number of design variables is less than that of the index functions, the direct method is computationally cost effective, otherwise is the adjoint method. Of course, prior to performing the sensitivity analysis, Eqs. (1) and (2) should be solved firstly and then the main task is to compute dK/dx, dR/dx, dKm/dx.

3. Sensitivity analysis of quasi-static problem Here the quasi-static problem means that the heat conduction concerned is now a transient problem M T_ þ KT ¼ R;

ð12Þ

where M is the heat capacity matrix. The structural problem is still a static one as defined in Eq. (2). It implies that the structural responses caused by the time dependent temperature vary so slow that the inertial terms can be neglected in the structural problem. The index function has the expression Z tf gðx; T; uÞ ¼ pðx; T; u; tÞ dt ð13Þ 0

and its sensitivity with respect to design variable x is given by  Z tf
dg op op dT op du ¼ þ þ dt: dx ox oT dx ou dx 0

ð14Þ

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

1895

The direct method and the adjoint method can be employed as follows: (1) Direct method: By differentiating Eq. (13), we obtain M

dT_ dT dR dM _ dK þK ¼  T: T dx dx dx dx dx

ð15Þ

The solutions to Eqs. (15) and (6) can be substituted into Eq. (14) to obtain the sensitivity of the index function. (2) Adjoint method: Following the same procedures in Section 2, the adjoint vector of heat conduction Kh and the adjoint vector of structural problem Km are introduced firstly. By multiplying ðKTh Þ to Eq. (15) and ðKTm Þ to Eq. (6), and then combining them with Eq. (14), we obtain !  Z tf
Z tf _ dg op op dT op du d T dT dR dM dK ¼ þ þ þK  þ T dt KTh M dt  T_ þ dx ox oT dx ou dx dx dx dx dx dx 0 0
 Z tf of of dT dK m T m du   þ u dt: ð16Þ  Km K dx ox oT dx dx 0 Prior to the simplifying the above expression, it is noted that t Z tf Z tf   T dT_ dT  f T T _ dT dt: dt ¼ Kh M   ð17Þ Kh M K_ h M þ KTh M d dx 0 dx 0 0   ¼ 0 and M _ ¼ 0 for linear problem, we Let the adjoint vector satisfy Kh(tf) = 0, and noting that dT dx t¼0 obtain Z tf Z tf dT_ T dT T dt ¼  dt: ð18Þ Kh M K_ h M dx dx 0 0 By substituting Eq. (18) into Eq. (16) and through some algebraic manipulations, Eq. (16) can be expressed as

 Z tf
dg op dM _ dK dK m T dR T of ¼ þ Kh  T þ Km  u dt T dx ox dx dx dx ox dx 0   Z tf
Z tf
T op op dT T m du T T of _  Km K dt þ þ Kh M  Kh K þ Km dt: ð19Þ þ ou dx oT oT dx 0 0 To avoid solving for dT/dx and du/dx, let their coefficient terms be zero and thus we have the adjoint equation of the structural problem
T op K m Km ¼ ð20Þ ou and the adjoint equation for the heat conduction
T
T op of _ M Kh þ KKh ¼ þ Km oT oT

ð21Þ

with the boundary condition at final time Kh ðtf Þ ¼ 0: Then Eq. (21) can be solved by using the backward time marching method.

ð22Þ

1896

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

The sensitivity of index function in Eq. (19) can be simplified as

 Z tf
dg op dR dM _ dK of dK m ¼ þ KTh  T þ KTm  u dt: T dx ox dx dx dx ox dx 0

ð23Þ

The solution sequence of adjoint method here is Eq. (20) first and then Eq. (21). The difference between the direct method and the adjoint method in computational expense for quasi-static problem is the same as that of static problem. It is noted that the sensitivity equations and the adjoint equations for static and quasi-static problems are identical to their original equations except for the vector on the right hand side. Therefore same solution method can be applied. The precise time integration methods in Ref. [11] are employed to perform the sensitivity analysis of transient heat conduction.

4. Graded finite element method and sensitivity analysis of FGMs The primary characteristic of FGMs is its spatial variance of material properties. Therefore, the generalized isoparametric formulation [14] method is employed to construct the finite element model. The graded finite element uses the isoparametric shape functions, which are the same for spatial coordinates, to interpolate the element material properties with nodal values. That is P¼

ne X

N iP i;

ð24Þ

i¼1

where the property P can be the conductivity k, or capacity qc, or the YoungÕs modulus E, or PoissonÕs ratio m, or thermal expansion coefficient a. The graded finite element method indicates that the material properties should be assigned to each node and then the corresponding elemental matrix and vector can be calculated using Gaussian quadrature. The above isoparametric formulations are valid for any isoparametric elements, such as membrane, axisymmetric and three-dimensional problem. In the process of sensitivity analysis, the main work for both the direct method and the adjoint method is to compute the derivatives of dK/dx, dM/dx, dR/dx for heat conduction and dKm/dx, of/ox, of/oT for structural problem. As for the design variable of volume fraction, we can derive their analytical expressions for each element and then assembly them for the global solution. In terms of the element heat conductance, we have K e ¼ K e1 þ K e2 ; ð25Þ where K e1 is related to the element conductivity and K e2 is the contribution of surface convection. Here, just K e1 depends on design variable and its sensitivity is
Z  Z d d dk ðK e1 Þ ¼ ð26Þ B T kBdX ¼ B T BdX; dx dx Xe dx Xe where B is the derivatives of element shape function with respect to spatial coordinates. As for the element heat capacity matrix, terms of the lumped mass matrix are Z M eii ¼ qc dX=ne ; ð27Þ Xe

where ne is the number of element nodes. Eq. (27) means that the total element heat capacity is assigned to each node evenly. Then its sensitivity expression is given by

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

d d ðM eii Þ ¼ dx dx


Z

qc dX=ne

 ¼

Xe

Z
Xe

 dðqcÞ dX=ne : dx

1897

ð28Þ

The element heat load can be expressed as Re ¼ ReQ þ ReC  ReB  Req ;

ð29Þ

where ReQ , ReC , ReB and Req are the element heat load vectors related to the internal heat generation, surface convection, specific nodal temperature and the surface heat flux. Here only the term ReB depends on the design variables and thus its derivative takes the form of d d dK eB e ðReB Þ ¼ ðK eB T eB Þ ¼ T : dx dx dx B In the structural problem, the element elastic matrix is Z e BT DB dX; ðK m Þ ¼

ð30Þ

ð31Þ

X

where D is the elastic matrix. The sensitivity of the element elastic matrix is Z dðK m Þe dD B dX: ¼ BT dx dx X

ð32Þ

In order to clarify the calculation of of/ox, we firstly express the external load vector f as f ¼ f th þ f m ;

ð33Þ

where f th is the thermal load and f m is independent of temperature change. Here f m may be the contribution of specific displacement, external mechanical load, gravity load or centrifugal load etc. If f m depends on the design variable, its sensitivity must be considered. Here just f th is considered Z e ðf th Þ ¼ B T DaDT dX; ð34Þ X

then its partial derivative with respect to design variable is
 Z e oðf th Þ oa T oD aþD ¼ B DT dX: ox ox ox X The partial derivatives with respect to the nodal temperature are Z e oðf th Þ oðDT Þ dX; ¼ B T Da oT oT X

ð35Þ

ð36Þ

where DT can be expressed as DT ¼

ne X

N i ðT i  T r Þ;

ð37Þ

i¼1

where Tr is the nodal reference temperature for zero thermal stress. The calculation of Eq. (36) can be conducted one by one for each nodal temperature. In the research of FGMs, there are two methods to describe the variance of the material properties. The first method is to use the specific functions for all kinds of material properties. In general, this method is proposed for the theoretical analysis. For example, one kind of the material properties and its sensitivity can be expressed as

1898

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

dP dP ðx; rÞ ¼ ; dx dx

P ¼ P ðx; rÞ;

ð38Þ

where x is design variable and r denotes the spatial coordinates. The expressions for all the specific functions of material properties can be derived analytically. The second method is to employ the specific functions of volume fraction of FGMs. In conjunction with the estimation methods of material properties, we can obtain the mixture properties according to volume fraction. The second method is of much more importance for engineering application. The specific function of volume fractions and its sensitivity can be expressed as dV dV ðx; rÞ ¼ ; dx dx

V ¼ V ðx; rÞ;

V 2 ½0; 1;

ð39Þ

where V is the volume fraction. If the function form is defined, its analytical sensitivity can be obtained. FGM usually consists of two kinds of materials. One is metal and the other is ceramic. The volume fraction relationship is V c þ V m ¼ 1;

dV c dV m ¼ ; dx d

V c ; V m 2 ½0; 1;

ð40Þ

where the subscript ‘‘c’’ and ‘‘m’’ denote the ceramic and metal, respectively. There are several models to predict the mixture properties of FGMs. In this paper the Voigt-type model and the Wakashima–Tsukamoto (W–K) model (provided in Ref. [5]) are considered. The materials are assumed to be macroscopically isotropic. The formulations and their derivatives for Voigt-type model and W–K model are provided in Appendix A. After obtaining the sensitivity of all the mixture properties, the next step is to compute sensitivity of the elastic matrix. The 2D plane problem (plane stress and plane strain) and axisymmetric problem are considered in this paper and their derivatives are detailed in Appendix B. And then all the sensitivity of element matrix and vector can be carried out by using the Gaussian quadrature. It is noted that the above equations are algebraic manipulations only. If other kind of problem (such as 3D problem) is considered, we can follow the above steps to perform the sensitivity analysis. As for the shape design variables, the derivatives of element matrix have the expressions d d ðK e Þ ¼ dx 1 dx dðK m Þe ¼ dx


Z

T



B kB dX Xe

 Z
T dB T dk T dB BþB k kB þ B ¼ dX; d dx dx Xe

 Z
T dB dB T dD T BþB D DB þ B dX: dx dx dx X

ð41Þ

ð42Þ

Then the derivatives of material properties are expressed as  ne
dk X dN i dki ¼ ki þ N i ; dx dx dx i¼1

ð43Þ

 ne
dðqcÞ X dN i dðqcÞi ¼ ðqcÞi þ N i ; dx dx dx i¼1

ð44Þ

 ne
dD X dN i dDi ¼ Di þ N i : dx dx dx i¼1

ð45Þ

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

1899

Other kinds of derivative expressions can be deduced by using the same procedures. The above expression results in the derivatives of B and N with respect to the design variable. In practice, the shape design variables will have many forms and the function relationships with B and N are implicit, hence it is difficult to derive the explicit expressions of their derivatives. Here, the semi-analytical technique [2,9] is employed. The formulations for the heat conduction are dM e M e ðx þ dxÞ  M e ðxÞ ; Dx dx

ð46Þ

dK e K e ðx þ dxÞ  K e ðxÞ : Dx dx

ð47Þ

Because the heat conduction is a scalar field problem, the above expressions can possess high efficiency and precision. As for the structural sensitivity analysis, it has been discovered that the semi-analytical method could cause significant inaccuracy in some cases. Detailed discussion on the subject can be found in Refs. [2,9] and references therein. In this paper, the alternative forward–backward finite difference method [4] is used to improve the numerical accuracy. The method is implemented by using the forward difference scheme for one element and the backward difference scheme for its neighboring element and so forth. Because of the employment of mapping method in this investigation to generate the finite element mesh, it is easy to control the ordering of element number and then the one easy way to identify the neighboring relationship in program implementation is to use the even number and odd number of the elements. In fact, the software platform JIFEX [10] was developed by the authors of Ref. [4] and the even–odd number method was implemented therein. That is in the structural finite element model, the derivative of element with even number is expressed as dðK m Þe K m ðx þ DxÞ  K m ðxÞ Dx dx

ð48Þ

and that of element with odd number is dðK m Þe K m ðxÞ  K m ðx  DxÞ : Dx dx

ð49Þ

Other kinds of derivatives of element matrix or vector with respect to shape design variables can be deduced with the same approach. After assembling the element derivatives, the sensitivity equations of the adjoint equations can be solved. The semi-analytical method is adaptive to general problems. And because it is used just at the element level, it is computationally cost effective.

5. Program implementation Both the direct method and the adjoint method can be used to compute the sensitivity of the problem by following the above equations. However in the program implementation, some practical problems should be considered. One of those is the software complexity of modules integration. Here the heat conduction module and the structural module should be integrated and then the amount of information exchange between two modules is the key problem. In the direct method, we just need to transfer the results of temperature and their sensitivity from heat conduction module to the structural module. The direction of information exchange is one way. However, in the adjoint method, the direction of information exchange is two way, namely we not only need to transfer the results of temperature from heat conduction module to the structural module in the analysis phase, but also need to transfer the results of adjoint vector Km from the structural module to the heat conduction module in the sensitivity analysis phase. The program

1900

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

complexity of the adjoint method is greater than that of the direct method. Therefore in this research, the direct method is employed. Another reason to choose the direct method is that some approximate techniques can be applied. In Eqs. (35) and (36), it can be seen that the computational cost of the sensitivity of thermal load will be great, especially for Eq. (36), the partial derivatives with respect to nodal temperature. In addition, if other kinds of loads (such as gravity load, centrifugal load) which depend on design variables are considered, their sensitivity must be taken into account and then the job will become more complex. In the program implementation, the matured and reliable method to manipulate the shape design variables is to use the mapping method [9], which is a structured method of mesh generation, to manipulate the shape design variables. The control sequence is: shape design variables ! design nodes and design elements ! finite element model. The requirement is that a design element model should be provided in advance. In the direct method, based on the temperature sensitivity, the finite difference method at the element level can be applied to compute the derivatives of thermal load in structural analysis, i.e., df e f e ðx þ Dx; T þ DTÞ  f e ðx; TÞ : Dx dx

ð50Þ

In order to implement the semi-analytical method, our job is to perturb the finite element model for each design variable and then to carry out the sensitivity analysis of heat conduction firstly. Based on the results of temperature sensitivity, the additional work is to perturb the nodal temperature and to compute the element stiffness matrix and thermal load vector. After obtaining the sensitivity of structural displacement, the sensitivity of elemental stress can be calculated as dre re ðx þ Dx; T þ DTÞ  re ðx; TÞ : Dx dx

ð51Þ

Based on the results of analysis and sensitivity analysis, the optimization procedures are outlined below: Step Step Step Step Step Step

1. 2. 3. 4. 5. 6.

Analysis of heat conduction; Analysis of structural problem induced by thermal load; Sensitivity analysis of heat conduction; Sensitivity analysis of structural problem; Solution of the optimal model by sequential linear programming (SLP); If converge, then stop, otherwise return to step 1.

The present methods have been implemented in the platform of JIFEX [10].

6. Numerical examples Example 1. A long slab L · H = 1.0 m · 0.1 m of FGMs is considered. Two ends of the slab are fixed and the upper surface is subjected to a heat flux of q = 1.0 · 105 W/m2 while the lower surface is subjected to specific temperature Tb = 298.0 K. The reference temperature for zero thermal stress is Ti = 273.0 K. The computational model is illustrated in Fig. 1 by considering the symmetry of the problem. The YoungÕs modulus and the thermal conductivity are defined by specific functions E ¼ E0 eky

ðE0 ¼ 10:0 GPa; k ¼ 20:0Þ;

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

1901

y Upper Surface: Heat Flux

H

x Lower Surface: Specific Temperature L/2

Fig. 1. Model of the slab.

k ¼ k 0 ejy

ðk 0 ¼ 100:0 W=m K; j ¼ 10:0Þ:

Other material properties are constant: PoissonÕs ratio m = 0.3 and the coefficient of thermal expansion a = 1.0 · 104. The design variables are k and j. Both plane stress model and plane strain model are considered. The sensitivity results of nodal temperature, displacement and stress are listed in Tables 1–4 where FDM means global finite difference method and SA means sensitivity analysis. The temperature is independent of k, then the corresponding results of sensitivity are zero and are not given in Tables 2 and 4. Firstly the above tables show that for the design variables of material parameters the FDM results are not so sensitive to the step size. And the results by the sensitivity analysis correlate well with those by FDM, which validates the correctness of the methods proposed in this investigation. An important information revealed here is that even though the elastic property is not directly related to j, the results of sensitivity in Tables 1 and 3 exhibit the effects of temperature on the structural response, especially for the structural stresses. Such strong coupling effects emphasize the necessity to carry out the entire sensitivity analysis for both the heat conduction and the structural analyses when performing the design optimization for thermomechanically structural system. Example 2. The sensitivity analysis of a circular disk with respect to shape design variables is presented. The axisymmetry finite element model is employed. The configuration and the finite element model are shown in Fig. 2. The geometry parameters are disk thickness = 0.02 m, r1 = 0.2 m, r2 = 0.4 m. Ni and Al2O3 are selected as the constituents of FGM. Their material properties are summarized in Table 5. Ni acting as the metal layer is placed on the outer surface while Al2O3 acting as the ceramic layer is placed on the inner surface. The inner surface is subjected to a heat flux of q = 1.0 · 105 W/m2 infusing into the disk. The outer surface is subjected to a convection which has the b = 5 · 104 W/(m2 K) and the ambient temperature T1 = 298.0 K. The initial temperature for all nodes is Ti = 298.0 K. Vc denotes the ceramic

Table 1 Sensitivity results with respect to j for plane stress problem Location x = 5.0

y = 0.025

y = 0.05

y = 0.075

y = 0.1

T

FDM Dk = 0.02 FDM Dk = 0.10 SA

0.30768 0.30772 0.30767

1.60463 1.60303 1.60503

4.42742 4.41923 4.42947

9.53743 9.51208 9.54378

uy

FDM Dk = 0.02 FDM Dk = 0.10 SA

2.30067E05 2.29557E05 2.30194E05

1.96986E05 1.96504E05 1.97107E05

9.47499E06 9.44371E06 9.48284E06

1.39614E05 1.39357E05 1.39678E05

rxx

FDM Dk = 0.02 FDM Dk = 0.10 SA

3.53463E+05 3.53784E+05 3.53382E+05

3.95938E+06 3.95637E+06 3.96013E+06

1.87745E+07 1.87431E+07 1.87824E+07

7.86578E+07 7.84437E+07 7.87114E+07

1902

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

Table 2 Sensitivity results with respect to k for plane stress problem Location x = 5.0

y = 0.025

y = 0.05

y = 0.075

y = 0.1

uy

FDM Dk = 0.02 FDM Dk = 0.10 SA

3.17052E06 3.16995E06 3.23705E06

2.93468E06 2.93417E06 2.99452E06

2.63850E06 2.63806E06 2.67286E06

2.26347E06 2.26312E06 2.22829E06

rxx

FDM Dk = 0.02 FDM Dk = 0.10 SA

2.16428E+06 2.16643E+06 2.16383E+06

1.20407E+07 1.20647E+07 1.20383E+07

4.53082E+07 4.54439E+07 4.52981E+07

1.54317E+08 1.54935E+08 1.54328E+08

Table 3 Sensitivity results with respect to j for plane strain problem Location x = 5.0

y = 0.025

y = 0.05

y = 0.075

y = 0.1

T

FDM Dk = 0.02 FDM Dk = 0.10 SA

0.30768 0.30772 0.30767

1.60463 1.60303 1.60503

4.42742 4.41923 4.42947

9.53743 9.51208 9.54378

uy

FDM Dk = 0.02 FDM Dk = 0.10 SA

3.20181E05 3.19442E05 3.20366E05

2.83504E05 2.82797E05 2.83681E05

1.70592E05 1.70072E05 1.70723E05

8.77550E06 8.76481E06 8.77818E06

rxx

FDM Dk = 0.02 FDM Dk = 0.10 SA

2.50969E+05 2.51553E+05 2.50824E+05

3.99612E+06 3.99397E+06 3.99665E+06

1.97040E+07 1.96738E+07 1.97115E+07

9.31105E+07 9.28535E+07 9.31749E+07

Table 4 Sensitivity results with respect to k for plane strain problem Location x = 5.0

y = 0.025

y = 0.05

y = 0.075

y = 0.1

uy

FDM Dk = 0.02 FDM Dk = 0.10 SA

2.49877E06 2.49809E06 2.56247E06

2.23306E06 2.23244E06 2.29180E06

1.90581E06 1.90527E06 1.94577E06

1.49665E06 1.49620E06 1.48552E06

rxx

FDM Dk = 0.02 FDM Dk = 0.10 SA

2.34059E+06 2.34289E+06 2.33964E+06

1.31149E+07 1.31409E+07 1.31080E+07

4.94542E+07 4.96018E+07 4.94194E+07

1.77672E+08 1.78382E+08 1.77588E+08

r2 12 13 14 15 16 17 18 19 20 21 22 Heat flux

Heat convection 1

r1 axis

2

3

4

5

6

7

8

9

10

11

z r Fig. 2. Axisymmetry model of circular disk.

volume fraction which varies linearly along the radius. The nodal volume fractions are listed in Table 6. The precise time integration method [11] is employed to solve the transient heat conduction and to perform its

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

1903

Table 5 Thermal properties of Ni and Al2O3 Properties

Constituents

Density (kg/m3) Specific heat (J/kg K) Thermal conductivity (W/m K) YoungÕs modulus (Gpa) PoissonÕs ratio Thermal expansion coefficient (K · 106)

Ni (metal)

Al2O3 (ceramic)

8900.0 444.0 90.7 199.5 0.3 13.3

3970.0 775.0 30.7 393.0 0.25 8.8

Table 6 Nodal volume fractions Node

1, 12

2, 13

3, 14

4,15

5, 16

6, 17

7, 18

8, 19

9, 20

10, 21

11, 22

Vc

1.0

0.9

0.8

0.7

0.6

0.5

0.4

0.3

0.2

0.1

0.0

sensitivity analysis. The time step length is set to be Dt = 10 s. The design variables are the inner radius r1 and outer radius r2. The purpose of using the nodal volume fractions instead of a specific analytical function for the variances of volume fraction along the radius is to make the design variables and the volume fraction independent of each other. Both the Voigt-type model and the Wakashima–Tsukamoto (W–K) model are employed to predict the properties of FGMs. The sensitivity of thermal displacement and stress at t = 15 min and t = 30 min are considered for three nodes. The sensitivities are listed in Table 7 (for r1) and Table 8 (for r2). It can be seen that all the sensitivity results (SA) correlate well with those obtained by the global finite difference method (FDM), which confirms the correctness of the methods in this paper for the shape design variables. It is worth pointing out that the estimation models have some effects on the sensitivity results and the results by Voigt and W–K models for the same node at the same time are different. Some extreme cases are those in Table 7, for example the results for node 12 at 30 min, where the signs of the number are opposite. This means that the estimation models should have some effects on the design optimization via the sensitivity analysis. Another valuable information is that the sensitivity results will be different for different time and then such time history should be taken into account for time dependent problem. Table 7 Sensitivity results with respect to r1 Node no.

ur

rhh

12

18

22

12

18

22

4.421E04 4.412E04

2.484E04 2.479E04

2.030E04 2.025E04

9.313E+08 9.316E+08

8.667E+08 8.656E+08

8.784E+08 8.784E+08

SA FDM

2.511E04 2.506E04

1.028E05 9.820E06

9.758E05 9.703E05

7.707E+08 7.703E+08

7.030E+08 7.034E+08

5.192E+08 5.193E+08

t = 30 min Voit SA FDM

8.556E05 8.475E05

5.575E04 5.586E04

7.651E04 7.667E04

1.452E+09 1.453E+09

1.198E+09 1.198E+09

3.085E+08 3.076E+08

1.136E04 1.143E04

7.740E04 7.753E04

9.794E04 9.813E04

1.211E+09 1.212E+09

8.240E+08 8.247E+08

2.513E+07 2.513E+07

t = 15 min Voit SA FDM W–K

W–K

SA FDM

1904

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

Table 8 Sensitivity results with respect to r2 Node no.

t = 15 min Voit W–K t = 30 min Voit W–K

ur

rhh

12

18

22

12

18

22

SA FDM

1.376E03 1.384E03

1.132E03 1.134E03

1.394E03 1.396E03

1.724E+09 1.723E+09

1.218E+09 1.216E+09

1.256E+09 1.258E+09

SA FDM

1.092E03 1.097E03

7.480E04 7.491E04

9.748E04 9.759E04

1.514E+09 1.513E+09

9.822E+08 9.810E+08

8.255E+08 8.264E+08

SA FDM

1.396E03 1.402E03

5.721E04 5.728E04

8.966E04 8.973E04

2.320E+09 2.315E+09

1.834E+09 1.833E+09

8.704E+08 8.711E+08

SA FDM

9.726E04 9.762E04

1.531E04 1.532E04

4.046E04 4.047E04

2.011E+09 2.007E+09

1.384E+09 1.383E+09

3.956E+08 3.957E+08

Example 3. In order to illustrate the application scope of the present method, the design optimization of volume fraction of a turbine disk are studied. The geometrical configuration is shown in Fig. 3. Its design element model is shown in Fig. 4. The material constituents are Ni and Al2O3 whose properties are listed in Table 5. The outer surface of the disk is subjected to high specific temperature while the inner surface is subjected to low specific temperature. The upper and lower surfaces are subjected to convection condition. In addition to the thermal loads, the turbine disk is subjected to centrifugal load and distributed load parallel to the radius direction acting on the outer surface. The inner part of the disk (low temperature area) is made of Ni while the outer part (high temperature area) is made of Al2O3. In order to maintain the manufacturability, special geometrical parameters are selected as design variables (shown in Table 9). The function of volume fraction of ceramic is defined by an specific analytical function as

Fig. 3. Geometrical configuration of turbine disk (dimensions in millimeters).

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

1905

Fig. 4. Design element model of turbine disk.

Table 9 Optimal results of design variables for turbine disk Design variables

Initial

Upper bound

Lower bound

Voit optimal

W–K optimal

L2 L3 L4 L5 L6 H1 H2 H3 H4 R1 R2 R3 R4 B3 x

95.00 106.00 115.00 221.00 225.00 105.00 83.00 87.00 43.00 60.00 57.00 6.00 10.00 77.00 1.00

97.00 108.00 117.00 224.00 229.00 105.00 83.00 87.00 56.00 63.00 65.00 9.00 14.00 81.00 3.00

92.00 105.00 112.00 221.00 225.00 97.00 75.00 79.00 40.00 52.00 45.00 3.00 6.00 70.00 0.05

94.93 107.32 112.98 223.48 227.85 103.50 78.91 83.60 41.90 59.10 64.17 6.55 9.40 74.48 1.22

96.19 105.45 115.48 221.92 225.21 104.15 82.03 83.60 44.84 62.97 63.28 4.79 7.20 78.88 0.85

8 r < 115; > < 0; x V c ðrÞ ¼ ððr  115Þ=100Þ ; 115 6 r 6 215; > : 1; r > 215; where x is the design variable of volume fraction. The mathematic model of design optimization is formulated as Min s:t:

W maxðrceramic Þ 6 250 MPa; maxðrmetal Þ 6 300 Mpa; maxður Þ 6 1:2 mm;

where W, rceramic, rmetal, ur are the structural weigh, the VonMisses stress in the ceramic region, the VonMisses stress in the metal region and the radial displacement. The Voigt-type model and the

1906

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

Table 10 Optimal results of structural responses for turbine disk Voit

Initial Optimal

W–K

W (kg)

rmetal max

36.40 30.40

302.41 300.66

(MPa)

rceramic max 222.07 249.89

(MPa)

ur (mm)

W (kg)

rmetal max (MPa)

rceramic (MPa) max

ur (mm)

1.16 1.20

36.40 30.57

304.50 296.38

252.56 194.33

1.17 1.17

Wakashima–Tsukamoto (W–K) model are employed in the computation process. The results of design variables are listed in Table 9. The results of structural responses are listed in Table 10. The variances of finite element meshes are shown in Fig. 5. In Fig. 5(a), the finite element meshes of initial design are shown together with the element number in order to present the relationship of elements when performing the alternative forward–backward finite difference method. Firstly after design optimization, the structural responses arrive at the desirable stage according to the design optimization model where the objective, the structural weight, is reduced by over 10% while keeping the stress and displacement below the initial given values, respectively. The example demonstrate the correctness and applicability of the optimization procedures. Another important information worth noting is that for different estimation models, the optimal results will be different. However it is not easy to determine which estimation models is better under the

Fig. 5. (a) Finite element meshes of initial design, (b) finite element meshes of optimal design with Voit-type model, (c) finite element meshes of optimal design with Wakashima–Tsukamoto model.

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

1907

circumstances. Therefore the selection of estimation models is of great importance for the design optimization of the FGMs structures in practice.

7. Conclusions In this paper, the thermomechanically coupled sensitivity analysis of heat conduction and structural responses for functionally graded materials is presented based on finite element model. The formulations of both the direct method and the adjoint method are developed for static and quasi-static problems. In the modeling of FGMs, the graded finite element method is applied. As for the design variable of volume fraction, the derivatives of element matrix (conductance matrix, elastic matrix) are deduced analytically. As for the shape design variables, the semi-analytical method is employed to conduct the sensitivity analysis. The applications of Voigt and Wakashima–Tsukamoto estimation methods and their sensitivity analysis are discussed in detail with plane stress, plane strain and axisymmetry problems. Based on the sensitivity results, an optimal model is constructed and solved by using the sequential linear programming. The methods are adaptive to the design optimization of general FGMs structures. Numerical examples are presented to validate the accuracy and efficiency of the present methods of sensitivity analysis and design optimization. At the same time, numerical results reflect the influences of FGMs estimation methods on the optimal results.

Acknowledgement The authors are grateful to the support of Australian Research Council.

Appendix A In terms of Voigt-type model which is simply an arithmetic mean, the mixture property can be P ¼ P mV m þ P cV c;

ðA:1Þ

where the over-bar denotes the mixture property. P may be E, m, q, a, qc and k. Then the sensitivity is given by dP dV c ¼ ðP c  P m Þ : dx dx

ðA:2Þ

The derivative formulations for W–K model and As for the W–K model, the mixture properties of q and qc are defined by Eqs. (A.1) and (A.2), whereas other material properties can be calculated using more complex formulas. The thermal conductivity is k ¼ k m þ

k m V c ðk c  k m Þ : k m þ ðk m  k c ÞV m =3

ðA:3Þ

As for the mechanical properties, if we start from the properties of E and m, we should firstly compute the bulk modulus K and the shear modulus G K c;m ¼

Ec;m ; 3ð1  2mc;m Þ

Gc;m ¼

Ec;m : 2ð1 þ mc;m Þ

ðA:4Þ

1908

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

Then the mixture K and G are K ¼ Km þ

aV c K m ðK c  K m Þ ; V m K c þ aV c K m

G ¼ Gm þ

bV c Gm ðGc  Gm Þ ; V m Gc þ bV c Gm

ðA:5Þ

where a and b are a¼

K c ð3K m þ 4Gm Þ ; K m ð3K c þ 4Gc Þ

b¼

ð1 þ eÞGc ; Gm þ eGc

e¼

9K m þ 8Gm : 6K m þ 12Gm

ðA:6Þ

The mixture E and m are E¼

9K G ; 3K þ G

m ¼

3K  2G 2ð3K þ GÞ

ðA:7Þ

and the coefficient of thermal expansion  a is  a ¼ am 

ð1=K  1=K m Þðac  am Þ : 1=K m  1=K c

ðA:8Þ

Then the sensitivity of macroscopic mixture properties should have the following expressions where the derivatives with respect to design variable dP/dx is simplified as P 0 and so on k 0 ¼ k m V c ðk c  k m Þðk m þ ðk m  k c Þ=3Þ V 0 ; c ðk m þ ðk m  k c ÞV m =3Þ2 0

K ¼ 0

G ¼

aV 0c K m ðK c  K m ÞðV m K c þ aV c K m Þ  aV c K m ðK c  K m ÞðV 0m K c þ aV 0c K m Þ 2

ðV m K c þ aV c K m Þ

bV 0c Gm ðGc  Gm ÞðV m Gc þ bV c Gm Þ  bV c Gm ðGc  Gm ÞðV 0m Gc þ bV 0c Gm Þ ðV m Gc þ bV c Gm Þ 0

0

0

E ¼

ð3K þ GÞ2 0

2

0

2

ðA:10Þ

;

ðA:11Þ

ðA:12Þ

; 0

ð3K  2G Þð3K þ GÞ  ð3K  2GÞð3K þ G Þ 2ð3K þ GÞ

;

0

0

9ðK G þ K G Þð3K þ GÞ  9K Gð3K þ G Þ

0

m0 ¼

ðA:9Þ

;

ðA:13Þ

0

 a0 ¼

ðac  am Þ K : 1=K m  1=K c K 2

ðA:14Þ

Appendix B As for the plane stress problem, we have 2 3 1 m 0 E 6 0 7 6m 1 7; D¼ 4 1  m2 1 m5 0 0 2

ðB:1Þ

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

2 D0 ¼

1

E ð1  m Þ þ 2mm E 6 6m 4 ð1  m2 Þ2 0 0

2

0

0

m

3

2

m0

1

6 0 0 7 7 þ E 6m 5 2 1m 4 1m 0 2

1 0

1 0

0

1909

3

0 7 7: m0 5 2

ðB:2Þ

As for the plane strain problem, we have 2 3 1m m 0 6 m E 1m 0 7 6 7; D¼ ð1 þ mÞð1  2mÞ 4 1  2m 5 0 0 2

ðB:3Þ

2 D¼

E0 ð1 þ mÞð1  2mÞ  E½m0 ð1  2mÞ  2m0 ð1 þ mÞ 6 6 4 ½ð1 þ mÞð1  2mÞ2 2 þ

m0

E 6 0 4 m ð1 þ mÞð1  2mÞ 0

m0 m0 0

0

1m

m

m

1m

0

0

3

0

3

0 7 7 1  2m 5 2

ðB:4Þ

7 0 5: m0

As for the axisymmetry problem, we have 2 m m 1 1m 1m 6 6 m m 6 1 1  m 1  m Eð1  mÞð1 þ mÞ 6 6 D¼ 6 m m 6 ð1  2mÞ 1 61  m 1  m 6 4 0 0 0

0

3

7 7 7 0 7 7 7; 7 0 7 7 1  2m 5 2ð1  mÞ

ðB:5Þ

2

1

6 6 m 6 0 0 0 0 ½E ð1  mÞ  Em ð1 þ mÞð1  2mÞ þ Eð1  mÞ½m þ 4mm  6 61  m 0 D ¼

6 m 2 6 ½ð1 þ mÞð1  2mÞ 61  m 6 4 0 3 2 m0 m0 0 0 7 6 ð1  mÞ2 ð1  mÞ2 7 6 7 6 7 6 m0 m0 7 6 0 0 2 2 7 6 ð1  mÞ Eð1  mÞ 7 6 ð1  mÞ þ 7: 6 0 0 7 m ð1 þ mÞð1  2mÞ 6 m 7 6 0 0 7 6 ð1  mÞ2 ð1  mÞ2 7 6 7 6 m0 5 4 0 0 0 2 2ð1  mÞ

m m 1m 1m m 1 1m m 1 1m 0

0

0

3

7 7 7 0 7 7 7 7 0 7 7 1  2m 5 2ð1  mÞ

ðB:6Þ

1910

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

References [1] M. Autio, Coupled thermal-structural problems in the optimization of laminated plates, Struct. Optimization 15 (1998) 49–56. [2] H. Boer, F. Keulen, Error analysis of refined semi-analytical design sensitivities, Struct. Optimization 14 (1997) 242–247. [3] B.S. Chen, H.B. Zhao, Y.X. Gu, Sensitivity analysis and design optimization of the coupling thermo-structures, in: Proceedings of the Fourth World Congress on Structural and Multidisciplinary Optimization, 2001, pp. 406–413. [4] G.D. Cheng, Y.X. Gu, Y.Y. Zhou, Accuracy of semi-analytic sensitivity analysis, Finite Elements Anal. Des. 6 (1989) 113–128. [5] J.R. Cho, J.T. Oden, Functionally graded material: a parametric sturdy on thermal-stress characteristics using the Crank– Nicolson–Galerkin scheme, Comput. Methods Appl. Mech. Engrg. 188 (2000) 17–38. [6] J.R. Cho, D.Y. Ha, Volume fraction optimization for minimizing thermal stress in Ni–Al2O3 functionally graded materials, Mater. Sci. Engrg. A 334 (2002) 147–155. [7] J.R. Cho, D.Y. Ha, Optimal tailoring of 2D volume fraction distributions for heat resisting functionally graded material using FDM, Comput. Methods Appl. Mech. Engrg. 191 (2002) 3195–3211. [8] K. Dems, Z. Mroz, Variational approach to sensitivity analysis in thermoelasticity, J. Therm. Stress. 10 (1987) 283–306. [9] Y.X. Gu, G.D. Cheng, Structural modeling and sensitivity analysis of shape optimization, Struct. Optimization 6 (1) (1993) 29–37. [10] Y.X. Gu, H.W. Zhang, Z.Q. Guan, Z. Kang, Y.P. Li, W.X. Zhong, New generation software of structural analysis and design optimization—JIFEX, Struct. Engrg. Mech. 7 (6) (1999) 589–599. [11] Y.X. Gu, B.S. Chen, H.W. Zhang, R.V. Grandhi, A sensitivity analysis method for linear and nonlinear transient heat conduction with precise time integration, Struct. Multidisciplinary Optimization 24 (1) (2002) 23–37. [12] G.J.W. Hou, J.S. Sheen et al., Shape-sensitivity analysis and design optimization of linear thermoelastic solids, AIAA J. 30 (2) (1993) 528–537. [13] J.H. Huang, G.M. Fadel, V.Y. Blouin, M. Grujicic, Bi-objective optimization design of functionally gradient materials, Mater. Des. 23 (2002) 657–666. [14] J.H. Kim, G.H. Paulino, Isoparometric graded finite elements for nonhomogeneous isotropic and orthotropic materials, J. Appl. Mech. 69 (2002) 502–514. [15] M. Koizumi, FGM activities in Japan, Composites, Part B 128B (1997) 1–4. [16] S. Kok, N. Stander, W. Roux, Thermal optimization in transient thermoelasticity using response surface approximations, Int. J. Numer. Methods Engrg. 43 (1) (1998) 1–21. [17] B.Y. Lee, B.M. Kwak, Shape optimization of two-dimensional thermoelastic structures using boundary integral equation formulation, Comput. Struct. 41 (4) (1991) 709–722. [18] B.Y. Lee, B.M. Kwak, Axisymmetric thermoelastic shape sensitivity analysis and its application to turbine disc design, Int. J. Numer. Methods Engrg. 33 (10) (1992) 2073–2089. [19] Q. Li, G.P. Steven, Y.M. Xie, Displacement minimization of thermoelastic structures by evolutionary thickness design, Comput. Methods Appl. Mech. Engrg. 179 (1999) 361–378. [20] T. Liu, Design optimization and sensitivity analysis of thermal structure, Master thesis, Department of Engineering Mechanics, Dalian University of Technology, China, 2003. [21] A.J. Markworth, K.S. Ramesh, W.P. Parks Jr., Review modeling studies applied to functionally graded materials, J. Mater. Sci. 30 (1995) 2183–2193. [22] A.J. Markworth, J. Saunders, A model of structure optimization for a functionally graded material, Mater. Lett. 22 (1995) 103– 107. [23] R.A. Meric, Material and load optimization of thermoelastic solids Part I: sensitivity anslysis, J. Therm. Stress. 9 (1986) 359–372. [24] R.A. Meric, Material and load optimization of thermoelastic solids Part II: numerical results, J. Therm. Stress. 9 (1986) 373–389. [25] R.A. Meric, Shape optimization of thermoelastic solids, J. Therm. Stress. 11 (1986) 187–206. [26] P. Michaleris, D.A. Tortorelli, C.A. Vidal, Analysis and optimization of weakly coupled thermoelastoplastic systems with applications to weldment design, Int. J. Numer. Methods Engrg. 38 (1995) 1259–1285. [27] N. Noda, Thermal stresses in functionally graded materials, J. Therm. Stress. 22 (1999) 477–512. [28] Y. Ootao, R. Kawamra, Y. Tanigawa, R. Imamura, Optimization of material composition of nonhomogeneous hollow sphere for thermal stress relaxation making use of neural network, Comput. Methods Appl. Mech. Engrg. 180 (1999) 185–201. [29] Y. Ootao, Y. Tanigawa, T. Nakamura, Optimization of material composition of FGM hollow circular cylinder under thermal loading: a neural network approach, Composites, Part B 30 (1999) 415–422. [30] Y. Ootao, Y. Tanigawa, O. Ishimaru, Optimization of material composition of functionally graded plate for thermal stress relaxation using a genetic algorithm, J. Therm. Stress. 23 (2000) 257–271. [31] S. Suresh, A. Mortensen, Fundamental of functionally graded material: processing and thermomechnical behaviour of graded metal and metal–ceramic composites, IOM Communication, London, 1998. [32] K. Tanaka, Y. Tanaka, K. Enomoto, V.F. Poterasu, Y. Sugano, Design of thermoelastic materials using direct sensitivity and optimization methods reduction of thermal stresses in functionally gradient materials, Comput. Methods Appl. Mech. Engrg. 106 (1993) 271–284.

B. Chen, L. Tong / Comput. Methods Appl. Mech. Engrg. 194 (2005) 1891–1911

1911

[33] K. Tanaka, Y. Tanaka, H. Watanabe, V.F. Poterasu, Y. Sugano, An improved solution to thermoelastic design in functionally gradient materials: scheme to reduce thermal stresses, Comput. Methods Appl. Mech. Engrg. 109 (1993) 377–389. [34] K. Tanaka, H. Watanabe, Y. Sugano, V.F. Poterasu, A multicriterial material tailoring of a hollow cylinder in functionally gradient materials: scheme to global reduction of thermoelastic stresses, Comput. Methods Appl. Mech. Engrg. 135 (1996) 369– 380. [35] Y. Tanigawa, Some basic thermoelastic problems for nonhomogeneous structural materials, Appl. Mech. Rev. 48 (6) (1995) 287– 300. [36] D.A. Tortorelli, G. Subramani, S. Lu, R.B. Haber, Sensitivity analysis for coupled thermoelastic systems, Int. J. Solids Struct. 27 (12) (1991) 1477–1497. [37] S. Turteltaub, Optimal control and optimization of functionally graded materials for thermomechanical processes, Int. J. Solid Struct. 39 (2002) 3175–3197. [38] S. Xu, R.V. Grandhi, Multipoint approximation development: thermal structural optimization case study, Int. J. Numer. Methods Eng. 48 (8) (2000) 1151–1164. [39] R.J. Yang, Shape design sensitivity analysis of thermoelasticity problems, Comput. Methods Appl. Mech. Engrg. 102 (1993) 41– 60.