Microstructural optimization of a functionally graded transversely isotropic layer

Mechanics of Materials 31 (1999) 637±651

www.elsevier.com/locate/mechmat

Microstructural optimization of a functionally graded transversely isotropic layer J.C. Nadeau a,*, M. Ferrari b,1 a b

Department of Civil and Environmental Engineering, Box 90287, Duke University, Durham, NC 27708-0287, USA Department of Materials Science and Mineral Engineering, Department of Civil and Environmental Engineering, Biomedical Microdevices Center, University of California at Berkeley, Berkeley, CA 94720±1710, USA Received 30 July 1998; received in revised form 19 May 1999

Abstract Motivated by current problems in coating technology, this paper addresses the microstructural optimization of a layer which is free of tractions, transversely isotropic, in®nite and subjected to a prescribed thermal gradient. The layer's microstructure is characterized as a bi-constituent composite in the form of a continuous matrix perfectly bonded to embedded spheroidal short ®bers. Both constituents are assumed to be isotropic in both mechanical and thermal properties. The microstructural parameters are taken to be the volume fraction, aspect ratio, and orientation distribution of the ®bers. The composite layer is made functionally graded by assuming that the microstructural parameters vary through the thickness of the layer. The e€ective properties of the bi-constituent composite are given by the Mori± Tanaka, Hatta±Taya and Rosen±Hashin homogenization theories. The compositional and microstructural properties are determined such that an objective function de®ned in terms of strain energy and curvature is minimized. Speci®c results are presented for an aluminum (Al) layer reinforced with silicon carbide (SiC). Comparisons are made to conventional coating technology. Ó 1999 Elsevier Science Ltd. All rights reserved. Keywords: Functionally graded/gradient material (FGM); Spatially optimal microstructure; Bi-constituent composite material; Optimization; Thermomechanical layer/plate

1. Introduction Materials with a microstructure that varies on a macroscopic length scale are termed functionally gradient materials (FGMs) (Neubrand and R oedel, 1997; Suresh and Mortensen, 1997; Mortensen and Suresh, 1995; Rabin and Shiota,

* Corresponding author. Tel.: +1-919-660-5216; fax: +1-919660-5219; e-mail: [email protected] 1 Currently address: Biomedical Engineering Center, The Ohio State University, Columbus, OH 43210-1002, USA.

1995; Ilscher and Cherradi, 1995; Holt et al., 1993; Yamanouchi et al., 1990). Interest in FGMs stems from their applicability as structural materials operating in extreme thermal environments. For example, the High-Speed Civil Transport (HSCT) is a proposed Mach 2.4 space-plane which will incur leading edge and skin temperatures in the vicinity of 175 C and engine temperatures in the neighborhood of 1600 C. To address the structural and thermal issues associated with the HSCT it has been proposed to combine the strength and ductility properties of a metal with the excellent thermal properties of a ceramic.

0167-6636/99/$ - see front matter Ó 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 1 6 7 - 6 6 3 6 ( 9 9 ) 0 0 0 2 3 - X

638

J.C. Nadeau, M. Ferrari / Mechanics of Materials 31 (1999) 637±651

Of particular interest in this paper is, given a characterization of the microstructure, to determine the optimal spatial distribution of that microstructure in order to achieve an optimal response of a speci®c system for a given loading condition. Much of the work, however, in the area of FGM optimization has been in the area of determining the optimal distribution of the e€ective material properties (Bendsùe et al., 1995; Bendsùe et al., 1994; Thomsen et al., 1994; Salzar and Barton, 1994; Carman et al., 1993; Kheinloo, 1987a; Kheinloo, 1987b). Then, when such optimal macroscopic property variations are determined, it remains to determine the microstructural variations which correspond to the optimal macroscopic properties (Sigmund, 1994; Zuiker, 1995). A more fundamental approach is to consider a characterization of the microstructure of a material and optimize it directly since, ultimately, in the manufacturing process it is the microstructure that is controlled. By relating the e€ective properties to the underlying microstructure it is possible to optimize the microstructure. Further justi®cation for optimizing a characterization of a microstructure rather than the e€ective material properties themselves is the correlation between di€erent e€ective material properties a€ected by the microstructure. For example, the e€ective thermal expansion is a function of the e€ective compliance (Rosen and Hashin, 1970). When optimizing material property distributions directly this underlying connection is neglected. With respect to microstructural optimization of FGMs, Obata and Noda (1996) investigated the optimal discrete gradation of volume fraction in order to maximize the minimum ratio of material strength to induced stress while considering temperature dependent material properties. Similarly, Tanigawa et al. (1997) minimized the maximum thermal stress in a transiently heated cylinder. In both cases the volume fraction distribution was characterized by a single parameter. Tanaka et al. (1993) have investigated the optimal volume fraction distribution in a hollow circular cylinder to reduce the maximum stress. The volume fraction distribution through the wall thickness was completely characterized by two parameters and the e€ective properties were predicted by the rule of

mixtures. A similar analysis optimizing heat ¯ow has recently been performed by Markworth and Saunders (1995). Allaire and Kohn (1993) and Jog et al. (1994) optimize microstructural features in the form of periodic arrays of rectangular perforations. For this investigation, motivated by the application to the HSCT, we consider a thermomechanical system in the form of an in®nite, transversely isotropic layer subjected to a prescribed thermal gradient. Motivated by natural FGMs such as bamboo, we model the functionally graded (FG) layer as a bi-constituent composite in the form of a matrix with embedded spheroidal ®bers. The layer is made functionally graded by taking the microstructural parameters: ®ber volume fraction, ®ber aspect ratio (i.e., length-to-diameter ratio) and ®ber orientation distribution to be functions of the through thickness coordinate of the layer. At each point through the thickness we assume that the spatial gradient of the microstructural ®elds is small and that homogenization theories for a macroscopically homogeneous medium can be employed. We term this pointwise homogenization. For comparison purposes we also consider an optimally designed bi-coated layer in the form of a substrate of matrix material coated, in general, on the top and bottom surfaces with ®ber material. The bi-coated layer is in some ways an indicator of performance achievable by current technology and therefore serves as a benchmark for the FGM design. The performance of the layer, either FG or bicoated, is quanti®ed by a scalar-valued function termed the objective function. The optimal layer is de®ned by the microstructural ®elds which minimize the objective function. In this paper a function is considered that comprises of a linear combination of the square of the layer's curvature and the layer's macroscopic strain energy density. 2. A non-homogeneous, transversely isotropic layer The non-homogeneous, transversely isotropic layer under consideration is illustrated in Fig. 1.

