Application of HFS-FEM to Functionally Graded Materials

I. INTRODUCTION

FGMs are a class of relatively new and promising composite materials that have optimized material properties by combining different material components following a predetermined law [1][2][3][4]. They are heterogeneous composite materials with graded variation of constituents from one material phase to another, which results in continuously varying material properties. FGMs thus have superior thermal and mechanical performance to conventional homogeneous materials, and have a wide variety of engineering applications especially for the purpose of removing mismatches of thermo-mechanical properties between coating and substrate and reducing stress level in structures.

Recently, two effective numerical methods have beed developed for analysing mechanical performance of FGMs [5,6]. The first is the so-called hybrid Trefftz FEM (or T-Trefftz method) [7][8][9]. Unlike in the conventional FEM, the T-Trefftz method couples the advantages of conventional FEM [10,11] and BEM [12,13]. In contrast to the standard FEM, the T-Trefftz method is based on a hybrid method which includes the use of an independent auxiliary inter-element frame field defined on each element boundary and an independent internal field chosen so as to a prior satisfy the homogeneous governing differential equations by means of a suitable truncated T-complete function set of homogeneous solutions. Since 1970s, T-Trefftz model has been considerably improved and has now become a highly efficient computational tool for the solution of complex boundary value problems. It has been applied to potential problems [14][15][16][17], two-dimensional elastics [18,19], elastoplasticity [20,21], fracture mechanics [22][23][24], micromechanics analysis [25], problem with holes [26,27], heat conduction [6,[28][29][30], thin plate bending [31][32][33][34], thick or moderately thick plates [35][36][37][38][39], threedimensional problems [40], piezoelectric materials [41][42][43][44][45], and contact problems [46][47][48].

On the other hand, the hybrid FEM based on the fundamental solution (F-Trefftz method for short) was initiated in 2008 [7,49] and has now become a very popular and powerful computational methods in mechanical engineering. The F-Trefftz method is significantly different from the T-Trefftz method discussed above. In this method, a linear combination of the fundamental solution at different points is used to approximate the field variable within the element. The independent frame field defined along the element boundary and the newly developed variational functional are employed to guarantee the inter-element continuity, generate the final stiffness equation and establish linkage between the boundary frame field and internal

II. Steady-state heat transfer in FGM

This section is concerned with the application of the T-Trefftz to the solution of Steady-state heat transfer in FGMs. A hybrid graded element model is described and used to analyse two-dimensional heat conduction problems in both isotropic and anisotropic exponentially graded materials.

II.1 Basic formulations

Consider a 2D heat-conduction problem defined in an anisotropic inhomogeneous media: (1) with the following boundary conditions:

-Specified temperature boundary condition (2) -Specified heat flux boundary condition (3) where denotes the thermal conductivity in terms of spatial variable and is assumed to be symmetric and positive-definite ( ). is the sought field variable and represents the boundary heat flux.

is the direction cosine of the unit outward normal vector to the boundary , and and are specified functions on the related boundaries, respectively. For convenience, the space derivatives are indicated by a comma, i.e. , and the subscript index takes values 1 and 2 in our analysis.

Moreover, the repeated subscript indices stand for summation convention.

II.2 Fundamental solution in FGMs

For simplicity, we assume the thermal conductivity varies exponentially with position vector, for example (4) where vector is a graded parameter and matrix is symmetric and positive-definite with constant entries. Substituting Eq (4) into Eq (81) yields (5) whose fundamental function defined in the infinite domain necessarily satisfies following equation ( (6) in which and denote arbitrary field point and source point in the infinite domain, respectively. is the Dirac delta function. The closed-form solution to Eq (6) in two dimensions can be expressed as [81]

where , is the geodesic distance defined as and .

is the modified Bessel function of the second kind of zero order. For isotropic materials, , , (5) recasts as (8) Then the fundamental solution given by (7) reduces to, 12 21 11 22 12 ,det 0

286

which agrees with the result in [82].

II.3 Generation of graded element

In this section, the procedure for developing a hybrid graded element model is described based on the boundary value problem (BVP) defined in Eqs (1)-(4). The focus is to fully introduce the smooth variation of material properties into element formulation, instead of stepwise constant approximation frequently used in the conventional FEM Similar to T-Trefftz FEM, the main idea of the F-Trefftz approach is to establish an appropriate hybrid FE formulation whereby intra-element continuity is enforced on a nonconforming internal displacement field formed by a linear combination of fundamental solutions at points outside the element domain under consideration, while an auxiliary frame field is independently defined on the element boundary to enforce the field continuity across inter-element boundaries. But unlike in the HT FEM, the intra-element fields are constructed based on the fundamental solution defined in Eq (7) (11) where is undetermined coefficients and is the number of virtual sources outside the element .

is the required fundamental solution expressed in local element coordinates , instead of global coordinates (see Fig. 2).

Figure 2

Intra-element field, frame field in a particular element in HFS-FEM, and generation of source pointsThe corresponding normal heat flux on is given

Evidently, Eq (11) analytically satisfies the heat conduction equation (5) due to the inherent property of .

In practice, the generation of virtual source points is usually done by means of the following formulation employed in the MFS [83][84][85]

where is a dimensionless coefficient, is the elementary boundary point and is the geometrical centroid of the element. For a particular element shown in Fig. 2, we can use the nodes of element to generate related source points for simplicity.

International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 287 by (13) where (14) with (15) II.3.2 Auxiliary conforming frame field In order to enforce the conformity on the field variable , for instance, on of any two

neighboring elements e and f, an auxiliary inter-element frame field is used and expressed in terms of the same degrees of freedom (DOF), , as used in the conventional finite elements. In this case, is confined to the whole element boundary (16) which is independently assumed along the element boundary in terms of nodal DOF , where represents the conventional FE interpolating functions. For example, a simple interpolation of the frame field on a side with three nodes of a particular element can be given in the form (17) where ( ) stands for shape functions in terms of natural coordinate defined in Fig. 3. The fundamental solution for FGM is used as in Eq (11) to approximate the intra-element field. It can be seen from Eq (9) that varied throughout each element due to different geodesic distance for each field point, so the smooth variation of material properties can be achieved by this inherent property, instead of stepwise constant approximation frequently used in the conventional FEM, for example, Fig. 4 illustrates the two models when the thermal conductivity varies along direction X 2 in isotropic material. It should be mentioned here that Eq (4) which describes variation of the thermal conductivity is defined under global coordinate system. When contriving the intraelement field for each element, this formulation has to be transferred into local element coordinate defined at the center of the element, the graded matrix in Eq (4) can, then, be expressed by (18) for a particular element e, where denotes the value of conductivity at the centroid of each element and can be calculated as follow: (19) where is the global coordinate of the element centroid. Accordingly, the matrix is used to replace (see Eq (7)) in the formulation of fundamental solution for FGM and the construction of intra-element field under local element coordinate for each element. (20) in which the governing equation (81) is assumed to be satisfied, a priori, in deriving the HFS-FE model. The boundary of a particular element consists of the following parts (21) where represents the inter-element boundary of the element 'e' shown in Fig. 1.

Figure 3

Typical quadratic interpolation for the frame field II.3.2 Graded element

Figure 1

II.4 Variational principle and stiffness equation

The stationary condition of the functional (20) can lead to the governing equation, boundary conditions and continuity conditions, which is shown here briefly. Eq (20) gives the following functional defined in a particular element:

whose first-order variational yields

From the notation and the Gauss theorem (24) for any smooth function in the domain, we have

For the displacement-based method, the potential conformity should be satisfied in advance, that is, (26) then, Eq (25) can be rewritten as (27) from which the Euler equation in the domain and boundary conditions on can be obtained (28) using the stationary condition .

II.4.2 Stiffness equation

Having independently defined the intra-element field and frame field in a particular element (see Fig. 2), the next step is to generate the element stiffness equation through a variational approach.

The variational functional corresponding to a particular element of the present problem can be written as

Appling the Gauss theorem (24) again to the above functional, we have the following functional for the HFS-FE model

Considering the governing equation (8), we finally have the functional defined on the element boundary only (31) which yields by substituting Eqs (11), (13) and (16) (34) from which the optional relationship between and , and the stiffness equation can be produced

III. Transient heat conduction in FGMs

III.1 Statement of heat conduction problems in FGMs

Consider a two-dimensional (2D) transient heat conduction problem: (36) with the boundary conditions: -Dirichlet boundary condition (37) -Neumann boundary condition (38) where denotes the time variable ( ). is the thermal conductivity dependent on the special variables .

is the number of dimensions of the solution domain ( in this study). is the mass density. is the specific heat, and is the unknown temperature field. represents the boundary heat flux defined by , where is the unit outward normal to the boundary . and are specified temperature and heat flow, respectively, on the related boundaries. In addition, an initial condition must be given for the time dependent problem. In this paper, a zero initial temperature distribution is considered, i.e.

The composition and the volume fraction of FGM constituents vary gradually with the coordinate X, giving a non-uniform microstructure with continuously graded macro-properties (conductivity, specific heat, density).

In the present discussion, to make the derivation is tractable, the mass density is assumed to be constant within each element and taken the value of  at the centroid of the element. The thermal conductivity and specific heat have been chosen to have the same functional variation so that the thermal diffusivity is constant, that is

and (42) It should be mentioned that the above assumption in FGMs leads to a class of solvable problems and can provide benchmark solutions to other numerical methods, such as FEM, meshless and BEM. Moreover, it can provide valuable insight into the thermal behavior of FGMs [86]. So this assumption has been followed by a lot of researchers in solving transient thermal problems in FGMs [4,86].

III.2 LT and fundamental solution in Laplace space

The LT of a function is defined by (43) where is the Laplace parameter. By integration by parts, one can show that: (44) The boundary conditions (37) and (38)

where is the non-homogeneity graded parameter.

Substituting Eq. (47) (51) In this case, the differential Eq. (50) in Laplace space becomes (52) Obviously, Eq.(52) is the modified Helmholz equation, whose fundamental solution is (53) Making use of Eq.(51), we obtain the fundamental solution of Eq. (50) in Laplace space (54) where , and denote arbitrary field point and source point in the infinite domain, respectively. is the modified Bessel function of the second kind of zero order.

III.2.2 General method for FGMs with different variation of properties

The method can be extended to a broader range of FGMs, not only exponential but also quadratic and trigonometric material variation, by variable transformations [86]. By defining a variable [86] (55)

Eq.(36) can be rewritten as (56) For simplicity, define (57) Then, Eq. (56) can be rewritten as (58) After performing the LT, the differential equation (58) becomes (59) When is a constant, Eq.(59) is a modified Helmholz equation whose fundamental solution is known. Then the fundamental solution of Eq. (36) in Laplace space can be obtained by inverse transformation: (60) For quadratic material, (61) In this case, in Eq. (59). For trigonometric material, (62) In this case, in Eq. (59).

For exponential material,

In this case, in Eq. (59). Substituting into Eq.(60) and using the exponential law, the fundamental solution given by Eq.(60) reduces to Eq. (54). Note that for quadratic, trigonometric and exponential variations of both heat conductivity and specific heat, the FGM transient problem can be transformed into the same differential equation which has a simple and standard form (Eq.(58)) [86].

III.3 Generation of graded element

As in the conventional FEM, the solution domain is divided into sub-domains or elements. For a particular element, say element e, its domain is denoted by and bounded by . Since a nonconforming function is used for modeling the internal fields, additional continuities are usually required over the common boundary between any two adjacent elements 'e' and 'f' (see Fig. 1) [41]: (64) in the proposed hybrid FE approach.

III.3.1 Non-conforming intra-element field

For a particular element, say element e, which occupies the sub-domain , the field variable within the element is extracted from a linear combination of fundamental solutions centered at different source points (see Fig. 2

where is undetermined coefficients and is the (65)) is used to approximate the intra-element field for a FGM. The smooth variation of material properties throughout an element can be achieved by using the fundamental solution which can reflect the impact of a concentrated unit source acting at a point on any other points. The model inherits the essence of a FGM, so it can simulate FGMs more naturally than the stepwise constant approximation, which has been frequently used in conventional FEM. Fig. 3 illustrates the difference between the two models when the thermal conductivity varies along direction X 2 .

Note that the thermal conductivity in Eq.(36) is defined in the global coordinate system. When contriving the intra-element field for each element, this formulation must be transferred into the local element coordinate system defined at the center of the element, and the graded heat conductivity in Eq.(40) can then be expressed by (66) for a particular element e, where denotes the value of conductivity at the centroid of each element and can be calculated as follows:

where is the global coordinates of the element centroid.

Accordingly, is used to replace (see Eq. (60)) in the formulation of the fundamental solution for the FGM and to construct the intra-element field in the local element coordinate system for each element. In practice, the generation of virtual sources is usually achieved by means of the following formulation employed in the MFS [4]

where is a dimensionless coefficient, and are, respectively, boundary point and geometrical centroid of the element. For a particular element as shown in Fig. 2, we can use the nodes of the element to generate related source points for simplicity.

The corresponding normal heat flux on is given by (69) where (70) with (71) III.3.2 Auxiliary conforming frame field

In order to enforce conformity on the field variable , for instance, on of any two

neighboring elements e and f, an auxiliary inter-element frame field is used and expressed in terms of nodal degrees of freedom (DOF), , as used in conventional FEM as (72) which is independently assumed along the element boundary, where represents the conventional FE interpolating functions. For example, a simple interpolation of the frame field on the side with three nodes of a particular element can be given in the form (73) where ( ) stands for shape functions which are the same as those in conventional FEM.

III.4 Modified variational and stiffness equation

Having independently defined the intra-element field and frame field in a particular element (see Fig. 2 (77) from which the optional relationship between and , and the stiffness equation can be produced in the form (78) where stands for the element stiffness matrix.

III.5. Numerical inversion of LT

In this section, we present a brief review of the inversion of the LT used in this work. In general, once the solution for in the Laplace space is found numerically by the method proposed above, inversion of the LT is needed to obtain the solution for in the original physical domain. There are many inversion approaches for LT algorithms available in the literature [87]. A comprehensive review on those approaches can be found in [88]. In terms of numerical accuracy, computational efficiency and ease of implementation, Davies and Martin showed that Stehfest's algorithm gives good accuracy with a fairly wide range of functions [89]. Therefore, Stehfest's algorithm is chosen in our study. If is the LT of , an approximate value of the function for a specific time is given by (79) where /2 min( , / 2) /2

in which must be taken as an even number. In implementation, one should compare the results for different N to investigate whether the function is smooth enough, and determine an optimum N [87]. Stehfest suggested and other researchers have found no significant change for [89]. Therefore, is adopted here. That means that for each specific time it is necessary to solve different boundary value problems with different corresponding Laplace parameters in Laplace space 10 times, then weight and sum the solutions obtained in Laplace space.

IV. F-Trefftz method for Nonlinear FGMs

IV.1 Basic formulations

Consider a two-dimensional (2D) heat conduction problem defined in an anisotropic inhomogeneous media: (81) For an inhomogeneous nonlinear functionally graded material, we assume the thermal conductivity varies exponentially with position vector and also be a function of temperature, that is (2) Iterative method Since the heat conductivity depends on the unknown function u , an iterative procedure is employed for determining the temperature distribution. The algorithm is given as follows:

1. Assume an initial temperature 0 u .

Calculate the heat conductivity in

IV.3 Generation of graded element

In this section, an element formulation is presented to deal with materials with continuous variation of physical properties. Such an element model is usually known as a hybrid graded element which can be used for solving the boundary value problem (BVP) defined in Eqs. (86) 95)) is used to approximate the intra-element field in FGM. It is well known that the fundamental solution represents the filed generated by a concentrated unit source acting at a point, so the smooth variation of material properties throughout an element can be achieved by this inherent property, instead of the stepwise constant approximation, which has been frequently used in the conventional FEM. For example, Fig. 4 illustrates the difference between the two models when the thermal conductivity varies along direction X 2 in isotropic material.

Note that the thermal conductivity in Eq. (87) is defined in the global coordinate system. When contriving the intra-element field for each element, this formulation has to be transferred into local element coordinate system defined at the center of the element, the graded matrix in Eq. (87) can, then, be expressed by

where is the global coordinates of the element centroid. Accordingly, the matrix is used to replace K (see Eq. (90)) in the formulation of fundamental solution for FGM and to construct intra-element field in the coordinate system local to element. In practice, the generation of virtual sources is usually done by means of the following formulation employed in the MFS [4] (98)

where is a dimensionless coefficient ( =2.5 in our analysis [4]), and are, respectively, boundary point and geometrical centroid of the element. For a particular element shown in Fig. 2, we can use the nodes of element to generate related source points. The corresponding normal heat flux on is given by where ( ) stands for shape functions in terms of natural coordinate defined in Fig. 3. where represents the inter-element boundary of the element 'e' shown in Fig. 1.

IV.4 Modified variational principle and stiffness equation

The stationary condition of the functional (104) can lead to the governing equation (Euler equation), boundary conditions and continuity conditions, details of the derivation can refer to Ref. [7].

IV.4.2 Stiffness equation

Having independently defined the intra-element field and frame field in a particular element (see Fig. 2), the next step is to generate the element stiffness equation through a variational approach and to establish a linkage between the two independent fields.

If the material is homogeneous, i.e., the Lame parameters are independent of the spatial variable , Eq. constants, is a graded parameter, which has a dimension of . In particular, if the graded parameter is equal to zero, the Lame constants in Eq. (121) will be reduced to two constants, and then the system of partial differential equations (118) will be the standard Navier-Cauchy equations for homogeneous isotropic elastic materials.

According to the work of [92], when the Poisson ratio is equal to 0.25 (a rather common value for rock materials) and the plane strain state is considered, one obtains (122) which can significantly simplify the derivation of fundamental solutions.

Generally, the free space fundamental displacement solution for an isotropic inhomogeneous elastic continuum must satisfy the following equation system In the present hybrid formulation, in order to obtain the element stiffness equation involving element boundary integrals only, the element interior displacement field is approximated by the linear

being the element stiffness matrix, which is sparse and symmetric.

VI. CONCLUSIONS

On the basis of the preceding discussion, the following conclusions can be drawn. This review reports recent developments on applications of HFS-FEM to functionally graded materials and structures. It proved to be a powerful computational tool in modeling materials and structures with inhomogeneous properties. However, there are still many possible extensions and areas in need of further development in the future. Among those developments one could list the following: