Heat conduction analysis of heterogeneous objects based on multi-color distance field

Materials and Design 31 (2010) 3331–3338

Contents lists available at ScienceDirect

Materials and Design journal homepage: www.elsevier.com/locate/matdes

Heat conduction analysis of heterogeneous objects based on multi-color distance field H.M. Zhou, Z.G. Liu *, B.H. Lu State Key Laboratory for Manufacturing Systems Engineering, Xi’an Jiaotong University, Xi’an, China

a r t i c l e

i n f o

Article history: Received 27 October 2009 Accepted 1 February 2010 Available online 6 February 2010 Keywords: Heat conduction Heterogeneous objects Multi-color distance field Meshless weighted least-square

a b s t r a c t The separation of computer aided design (CAD) modeling and finite element analysis (FEA) of heterogeneous objects makes it cumbersome for designers to exchange the data information. In this paper a systematic modeling and analysis for heterogeneous objects based on multi-color distance field was proposed. Firstly, the geometric model of an object was represented by the signed distance field; next, each grid was assigned with material property and material distance field which had the same resolution with geometric distance field was calculated according to material feature; the geometric and material distance fields were then unified, named as multi-color distance field. Then, through sampling the multi-color distance field, boundary and interior nodes were directly as nodes for meshless weighted least-square (MWLS) analysis. Using the moving least-squares approximation scheme as shape function, the MWLS analysis leads to a pure meshless analysis for heterogeneous objects. In the end, the effectiveness of the approach was illustrated by several numerical examples. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Heterogeneous objects can be classified into three types: multimaterial objects with clear material boundary, functionally graded materials (FGMs) and composite materials. In order to systemically design and analyze the heterogeneous objects, a lot of former researchers introduced many solutions. Voxel based model [1] and volume mesh based model [2] describe heterogeneous material distribution through intensive space decompositions with a collection of voxels and sub-volumes/polyhedrons, respectively. Kumar and Dutta [3] employed the extended r-sets to model and represent heterogeneous objects and Shin et al. [4] presented a representation scheme for heterogeneous objects based on the constructive representation. Kou and Tan [5] presented Heterogeneous Feature Tree structure to construct complex FGM distributions. Yang and Qian [6] presented the heterogeneous design and analysis method based on B-spline. Shapiro et al. [7] proposed the construction of field modeling using sampling distance field and interpolated physics field by way of extended R-functions. Toward reverse engineering, Zhou et al. [8] proposed the implicit modeling of heterogeneous object of multi-color distance field. Due to excellent thermal properties of heterogeneous objects, heat conduction problems in heterogeneous objects exist in many engineering processes, such as electronic cooling, encapsulation and cryogenics, and have received a considerable attention from * Corresponding author. Tel.: +86 13484461026; fax: +86 029 82660114. E-mail address: [email protected] (Z.G. Liu). 0261-3069/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.matdes.2010.02.001

the researchers. Various numerical techniques, such as the finite difference method (FDM) [9–11], finite element method (FEM) [12], boundary element method (BEM) [13] or more recently developed meshless methods, have been developed for analyzing the heat conduction problems because of the complexity of the corresponding governing equation, analytical solutions are usually difficult to obtain for those problems with arbitrary geometry and complex boundary conditions. Compared to FEM, FDM and BEM, the meshless methods is associated with a class of numerical techniques that approximate a given differential equation or a set of differential equations using global interpolations on an unstructured distribution of nodes, exhibiting the advantages of avoiding mesh generation, simple data preparation, easy post-processing and so on. Sutradhar and Paulino [13] presented a BEM for transient heat conduction in functionally graded materials without any domain discretization. However, the use of BEM usually results in domain integrals which may increase computing time and even cause some numerical problems. Liu and Gu [14] introduced meshless methods and their programming, such as the element-free Galerkin (EFG) method, the hp-clouds method, the meshless local Petrov–Galerkin (MLPG) method, meshless Galerkin method using radial basis functions, the least-square method and meshless point collocation method. The main advantage of the MLPG method [12,15,16] compared with regular Galerkin-based methods is that no background mesh is used to evaluate various integrals appearing in the local weak formulation of problem, but it requires a highorder quadrature rule to obtain converged results and thus needs

3332

H.M. Zhou et al. / Materials and Design 31 (2010) 3331–3338