J.C. Nadeau, M. Ferrari / Mechanics of Materials 31 (1999) 637±651

639

rium equations are now exactly satis®ed and it remains to satisfy compatibility and the boundary conditions on the lateral boundary of the layer. The compatibility equations in indicial notation read eijk elmn jm;kn ˆ 0; Fig. 1. A non-homogeneous, transversely isotropic, layer.

It is of ®nite thickness and of in®nite extent in the x1 ±x2 plane. The through thickness coordinate of the layer is x3 . The layer is subjected to traction free boundary conditions. The top (x3 ˆ a) and bottom (x3 ˆ b) surfaces of the layer are subjected to uniform temperature boundary conditions: …1† h…a† ˆ ha ; h…b† ˆ hb :

…2†

The lateral boundaries may be considered to be insulated. The layer is inhomogeneous through its thickness only. The steady-state temperature distribution h…x3 † is governed by the one-dimensional heat equation (with no internal heat generation) …k33 h;3 †;3 ˆ 0;

…3†

where k33 is the thermal conductivity in the direction of the x3 -axis. We consider only the steady-state solution in the formulation of the objective function. As a result it is possible to solve the heat equation independently of the mechanical problem. We hereby assume that in the foregoing the heat equation has been solved and the temperature distribution through the thickness of the layer is known. The equilibrium equations with zero body forces are rij;j ˆ 0;

…4†

where r is the stress tensor. Assuming that all ®eld quantities are independent of the in-plane coordinates, the equilibrium equations, in conjunction with the traction-free boundary conditions on the top and bottom surfaces of the layer, yield r13 …x3 † ˆ r23 …x3 † ˆ r33 …x3 † ˆ 0. The equilib-

…5†

where eijk is the permutation symbol. For a thermoelastic problem the total strain is given by ij ˆ Sijkl rkl ‡ aij dh;

…6†

where dh is the di€erence between the steady-state temperature distribution h and the stress-free reference temperature distribution hr . That is, dh :ˆ h ÿ hr . Substituting Eq. (6) into the compatibility equations (5) and simplifying for the problem at hand, yields a single, non-trivial compatibility equation: …Sabcd rcd ‡ aab dh†;33 ˆ 0:

…7†

Integration yields Sabcd rcd ˆ Pab x3 ‡ Qab ÿ aab dh;

…8†

where Pab and Qab are integration constants which shall be determined from the boundary conditions on the lateral surfaces of the layer. The physical interpretation of these integration constants are as follows: Pab are curvatures and Qab are in-plane strains at x3 ˆ 0. In general, pointwise satisfaction of the lateral boundary condition is impossible and so it will be satis®ed in an integral sense: Z b rab dx3 ˆ 0 …9† a

and Z b rab x3 dx3 ˆ 0: a

…10†

Due to rotational symmetry about the x3 -axis we conclude that P11 ˆ P22 ; P12 ˆ P21 ˆ 0; Q11 ˆ Q22 ; and Q12 ˆ Q21 ˆ 0. For notational convenience, let P :ˆ P11 ˆ P22 and Q :ˆ Q11 ˆ Q22 . Substituting the stresses from 8 into the integral boundary conditions (9) and (10) yields a system of two equations for the solution of P and Q:

640

J.C. Nadeau, M. Ferrari / Mechanics of Materials 31 (1999) 637±651

Fig. 2. A bi-coated layer.

38 9

C x3 dx3

Rb a

C dx3

…11†

where a :ˆ a11 ˆ a22 and C :ˆ C1111 ‡ C1122 2 ÿ2…C1133 † =C3333 . This completes the formulation of the mechanical problem. A more general solution which does not make the assumption of transverse isotropy has been given by Ferrari (1992). As a means of comparison we consider a bicoated layer (see Fig. 2) in the form of a substrate with perfectly bonded coatings on its top and bottom surfaces. The x3 ˆ 0 plane is taken to be coincident with the top surface of the coated layer (i.e., a ˆ 0). The layer has a total thickness h ˆ b ÿ a. The top and bottom coatings are of thickness vmh and …1 ÿ v†mh, respectively, where m is the volume fraction of coating material and v is the fraction of coating material on the top surface. Thus, 0 6 m; v 6 1. The material properties of the substrate and coatings will be distinguished by a superscript s and c, respectively. Note that the bi-coated layer is a special case of the FG layer. As a result, for a given choice of a performance quantifying function, the FG layer is guaranteed to perform at least as well as the coated layer.

3. Microstructural parameters: a, r, f …g† The microstructure of the FG layer is modeled by taking the layer to be a bi-constituent com-

