# Thermal Analysis of 2D FGM Beam Subjected to Thermal Loading Using Meshless Weighted Least-Square Method

Authors: H. M. Zhou; X. M. Zhang; Z. Y. Wang
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The paper analyzed the thermal problem of the 2D FGM beam using meshless weighted least-square (MWLS) method. The MWLS as a meshless method is fully independent of mesh, and an approximate function was used to construct a series of linear equations to solve the unknown field variable, which avoided the troublesome task of numerical integration. The effectiveness and accuracy of the approach were illustrated by a clamped-clamped FGM beam which was subjected with interior heat source. The volume fraction of FGM beam was assumed to be given by a simple power law distribution. The effective material properties of the FGM beam were assumed to be temperature independent and calculated by Mori-Tanaka method. The results showed that a good agreement was achieved between the proposed meshless method and commercial COMSOL Multiphysics.


---

## Body

## 1. Introduction

FGMs can resist high temperatures and are proficient in reducing the thermal stress and have received more attention from the researchers [1]. Most of these researches on FGMs have been restricted to heat conduction analyses, thermal stress analyses, thermal buckling analyses, thermal vibration, and optimization problem. Various numerical techniques, such as the finite difference method (FDM) [2, 3], finite element method (FEM) [4–7], boundary element method (BEM) [8], or more recently developed meshless methods [9–23], have been developed for analyzing thermal related problems and other problems. Because of the complexity of the relevant governing equation, analytical solutions are usually difficult to obtain for those arbitrary geometry and complex boundary conditions, and the exact solutions are usually obtained based on classical plate theory, first-order shear deformation theory, high-order shear deformation theory, and so on [24, 25]. Compared with FEM, FDM, and BEM, the meshless methods are associated with a class of numerical techniques that approximate a given differential equation or a set of differential equations using global interpolations on the discrete nodes or background mesh, exhibiting the advantages of avoiding mesh generation, simple data preparation, easy postprocessing, and so forth.Liu and Gu [15] introduced meshless methods and their programming, such as the element-free Galerkin (EFG) method [12, 16], the hp-clouds method, the meshless local Petrov-Galerkin (MLPG) method [9, 10], meshless Galerkin method using radial basis functions [11, 14, 21], the least-square method [20], and meshless point interpolation method [19]. The main advantage of the MLPG method compared with regular Galerkin-based methods is that no background mesh is used to evaluate various integrals appearing in the local weak formulation of problem, but it requires a high-order quadrature rule to obtain converged results and thus needs much more computational effort in terms of CPU time than that for the FEM. Ching and Chen [9] proposed the thermoelastic analysis of a functionally graded composite considering temperature-dependent thermomechanical properties by the MLPG method. Andrew and Senthil [12] proposed a methodology for the two-dimensional simulation and optimization of material distribution of functionally graded materials for thermomechanical processes using a genetic algorithm. Katsikadelis [13] employed the meshless analog equation method to solve the 2D elastostatic problem for inhomogeneous anisotropic. A meshless algorithm of fundamental solution coupling with radial basis functions based on analog equation theory was proposed to simulate the static thermal stress distribution in 2D FGMs [14]. Mierzwiczak et al. [17] presented the singular boundary method for steady-state nonlinear heat conduction problems. Bhavani et al. [26] solved thermoelastic equilibrium equations for a functionally graded beam to obtain the axial stress distribution. Sohn and Kim [27] analyzed static and dynamic stabilities of FG panels which are subjected to thermal and aerodynamic loads. Hiroyuki [28] presented a two-dimensional (2D) higher-order deformation theory presented for the evaluation of displacements and stresses in functionally graded (FG) plates subjected to thermal and mechanical loading.Zhou et al. presented steady-state [29] and transient-state [30] heat conduction analysis of heterogeneous material using the meshless weighted least-square method. In this paper the pure meshless method (MWLS) was then extended to solve problems of thermoelastic analysis for the FGM beam with interior heat source. The volume fractions of constituent materials composing the FGM beam are assumed to be given by a simple power law distribution. Material properties of the FGM are obtained by Mori-Tanaka method. The paper is divided as follows: firstly, we give problem description and MWLS analysis about the thermal problem. Then, in order to demonstrate the efficiency and accuracy of the proposed method, numerical implementation is given in the next section. The last section includes some conclusions.


## 2. MWLS Analysis of the Thermal Problem

The solution of MWLS analysis to the thermal problem is described in this section. The shape functions in MWLS analysis is a moving least-squares approximation scheme which is originally developed for the smooth interpolation of irregularly distributed data.