much more computational effort in terms of CPU time than that for the FEM. A meshless algorithm of fundamental solution coupling with radial basis functions based on analog equation theory was proposed to conduct steady-state and transient thermal analysis of FGMs [17–18] and to simulate the static thermal stress distribution in 2D FGMs [19]. Katsikadelis [20] employed the meshless analog equation method to solve the 2D elastostatic problem for inhomogeneous anisotropic. Sladek et al. [16] presented the local boundary integral equation method with the MLS approximation for spatial variations of physical fields for stationary and transient heat conduction inverse problems in 2D and 3D axisymmetric body. Liu et al. [21,22] presented the meshless weighted leastsquares method (MWLS) to solve steady-state and transient-state heat transfer problems for homogeneous objects. The problem of optimization of FGMs can be formulated using this solution yielding the rate of change in the response with respect to design variables, i.e., performing a sensitivity analysis [23]. The element-free Galerkin analysis and genetic algorithm were applied to do single-objective [24] or multi-objective [25] optimization of the material distribution for thermo-mechanical processes. Chen and Tong [26] presented a systematic numerical technique for performing sensitivity analysis of coupled thermomechanical problem of FGMs. In the modeling of spatial variances of material properties, the graded FEM was employed to conduct the heat transfer analysis and structural analysis and their sensitivity analysis with known material property function. Na and Kim [27] employed FEM to optimize the volume fraction of FGMs considering stress and critical temperature. Gao et al. [28] studied topology optimization of steady heat conduction problem under both design-independent and design-dependent heat loads by means of a modified bidirectional evolutionary structural optimization method. These methods have mainly focused on simple topology heterogeneous objects with simple gradient material. It is difficult to use these methods to process arbitrary-shaped complex objects. Moreover, motivated to bridge the gap between CAD modeling and analysis of heterogeneous objects, this paper proposed a systematic approach to model and analyze the heterogeneous objects based on multi-color distance field. Simultaneously, the MWLS analysis was then extended to solve problems of heat conduction problem for heterogeneous objects combined with multi-color distance field modeling. The main contributions of this paper are: (1) The method allows the designers to design and analysis not only the geometry but also the material composition of an object. Objects with complex material distributions, which are hardly modeled with traditional approaches, can be modeled and analyzed. (2) Through sampling the multi-color distance field, boundary and interior nodes were directly as nodes for meshless weighted least-square (MWLS) analysis and no modification needs to be made. Compared with other meshless methods, the computing process is relatively simple and efficient.

Geometry scan

Design intent

Uniform subdividing

Material feature

Signed distance field

Material distance field Multi-color distance field

Node(coordinate and material properties) MWLS analysis Temperature and stress field Fig. 1. Modeling and MWLS analysis of multi-color distance field.

scanned or sampled to produce a 2D or 3D point cloud. The point cloud is then enclosed by a box which is uniformly subdivided to the prescribed resolution. The signed distance field of the model is then calculated by solving Eikonal equation [29] and zero isosurface of the model is extracted by Marching cube algorithm [30] for representing the model. At the same time, the material feature is collected by design intent or known material properties and then material distance field is calculated with the same resolution as the above signed distance field. The multi-color distance field is then formed by unifying the geometry field with material field. Sampling the zero iso-surface of the implicit model the boundary nodes can be obtained. The subdivided grid points inside the model are as the interior nodes. The temperature or stress field can be obtained by the MWLS analysis for the above boundary and interior nodes. 2.1. Geometric distance field Compute the distance field within a bounding box represented as an axis-aligned uniform 3D grid. Let the size of each voxel in the 3D Cartesian grid be Nx  Ny  Nz. The optimal grid size is determined by a specified reconstruction tolerance. For a given point cloud S2R3, an unsigned distance function is defined as the function that calculate the distance from point p (grid point) to the closest point x in S (see Fig. 2). The signed distance function is the distance to the boundary oS of the model, and the sign is used to denote inside or outside the model.

dS ðpÞ ¼ signðpÞ inf kx  pk

where signðpÞ ¼ 2. Multi-color distance field Modeling heterogeneous objects is explained to integrate material information into geometry and topology model, the corresponding material modeling mode should be selected for different modeling methods. Using implicit distance field to set up geometry model for those complex topology structure models, therefore the modeling of material information is also described by distance field in order to reduce data storage and convenient analysis. The modeling and analysis process of multi-color distance field is depicted in the flowchart in Fig. 1. At the outset, a model is

ð1Þ

x2@S



1

if p inside the model

1

otherwise