posite in the form of a matrix embedded with spheroidal ®bers. The constitution of the matrix and ®ber materials is taken to be isotropic. The ®bers are assumed to be perfectly bonded to the matrix. The parameters which characterize the microstructure of this composite are: (i) the volume fraction a of the ®ber material, (ii) the aspect ratio r (i.e., length-to-diameter ratio) of the spheroidal ®bers and (iii) the orientation distribution f …g† of the ®bers. The orientation distribution function (ODF) is denoted by f …g† where g :ˆ …w1 ; /; w2 † is the Euler triad and w1 , / and w2 are the Euler angles. The triad g represents the relative orientation of a given ®ber relative to the global coordinate frame, x1 x2 x3 . The quantity f …g† dg is thus the probability of a ®ber orientation occurring within the ``interval'' ‰g; g ‡ dg†. There are two possible ways to proceed with quantifying the ODF. The ®rst way is to expand the arbitrary ODF f …g† in terms of generalized spherical harmonics Tlmn …g† (Viglin, 1961) f …g† ˆ

1 X l l X X lˆ0 mˆÿl nˆÿl

Clmn Tlmn …g†:

…12†

The coecients Clmn of this expansion are termed texture coecients. Since the ODF is a probability density function it is a real-valued, strictly positive function which is normalized to unity. These properties imply conditions on the texture coecients. For example, the texture coecients are bounded (Nadeau and Ferrari, 1998a) jClmn j 6 2l ‡ 1

…13†

and C000 ˆ 1. These bounds are necessary below in the context of constrained optimization. The truncation theorem (Ferrari and Johnson, 1988) states that the orientational average of a tensor T 2 Tr (Tr is the set of tensors of order r) is not a function of texture coecients Clmn for l > r. Since four (4) is the largest order of any tensor to be considered it follows from the truncation theorem that only Clmn for l 6 4 need to be considered. Since only even-ordered properties are being considered (i.e., thermal conductivity, thermal linear expansion, and elastic sti€ness) it

J.C. Nadeau, M. Ferrari / Mechanics of Materials 31 (1999) 637±651

follows that Clmn for l ˆ 1; 3; 5; . . . do not contribute to orientational averaging. As a result of the transverse isotropy of the layer and the material and geometric symmetry of the spheroidal ®bers it can be concluded that Clmn ˆ 0 for n 6ˆ 0 or m 6ˆ 0 (Ferrari and Johnson, 1989). It follows that the only texture coecients which contribute to the orientational averaging of tensors for the layer are C200 and C400 . In what is to follow, the expressions for the e€ective material properties involve orientational averaging as a result of the ODF. As a result of the above discussion, the portion of the ODF relevant to the orientational averaging is thus given by   1 3 f~…g†  f~…/† ˆ 1 ‡ C200 ‡ cos…2/† 4 4   9 5 35 00 ‡ cos…2/† ‡ cos…4/† : …14† ‡ C4 64 16 64 Depending on the values of C200 and C400 the function (14) may or may not be an ODF. The remaining texture coecients (Clmn , l ˆ 1; 3; 5; 6; 7; 8; . . .), which are not relevant to the orientational averaging, could then be determined such that the sum of this function and 14 is an ODF. It has not been proven, however, that given a choice of C200 and C400 (satisfying the known bounds) that there exists a choice of the remaining texture coecients such that the sum is an ODF. Second, rather than characterizing the orientation distribution of the ®bers by the two texture coecients C200 and C400 it is possible to use /, which is the angle that the axis of the ®ber makes with the normal to the plane of the layer. In this case, / ˆ /0 …x3 †. Because of the transverse isotropy of the layer, and the in®nite-fold rotation axis of the ®ber, the ODF is independent of the other two Euler angles, w1 and w2 . This particular ODF is given by f …g† ˆ

2 d…/ ÿ /0 †; sin…/0 †

…15†

where d is the Dirac delta function. The C200 and C400 texture coecients for the ODF (15) are given by

641

  1 3 C200 ˆ 5 ‡ cos…2/0 † ; 4 4   9 5 35 ‡ cos…2/0 † ‡ cos…4/0 † : C400 ˆ 9 64 16 64

…16† …17†

This choice of an ODF is less general than that using the two independent texture coecients, however, it better lends itself to physical interpretation. The layer is made functionally gradient by taking the microstructural parameter ®elds (i.e., a, r, C200 and C400 (or /)) to be functions of the layer's through-thickness coordinate x3 . Each point through the thickness is assumed to possess a representative volume element (RVE) which is effectively homogeneous. This is a valid assumption provided that the microstructural parameters do not vary signi®cantly over a length scale comparable to the sizes of the ®bers. The RVE at each point is then homogenized to yield the e€ective material properties at that point. The e€ective elastic properties are predicted by the Mori-Tanaka e€ective medium theory (Benveniste, 1987) while recognizing, however, that this theory has its limitations (Reiter et al., 1997; Zuiker and Dvorak, 1994; Ferrari, 1991; Christensen, 1990). The e€ective thermal conductivity is predicted by the Hatta±Taya theory (Hatta and Taya, 1985, 1986; Nadeau and Ferrari, 1995a) which is form equivalent to the Mori-Tanaka theory of e€ective elastic properties. The e€ective coecient of thermal expansion is obtained from an extension of the Rosen-Hashin result (Nadeau and Ferrari, 1995b; Rosen and Hashin, 1970), which relates the e€ective CTLE to the e€ective compliance. These homogenization theories are summarized in Appendix A. For notational convenience we make the following de®nition: z :ˆ x3 . For numerical implementation we approximate each of the characteristic functions, say, a…z†, r…z†, C200 …z†, C400 …z†, by continuous, piecewise linear functions. To achieve this we partition the interval ‰a; bŠ through the thickness of the layer (see Fig. 3): …a ˆ zp1 ; zp2 ; . . . ; zpi ; zpi‡1 ; . . . ; zpn ˆ b†; zpi

zpi‡1 .

where < al variables,

At each zpi ai , ri , …C200 †i

…18†

we assign the nodand …C400 †i , which

642

J.C. Nadeau, M. Ferrari / Mechanics of Materials 31 (1999) 637±651

Fig. 3. Layer discretization for numerical solution.