### 2.1. The Moving Least-Square (MLS) Approximation Scheme

Construct the local approximate functionfh(x) of an unknown field variable function f(x) expressed as(1)fx≈fhx=pTxaxwhere pT(x) is the basis function and the quadratic basis pT(x)={1,x,y,x2,xy,y2} is used in this paper. In 2D space, the number of basis function m=6; a(x) are unknown coefficients, which are solved by minimizing a weighted discrete residual given as(2)JJ=∑I=1NωIxfhx,xI-fxI2=∑I=1NωIxpTxax-fxI2The minimum value ofJJ may be achieved through differentiating with respect to a(x)(3)∂JJ∂ajx=2∑I=1NωIx∑i=1mpixIaix-fIpjxI=0j=1,2,…,min which xI are the positions of theN nodes, fI refers to the nodal parameter of the field variable at node I, fI=f(xI). ωI(x) is the weighting function and usually a compactly supported function that is only nonzero in a small neighborhood called the “support domain” of node xI. The exponential function is used in this study. (4)ωr=exp-r2β2-exp-β21-exp-β2r≤10r>1where r=x-xI/dmI=(x-xI)2+(y-yI)2/dmI for circular support domain and β is a constant. dmI denotes the radius of the circular support domain.To obtaina(x), (3) can be rewritten in the matrix form(5)Axax=Bxfwhere Ax=∑I=1NωI(x)p(xI)pT(xI), B(x)=ω1xpx1ω2(x)p(x2)…ωN(x)p(xN), f=[fx1,fx2,…,fxN]T.d m I is chosen to make A(x) be the nonsingular matrix everywhere in the whole domain; however, the circular support domain must have enough neighborhood nodes. Through finding out thekkth nearest points of the evaluation point x, the smallest support domain radius including these points can be obtained. The value ofkk is gained by comparing some numerical examples with their analytical solutions in Zhou et al. [29, 30].Solvinga(x) from (5), coefficients a(x) can be obtained(6)ax=A-1xBxfSubstituting (6) back into (1) and removing unknown a(x), we have(7)fhx=pTxA-1xBxf=SxfSet pTxA-1xBx=S(x), where S(x) is the shape function.


### 2.2. Heat Conduction Analysis of FGM Object

The steady-state heat conduction equation and the thermal boundary conditions of FGM objects are as follows, respectively,(8)kx∇2Tx+∇kx·∇Tx+Q=0(9)Dirichletboundary:Tx=T-x∈Γ1Neumannboundary:n·kx∇Tx=q-x∈Γ2Mixedboundary:n·kx∇Tx=hT∞-Txx∈Γ3where T(x) is the temperature field on a fixed domain Ω surrounded by a closed boundary Γ=Γ1+Γ2+Γ3. The variable x denotes the physical dimensions expressed in Cartesian coordinates, x: (x,y). n is the outward surface normal. The parameters k,T∞,h and Q are the thermal conductivity, ambient temperature, the heat transfer coefficient, and the heat resource, respectively.Substituting the unknown temperature field variableT(x) of (8) with the approximate function of (7), the residuals are minimized in a least-squares manner,(10)δ∏=∫ΩδfJk∇2SJx+∇k∇SJx·k∇2SIxfI+∇k∇SIxfI+QdΩ+∫Γ1δfJSJxSIxfI-f¯dΓ+∫Γ2δfJn·k∇SJx·n·k∇SIxfI-q-dΓ+∫Γ3δfJn·k∇SJx+hSJx·n·k∇SIxfI-hT∞-SIxfIdΓWe use an alternative discrete equation to avoid integration and (10) can be rewritten as(11)δ∏=δfJ∑s=1Nks∇2SJxs+∇ks∇SJxs·ks∇2SIxsfI+∇ks∇SIxsfI+Q+δfJ∑s=1N1SJxsSIxsfI-T¯+δfJ∑s=1N2n·ks∇SJxs·n·ks∇SIxsfI-q-+δfJ∑s=1N3n·ks∇SJxs+hSJxs·n·ks∇SIxfI-hT∞-SIxsfIThe system equations of the MWLS method for solving steady-state heat conduction equations are written in a matrix form(12)KT=P(13)K=∑s=1Nks∇2SJxs+∇ks∇SJxs·ks∇2SIxs+∇ks∇SIxs+∑s=1N1SJxsSIxs+∑s=1N2n·ks∇SJxsn·ks∇SIxs+∑s=1N3n·ks∇SJxs+hSJxs·n·ks∇SIx+hSIx(14)P=-∑s=1Nks∇2SJxs+∇ks∇SJxs]Q+∑s=1N1SJxsT¯+∑s=1N2n·ks∇SJxsq-+∑s=1N3n·ks∇SJxs+hSJxshT∞where I and J refer to node indices,N1,N2, andN3 is the number of interior nodes in the boundary Γ1,Γ2, and Γ3, respectively.