The calculation of unsigned distance is easy to implement but the sign assignment is very difficult for unorganized point cloud. However, the sign assignment is important for extracting iso-surface of model exactly. There are two commonly used methods to estimate the signs of distance field, that is, oriented tangent plane method [31] and solving Eikonal equation method [29]. Although the former is the most direct method to determine the signs of distance field, the computing process will be problematic and errorprone especially for regions with high curvature. In this paper, the latter is used to estimate the signs. The bounding cube of point

3333

H.M. Zhou et al. / Materials and Design 31 (2010) 3331–3338

where pT(x) is the basis function and the quadratic basis pT(x) = {1, x, y, x2, xy, y2} is used in this paper; a(x) is the coefficient, which is determined by minimizing a functional of weighted residual

J¼

N X



xI ðxÞ uh ðx; xI Þ  uðxI Þ

2

I¼1

¼

N X

xI ðxÞ½pT ðxÞaðxÞ  uðxI Þ2

ð5Þ

I¼1

The minimum value of J may be obtained through differentiating with respect to a(x)

" # N m X X @J xI ðxÞ pi ðxI Þai ðxÞ  uI pj ðxI Þ ¼ 0 j ¼ 1; 2; . . . ; m ¼2 @aj ðxÞ I¼1 i¼1

Fig. 2. 2D turbine blade (32  34 Cartesian grid).

cloud S is regularly subdivided into the desired accuracy sub-cubes and the distance function d(x) to S satisfies the Eikonal equation with boundary condition.

jrdðxÞj ¼ 1

ð2Þ

where dðx 2 SÞ ¼ 0. The Eikonal equation is solved using upwind schemes and Gauss–Seidel iterations with various orderings of the sweepings. Details of solving these PDEs can be found in Zeng et al. [32].

where xI are the positions of the N nodes, uI is the nodal parameter of the field variable at node I. xI(x) is the weighting function and usually a compactly supported function that is only nonzero in a small neighborhood called the ‘‘support domain” of node xI. The gauss function is used in this paper.

xðrÞ ¼