correspond to the nodal values of the volume fraction, aspect ratio, and the two texture coeand C400 , respectively. Let cients, C200 00 00 u 2 fa; r; C2 ; C4 g. The function u…z† is approximated as a piecewise linear function as follows, u…z†  ui ‡

z ÿ zpi …u ÿ ui †; ÿ zpi i‡1

zpi‡1

zpi 6 z 6 zpi‡1 :

…19†

The four characteristic functions are now de®ned by 4n scalar quantities. For notational convenience denote this set of 4n quantities by n. That is,

We consider a quadratic variation in the thermal conductivity within a temperature ®nite element. A typical ®nite element is depicted in Fig. 4 where the two solid circles () denote the terminal ends of the element and the hollow circle () denotes the geometric center of the element. Nodal temperatures are associated only with the solid circles while thermal conductivities: k~1 , k~2 and k~3 , are assumed to be known at the three circle locations. When the three (3) conductivities are known

n :ˆ fa1 ; . . . ; an ; r1 ; . . . ; rn ; …C200 †1 ; . . . ; …C200 †n ; …C400 †1 ; . . . ; …C400 †n g:

…20†

When each element of n takes on a speci®c, admissible value, each of the four characteristic functions are uniquely de®ned. The interval ‰zpi ; zpi‡1 Š will at times be referred to as a ``parameter element''.

4. Numerical solution methods In this section, we present some details concerning the numerical solution of the heat equation and the mechanical problem for the layer model discussed in Section 2.

Fig. 4. A 1-D thermal element with quadratic variation of the thermal conductivity.

J.C. Nadeau, M. Ferrari / Mechanics of Materials 31 (1999) 637±651

the quadratic form of the thermal conductivity is given by k~33 …~z† ˆ a ‡ b~z ‡ c~z2 ;

…21†

where a ˆ k~3 ; b ˆ …k~2 ÿ k~1 †=l;

…22†

c ˆ 2…k~1 ‡ k~2 ÿ 2 k~3 †=l2 :

…24†

…23†

Given the evaluation of the e€ective thermal conductivities at the three locations depicted in Fig. 4, the thermal conductivity ®eld, at the element level, is approximated by the quadratic function (21). The temperature interpolation function is taken as the exact thermal distribution for the speci®c quadratic variation of the element thermal conductivity. At the level of the layer, the discretized form of the heat equation is a tridiagonal system of equations which is then solved for the nodal temperatures. In order to solve the mechanical layer problem of Section 2, ®ve (5) integrals are required through the thickness of the layer (see 11). Romberg integration of order 2k is implemented to perform the required integrations. In order to utilize Romberg integration with k > 1 it is required that the integrand be suciently smooth. To assure this, integration through the thickness is partitioned into a sum of integrations over each thermal element. We note that C and a are very computationally intensive to compute and therefore it is desired to keep the number of evaluations of these functions to a minimum. 5. Microstructural optimization Optimization algorithms search an admissible space for the extrema of a scalar-valued function. In order to implement an algorithm to optimize the microstructure of the FG layer we quantify the ``performance'', or response, of the layer by means of a scalar-valued function. This function is typically referred to as an objective function, or a cost function in the optimization literature. It is assumed that the optimal performance of the layer is achieved when the objective function is

643

minimized. Optimization algorithms can only be certain to have found a local minimum within a speci®ed tolerance. Whether that local minimum is the global minimum is generally unknown. For this reason, when feasible, a given optimization should be initiated from a number of di€erent starting points. It is not our intention to put forth an extensive search for the global minimum of the objective function. A local minimum which achieves a substantial improvement over a bicoated layer is sucient. Since it may be desired, for example, to require the aspect ratio r to be uniform through the thickness the nodal quantities n are not necessarily independent. There exists, then, an independent set of quantities v  n. Optimization can be a myopic pursuit with detrimental consequences. The desire to minimize, or maximize, a single speci®c quantity can often lead to unfavorable consequences in other quantities which prior to optimization were probably deemed satisfactory. It is therefore essential to be aware of potential adverse consequences of the optimization procedure. This requires a comprehensive review of the optimization results to assure that the design of the system is not adversely affected by the optimization procedure. For eciency we restrict attention to continuous objective functions. Further eciency would require the evaluation of gradients of the objective function. Due to the complexity of the expressions for the e€ective properties it is intractable to analytically determine the gradient of any objective function that we might select. Numerical approximation to the gradients could be utilized but this requires additional function evaluations which we are trying to minimize. Therefore we restrict attention to algorithms which require only objective function evaluations, e.g., Powell's Direction Set method and the Downhill Simplex method. A problem with Powell's method is that it will fail when operating on an objective function with an asymptotic minimum. It is quite likely that any objective function for the layer will have a asymptotic minima due to the fact that when the e€ective properties are evaluated at a volume fraction of either zero (0) or one (1) the values for the aspect ratio, and two texture coecients are

644

J.C. Nadeau, M. Ferrari / Mechanics of Materials 31 (1999) 637±651

irrelevant. This rules out the use of Powell's Direction Set method. The optimization algorithm which has thus been implemented is the Downhill Simplex method. The Downhill Simplex method is for unconstrained optimization. Since there are of n (see, e.g., (13)): 0 6 ai 6 1;

…25†

0 6 ri < 1;

…26†

ÿ5 6 …C200 †i 6 5;

…27†

ÿ9 6 …C400 †i 6 9;

…28†

it is necessary to make use of the following transformation. For simplicity let ni 2 n only be a function of vj 2 v. It follows that ni will be bounded by ‡ nÿ i 6 ni 6 ni

…29†

for all vj 2 R if the transformation from v to n for ni is given, for example, by 2 ‡ ÿ n i ˆ nÿ i ‡ …ni ÿ ni † sin …vj †:

