- Email: [email protected]
A local BIEM for analysis of transient heat conduction with nonlinear source terms in FGMs
Engineering Analysis with Boundary Elements 28 (2004) 1–11 www.elsevier.com/locate/enganabound
A local BIEM for analysis of transient heat conduction with nonlinear source terms in FGMs Jan Sladeka,*, Vladimir Sladeka, Ch. Zhangb b
a Institute of Construction and Architecture, Slovak Academy of Sciences, 84220 Bratislava, Slovakia Department of Civil Engineering, University of Applied Sciences Zittau/Go¨rlitz, D-02763 Zittau, Germany
Received 6 February 2003; revised 13 May 2003; accepted 14 May 2003
Abstract The diffusion equation with nonlinear heat source intensity in functionally graded materials (FGMs) is considered. In FGMs the thermal material properties are dependent on spatial coordinates. For transient or steady-state heat problems in FGMs the conventional boundary integral equation method or boundary element method cannot be applied due to the lack of a fundamental solution. In this paper, a local boundary integral equation method is proposed to analyse a temperature distribution in a nonhomogeneous body under a microwave heating. To eliminate time variable in the heat equation, the Laplace transform technique is used. The boundary-domain integral formulation with a simple fundamental solution corresponding to the Laplace operator is related to all subdomains which cover the analysed domain. If such integral equations are considered on small subdomains with a simple geometry (circle), domain integrals can be easily evaluated. Physical fields (temperature, heat flux) on the local boundary and in the interior of the subdomain are approximated by the moving least-square. The method is completely element free. q 2003 Elsevier Ltd. All rights reserved. Keywords: Hot spots; Transient heat conduction analysis; Local boundary integral equations; Functionally graded materials; Laplace transform; Moving leastsquare approximation; Absorbtivity
1. Introduction Transient heat conduction equations with nonlinear source terms arise as governing equations in many different areas of mathematical physics. Examples in the literature include microwave heating process, mass transport in groundwater, and ignition. The use of microwave heating has attracted much attention in industrial processes due to its high efficiency. A crucial point in microwave heating is to analyse hot spots, regions where temperature is enormously running. This phenomenon is occurred due to dependence of material properties on temperature. Namely, thermal absorptivity increases with temperature. Hot spots out of control can damage the sample in an industrial process. The first attempt to analyse hot spots in a homogeneous body is given by Zhu et al. [1]. The dual reciprocity boundary element method (BEM) is proposed there. In the DRBEM, some interior nodes have to be supplemented to boundary * Corresponding author. Tel.: þ 421-7-54788662; fax: þ 421-754772494. E-mail address: [email protected] (J. Sladek). 0955-7997/04/$ - see front matter q 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0955-7997(03)00093-6
nodes for an adequate spatial approximation of physical fields in transforming domain integrals into boundary integrals. Functionally graded materials (FGMs) are preferred and favoured in many engineering structures and components due to their excellent thermal properties [2,3]. In FGMs the volume fractions of composite constituents are varying continuously in space. Solution of boundary or initial boundary value problems for FGMs requires special numerical methods due to the high mathematical complexity caused by the material nonhomogeneity. Beside the wellestablished finite element method (FEM), the boundary integral equation method (BIEM) or BEM provides an efficient and popular alternative to the FEM for solving certain class of boundary or initial boundary value problems. However, the pure BIEM/BEM formulation can be efficiently applied only to problems where fundamental solutions or Green’s functions are available in simple forms. Unfortunately, for general FGMs the corresponding fundamental solutions are either not available or they are mathematically too complicated. To avoid this difficulty,
2
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
one possibility to obtain a BEM formulation is the use of fundamental solutions for a fictitious homogeneous medium, as has been suggested for the first time by Butterfield [4] for potential flow problems. This approach, however, leads to a boundary-domain integral formulation with additional domain integrals representing the effects of the gradients of primary fields to obtain a unique formulation. Some computational difficulties in the numerical implementation of the boundary-domain integral method have been reported. A BEM based on Green’s function method for three-dimensional transient heat conduction analysis in FGMs is presented in Ref. [5], where the method is restricted to exponentially varying material properties. An alternative way is recently suggested in Refs. [6 – 9] by using the local BIEM. In this method, the analysed domain is divided into small subdomains with a simple geometry such as circles and spheres. In each sufficiently small subdomain one can apply the fundamental solutions corresponding to a homogeneous material with adequate local material constants [6 – 9]. On the boundary of subdomains the local boundary integral equations (LBIEs) are applied. The temperature and the heat flux on the local boundary and in the interior of the subdomain are approximated by the moving least-square (MLS) scheme. The key idea of MLS interpolants is that it is only necessary to construct an array of nodes in the domain under consideration [10]. Thus, the method is completely element-free. In this paper, a local BIEM for microwave heating analysis in FGMs is presented. The transient heat conduction problem is described by a parabolic partial differential equation, while its steady-state is governed by an elliptic one. To eliminate the time dependence of the differential equation for transient heat conduction analysis, 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 yet not known to our knowledge for general FGMs. 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 integral formulation in the global BIEM/BEM is numerically not efficient, we apply these integral equations to small subdomains, which result in a LBIE formulation. By choosing a simple geometric shape such as a circle for the subdomains, the arising domain integrals in the local BIEM/BEM can be easily evaluated. This means that the boundary-domain integral formulation for simple subdomains does not bring any numerical difficulties. To describe the material gradation of FGMs, an exponential spatial variation of thermal material properties is applied in numerical calculations. The temperature dependence of the thermal absorptivity is approximated by a power law. Numerical examples are presented for both steady-state and transient microwave heating problems.
2. Local boundary integral equations Microwave heating of a continuously nonhomogeneous body can be modelled by the nonstationary forced heat equation [1]
u;ii ðx; tÞ 2 ¼2
k 1 ›u ðx; tÞ þ ;i ðxÞu;i ðx; tÞ kðxÞ ›t k
1 bðuÞlEl2 kðxÞ
ð1Þ
and Maxwell’s equations, as governing equations for diffusion of heat and propagation of electromagnetic field including its absorption. In Eq. (1) t is the temporal variable; uðx; tÞ; the temperature field; k and k; the thermal conductivity and diffusivity, respectively; bðuÞ; the thermal absorptivity; lEl is the amplitude of the electric field. For simplicity the electric conductivity, magnetic and electric permeabilities are considered to be constant. Only thermal conductivity and absorptivity are considered to be dependent on temperature. The amplitude of the electric field is exponentially dependent on spatial variables [1] lEl ¼ expð2gz=2Þ;
ð2Þ
for decay from an incident boundary at z ¼ 0; where g is the decay constant, and thus Eq. (1) is uncoupled from Maxwell’s equations. Generally g should be complex propagation constant. In such a case the coupled system of both equations has to be solved. The simplified model for microwave heating is then determined by
u;ii ðx; tÞ 2 ¼2
1 ›u k ðx; tÞ þ ;i ðxÞu;i ðx; tÞ kðxÞ ›t k
1 bðuÞexpð2gzÞ: kðxÞ
ð3Þ
The microwave heating problem is described by a nonlinear parabolic partial differential equation (3). Analytical solutions exist only for simple geometrical and loading conditions, with uniform material properties and g ¼ 0: For general cases, Eq. (3) has to be solved numerically. A dual reciprocity BEM has been presented in Ref. [1], where numerical results are given for uniform material properties. If the thermal absorptivity is a nonlinear function of temperature, a linearization scheme has to be applied to Eq. (3). This can be done by using the following first-order Taylor-series expansion [11,12] ›b ~ ~ bðuÞ ¼ bðuÞ þ ðu 2 uÞ ; ð4Þ ›u u¼u~ where u~ is the previously iterated solution at a particular time. To eliminate the time variable in the differential equation (3), the Laplace transform technique is used.
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
Applying the Laplace transform ð1 L½uðx; tÞ ¼ uðx; pÞ ¼ uðx; tÞe2pt dt;
the evaluation of hypersingular integrals make the global boundary-domain formulation somewhat cumbersome. Another approach is the application of the LBIEs to the boundaries of small regular no-overlapping subdomains randomly spread in the analysed domain. If we consider a subdomain Vs instead of the entire domain V; the following LBIE should also valid over the subdomain
0
to the linearized governing equation (3) we obtain k p u;ii ðx; pÞ 2 uðx; pÞ þ ;i ðxÞu;i ðx; pÞ þ f1 ðx; z; u~Þuðx; pÞ k k pÞ; ¼ 2Fðx; where p is the Laplace transform parameter 1 ›b expð2gzÞ f1 ðx; z; u~Þ ¼ kðxÞ ›u u¼u~
ð5Þ
ð6Þ
and z; u~Þ ¼ 1 bðu~Þ 2 u~ ›b u¼u~ expð2gzÞ þ 1 uðx;0Þ; Fðx; kðxÞ ›u kðxÞ ð7Þ is the redefined body heat source in the Laplace transform domain with the initial condition for the temperature field. The solution of the governing equation (5) can be found in a weak form with a fundamental solution or Green’s function Gðx; yÞ ¼
1 1 ln ; 2p r
ð8Þ
corresponding to the Poisson’s equation G;ii ðx; yÞ þ dðx; yÞ ¼ 0;
3
ð9Þ
where dðx; yÞ is the Dirac delta function and r is the distance between the field and the source points, i.e. r ¼ lx 2 yl: By the use of the Gauss divergence theorem to the weak form of Eq. (5) we obtain an integral representation for the temperature field in the Laplace transform domain ð ›u ð ›G ðx; pÞGðx; yÞdG 2 ðx; yÞdG uðy; pÞ ¼ uðx; pÞ ›n G ›n G ð p 2 f1 ðx; z; u~Þ uðx; pÞGðx; yÞdV 2 V kðxÞ ð k ;i ðxÞu;i ðx; pÞGðx; yÞdV þ V k ð z; u~ÞGðx; yÞdV; þ ð10Þ Fðx; V
where G is the boundary of the analysed domain V: Since the boundary-domain integral representation (10) for the Laplace transform of the temperature contains temperature gradients, an 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 a unique global boundary-domain formulation is obtained. Domain discretization and
uðy; pÞ ¼
ð
ð ›u ›G ðx;pÞGðx; yÞdG 2 uðx;pÞ ðx; yÞdG ›n ›Vs ›n ›Vs ð p 2 f1 ðx; z; u~Þ uðx; pÞGðx;yÞdV 2 Vs kðxÞ ð k ;i ðxÞu;i ðx; pÞGðx;yÞdV þ Vs k ð þ u~ÞGðx;yÞdV: ð11Þ Fðx;z; Vs
The integral equation (11) is considered on small subdomains Vs , V: Hence, no of the boundary densities in the Laplace transformed domain are prescribed on the local boundary ›Vs as long as it lies entirely inside Vs : Both quantities (temperature and heat flux) are unknowns on ›Vs : To reduce the number of unknowns on ›Vs the concept of a companion solution can be utilized successfully [6]. The companion solution is associated with the fundamental solution Gðx; yÞ and is defined as the solution to the following equations on V0s ;
~ ;ii ðx;yÞ ¼ 0 G
~ ;ii ðx;yÞ ¼ Gðx; yÞ G
on ›V0s ;
ð12Þ
where V0s and ›V0s are depicted in Fig. 1. For simplicity, V0s is taken as a circle in the present implementation. The modified test function ~ p ðx; yÞ ¼ Gðx; yÞ 2 Gðx; ~ yÞ ¼ 1 ln r0 G 2p r has to satisfy the governing equation (9), and it is the modified fundamental solution or Green’s function vanishing on the boundary of the circular subdomain of radius r0 : Since the integral equation (11) is also valid for the modified ~ p ðx; yÞ; we obtain fundamental solution G ð
~p ›G ðx; yÞdG uðx; pÞ ›n ›Vs ð p ~ p ðx; yÞdV 2 f1 ðx; z; u~Þ uðx; pÞG 2 Vs kðxÞ ð k ;i ~ p ðx; yÞdV ðxÞu;i ðx; pÞG þ Vs k ð ~ p ðx; yÞdV; z; u~ÞG þ ð13Þ Fðx;
uðy; pÞ ¼ 2
Vs
4
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
Fig. 1. Local boundaries and the domain of definition of the MLS approximation for the trial function.
for interior source points y [ V and
uðz; pÞ ¼ 2
ð
described by a power law of the form
ð ~p ›G ðx; zÞdG 2 lim uðx; pÞ uðx; pÞ y!z Gs ›n Ls
bðuÞ ¼ b0 un :
ð ›u ~p ›G ~ p ðx; zÞdG ðx; yÞdG þ ðx; pÞG ›n G s ›n ð p ~ p ðx; zÞdV 2 f1 ðx; z; u~Þ uðx; pÞG 2 Vs kðxÞ ð k ð ;i ~ p ðx; zÞdV þ z; u~Þ ðxÞu;i ðx; pÞG þ Fðx; Vs k Vs ~ p ðx; zÞdV; G
ð14Þ
for source points located on the global boundary z [ Gs , G: In the latter case, ›Vs ¼ Ls < Gs ; where Ls is the circular part of ›Vs : The strongly singular integral in the integral equation (14) is written in the limit form. One can see that only temperature is unknown in the local boundary integral equation (13) at the interior source point. The heat flux is eliminated there due to convenient behaviour of the Green ~ p ðx; yÞ on the circle. For the source point on the function G global boundary the LBIE (14) contains both quantities (temperature and heat flux), however one of them is prescribed in a well posed boundary value problem. Next, we consider the steady-state or stationary heat conduction problem without time-dependence in the nonlinear parabolic partial differential equation (3). In this case, the governing equation (3) can be simplified to k 1 u;ii ðxÞ þ ;i ðxÞu;i ðxÞ ¼ 2 bðuÞexpð2gzÞ: k kðxÞ
ð15Þ
The steady-state heat conduction equation (3) is an elliptic partial differential equation and it is suitable for describing hot spots problems. Hot spots are the localized areas of high temperature caused by the graduation of the thermal absorptivity with temperature. To simplify the analysis, we assume that the thermal absorptivity is
ð16Þ
Note here that many engineering materials can be adequately approximated by the power law (16). According to Eq. (16), a thermal runway is possible only if n $ 1 [1]. Substitution of Eq. (16) into Eq. (15) leads to a eigen-value problem for the determination of the critical value of the absorptivity coefficient b0 ; which can be written as
u;ii ðxÞ þ
k;i 1 ðxÞu;i ðxÞ þ b ½uðxÞn expð2gzÞ ¼ 0: k kðxÞ 0
ð17Þ
By using the first-order Taylor-series expansion (4), Eq. (17) can be written in a linearized form as
u;ii ðxÞ þ
k;i ðxÞu;i ðxÞ þ f1 ðx; z; u~ÞuðxÞ ¼ 2Fðx; z; u~Þ; k
ð18Þ
where f1 ðx; z; u~Þ ¼
1 ›b 1 expð2gzÞ ¼ expð2gzÞb0 nu~n21 ; ~ kðxÞ ›u u¼u kðxÞ
Fðx; z; u~Þ ¼
1 ~ ›b bðuÞ 2 u~ u¼u~ expð2gzÞ kðxÞ ›u
¼
1 expð2gzÞb0 ð1 2 nÞu~n ; kðxÞ
in which u~ is the previously iterated solution for a nonlinear variation of the thermal absorptivity with temperature. Note here that no iteration is needed if the coefficient n ¼ 0 and 1, i.e. when the thermal absorptivity is independent of temperature or linearly dependent on temperature. The solution of the eigen-value problem (18) can be given in an integral form through the LBIEs. The procedure to obtain the LBIEs is the same as in the nonstationary case. By omitting integrals representing the nonstationary term in the BIEs (13) and (14), a system of LBIEs at interior and
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
boundary nodes can be obtained as ð ð ~p ›G ~ p ðx;yÞdV ðx;yÞdG þ uðyÞ¼2 uðxÞ f1 ðx;z; u~ÞuðxÞG ›n ›Vs Vs ð k ð ;i ~ p ðx;yÞdV þ ~ p ðx;yÞdV ðxÞu;i ðxÞG þ Fðx;z; u~ÞG Vs k Vs ð19Þ and
given by a fitting curve obtained by the least square method from the nodal values. The stationarity of J with respect to aðxÞ leads to the following linear relation
›J=›a ¼ AðxÞaðxÞ 2 BðxÞu^ ¼ 0; where matrices AðxÞ and BðxÞ are defined AðxÞ ¼ PT WP ¼
n X
wi ðxÞpðxi ÞpT ðxi Þ;
i¼1
ð
ð ~p ~p ›G ›G ðx; zÞdG 2lim ðx;yÞdG uðxÞ y!z Gs ›n ›n Ls ð ›u ð ~ p ðx; zÞdG 2 ~ p ðx; zÞdV ðxÞG þ f1 ðx;z; u~ÞuðxÞG G s ›n Vs ð k ð ;i ~ p ðx; zÞdV þ ~ p ðx; zÞdV: ðxÞu;i ðxÞG þ Fðx;z; u~ÞG Vs k Vs ð20Þ
uðzÞ¼2
uðxÞ
Here again, the strongly singular integral in Eq. (20) is written in the limit form, which is convenient for the numerical evaluation of the strongly singular integral when the approximation of the boundary density is not known in a closed form [13]. The remaining problem is how to represent the values of the unknown function on the local boundary and in the interior of a subdomain. This will be explained and discussed in Section 3.
3. The MLS approximation 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. To explain the key idea of the MLS approximation, let us consider now a subdomain Vx of the analysed domain V: The subdomain is the neighbourhood of a point x and denoted as the domain of definition of the MLS approximation for the trial function at x: To approximate the distribution of the temperature field uðxÞ in Vs over a number of randomly located nodes {xi } ði ¼ 1; 2; …; nÞ; the moving least square approximant uh ðxÞ of u for ;x [ Vs is defined as
uh ðxÞ ¼ pT ðxÞaðxÞ; T
1
;x [ Vs ; 2
ð21Þ
n X
BðxÞ ¼ PT W ¼ ½w1 ðxÞpðx1 Þ; w2 ðxÞpðx2 Þ; …; wn ðxÞpðxn Þ: The approximated function can be written as [10]
uh ðxÞ ¼ FT ðxÞ·u^ ¼
wi ðxÞ½pT ðxi ÞaðxÞ 2 u^i 2 ;
i¼1
where wi ðxÞ is the weight function associated with the node i and u^i are certain nodal values. Approximation values are
n X
fi ðxÞu^i ;
ð22Þ
i¼1
where
fi ðxÞ ¼
m X
pj ðxÞ½A21 ðxÞBðxÞji ;
j¼1
is the shape function associated with node i: Since the shape functions do not exhibit the Kronecker-delta property fi ðxj Þ dij ; the nodal values are fictitious nodal values which are different from real nodal values. The number of nodes n used for the approximation of uðxÞ is prescribed by the weight function wi ðxÞ: The Gaussian weight function is applied in the present work 8 > exp½2ðdi =ci Þ2 2 exp½2ðr i =ci Þ2 > < ; 0 # di # ri 1 2 exp½2ðr i =ci Þ2 wi ðxÞ ¼ ; > > : 0; di $ ri ð23Þ where di ¼ kx 2 xi k is the distance between the integration point and the source point, ci is a constant controlling the shape of the weight function wi ðxÞ; and r i is the size of the support domain. For a sharper Gaussian weight function the number of nodes used for approximation is smaller than in an opposite case. A necessary condition for a well-defined MLS approximation is that at least m weight functions are nonzero (i.e. n $ m) due to the regularity of the matrix A. Finally, the approximation formulae for the Laplace transform of the temperature and its normal derivative can be written as
m
where p ðxÞ ¼ ½p ðxÞ; p ðxÞ; …; p ð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 weighted discrete L2 norm, defined as JðxÞ ¼
5
uðx; pÞ ¼
n X
fi ðxÞu^i ðpÞ;
ð24Þ
i¼1 n X ›u ðx; pÞ ¼ nk ðxÞ u^i ðpÞfi;k ðxÞ: ›n i¼1
ð25Þ
Eqs. (24) and (25) imply that both quantities u and ›u=›n can be approximated by the fictitious nodal values u^i ðpÞ: Substituting the approximation formulae (24) and (25) into the local boundary integral Eqs. (13) and (14) we obtain a set
6
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
of linear algebraic equations ( ) n ð X ~p ›G i b b i ðx; y Þf ðxÞdG u^i ðpÞ f ðy Þ þ › n › V s i¼1 ¼2
n ð X i¼1
~ p ðx; yb Þ p fi ðxÞ 2 f1 ðx; z; u~Þfi ðxÞ G kðxÞ Vs
ð k;j ~ p ðx; yb ÞdV; z; u~ÞG ðxÞfi;j ðxÞ dVu^i ðpÞ þ Fðx; k Vs
2
for yb [ V and n X i¼1
fa ðtÞ ¼
fi ðzb Þ þ limb
ð
y!z
i
f ðxÞdG 2
Gs
ð
i¼1
2
ð ›G ~p ~p ›G ðx; yÞfi ðxÞdG þ ðx; zb Þ ›n Ls › n )
ð Gs
¼2
~ p ðx; zb ÞdG nk fi;k ðxÞG
~ p ðx; zb Þ G Vs
u ðpÞ
p i f ðxÞ 2 f1 ðx; z; u~Þfi ðxÞ kðxÞ
ð k;j ~ p ðx; zb ÞdV; z; u~ÞG ðxÞfi;j ðxÞ dVu^i ðpÞ þ Fðx; k Vs ð27Þ
where Gsq is a part of the global boundary G; over which the heat flux is prescribed. If the temperature is prescribed on a part Gsu we propose to collocate the approximation formula (24) by using
fi ðxÞu^i ðpÞ ¼ u~ðzb ; pÞ;
where vi ¼ ð21ÞN=2þi
^i
maxði;N=2Þ X
kN=2 ð2kÞ! : ðN=2 2 kÞ!k!ðk 2 1Þ!ði 2 kÞ!ð2k 2 iÞ! k¼½ðiþ1Þ=2
for zb [ Gsu ;
ð28Þ
Sutradhar et al. [5] 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 j Þ is total number of the Laplace transform solutions fðp then M £ N: By following essentially the same procedure, the local BIEM can be applied to steady-state or stationary microwave heating problems in FGMs. In this case, the temperature and the its normal derivative are approximated by
i¼1
where u~ðz ; pÞ is the Laplace transform of the prescribed temperature. The numerical integration is a crucial point in meshless methods, because shape functions are rather complex with steep variation. Insufficiently accurate numerical integrations may cause a deterioration and a rank-deficiency in the numerical solutions. To overcome these difficulties the integration domain is divided into smaller partitions. On each partition the standard Gaussian quadrature is applied. The sensitivity of the ‘final’ solution with respect to the number of integration cells and the convergence properties in potential and elasticity theory for the boundary node method are given in Ref. [14]. All integrands in the LBIE method are regular since the collocation point is at the centre of circle on which the integration is applied. The size of subdomain (circle) has not influence the accuracy of results if an accurate numerical integration is used. A smaller size of subdomains may induce larger oscillations in the shape functions derived from meshless interpolations. b
ð29Þ
ð30Þ
for zb [ Gsq ;
n X
N ln 2 X ln 2 i ; vi f t i¼1 t
ð26Þ
(
n X
The time dependent values of the transformed variables can be obtained by an inverse Laplace transform. There are many inversion methods available for Laplace transform. As the inversion of the Laplace transform 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, the Stehfest’s inversion algorithm [15] as a very stable and accurate method [5] is used. Accordingly, an approximate value fa of the inverse f ðtÞ for a specific time t is given by
n X
uðxÞ ¼
fi ðxÞu^i ;
ð31Þ
i¼1 n X ›u ðxÞ ¼ nk ðxÞ fi;k ðxÞu^i : ›n i¼1
ð32Þ
Substitution of Eqs. (31) and (32) into the LBIEs (19) and (20) results in a system of linear algebraic equations for the nodal values of the temperature as n X
( i
b
f ðy Þ þ
ð ›Vs
i¼1
¼
) ~p ›G b i ðx; y Þf ðxÞdG u^i ›n
k;j p b i i ~ ~ G ðx; y Þ f1 ðx; z; uÞf ðxÞ þ ðxÞf;j ðxÞ dV u^i k Vs
n ð X i¼1
þ
ð Vs
~ p ðx;yb ÞdV; u~ÞG Fðx;z;
ð33Þ
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
for yb [ V; and n ð ›G ð ›G X ~p ~p ðx; yÞfi ðxÞdG þ ðx; zb Þ fi ðzb Þ þ limb › n › n y! z G L s s i¼1 ð ~ p ðx; zb ÞdG u^i fi ðxÞdG 2 nk fi;k ðxÞG
Table 1 Critical values of b0 for linear thermal absorptivity in a homogeneous material
Gs
¼
n ð X i¼1
þ
ð Vs
k ~ p ðx; zb Þ f1 ðx; z; u~Þfi ðxÞ þ ;j ðxÞfi;j ðxÞ dV u^i G k Vs
~ p ðx; zb ÞdV; z; u~ÞG Fðx;
ð34Þ
for zb [ G: Here again, for prescribed temperature on Gsu the following collocation formula is applied n X
fi ðxÞu^i ¼ u~ðzb Þ;
for zb [ Gsu ;
7
ð35Þ
i
where u~ðzb Þ is the prescribed temperature on zb [ Gsu :
4. Numerical examples To test the efficiency and the accuracy of the present local BIEM, let us first consider the stationary heat conduction problem. As first example, we consider a square slab being heated with microwave energy [1,16] as shown in Fig. 2. On the whole boundary of the square slab a uniform temperature u ¼ 1 deg is prescribed to analyse hot spots in the slab. In our numerical calculations, a side-length a ¼ 1:0 m and a regular node distribution with 24 boundary and 23 interior nodes are used for the MLS approximation, see Fig. 2. Firstly, a homogeneous material with a linear dependence of the thermal absorptivity on the temperature is considered, i.e. n ¼ 1 in Eq. (16). To obtain the critical value of the thermal absorptivity b0 analytically, the method of Frank-Kamenetskii [17] is used. A comparison of the present numerical results with analytical ones and those numerical results based on the DRBEM [1] is given in Table 1 for three different values of the decay constant.
Fig. 2. Boundary conditions and node distribution for a square slab.
Decay constant
Analytical
Numerical [1]
Present results
g¼0 g¼2 g¼4
19.7 46.4 86.1
19.6 46.0 88.0
19.87 46.7 86.7
Table 1 shows a good agreement of our results with analytical and other numerical results, which confirms the accuracy of the present method. To visualize the dependence of the temperature distribution on the parameters b0 and g; the steady-state temperature profiles on the cross-line x2 ¼ 0:5 are graphically presented in Figs. 3– 6 for several values of b0 and g: Fig. 3 shows that the temperature profile for decay parameter g ¼ 0: is symmetric with respect to the line x1 ¼ 0:5: The dependence of the temperature field on the decay parameter g is shown in Figs. 4 and 5, which imply that with increasing value of the decay parameter the maximum of the temperature profile is shifted to the left side of the square slab, where a higher value of microwave energy is expected. It is caused by the exponential variation of heat sources in the governing equation (3). From Fig. 6 it can be concluded that for a fixed value of the absorptivity parameter b0 ¼ 4 the maximum value of the temperature profile decreases with increasing decay parameter g: To explore the effect of the temperature dependence of the thermal absorptivity on the critical value b0 ; a homogeneous material with a quadratic temperature dependence of the thermal absorptivity is considered in the second example, i.e. n ¼ 2 in Eq. (16). The linearized LBIEs (19) and (20) are solved by an iteration technique where the first iteration value of u~ is used as the prescribed uniform boundary quantity. The tolerance for stopping the iterations was chosen to be 1023. The total number of iterations was less than 7 in all numerical calculations. The present
Fig. 3. Steady-state temperature profiles along x1 -axis at x2 =a ¼ 0:5 for n ¼ 1 and g ¼ 0:
8
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
Fig. 4. Steady-state temperature profiles along x1 -axis at x2 =a ¼ 0:5 for n ¼ 1 and g ¼ 2:
numerical results are compared with analytical results and those obtained by a time marching dual reciprocity method [1] and a Laplace transform dual reciprocity method [16] in Table 3. Analytical results are obtained again by the FrankKamenetskii method [17]. Table 3 shows a quite good agreement of our results with analytical and other numerical results. A comparison of Table 2 with Table 1 reveals that the critical values of b0 are significantly lower for quadratic variation than for linear variation of the thermal absorptivity as presented in Table 1. This means that hot spots for a quadratic temperature dependence of the thermal absorptivity are occurred at lower values of b0 : In the third example, a FGM with a spatial variation of the thermal conductivity is considered. Though the proposed local BIEM in this paper has no restrictions on the spatial variation of the material parameters of the FGM,
Fig. 6. Steady-state temperature profiles along x1 -axis for fixed absorptivity parameter b0 ¼ 4:
the numerical example presented here uses an exponential variation of the material properties with Cartesian coordinates. A uni-directional variation of the thermal conductivity kðxÞ ¼ k0 edx1 ;
ð36Þ
is applied, where k0 ¼ 1:0 W m21 deg21 : Two values of the gradient parameter d ¼ ^1:0 m21 are selected in numerical calculations. In Table 3, the corresponding numerical results for the critical value of the thermal absorptivity parameter b0 in the FGM are compared with critical vales corresponding to a homogeneous material. One can see that the critical value of the absorptivity is substantially increased with increasing gradient parameter for both decay constants considered. For a fixed gradient parameter d; the critical value of the thermal absorptivity increases with increasing decay constant g; which is caused by the fact that an increase of the decay constant leads to a drop of the intensity of the heat source. The temperature profiles at the section x2 ¼ 0:5 are presented in Figs. 7 and 8 for g ¼ 0 and 2, respectively. The gradient parameter d has a similar influence on the temperature profile like the decay constant d: For an exponential graduation of the thermal conductivity along the x1-axis, the maximum value of the temperature is lower than the reference value corresponding to a homogeneous Table 2 Critical values of b0 for quadratic thermal absorptivity in a homogeneous material Decay constant Analytical DRBEM [1] LTDRM [16] Present results
Fig. 5. Steady-state temperature profiles along x1 -axis at x2 =a ¼ 0:5 for n ¼ 1 and g ¼ 4:
g¼0 g¼2 g¼4
4.9 11.4 20.1
4.7 11.0 21.0
4.78 11.27 21.09
4.77 11.2 20.3
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
9
Table 3 Critical values of b0 for thermal absorptivity in FGMs Decay constant
g¼0 g¼2
Gradient parameter
d¼0
d¼1
d ¼ 21
19.99 46.7
39.5 86.2
10.85 25.15
material (i.e. d ¼ 0) as shown in Fig. 7 for a fixed absorptivity parameter b0 ¼ 14 and in Fig. 8 for b0 ¼ 35: In the fourth numerical example, an infinitely long thickwalled hollow cylinder is considered, where the radii R1 ¼ 8:0 m and R2 ¼ 10:0 m are selected. Stationary boundary conditions with a uniform temperature u ¼ 1 deg are prescribed on the whole boundary. Due to the symmetry in geometry and boundary conditions it is sufficient to analyse only a half of the cross-section of the hollow cylinder.
Fig. 9. Boundary conditions and node distribution for a hollow cylinder.
The total number of nodes used for the MLS approximation is 161 with 72 nodes lying on the global boundary as depicted in Fig. 9. A linear temperature dependence of the thermal absorptivity is considered, and an exponential variation of the thermal conductivity on the radius, r; i.e. kðrÞ ¼ k0 edr ;
Fig. 7. Comparison of the temperature profiles in a homogeneous and a FGM square slab for b0 ¼ 14 and g ¼ 0:
ð37Þ
is applied, where k0 ¼ 1:0 W m21 deg21 and d is the gradient parameter. Numerical results for the critical value of the thermal absorptivity b0 in a homogeneous and a FGM cylinder are given in Table 4 for several values of the decay constant. One can observe from Table 4 that a negative value of the gradient parameter leads to a decrease of the critical value of the thermal absorptivity in comparison with the corresponding homogeneous case. For the considered geometry of the hollow cylinder, the value of the thermal conductivity is 0:135k0 on the outer surface while k0 on the inner surface of the hollow cylinder. In the governing equation (3) for the forced heating, the thermal conductivity appears in the dominator, which means that a lower k-value increases the intensity of the heating. Therefore, the critical value of the absorptivity is lower than in the homogeneous case. On the other hand, an increasing value of the decay constant decreases the intensity of the forced heating according to Eq. (3). Therefore, the critical value of the absorptivity b0 increases with increasing decay constant g: The present local BIEM can be also successfully applied to nonstationary heat conduction analysis in a solid under a microwave heating. To demonstrate this, let us consider a homogeneous square slab considered in example 1. On Table 4 Critical values of b0 for thermal absorptivity in a hollow cylinder Decay constant
Fig. 8. Comparison of the temperature profiles in a homogeneous and a FGM square slab for b0 ¼ 35 and g ¼ 2:
g¼0 g ¼ 0:5 g¼1 g¼2
Gradient parameter
d¼0
d ¼ 21
2.5 4.3 7.7 17.1
0.95 1.62 2.87 6.2
10
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
at the centre of the slab (i.e. x1 ¼ 0:5; x2 ¼ 0:5;) and for a fixed value b0 ¼ 14: It shows here that the maximum or the stationary temperature decreases with increasing decay constant g as expected.
5. Conclusions
Fig. 10. Time variations of the temperature at the section x2 =a ¼ 0:5 and for b0 ¼ 64 and g ¼ 4:
the whole boundary of the square slab a thermal shock of the form uðx; tÞ ¼ T·HðtÞ is applied, where T ¼ 1 deg is the amplitude of the thermal shock, and HðtÞ is the Heaviside step function. The initial condition u ¼ 1 deg is prescribed at t ¼ 0: A linear temperature dependence of the thermal absorptivity is assumed, and the Stehfest’s numerical inversion method [15] of the Laplace transform is applied. For fixed values b0 ¼ 64; g ¼ 4 and at three different positions x1 =a ¼ 0:25; 0.5, 0.75 at the middle section x2 ¼ 0:5 of the square slab, numerical results for the temperature are presented in Fig. 10 versus time. Fig. 10 shows that the temperature increases sharply with increasing time in the initial stage, and after a short time it tends to its stationary or steady-state value. The effect of the decay constant g on the temporal variation of the temperature is presented in Fig. 11. Here, numerical results are computed
A local BIEM in conjunction with the MLS approximation of physical fields is presented in this paper for steady-state and transient heat conduction problems with a nonlinear heat source in FGMs. A power-law and an exponential decay of the heat source intensity are used for modelling the microwave heating. The method can be applied to a general spatial variation of the material constants for describing FGMs. The initial-boundary value problem for the transient heat analysis is solved in the Laplace transform domain with a subsequent numerical Laplace inversion to obtain time-dependent solutions. The LBIEs are formulated in a boundary-domain integral form with a simple fundamental solution corresponding to the Poisson’s equation. The boundary-domain integral formulation can be easily implemented by choosing a simple circular subdomain to which the LBIE is applied. Contrary to the conventional BIEM or BEM, all integrands in the present LBIEM formulation are regular. Thus, no special integration techniques are required to evaluate the arising integrals. The present LBIEM or LBEM removes the wellknown restriction of the conventional BIEM/BEM to problems with homogeneous material properties. Numerical examples for both steady-state and transient microwave heating problems in FGMs are given to demonstrate the efficiency and the accuracy of the present LBIEM. Numerical results show the effects of the temperature dependence of the thermal absorptivity, the decay constant of the microwave heating source, and the material gradient on the occurrence of hot spots and the temperature distribution in FGMs.
Acknowledgements The authors acknowledge the support by the Slovak Science and Technology Assistance Agency registered under number APVT-51-003702, as well as by the Slovak Grant Agency VEGA-2303823, 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. 11. Time variations of the temperature at the centre of a square slab and for b0 ¼ 14:
[1] Zhu S, Zhang Y, Marchant TR. A DRBEM model for microwave heating problems. Appl Math Modelling 1995;19:287 –97.
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11 [2] Miyamoto Y, Kaysser WA, Rabin BH, Kawasaki A, Ford RG. Functionally graded materials: design processing and applications. Dordrecht: Kluwer Academic Publishers; 1999. [3] Suresh S, Mortensen A. Fundamentals of functionally graded materials. London: Institute of Materials; 1998. [4] Butterfield R. An application of the boundary element method to potential flow problems in generally inhomogeneous bodies. In: Brebbia CA, editor. Recent advances in boundary element method. London: Pentech Press; 1978. [5] Sutradhar A, Paulino GH, Gray LJ. Transient heat conduction in homogeneous and non-homogeneous materials by the Laplace transform Galerkin boundary element method. Engng Anal Bound Elem 2002;26:119– 32. [6] Zhu T, Zhang JD, Atluri SN. A local boundary integral equation (LBIE) method in computational mechanics, and a meshless discretization approach. Comput Mech 1998;21:223– 35. [7] De S, Bathe KJ. The method of finite spheres. Comput Mech 2000;25: 329–45. [8] Sladek J, Sladek V, Atluri SN. Local boundary integral equation (LBIE) method for solving problems of elasticity with nonhomogeneous material properties. Comput Mech 2000;24: 456 – 62. [9] Sladek J, Sladek V, Zhang Ch. Transient heat conduction analysis in functionally graded materials by the meshless local
[10]
[11] [12]
[13]
[14]
[15] [16]
[17]
11
boundary integral equation method. Submitted for publication; 2002. Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: an overview and recent developments. Comput Meth Appl Mech Engng 1996;139:3–47. Chen HT, Lin JY. Application of Laplace transform to nonlinear transient problems. Appl Math Modelling 1992;15:144– 51. Ramachandran PA. Application of the boundary element method to nonlinear diffusion with reaction problems. Int J Numer Meth Engng 1990;29:1021 –31. Sladek V, Sladek J, Atluri SN, Van Keer R. Numerical integration of singularities in meshless implementation of local boundary integral equations. Comput Mech 2000;25:394 –403. Chati MK, Paulino GH, Mukherjee S. The meshless standard and hypersingular boundary node methods: applications to error estimation and adaptivity in three-dimensional problems. Int J Numer Meth Engng 2001;50:2233–69. Stehfest H. Algorithm 368: numerical inversion of Laplace transform. Comm Assoc Comput Mach 1970;13:47–9. Zhu S, Satravaha P. An efficient computational method for modelling transient heat conduction with nonlinear source terms. Appl Math Modelling 1996;20:513– 22. Frank-Kamenetskii DA. Diffusion and heat exchange in chemical kinetics. Princeton: Princeton University Press; 1955.
A local BIEM for analysis of transient heat conduction with nonlinear source terms in FGMs Jan Sladeka,*, Vladimir Sladeka, Ch. Zhangb b
a Institute of Construction and Architecture, Slovak Academy of Sciences, 84220 Bratislava, Slovakia Department of Civil Engineering, University of Applied Sciences Zittau/Go¨rlitz, D-02763 Zittau, Germany
Received 6 February 2003; revised 13 May 2003; accepted 14 May 2003
Abstract The diffusion equation with nonlinear heat source intensity in functionally graded materials (FGMs) is considered. In FGMs the thermal material properties are dependent on spatial coordinates. For transient or steady-state heat problems in FGMs the conventional boundary integral equation method or boundary element method cannot be applied due to the lack of a fundamental solution. In this paper, a local boundary integral equation method is proposed to analyse a temperature distribution in a nonhomogeneous body under a microwave heating. To eliminate time variable in the heat equation, the Laplace transform technique is used. The boundary-domain integral formulation with a simple fundamental solution corresponding to the Laplace operator is related to all subdomains which cover the analysed domain. If such integral equations are considered on small subdomains with a simple geometry (circle), domain integrals can be easily evaluated. Physical fields (temperature, heat flux) on the local boundary and in the interior of the subdomain are approximated by the moving least-square. The method is completely element free. q 2003 Elsevier Ltd. All rights reserved. Keywords: Hot spots; Transient heat conduction analysis; Local boundary integral equations; Functionally graded materials; Laplace transform; Moving leastsquare approximation; Absorbtivity
1. Introduction Transient heat conduction equations with nonlinear source terms arise as governing equations in many different areas of mathematical physics. Examples in the literature include microwave heating process, mass transport in groundwater, and ignition. The use of microwave heating has attracted much attention in industrial processes due to its high efficiency. A crucial point in microwave heating is to analyse hot spots, regions where temperature is enormously running. This phenomenon is occurred due to dependence of material properties on temperature. Namely, thermal absorptivity increases with temperature. Hot spots out of control can damage the sample in an industrial process. The first attempt to analyse hot spots in a homogeneous body is given by Zhu et al. [1]. The dual reciprocity boundary element method (BEM) is proposed there. In the DRBEM, some interior nodes have to be supplemented to boundary * Corresponding author. Tel.: þ 421-7-54788662; fax: þ 421-754772494. E-mail address: [email protected] (J. Sladek). 0955-7997/04/$ - see front matter q 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0955-7997(03)00093-6
nodes for an adequate spatial approximation of physical fields in transforming domain integrals into boundary integrals. Functionally graded materials (FGMs) are preferred and favoured in many engineering structures and components due to their excellent thermal properties [2,3]. In FGMs the volume fractions of composite constituents are varying continuously in space. Solution of boundary or initial boundary value problems for FGMs requires special numerical methods due to the high mathematical complexity caused by the material nonhomogeneity. Beside the wellestablished finite element method (FEM), the boundary integral equation method (BIEM) or BEM provides an efficient and popular alternative to the FEM for solving certain class of boundary or initial boundary value problems. However, the pure BIEM/BEM formulation can be efficiently applied only to problems where fundamental solutions or Green’s functions are available in simple forms. Unfortunately, for general FGMs the corresponding fundamental solutions are either not available or they are mathematically too complicated. To avoid this difficulty,
2
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
one possibility to obtain a BEM formulation is the use of fundamental solutions for a fictitious homogeneous medium, as has been suggested for the first time by Butterfield [4] for potential flow problems. This approach, however, leads to a boundary-domain integral formulation with additional domain integrals representing the effects of the gradients of primary fields to obtain a unique formulation. Some computational difficulties in the numerical implementation of the boundary-domain integral method have been reported. A BEM based on Green’s function method for three-dimensional transient heat conduction analysis in FGMs is presented in Ref. [5], where the method is restricted to exponentially varying material properties. An alternative way is recently suggested in Refs. [6 – 9] by using the local BIEM. In this method, the analysed domain is divided into small subdomains with a simple geometry such as circles and spheres. In each sufficiently small subdomain one can apply the fundamental solutions corresponding to a homogeneous material with adequate local material constants [6 – 9]. On the boundary of subdomains the local boundary integral equations (LBIEs) are applied. The temperature and the heat flux on the local boundary and in the interior of the subdomain are approximated by the moving least-square (MLS) scheme. The key idea of MLS interpolants is that it is only necessary to construct an array of nodes in the domain under consideration [10]. Thus, the method is completely element-free. In this paper, a local BIEM for microwave heating analysis in FGMs is presented. The transient heat conduction problem is described by a parabolic partial differential equation, while its steady-state is governed by an elliptic one. To eliminate the time dependence of the differential equation for transient heat conduction analysis, 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 yet not known to our knowledge for general FGMs. 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 integral formulation in the global BIEM/BEM is numerically not efficient, we apply these integral equations to small subdomains, which result in a LBIE formulation. By choosing a simple geometric shape such as a circle for the subdomains, the arising domain integrals in the local BIEM/BEM can be easily evaluated. This means that the boundary-domain integral formulation for simple subdomains does not bring any numerical difficulties. To describe the material gradation of FGMs, an exponential spatial variation of thermal material properties is applied in numerical calculations. The temperature dependence of the thermal absorptivity is approximated by a power law. Numerical examples are presented for both steady-state and transient microwave heating problems.
2. Local boundary integral equations Microwave heating of a continuously nonhomogeneous body can be modelled by the nonstationary forced heat equation [1]
u;ii ðx; tÞ 2 ¼2
k 1 ›u ðx; tÞ þ ;i ðxÞu;i ðx; tÞ kðxÞ ›t k
1 bðuÞlEl2 kðxÞ
ð1Þ
and Maxwell’s equations, as governing equations for diffusion of heat and propagation of electromagnetic field including its absorption. In Eq. (1) t is the temporal variable; uðx; tÞ; the temperature field; k and k; the thermal conductivity and diffusivity, respectively; bðuÞ; the thermal absorptivity; lEl is the amplitude of the electric field. For simplicity the electric conductivity, magnetic and electric permeabilities are considered to be constant. Only thermal conductivity and absorptivity are considered to be dependent on temperature. The amplitude of the electric field is exponentially dependent on spatial variables [1] lEl ¼ expð2gz=2Þ;
ð2Þ
for decay from an incident boundary at z ¼ 0; where g is the decay constant, and thus Eq. (1) is uncoupled from Maxwell’s equations. Generally g should be complex propagation constant. In such a case the coupled system of both equations has to be solved. The simplified model for microwave heating is then determined by
u;ii ðx; tÞ 2 ¼2
1 ›u k ðx; tÞ þ ;i ðxÞu;i ðx; tÞ kðxÞ ›t k
1 bðuÞexpð2gzÞ: kðxÞ
ð3Þ
The microwave heating problem is described by a nonlinear parabolic partial differential equation (3). Analytical solutions exist only for simple geometrical and loading conditions, with uniform material properties and g ¼ 0: For general cases, Eq. (3) has to be solved numerically. A dual reciprocity BEM has been presented in Ref. [1], where numerical results are given for uniform material properties. If the thermal absorptivity is a nonlinear function of temperature, a linearization scheme has to be applied to Eq. (3). This can be done by using the following first-order Taylor-series expansion [11,12] ›b ~ ~ bðuÞ ¼ bðuÞ þ ðu 2 uÞ ; ð4Þ ›u u¼u~ where u~ is the previously iterated solution at a particular time. To eliminate the time variable in the differential equation (3), the Laplace transform technique is used.
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
Applying the Laplace transform ð1 L½uðx; tÞ ¼ uðx; pÞ ¼ uðx; tÞe2pt dt;
the evaluation of hypersingular integrals make the global boundary-domain formulation somewhat cumbersome. Another approach is the application of the LBIEs to the boundaries of small regular no-overlapping subdomains randomly spread in the analysed domain. If we consider a subdomain Vs instead of the entire domain V; the following LBIE should also valid over the subdomain
0
to the linearized governing equation (3) we obtain k p u;ii ðx; pÞ 2 uðx; pÞ þ ;i ðxÞu;i ðx; pÞ þ f1 ðx; z; u~Þuðx; pÞ k k pÞ; ¼ 2Fðx; where p is the Laplace transform parameter 1 ›b expð2gzÞ f1 ðx; z; u~Þ ¼ kðxÞ ›u u¼u~
ð5Þ
ð6Þ
and z; u~Þ ¼ 1 bðu~Þ 2 u~ ›b u¼u~ expð2gzÞ þ 1 uðx;0Þ; Fðx; kðxÞ ›u kðxÞ ð7Þ is the redefined body heat source in the Laplace transform domain with the initial condition for the temperature field. The solution of the governing equation (5) can be found in a weak form with a fundamental solution or Green’s function Gðx; yÞ ¼
1 1 ln ; 2p r
ð8Þ
corresponding to the Poisson’s equation G;ii ðx; yÞ þ dðx; yÞ ¼ 0;
3
ð9Þ
where dðx; yÞ is the Dirac delta function and r is the distance between the field and the source points, i.e. r ¼ lx 2 yl: By the use of the Gauss divergence theorem to the weak form of Eq. (5) we obtain an integral representation for the temperature field in the Laplace transform domain ð ›u ð ›G ðx; pÞGðx; yÞdG 2 ðx; yÞdG uðy; pÞ ¼ uðx; pÞ ›n G ›n G ð p 2 f1 ðx; z; u~Þ uðx; pÞGðx; yÞdV 2 V kðxÞ ð k ;i ðxÞu;i ðx; pÞGðx; yÞdV þ V k ð z; u~ÞGðx; yÞdV; þ ð10Þ Fðx; V
where G is the boundary of the analysed domain V: Since the boundary-domain integral representation (10) for the Laplace transform of the temperature contains temperature gradients, an 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 a unique global boundary-domain formulation is obtained. Domain discretization and
uðy; pÞ ¼
ð
ð ›u ›G ðx;pÞGðx; yÞdG 2 uðx;pÞ ðx; yÞdG ›n ›Vs ›n ›Vs ð p 2 f1 ðx; z; u~Þ uðx; pÞGðx;yÞdV 2 Vs kðxÞ ð k ;i ðxÞu;i ðx; pÞGðx;yÞdV þ Vs k ð þ u~ÞGðx;yÞdV: ð11Þ Fðx;z; Vs
The integral equation (11) is considered on small subdomains Vs , V: Hence, no of the boundary densities in the Laplace transformed domain are prescribed on the local boundary ›Vs as long as it lies entirely inside Vs : Both quantities (temperature and heat flux) are unknowns on ›Vs : To reduce the number of unknowns on ›Vs the concept of a companion solution can be utilized successfully [6]. The companion solution is associated with the fundamental solution Gðx; yÞ and is defined as the solution to the following equations on V0s ;
~ ;ii ðx;yÞ ¼ 0 G
~ ;ii ðx;yÞ ¼ Gðx; yÞ G
on ›V0s ;
ð12Þ
where V0s and ›V0s are depicted in Fig. 1. For simplicity, V0s is taken as a circle in the present implementation. The modified test function ~ p ðx; yÞ ¼ Gðx; yÞ 2 Gðx; ~ yÞ ¼ 1 ln r0 G 2p r has to satisfy the governing equation (9), and it is the modified fundamental solution or Green’s function vanishing on the boundary of the circular subdomain of radius r0 : Since the integral equation (11) is also valid for the modified ~ p ðx; yÞ; we obtain fundamental solution G ð
~p ›G ðx; yÞdG uðx; pÞ ›n ›Vs ð p ~ p ðx; yÞdV 2 f1 ðx; z; u~Þ uðx; pÞG 2 Vs kðxÞ ð k ;i ~ p ðx; yÞdV ðxÞu;i ðx; pÞG þ Vs k ð ~ p ðx; yÞdV; z; u~ÞG þ ð13Þ Fðx;
uðy; pÞ ¼ 2
Vs
4
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
Fig. 1. Local boundaries and the domain of definition of the MLS approximation for the trial function.
for interior source points y [ V and
uðz; pÞ ¼ 2
ð
described by a power law of the form
ð ~p ›G ðx; zÞdG 2 lim uðx; pÞ uðx; pÞ y!z Gs ›n Ls
bðuÞ ¼ b0 un :
ð ›u ~p ›G ~ p ðx; zÞdG ðx; yÞdG þ ðx; pÞG ›n G s ›n ð p ~ p ðx; zÞdV 2 f1 ðx; z; u~Þ uðx; pÞG 2 Vs kðxÞ ð k ð ;i ~ p ðx; zÞdV þ z; u~Þ ðxÞu;i ðx; pÞG þ Fðx; Vs k Vs ~ p ðx; zÞdV; G
ð14Þ
for source points located on the global boundary z [ Gs , G: In the latter case, ›Vs ¼ Ls < Gs ; where Ls is the circular part of ›Vs : The strongly singular integral in the integral equation (14) is written in the limit form. One can see that only temperature is unknown in the local boundary integral equation (13) at the interior source point. The heat flux is eliminated there due to convenient behaviour of the Green ~ p ðx; yÞ on the circle. For the source point on the function G global boundary the LBIE (14) contains both quantities (temperature and heat flux), however one of them is prescribed in a well posed boundary value problem. Next, we consider the steady-state or stationary heat conduction problem without time-dependence in the nonlinear parabolic partial differential equation (3). In this case, the governing equation (3) can be simplified to k 1 u;ii ðxÞ þ ;i ðxÞu;i ðxÞ ¼ 2 bðuÞexpð2gzÞ: k kðxÞ
ð15Þ
The steady-state heat conduction equation (3) is an elliptic partial differential equation and it is suitable for describing hot spots problems. Hot spots are the localized areas of high temperature caused by the graduation of the thermal absorptivity with temperature. To simplify the analysis, we assume that the thermal absorptivity is
ð16Þ
Note here that many engineering materials can be adequately approximated by the power law (16). According to Eq. (16), a thermal runway is possible only if n $ 1 [1]. Substitution of Eq. (16) into Eq. (15) leads to a eigen-value problem for the determination of the critical value of the absorptivity coefficient b0 ; which can be written as
u;ii ðxÞ þ
k;i 1 ðxÞu;i ðxÞ þ b ½uðxÞn expð2gzÞ ¼ 0: k kðxÞ 0
ð17Þ
By using the first-order Taylor-series expansion (4), Eq. (17) can be written in a linearized form as
u;ii ðxÞ þ
k;i ðxÞu;i ðxÞ þ f1 ðx; z; u~ÞuðxÞ ¼ 2Fðx; z; u~Þ; k
ð18Þ
where f1 ðx; z; u~Þ ¼
1 ›b 1 expð2gzÞ ¼ expð2gzÞb0 nu~n21 ; ~ kðxÞ ›u u¼u kðxÞ
Fðx; z; u~Þ ¼
1 ~ ›b bðuÞ 2 u~ u¼u~ expð2gzÞ kðxÞ ›u
¼
1 expð2gzÞb0 ð1 2 nÞu~n ; kðxÞ
in which u~ is the previously iterated solution for a nonlinear variation of the thermal absorptivity with temperature. Note here that no iteration is needed if the coefficient n ¼ 0 and 1, i.e. when the thermal absorptivity is independent of temperature or linearly dependent on temperature. The solution of the eigen-value problem (18) can be given in an integral form through the LBIEs. The procedure to obtain the LBIEs is the same as in the nonstationary case. By omitting integrals representing the nonstationary term in the BIEs (13) and (14), a system of LBIEs at interior and
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
boundary nodes can be obtained as ð ð ~p ›G ~ p ðx;yÞdV ðx;yÞdG þ uðyÞ¼2 uðxÞ f1 ðx;z; u~ÞuðxÞG ›n ›Vs Vs ð k ð ;i ~ p ðx;yÞdV þ ~ p ðx;yÞdV ðxÞu;i ðxÞG þ Fðx;z; u~ÞG Vs k Vs ð19Þ and
given by a fitting curve obtained by the least square method from the nodal values. The stationarity of J with respect to aðxÞ leads to the following linear relation
›J=›a ¼ AðxÞaðxÞ 2 BðxÞu^ ¼ 0; where matrices AðxÞ and BðxÞ are defined AðxÞ ¼ PT WP ¼
n X
wi ðxÞpðxi ÞpT ðxi Þ;
i¼1
ð
ð ~p ~p ›G ›G ðx; zÞdG 2lim ðx;yÞdG uðxÞ y!z Gs ›n ›n Ls ð ›u ð ~ p ðx; zÞdG 2 ~ p ðx; zÞdV ðxÞG þ f1 ðx;z; u~ÞuðxÞG G s ›n Vs ð k ð ;i ~ p ðx; zÞdV þ ~ p ðx; zÞdV: ðxÞu;i ðxÞG þ Fðx;z; u~ÞG Vs k Vs ð20Þ
uðzÞ¼2
uðxÞ
Here again, the strongly singular integral in Eq. (20) is written in the limit form, which is convenient for the numerical evaluation of the strongly singular integral when the approximation of the boundary density is not known in a closed form [13]. The remaining problem is how to represent the values of the unknown function on the local boundary and in the interior of a subdomain. This will be explained and discussed in Section 3.
3. The MLS approximation 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. To explain the key idea of the MLS approximation, let us consider now a subdomain Vx of the analysed domain V: The subdomain is the neighbourhood of a point x and denoted as the domain of definition of the MLS approximation for the trial function at x: To approximate the distribution of the temperature field uðxÞ in Vs over a number of randomly located nodes {xi } ði ¼ 1; 2; …; nÞ; the moving least square approximant uh ðxÞ of u for ;x [ Vs is defined as
uh ðxÞ ¼ pT ðxÞaðxÞ; T
1
;x [ Vs ; 2
ð21Þ
n X
BðxÞ ¼ PT W ¼ ½w1 ðxÞpðx1 Þ; w2 ðxÞpðx2 Þ; …; wn ðxÞpðxn Þ: The approximated function can be written as [10]
uh ðxÞ ¼ FT ðxÞ·u^ ¼
wi ðxÞ½pT ðxi ÞaðxÞ 2 u^i 2 ;
i¼1
where wi ðxÞ is the weight function associated with the node i and u^i are certain nodal values. Approximation values are
n X
fi ðxÞu^i ;
ð22Þ
i¼1
where
fi ðxÞ ¼
m X
pj ðxÞ½A21 ðxÞBðxÞji ;
j¼1
is the shape function associated with node i: Since the shape functions do not exhibit the Kronecker-delta property fi ðxj Þ dij ; the nodal values are fictitious nodal values which are different from real nodal values. The number of nodes n used for the approximation of uðxÞ is prescribed by the weight function wi ðxÞ: The Gaussian weight function is applied in the present work 8 > exp½2ðdi =ci Þ2 2 exp½2ðr i =ci Þ2 > < ; 0 # di # ri 1 2 exp½2ðr i =ci Þ2 wi ðxÞ ¼ ; > > : 0; di $ ri ð23Þ where di ¼ kx 2 xi k is the distance between the integration point and the source point, ci is a constant controlling the shape of the weight function wi ðxÞ; and r i is the size of the support domain. For a sharper Gaussian weight function the number of nodes used for approximation is smaller than in an opposite case. A necessary condition for a well-defined MLS approximation is that at least m weight functions are nonzero (i.e. n $ m) due to the regularity of the matrix A. Finally, the approximation formulae for the Laplace transform of the temperature and its normal derivative can be written as
m
where p ðxÞ ¼ ½p ðxÞ; p ðxÞ; …; p ð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 weighted discrete L2 norm, defined as JðxÞ ¼
5
uðx; pÞ ¼
n X
fi ðxÞu^i ðpÞ;
ð24Þ
i¼1 n X ›u ðx; pÞ ¼ nk ðxÞ u^i ðpÞfi;k ðxÞ: ›n i¼1
ð25Þ
Eqs. (24) and (25) imply that both quantities u and ›u=›n can be approximated by the fictitious nodal values u^i ðpÞ: Substituting the approximation formulae (24) and (25) into the local boundary integral Eqs. (13) and (14) we obtain a set
6
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
of linear algebraic equations ( ) n ð X ~p ›G i b b i ðx; y Þf ðxÞdG u^i ðpÞ f ðy Þ þ › n › V s i¼1 ¼2
n ð X i¼1
~ p ðx; yb Þ p fi ðxÞ 2 f1 ðx; z; u~Þfi ðxÞ G kðxÞ Vs
ð k;j ~ p ðx; yb ÞdV; z; u~ÞG ðxÞfi;j ðxÞ dVu^i ðpÞ þ Fðx; k Vs
2
for yb [ V and n X i¼1
fa ðtÞ ¼
fi ðzb Þ þ limb
ð
y!z
i
f ðxÞdG 2
Gs
ð
i¼1
2
ð ›G ~p ~p ›G ðx; yÞfi ðxÞdG þ ðx; zb Þ ›n Ls › n )
ð Gs
¼2
~ p ðx; zb ÞdG nk fi;k ðxÞG
~ p ðx; zb Þ G Vs
u ðpÞ
p i f ðxÞ 2 f1 ðx; z; u~Þfi ðxÞ kðxÞ
ð k;j ~ p ðx; zb ÞdV; z; u~ÞG ðxÞfi;j ðxÞ dVu^i ðpÞ þ Fðx; k Vs ð27Þ
where Gsq is a part of the global boundary G; over which the heat flux is prescribed. If the temperature is prescribed on a part Gsu we propose to collocate the approximation formula (24) by using
fi ðxÞu^i ðpÞ ¼ u~ðzb ; pÞ;
where vi ¼ ð21ÞN=2þi
^i
maxði;N=2Þ X
kN=2 ð2kÞ! : ðN=2 2 kÞ!k!ðk 2 1Þ!ði 2 kÞ!ð2k 2 iÞ! k¼½ðiþ1Þ=2
for zb [ Gsu ;
ð28Þ
Sutradhar et al. [5] 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 j Þ is total number of the Laplace transform solutions fðp then M £ N: By following essentially the same procedure, the local BIEM can be applied to steady-state or stationary microwave heating problems in FGMs. In this case, the temperature and the its normal derivative are approximated by
i¼1
where u~ðz ; pÞ is the Laplace transform of the prescribed temperature. The numerical integration is a crucial point in meshless methods, because shape functions are rather complex with steep variation. Insufficiently accurate numerical integrations may cause a deterioration and a rank-deficiency in the numerical solutions. To overcome these difficulties the integration domain is divided into smaller partitions. On each partition the standard Gaussian quadrature is applied. The sensitivity of the ‘final’ solution with respect to the number of integration cells and the convergence properties in potential and elasticity theory for the boundary node method are given in Ref. [14]. All integrands in the LBIE method are regular since the collocation point is at the centre of circle on which the integration is applied. The size of subdomain (circle) has not influence the accuracy of results if an accurate numerical integration is used. A smaller size of subdomains may induce larger oscillations in the shape functions derived from meshless interpolations. b
ð29Þ
ð30Þ
for zb [ Gsq ;
n X
N ln 2 X ln 2 i ; vi f t i¼1 t
ð26Þ
(
n X
The time dependent values of the transformed variables can be obtained by an inverse Laplace transform. There are many inversion methods available for Laplace transform. As the inversion of the Laplace transform 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, the Stehfest’s inversion algorithm [15] as a very stable and accurate method [5] is used. Accordingly, an approximate value fa of the inverse f ðtÞ for a specific time t is given by
n X
uðxÞ ¼
fi ðxÞu^i ;
ð31Þ
i¼1 n X ›u ðxÞ ¼ nk ðxÞ fi;k ðxÞu^i : ›n i¼1
ð32Þ
Substitution of Eqs. (31) and (32) into the LBIEs (19) and (20) results in a system of linear algebraic equations for the nodal values of the temperature as n X
( i
b
f ðy Þ þ
ð ›Vs
i¼1
¼
) ~p ›G b i ðx; y Þf ðxÞdG u^i ›n
k;j p b i i ~ ~ G ðx; y Þ f1 ðx; z; uÞf ðxÞ þ ðxÞf;j ðxÞ dV u^i k Vs
n ð X i¼1
þ
ð Vs
~ p ðx;yb ÞdV; u~ÞG Fðx;z;
ð33Þ
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
for yb [ V; and n ð ›G ð ›G X ~p ~p ðx; yÞfi ðxÞdG þ ðx; zb Þ fi ðzb Þ þ limb › n › n y! z G L s s i¼1 ð ~ p ðx; zb ÞdG u^i fi ðxÞdG 2 nk fi;k ðxÞG
Table 1 Critical values of b0 for linear thermal absorptivity in a homogeneous material
Gs
¼
n ð X i¼1
þ
ð Vs
k ~ p ðx; zb Þ f1 ðx; z; u~Þfi ðxÞ þ ;j ðxÞfi;j ðxÞ dV u^i G k Vs
~ p ðx; zb ÞdV; z; u~ÞG Fðx;
ð34Þ
for zb [ G: Here again, for prescribed temperature on Gsu the following collocation formula is applied n X
fi ðxÞu^i ¼ u~ðzb Þ;
for zb [ Gsu ;
7
ð35Þ
i
where u~ðzb Þ is the prescribed temperature on zb [ Gsu :
4. Numerical examples To test the efficiency and the accuracy of the present local BIEM, let us first consider the stationary heat conduction problem. As first example, we consider a square slab being heated with microwave energy [1,16] as shown in Fig. 2. On the whole boundary of the square slab a uniform temperature u ¼ 1 deg is prescribed to analyse hot spots in the slab. In our numerical calculations, a side-length a ¼ 1:0 m and a regular node distribution with 24 boundary and 23 interior nodes are used for the MLS approximation, see Fig. 2. Firstly, a homogeneous material with a linear dependence of the thermal absorptivity on the temperature is considered, i.e. n ¼ 1 in Eq. (16). To obtain the critical value of the thermal absorptivity b0 analytically, the method of Frank-Kamenetskii [17] is used. A comparison of the present numerical results with analytical ones and those numerical results based on the DRBEM [1] is given in Table 1 for three different values of the decay constant.
Fig. 2. Boundary conditions and node distribution for a square slab.
Decay constant
Analytical
Numerical [1]
Present results
g¼0 g¼2 g¼4
19.7 46.4 86.1
19.6 46.0 88.0
19.87 46.7 86.7
Table 1 shows a good agreement of our results with analytical and other numerical results, which confirms the accuracy of the present method. To visualize the dependence of the temperature distribution on the parameters b0 and g; the steady-state temperature profiles on the cross-line x2 ¼ 0:5 are graphically presented in Figs. 3– 6 for several values of b0 and g: Fig. 3 shows that the temperature profile for decay parameter g ¼ 0: is symmetric with respect to the line x1 ¼ 0:5: The dependence of the temperature field on the decay parameter g is shown in Figs. 4 and 5, which imply that with increasing value of the decay parameter the maximum of the temperature profile is shifted to the left side of the square slab, where a higher value of microwave energy is expected. It is caused by the exponential variation of heat sources in the governing equation (3). From Fig. 6 it can be concluded that for a fixed value of the absorptivity parameter b0 ¼ 4 the maximum value of the temperature profile decreases with increasing decay parameter g: To explore the effect of the temperature dependence of the thermal absorptivity on the critical value b0 ; a homogeneous material with a quadratic temperature dependence of the thermal absorptivity is considered in the second example, i.e. n ¼ 2 in Eq. (16). The linearized LBIEs (19) and (20) are solved by an iteration technique where the first iteration value of u~ is used as the prescribed uniform boundary quantity. The tolerance for stopping the iterations was chosen to be 1023. The total number of iterations was less than 7 in all numerical calculations. The present
Fig. 3. Steady-state temperature profiles along x1 -axis at x2 =a ¼ 0:5 for n ¼ 1 and g ¼ 0:
8
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
Fig. 4. Steady-state temperature profiles along x1 -axis at x2 =a ¼ 0:5 for n ¼ 1 and g ¼ 2:
numerical results are compared with analytical results and those obtained by a time marching dual reciprocity method [1] and a Laplace transform dual reciprocity method [16] in Table 3. Analytical results are obtained again by the FrankKamenetskii method [17]. Table 3 shows a quite good agreement of our results with analytical and other numerical results. A comparison of Table 2 with Table 1 reveals that the critical values of b0 are significantly lower for quadratic variation than for linear variation of the thermal absorptivity as presented in Table 1. This means that hot spots for a quadratic temperature dependence of the thermal absorptivity are occurred at lower values of b0 : In the third example, a FGM with a spatial variation of the thermal conductivity is considered. Though the proposed local BIEM in this paper has no restrictions on the spatial variation of the material parameters of the FGM,
Fig. 6. Steady-state temperature profiles along x1 -axis for fixed absorptivity parameter b0 ¼ 4:
the numerical example presented here uses an exponential variation of the material properties with Cartesian coordinates. A uni-directional variation of the thermal conductivity kðxÞ ¼ k0 edx1 ;
ð36Þ
is applied, where k0 ¼ 1:0 W m21 deg21 : Two values of the gradient parameter d ¼ ^1:0 m21 are selected in numerical calculations. In Table 3, the corresponding numerical results for the critical value of the thermal absorptivity parameter b0 in the FGM are compared with critical vales corresponding to a homogeneous material. One can see that the critical value of the absorptivity is substantially increased with increasing gradient parameter for both decay constants considered. For a fixed gradient parameter d; the critical value of the thermal absorptivity increases with increasing decay constant g; which is caused by the fact that an increase of the decay constant leads to a drop of the intensity of the heat source. The temperature profiles at the section x2 ¼ 0:5 are presented in Figs. 7 and 8 for g ¼ 0 and 2, respectively. The gradient parameter d has a similar influence on the temperature profile like the decay constant d: For an exponential graduation of the thermal conductivity along the x1-axis, the maximum value of the temperature is lower than the reference value corresponding to a homogeneous Table 2 Critical values of b0 for quadratic thermal absorptivity in a homogeneous material Decay constant Analytical DRBEM [1] LTDRM [16] Present results
Fig. 5. Steady-state temperature profiles along x1 -axis at x2 =a ¼ 0:5 for n ¼ 1 and g ¼ 4:
g¼0 g¼2 g¼4
4.9 11.4 20.1
4.7 11.0 21.0
4.78 11.27 21.09
4.77 11.2 20.3
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
9
Table 3 Critical values of b0 for thermal absorptivity in FGMs Decay constant
g¼0 g¼2
Gradient parameter
d¼0
d¼1
d ¼ 21
19.99 46.7
39.5 86.2
10.85 25.15
material (i.e. d ¼ 0) as shown in Fig. 7 for a fixed absorptivity parameter b0 ¼ 14 and in Fig. 8 for b0 ¼ 35: In the fourth numerical example, an infinitely long thickwalled hollow cylinder is considered, where the radii R1 ¼ 8:0 m and R2 ¼ 10:0 m are selected. Stationary boundary conditions with a uniform temperature u ¼ 1 deg are prescribed on the whole boundary. Due to the symmetry in geometry and boundary conditions it is sufficient to analyse only a half of the cross-section of the hollow cylinder.
Fig. 9. Boundary conditions and node distribution for a hollow cylinder.
The total number of nodes used for the MLS approximation is 161 with 72 nodes lying on the global boundary as depicted in Fig. 9. A linear temperature dependence of the thermal absorptivity is considered, and an exponential variation of the thermal conductivity on the radius, r; i.e. kðrÞ ¼ k0 edr ;
Fig. 7. Comparison of the temperature profiles in a homogeneous and a FGM square slab for b0 ¼ 14 and g ¼ 0:
ð37Þ
is applied, where k0 ¼ 1:0 W m21 deg21 and d is the gradient parameter. Numerical results for the critical value of the thermal absorptivity b0 in a homogeneous and a FGM cylinder are given in Table 4 for several values of the decay constant. One can observe from Table 4 that a negative value of the gradient parameter leads to a decrease of the critical value of the thermal absorptivity in comparison with the corresponding homogeneous case. For the considered geometry of the hollow cylinder, the value of the thermal conductivity is 0:135k0 on the outer surface while k0 on the inner surface of the hollow cylinder. In the governing equation (3) for the forced heating, the thermal conductivity appears in the dominator, which means that a lower k-value increases the intensity of the heating. Therefore, the critical value of the absorptivity is lower than in the homogeneous case. On the other hand, an increasing value of the decay constant decreases the intensity of the forced heating according to Eq. (3). Therefore, the critical value of the absorptivity b0 increases with increasing decay constant g: The present local BIEM can be also successfully applied to nonstationary heat conduction analysis in a solid under a microwave heating. To demonstrate this, let us consider a homogeneous square slab considered in example 1. On Table 4 Critical values of b0 for thermal absorptivity in a hollow cylinder Decay constant
Fig. 8. Comparison of the temperature profiles in a homogeneous and a FGM square slab for b0 ¼ 35 and g ¼ 2:
g¼0 g ¼ 0:5 g¼1 g¼2
Gradient parameter
d¼0
d ¼ 21
2.5 4.3 7.7 17.1
0.95 1.62 2.87 6.2
10
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11
at the centre of the slab (i.e. x1 ¼ 0:5; x2 ¼ 0:5;) and for a fixed value b0 ¼ 14: It shows here that the maximum or the stationary temperature decreases with increasing decay constant g as expected.
5. Conclusions
Fig. 10. Time variations of the temperature at the section x2 =a ¼ 0:5 and for b0 ¼ 64 and g ¼ 4:
the whole boundary of the square slab a thermal shock of the form uðx; tÞ ¼ T·HðtÞ is applied, where T ¼ 1 deg is the amplitude of the thermal shock, and HðtÞ is the Heaviside step function. The initial condition u ¼ 1 deg is prescribed at t ¼ 0: A linear temperature dependence of the thermal absorptivity is assumed, and the Stehfest’s numerical inversion method [15] of the Laplace transform is applied. For fixed values b0 ¼ 64; g ¼ 4 and at three different positions x1 =a ¼ 0:25; 0.5, 0.75 at the middle section x2 ¼ 0:5 of the square slab, numerical results for the temperature are presented in Fig. 10 versus time. Fig. 10 shows that the temperature increases sharply with increasing time in the initial stage, and after a short time it tends to its stationary or steady-state value. The effect of the decay constant g on the temporal variation of the temperature is presented in Fig. 11. Here, numerical results are computed
A local BIEM in conjunction with the MLS approximation of physical fields is presented in this paper for steady-state and transient heat conduction problems with a nonlinear heat source in FGMs. A power-law and an exponential decay of the heat source intensity are used for modelling the microwave heating. The method can be applied to a general spatial variation of the material constants for describing FGMs. The initial-boundary value problem for the transient heat analysis is solved in the Laplace transform domain with a subsequent numerical Laplace inversion to obtain time-dependent solutions. The LBIEs are formulated in a boundary-domain integral form with a simple fundamental solution corresponding to the Poisson’s equation. The boundary-domain integral formulation can be easily implemented by choosing a simple circular subdomain to which the LBIE is applied. Contrary to the conventional BIEM or BEM, all integrands in the present LBIEM formulation are regular. Thus, no special integration techniques are required to evaluate the arising integrals. The present LBIEM or LBEM removes the wellknown restriction of the conventional BIEM/BEM to problems with homogeneous material properties. Numerical examples for both steady-state and transient microwave heating problems in FGMs are given to demonstrate the efficiency and the accuracy of the present LBIEM. Numerical results show the effects of the temperature dependence of the thermal absorptivity, the decay constant of the microwave heating source, and the material gradient on the occurrence of hot spots and the temperature distribution in FGMs.
Acknowledgements The authors acknowledge the support by the Slovak Science and Technology Assistance Agency registered under number APVT-51-003702, as well as by the Slovak Grant Agency VEGA-2303823, 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. 11. Time variations of the temperature at the centre of a square slab and for b0 ¼ 14:
[1] Zhu S, Zhang Y, Marchant TR. A DRBEM model for microwave heating problems. Appl Math Modelling 1995;19:287 –97.
J. Sladek et al. / Engineering Analysis with Boundary Elements 28 (2004) 1–11 [2] Miyamoto Y, Kaysser WA, Rabin BH, Kawasaki A, Ford RG. Functionally graded materials: design processing and applications. Dordrecht: Kluwer Academic Publishers; 1999. [3] Suresh S, Mortensen A. Fundamentals of functionally graded materials. London: Institute of Materials; 1998. [4] Butterfield R. An application of the boundary element method to potential flow problems in generally inhomogeneous bodies. In: Brebbia CA, editor. Recent advances in boundary element method. London: Pentech Press; 1978. [5] Sutradhar A, Paulino GH, Gray LJ. Transient heat conduction in homogeneous and non-homogeneous materials by the Laplace transform Galerkin boundary element method. Engng Anal Bound Elem 2002;26:119– 32. [6] Zhu T, Zhang JD, Atluri SN. A local boundary integral equation (LBIE) method in computational mechanics, and a meshless discretization approach. Comput Mech 1998;21:223– 35. [7] De S, Bathe KJ. The method of finite spheres. Comput Mech 2000;25: 329–45. [8] Sladek J, Sladek V, Atluri SN. Local boundary integral equation (LBIE) method for solving problems of elasticity with nonhomogeneous material properties. Comput Mech 2000;24: 456 – 62. [9] Sladek J, Sladek V, Zhang Ch. Transient heat conduction analysis in functionally graded materials by the meshless local
[10]
[11] [12]
[13]
[14]
[15] [16]
[17]
11
boundary integral equation method. Submitted for publication; 2002. Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: an overview and recent developments. Comput Meth Appl Mech Engng 1996;139:3–47. Chen HT, Lin JY. Application of Laplace transform to nonlinear transient problems. Appl Math Modelling 1992;15:144– 51. Ramachandran PA. Application of the boundary element method to nonlinear diffusion with reaction problems. Int J Numer Meth Engng 1990;29:1021 –31. Sladek V, Sladek J, Atluri SN, Van Keer R. Numerical integration of singularities in meshless implementation of local boundary integral equations. Comput Mech 2000;25:394 –403. Chati MK, Paulino GH, Mukherjee S. The meshless standard and hypersingular boundary node methods: applications to error estimation and adaptivity in three-dimensional problems. Int J Numer Meth Engng 2001;50:2233–69. Stehfest H. Algorithm 368: numerical inversion of Laplace transform. Comm Assoc Comput Mach 1970;13:47–9. Zhu S, Satravaha P. An efficient computational method for modelling transient heat conduction with nonlinear source terms. Appl Math Modelling 1996;20:513– 22. Frank-Kamenetskii DA. Diffusion and heat exchange in chemical kinetics. Princeton: Princeton University Press; 1955.