(

   expðr2 b2 Þ  expðb2 Þ = 1  expðb2 Þ

0

To represent the heterogeneous objects, the traditional Euclidean space R3 is extended to a tensor space: T = R3  Mn, where Mn denotes n-dim material space. The vector of material composition: M = [m1, m2, . . . ,mn], where, mi denotes primary material and P satisfies ni¼1 mi ¼ 1, that is, the total volume fraction of all material composition is summed to one. In this paper, the material distributions are defined in terms of the distances from the grid points to the material features. The total distance values of each grid point make up of material distance field. There are three major steps involved in calculating material distance field. In the first step, the suitable primary materials and material features are selected for the model to improve its performance. The second step involves calculating distance values from each grid point to material features. In the third step, for the multiple-material object, the several distance fields are blended. In other words, the distance field of material feature is selected as material distribution function, which is added to the above signed distance field of geometry model, then the multi-color distance field is formed. It is described as:

r>1

NðxÞ ¼ pT ðxÞA1 ðxÞBðxÞ 8 AðxÞaðxÞ ¼ BðxÞu > > > < N P xI ðxÞpðxI ÞpT ðxI Þ AðxÞ ¼ > I¼1 > > : BðxÞ ¼ ½x1 ðxÞpðx1 Þ x2 ðxÞpðx2 Þ:::xN ðxÞpðxN Þ

The shape functions used in MWLS analysis are based on a moving least-squares approximation scheme which is originally developed for the smooth interpolation of irregularly distributed function-value data. 3.1. The moving least-square (MLS) approximation The local approximation of the field variable u(x) is expressed as

ð4Þ

ð9Þ

In Eq. (7) the radius of the circular support domain, dmI, is chosen to make the matrix A(x) nonsingular everywhere in the domain, i.e., the support domain must have enough neighborhood nodes. 3.2. Steady heat conduction equation of heterogeneous objects The steady heat conduction equation and the thermal boundary conditions of heterogeneous objects, respectively, are

kðxÞr2 uðxÞ þ rkðxÞ  ruðxÞ þ Q ¼ 0

Prescribed flux :

3. Problem description and MWLS analysis

ð8Þ

where the matrices A(x) and B(x) are defined as

ð3Þ

where d(S) is the signed distance to model surface and F(r) is the distance to material feature in grid size Nx  Ny  Nz.

uðxÞ  uh ðxÞ ¼ NI ðxÞuI ¼ pT ðxÞaðxÞ

ð7Þ

Solving N(x) from Eqs. (4) and (6), the shape function is given by:

Prescribed temperature : ðNx  Ny  NzÞ

r61

r ¼ kx  x1 k=dml

2.2. Material distance field

8 3 > < R ! dðSÞ M n ! FðrÞ ; > :

ð6Þ

Convection :

 x 2 C1 uðxÞ ¼ u

 x 2 C2 n  kruðxÞ ¼ q

n  kruðxÞ ¼ hðu1  uðxÞÞ x 2 C3

ð10Þ ð11Þ ð12Þ ð13Þ

where u(x) is the temperature field on a fixed domain X surrounded by a closed boundary C = C1 + C2 + C3. The variable x denotes the physical dimensions expressed in Cartesian coordinates, x: (x,y). n is the outward surface normal. The parameters k, u1, h and Q are the thermal conductivity, ambient temperature, the heat transfer coefficient and the volumetric heat resource, respectively. Eq. (10) will degenerate into the heat conduction equation of homogeneous objects when thermal conductivity is a constant (independent of spatial variables) within the solution domain.

kr2 uðxÞ þ Q ¼ 0

ð14Þ

Substituting the approximate shape function u of Eq. (10) into Eq. (4), the residuals are minimized in a least-squares manner,

3334

d

Y

H.M. Zhou et al. / Materials and Design 31 (2010) 3331–3338

¼

Z

The system equations of the MWLS method for solving steadystate heat conduction equations are obtained as

duJ ½kr2 NJ ðxÞ þ rkrNJ ðxÞ  ½kr2 NI ðxÞuI Z  dC þ rkrN I ðxÞuI þ Q dX þ duJ N J ðxÞ½N I ðxÞuI  u X

þ

KU ¼ P

C1

Z

ÞdC duJ n  krNJ ðxÞ  ðn  krNI ðxÞuI  q

s¼1

C2

þ

Z

duJ ðn  krNJ ðxÞ þ hN J ðxÞÞ  ½n  krNI ðxÞuI

þ

C3

ð15Þ

We use an alternative discrete equation to avoid integration:

d

¼ duJ

N X

½ks r2 NJ ðxs Þ þ rks rNJ ðxs Þ  ½ks r2 NI ðxs ÞuI N1 X

N2 X

n  ks rNJ ðxs Þn  ks rNI ðxs Þ

s¼1

ks r2 NJ ðxs Þ þ rks rNJ ðxs ÞQ þ

N2 X

 NJ ðxs Þu

þ n  ks rNJ ðxs Þq

N3 X ½n  ks rNJ ðxs Þ þ hN J ðxs Þhu1

ð19Þ

s¼1

where N1, N2 and N3 is the number of interior nodes in the boundary C1, C2 and C3, respectively.

Þ n  ks rNJ ðxs Þ  ðn  ks rNI ðxs ÞuI  q

s¼1 N3 X ½n  ks rNJ ðxs Þ þ hN J ðxs Þ  ½n  ks rNI ðxÞuI

3.3. Material properties

s¼1

ð16Þ

There are two methods to describe the variance of the material properties. One is to use the specific functions for all kinds of material properties like Section 4.1. The other is to employ the specific functions of volume fraction of heterogeneous objects such as Sections 4.2 and 4.3. For those heterogeneous objects whose material composition is linear or non-linear gradients along the object’ thickness direction, so the material volume fraction is the function of distance field that is composed of two materials.

material1 : V 1 ðdÞ ¼ ðadÞk material2 : V 2 ðdÞ ¼ 1  ðadÞk

Fig. 3. Geometry of square region and boundary conditions.

N1 X s¼1

s¼1

 hðu1  NI ðxs ÞuI Þ

ð18Þ

s¼1

N3 X ½n  ks rN J ðxs Þ þ hNJ ðxs Þ  ½n  ks rNI ðxÞ þ hN I ðxÞ

N X

þ

 NJ ðxs Þ½NI ðxs ÞuI  u

s¼1

þ duJ

NJ ðxs ÞNI ðxs Þ þ

s¼1

þ rks rNI ðxs ÞuI þ Q  þ duJ N2 X

þ P¼

s¼1

þ duJ

N1 X s¼1

 hðu1  NI ðxÞuI ÞdC

Y

ð17Þ

N h i X ½ks r2 NJ ðxs Þ þ rks rN J ðxs Þ  ks r2 NI ðxs Þ þ rks rNI ðxs Þ K¼

) ð20Þ