…30†

As a benchmark the optimal FG layer is compared with an optimally designed bi-coated layer. There are two (2) parameters which characterize the bicoated layer: the volume fraction m of the coating material and the fraction of coating material v on the top surface of the layer. It is these parameters: m and v, which are optimized for the bi-coated layer.

6. Results In this section optimization results are presented for a layer with matrix/substrate material of aluminum and the ®ber/coating material of silicon carbide. The isotropic material properties for Al and SiC are given in Table 1. Results are presented for the FG and bi-coated layers for the objective function N ˆ KP 2 ‡ W where P is the curvature of the layer, W is the macroscopic strain energy and K is a parameter which serves only to weight the relative impor-

Table 1 Aluminum (Al) and silicon carbide (SiC) isotropic material properties Material Bulk modulus (GPa) Shear modulus (GPa) CTLE (10ÿ6 /K) Thermal conductivity (W/m K)

Matrix

Fiber

Al

SiC

73 26 23.0 237

226 176 4.8 490

tance of P 2 and W. For example, for smaller values of K greater emphasis is placed on reducing the strain energy W while for larger values of K greater emphasis is placed on reducing the curvature. The applied temperature boundary conditions are ha ˆ 1 C and hb ˆ 0 C with a stress-free reference temperature of hr ˆ hb ˆ 0 C. The top of the layer is taken at z ˆ a ˆ 0, and the bottom of the layer is taken at z ˆ b ˆ h ˆ 0:01 m. The partitioning of the interval ‰a; bŠ for the discretization of n corresponds to 10 equally sized subintervals The partition used for the solution of the heat equation is a subpartition of the partition used for n. In particular, each subinterval of the partition used for n is partitioned into ®ve (5) equally sized subintervals. In other words, there are 10 equally sized parameter elements and 50 equally sized thermal elements through the thickness of the layer. With respect to the numerical evaluation of the objective function the following parameters were used. For the integrations Romberg integration of order 2k ˆ 10 with a convergence tolerance of 10ÿ5 was utilized. The convergence tolerance for the optimization algorithm was 10ÿ4 . In order to arrive at a problem statement of reasonable size the aspect ratio r is taken to be uniform through the thickness of the layer. In addition to reducing the size of the admissible domain, the physical understanding of the problem is enhanced by considering optimization of the Euler angle / rather than the two texture coecients, C200 and C400 . As a result the number of optimization parameters have been reduced from 4…11† ˆ 44 in n to 23 in v. The optimization results are presented in Fig. 5 in the form of the minimum objective function value N versus the parameter K. We note the

J.C. Nadeau, M. Ferrari / Mechanics of Materials 31 (1999) 637±651

645

Fig. 5. Objective function N ˆ KP 2 ‡ W versus penalty parameter K for the optimal FG and bi-coated layers.

Fig. 6. Optimal volume fractions versus penalty parameter K for the objective function N ˆ KP 2 ‡ W .

following two (2) limiting conditions on the optimal bi-coated layer. As K ! 0 the optimal bicoated layer converges to a homogeneous layer of SiC. (We also note that the optimal FG layer also tends toward a homogeneous layer of SiC). As the penalty parameter approaches zero the objective function strives to minimize the strain energy W. The strain energy associated with a homogeneous, traction-free layer subjected to a thermal gradient is zero. Since the curvature of a homogeneous  ÿ hb †=h, and since layer is given by P ˆ ÿa……h† a aSiC < aAl it follows that as K ! 0 the optimal FG and bi-coated layers converge to a homogeneous layer of SiC. As K ! 1 the optimal bi-coated layer is an Al substrate with SiC only on its top surface and in a quantity such that the curvature is zero (0). Using a root ®nding algorithm, it was determined that this corresponds to a SiC volume fraction m ˆ 0:61005 with v ˆ 1. As can be observed in Fig. 6 this is what the bi-coated layer approaches as K ! 1. From Fig. 5 it is uncertain as to what the FG layer converges to as K ! 1. Since the bi-coated layer is a special case of the FG layer it is known that the FG layer can achieve P ˆ 0 since it has been demonstrated that the bicoated layer can achieve such a response. It is quite possible, however, that the FG layer will have an asymptote lower than that for the bi-coated layer as a result of potentially being able to achieve a zero curvature with a reduced strain energy over

that for the optimal bi-coated layer. Indications from the optimal microstructural ®elds for large values of K are that this is the case; The FG layer (most likely) has a lower asymptote than that for a bi-coated layer as K ! 1. The optimal FG layer achieves the greatest percentage increase in performance over the optimal bi-coated layer at approximately K ˆ 105 . From Fig. 5 it can be deduced that the optimal bi-coated layer is sensitive to the value of the penalty parameter K only over approximately three orders of magnitude: …103 ; 106 †. This zone corresponds to the transition region between the two straight line portions of the objective function. As K ! 0 this forces the strain energy toward zero faster than P 2 and thus the objective function becomes a linear function of K with a slope equal to P 2 . As K ! 1 this forces P 2 to zero faster than W and thus the objective function becomes a constant equal to W. This point is also illustrated by considering Fig. 6 which contains a plot of the two ``microstructural'' parameters for the bi-coated layer (i.e., m, v). Note that for all values of K the parameter v is equal to unity. Recall that v ˆ 1 implies that all of the ®ber material is located on the top surface of the layer. The fact that the gradient dm=dK 6ˆ 0 primarily only over this zone indicates that the optimal bi-coated layer is sensitive to K only over this domain. The results presented in Fig. 6 for the optimal FG layer are for

646

J.C. Nadeau, M. Ferrari / Mechanics of Materials 31 (1999) 637±651

the volume fraction of ®ber material per unit surRb face area, i.e., a a dz. On the other hand, as seen from Fig. 5, the optimal FG layer is sensitive to the penalty parameter well over nine (9) orders of magnitude. The relative sizes of these domains of K for which the objective function exhibits a sensitivity is another indicator as to the ability of a FG layer to improve performance over that of a bi-coated layer. From the objective function itself it is not possible to determine the speci®cs of how the layer is responding. For this reason the curvature P, strain Q and strain energy W, corresponding to each of the optimal solutions presented in Fig. 5, are presented in Figs. 7±9, respectively. Note from Fig. 7 that the absolute value of the curvature jP j of the optimal FG layer is always less than that of the optimal bi-coated layer. From Fig. 9 is seen that for 10 < K < 104 the optimal FG layer can achieve a smaller objective function value than that of the bi-coated layer by incurring additional strain energy to achieve a smaller curvature. Over this range of K the optimal bi-coated layer is essentially a homogeneous layer of SiC with a strain energy of zero (0). For K > 104 the optimal FG layer has a smaller curvature, smaller strain energy and, obviously, a smaller objective function value.

Fig. 8. Optimal strain Q versus penalty parameter K for objective function N ˆ KP 2 ‡ W .

Fig. 9. Optimal strain energy per unit area W versus penalty parameter K for objective function N ˆ KP 2 ‡ W .

Fig. 7. Optimal curvature P versus penalty parameter K for objective function N ˆ KP 2 ‡ W .

Arriving at the optimal FG layer results indicated by the circles () was not trivial. Since one can only say with certainty that a numerical algorithm like the Downhill Simplex method has found a local minimum we performed three (3) optimization runs for each data point () corresponding to three (3) di€erent starting conditions for the optimization algorithm. The smallest of the three (3) minima was selected. This procedure resulted in each data point requiring 10,000 to 15,000 objective function evaluations consuming 3 to 5 hours of computation time on a SPARC 20.

J.C. Nadeau, M. Ferrari / Mechanics of Materials 31 (1999) 637±651

As was remarked above, it was not our intention to extensively search for the global minimum of the objective function but, more importantly, to demonstrate that at least local minima exist for which an FG layer can adequately out perform a bi-coated layer. More extensive global minimum searches should be conducted in speci®c design applications of microstructural optimization. We now present the microstructural and thermomechanical ®elds for a typical optimization. We choose the optimal solutions corresponding to K ˆ 105 . The optimal bi-coated layer was determined by a di€erent program. The distribution of the microstructural parameters: a and /, through the thickness of the layer for the optimal FG solution are presented in Figs. 10 and 11, respectively. The limits on the vertical axes of these plots correspond to the bounds on the corresponding microstructural parameter with the exception of the vertical axis for the aspect ratio. The optimal ®ber volume fraction distribution (see Fig. 10) is a non-trivial, gradual gradation from pure SiC at the top surface to pure Al at the bottom. The aspect ratio, taken to be constant through the thickness, is optimal at a value of approximately 1630 ± very slender ®bers with a length 1630 times its diameter. We note that the value of the objective function is not appreciably changed when the aspect ratio is reduced to a value of 20. The optimal ®ber orientation distribution

Fig. 10. Optimal volume fraction distribution for FG layer: N ˆ 105 P 2 ‡ W .

647

Fig. 11. Optimal Euler angle / distribution for FG layer: N ˆ 105 P 2 ‡ W .

(see Fig. 11) is also non-trivial. At the top surface of the layer the ®bers' rotation axis is perpendicular to the normal of the layer. The ®bers gradually rotate through the thickness so that their rotation axis is in the plane of the layer at approximately one quarter the way through the thickness (i.e., z  0:25). Toward the mid-thickness of the layer the ®bers are once again oriented perpendicular to the layer. Through the bottom half of the layer the ®bers change orientation gradually until, at the bottom of the layer, they are oriented in the plane of the layer. The thermomechanical ®elds: h, ,  , m , r and W, for the same optimal FG layer are presented in Figs. 12±17, respectively. Also presented in these ®gures, for comparison purposes, are the thermomechanical ®elds for the optimal bi-coated layer for K ˆ 105 . From Fig. 13 it is seen that the optimal FG layer has a greater total strain than the optimal bicoated layer. One might want to view this as an adverse, unintended consequence of the optimization. For further comparison, however, the total strain distributions for homogeneous Al and SiC layers are also presented. The discontinuities of  , m , r and W at z  0:0069 corresponds to the substrate/coating interface of the optimal bi-coated layer. Note that the maximum macroscopic stress of the FG layer is approximately half that of the maximum stress of

648

J.C. Nadeau, M. Ferrari / Mechanics of Materials 31 (1999) 637±651

Fig. 12. Optimal temperature distribution: N ˆ 105 P 2 ‡ W .

Fig. 14. Optimal N ˆ 105 P 2 ‡ W .

mechanical

strain

distribution:

Fig. 13. Optimal total strain distribution: N ˆ 105 P 2 ‡ W .

the bi-coated layer. Of great interest is Fig. 17 which illustrates that the optimal FG layer has a signi®cantly smaller strain energy distribution through the thickness of the layer (except over the interval (0.009,0.010)) and without a discontinuity. FGMs could thus alleviate the debonding problems associated with coatings while at the same time achieving better performance (debonding aside). 7. Closure This paper has addressed microstructural optimization of a FG layer. The microstructure was

Fig. 15. Optimal thermal strain distribution: N ˆ 105 P 2 ‡ W .

characterized by volume fraction, aspect ratio and orientation distribution. Results were presented for the objective function N ˆ KP 2 ‡ W . For comparison purposes, the optimal FG layer was compared with an equivalently optimized bi-coated layer. It was found that the FG layer could substantially better the performance of a bi-coated layer and over a broader range of the parameter K. In particular, for K ˆ 105 the optimal FG layer bettered the optimal bi-coated layer by 86.7%. More speci®cally, the optimal FG layer exhibited a non-

J.C. Nadeau, M. Ferrari / Mechanics of Materials 31 (1999) 637±651

Fig. 16. Optimal stress distribution: N ˆ 105 P 2 ‡ W .

649

This has been the ®rst investigation into the microstructural optimization of a FGM system characterized by volume fraction, aspect ratio and orientation distribution. It can be concluded, in particular, that even in the elastic regime it is possible for a FG layer to signi®cantly out perform a bi-coated layer. This is signi®cant since many investigations to date have dealt with FGMs in non-linear regimes since it is in this realm that FGMs are expected to perform well. For example, Qui and Weng (1991) have demonstrated signi®cant sensitivity of the plastic response of a biconstituent composite with randomly oriented ®bers to the shape of the ®bers. This indicates that maybe optimally designed FGMs in the inelastic regime can achieve even greater increases in performance than shown here for the elastic regime.

Acknowledgements JCN was supported by the National Science Foundation (NSF) grant MSS±9215671 and MF gratefully acknowledges the support of the NSF through the National Young Investigator Award in the Mechanics and Materials program. Appendix A. Homogenization The Mori±Tanaka theory (Benveniste, 1987) for the e€ective elastic sti€ness C  is given by Fig. 17. Optimal strain energy per unit surface area distribution: N ˆ 105 P 2 ‡ W .

trivial volume fraction distribution and a nontrivial ®ber orientation distribution through the thickness of the layer. With respect to the thermomechanical ®elds a signi®cant result is that the macroscopic ®elds:  , m , r and W, are continuous. The same ®elds for the optimal bi-coated layer exhibit discontinuities at the substrate/coating interface. In addition the maximum stress in the optimal FG layer (for K ˆ 105 ) is half that of the maximum stress in the bi-coated layer.

C  ˆ C m ‡ ah…C f ÿ C m † : Aig ;

…A:1†

where h
ig denotes orientational averaging (Ferrari and Johnson, 1989) over the ODF f …g† and iÿ1 h …A:2† A ˆ T : …1 ÿ a†I s ‡ ahTig ;  s  ÿ1 …A:3† T ˆ I ‡ E : S m : …C f ÿ C m † : The tensor I s is the fourth-rank, symmetric idens ˆ …dik djl ‡ dil djk †=2, dij is the tity tensor (i.e., Iijkl Kronecker delta). E is the fourth-rank Eshelby tensor and S is a compliance tensor. The Eshelby tensor E is a function of the shape of the ®ber (e.g., the spheroidal aspect ratio r). The numerical evaluation of the e€ective sti€ness tensor is greatly

650

J.C. Nadeau, M. Ferrari / Mechanics of Materials 31 (1999) 637±651

aided by invariant tensor-to-matrix mappings which map tensorial expressions like those above to matrix expressions involving standard matrix operations (Nadeau and Ferrari, 1998b). The Hatta±Taya theory (Hatta and Taya, 1985, 1986; Nadeau and Ferrari, 1995a) for the e€ective thermal conductivity k is given by k ˆ km ‡ ah…kf ÿ km †aig ;

…A:4†

where a ˆ t‰…1 ÿ a†i ‡ ahtig Šÿ1 ; ÿ1  t ˆ i ‡ eqm …kf ÿ km † :

…A:5† …A:6†

i is the second-rank identity tensor (i.e., iij ˆ dij ) and e is the second-rank Eshelby tensor (Hatta and Taya, 1986) and q is a thermal resistivity tensor. The Eshelby tensor e is a function of the shape of the ®ber. A convenient presentation of both the second- and fourth-rank Eshelby tensors for spheroidal inclusions can be found in Nadeau (1996). The generalized result of Rosen±Hashin (Nadeau and Ferrari, 1995b; Rosen and Hashin, 1970) for the e€ective coecient of thermal linear expansion a is given by a ˆ am ‡ …af ÿ am † : …S f ÿ S m †

ÿ1

ÿ1

: …S  ÿ S m † : …A:7†

This result is applicable to the class of composites utilized in this paper. This result, however, is not generally valid for all bi-constituent composites (Nadeau and Ferrari, 1995b). References Allaire, G., Kohn, R.V., 1993. Optimal design for minimum weight and compliance in plane stress using extremal microstructures. European Journal of Mechanics A/Solids 12 (6), 839±878. Bendsùe, M.P., Guedes, J.M., Haber, R.B., Petersen, P., Taylor, J.E., 1994. An analytical model to predict optimal material properties in the context of optimal structural design. Journal of Applied Mechanics 61, 930±937. Bendsùe, M.P., Di az, A.R., Lipton, R., Taylor, J.E., 1995. Optimal design of material properties and material distribution for multiple loading conditions. International Journal for Numerical Methods in Engineering 38, 1149±1170.

