- Email: [email protected]
Transient heat conduction analysis in functionally graded materials by the meshless local boundary integral equation method
Computational Materials Science 28 (2003) 494–504 www.elsevier.com/locate/commatsci
Transient heat conduction analysis in functionally graded materials by the meshless local boundary integral equation method J. Sladek b
a,*
, V. Sladek a, Ch. Zhang
b
a Institute of Construction and Architecture, Slovak Academy of Sciences, 84220 Bratislava, Slovakia Department of Civil Engineering, University of Applied Sciences Zittau/G€orlitz, D-02763 Zittau, Germany
Abstract Advanced computational method for transient heat conduction analysis in continuously nonhomogeneous functionally graded materials (FGM) is proposed. The method is based on the local boundary integral equations with moving least square approximation of the temperature and heat flux. The initial-boundary value problem is solved by the Laplace transform technique. Both Papoulis and Stehfest algorithms are applied for the numerical Laplace inversion to obtain the time-dependent solutions. Numerical results are presented for a finite strip and a hollow cylinder with an exponential spatial variation of material parameters. Ó 2003 Elsevier B.V. All rights reserved. Keywords: Transient heat conduction analysis; Local boundary integral equations; Functionally graded materials (FGM); Laplace transform; MLS approximation
1. Introduction Due to their excellent thermal properties at high temperatures and their superior corrosion and wear resistances, ceramic materials are preferred and favored in many engineering structures and components. One major disadvantage of ceramics is their brittleness and low fracture toughness. To eliminate this disadvantage and to increase their toughness, ceramic materials are frequently reinforced by metallic particles or fibres, which yields
*
Corresponding author. Tel.: +421-2-54788662; fax: +421-254773548. E-mail address: [email protected] (J. Sladek).
new composite materials. To fabricate a structural component with a constant mixture ratio of ceramics and metals is not an efficient way due to specific requirements in material processing on individual parts of the material surfaces. A more sophisticated way is to have volume fractions of composite constituents varying continuously in space. This class of composite materials is referred to as functionally graded materials (FGMs) [1,2]. Since functionally graded materials are frequently used for structures under thermal load, it is indispensable to analyze their thermal properties. Due to the complexity of the problem, only a limited number of investigations dealing with heat conduction analysis in FGMs can be found in literature. Most of them used an exponential
0927-0256/$ - see front matter Ó 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2003.08.006
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
variation of material constants with Cartesian coordinates under stationary thermal loading conditions [3–5]. Transient heat transfer in FGMs with an exponential spatial variation of material constants has been examined by Jin and Batra [6], Noda and Jin [7] and Jin and Paulino [8]. Recently, Sutradhar et al. [9] presented a Laplace transform Galerkin BEM for 3-D transient heat conduction analysis by using the GreenÕs function approach where an exponential law for the FGMs was used. Since analytical methods are usually restricted to simple geometry and boundary conditions, solution of the boundary value problems for nonhomogeneous bodies requires special numerical methods due to the high mathematical complexity. Beside the well established finite element method (FEM), the boundary element method (BEM) provides an efficient and popular alternative to the FEM for solving certain class of boundary value problems. However, the pure BEM formulation can be applied only to problems where fundamental solution is available. For a general nonhomogeneous body the fundamental solution is unfortunately not available. One possibility to obtain a BEM formulation is based on the use of fundamental solutions for a fictitious homogeneous medium, as has been suggested for the first time by Butterfield [10] for potential flow problems. This approach, however, leads to a boundary-domain integral formulation with additional domain integrals for the gradients of primary fields to obtain a unique formulation. This approach has been applied to the heat conduction analysis in anisotropic nonhomogeneous media by Tanaka and Tanaka [11]. In this approach, some computational difficulties in the numerical implementation have been reported. Alternatively, the analyzed domain can be divided into small subdomains with a simple geometry such as circle and sphere. To each subdomain one can relate the fundamental solution to the corresponding material constants [12,13]. On the surface of subdomains the local boundary integral equations (LBIE) are applied. The temperature and heat flux on the local boundary and in the interior of the subdomain are approximated by the moving least square (MLS). The key idea of MLS interpolants is that it is only necessary to construct an array of nodes in the domain under consideration [14]. Thus, the method is completely element-free.
495
In this paper, a local boundary element method for transient heat conduction analysis is presented. The transient heat conduction problem is described by a parabolic partial differential equation. To eliminate the time dependence of the differential equation, the Laplace transform technique is used. A pure boundary integral formulation for the global boundary requires the application of the fundamental solution in the Laplace transform domain, which is not known for general FGM. In the present work, a simpler fundamental solution corresponding to the PoissonÕs equation is adopted which leads to a boundary-domain integral formulation. Since the boundary-domain formulation in the global BEM is numerically not efficient, we apply these integral equations to small subdomains, which result in a local boundary integral equation formulation. By choosing a simple geometric shape such as a circle for the subdomains, the arising domain integrals in the local BEM can be easily evaluated. It means that the boundarydomain integral formulation for simple subdomains does not bring any difficulties. Several quasi-static boundary value problems are solved for various values of the Laplace transform parameter. The Papoulis [15] and Stehfest [16] numerical inversion methods are applied to obtain the time-dependent solutions. The exponential spatial variation of thermal material properties is considered in numerical examples. Results are given for a finite strip and a hollow cylinder. Both stationary and transient loading (thermal shock) conditions are considered.
2. Local boundary integral equations Let us consider a continuously nonhomogeneous FGM with zero initial condition hðx; tÞ ¼ 0 at t ¼ 0, which is subjected to a thermal loading on its boundary. The thermal heat conduction problem in the FGM is described by the following governing equation h;ii ðx; tÞ ¼
1 oh k;i ðx; tÞ þ ðxÞh;i ðx; tÞ jðxÞ ot k
1 Qðx; tÞ; jðxÞ
ð1Þ
496
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
where hðx; tÞ is the temperature field, k and j stand for the thermal conductivity and diffusivity, and Qðx; tÞ is the density of the body heat sources, respectively. Applying the Laplace transform Z 1 L½hðx; tÞ ¼ hðx; pÞ ¼ hðx; tÞept dt 0
to the governing equation (1) we obtain p k;i h;ii ðx; pÞ hðx; pÞ þ ðxÞh;i ðx; pÞ j k 1 ¼ F ðx; pÞ; j
ð2Þ
where p is the Laplace transform parameter and F ðx; pÞ ¼ Qðx; pÞ þ hðx; 0Þ is the redefined body heat source in the Laplace transform domain with the initial-boundary condition for the temperature field. The solution of the governing equation (2) can be found in a weak form with a fundamental solution h ðx; yÞ ¼
1 1 ln 2p r
ð3Þ hðy; pÞ ¼
corresponding to the Poisson equation h;ii ðx; yÞ
þ dðx; yÞ ¼ 0;
where C is the boundary of the analysed domain X. Since the boundary-domain integral representation (5) for the Laplace transform of the temperature contains temperature gradients, the integral representation for the temperature gradients at interior points has to be added to obtain a unique integral formulation. It should be noted that the singularity of the kernel in the integral representation for the temperature gradients is increased to a hypersingularity. By applying the representation formulae for the temperature and its gradients to the global boundary of the considered domain, a unique global BEM formulation is obtained. Domain discretization and evaluation of hypersingular integrals make the global BEM cumbersome. Another approach is the use of the local boundary integral equations, related to the boundaries of small regular no-overlapping subdomains randomly spread in the analyzed domain. If, instead of the entire domain X, we consider a subdomain Xs , the following local boundary integral equation should also hold over the subdomain Z oXs
ð4Þ
where dðx; yÞ is the Dirac delta function and r is the distance of the field and source points, i.e., r ¼ jx yj. By the use of the Gauss divergence theorem to the weak form of Eq. (2) we obtain an integral representation for the temperature field in the Laplace transform domain Z oh hðy; pÞ ¼ ðx; pÞh ðx; yÞ dC C on Z oh ðx; yÞ dC hðx; pÞ on ZC p hðx; pÞh ðx; yÞ dX jðxÞ ZX k;i ðxÞh;i ðx; pÞh ðx; yÞ dX þ X k Z 1 F ðx; pÞh ðx; yÞ dX; ð5Þ þ X jðxÞ
þ þ
Z
Z Z Z
oh ðx; pÞh ðx; yÞ dC on
oh hðx; pÞ ðx; yÞ dC on oXs Xs
p hðx; pÞh ðx; yÞ dX jðxÞ
Xs
k;i ðxÞh;i ðx; pÞh ðx; yÞ dX k
Xs
1 F ðx; pÞh ðx; yÞ dX: jðxÞ
ð6Þ
The integral equation (6) is considered for small subdomains Xs X. Hence, no of the boundary densities in the Laplace transformed domain are required on the local boundary oXs as long as it lies entirely inside Xs . To reduce the number of unknowns on oXs the concept of a companion solution can be utilized successfully [12]. The companion solution is associated with the fundamental solution h ðx; yÞ and is defined as the solution to the following equations:
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
h~;ii ðx; yÞ ¼ 0 on X0s ; h~;ii ðx; yÞ ¼ h ðx; yÞ on oX0s ;
local boundary ∂Ωs=∂Ω's=Ls subdomain Ωs=Ω's
ð7Þ
where X0s and oX0s are defined in Fig. 1. For simplicity, X0s is taken as a circle in the present implementation. The modified test function
ri
Ωx x
1 r0 ln h~ ðx; yÞ ¼ h ðx; yÞ h~ðx; yÞ ¼ 2p r
Ls
node xi
has to satisfy the governing equation (4). It is the GreenÕs function vanishing on the boundary of the circular subdomain of radius r0 . Since the integral equation (6) is also valid for the modified fundamental solution h~ ðx; yÞ, we can obtain Z oh~ hðy; pÞ ¼ hðx; pÞ ðx; yÞ dC on oX Z s p hðx; pÞh~ ðx; yÞ dX jðxÞ Xs Z k;i ðxÞh;i ðx; pÞh~ ðx; yÞ dX þ k X Z s 1 þ F ðx; pÞh~ ðx; yÞ dX ð8Þ jðxÞ Xs for the interior source point y 2 X and Z oh~ hðf; pÞ ¼ hðx; pÞ ðx; fÞ dC on Ls Z oh~ lim hðx; pÞ ðx; yÞ dC y!f C on s Z oh ðx; pÞh~ ðx; fÞ dC þ Cs on Z p hðx; pÞh~ ðx; fÞ dX jðxÞ Xs Z k;i ðxÞh;i ðx; pÞh~ ðx; fÞ dX þ k X Z s 1 þ F ðx; pÞh~ ðx; fÞ dX jðxÞ Xs
497
Γs support of node xi
for the source point located on the global boundary f 2 Cs C. Note that oXs ¼ Ls [ Cs , where Ls is the circular part of oXs . The integral equation (9) is written in the limit form. Such an expression is appropriate [17] if the physical fields such as the temperature are given in a nonpolynomial approximation form.
Ωs'
∂Ωs'
Fig. 1. Local boundaries and the domain of definition of the MLS approximation for the trial function.
3. The MLS approximation scheme In general, a meshless method uses a local interpolation to represent the trial function with the values (or the fictitious values) of the unknown variable at some randomly located nodes. The moving least squares (MLS) approximation used in the present analysis may be considered as one of such schemes. Consider a subdomain Xx , the neighbourhood of a point x and denoted as the domain of definition of the MLS approximation for the trial function at x, which is located in the analyzed domain X. To approximate the distribution of the field hðxÞ in Xs over a number of randomly located nodes fxi g; i ¼ 1; 2; . . . ; n, the moving least square approximant hh ðxÞ of h for 8x 2 Xs is defined as hh ðxÞ ¼ pT ðxÞaðxÞ;
ð9Þ
∂Ωs=Ls ∪ Γs
8x 2 Xs ;
ð10Þ
where pT ðxÞ ¼ ½p1 ðxÞ; p2 ðxÞ; . . . ; pm ðxÞ is a complete monomial basis of order m and aðxÞ is a vector with components aa ðxÞ; a ¼ 1; 2; . . . ; m. The coefficient vector aðxÞ is determined by minimizing a weight discrete L2 norm which is defined as n X 2 wi ðxÞ½pT ðxi ÞaðxÞ h^i ; ð11Þ J ðxÞ ¼ i¼1
where wi ðxÞ is the weight function associated with the node i. Recall that n is the number of nodes in Xx for which the weight function wi ðxÞ > 0 and h^i are the fictitious nodal values, but not the nodal
498
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
values of the unknown trial function hh ðxÞ in general [18]. The stationarity of J ðxÞ in Eq. (11) with respect to aðxÞ leads to the following linear relation between aðxÞ and ^ h AðxÞaðxÞ ¼ BðxÞ^ h;
ð12Þ
where AðxÞ ¼
n X wi ðxÞpðxi ÞpT ðxi Þ; i¼1
) oh~ i b ðx; y Þ/ ðxÞ dC h^i ðpÞ / ðy Þ þ on oXs i¼1 n Z X p b ~ /i ðxÞ h ðx; y Þ ¼ jðxÞ Xs i¼1 k;j i ðxÞ/;j ðxÞ dX hi ðpÞ k Z 1 þ F ðx; pÞh~ ðx; yb Þ dX for yb 2 X Xs jðxÞ
n X
(
i
b
Z
BðxÞ ¼ ½w1 ðxÞpðx1 Þ; w2 ðxÞpðx2 Þ; . . . ; wn ðxÞpðxn Þ: Solving for aðxÞ from Eq. (12) and substituting it into Eq. (10) give a relation n X h¼ /i ðxÞh^i ; ð13Þ hh ðxÞ ¼ UT ðxÞ ^
ð17Þ and n X
(
Z Ls
UT ðxÞ ¼ pT ðxÞA1 ðxÞBðxÞ:
where d i ¼ kx xi k, ci is a constant controlling the shape of the weight function wi ðxÞ and ri is the size of the support domain. The size of the support should be large enough in the domain of definition to ensure the regularity of the matrix A. Finally, one can write the approximation formula for the Laplace transform of the temperature and its normal derivative as n X hðx; pÞ ¼ /i ðxÞh^i ðpÞ; ð15Þ i¼1
Z
Z Cs
oh~ ðx; yÞ/i ðxÞ dC on
oh~ ðx; fb Þ/i ðxÞ dC on ) nk /i;k ðxÞh~ ðx
f Þ dC h^i ðpÞ b
Cs
p b ~ /i ðxÞ ¼ h ðx f Þ jðxÞ Xs i¼1 k;j ðxÞ/i;j ðxÞ dX h^i ðpÞ k Z 1 F ðx; pÞh~ ðx; fb Þ dX for fb 2 Csq ; þ jðxÞ Xs n Z X
ð14Þ
n X oh ðx; pÞ ¼ nk ðxÞ h^i ðpÞ/i;k ðxÞ: on i¼1
/ ðf Þ þ lim
y!fb
þ
The Gaussian weight function is applied in the present work as: ( exp½ðd i =ci Þ2 exp½ðri =ci Þ2 ; 0 6 d i 6 ri i 1exp½ðri =ci Þ2 w ðxÞ ¼ 0; d i P ri ;
b
i¼1
i¼1
where
i
ð18Þ Csq is a part of the global boundary C, over which the heat flux is prescribed. If the temperature is prescribed on a part Csh we propose to collocate the approximation formula (15) by using n X
/i ðfb Þh^i ðpÞ ¼ ~hðfb ; pÞ for fb 2 Ch ;
ð19Þ
i¼1
ð16Þ
Eqs. (15) and (16) imply that both quantities h and oh=on can be approximated by the fictitious nodal values h^i ðpÞ. Substituting the approximation formulae (15) and (16) into the local boundary integral equations (8) and (9) we obtain a set of linear algebraic equations:
where ~hðfb ; pÞ is the Laplace transform of the prescribed temperature. 4. Numerical inversion of Laplace transform The time-dependent values of the transformed variables in the previous sections can be obtained
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
by an inverse transform. There are many inversion methods available for Laplace transform. As the Laplace transform inversion is an ill-posed problem, small truncation errors can be greatly magnified in the inversion process and lead to poor numerical results. In the present analysis, two numerical inversion methods are used. The first one is the most familiar Papoulis inversion method [15] and the second one is the sophisticated StehfestÕs algorithm [16], which are explained briefly in the following. If f ðpÞ is the Laplace transform of f ðtÞ, an approximate value fa of f ðtÞ in the Papoulis method is given by fa ðsÞ ¼
L X
Cm sinð2m 1Þs;
ð20Þ
m¼1
where 1 t ¼ ln coss b with a transformation of the original time interval tð0; 1Þ into the transformed interval sð0; p=2Þ. The coefficients Cm in Eq. (20) are obtained from the system of algebraic equations L X
Bnk Ck ¼ An ;
ð21Þ
499
where vi ¼ ð1ÞN =2þi
maxði;N X=2Þ
k N =2 ð2kÞ! : ðN =2 kÞ!k!ðk 1Þ!ði kÞ!ð2k iÞ! k¼½ðiþ1Þ=2 ð23Þ
Sutradhar et al. [9] have suggested to use N ¼ 10 for single precision arithmetic. It means that for each time t it is needed to solve N boundary value problems for the corresponding Laplace parameters, p ¼ i ln 2=t, with i ¼ 1; 2; . . . ; N . If M denotes the number of the time instants in which we are interested to know f ðtÞ, the number of the Laplace transform solutions f ðpj Þ is then M N , which can be much higher than that in the Papoulis method, where for selected L-values of the Laplace transform parameter we obtain fa in the entire time interval. Most of the methods for the numerical inversion of the Laplace transform require the use of complex valued Laplace transform parameter, and as a result, the application of complex arithmetic may lead to additional storage requirement and an increase in computing time. A critical study of various inversion algorithms is given by Davies and Martin [21].
k¼1
where
5. Numerical examples
4 An ¼ 22ðn1Þ bf fð2n 1Þbg; p
Though the proposed local boundary integral equation method has no restrictions on the spatial variation of the material parameters of the FGM, the numerical examples presented here use an exponential variation of the material properties with Cartesian coordinates. In such a special case the equations for the heat conduction problems are substantially simplified in comparison with a general nonhomogeneous case. For comparison purposes it is desirable to have reliable accurate results obtained by another independent method. Analytical solutions are of course the best choice if they are available. First, we consider a finite strip with a unidirectional variation of the thermal conductivity and diffusivity. The same exponential spatial variation for both parameters is taken
and the coefficients of triangular matrix Bnk can be found in Papoulis [15] and Balas et al. [19]. The number of terms L in the expansion series (20) is selected in an interval ð12; 16Þ to obtain a high accuracy of the Laplace inversion. The parameter b is determined by the maximum value of the time interval and the practical limit of smax used is smax ffi 70 deg [19]. The StehfestÕs algorithm originates from Gaver [20]. Accordingly, an approximate value fa of the inverse f ðtÞ for a specific time t is given by N ln 2 X ln 2 i ; fa ðtÞ ¼ vi f t i¼1 t
ð22Þ
500
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
kðxÞ ¼ k0 ecx1 ; jðxÞ ¼ j0 ecx1
ð24Þ
with j0 ¼ 0:17 104 m2 s1 and k0 ¼ 17 W m1 deg1 . Two different exponential parameters c ¼ 0:2 and 0.5 cm1 are selected in numerical calculation. On both opposite sides parallel to the x2 -axis two different temperatures are prescribed. One side is kept to zero temperature and the other has the Heaviside step time variation hðx2 ; tÞ ¼ T H ðtÞ with T ¼ 1 deg. On the lateral sides of the strip the heat flux vanishes. In numerical calculations, a square with a side-length a ¼ 0:04 m and a regular node distribution with equal 16 boundary and interior nodes are used for the MLS approximation (see Fig. 2). The special case with an exponential parameter c ¼ 0 corresponds to a homogeneous material. In such a case an analytical solution is available 1 x1 2 X T cos np npx1 hðx1 ; tÞ ¼ T þ sin n a p n a jn2 p2 t exp ; ð25Þ a2 which can be used to check the accuracy of the present numerical method. Numerical results are computed at three locations along the x1 -axis with x1 =a ¼ 0:25, 0.5 and 0.75. Both Papoulis and StehfestÕs numerical inversion methods of the Laplace transform are applied. In the Papoulis method 15 different Laplace transform parameters and b ¼ 0:01 are used. In the StehfestÕs method
x2 q=0
9
13
θ=0
θ =ΤΗ(t)
x1
1
q=0
5
Fig. 2. Boundary conditions and node distribution for a finite strip.
N ¼ 10 is chosen. For the present case, both methods give nearly identical results. Numerical results for the temperature distribution are presented in Fig. 3. An excellent agreement between numerical and analytical results is obtained. Next, we continue to analyze the heat conduction problem in the FGM with c ¼ 0:2 and 0.5 cm1 as mentioned above. The time variation of the temperature at the position x1 =a ¼ 0:25 for three c-values is presented in Fig. 4. Here again, both inversion methods give almost identical results. With increasing c-values the thermal conductivity increases and a higher level of temperature at the considered position is obtained in the large time range corresponding to a steady state. For steady state an analytical solution can be obtained as ecx1 1 hðx1 Þ ¼ T ca : ð26Þ e 1 Analytical and numerical results corresponding to stationary or static loading conditions are presented in Fig. 5. The numerical results agree very well with the analytical results for the steady state case. In the second numerical example, an infinitely long and functionally graded thick-walled hollow cylinder is considered, where the following radii R1 ¼ 8 102 m and R2 ¼ 10 102 m are selected. Heaviside step boundary condition is prescribed on the external surface of the hollow cylinder for the time-dependent thermal loading as a thermal shock. The inner surface is kept at zero temperature. Due to the symmetry in geometry and boundary conditions it is sufficient to analyze only a quarter of the cross section of the hollow cylinder. The total number of nodes used for the MLS approximation is 84 with 42 nodes lying on the global boundary as depicted in Fig. 6. For comparison purpose, a homogeneous material is first considered. The value of the thermal diffusivity used here is the same as in the previous numerical example. Also in this homogeneous case, an analytical solution [22] is known as: 1 X lnðr=R1 Þ J 2 ðR1 an ÞU0 ðran Þ hðr; tÞ ¼ T p T 20 lnðR2 =R1 Þ J0 ðR1 an Þ J02 ðR2 an Þ n¼1 expðja2n tÞ;
ð27Þ
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
501
0.8
0.7
0.6
Temperature θ
0.5
0.4
0.3
0.2
analytical: x1/a = 0.25 x1/a = 0.5
0.1
x1/a = 0.75 LBIE:
x1/a = 0.25 x1/a = 0.5
0
x1/a = 0.75 -0.1 0
10
20
30
40
50
60
70
Time t[sec] Fig. 3. Time variation of the temperature in a finite strip at three positions.
Fig. 4. Time variation of the temperature at position x1 =a ¼ 0:25 of the functionally graded finite strip.
Fig. 5. Temperature variation with x1 -coordinate for a functionally graded finite strip under steady state loading conditions.
where U0 ðran Þ ¼ J0 ðran ÞY0 ðan R2 Þ J0 ðan R2 ÞY0 ðran Þ;
J0 ðrÞY0 ðrR2 =R1 Þ J0 ðrR2 =R1 ÞY0 ðrÞ ¼ 0
and an are roots of the following transcendental equation
with J0 ðxÞ and Y0 ðxÞ being Bessel functions of first and second kind and zeroth order.
502
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504 22 R2
q=0 θ =ΤΗ(t)
28 R1 x2
θ =0 42 x1
1
q=0
7
Fig. 6. Boundary conditions and node distribution for a hollow cylinder.
In the numerical calculations, both Papoulis and StehfestÕs inversion algorithms are applied to obtain the temporal variation of the temperature. In the Papoulis method b ¼ 0:0025 is chosen. In the Papoulis method, some small discrepancies between analytical and numerical results are noted as can be seen in Fig. 7. In contrast, the StehfestÕs method yields an excellent agreement between
Fig. 8. Temperature variation with radial coordinate r at different time levels for a hollow cylinder with homogeneous material properties.
analytical and numerical results. Fig. 8 shows the variation of the temperature with radial coordinate at four different time instants. Finally, we consider the functionally graded hollow cylinder with the thermal diffusivity and conductivity being graded in the radial direction r. The parameters j0 , k0 and c required in the exponential spatial variation of the material constants are selected as the same as in the first numerical example. Numerical results for the time variation of the temperature are shown in Fig. 9. Similar to the results for a finite strip in the first example, the temperature level at interior points in the steady state increases with increasing c-value. Fig. 10 illustrates the variation of the temperature with the normalized radial coordinate (r=R1 ) in such a case.
6. Conclusions
Fig. 7. Time variation of the temperature at r ¼ 8:3 cm in the hollow cylinder with homogeneous material properties.
A local boundary integral equation method with the MLS approximation of physical fields is applied to transient heat conduction analysis in functionally graded materials. In this method, a general spatial variation of the material constants can be used for modelling the functionally graded
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
503
plemented for a simple circular subdomain to which the local boundary integral equation is related. Contrary to the conventional boundary integral equation methods, all integrands in the present formulation are regular. Thus, no special integration techniques are required to evaluate the integrals. The present local boundary integral equation method removes the well-known restriction of conventional BEM to problems with homogeneous material properties.
Acknowledgements
Fig. 9. Time variation of the temperature at r ¼ 9:0 cm in a functionally graded hollow cylinder.
The authors acknowledge the support by the Slovak Science and Technology Assistance Agency registered under number APVT-51-003702, and the Project for Bilateral Cooperation in Science and Technology supported jointly by the International Bureau of the German BMBF and the Ministry of Education of Slovak Republic under the project number SVK 01/020.
References
Fig. 10. Temperature variation with radial coordinate r for a functionally graded hollow cylinder.
materials. The initial-boundary value problem is solved in the Laplace transform domain with a subsequent numerical Laplace inversion to obtain time-dependent solutions. The local boundary integral equations are formulated in a boundarydomain form with a simple fundamental solution corresponding to the PoissonÕs equation. The boundary-domain formulation can be easily im-
[1] S. Suresh, A. Mortensen, Fundamentals of Functionally Graded Materials, Institute of Materials, London, 1998. [2] Y. Miyamoto, W.A. Kaysser, B.H. Rabin, A. Kawasaki, R.G. Ford, Functionally Graded Materials: Design, Processing and Applications, Kluwer Academic Publishers, Dordrecht, 1999. [3] N. Noda, Z.H. Jin, Thermal stress intensity factors for a crack in a strip of a functionally gradient material, Int. J. Solids Struct. 30 (1993) 1039–1056. [4] F. Erdogan, B.H. Wu, Crack problems in FGM layers under thermal stresses, J. Therm. Stresses 19 (1996) 237– 265. [5] Z.H. Jin, N. Noda, An internal crack parallel to the boundary of a nonhomogeneous half plane under thermal loading, Int. J. Eng. Sci. 31 (1993) 793–806. [6] Z.H. Jin, R.C. Batra, Stress intensity relaxation at the tip of an edge crack in a functionally graded material subjected to a thermal shock, J. Therm. Stresses 19 (1996) 317–339. [7] N. Noda, Z.H. Jin, A crack in functionally gradient materials under thermal shock, Arch. Appl. Mech. 64 (1994) 99–110. [8] Z.H. Jin, G.H. Paulino, Transient thermal stress analysis of an edge crack in a functionally graded material, Int. J. Fract. 107 (2001) 73–98. [9] A. Sutradhar, G.H. Paulino, L.J. Gray, Transient heat conduction in homogeneous and non-homogeneous materials
504
[10]
[11]
[12]
[13] [14]
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504 by the Laplace transform Galerkin boundary element method, Eng. Anal. Boundary Elem. 26 (2002) 119–132. R. Butterfield, An application of the boundary element method to potential flow problems in generally inhomogeneous bodies, in: C.A. Brebbia (Ed.), Recent Advances in Boundary Element Method, Pentech Press, London, 1978. M. Tanaka, K. Tanaka, Transient heat conduction problems in inhomogeneous media discretized by means of boundary-volume element, Nucl. Eng. Des. 60 (1980) 381– 387. T. Zhu, J.D. Zhang, S.N. Atluri, A local boundary integral equation (LBIE) method in computational mechanics, and a meshless discretization approach, Comput. Mech. 21 (1998) 223–235. S. De, K.J. Bathe, The method of finite spheres, Comput. Mech. 25 (2000) 329–345. T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, P. Krysl, Meshless methods: an overview and recent developments, Comput. Meth. Appl. Mech. Eng. 139 (1996) 3– 47.
[15] A. Papoulis, A new method of inversion of the Laplace transform, Quart. Appl. Math. 14 (1957) 405–414. [16] H. Stehfest, Algorithm 368: numerical inversion of Laplace transform, Comm. Assoc. Comput. Mach. 13 (1970) 47–49. [17] V. Sladek, J. Sladek, S.N. Atluri, R. Van Keer, Numerical integration of singularities in meshless implementation of local boundary integral equations, Comput. Mech. 25 (2000) 394–403. [18] Y.X. Mukherjee, S. Mukherjee, On boundary conditions in the element-free Galerkin method, Comput. Mech. 19 (1997) 264–270. [19] J. Balas, J. Sladek, V. Sladek, Stress Analysis by Boundary Element Methods, Elsevier, Amsterdam, 1989. [20] D.P. Gaver, Observing stochastic processes, and approximate transform inversion, Oper. Res. 14 (1966) 444–459. [21] B. Davies, B. Martin, Numerical inversion of the Laplace transform: a survey and comparison of methods, J. Comput. Phys. 33 (1979) 1–32. [22] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, Clarendon Press, Oxford, 1959.
Transient heat conduction analysis in functionally graded materials by the meshless local boundary integral equation method J. Sladek b
a,*
, V. Sladek a, Ch. Zhang
b
a Institute of Construction and Architecture, Slovak Academy of Sciences, 84220 Bratislava, Slovakia Department of Civil Engineering, University of Applied Sciences Zittau/G€orlitz, D-02763 Zittau, Germany
Abstract Advanced computational method for transient heat conduction analysis in continuously nonhomogeneous functionally graded materials (FGM) is proposed. The method is based on the local boundary integral equations with moving least square approximation of the temperature and heat flux. The initial-boundary value problem is solved by the Laplace transform technique. Both Papoulis and Stehfest algorithms are applied for the numerical Laplace inversion to obtain the time-dependent solutions. Numerical results are presented for a finite strip and a hollow cylinder with an exponential spatial variation of material parameters. Ó 2003 Elsevier B.V. All rights reserved. Keywords: Transient heat conduction analysis; Local boundary integral equations; Functionally graded materials (FGM); Laplace transform; MLS approximation
1. Introduction Due to their excellent thermal properties at high temperatures and their superior corrosion and wear resistances, ceramic materials are preferred and favored in many engineering structures and components. One major disadvantage of ceramics is their brittleness and low fracture toughness. To eliminate this disadvantage and to increase their toughness, ceramic materials are frequently reinforced by metallic particles or fibres, which yields
*
Corresponding author. Tel.: +421-2-54788662; fax: +421-254773548. E-mail address: [email protected] (J. Sladek).
new composite materials. To fabricate a structural component with a constant mixture ratio of ceramics and metals is not an efficient way due to specific requirements in material processing on individual parts of the material surfaces. A more sophisticated way is to have volume fractions of composite constituents varying continuously in space. This class of composite materials is referred to as functionally graded materials (FGMs) [1,2]. Since functionally graded materials are frequently used for structures under thermal load, it is indispensable to analyze their thermal properties. Due to the complexity of the problem, only a limited number of investigations dealing with heat conduction analysis in FGMs can be found in literature. Most of them used an exponential
0927-0256/$ - see front matter Ó 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2003.08.006
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
variation of material constants with Cartesian coordinates under stationary thermal loading conditions [3–5]. Transient heat transfer in FGMs with an exponential spatial variation of material constants has been examined by Jin and Batra [6], Noda and Jin [7] and Jin and Paulino [8]. Recently, Sutradhar et al. [9] presented a Laplace transform Galerkin BEM for 3-D transient heat conduction analysis by using the GreenÕs function approach where an exponential law for the FGMs was used. Since analytical methods are usually restricted to simple geometry and boundary conditions, solution of the boundary value problems for nonhomogeneous bodies requires special numerical methods due to the high mathematical complexity. Beside the well established finite element method (FEM), the boundary element method (BEM) provides an efficient and popular alternative to the FEM for solving certain class of boundary value problems. However, the pure BEM formulation can be applied only to problems where fundamental solution is available. For a general nonhomogeneous body the fundamental solution is unfortunately not available. One possibility to obtain a BEM formulation is based on the use of fundamental solutions for a fictitious homogeneous medium, as has been suggested for the first time by Butterfield [10] for potential flow problems. This approach, however, leads to a boundary-domain integral formulation with additional domain integrals for the gradients of primary fields to obtain a unique formulation. This approach has been applied to the heat conduction analysis in anisotropic nonhomogeneous media by Tanaka and Tanaka [11]. In this approach, some computational difficulties in the numerical implementation have been reported. Alternatively, the analyzed domain can be divided into small subdomains with a simple geometry such as circle and sphere. To each subdomain one can relate the fundamental solution to the corresponding material constants [12,13]. On the surface of subdomains the local boundary integral equations (LBIE) are applied. The temperature and heat flux on the local boundary and in the interior of the subdomain are approximated by the moving least square (MLS). The key idea of MLS interpolants is that it is only necessary to construct an array of nodes in the domain under consideration [14]. Thus, the method is completely element-free.
495
In this paper, a local boundary element method for transient heat conduction analysis is presented. The transient heat conduction problem is described by a parabolic partial differential equation. To eliminate the time dependence of the differential equation, the Laplace transform technique is used. A pure boundary integral formulation for the global boundary requires the application of the fundamental solution in the Laplace transform domain, which is not known for general FGM. In the present work, a simpler fundamental solution corresponding to the PoissonÕs equation is adopted which leads to a boundary-domain integral formulation. Since the boundary-domain formulation in the global BEM is numerically not efficient, we apply these integral equations to small subdomains, which result in a local boundary integral equation formulation. By choosing a simple geometric shape such as a circle for the subdomains, the arising domain integrals in the local BEM can be easily evaluated. It means that the boundarydomain integral formulation for simple subdomains does not bring any difficulties. Several quasi-static boundary value problems are solved for various values of the Laplace transform parameter. The Papoulis [15] and Stehfest [16] numerical inversion methods are applied to obtain the time-dependent solutions. The exponential spatial variation of thermal material properties is considered in numerical examples. Results are given for a finite strip and a hollow cylinder. Both stationary and transient loading (thermal shock) conditions are considered.
2. Local boundary integral equations Let us consider a continuously nonhomogeneous FGM with zero initial condition hðx; tÞ ¼ 0 at t ¼ 0, which is subjected to a thermal loading on its boundary. The thermal heat conduction problem in the FGM is described by the following governing equation h;ii ðx; tÞ ¼
1 oh k;i ðx; tÞ þ ðxÞh;i ðx; tÞ jðxÞ ot k
1 Qðx; tÞ; jðxÞ
ð1Þ
496
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
where hðx; tÞ is the temperature field, k and j stand for the thermal conductivity and diffusivity, and Qðx; tÞ is the density of the body heat sources, respectively. Applying the Laplace transform Z 1 L½hðx; tÞ ¼ hðx; pÞ ¼ hðx; tÞept dt 0
to the governing equation (1) we obtain p k;i h;ii ðx; pÞ hðx; pÞ þ ðxÞh;i ðx; pÞ j k 1 ¼ F ðx; pÞ; j
ð2Þ
where p is the Laplace transform parameter and F ðx; pÞ ¼ Qðx; pÞ þ hðx; 0Þ is the redefined body heat source in the Laplace transform domain with the initial-boundary condition for the temperature field. The solution of the governing equation (2) can be found in a weak form with a fundamental solution h ðx; yÞ ¼
1 1 ln 2p r
ð3Þ hðy; pÞ ¼
corresponding to the Poisson equation h;ii ðx; yÞ
þ dðx; yÞ ¼ 0;
where C is the boundary of the analysed domain X. Since the boundary-domain integral representation (5) for the Laplace transform of the temperature contains temperature gradients, the integral representation for the temperature gradients at interior points has to be added to obtain a unique integral formulation. It should be noted that the singularity of the kernel in the integral representation for the temperature gradients is increased to a hypersingularity. By applying the representation formulae for the temperature and its gradients to the global boundary of the considered domain, a unique global BEM formulation is obtained. Domain discretization and evaluation of hypersingular integrals make the global BEM cumbersome. Another approach is the use of the local boundary integral equations, related to the boundaries of small regular no-overlapping subdomains randomly spread in the analyzed domain. If, instead of the entire domain X, we consider a subdomain Xs , the following local boundary integral equation should also hold over the subdomain Z oXs
ð4Þ
where dðx; yÞ is the Dirac delta function and r is the distance of the field and source points, i.e., r ¼ jx yj. By the use of the Gauss divergence theorem to the weak form of Eq. (2) we obtain an integral representation for the temperature field in the Laplace transform domain Z oh hðy; pÞ ¼ ðx; pÞh ðx; yÞ dC C on Z oh ðx; yÞ dC hðx; pÞ on ZC p hðx; pÞh ðx; yÞ dX jðxÞ ZX k;i ðxÞh;i ðx; pÞh ðx; yÞ dX þ X k Z 1 F ðx; pÞh ðx; yÞ dX; ð5Þ þ X jðxÞ
þ þ
Z
Z Z Z
oh ðx; pÞh ðx; yÞ dC on
oh hðx; pÞ ðx; yÞ dC on oXs Xs
p hðx; pÞh ðx; yÞ dX jðxÞ
Xs
k;i ðxÞh;i ðx; pÞh ðx; yÞ dX k
Xs
1 F ðx; pÞh ðx; yÞ dX: jðxÞ
ð6Þ
The integral equation (6) is considered for small subdomains Xs X. Hence, no of the boundary densities in the Laplace transformed domain are required on the local boundary oXs as long as it lies entirely inside Xs . To reduce the number of unknowns on oXs the concept of a companion solution can be utilized successfully [12]. The companion solution is associated with the fundamental solution h ðx; yÞ and is defined as the solution to the following equations:
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
h~;ii ðx; yÞ ¼ 0 on X0s ; h~;ii ðx; yÞ ¼ h ðx; yÞ on oX0s ;
local boundary ∂Ωs=∂Ω's=Ls subdomain Ωs=Ω's
ð7Þ
where X0s and oX0s are defined in Fig. 1. For simplicity, X0s is taken as a circle in the present implementation. The modified test function
ri
Ωx x
1 r0 ln h~ ðx; yÞ ¼ h ðx; yÞ h~ðx; yÞ ¼ 2p r
Ls
node xi
has to satisfy the governing equation (4). It is the GreenÕs function vanishing on the boundary of the circular subdomain of radius r0 . Since the integral equation (6) is also valid for the modified fundamental solution h~ ðx; yÞ, we can obtain Z oh~ hðy; pÞ ¼ hðx; pÞ ðx; yÞ dC on oX Z s p hðx; pÞh~ ðx; yÞ dX jðxÞ Xs Z k;i ðxÞh;i ðx; pÞh~ ðx; yÞ dX þ k X Z s 1 þ F ðx; pÞh~ ðx; yÞ dX ð8Þ jðxÞ Xs for the interior source point y 2 X and Z oh~ hðf; pÞ ¼ hðx; pÞ ðx; fÞ dC on Ls Z oh~ lim hðx; pÞ ðx; yÞ dC y!f C on s Z oh ðx; pÞh~ ðx; fÞ dC þ Cs on Z p hðx; pÞh~ ðx; fÞ dX jðxÞ Xs Z k;i ðxÞh;i ðx; pÞh~ ðx; fÞ dX þ k X Z s 1 þ F ðx; pÞh~ ðx; fÞ dX jðxÞ Xs
497
Γs support of node xi
for the source point located on the global boundary f 2 Cs C. Note that oXs ¼ Ls [ Cs , where Ls is the circular part of oXs . The integral equation (9) is written in the limit form. Such an expression is appropriate [17] if the physical fields such as the temperature are given in a nonpolynomial approximation form.
Ωs'
∂Ωs'
Fig. 1. Local boundaries and the domain of definition of the MLS approximation for the trial function.
3. The MLS approximation scheme In general, a meshless method uses a local interpolation to represent the trial function with the values (or the fictitious values) of the unknown variable at some randomly located nodes. The moving least squares (MLS) approximation used in the present analysis may be considered as one of such schemes. Consider a subdomain Xx , the neighbourhood of a point x and denoted as the domain of definition of the MLS approximation for the trial function at x, which is located in the analyzed domain X. To approximate the distribution of the field hðxÞ in Xs over a number of randomly located nodes fxi g; i ¼ 1; 2; . . . ; n, the moving least square approximant hh ðxÞ of h for 8x 2 Xs is defined as hh ðxÞ ¼ pT ðxÞaðxÞ;
ð9Þ
∂Ωs=Ls ∪ Γs
8x 2 Xs ;
ð10Þ
where pT ðxÞ ¼ ½p1 ðxÞ; p2 ðxÞ; . . . ; pm ðxÞ is a complete monomial basis of order m and aðxÞ is a vector with components aa ðxÞ; a ¼ 1; 2; . . . ; m. The coefficient vector aðxÞ is determined by minimizing a weight discrete L2 norm which is defined as n X 2 wi ðxÞ½pT ðxi ÞaðxÞ h^i ; ð11Þ J ðxÞ ¼ i¼1
where wi ðxÞ is the weight function associated with the node i. Recall that n is the number of nodes in Xx for which the weight function wi ðxÞ > 0 and h^i are the fictitious nodal values, but not the nodal
498
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
values of the unknown trial function hh ðxÞ in general [18]. The stationarity of J ðxÞ in Eq. (11) with respect to aðxÞ leads to the following linear relation between aðxÞ and ^ h AðxÞaðxÞ ¼ BðxÞ^ h;
ð12Þ
where AðxÞ ¼
n X wi ðxÞpðxi ÞpT ðxi Þ; i¼1
) oh~ i b ðx; y Þ/ ðxÞ dC h^i ðpÞ / ðy Þ þ on oXs i¼1 n Z X p b ~ /i ðxÞ h ðx; y Þ ¼ jðxÞ Xs i¼1 k;j i ðxÞ/;j ðxÞ dX hi ðpÞ k Z 1 þ F ðx; pÞh~ ðx; yb Þ dX for yb 2 X Xs jðxÞ
n X
(
i
b
Z
BðxÞ ¼ ½w1 ðxÞpðx1 Þ; w2 ðxÞpðx2 Þ; . . . ; wn ðxÞpðxn Þ: Solving for aðxÞ from Eq. (12) and substituting it into Eq. (10) give a relation n X h¼ /i ðxÞh^i ; ð13Þ hh ðxÞ ¼ UT ðxÞ ^
ð17Þ and n X
(
Z Ls
UT ðxÞ ¼ pT ðxÞA1 ðxÞBðxÞ:
where d i ¼ kx xi k, ci is a constant controlling the shape of the weight function wi ðxÞ and ri is the size of the support domain. The size of the support should be large enough in the domain of definition to ensure the regularity of the matrix A. Finally, one can write the approximation formula for the Laplace transform of the temperature and its normal derivative as n X hðx; pÞ ¼ /i ðxÞh^i ðpÞ; ð15Þ i¼1
Z
Z Cs
oh~ ðx; yÞ/i ðxÞ dC on
oh~ ðx; fb Þ/i ðxÞ dC on ) nk /i;k ðxÞh~ ðx
f Þ dC h^i ðpÞ b
Cs
p b ~ /i ðxÞ ¼ h ðx f Þ jðxÞ Xs i¼1 k;j ðxÞ/i;j ðxÞ dX h^i ðpÞ k Z 1 F ðx; pÞh~ ðx; fb Þ dX for fb 2 Csq ; þ jðxÞ Xs n Z X
ð14Þ
n X oh ðx; pÞ ¼ nk ðxÞ h^i ðpÞ/i;k ðxÞ: on i¼1
/ ðf Þ þ lim
y!fb
þ
The Gaussian weight function is applied in the present work as: ( exp½ðd i =ci Þ2 exp½ðri =ci Þ2 ; 0 6 d i 6 ri i 1exp½ðri =ci Þ2 w ðxÞ ¼ 0; d i P ri ;
b
i¼1
i¼1
where
i
ð18Þ Csq is a part of the global boundary C, over which the heat flux is prescribed. If the temperature is prescribed on a part Csh we propose to collocate the approximation formula (15) by using n X
/i ðfb Þh^i ðpÞ ¼ ~hðfb ; pÞ for fb 2 Ch ;
ð19Þ
i¼1
ð16Þ
Eqs. (15) and (16) imply that both quantities h and oh=on can be approximated by the fictitious nodal values h^i ðpÞ. Substituting the approximation formulae (15) and (16) into the local boundary integral equations (8) and (9) we obtain a set of linear algebraic equations:
where ~hðfb ; pÞ is the Laplace transform of the prescribed temperature. 4. Numerical inversion of Laplace transform The time-dependent values of the transformed variables in the previous sections can be obtained
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
by an inverse transform. There are many inversion methods available for Laplace transform. As the Laplace transform inversion is an ill-posed problem, small truncation errors can be greatly magnified in the inversion process and lead to poor numerical results. In the present analysis, two numerical inversion methods are used. The first one is the most familiar Papoulis inversion method [15] and the second one is the sophisticated StehfestÕs algorithm [16], which are explained briefly in the following. If f ðpÞ is the Laplace transform of f ðtÞ, an approximate value fa of f ðtÞ in the Papoulis method is given by fa ðsÞ ¼
L X
Cm sinð2m 1Þs;
ð20Þ
m¼1
where 1 t ¼ ln coss b with a transformation of the original time interval tð0; 1Þ into the transformed interval sð0; p=2Þ. The coefficients Cm in Eq. (20) are obtained from the system of algebraic equations L X
Bnk Ck ¼ An ;
ð21Þ
499
where vi ¼ ð1ÞN =2þi
maxði;N X=2Þ
k N =2 ð2kÞ! : ðN =2 kÞ!k!ðk 1Þ!ði kÞ!ð2k iÞ! k¼½ðiþ1Þ=2 ð23Þ
Sutradhar et al. [9] have suggested to use N ¼ 10 for single precision arithmetic. It means that for each time t it is needed to solve N boundary value problems for the corresponding Laplace parameters, p ¼ i ln 2=t, with i ¼ 1; 2; . . . ; N . If M denotes the number of the time instants in which we are interested to know f ðtÞ, the number of the Laplace transform solutions f ðpj Þ is then M N , which can be much higher than that in the Papoulis method, where for selected L-values of the Laplace transform parameter we obtain fa in the entire time interval. Most of the methods for the numerical inversion of the Laplace transform require the use of complex valued Laplace transform parameter, and as a result, the application of complex arithmetic may lead to additional storage requirement and an increase in computing time. A critical study of various inversion algorithms is given by Davies and Martin [21].
k¼1
where
5. Numerical examples
4 An ¼ 22ðn1Þ bf fð2n 1Þbg; p
Though the proposed local boundary integral equation method has no restrictions on the spatial variation of the material parameters of the FGM, the numerical examples presented here use an exponential variation of the material properties with Cartesian coordinates. In such a special case the equations for the heat conduction problems are substantially simplified in comparison with a general nonhomogeneous case. For comparison purposes it is desirable to have reliable accurate results obtained by another independent method. Analytical solutions are of course the best choice if they are available. First, we consider a finite strip with a unidirectional variation of the thermal conductivity and diffusivity. The same exponential spatial variation for both parameters is taken
and the coefficients of triangular matrix Bnk can be found in Papoulis [15] and Balas et al. [19]. The number of terms L in the expansion series (20) is selected in an interval ð12; 16Þ to obtain a high accuracy of the Laplace inversion. The parameter b is determined by the maximum value of the time interval and the practical limit of smax used is smax ffi 70 deg [19]. The StehfestÕs algorithm originates from Gaver [20]. Accordingly, an approximate value fa of the inverse f ðtÞ for a specific time t is given by N ln 2 X ln 2 i ; fa ðtÞ ¼ vi f t i¼1 t
ð22Þ
500
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
kðxÞ ¼ k0 ecx1 ; jðxÞ ¼ j0 ecx1
ð24Þ
with j0 ¼ 0:17 104 m2 s1 and k0 ¼ 17 W m1 deg1 . Two different exponential parameters c ¼ 0:2 and 0.5 cm1 are selected in numerical calculation. On both opposite sides parallel to the x2 -axis two different temperatures are prescribed. One side is kept to zero temperature and the other has the Heaviside step time variation hðx2 ; tÞ ¼ T H ðtÞ with T ¼ 1 deg. On the lateral sides of the strip the heat flux vanishes. In numerical calculations, a square with a side-length a ¼ 0:04 m and a regular node distribution with equal 16 boundary and interior nodes are used for the MLS approximation (see Fig. 2). The special case with an exponential parameter c ¼ 0 corresponds to a homogeneous material. In such a case an analytical solution is available 1 x1 2 X T cos np npx1 hðx1 ; tÞ ¼ T þ sin n a p n a jn2 p2 t exp ; ð25Þ a2 which can be used to check the accuracy of the present numerical method. Numerical results are computed at three locations along the x1 -axis with x1 =a ¼ 0:25, 0.5 and 0.75. Both Papoulis and StehfestÕs numerical inversion methods of the Laplace transform are applied. In the Papoulis method 15 different Laplace transform parameters and b ¼ 0:01 are used. In the StehfestÕs method
x2 q=0
9
13
θ=0
θ =ΤΗ(t)
x1
1
q=0
5
Fig. 2. Boundary conditions and node distribution for a finite strip.
N ¼ 10 is chosen. For the present case, both methods give nearly identical results. Numerical results for the temperature distribution are presented in Fig. 3. An excellent agreement between numerical and analytical results is obtained. Next, we continue to analyze the heat conduction problem in the FGM with c ¼ 0:2 and 0.5 cm1 as mentioned above. The time variation of the temperature at the position x1 =a ¼ 0:25 for three c-values is presented in Fig. 4. Here again, both inversion methods give almost identical results. With increasing c-values the thermal conductivity increases and a higher level of temperature at the considered position is obtained in the large time range corresponding to a steady state. For steady state an analytical solution can be obtained as ecx1 1 hðx1 Þ ¼ T ca : ð26Þ e 1 Analytical and numerical results corresponding to stationary or static loading conditions are presented in Fig. 5. The numerical results agree very well with the analytical results for the steady state case. In the second numerical example, an infinitely long and functionally graded thick-walled hollow cylinder is considered, where the following radii R1 ¼ 8 102 m and R2 ¼ 10 102 m are selected. Heaviside step boundary condition is prescribed on the external surface of the hollow cylinder for the time-dependent thermal loading as a thermal shock. The inner surface is kept at zero temperature. Due to the symmetry in geometry and boundary conditions it is sufficient to analyze only a quarter of the cross section of the hollow cylinder. The total number of nodes used for the MLS approximation is 84 with 42 nodes lying on the global boundary as depicted in Fig. 6. For comparison purpose, a homogeneous material is first considered. The value of the thermal diffusivity used here is the same as in the previous numerical example. Also in this homogeneous case, an analytical solution [22] is known as: 1 X lnðr=R1 Þ J 2 ðR1 an ÞU0 ðran Þ hðr; tÞ ¼ T p T 20 lnðR2 =R1 Þ J0 ðR1 an Þ J02 ðR2 an Þ n¼1 expðja2n tÞ;
ð27Þ
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
501
0.8
0.7
0.6
Temperature θ
0.5
0.4
0.3
0.2
analytical: x1/a = 0.25 x1/a = 0.5
0.1
x1/a = 0.75 LBIE:
x1/a = 0.25 x1/a = 0.5
0
x1/a = 0.75 -0.1 0
10
20
30
40
50
60
70
Time t[sec] Fig. 3. Time variation of the temperature in a finite strip at three positions.
Fig. 4. Time variation of the temperature at position x1 =a ¼ 0:25 of the functionally graded finite strip.
Fig. 5. Temperature variation with x1 -coordinate for a functionally graded finite strip under steady state loading conditions.
where U0 ðran Þ ¼ J0 ðran ÞY0 ðan R2 Þ J0 ðan R2 ÞY0 ðran Þ;
J0 ðrÞY0 ðrR2 =R1 Þ J0 ðrR2 =R1 ÞY0 ðrÞ ¼ 0
and an are roots of the following transcendental equation
with J0 ðxÞ and Y0 ðxÞ being Bessel functions of first and second kind and zeroth order.
502
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504 22 R2
q=0 θ =ΤΗ(t)
28 R1 x2
θ =0 42 x1
1
q=0
7
Fig. 6. Boundary conditions and node distribution for a hollow cylinder.
In the numerical calculations, both Papoulis and StehfestÕs inversion algorithms are applied to obtain the temporal variation of the temperature. In the Papoulis method b ¼ 0:0025 is chosen. In the Papoulis method, some small discrepancies between analytical and numerical results are noted as can be seen in Fig. 7. In contrast, the StehfestÕs method yields an excellent agreement between
Fig. 8. Temperature variation with radial coordinate r at different time levels for a hollow cylinder with homogeneous material properties.
analytical and numerical results. Fig. 8 shows the variation of the temperature with radial coordinate at four different time instants. Finally, we consider the functionally graded hollow cylinder with the thermal diffusivity and conductivity being graded in the radial direction r. The parameters j0 , k0 and c required in the exponential spatial variation of the material constants are selected as the same as in the first numerical example. Numerical results for the time variation of the temperature are shown in Fig. 9. Similar to the results for a finite strip in the first example, the temperature level at interior points in the steady state increases with increasing c-value. Fig. 10 illustrates the variation of the temperature with the normalized radial coordinate (r=R1 ) in such a case.
6. Conclusions
Fig. 7. Time variation of the temperature at r ¼ 8:3 cm in the hollow cylinder with homogeneous material properties.
A local boundary integral equation method with the MLS approximation of physical fields is applied to transient heat conduction analysis in functionally graded materials. In this method, a general spatial variation of the material constants can be used for modelling the functionally graded
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504
503
plemented for a simple circular subdomain to which the local boundary integral equation is related. Contrary to the conventional boundary integral equation methods, all integrands in the present formulation are regular. Thus, no special integration techniques are required to evaluate the integrals. The present local boundary integral equation method removes the well-known restriction of conventional BEM to problems with homogeneous material properties.
Acknowledgements
Fig. 9. Time variation of the temperature at r ¼ 9:0 cm in a functionally graded hollow cylinder.
The authors acknowledge the support by the Slovak Science and Technology Assistance Agency registered under number APVT-51-003702, and the Project for Bilateral Cooperation in Science and Technology supported jointly by the International Bureau of the German BMBF and the Ministry of Education of Slovak Republic under the project number SVK 01/020.
References
Fig. 10. Temperature variation with radial coordinate r for a functionally graded hollow cylinder.
materials. The initial-boundary value problem is solved in the Laplace transform domain with a subsequent numerical Laplace inversion to obtain time-dependent solutions. The local boundary integral equations are formulated in a boundarydomain form with a simple fundamental solution corresponding to the PoissonÕs equation. The boundary-domain formulation can be easily im-
[1] S. Suresh, A. Mortensen, Fundamentals of Functionally Graded Materials, Institute of Materials, London, 1998. [2] Y. Miyamoto, W.A. Kaysser, B.H. Rabin, A. Kawasaki, R.G. Ford, Functionally Graded Materials: Design, Processing and Applications, Kluwer Academic Publishers, Dordrecht, 1999. [3] N. Noda, Z.H. Jin, Thermal stress intensity factors for a crack in a strip of a functionally gradient material, Int. J. Solids Struct. 30 (1993) 1039–1056. [4] F. Erdogan, B.H. Wu, Crack problems in FGM layers under thermal stresses, J. Therm. Stresses 19 (1996) 237– 265. [5] Z.H. Jin, N. Noda, An internal crack parallel to the boundary of a nonhomogeneous half plane under thermal loading, Int. J. Eng. Sci. 31 (1993) 793–806. [6] Z.H. Jin, R.C. Batra, Stress intensity relaxation at the tip of an edge crack in a functionally graded material subjected to a thermal shock, J. Therm. Stresses 19 (1996) 317–339. [7] N. Noda, Z.H. Jin, A crack in functionally gradient materials under thermal shock, Arch. Appl. Mech. 64 (1994) 99–110. [8] Z.H. Jin, G.H. Paulino, Transient thermal stress analysis of an edge crack in a functionally graded material, Int. J. Fract. 107 (2001) 73–98. [9] A. Sutradhar, G.H. Paulino, L.J. Gray, Transient heat conduction in homogeneous and non-homogeneous materials
504
[10]
[11]
[12]
[13] [14]
J. Sladek et al. / Computational Materials Science 28 (2003) 494–504 by the Laplace transform Galerkin boundary element method, Eng. Anal. Boundary Elem. 26 (2002) 119–132. R. Butterfield, An application of the boundary element method to potential flow problems in generally inhomogeneous bodies, in: C.A. Brebbia (Ed.), Recent Advances in Boundary Element Method, Pentech Press, London, 1978. M. Tanaka, K. Tanaka, Transient heat conduction problems in inhomogeneous media discretized by means of boundary-volume element, Nucl. Eng. Des. 60 (1980) 381– 387. T. Zhu, J.D. Zhang, S.N. Atluri, A local boundary integral equation (LBIE) method in computational mechanics, and a meshless discretization approach, Comput. Mech. 21 (1998) 223–235. S. De, K.J. Bathe, The method of finite spheres, Comput. Mech. 25 (2000) 329–345. T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, P. Krysl, Meshless methods: an overview and recent developments, Comput. Meth. Appl. Mech. Eng. 139 (1996) 3– 47.
[15] A. Papoulis, A new method of inversion of the Laplace transform, Quart. Appl. Math. 14 (1957) 405–414. [16] H. Stehfest, Algorithm 368: numerical inversion of Laplace transform, Comm. Assoc. Comput. Mach. 13 (1970) 47–49. [17] V. Sladek, J. Sladek, S.N. Atluri, R. Van Keer, Numerical integration of singularities in meshless implementation of local boundary integral equations, Comput. Mech. 25 (2000) 394–403. [18] Y.X. Mukherjee, S. Mukherjee, On boundary conditions in the element-free Galerkin method, Comput. Mech. 19 (1997) 264–270. [19] J. Balas, J. Sladek, V. Sladek, Stress Analysis by Boundary Element Methods, Elsevier, Amsterdam, 1989. [20] D.P. Gaver, Observing stochastic processes, and approximate transform inversion, Oper. Res. 14 (1966) 444–459. [21] B. Davies, B. Martin, Numerical inversion of the Laplace transform: a survey and comparison of methods, J. Comput. Phys. 33 (1979) 1–32. [22] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, Clarendon Press, Oxford, 1959.