Application of HFS-FEM to Functionally Graded Materials

© 2015 IJSRSET | Volume 1 | Issue 3 | Print ISSN : 2395-1990 | Online ISSN : 2394-4099 Themed Section: Engineering and Technology Application of HFS-FEM to Functionally Graded Materials Yi Xiao Research School of Engineering, Australian National University, Acton, ACT 2601, Australia ABSTRACT This paper presents an overview on applications of HFS-FEM to functionally graded materials. Recent developments on the hybrid fundamental solution (HFS) based finite element model (FEM) of steady-state heat transfer, transient heat conduction, nonlinear heat transfer, and elastic problems of functionally graded materials (FGMs) are described. Formulations for all cases are derived by means of modified variational functional and fundamental solutions. Generation of elemental stiffness equations from the modified variational principle is also discussed. Finally, a brief summary of the approach is provided. Keywords: Finite Element Method, Fundamental Solution, Functionally graded material 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-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-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-17], two-dimensional elastics [18, 19], elastoplasticity [20, 21], fracture mechanics [22-24], micromechanics analysis [25], problem with holes [26, 27], heat conduction [6, 28-30], thin plate bending [3134], thick or moderately thick plates [35-39], threedimensional problems [40], piezoelectric materials [4145], and contact problems [46-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 IJSRSET151364 | Received: 13 June 2015 | Accepted: 16 June 2015 | May-June 2015 [(1)3: 284-301] 284 field in the element. This review will focus on the FTrefftz method. spatial variable X and is assumed to be symmetric and positive-definite ( K12  K21 ,det K  K11 K22  K122  0 ). u is the sought ~ The F-Trefftz method, newly developed recently [7, 49], has gradually become popular in the field of mechanical and physical engineering since it is initiated in 2008 [7, 50, 51]. It has been applied to potential problems [16, 52-54], plane elasticity [19, 55, 56], composites [57-60], piezoelectric materials [61-63], three-dimensional problems [64], functionally graded materials [5, 65-67], bioheat transfer problems [68-72], thermal elastic problems [73], hole problems [74, 75], heat conduction problems [49, 76], micromechanics problems [25, 77], and anisotropic elastic problems [78-80]. Following this introduction, the present review consists of six sections. T-Trefftz FEM steady-state heat transfer in FGMs is described in Section 2. It describes in detail the method of deriving element stiffness equations. Section 3 focuses on the essentials of F-Trefftz elements for transient heat conduction in FGMs. The applications of F-Trefftz elements to heat transfer in nonlinear FGMs and elastic anaysis are discussed in Sections 4-5. Finally, a brief summary of the developments of the Treffz methods is provided. II. Steady-state heat transfer in FGM This section is concerned with the application of the TTrefftz 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: Kij u( X)  2u( X)  Kij 0 X i X j X i X j X  on u ~ on q ~ ~ ~ field variable and q represents the boundary heat flux. n j is the direction cosine of the unit outward normal vector n to the boundary   u   q , and u and q are specified functions on the related boundaries, respectively. For convenience, the space derivatives are indicated by a comma, i.e. u, j  u / X j , and the subscript index i, j 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 K ( X)  K exp(2β  X) (4) where vector β  ( 1 ,  2 ) is a graded parameter and matrix K is symmetric and positive-definite with constant entries. Substituting Eq (4) into Eq (81) yields Kij  i  j u ( X)  2i Kij  j u ( X)  0 (5) whose fundamental function defined in the infinite domain necessarily satisfies following equation Kij i  j N ( X, Xs )  2i Kij  j N ( X, Xs )   ( X, Xs )  0 (6) in which X and X s 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] N ( X, Xs )   K0 ( R) exp{β  ( X  Xs )} (7) 2 det K where   β  Kβ , R is the geodesic distance defined K 0 is the modified Bessel function of the second kind of zero order. For isotropic materials, K12  K21  0 , (2) －Specified heat flux boundary condition q  Kij u, j n j  q ~ (1) as R  R(X, X )  r  K 1r and r  r(X, X )  X  X . s s s with the following boundary conditions: －Specified temperature boundary condition u u ~ K11  K 22  k0  0 , (5) recasts as k0 2u ( X)  2k0 i  i u ( X)  0 (8) (3) Then the fundamental solution given by (7) reduces to, where Kij denotes the thermal conductivity in terms of N ( X, Xs )   K0 ( R) exp{β  ( X  Xs )} 2 k0 International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 285 (9) 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), rather than T-functions. Consequently, a variational functional corresponding to the new trial function is required to derive the related stiffness matrix equation. With the problem domain divided into some sub-domains or elements denoted by  e with the element boundary  e , additional continuities are usually required on the common boundary  Ief between any two adjacent elements ‘e’ and ‘f’ (see Fig. 1): ue  u f (conformity)   on  Ief  e   f (10) qe  q f  0 (reciprocity)  solution (MFS) [4] to remove the singularity of fundamental solution, for a particular element, say element e , which occupies sub-domain  e , we first assume that the field variable within an element is extracted from a linear combination of fundamental solutions centered at different source points (see Fig. 2), that is, ue  x    N e  x, y j  cej  N e  x  c e ns j 1 (11) where cej is undetermined coefficients and ns is the Ne  x, y j  is the required fundamental solution number of virtual sources outside the element e . expressed in local element coordinates ( x1 , x2 ) , instead of global coordinates ( X1 , X 2 ) (see Fig. 2). Evidently, Eq (11) analytically satisfies the heat conduction equation (5) due to the inherent property of Ne  x, y j  . In practice, the generation of virtual source points is usually done by means of the following formulation employed in the MFS [83-85] y  xb    xb  xc  (12) where  is a dimensionless coefficient, xb is the elementary boundary point and xc 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. in the proposed hybrid FE approach. e f  Ief Fig. 1 Illustration of continuity between two adjacent elements ‘e’ and ‘f’ II.3.1 Non-conforming intra-element field Activating by the idea of method of fundamental Fig. 2 Intra-element field, frame field in a particular element in HFS-FEM, and generation of source points The corresponding normal heat flux on  e is given International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 286 by where with qe  K e Q e  K e ue  Q ece n N e  K e ATe n Te   Ne,1 Ne,2  T A   n1 n2  (13) (14) (15) (11) to approximate the intra-element field. It can be seen from Eq (9) that N e 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 X2 in isotropic material. II.3.2 Auxiliary conforming frame field In order to enforce the conformity on the field variable u , for instance, ue  u f on  e   f of any two neighboring elements e and f, an auxiliary inter-element frame field u is used and expressed in terms of the same degrees of freedom (DOF), d , as used in the conventional finite elements. In this case, u is confined to the whole element boundary ue  x  Ne  x de (16) which is independently assumed along the element boundary in terms of nodal DOF de , where N e 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 u  N1u1  N 2u2  N 3u3 (17) where N i ( i  1, 2,3 ) stands for shape functions in terms of natural coordinate  defined in Fig. 3. Fig. 4 Comparison of computational cell in the conventional FEM and the proposed HFS-FEM 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 K in Eq (4) can, then, be expressed by K e (x)  K C exp(2β  x) (18) for a particular element e, where K C denotes the value of conductivity at the centroid of each element and can be calculated as follow: K C  K exp(2β  Xc ) (19) where X c is the global coordinate of the element centroid. Accordingly, the matrix K C is used to replace K (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. II.4 Variational principle and stiffness equation Fig. 3 Typical quadratic interpolation for the frame field II.3.2 Graded element The fundamental solution for FGM is used as N e in Eq II.4.1 Modified functional For the boundary value problem defined in Eqs (1)-(4), since the stationary conditions of the traditional potential or complementary variational functional can’t guarantee the satisfaction of inter-element continuity condition International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 287 required in the proposed HFS-FE model, a modified potential functional is developed as follows [49] 1    m     ku,i u,i d   qud    u  u  qd  e 2 qe e  e  (20) in which the governing equation (81) is assumed to be satisfied, a priori, in deriving the HFS-FE model. The boundary  e of a particular element consists of the following parts  e  ue   qe   Ie (21) where  Ie represents the inter-element boundary of the element ‘e’ shown in Fig. 1. 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:  m  [  e 1 ku,i u,i d   qud   q  u  u  d] qe e 2 e  m  [   ku,i u,i d   q ud (22) whose first-order variational yields  e qe  u   u  qd    u  u  qd]  e (23) e  h d   hni d  ,i   m  [  (ku,i ),i  ud   (24) for any smooth function h in the domain, we have e qe   q ud   q ud   e ue  Ie e  q  q   ud  u  u  qd] (25)  ue   u f on ue on  Ief ( u  u) ( ue  u f )  m  [  (ku,i ),i  ud   then, Eq (25) can be rewritten as e   q ud   e  Ie qe  q  q   ud  u  u   qd]  (26) (27) e boundary conditions on  e can be obtained u u in e on  qe on  e particular element e of the present problem can be written as 1  me    ku,i u,i d   qud   q  u  u  d qe e 2 e (29) Appling the Gauss theorem (24) again to the above functional, we have the following functional for the HFS-FE model 1  me    qud   u (ku,i ),i d  e  2  e (30)   qud   q  u  u  d e Considering the governing equation (8), we finally have the functional defined on the element boundary only  me   1 qud   qud   qud  qe e 2 e (31) which yields by substituting Eqs (11), (13) and (16) into the functional (31) 1  e   ceT H ece  d eT g e  ceT G ed e 2 (32) H e   QeT Ne d, G e   QeT Ne d, g e   NeT qd e e  qe (33) Next, to enforce inter-element continuity on the common element boundary, the unknown vector c e should be expressed in terms of nodal DOF de . The minimization from which the Euler equation in the domain  e and (ku,i ),i  0 qq The variational functional  e corresponding to a with For the displacement-based method, the potential conformity should be satisfied in advance, that is, u  0 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. qe From the notation q  ku,i ni and the Gauss theorem e using the stationary condition  me  0 . of the functional  e with respect to c e and de , respectively, yields e e  Hece  Gede  0,  G eTce  ge  0 (34) T ce dTe from which the optional relationship between c e and (28) de , and the stiffness equation can be produced ce = H e1G ed e and International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) K ed e = g e 288 (35) T 1 where K e = G e H e G e stands for the element stiffness and  matrix. III. Transient heat conduction in FGMs III.1 Statement of heat conduction problems in FGMs Consider a two-dimensional conduction problem: (2D) transient   (k (X )u ( X , t ))   ( X )c( X ) with the boundary conditions: －Dirichlet boundary condition u ( X , t ) (36) t u ( X , t )  u (X , t ) on u －Neumann boundary condition on  q q( X , t )  q (X , t ) heat k0 c0  (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 (37) (38) where t denotes the time variable ( t  0 ). k is the thermal conductivity dependent on the special variables The LT of a function u( X , t ) is defined by L(u( X , t ))  U ( X , s)   u( X , t )e st dt  0 where s is the Laplace parameter. By integration by parts, one can show that: u ( X , t ) ]  sU ( X , s )  u0 ( X ) t X   R . d is the number of dimensions of the solution domain  ( d  2 in this study).  is the The boundary conditions (37) and (38) become u ( X , t) mass density. c is the specific heat, and u is the U ( X , s)  on u s unknown temperature field. q represents the boundary q heat flux defined by q  ku / n , where n is the unit P( X , s)  (X , t ) on  q s outward normal to the boundary      . u and L[ d u q (43) (44) (45) (46) III.2.1 Exponentially graded material q are specified temperature and heat flow, respectively, First, we consider a FGM with thermal conductivity and 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. u ( X , 0)  u0 ( X )  0 (39) 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 k ( X )  k0 f ( X ) c( X )  c0 f ( X ) (40) (41) specific heat varying exponentially in one Cartesian coordinate, direction X 2 only, k ( X 2 )  k0 e 2  X 2 (47) c( X 2 )  c0e 2  X 2 where  is the non-homogeneity graded parameter. (48) Substituting Eq. (47) and Eq. (48) into Eq. (36) yields 1 u (49)  t denotes the derivative of u with respect to  2u  2  u X 2  where u X 2 X 2 ( u X 2  u / X 2 ) After performing the LT, Eq. (49) becomes  2U  2 U X 2   s U 0 (50) in LT space ,where u0 ( X )  0 (at t  0 ) is considered (see Eq.(39)). To obtain the fundamental solution of Eq. (50), following substitution is used here: International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 289 U  e  X 2 G (51) In this case, the differential Eq.(50) in Laplace space becomes  2G  (  2   s )G  0 (52) Obviously, Eq.(52) is the modified Helmholz equation, whose fundamental solution is G s 1 K0 (  2  r )  2 (53) Making use of Eq.(51), we obtain the fundamental solution of Eq. (50) in Laplace space N(X , XS )  1   ( X 2  X S2 ) s e K0 (  2  r ) 2  K 0 ( k '  r) 1  N(X , XS )  2 k ( X )1/2 k ( X s )1/2 s (60) For quadratic material, k ( X )  k0 (a1   X 2 ) 2 (61) k ( X )  k0 (a1 cos  X 2  a2 sin  X 2 ) 2 (62) In this case, k  0 in Eq.(59). For trigonometric material, ' In this case, k '   2 in Eq.(59). For exponential material, k ( X )  k0 (a1e  X 2  a2e   X 2 ) 2 (54) k '   2 in Eq.(59). Substituting where r  X  X S , X and X S denote arbitrary field In this point and source point in the infinite domain, respectively. K 0 is the modified Bessel function of the 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]. 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] v( X , t )  k ( X )u( X , t ) (55) case, (63) k '   2 into Eq.(60) and using the exponential law, 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  e and bounded by  e . Since a nonconforming function is used k ( X ) k ( X )  2 k ( X )  c( X ) v  v(  )v  (56) for modeling the internal fields, additional continuities k ( X ) t 4k 2 ( X ) 2k ( X ) are usually required over the common boundary  Ief Eq.(36) can be rewritten as 2 For simplicity, define k ( X ) k ( X )  2 k ( X )  k (X )  4k 2 ( X ) 2k ( X ) ' (57) between any two adjacent elements ‘e’ and ‘f’ (see Fig. 1)[41]: (58) in the proposed hybrid FE approach. Then, Eq. (56) can be rewritten as  2 v  k ' ( X )v  1 v  t After performing the LT, the differential equation (58) becomes  2V  k 'V   s V 0 (59) When k 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: ' U e  U f (conformity)   on  Ief  e   f (64) Pe  Pf  0 (reciprocity)  III.3.1 Non-conforming intra-element field For a particular element, say element e, which occupies the sub-domain  e , the field variable within the element is extracted from a linear combination of fundamental solutions centered at different source points (see Fig. 2), that is, U e  x    N e  x,xS  cej  N e  x  ce ns (65) j 1 where cej is undetermined coefficients and ns is the International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 290 Ne  x,xS  is the required fundamental solution number of virtual sources outside the element e. expressed in local element coordinates ( x1 , x2 ) , rather than global coordinates ( X1 , X 2 ) (see Fig. 2). Clearly, 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  e is given by Pe  ke Eq.(55) analytically satisfies the transformed governing equation of Eq.(36) in Laplace space due to the inherent property of Ne  x,xS  . where Q e   ke The fundamental solution for FGMs ( N e in Eq.(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 X2. 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 k ( X ) in Eq.(40) can then be with U e ni  Qece X j (69) N e ni  ke ATe X j Te   Ne,1 Ne,2  T (70) A   n1 n2  (71) III.3.2 Auxiliary conforming frame field In order to enforce conformity on the field variable U , for instance, U e  U f on  e   f of any two neighboring elements e and f, an auxiliary inter-element frame field U is used and expressed in terms of nodal degrees of freedom (DOF), d , as used in conventional FEM as Ue  x   Ne  x  de (72) which is independently assumed along the element boundary, where N e represents the conventional FE (66) 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 for a particular element e, where kC ( X ) denotes the (73) expressed by ke ( X )  kC ( X ) f ( X ) value of conductivity at the centroid of each element and can be calculated as follows: kC ( X )  k0 f ( X C ) (67) where X C is the global coordinates of the element centroid. Accordingly, kC is used to replace k (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] y  xb    xb  xc  (68) where  is a dimensionless coefficient, xb and xc are, respectively, boundary point and geometrical centroid of U  N11  N 2  2  N 3 3 where N i ( i  1, 2,3 ) 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), the element stiffness equation can be generated through a variational approach. The final functional defined only on the element boundary is  me   1 q PUd   Ud   PUd (74)   qe s e 2 e Substituting Eqs.(65), (69) and (72) into the functional (74), yields 1  e   ceT H ece  d eT g e  ceT G ed e 2 (75) where International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 291 He   QeT Ne d, G e   QeT Ne d, g e   NeT e e qe q d s (76) Next, to enforce inter-element continuity on the common element boundary, the unknown vector c e must be expressed in terms of nodal DOF de . The minimization of the functional  e with respect to c e and de , respectively, yields e e  Hece  Gede  0,  G eTce  ge  0 (77) T T ce de suggested N  10 and other researchers have found no significant change for 6  N  10 [89]. Therefore, N  10 is adopted here. That means that for each specific time T it is necessary to solve different boundary value problems with different corresponding Laplace parameters s  ln 2 i, i  1, 2 10 in Laplace T space 10 times, then weight and sum the solutions obtained in Laplace space. from which the optional relationship between c e and IV. F-Trefftz method for Nonlinear FGMs de , and the stiffness equation can be produced in the IV.1 Basic formulations Consider a two-dimensional (2D) heat conduction problem defined in an anisotropic inhomogeneous media: form ce = H e1G ed e and K ed e = g e (78) T 1 where K e = G e H e G e stands for the element stiffness  X 2 matrix. i , j 1 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 U ( X , s) in the Laplace space is found numerically by the method proposed above, inversion of the LT is needed to obtain the solution for u( X , t ) 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 F (s) is the LT of f (t ) , an approximate value f a of the function f (t ) for a specific time t  T is given by where Vi  (1) N / 2i fa   min( i , N / 2) k i 1 2 ln 2 N ln 2 Vi F ( i)  T i 1 T (79) N /2 k (2k )! (80) ( N / 2  k )!k !(k  1)!(i  k )!(2k  i)! in which N 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  ( Kij (X,u) i u ( X) )=0 X j X   (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 Kij (X, u)   (u) Kij exp(2β  X) ~ (82) where  (u)  0 is a function of temperature which may be different for different materials, the vector β  ( 1 ,  2 ) is a dimensionless graded parameter and matrix K  [ K ij ]1i , j  2 is a symmetric, positive-definite constant matrix ( K12  K 21 , det K  K11K 22  K12 2  0 ). The boundary conditions are as follows: －Dirichlet boundary condition u u －Neumann boundary condition q    Kij 2 i , j 1 u ni  q X j on u (83) on  q (84) ~ where Kij denotes the thermal conductivity which is the function of spatial variable X and unknown temperature field u . q represents the boundary heat flux. n j is the direction cosine of the unit outward normal vector n to the boundary   u   q . u and q are specified functions on the related boundaries, respectively. IV.2 Kirchhoff transformation and iterative method International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 292 Two methods are employed here to deal with the nonlinear term  (u) , one is Kirchhoff transformation [90] and another is the iterative method. (1) Kirchhoff transformation (u)   (u(X))    (u)du (85) Making use of Eq.(85), Eq.(81) reduces to  X 2 i , j 1 where  ( Kij* (X) i  ( X) )=0 X   X j Kij* (X)  Kij exp(2β  X) (86) (87) Substituting Eq.(87) into Eq.(86) yields  2   2 (X) K  2β  (K(X))  exp(2β  X)  0 (88)   ij X i X j i , j 1  where u   1 ( ) (89) It should be mentioned that the inverse of  in Eq.(89) exists since  (u)  0 . The fundamental solution to Eq.(88) in two dimensions can be expressed as [90] N ( X, Xs )   K0 ( R) exp{β  ( X  Xs )} (90) 2 det K where   β  Kβ , R is the geodesic distance defined as R  R(X,Xs )  r  K r and r = X - Xs in which X -1 and X s denote observing field point and source point in the infinite domain, respectively. K 0 is the modified Bessel function of the second kind of zero order. For K12  K21  0 , isotropic materials, K11  K 22  k0  0 , then the fundamental solution given by (90) reduces to N ( X, Xs )   K 0 ( R) exp{β  ( X + Xs )} (91) 2 k0 which agrees with the result in [82]. Under the Kirchhoff transformation, the boundary conditions (83)-(84) are transformed into the corresponding boundary conditions in terms of  .    (u ) on u (92) 2  u p   K ni    Kij ni  q  q on  q (93) X j X j i , j 1 i , j 1 2 * ij Therefore, by Kirchhoff transformation, the original nonlinear heat conduction equation (81), in which the heat conductivity is a function of coordinate X and unknown function u , can be transformed into the linear equation (86) in which the heat conductivity is just a function of coordinate X . At the same time, the field variable becomes  in Eq.(86), rather than u in Eq.(81). The boundary conditions (83)-(84) are correspondingly transformed into Eqs.(92)-(93). Once  is determined, the temperature solution u can be found by the reversion of transformation (89), i.e. u   1 ( ) . (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 u 0 . 2. Calculate the heat conductivity in Eq.(82) using u0 . 3. Solve the boundary value problem defined by Eqs.(81)-(84) for the temperature u 4. Define the convergent criterion u  u 0   (=10-6 in our analysis). If the criterion is satisfied, output the result and terminate the process. If not satisfied, go to next step. 5. Update u 0 with u 6. Go to step 2. 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) and (92)-(93). As was done in conventional FEM, the solution domain is divided into sub-domains or elements. For a particular element, say element e, its domain is denoted by  e and bounded by  e . Since a nonconforming function is used for modeling intraelement field, additional continuities are usually required over the common boundary  Ief between any two adjacent elements ‘e’ and ‘f’ (see Fig. 1):  e   f (conformity)   on  Ief  e   f pe  p f  0 (reciprocity)  International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 293 (94) in the proposed hybrid FE approach. IV.3.1 Non-conforming intra-element field For a particular element, say element e, which occupies sub-domain  e , the field variable within the element is extracted from a linear combination of fundamental solutions centered at different source points (see Fig. 2) that  e  x    Ne  x, y j  cej  Ne  x  ce x e , y j e ns j 1 (95) where cej is undetermined coefficients and ns is the centroid. Accordingly, the matrix K C 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] y  xb    xb  xc  (98) where  is a dimensionless coefficient (  =2.5 in our analysis[4]), xb and xc are, respectively, boundary expressed in terms of local element coordinates ( x1 , x2 ) , 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  e is given by instead of global coordinates ( X1 , X 2 ) (see Fig. 2). (99) Ne  x, y j  is the required fundamental solution number of virtual sources outside the element e. Obviously, Eq. (95) analytically satisfies the heat conduction equation (88) due to the inherent property of Ne  x, y j  . The fundamental solution for FGM ( N e in Eq.(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 X2 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 K* in Eq. (87) can, then, be expressed by K *e (x)  K C exp(2β  x) (96) for a particular element e, where K C denotes the value of conductivity at the centroid of each element and can be calculated as follows: K C  K exp(2β  Xc ) (97) pe  K *e where with Qe  K *e  e ni  Qece X j N e ni   AK *e Te X j Te   Ne,1 Ne,2  T (100) A   n1 n2  (101) IV.3.2 Auxiliary conforming frame field In order to enforce the conformity on the field variable u , for instance,  e   f on  e   f of any two frame field  is used and expressed in terms of nodal degrees of freedom (DOF), d , as used in the conventional finite elements as neighboring elements e and f, an auxiliary inter-element e  x  Ne  x de (102) which is independently assumed along the element boundary, where N e 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   N11  N 2  2  N3 3 (103) where N i ( i  1, 2,3 ) stands for shape functions in terms of natural coordinate  defined in Fig. 3. IV.4 Modified variational principle and stiffness equation where X c is the global coordinates of the element International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 294 IV.4.1 Modified variational functional For the boundary value problem defined in Eqs.(86) and (92)-(93), since the stationary conditions of the traditional potential or complementary variational functional can’t guarantee the satisfaction of interelement continuity condition required in the proposed HFS-FE model, a modified potential functional is developed as follows [7] 1 *  m    me  [   Kij  ,i  , j d e 2 e e (104)   q d   qe e      pd] in which the governing equation (86) is assumed to be satisfied, a priori, in deriving the HFS-FE model (For convenience, the repeated subscript indices stand for summation convention). The boundary  e of a particular element consists of the following parts  e  ue   qe   Ie (105) where  Ie represents the inter-element boundary of the element ‘e’ shown in Fig. 1. 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. The variational functional  e corresponding to a particular element e of the present problem can be written as 1 Kij*  ,i  , j d   2 e   q d   p    d  me    qe  e  (106) 1 pd    ( K ij*u,i ), j d      e 2 e (107)   q d   p    d  qe e   me   1 pd   q d   pd  qe e 2 e (108) which yields by substituting Eqs (95), (99) and (102) into the functional (108) 1  e   ceT H ece  d eT g e  ceT G ed e 2 (109) H e   QeT Ne d, G e   QeT Ne d, g e   NeT qd with e e qe (110) V. Elastic problems in FGMs V.1 Formulation of the problem In this section, basic equations and the corresponding fundamental solutions for FGMs presented in [91] are briefly reviewed to provide notations and references for the subsequent sections. V.1.1 Basic equations For a 2D linear elastic problem, the governing equations of force equilibrium in the absence of body forces are given by  ij , j  0 (111) where  ij are the components of the Cauchy stress tensor. For plane problems, all indices range from 1 to 2 and an index followed by a comma stands for partial differentiation with respect to the spatial coordinate. The summation convention is implied for repeated indices. For the functionally graded materials considered in this study, the elastic stiffness tensor Cijkl is associated with the spatial variable x  ( x1 , x2 ) ; that is, Cijkl  Cijkl (x) . Therefore, the linear elastic strain-stress  ij  Cijkl (x) kl relation is written as Appling the Gauss theorem to the above functional, we have the following functional for the HFS-FE model  me  only  Considering the governing equation (86), we finally have the functional defined on the element boundary (112) The components of stiffness tensor Cijkl must satisfy the Cijkl  Cklij  Cijlk  C jikl usual symmetric condition (113) Specially, for isotropic inhomogeneous elastic media, Cijkl (x)   (x) ij kl   (x)  ik jl   il jk  (114) the elastic stiffness tensor Cijkl is written as International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 295 where  ij is Kronecker’s delta, the Lame elastic parameters  (x) and  (x) are the functions of spatial coordinate variable x and can be expressed in terms of elastic modulus E , and Poisson ratio  as  ( x)  3  (x),  1  ( x)  In this work, the Lame constants   x   0  c  i xi  ,  and  are assumed   x   0  c  i xi  to be quadratic variation of the spatial variable x , that is 2 2 (121) , 0 and  0 are the corresponding material E ( x) where c (115) 2 1   constants,  i is a graded parameter, which has a with   3  4 for plane strain and   (3  )/(1  ) dimension of m1 . In particular, if the graded parameter i is equal to zero, the Lame constants in Eq. (121) will for plane stress. Therefore, the constitutive law Eq. (112) can be simplified as  ij   (x) ij kk  2 (x) ij (116) As well, the infinitesimal strain tensor  ij related to the 1  ui , j  u j , i  2 displacement field is expressed as  ij  (117) 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 0  0 (122) which can significantly simplify the derivation of Substituting Eq. (117) into the constitutive equation (116) fundamental solutions. Generally, the free space fundamental displacement and then into the equilibrium equation (111) we have solution for an isotropic inhomogeneous elastic ,i  x  u j , j  x     x     x  u j , ji  x  (118) continuum must satisfy the following equation system ,i  x  u j , j  x     x     x  u j , ji  x     x  ui , jj  x     x  ui , jj  x   , j  x  ui , j  x   u j ,i  x   0 (123)  , j  x  ui , j  x   u j ,i  x     x  xs  ei  0 If the material is homogeneous, i.e., the Lame where x is a field point in the infinite plane, xs is a parameters are independent of the spatial variable x , Eq. source point at which the unit force ei along the i(118) becomes (119) direction is applied, and  (x) is the Dirac delta function.        u j, ji  x  ui, jj  x  0 which is the classic Navier-Cauchy equation with respect to displacements. The boundary conditions have the same form as those of homogeneous materials: ui  ui on u ti  ti on t (120) where ti   ij n j represents the ith component of the boundary traction, and ni is the ith component of outward normal to the boundary.  u and t are the boundaries on which the displacement and the traction are prescribed respectively. An overbar denotes that the variable is specified. V.1.2 Fundamental solutions for quadratic variation of elasticity To obtain the fundamental displacement solution for the equilibrium equations (123) the following transformation is established for the displacement vector [91] wi   c  k xk  ui from which we have ui , j  wi , j c   k xk  (124) wi  j  c  k xk  (125) 2  ij   c  l xl  0 ij wk ,k  0  wi , j  w j ,i   and then the stress component can be given by  0 ij wk  k  0  wi  j  w j i  (126) Substituting Eqs. (124) and (121) into Eq. (123), we obtain 1 20 w j , ji   0 wi , jj   x  xs ei  0 (127)  c  k xk  International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com)   296 uli*  If the concentrated force acts at the origin, using the logarithmic potential, Yuan and Yin [91] obtain wli*  xi xl  1  2 li ln r  2  60  r  (128) * li where w denote the generalized displacement solutions at the field point x along the i-direction when a unit point force is applied at the origin along the l-direction. After this, with the inverse transformation of Eq. (124), and at the same time moving the point force from the origin to an arbitrary source point x s , the displacement components can be written as rr  c  uli*  2 li ln r  i 2l  s  r  60  c   k xk   c   k xk   (129) where some useful quantities related to the distance r are r   ri ri  , ri  xi  xis 1/2 r,i  ri , r,i r,i  1, ri , j   ij r (130) Based on the displacement formulation (129) [91], we obtain the strain components by differentiating the solution (129) with respect to the spatial variable x j  lij*   120  c   k xk   c   k xks  c 2 rj rl    ri rl      j  2 il ln r  2   i  2 jl ln r  2   r  r       2 ij rl   il rj   jl ri 4ri rl rj  c   4  s  r2 r  120  c   k xk   c   k xk   (131) and then, the stress components are given by    rr  r  c   2l ln r  k 2k l    c   k xk  l2   ij s   r r 6  c   k xk      * lij rj rl   rr     j  2 il ln r  i 2l   i  2 jl ln r  2  r  r     2 ij rl   il rj   jl ri 4ri rl rj   c   k xk    4 r2 r       (132) It is obvious that the fundamental solutions (129) and (132) can easily be reduced to the homogeneous fundamental solutions, when the graded parameters i  0 (i  1, 2) and c  1 . For example, for   0.25 ,   0    0 , we have homogeneous isotropic materials with Poisson ratio and  lij*  rr  1  2 il ln r  l 2i   r  60  (133) r 1  rl ij  rj il  ri jl 4rr  i 4j l   2 r r  6  (134) which are same as the formulations used in BEM for homogeneous materials. V.2 Hybrid finite element formulation V.2.1 Hybrid functional and element stiffness equation The initial concept of the hybrid finite element method features two independent fields (interior and frame fields) assumed in a given element. In the present work, the variables ui and ui respectively represent the interior and frame field variables. In the absence of body forces, the variational functional for any given element, say element e, used in the present model can be constructed as[7] 1  me    ij  ij d   ti ui d   ti  ui  ui  d (135) te e 2 e where  e is the domain of element e,  te and  ue are boundaries, where the traction and displacement are respectively specified, and  e denotes the whole element boundary. The inter-element boundary is denoted by  Ie . Clearly, for the hybrid element shown e  ue  te   Ie in Fig. 2 we have (136) Making use of Gauss theorem, the first-order variational of the functional can be further written as  me    ij , j ui d   ti ui d  te e  Ie  ti  ti  ui d    ti ui  ui  d (137) e  ij , j  0 , and the second integral enforces the in which the first integral gives the equilibrium equation reciprocity condition by co-considering those from neighboring elements. The traction boundary condition can be enforced by the third integral, and the final integral enforces equality of ui and ui along the elemental frame boundary  e . 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 International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 297 combination of the fundamental solutions at a series of de  , the application of Gauss theorem to the functional (135) gives 1 1 ui  x   c u  x, x    Ni ce  ,i, l  1, 2, m  1, 2, , M (138)  me     ij , j ui d   ti ui d e 2 2 e (146) where M is the number of virtual sources outside the  t u   s u  d d T te i i e i i c1M c2 M  is element domain, ce   c11 c21 Because the assumed displacement field (138) and an unknown coefficient vector (not nodal displacement), stress field (140) analytically satisfy the governing and the interpolation matrix equation (123), we have u1*i  x, xMs  u2*i  x, xMs  Ni   u1*i  x, x1s  u2*i  x, x1s  1  me    ti ui d   ti ui d   ti ui d (147) (139) te e 2 e s m source points x located outside the element domain as * lm li s m  consists of the fundamental solutions uli* x, x s  at M source points. It is noted that the constructed displacement field (138) can analytically satisfy the inhomogeneous elastic governing equation (123), since the fundamental solutions (129) of the problem are used as the interpolation functions. Making use of the strain-displacement equation (117) and the stress-strain relationship (112), the corresponding stress and traction components are expressed as  ij  x   clm lij*  x, x ms   Sij  ce  , i, j , l  1, 2, m  1, 2, , M (140) ti  x   clm tli*  x, xms   Qi ce  ,i, l  1, 2, m  1, 2, and Sij   1*ij  x, x1s   2*ij  x, x1s   in which Qi   t  x, x  * 1i s 1 t * 2i  x, x  s 1 t  x, x s M  t where 1 T T T ce  He ce   de  ge   ce  Ge de  (148) 2  H e    Qi   Ni  d G e    Qi  T e T g e    e  N i  d (149) ti  N i  d T te with respect to ce  and de  yields, respectively, the The stationary condition of the functional (148) optional relationship between ce  and de  and the ce  = He  Ge de  element stiffness equation as  x, x  (142) * 2i  me   1 , M (141) 1*ij  x, xMs   2*ij  x, xMs  * 1i Substituting Eqs. (138), (141) and (145) into the functional (147) yields tli*  x, x1s    lij*  x, xms  n j s M (143) with the traction kernels being defined by with Ke de  = ge  (151)  K e  =  G e   H e  G e  T 1 (152) being the element stiffness matrix, which is sparse and symmetric. (144) To enforce conformity of the displacement field on the common interface of any two neighboring elements, frame displacement fields ui are separately assumed on ui  x   dlk Nli  x, xk   [Ni ]de  ,i, l  1, 2, k  1, 2, and (150) the element boundary as , K (145) where [ N i ] denotes the interpolation vector relating the de  . boundary displacement to the nodal displacement vector optional relationship of unknown coefficient ce  and To obtain the element stiffness equation and the VI. CONCLUSIONS On the basis of the preceding discussion, the following conclusions can be drawn. This review reports recent developments on applications of HFSFEM 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: International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 298 [12] Q.H. Qin, Y. Mai, BEM for crack-hole problems in thermopiezoelectric materials, Engineering Fracture Mechanics, 69(5) (2002) 577-588. [13] Q.H. Qin, Y. Huang, BEM of postbuckling analysis of thin plates, Applied Mathematical Modelling, 14(10) (1990) 544-548. [14] W. Chen, Z. Fu, Q.H. Qin, Boundary particle method with high-order Trefftz functions, Computers, Materials & Continua (CMC), 13(3) (2010) 201-217. [15] H. Wang, Q.H. Qin, D. Arounsavat, Application of hybrid Trefftz finite element method to non ‐ linear problems of minimal surface, International Journal for Numerical Methods in Engineering, 69(6) (2007) 12621277. [16] H. Wang, Q.H. Qin, X. Liang, Solving the nonlinear Poisson-type problems with F-Trefftz hybrid finite element model, Engineering Analysis with Boundary Elements, 36(1) (2012) 39-46. [17] L. Cao, Q.H. Qin, N. Zhao, A new RBF-Trefftz meshless method for partial differential equations, IOP Conference Series: Materials Science and Engineering, 10 (2010) 012217. [18] Q.H. Qin, Dual variational formulation for Trefftz finite VII. REFERENCES element method of elastic materials, Mechanics Research Communications, 31(3) (2004) 321-330. [1] J.N. Reddy, Analysis of functionally graded plates, [19] H. Wang, Q.H. Qin, Numerical implementation of local International Journal for Numerical Methods in effects due to two-dimensional discontinuous loads using Engineering, 47(1-3) (2000) 663-684. special elements based on boundary integrals, [2] E. Carrera, S. Brischetto, A. Robaldo, Variable kinematic Engineering Analysis with Boundary Elements, 36(12) model for the analysis of functionally graded material (2012) 1733-1745. plates, AIAA journal, 46(1) (2008) 194-203. [20] Q.H. Qin, Formulation of hybrid Trefftz finite element [3] J.H. Kim, G.H. Paulino, Isoparametric Graded Finite method for elastoplasticity, Applied mathematical Elements for Nonhomogeneous Isotropic and modelling, 29(3) (2005) 235-252. Orthotropic Materials, Journal of Applied Mechanics, 69 [21] Q.H. Qin, Trefftz plane elements of elastoplasticity with (2002) 502-514. p-extension capabilities, Journal of Mechanical [4] H. Wang, Q.H. Qin, Y. Kang, A meshless model for Engineering, 56(1) (2005) 40-59. transient heat conduction in functionally graded [22] Y. Cui, Q.H. Qin, Fracture analysis of mode III problems materials, Computational Mechanics, 38(1) (2006) 51-60. by Trefftz finite element approach, in: WCCM VI in [5] H. Wang, L. Cao, Q.H. Qin, Hybrid graded element model conjunction with APCOM, 2004, pp. v1-p380. for nonlinear functionally graded materials, Mechanics [23] Y. Cui, Q.H. Qin, J.-S. WANG, Application of HT finite of Advanced Materials and Structures, 19(8) (2012) 590element method to multiple crack problems of Mode I, II 602. and III, Chinese Journal of Engineering Mechanics, 23(3) [6] Z.J. Fu, Q.H. Qin, W. Chen, Hybrid-Trefftz finite element (2006) 104-110. method for heat conduction in nonlinear functionally [24] Y. Cui, J. Wang, M. Dhanasek, Q.H. Qin, Mode III graded materials, Engineering Computations, 28(5) fracture analysis by Trefftz boundary element method, (2011) 578-599. Acta Mechanica Sinica, 23(2) (2007) 173-181. [7] Q.H. Qin, H. Wang, Matlab and C programming for [25] C. Cao, Q.H. Qin, A. Yu, Micromechanical Analysis of Trefftz finite element methods, New York: CRC Press, Heterogeneous Composites using Hybrid Trefftz FEM 2008. and Hybrid Fundamental Solution Based FEM, Journal [8] Q.H. Qin, The Trefftz finite and boundary element method, of Mechanics, 29(4) (2013) 661-674. WIT Press, Southampton, 2000. [26] M. Dhanasekar, J. Han, Q.H. Qin, A hybrid-Trefftz [9] Q.H. Qin, Trefftz finite element method and its element containing an elliptic hole, Finite Elements in applications, Applied Mechanics Reviews, 58(5) (2005) Analysis and Design, 42(14) (2006) 1314-1323. 316-337. [27] Q.H. Qin, X.Q. He, Special elliptic hole elements of [10] H.C. Martin, G.F. Carey, Introduction to Finite Element Trefftz FEM in stress concentration analysis, Journal of Analysis: Theory and Applications, McGraw-Hill Book Mechanics and MEMS, 1(2) (2009) 335-348. Company, New York 1973. [28] J. Jirousek, Q.H. Qin, Application of hybrid-Trefftz [11] Q.H. Qin, C.X. Mao, Coupled torsional-flexural vibration element approach to transient heat conduction analysis, of shaft systems in mechanical engineering—I. Finite Computers & Structures, 58(1) (1996) 195-201. element model, Computers & Structures, 58(4) (1996) [29] L. Cao, Q.H. Qin, N. Zhao, Application of DRM-Trefftz 835-843. and DRM-MFS to Transient Heat Conduction Analysis, 1 Development of efficient F-Trefftz FE-BEM schemes for complex engineering structures containing heterogeneous materials and the related general purpose computer codes with preprocessing and postprocessing capabilities. 2 Generation of various special-purpose elements to effectively handle singularities attributable to local geometrical or load effects (holes, cracks, inclusions, interface, corner and load singularities). The specialpurpose functions warrant that excellent results are obtained at minimal computational cost and without local mesh refinement. 3 Development of F-Trefftz FE in conjunction with a topology optimization scheme to contribute to microstructure design. 4 Extension of the F-Trefftz FEM to elastodynamics and fracture mechanics of FGMs. International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 299 Recent Patents on Space Technology (Open access), 2 (2010) 41-50. [30] N. Zhao, L. Cao, Q.H. Qin, Application of Trefftz Method to Heat Conduction Problem in Functionally Graded Materials, Recent Patents on Space Technology, 1(2) (2011) 158-166. [31] Q.H. Qin, Hybrid Trefftz finite-element approach for plate bending on an elastic foundation, Applied Mathematical Modelling, 18(6) (1994) 334-339. [32] Q.H. Qin, Postbuckling analysis of thin plates by a hybrid Trefftz finite element method, Computer Methods in Applied Mechanics and Engineering, 128(1) (1995) 123136. [33] Q.H. Qin, Transient plate bending analysis by hybrid Trefftz element approach, Communications in Numerical Methods in Engineering, 12(10) (1996) 609-616. [34] Q.H. Qin, Postbuckling analysis of thin plates on an elastic foundation by HT FE approach, Applied Mathematical Modelling, 21(9) (1997) 547-556. [35] F. Jin, Q.H. Qin, A variational principle and hybrid Trefftz finite element for the analysis of Reissner plates, Computers & structures, 56(4) (1995) 697-701. [36] J. Jirousek, A. Wroblewski, Q.H. Qin, X. He, A family of quadrilateral hybrid-Trefftz p-elements for thick plate analysis, Computer Methods in Applied Mechanics and Engineering, 127(1) (1995) 315-344. [37] Q.H. Qin, Hybrid-Trefftz finite element method for Reissner plates on an elastic foundation, Computer Methods in Applied Mechanics and Engineering, 122(3) (1995) 379-392. [38] Q.H. Qin, Nonlinear analysis of thick plates by HT FE approach, Computers & structures, 61(2) (1996) 271-281. [39] Q.H. Qin, S. Diao, Nonlinear analysis of thick plates on an elastic foundation by HT FE with p-extension capabilities, International Journal of Solids and Structures, 33(30) (1996) 4583-4604. [40] C.Y. Lee, Q.H. Qin, H. Wang, Trefftz functions and application to 3D elasticity, Computer Assisted Mechanics and Engineering Sciences, 15 (2008) 251-263. [41] Q.H. Qin, Variational formulations for TFEM of piezoelectricity, International Journal of Solids and Structures, 40(23) (2003) 6335-6346. [42] Q.H. Qin, Solving anti-plane problems of piezoelectric materials by the Trefftz finite element approach, Computational Mechanics, 31(6) (2003) 461-468. [43] Q.H. Qin, Mode III fracture analysis of piezoelectric materials by Trefftz BEM, Structural Engineering and Mechanics, 20(2) (2005) 225-240. [44] Q.H. Qin, Fracture Analysis of Piezoelectric Materials by Boundary and Trefftz Finite Element Methods, WCCM VI in conjunction with APCOM’04, Sept. 5-10, 2004, Beijing, China, (2004). [45] Q.H. Qin, Trefftz Plane Element of Piezoelectric Plate with p-Extension Capabilities, IUTAM Symposium on Mechanics and Reliability of Actuating Materials, (2006) 144-153. [46] Q.H. Qin, K.Y. Wang, Application of hybrid-Trefftz finite element method fractional contact problems, Computer Assisted Mechanics and Engineering Sciences, 15 (2008) 319-336. [47] K. Wang, Q.H. Qin, Y. Kang, J. Wang, C. Qu, A direct constraint ‐ Trefftz FEM for analysing elastic contact [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] problems, International journal for numerical methods in engineering, 63(12) (2005) 1694-1718. M. Dhanasekar, K.Y. Wang, Q.H. Qin, Y.L. Kang, Contact analysis using Trefftz and interface finite elements, Computer Assisted Mechanics and Engineering Sciences, 13(3) (2006) 457-471. H. Wang, Q.H. Qin, Hybrid FEM with fundamental solutions as trial functions for heat conduction simulation, Acta Mechanica Solida Sinica, 22(5) (2009) 487-498. C. Cao, Q.H. Qin, Hybrid fundamental solution based finite element method: theory and applications, Advances in Mathematical Physics, 2015 (2015) Article ID: 916029, 916038 pages. Q.H. Qin, Fundamental Solution Based Finite Element Method, J Appl Mech Eng, 2 (2013) e118. Z.J. Fu, W. Chen, Q.H. Qin, Hybrid Finite Element Method Based on Novel General Solutions for Helmholtz-Type Problems, Computers Materials and Continua, 21(3) (2011) 187-208. Y.T. Gao, H. Wang, Q.H. Qin, Orthotropic Seepage Analysis using Hybrid Finite Element Method, Journal of Advanced Mechanical Engineering, 2(1) (2015) 1-13. H. Wang, Q.H. Qin, Fundamental solution-based hybrid finite element analysis for non-linear minimal surface problems, in: E.J. Sapountzakis (Ed.) Recent Developments in Boundary Element Methods: A Volume to Honour Professor John T. Katsikadelis, WIT Press, Southampton, 2010, pp. 309-321. H. Wang, Q.H. Qin, Fundamental-solution-based hybrid FEM for plane elasticity with special elements, Computational Mechanics, 48(5) (2011) 515-528. H. Wang, Q.H. Qin, W. Yao, Improving accuracy of opening-mode stress intensity factor in two-dimensional media using fundamental solution based finite element model, Australian Journal of Mechanical Engineering, 10(1) (2012) 41-52. Q.H. Qin, H. Wang, Special Elements for Composites Containing Hexagonal and Circular Fibers, International Journal of Computational Methods, 12(4) (2015) 1540012. H. Wang, Q.H. Qin, Special fiber elements for thermal analysis of fiber-reinforced composites, Engineering Computations, 28(8) (2011) 1079-1097. H. Wang, Q.H. Qin, A fundamental solution based FE model for thermal analysis of nanocomposites, in: C.A. Brebbia, V. Popov (Eds.) Boundary Elements and Other Mesh Reduction Methods Xxxiii, WIT press, New Forest, 2011, pp. 191-202. H. Wang, Q.H. Qin, Implementation of fundamentalsolution based hybrid finite element model for elastic circular inclusions, in: Proceedings of the Asia-Pacific Congress for Computational Mechanics, 11th-14th Dec. 2013, Singapore. , 2013. C. Cao, Q.H. Qin, A. Yu, Hybrid fundamental-solutionbased FEM for piezoelectric materials, Computational Mechanics, 50(4) (2012) 397-412. C. Cao, A. Yu, Q.H. Qin, A new hybrid finite element approach for plane piezoelectricity with defects, Acta Mechanica, 224(1) (2013) 41-61. H. Wang, Q.H. Qin, Fracture analysis in plane piezoelectric media using hybrid finite element model, in: International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 300 S.W. Yu, X.Q. Feng (Eds.) Proceedings of the 13th International Conference of Fracture, China Science Literature Publishing House, Beijing, 2013. [64] C. Cao, Q.H. Qin, A. Yu, A new hybrid finite element approach for three-dimensional elastic problems, Archives of Mechanics, 64(3) (2012) 261–292. [65] L. Cao, Q.H. Qin, N. Zhao, Hybrid graded element model for transient heat conduction in functionally graded materials, Acta Mechanica Sinica, 28(1) (2012) 128-139. [66] L. Cao, H. Wang, Q.H. Qin, Fundamental solution based graded element model for steady-state heat transfer in FGM, Acta Mechanica Solida Sinica, 25(4) (2012) 377392. [67] H. Wang, Q.H. Qin, Boundary integral based graded element for elastic analysis of 2D functionally graded plates, European Journal of Mechanics-A/Solids, 33 (2012) 12-23. [68] H. Wang, Q.H. Qin, FE approach with Green's function as internal trial function for simulating bioheat transfer in the human eye, Archives of Mechanics, 62(6) (2010) 493-510. [69] H. Wang, Q.H. Qin, Computational bioheat modeling in human eye with local blood perfusion effect, in: J.H.T. E.Y.K. Ng, U.R. Acharya, J.S Suri (Ed.) Human Eye Imaging and Modeling, CRC Press, 2012, pp. 311-328. [70] H. Wang, Q.H. Qin, A fundamental solution-based finite element model for analyzing multi-layer skin burn injury, Journal of Mechanics in Medicine and Biology, 12(5) (2012) 1250027. [71] Z.W. Zhang, H. Wang, Q.H. Qin, Transient bioheat simulation of the laser-tissue interaction in human skin using hybrid finite element formulation, MCB: Molecular & Cellular Biomechanics, 9(1) (2012) 31-54. [72] Z.W. Zhang, H. Wang, Q.H. Qin, Analysis of transient bioheat transfer in the human eye using hybrid finite element model, Applied Mechanics and Materials, 553 (2014) 356-361. [73] C. Cao, Q.H. Qin, A. Yu, A novel boundary-integral based finite element method for 2D and 3D thermoelasticity problems, Journal of Thermal Stresses, 35(10) (2012) 849-876. [74] Q.H. Qin, H. Wang, Special circular hole elements for thermal analysis in cellular solids with multiple circular holes, International Journal of Computational Methods, 10(4) (2013) 1350008. [75] H. Wang, Q.H. Qin, A new special element for stress concentration analysis of a plate with elliptical holes, Acta Mechanica, 223(6) (2012) 1323-1340. [76] Q.H. Qin, H. Wang, Fundamental solution based FEM for nonlinear thermal radiation problem, in: 12th International Conference on Boundary Element and Meshless Techniques (BeTeq 2011), ed. E.L. Albuquerque, M.H. Aliabadi, EC Ltd, Eastleigh, UK, 2011, pp. 113-118. [77] C. Cao, A. Yu, Q.H. Qin, Evaluation of effective thermal conductivity of fiber-reinforced composites by boundary integral based finite element method, International Journal of Architecture, Engineering and Construction, 1(1) (2012) 14-29. [78] C. Cao, A. Yu, Q.H. Qin, A novel hybrid finite element model for modeling anisotropic composites, Finite Elements in Analysis and Design, 64 (2013) 36-47. [79] C. Cao, A. Yu, Q.H. Qin, Mesh reduction strategy: Special element for modelling anisotropic materials with defects, in: Proceedings of the 36th International Conference on Boundary Elements and Other Mesh Reduction Methods, 22 - 24 October, 2013, Dalian, China, 2013. [80] H. Wang, Q.H. Qin, Fundamental-solution-based finite element model for plane orthotropic elastic bodies, European Journal of Mechanics-A/Solids, 29(5) (2010) 801-809. [81] J. Berger, P. Martin, V. Mantič, L. Gray, Fundamental solutions for steady-state heat transfer in an exponentially graded anisotropic material, Zeitschrift für angewandte Mathematik und Physik ZAMP, 56(2) (2005) 293-303. [82] L. Gray, T. Kaplan, J. Richardson, G.H. Paulino, Green's functions and boundary integral analysis for exponentially graded materials: heat conduction, Journal of Applied Mechanics, 70(4) (2003) 543-549. [83] H. Wang, Q.H. Qin, Y. Kang, A new meshless method for steady-state heat conduction problems in anisotropic and inhomogeneous media, Archive of Applied Mechanics, 74(8) (2005) 563-579. [84] H. Wang, Q.H. Qin, Meshless approach for thermomechanical analysis of functionally graded materials, Engineering Analysis with Boundary Elements, 32(9) (2008) 704-712. [85] H. Wang, Q.H. Qin, A meshless method for generalized linear or nonlinear Poisson-type problems, Engineering Analysis with Boundary Elements, 30(6) (2006) 515-521. [86] A. Sutradhar, G.H. Paulino, The simple boundary element method for transient heat conduction in functionally graded materials, Computer Methods in Applied Mechanics and Engineering, 193(42) (2004) 4511-4539. [87] A. Sutradhar, G.H. Paulino, L. Gray, Transient heat conduction in homogeneous and non-homogeneous materials by the Laplace transform Galerkin boundary element method, Engineering Analysis with Boundary Elements, 26(2) (2002) 119-132. [88] B. Davies, B. Martin, Numerical inversion of the Laplace transform: a survey and comparison of methods, Journal of computational physics, 33(1) (1979) 1-32. [89] S. Zhu, P. Satravaha, X. Lu, Solving linear diffusion equations with the dual reciprocity method in Laplace space, Engineering Analysis with Boundary Elements, 13(1) (1994) 1-10. [90] L. Marin, D. Lesnic, The method of fundamental solutions for nonlinear functionally graded materials, International journal of solids and structures, 44(21) (2007) 6878-6890. [91] Z.F. Yuan, H.M. Yin, elastic green’s functions for a specific graded material with a quadratic variation of elasticity, Journal of Applied Mechanics, 78(2) (2011) 021021. [92] G.D. Manolis, R.P. Shaw, Green's function for the vector wave equation in a mildly heterogeneous continuum, Wave Motion, 24(1) (1996) 59-83. International Journal of Scientific Research in Science, Engineering and Technology (ijsrset.com) 301