Benveniste, Y., 1987. A new approach to the application of Mori-Tanaka's theory in composite materials. Mechanics of Materials 6, 147±157. Carman, G.P., Averill, R.C., Reifsnider, K.L., Reddy, J.N., 1993. Optimization of ®ber coatings to minimize stress concentrations in composite materials. Journal of Composite Materials 27 (6), 589±612. Christensen, R.M., 1990. A critical evaluation for a class of micromechanics models. Journal of the Mechanics and Physics of Solids 38 (3), 379±404. Ferrari, M., Johnson, G.C., 1988. The equilibrium properties of a 6 mm polycrystal exhibiting transverse isotropy. Journal of Applied Physics 63 (9), 4460±4468. Ferrari, M., Johnson, G.C., 1989. E€ective elasticities of short®ber composites with arbitrary orientation distribution. Mechanics of Materials 8, 67±73. Ferrari, M., 1991. Asymmetry and the high concentration limit of the Mori-Tanaka e€ective medium theory. Mechanics of Materials 11, 251±256. Ferrari, M., 1992. Anisotropic layers with through-thickness thermal and material variations. Thermal Stresses 15, 439± 445. Hatta, H., Taya, M., 1985. E€ective thermal conductivity of a misoriented short ®ber composite. Journal of Applied Physics 58 (7), 2478±2486. Hatta, H., Taya, M., 1986. Equivalent inclusion method for steady state heat conduction in composites. International Journal Engineering Science 24 (7), 1159±1172. Holt, J.B., Koizumi, M., Hirai, T., Munir, Z.A., 1993. Functionally Gradient Materials, volume 34 of Ceramic Transactions, American Ceramic Society, Westerville, OH. Proceedings of the Second International Symposium on Functionally Gradient Materials, presented at the Third International Ceramic Science and Technology Congress, held in San Francisco, CA, November 1992, pp. 1±4 . Ilschner, B., Cherradi, N. (Eds.), 1995. Proceedings of the third international symposium on structural and functionally graded materials. Lausanne, Switzerland. Presses Polytechniques et Universitaires Romandes. Jog, C.S., Haber, R.B., Bendsùe, M.P., 1994. Topology design with optimized, self-adaptive materials. International Journal for Numerical Methods in Engineering 37, 1323±1350. Kheinloo, M.L., 1987a. Optimization of material properties of inhomogeneous circular discs subjected to pressure. Prikladnaya Mekhanika 23 (8), 83±89. Kheinloo, M.L., 1987b. English translation: Soviet Applied Mechanics 23(8), 777±782. Markworth, A.J., Saunders, J.H., 1995. A model of structure optimization for a functionally graded material. Materials Letters 22, 103±107. Mortensen, A., Suresh, S., 1995. Functionally graded metals and metal-ceramic composites: Part 1 processing. International Materials Reviews 40 (6), 239±265. Nadeau, J.C., Ferrari, M., 1995a. Second-rank equilibrium and transport properties of ®brous composites: e€ective predictions and bounds. Composites Engineering 5 (7), 821± 838.

J.C. Nadeau, M. Ferrari / Mechanics of Materials 31 (1999) 637±651 Nadeau, J.C., Ferrari, M., 1995b. Thermal expansion of heterogeneous materials with application to low temperature environments. Report No. UCB/SEMM-95/08. Nadeau, J.C., Ferrari, M., 1998a. Bounds on texture coecients, Submitted for publication. Nadeau, J.C., Ferrari, M., 1998b. Invariant tensor-to-matrix mappings for evaluation of tensorial expressions. Journal of Elasticity 52 (1), 43±61. Nadeau, J.C., 1996. E€ective predictions of and bounds on properties of composite materials with application to microstructural optimization of a functionally gradient, textured, thermoelastic, composite layer. Ph.D. thesis, Univerity of California at Berkeley, Berkeley, CA. Neubrand, A., R oedel, J., 1997. Gradient materials an overview of a novel concept. Zeitschrift fur Metallkunde 88 (5), 358± 371. Obata, Y., Noda, N., 1996. Optimum material design for functionally gradient material plate. Archive of Applied Mechanics 66 (8), 581±589. Qui, Y.P., Weng, G.J., 1991. The in¯uence of inclusion shape on the overall elastoplastic behavior of a two-phase isotropic composite. International Journal of Solids and Structures 27 (12), 1537±1550. Rabin, B.H., Shiota, I., 1995. Functionally gradient materials. MRS Bulletin 20 (1), 14±55. Reiter, T., Dvorak, G.J., Tvergaard, V., 1997. Micromechanical models for graded composite materials. Journal of the Mechanics and Physics of Solids 45 (8), 1281±1302. Rosen, B.W., Hashin, Z., 1970. E€ective thermal expansion coecients and speci®c heats of composite materials. International Journal of Engineering Science 8, 157±173. Salzar, R.S., Barton, F.W., 1994. Residual stress optimization in metal-matrix composites using discretely graded interfaces. Composites Engineering 4 (1), 115±128.

651

Sigmund, O., 1994. Materials with prescribed constitutive parameters: an inverse homogenization problem. International Journal of Solids and Structures 31 (17), 2313± 2329. Suresh, S., Mortensen, A., 1997. Functionally graded metals and metal-ceramic composites: Part 2 thermomechanical behaviour. International Materials Reviews 42 (3), 85±116. Tanaka, K., Tanaka, Y., Enomoto, K., Poterasu, V.F., Sugano, Y., 1993. Design of thermoelastic materials using direct sensitivity and optimization methods. Reduction of thermal stresses in functionally gradient materials. Computer Methods in Applied Mechanics and Engineering 106, 271± 284. Tanigawa, Y., Oka, N., Akai, T., Kawamura, R., 1997. Onedimensional transient thermal stress problem for nonhomogeneous hollow circular cylinder and its optimization of material composition for thermal stress relaxation. JSME International Journal Series A 40 (2), 117±127. Thomsen, N.B., Wang, J., Karihaloo, B.L., 1994. Optimization±a tool in advanced material design. Structural Optimization 8 (1), 9±15. Viglin, A.S., 1961. A quantitative measure of the texture of a polycrystalline material-texture function. Soviet Physics: Solid State 2 (10), 2195±2207. Yamanouchi, M., Koizumi, M., Shiota, I. (Eds.), 1990. Proceedings of the ®rst international symposium on functionally gradient materials, Tokyo. Zuiker, J., Dvorak, G., 1994. The e€ective properties of functionally graded composites±I. Extension of the MoriTanaka method to linearly varying ®elds. Composites Engineering 4 (1), 19±35. Zuiker, J.R., 1995. Functionally graded materials: Choice of micromechanics model and limitations in property variation. Composites Engineering 5 (7), 807±819.