### 2.3. Thermoelastic Analysis of FGM Object

Consider the 2D FGM anisotropic linear elastic body defined in the domainΩ bounded by Γ. The governing equation and boundary condition can be written in the following form disregarding the body forces (15)σij,j=0inΩ(16)stressboundarycondition:σijnj-t-i=0onΓtdisplacementboundarycondition:ui=u-ionΓuin which σij is the components of the Cauchy stress tensor. A comma followed by index j denotes the partial differentiation with respect to coordinate xj of a material point. nj is the unit outward normal to Γ. ui are the displacement components, and u-i are the prescribed displacements on Γu. t-i are the prescribed traction on Γt. Γu and Γt are the complementary parts of the boundary Γ. MWLS analysis requires a discretization of the domain Ω; (15)-(16) becomes (17)σij,jxk=0xk∈Ω,i,j=1,2;k=1,2,…,NΩ(18)σijxknj=t-ixkxk∈Γt,i,j=1,2;k=1,2,…,Nt(19)uixk=u-ixkxk∈Γu,i=1,2;k=1,2,…,Nuwhere NΩ, Nt, and Nu is the number of interior nodes in the domain Ω, in the boundary Γt and Γu, respectively.Substitutingσij,ui of (17)~(19) with the approximate function of (7), (20)∑I=1NHIxkfI=0xk∈Ω,k=1,2,…,NΩ(21)∑I=1NQIxkfI=t-kxk∈Γt,k=1,2,…,Nt(22)∑I=1NNIxkfI=u-kxk∈Γu,k=1,2,…,Nuwhere H, Q, and N denote shape function similar to S of (7)(23)HIxk=E1-υ2∂2NIxk∂x2+1-υ2∂2NIxk∂y21+υ2∂2NIxk∂x∂y1+υ2∂2NIxk∂x∂y∂2NIxk∂y2+1-υ2∂2NIxk∂x2(24)QIxk=E1-υ2l∂NIxk∂x+m1-υ2∂NIxk∂ylυ∂NIxk∂y+m1-υ2∂NIxk∂xmυ∂NIxk∂x+l1-υ2∂NIxk∂ym∂NIxk∂y+l1-υ2∂NIxk∂x(25)NI=NIxk00NIxk,fI=f1If2I,t-k=t-1xkt-2xk,u-k=u-1xku-2xkwhere l=cos(No,x), m=cos(No,y);Nois the normal vector of any point.Substituting the unknown field variableσij of (15)~(16) into (1), the residuals are minimized in a least-squares manner,(26)Π=∫Ωσij,jσik,kdΩ+∫Γuui-u-iui-u-idΓu+∫Γtσijnj-t-iσijnj-t-idΓtSimilar to (12) the system equations of the MWLS method for solving thermoelastic problem is written in the following matrix form(27)Kd=Pwhere (28)K=∑s=1NHTxsHxs+∑s=1NuNTxsNxs+∑s=1NtQTxsQxs(29)P=∑s=1NuNTxsu-s+∑s=1NtQTxst-swhere H, Q, and N are obtained by (23)~(25) and d denotes the displacement of x, y, [u1,u2].The thermal stresses is written in the matrix form(30)σ=Dε-βTwhere D is the stiffness matrix for a linearly elastic, isotropic 2-D solid.ɛ is the infinitesimal strain vector. (31)ε=∂u1∂x1∂u2∂x2∂u2∂x1+∂u1∂x2T,D=E-x1-υ-x21υ-x0υ-x10001-υ-x2in which E-=E;υ-=υ;β=αE/1-υ110, for plane stress with E,υ and α denoting the Young’s modulus, Poisson’s ratio, and coefficient of thermal expansion, respectively, and E-=E/1-υ2;υ-=υ/1-υ;β=αE/1-2υ110 for plane stain.


### 2.4. Material Properties

Two homogenization methods are often used to evaluate the effective material properties for FGMs. One is the rule of mixtures, and the other is the micromechanical model. The former is simply a linear rule of mixtures and the effective value can be determined by(32)P=P1V1+P2V2where the volume fractions satisfy V1+V2=1, and P may be elastic modulus E, bulk modulus G,Poisson’s ratio v, coefficient of thermal expansiona, thermal conductivity k, and shear modulusμ.The most widely used micromechanical model is Mori-Tanaka model [9, 31]. It is the modified rule of mixtures and the effective material properties can be defined using the following relation,(33)k=k1+3k1V2k2-k13k1+V1k2-k1(34)E=E1+V23E1+4μ1E2-E131-V2E2-E1+3E1+4μ1(35)μ=μ1+V2μ1+ff1μ2-μ11-V2μ2-μ1+μ1+ff1(36)α=α1+E2E1-Eα2-α1EE1-E2where G1=E1/3(1-2υ1);G2=E2/3(1-2υ2);μ1=E1/2(1+υ1);μ2=E2/2(1+υ2);ff1=μ1(9E1+8μ1)/6(E1+2μ1).


## 2.1. The Moving Least-Square (MLS) Approximation Scheme

Construct the local approximate functionfh(x) of an unknown field variable function f(x) expressed as(1)fx≈fhx=pTxaxwhere pT(x) is the basis function and the quadratic basis pT(x)={1,x,y,x2,xy,y2} is used in this paper. In 2D space, the number of basis function m=6; a(x) are unknown coefficients, which are solved by minimizing a weighted discrete residual given as(2)JJ=∑I=1NωIxfhx,xI-fxI2=∑I=1NωIxpTxax-fxI2The minimum value ofJJ may be achieved through differentiating with respect to a(x)(3)∂JJ∂ajx=2∑I=1NωIx∑i=1mpixIaix-fIpjxI=0j=1,2,…,min which xI are the positions of theN nodes, fI refers to the nodal parameter of the field variable at node I, fI=f(xI). ωI(x) is the weighting function and usually a compactly supported function that is only nonzero in a small neighborhood called the “support domain” of node xI. The exponential function is used in this study. (4)ωr=exp-r2β2-exp-β21-exp-β2r≤10r>1where r=x-xI/dmI=(x-xI)2+(y-yI)2/dmI for circular support domain and β is a constant. dmI denotes the radius of the circular support domain.To obtaina(x), (3) can be rewritten in the matrix form(5)Axax=Bxfwhere Ax=∑I=1NωI(x)p(xI)pT(xI), B(x)=ω1xpx1ω2(x)p(x2)…ωN(x)p(xN), f=[fx1,fx2,…,fxN]T.d m I is chosen to make A(x) be the nonsingular matrix everywhere in the whole domain; however, the circular support domain must have enough neighborhood nodes. Through finding out thekkth nearest points of the evaluation point x, the smallest support domain radius including these points can be obtained. The value ofkk is gained by comparing some numerical examples with their analytical solutions in Zhou et al. [29, 30].Solvinga(x) from (5), coefficients a(x) can be obtained(6)ax=A-1xBxfSubstituting (6) back into (1) and removing unknown a(x), we have(7)fhx=pTxA-1xBxf=SxfSet pTxA-1xBx=S(x), where S(x) is the shape function.


## 2.2. Heat Conduction Analysis of FGM Object

The steady-state heat conduction equation and the thermal boundary conditions of FGM objects are as follows, respectively,(8)kx∇2Tx+∇kx·∇Tx+Q=0(9)Dirichletboundary:Tx=T-x∈Γ1Neumannboundary:n·kx∇Tx=q-x∈Γ2Mixedboundary:n·kx∇Tx=hT∞-Txx∈Γ3where T(x) is the temperature field on a fixed domain Ω surrounded by a closed boundary Γ=Γ1+Γ2+Γ3. The variable x denotes the physical dimensions expressed in Cartesian coordinates, x: (x,y). n is the outward surface normal. The parameters k,T∞,h and Q are the thermal conductivity, ambient temperature, the heat transfer coefficient, and the heat resource, respectively.Substituting the unknown temperature field variableT(x) of (8) with the approximate function of (7), the residuals are minimized in a least-squares manner,(10)δ∏=∫ΩδfJk∇2SJx+∇k∇SJx·k∇2SIxfI+∇k∇SIxfI+QdΩ+∫Γ1δfJSJxSIxfI-f¯dΓ+∫Γ2δfJn·k∇SJx·n·k∇SIxfI-q-dΓ+∫Γ3δfJn·k∇SJx+hSJx·n·k∇SIxfI-hT∞-SIxfIdΓWe use an alternative discrete equation to avoid integration and (10) can be rewritten as(11)δ∏=δfJ∑s=1Nks∇2SJxs+∇ks∇SJxs·ks∇2SIxsfI+∇ks∇SIxsfI+Q+δfJ∑s=1N1SJxsSIxsfI-T¯+δfJ∑s=1N2n·ks∇SJxs·n·ks∇SIxsfI-q-+δfJ∑s=1N3n·ks∇SJxs+hSJxs·n·ks∇SIxfI-hT∞-SIxsfIThe system equations of the MWLS method for solving steady-state heat conduction equations are written in a matrix form(12)KT=P(13)K=∑s=1Nks∇2SJxs+∇ks∇SJxs·ks∇2SIxs+∇ks∇SIxs+∑s=1N1SJxsSIxs+∑s=1N2n·ks∇SJxsn·ks∇SIxs+∑s=1N3n·ks∇SJxs+hSJxs·n·ks∇SIx+hSIx(14)P=-∑s=1Nks∇2SJxs+∇ks∇SJxs]Q+∑s=1N1SJxsT¯+∑s=1N2n·ks∇SJxsq-+∑s=1N3n·ks∇SJxs+hSJxshT∞where I and J refer to node indices,N1,N2, andN3 is the number of interior nodes in the boundary Γ1,Γ2, and Γ3, respectively.


## 2.3. Thermoelastic Analysis of FGM Object

Consider the 2D FGM anisotropic linear elastic body defined in the domainΩ bounded by Γ. The governing equation and boundary condition can be written in the following form disregarding the body forces (15)σij,j=0inΩ(16)stressboundarycondition:σijnj-t-i=0onΓtdisplacementboundarycondition:ui=u-ionΓuin which σij is the components of the Cauchy stress tensor. A comma followed by index j denotes the partial differentiation with respect to coordinate xj of a material point. nj is the unit outward normal to Γ. ui are the displacement components, and u-i are the prescribed displacements on Γu. t-i are the prescribed traction on Γt. Γu and Γt are the complementary parts of the boundary Γ. MWLS analysis requires a discretization of the domain Ω; (15)-(16) becomes (17)σij,jxk=0xk∈Ω,i,j=1,2;k=1,2,…,NΩ(18)σijxknj=t-ixkxk∈Γt,i,j=1,2;k=1,2,…,Nt(19)uixk=u-ixkxk∈Γu,i=1,2;k=1,2,…,Nuwhere NΩ, Nt, and Nu is the number of interior nodes in the domain Ω, in the boundary Γt and Γu, respectively.Substitutingσij,ui of (17)~(19) with the approximate function of (7), (20)∑I=1NHIxkfI=0xk∈Ω,k=1,2,…,NΩ(21)∑I=1NQIxkfI=t-kxk∈Γt,k=1,2,…,Nt(22)∑I=1NNIxkfI=u-kxk∈Γu,k=1,2,…,Nuwhere H, Q, and N denote shape function similar to S of (7)(23)HIxk=E1-υ2∂2NIxk∂x2+1-υ2∂2NIxk∂y21+υ2∂2NIxk∂x∂y1+υ2∂2NIxk∂x∂y∂2NIxk∂y2+1-υ2∂2NIxk∂x2(24)QIxk=E1-υ2l∂NIxk∂x+m1-υ2∂NIxk∂ylυ∂NIxk∂y+m1-υ2∂NIxk∂xmυ∂NIxk∂x+l1-υ2∂NIxk∂ym∂NIxk∂y+l1-υ2∂NIxk∂x(25)NI=NIxk00NIxk,fI=f1If2I,t-k=t-1xkt-2xk,u-k=u-1xku-2xkwhere l=cos(No,x), m=cos(No,y);Nois the normal vector of any point.Substituting the unknown field variableσij of (15)~(16) into (1), the residuals are minimized in a least-squares manner,(26)Π=∫Ωσij,jσik,kdΩ+∫Γuui-u-iui-u-idΓu+∫Γtσijnj-t-iσijnj-t-idΓtSimilar to (12) the system equations of the MWLS method for solving thermoelastic problem is written in the following matrix form(27)Kd=Pwhere (28)K=∑s=1NHTxsHxs+∑s=1NuNTxsNxs+∑s=1NtQTxsQxs(29)P=∑s=1NuNTxsu-s+∑s=1NtQTxst-swhere H, Q, and N are obtained by (23)~(25) and d denotes the displacement of x, y, [u1,u2].The thermal stresses is written in the matrix form(30)σ=Dε-βTwhere D is the stiffness matrix for a linearly elastic, isotropic 2-D solid.ɛ is the infinitesimal strain vector. (31)ε=∂u1∂x1∂u2∂x2∂u2∂x1+∂u1∂x2T,D=E-x1-υ-x21υ-x0υ-x10001-υ-x2in which E-=E;υ-=υ;β=αE/1-υ110, for plane stress with E,υ and α denoting the Young’s modulus, Poisson’s ratio, and coefficient of thermal expansion, respectively, and E-=E/1-υ2;υ-=υ/1-υ;β=αE/1-2υ110 for plane stain.


## 2.4. Material Properties

Two homogenization methods are often used to evaluate the effective material properties for FGMs. One is the rule of mixtures, and the other is the micromechanical model. The former is simply a linear rule of mixtures and the effective value can be determined by(32)P=P1V1+P2V2where the volume fractions satisfy V1+V2=1, and P may be elastic modulus E, bulk modulus G,Poisson’s ratio v, coefficient of thermal expansiona, thermal conductivity k, and shear modulusμ.The most widely used micromechanical model is Mori-Tanaka model [9, 31]. It is the modified rule of mixtures and the effective material properties can be defined using the following relation,(33)k=k1+3k1V2k2-k13k1+V1k2-k1(34)E=E1+V23E1+4μ1E2-E131-V2E2-E1+3E1+4μ1(35)μ=μ1+V2μ1+ff1μ2-μ11-V2μ2-μ1+μ1+ff1(36)α=α1+E2E1-Eα2-α1EE1-E2where G1=E1/3(1-2υ1);G2=E2/3(1-2υ2);μ1=E1/2(1+υ1);μ2=E2/2(1+υ2);ff1=μ1(9E1+8μ1)/6(E1+2μ1).


## 3. Numerical Results and Discussions

In order to demonstrate the efficiency and accuracy of the presented method, firstly, we choose a isotropic square region (Case 1) with defined boundary conditions; through heat conduction analysis the results are compared with the analytical solutions and FDM. Then a clamped-clamped FGM beam (Case 2) which was subjected with interior heat source is analyzed using MWLS method.


### 3.1. Case 1

A 100×100m isotropic square region is shown in Figure1. The top and bottom boundaries are insulated. The left and right boundaries are assigned a temperature of 200°C and 100°C, respectively. The spatial variation of the thermal conductivity is taken to be cubic in the x-direction as k(x,y)=(1+x/100)3.Figure 1
Geometry of square region and boundary conditions.An analytical solution is given as(37)Tx,y=800611+x/1002+12The design region is discretized as 20×20,β=4 in (4), and the material properties are a linear rule of mixtures (Eq. (32)). The temperature field distribution is computed in Eq. (12)-(14) and shown in Figure 2. According to Wang et al. [32] the estimation rule is maximum of relative error% = maxf-g/f×100%, in which f is the analytical solution and g is numerical solution. Then, our MWLS method is compared with FDM in different resolutions. The results are shown in Table 1. Obviously our method has a high precision compared with the analytical solution in spite of any resolution.Table 1
Maximum of relative error.
Cartesian grid MWLS FDM 10×10 0.15% 0.35% 20×20 0.05% 0.077% 30×30 0.024% 0.033%Figure 2
Temperature field distribution (20×20).
(a)
2D display (b)
Top boundary


### 3.2. Case 2

A clamped-clamped FGM beam is shown in Figure3, length L=1m, width D=0.5m, and thickness H=0.1m; material property is shown in Table 2, interior heat source Q=5e5W/m3; the spatial variation of the volume fraction of Al is taken to be a power law distribution in the y-direction as V=y/Da.Table 2
Material property.
Material property Al SiC k(W/mK) 233 65 E(Pa) 7e10 4.27e11 υ 0.3 0.17 α (/K) 2.34e-5 4.3e-6Figure 3
Initial condition.The beam is assumed to be in a state of plane strain normal to the xy plane, and the design region is discretized as 31×15 and 61×30. The effective material properties are determined by the Mori-Tanaka model (Eq. (33)~(36)). In order to verify the proposed computational method, we do some comparisons between the MWLS and the commercial COMSOL Multiphysics for a homogeneous material (a=0), relevant results are shown in Figure 4 and listed in Table 3. The results obtained with the two methods are in good agreement in temperature field aspect; however, the x-displacements and y-displacement have a little difference in 31×15 grid and in 61×30 the results are in good agreement in Table 3. From Figure 4 we also can know that our method and COMSOL Multiphysics have the same distribution trend. The maximum temperature 361.3K is in the center (0.5,0.25) of the beam.Table 3
Comparison of MWLS method and COMSOL Multiphysics.
method Temperature/K X displacement /mm Y displacement/mm maximum minimum maximum minimum maximum minimum MWLS (31×15) 361.8 300 .596 -.873 2.95 -2.8 MWLS (61×30) 361.3 300 .71 -.72 3.88 -3.84 COMSOL (948 triangle mesh) 361.1 300 .639 -.639 3.81 -3.81Figure 4
Comparison of our method with COMOSOL of a=0.
(a)
Temperature field of our method (b)
Temperature field of COMSOL (c)
x-displacement of our method (d)
x-displacement of COMSOL (e)
y-displacement of our method (f)
y-displacement of COMSOLFor a=2, we analyzed the heat conduction and thermoelastic problem using MWLS method. Temperature field distribution, x-displacement, and y-displacement are plotted in Figures5, 6, and 7, respectively. Figure 5 indicates that the maximum temperature 425.6K is higher than the homogeneous material of Figure 2(a), in the Cartesian coordinates (0.5,0.179) of the beam. Figures 6 and 7 show that when subjected to temperature rise, the beam expands and the maximum y-displacement is located at the top middle of the beam. Then we do heat conduction analysis in different material distribution in different heat source; the result is listed in Table 4. From Table 4, we can know that the maximum temperature of FGM model is higher than that of the fully metal model. Moreover, as the volume fraction index is increased, the maximum temperature increases. This is because for FGMs, when the volume fraction index is increased, the contained quantity of ceramic increases. Finally, to make a comparison, we do thermoelastic analysis and obtain thermal stresses in the neutral axis of the beam among a=0, a=2, and a=3, as shown in Figure 8. In Figures 8(a) and 8(c), the volume fraction of Al is gradually decreased from a=0 to a=3, the σx and σy stresses are in an upward trend. The maximum thermal stress always occurred in the vicinity of neutral axis of the beam from Figures 8(a), 8(c), and 8(d). The results also agreed well with the presented elasticity solutions of [26].Table 4
Comparison of different material distribution in different heat source.
MWLS (31×15) Temperature/K Temperature/K Q=5e5W/m3 Q=5e6W/m3 maximum minimum maximum minimum a=0 361.8 300 910.8 300 a=2 425.6 300 1558.2 300 a=3 436.1 300 1660.9 300Figure 5
Distribution of temperature field.Figure 6
x-displacement/m.Figure 7
y-displacement/m.Figure 8
Stress results comparison of thermoelastic analysis among a=0, a=2, and a=3.
(a)
σ x ( x , D / 2 ) (b)
σ x ( L / 2 , y ) (c)
σ y ( x , D / 2 ) (d)
σ y ( L / 2 , y )


## 3.1. Case 1

A 100×100m isotropic square region is shown in Figure1. The top and bottom boundaries are insulated. The left and right boundaries are assigned a temperature of 200°C and 100°C, respectively. The spatial variation of the thermal conductivity is taken to be cubic in the x-direction as k(x,y)=(1+x/100)3.Figure 1
Geometry of square region and boundary conditions.An analytical solution is given as(37)Tx,y=800611+x/1002+12The design region is discretized as 20×20,β=4 in (4), and the material properties are a linear rule of mixtures (Eq. (32)). The temperature field distribution is computed in Eq. (12)-(14) and shown in Figure 2. According to Wang et al. [32] the estimation rule is maximum of relative error% = maxf-g/f×100%, in which f is the analytical solution and g is numerical solution. Then, our MWLS method is compared with FDM in different resolutions. The results are shown in Table 1. Obviously our method has a high precision compared with the analytical solution in spite of any resolution.Table 1
Maximum of relative error.
Cartesian grid MWLS FDM 10×10 0.15% 0.35% 20×20 0.05% 0.077% 30×30 0.024% 0.033%Figure 2
Temperature field distribution (20×20).
(a)
2D display (b)
Top boundary


## 3.2. Case 2

A clamped-clamped FGM beam is shown in Figure3, length L=1m, width D=0.5m, and thickness H=0.1m; material property is shown in Table 2, interior heat source Q=5e5W/m3; the spatial variation of the volume fraction of Al is taken to be a power law distribution in the y-direction as V=y/Da.Table 2
Material property.
Material property Al SiC k(W/mK) 233 65 E(Pa) 7e10 4.27e11 υ 0.3 0.17 α (/K) 2.34e-5 4.3e-6Figure 3
Initial condition.The beam is assumed to be in a state of plane strain normal to the xy plane, and the design region is discretized as 31×15 and 61×30. The effective material properties are determined by the Mori-Tanaka model (Eq. (33)~(36)). In order to verify the proposed computational method, we do some comparisons between the MWLS and the commercial COMSOL Multiphysics for a homogeneous material (a=0), relevant results are shown in Figure 4 and listed in Table 3. The results obtained with the two methods are in good agreement in temperature field aspect; however, the x-displacements and y-displacement have a little difference in 31×15 grid and in 61×30 the results are in good agreement in Table 3. From Figure 4 we also can know that our method and COMSOL Multiphysics have the same distribution trend. The maximum temperature 361.3K is in the center (0.5,0.25) of the beam.Table 3
Comparison of MWLS method and COMSOL Multiphysics.
method Temperature/K X displacement /mm Y displacement/mm maximum minimum maximum minimum maximum minimum MWLS (31×15) 361.8 300 .596 -.873 2.95 -2.8 MWLS (61×30) 361.3 300 .71 -.72 3.88 -3.84 COMSOL (948 triangle mesh) 361.1 300 .639 -.639 3.81 -3.81Figure 4
Comparison of our method with COMOSOL of a=0.
(a)
Temperature field of our method (b)
Temperature field of COMSOL (c)
x-displacement of our method (d)
x-displacement of COMSOL (e)
y-displacement of our method (f)
y-displacement of COMSOLFor a=2, we analyzed the heat conduction and thermoelastic problem using MWLS method. Temperature field distribution, x-displacement, and y-displacement are plotted in Figures5, 6, and 7, respectively. Figure 5 indicates that the maximum temperature 425.6K is higher than the homogeneous material of Figure 2(a), in the Cartesian coordinates (0.5,0.179) of the beam. Figures 6 and 7 show that when subjected to temperature rise, the beam expands and the maximum y-displacement is located at the top middle of the beam. Then we do heat conduction analysis in different material distribution in different heat source; the result is listed in Table 4. From Table 4, we can know that the maximum temperature of FGM model is higher than that of the fully metal model. Moreover, as the volume fraction index is increased, the maximum temperature increases. This is because for FGMs, when the volume fraction index is increased, the contained quantity of ceramic increases. Finally, to make a comparison, we do thermoelastic analysis and obtain thermal stresses in the neutral axis of the beam among a=0, a=2, and a=3, as shown in Figure 8. In Figures 8(a) and 8(c), the volume fraction of Al is gradually decreased from a=0 to a=3, the σx and σy stresses are in an upward trend. The maximum thermal stress always occurred in the vicinity of neutral axis of the beam from Figures 8(a), 8(c), and 8(d). The results also agreed well with the presented elasticity solutions of [26].Table 4
Comparison of different material distribution in different heat source.
MWLS (31×15) Temperature/K Temperature/K Q=5e5W/m3 Q=5e6W/m3 maximum minimum maximum minimum a=0 361.8 300 910.8 300 a=2 425.6 300 1558.2 300 a=3 436.1 300 1660.9 300Figure 5
Distribution of temperature field.Figure 6
x-displacement/m.Figure 7
y-displacement/m.Figure 8
Stress results comparison of thermoelastic analysis among a=0, a=2, and a=3.
(a)
σ x ( x , D / 2 ) (b)
σ x ( L / 2 , y ) (c)
σ y ( x , D / 2 ) (d)
σ y ( L / 2 , y )


## 4. Conclusion

In this paper, a novel thermoelastic analysis of FGM beam based on MWLS method was presented. We do thermoelastic and heat conduction analysis aimed at a clamped-clamped thick beam which is subjected with interior heat source. The FGM beam is assumed to be given by a simple power law distribution. Material properties of the FGM beam are obtained by Mori-Tanaka method. Through being compared with analytical solution and the commercial software of COMSOL Multiphysics, the effectiveness and accuracy are verified. We also listed the comparison of thermal stresses with the variation of power law index. The present method of analysis will be also useful in the design and optimization of FGM objects.



---

*Source: 2541707-2019-04-30.xml*