where d is normalized distance from a query point to model surface (thickness of the heterogeneous object); a is the scaling factor; k is the distribution exponent to describe gradient property. If k is 1, material property is linear distribution. Usually a and k are determined by design intent or performance test. Details of defining material feature can be found in Zhou et al. [8]. In this paper, relevant material properties at some sampled points are then determined based on Voigt-type model [26] and Mori–Tanaka model [12].

Fig. 4. Temperature field distribution (a) 2D display (b) top boundary.

H.M. Zhou et al. / Materials and Design 31 (2010) 3331–3338

3.3.1. (1) Voigt-type model It is simply a linear rule of mixtures and the mixture property can be:

P ¼ P1 V 1 þ P2 V 2

ð21Þ

where P may be elastic modul E, Poisson’s ratio v, density q, coefficient of thermal expansion a and thermal conductivity k.

3.3.2. (2) Mori–Tanaka model It is modified rule of mixtures and the effective thermal conductivity k of the materials mixture can be determined as follows:

Table 1 Maximum of relative error. Cartesian grid

MWLS

FDM

10  10 20  20

0.15% 0.05%

0.35% 0.077%

k ¼ k1 þ

3k1 V 2 ðk2  k1 Þ 3k1 þ V 1 ðk2  k1 Þ

3335

ð22Þ

where k1 and k2 are the thermal conductivity of two materials.

4. Numerical results and discussions In order to demonstrate the efficiency and accuracy of the proposed method, we choose examples of relatively simple geometrical shape (in the Sections 4.1 and 4.2) and their results are compared with the analytical results and FDM. However, our approach could be applied to the modeling and analysis of more complicated case (in the Section 4.3), where the heat conduction in heterogeneous objects with arbitary-shape is discussed.

4.1. Case 1 A 100  100 m isotropic square region is illustrated in Fig. 3. The top and bottom boundaries are insulated. The left boundary is assigned a temperature of 200 °C and the right boundary is assigned a temperature of 100 °C. The spatial variation of the thermal conductivity is taken to be cubic in the x-direction as k(x, y) = (1 + x/100)3. An analytical solution is given as

" # 800 1 1 Tðx; yÞ ¼ þ 6 ð1 þ x=100Þ2 2

ð23Þ

The temperature field distribution is shown in Fig. 4. According as Wang et al. [17] the estimation rule is maximum of relative      100%. where f is the analytical solution, g is error% ¼ max f g f  ðx;yÞ2X

numerical solution and X is the domain within which f and g are defined. The results are shown in Table 1. From Table 1 it is obtained that the MWLS is better than FDM because of the smaller relative error. At the same time the maximum of relative error is also been influenced by the Cartesian grid size. Fig. 5. Geometry of square region and boundary conditions.

Fig. 6. Heterogeneous material temperature field analysis (top, right and bottom boundary).

3336

H.M. Zhou et al. / Materials and Design 31 (2010) 3331–3338

Fig. 7. Homogeneous material temperature field (top, right and bottom boundary).

Heat conduction results are shown in Figs. 6 and 7 for heterogeneous material and homogeneous material, respectively. It can be observed that a very good agreement between the MWLS and FDM in two homogenization models. This experiment also shows that using either Mori–Tanaka model (MT) or Voigt-type model (Voigt) to evaluate the effective thermal conductivities yields the similar results in calculation of the temperature fields. 4.3. Case 3

Fig. 8. Node and boundary condition.

4.2. Case 2 A 100  100 isotropic square region is illustrated in Fig. 5. The top and right boundaries are insulated. The left boundary is assigned a temperature of 0 °C and the bottom boundary is assigned a flux of 10. The spatial variation of the volume fraction is taken to be quadric in the x-direction as V(x, y) = (x/100)2.

In aerospace applications, in order to reduce the mechanical stress and thermal stress at the same time, metal on the low temperature side and ceramic on the high temperature side are used. Moreover, in order to reduce the high stress concentration at the interfaces of metal and ceramic, a mixture of metal and ceramic with varying graded material are designed [33]. Fig. 8 is node through sampling zero iso-surface and interior grid of Fig. 2. Using geometry distance field in the Section 2.1 the distance field distribution is obtained in Fig. 9. The material composition of turbine blade is varying along the thickness direction of model surface. Material field distribution is

Titanium : V 1 ðdÞ ¼ ðadÞk

) ð24Þ

Ceramics : V 2 ðdÞ ¼ 1  ðadÞk

In Fig. 10 material field distribution of titanium is computed in different a and k, where graded region is from d = 0 to d = 3 in Fig. 10(a) and from d = 0 to d = 5 in Fig. 10(b), Ceramic region is d > 3 in Fig. 10(a) d > 5 in Fig. 10(b). Heat conduction equation and boundary condition is

8 2 > < kðx; yÞr Tðx; yÞ þ rkðx; yÞrTðx; yÞ ¼ 0 k1 nrTðx; yÞ ¼ h1 ½T 1  Tðx; yÞ > : k2 nrTðx; yÞ ¼ h2 ½T 2  Tðx; yÞ

Fig. 9. Distance field distribution.

x; y 2 X x; y 2 C1

ð25Þ

x; y 2 C2

where, T1 = 1073 K, T2 = 1573 K. Material properties are shown in Table 2. The thermal conductivity distribution and temperature field distribution are shown in Figs. 11 and 12 through MWLS computation in Voigt-type model, respectively. As indicated in Table 2 ceramic has a very high thermal conductivity compared with titanium and it has the effect of promoting

H.M. Zhou et al. / Materials and Design 31 (2010) 3331–3338

3337

Fig. 10. Material field distribution of tianium (a) a = 1/3, k = 2; (b) a = 1/5, k = 2.

Table 2 Material heat properties.

5. Conclusion

Material

k (W/(m K))

h (W/(m2 K)

Titanium Ceramic

7 150

100 1500

the heats wide propagations. From Fig. 12a–b, it can be clearly seen the area with temperature below 1200 K is reduced and the area with temperature above 1200 K is greatly enlarged. This is reasonable because in Fig. 11b, the thermal conductivity in this region is designed to be higher and the heats are easier to propagate than the case in Fig. 11a.

In this paper, a meshless weighted least-square (MWLS) analysis for heterogeneous objects steady-state heat conduction analysis based on multi-color distance field was presented, in which the combination of modeling and heat conduction analysis using multi-color distance field provided a powerful numerical procedure. A discrete function was employed in the MWLS method to construct a set of linear equation, which avoided the burdensome task of numerical integration. Numerical cases not only showed that a good agreement was achieved between the results obtained from the proposed meshless method and available analytical solutions, but also showed that the method can take advantages of the de-

Fig. 11. Thermal conductivity distribution (a) a = 1/3, k = 2; (b) a = 1/5, k = 2.

Fig. 12. Temperature field in turbine blade (a) a = 1/3, k = 2; (b) a = 1/5, k = 2.

3338

H.M. Zhou et al. / Materials and Design 31 (2010) 3331–3338

signer’s experience in the formulation of complex material distribution. Finally, we list some issues that are suggested for further research. (1) The proposed method can be easily extended to a transientstate heat conduction problem through discretizing the time domain. (2) The proposed method can also be extended to some engineering problems of optimizing material distribution (such as inverse heat conduction problem and thermal stress problem) combined with a kind of optimization method. Acknowledgements The work described in this paper was supported by a Grant from the National Natural Science Foundation of China (Projects No. 50575177 and 50305027) and by Program for Changjiang Scholars and Innovative Research Team in University (IRT0646). The correlative members of the projects are hereby acknowledged. References [1] Chen M, Tucker JV. Constructive volume geometry. Comput Graph Forum 2000;19(4):281–93. [2] Jackson TR. Analysis of functionally graded material object representation methods. Ph.D. Thesis, Clemson Univerisity, South Carolina, USA; 2000. [3] Kumar V, Dutta D. An approach to modeling and representation of heterogeneous objects. J Mech Des 1998;120(4):659–67. [4] Shin KH, Natu H, Dutta D, Mazumder J. A method for the design and fabrication of heterogeneous objects. Mater Des 2003;24(5):339–53. [5] Kou XY, Tan ST. A systematic approach for integrated computer-aided design and finite element analysis of functionally-graded-material objects. Mater Des 2007;28:2549–65. [6] Yang PH, Qian XP. A B-spline-based approach to heterogeneous objects design and analysis. Comput-Aided Des 2007;39:95–111. [7] Biswas A, Shapiro V, Tsukanov I. Heterogeneous material modeling with distance fields. Comput-Aided Geom Des 2004;21:215–42. [8] Zhou HM, Liu ZG, Lu BH. Heterogeneous object modeling based on multi-color distance field. Mater Des 2009;30:939–46. [9] Turteltaub S. Optimal material properties for transient problems. Struct Multidisc Optim 2001;22:157–66. [10] Turteltaub S. Optimal control and optimization of functionally graded materials for thermomechanical processes. Int J Solids Struct 2002;39: 3175–97. [11] Turteltaub S. Functionally graded materials for prescribed field evolution. Comput Meth Appl Mech Eng 2002;191:2283–96. [12] Ching HK, Yen SC. Meshless local Petrov–Galerkin analysis for 2D functionally graded elastic solids under mechanical and thermal loads. Composites Part B 2005;36:223–40.

[13] Sutradhar A, Paulino GH. The simple boundary element method for transient heat conduction in functionally graded materials. Comput Meth Appl Mech Eng 2004;193:4511–39. [14] Liu GR, Gu YT. An introduction to meshfree methods and their programming. Berlin: Springer; 2005. [15] Gilhooley DF, Xiao JR, Batra RC, McCarthy MA, Gillespie Jr JW. Twodimensional stress analysis of functionally graded solids using the MLPG method with radial basis functions. Comput Mater Sci 2008;41:467–81. [16] Sladek J, Sladek V, Hon YC. Inverse heat conduction problems by meshless local Petrov–Galerkin method. Eng Anal Boundary Elem 2006;30:650–61. [17] Wang H, Qin QH, Kang YL. A new meshless method for steady-state heat conduction problems in anisotropic and inhomogeneous media. Arch Appl Mech 2005;74:563–79. [18] Wang H, Qin QH, Kang YL. A meshless model for transient heat conduction in functionally graded materials. Comput Mech 2006;38:51–60. [19] Wang H, Qin QH. Meshless approach for thermo-mechanical analysis of functionally graded materials. Eng Anal Boundary Elem 2008;32:704–12. [20] Katsikadelis JT. The 2D elastostatic problem in inhomogeneous anisotropic bodies by the meshless analog equation method (MAEM). Eng Anal Boundary Elem 2008;32:997–1005. [21] Liu Y, Zhang X, Lu MW. Meshless least-squares method for solving the steadystate heat conduction equation. Tsinghua Sci Technol 2005;10:61–6. [22] Liu Y, Zhang X, Lu MW. A meshless method based on least-squares approach for steady-and unsteady-state heat conduction problems. Numer Heat Transfer, Part B 2005;47:257–75. [23] Birman V, Byrd LW. Modeling and analysis of functionally graded materials and structures. Appl Mech Rev 2007;60:195–214. [24] Goupee AJ, Vel SS. Two-dimensional optimization of material composition of functionally graded materials using meshless analyses and a genetic algorithm. Comput Meth Appl Mech Eng 2006;195:5926–48. [25] Goupee AJ, Vel SS. Multi-objective optimization of functionally graded materials with temperature-dependent materials properties. Mater Des 2007;28:1861–79. [26] Chen BS, Tong LY. Thermomechanically coupled sensitivity analysis and design optimization of functionally graded materials. Comput Meth Appl Mech Eng 2005;194:1891–911. [27] Na KS, Kim JH. Optimization of volume fractions for functionally graded panels considering stress and critical temperature. Compos Struct 2009;89:509–16. [28] Gao T, Zhang WH, Zhu JH, Xu YJ, Bassir DH. Topology optimization of heat conduction problem involving design-dependent heat load effect. Finite Elem Anal Des 2008;44:805–13. [29] Zhao HK, Osher S, Merriman B, Kang M. Implicit and non-parametric shape reconstruction from unorganized data using a variational level set method. Comput Vis Image Underst 2001;80:295–314. [30] Theisel H. Exact isosurfaces for marching cubes. Comput Graph Forum 2002;21:19–31. [31] Hoppe H, DeRose T, Duchamp T, McDonald J, Stuetzle W. Surface reconstruction from unorganized points. In: SIGGRAPH’92 proceedings 1992; 26: 71–8. [32] Zeng HF, Liu ZG, Lin ZH. PDE-driven implicit reconstruction of 3D object. In: Proceedings of the computer graphics, imaging and vision: new trends; 2005. p. 251–6. [33] Qian XP, Dutta D. Design of heterogeneous turbine blade. Computer-Aided Des 2003;35:319–29.