Scaled boundary finite-element method for solving non-homogeneous anisotropic heat conduction problems

Scaled boundary finite-element method for solving non-homogeneous anisotropic heat conduction problems

Accepted Manuscript A Scaled Boundary Finite-Element Solution to Non-homogeneous Anisotropic Heat Conduction Problems Mohammad Hossein Bazyar, Abbas T...

2MB Sizes 0 Downloads 12 Views

Accepted Manuscript A Scaled Boundary Finite-Element Solution to Non-homogeneous Anisotropic Heat Conduction Problems Mohammad Hossein Bazyar, Abbas Talebi PII: DOI: Reference:

S0307-904X(15)00186-9 http://dx.doi.org/10.1016/j.apm.2015.03.024 APM 10499

To appear in:

Appl. Math. Modelling

Received Date: Revised Date: Accepted Date:

3 March 2013 16 March 2015 23 March 2015

Please cite this article as: M.H. Bazyar, A. Talebi, A Scaled Boundary Finite-Element Solution to Non-homogeneous Anisotropic Heat Conduction Problems, Appl. Math. Modelling (2015), doi: http://dx.doi.org/10.1016/j.apm. 2015.03.024

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A Scaled Boundary Finite-Element Solution to Non-homogeneous Anisotropic Heat Conduction Problems

Mohammad Hossein Bazyar* & Abbas Talebi Department of Civil Engineering, Faculty of Engineering, Yasouj University, Yasouj, Iran *

Corresponding author email: [email protected]

Abstract In this paper the scaled boundary finite-element method (SBFEM) is employed to solve the two-dimensional heat conduction problems. In the last two decades SBFEM has been successfully applied to different fields of engineering and science. It utilizes the advantages of both boundary element and finite-element methods. Same as in the boundary element method only the boundary is discretized reducing the special discretization by one. No fundamental solution is needed as in the finite-element method. Making use of a scaling center, the geometry of the problems is transformed into the scaled boundary coordinates including radial and circumferential coordinates. The boundary of the problem represents the whole computational domain. The finite-element approximation on the circumferential coordinates leads to the analytical equation in the radial coordinate. In the steady-state analysis, an eigen-value problem is solved to compute the stiffness matrix and to obtain the temperature field in the domain. In the transient analysis, the stiffness matrix determined from the eigen-value problem along with the mass matrix determined from the low frequency behavior forms a system of first-order ordinary differential equations to be solved for temperature field using time integration schemes. The formulation and solution procedure of the SBFEM for modeling heat conduction problems are outlined. A code in programming environment MATLAB has been developed to obtain the numerical solution. Several numerical examples are addressed to demonstrate the simplicity, versatility and applicability of the SBFEM in simulating heat conduction problems with non-homogeneous and anisotropic media. Complex irregular geometries are modeled. An example with flux singularity is numerically simulated to further demonstrate the advantages of the technique. The validation of method has been achieved by comparing the SBFEM results with those obtained by FEM or those available in the literature. Key words: scaled boundary finite-element method, heat transfer, stiffness matrix, flux

1.

INTRODUCTION

Two-dimensional parabolic equations play a very important role in many engineering applications. Analysis of seepage in porous media and the temperature distribution among solids are the most important applications. Because of its considerable importance in many engineering applications, the fundamentals of heat transfer in solid particles have long attracted the attention of researchers. Many practical heat conduction questions lead to problems not conveniently solvable by classical methods, such as separation of variables techniques or the use of Green's functions. However, problems with only simple geometry and material properties can be solved via analytical techniques. As a result, practice engineers have to use approximate numerical approaches for dealing with complex geometry and material properties. The finite difference method (FDM), is often used as a conventional tool for studying heat conduction problems. FDM is a mesh dependent method as the boundary fitted-mesh is needed. Consequently advanced and sophisticated techniques have to be employed for discretizing complex geometries. The most popular numerical tool in engineering practice is the finite-element method (FEM). FEM has been successfully applied to heat conduction problems. When dealing with infinite domains it needs to truncate the domain resulting in an accuracy reduction. However, infinite elements coupled with finite elements have been widely employed in the literature for capturing unbounded domains. Apart from the mentioned numerical schemes, the boundary element method (BEM) based on integral equations has also been attracted the attention of researchers for solving the partial differential governing equations. Only the boundary is discretized and the complex geometries can be easily discretized. However, the complexity of the fundamental solutions needed in the modeling increases considerably for general anisotropic materials limiting the application of the boundary element method. Many researchers have used the three mentioned popular numerical schemes for solving two-dimensional heat conduction problems. A short literature review is presented here. Wang et al. [1] used the finite-volume method for cylindrical heat conduction problems. A finite element/finite difference scheme is proposed by Wang et al. [2] for the non-classical heat conduction problems. Pulvirenti et al. [3] studied finite-difference solution of hyperbolic heat conduction with temperature-dependent properties. Wang [4] solved transient onedimensional heat conduction problems by finite-element method. Desilva [5] presented a coupled boundary element method and finite difference method for the heat conduction in laser processing. Kutanaei et al. [6] used mesh-free modeling of two-dimensional heat conduction between eccentric circular cylinders. Meshless method based on the local weakforms is applied by Wu and Tao [7] for steady-state heat conduction problems. Gao and Wang [8] used the interface integral BEM for solving multi-medium heat conduction problems. Zhang et al. [9] employed a novel boundary element approach for solving the anisotropic potential problems. Local radial basis function-based differential quadrature (RBF-DQ) method [10] is used for two-dimensional heat conduction problems.

It is the objective of this paper to apply the semi-analytical technique of SBFEM directly to heat conduction problems in either two or three space variables. Only formulation and solution procedures for 2D problems are addressed in this study. The scaled boundary finite-element method was developed in 90s based on similarity concept for dynamic stiffness of unbounded domains. The original formulation and its application in 2D and 3D vector waves and fracture mechanics were presented in the book by Wolf and Song [11]. Later a new derivation consistent with FEM based on Galerkin's weighted residual method was derived [12]. In the work by Deeks and Wolf [13] the principal of virtual work was utilized to re-derive the scaled boundary finite-element equations. Analytical solutions were derived in [14] for problems with certain body loads such as the concentrated loads and loads varying as power function of radial coordinates. Extension of the SBFEM to model time-harmonic and transient response of non-homogeneous elastic unbounded domains was performed in the works by Bazyar and Song [15, 16]. New developments were carried out to enhance the solution procedures for modeling unbounded domains in both frequency and time domains [17, 18]. Steady-state confined and unconfined seepage problems were dealt with in [19, 20]. Transient seepage was addressed in the work by authors [21]. The paper is organized as follows. Governing differential equations of the heat conduction problems are presented in Section 2. A brief summary of the concept of the SBFEM is provided in Section 3. Formulation of the scaled boundary finite-element method for heat conduction problems is detailed in Section 4. Section 5 deals with the solution procedures of the SBFEM for transient heat conduction problems. Five numerical examples are solved in Section 6 to evaluate performance of the proposed scheme. Finally the paper is closed by conclusions in Section 7 followed by references. 2.

GOVERNING EQUATIONS

The governing equations for two-dimensional heat conduction problems are described in this section. In a heat conduction problem, a solid  undergoes thermal loads and is tried to determine the temperature field in this solid, when it is in equilibrium with the external environment. The equilibrium state sought corresponds to transient state. Besides, our study is restricted to the linear case in which the thermal conductivity is independent of temperature. The governing equations for 2D transient heat conduction problems are written as:    +  = 0 (1)

where the material constant  is the product of the mass density and the specific heat  of the material.  is the heat flux vector related to the temperature field defined as

 = −

 with   =  

  as the anisotropic thermal conductivity matrix.  is the gradient operator 

     Eq. (1) is transformed into the frequency domain by performing the Fourier transform  = 

(2)

(3)

  +  ! " = 0

(4)

where symbols  and " are the Fourier transform of  and , respectively. is the frequency. By setting = 0, the steady state of the heat conduction equation is obtained. The thermal conduction equation (Eq. (4)) must be solved for prescribed boundary and initial conditions. Heat conduction boundary conditions take several forms. The frequently encountered conditions are specified surface temperature and specified surface heat flow. On the boundary of the domain, either the value of the temperature (Dirichlet boundary condition) or the heat flux (Neumann boundary condition) must be specified.

3.

METHODOLOGY

In this study the SBFEM is employed to model 2D heat conduction problems. The SBFEM is a semi-analytical computational scheme in which only the boundary of the domain is discretized. In contrast to BEM and similar to FEM, no fundamental solution is required. Its development can be found in the literature [11-12]. 3.1 Scaled Boundary Transformation of the Geometry In this approach which is based on similarity concept, the boundary represents the computational domain. Figure 1 demonstrates how a domain of analysis is modeled by this approach. A scaling center is selected in the domain (Figure 1a) or on the boundary (Figure 1b) from which the entire boundary must be visible. When it is not possible to find a point from which the whole boundary is visible, the total domain can be divided into subdomains in each sub-domain similarity is satisfied (Figure 1b). The local coordinate system of the domain ( ξ , s) is defined by a radial coordinate ξ and a circumferential coordinate s. ξ runs from 0 at the scaling center O to 1 on the boundary. For each element, s is from -1 to +1 (Figure 1c) as in isoparametric elements in the FEM. Figure 1c illustrates a 3-node quadratic element for 2D boundary discretization. As in FEM, higher order elements can also be used to improve accuracy. By scaling the boundary according to the radial coordinate ξ with respect to the scaling center O, the whole domain is covered. The straight surfaces forming the boundary and passing through the scaling center (side faces) and the interfaces between different materials are not discretized (Figure 1b). The scaled boundary coordinate system ( ξ , s) is related to the Cartesian coordinate system by the scaling equation as

x&ξ, s( = xˆ 0 + ξx&s( = ξN&s(x

 y&ξ, s( = y 0 + ξy&s( = ξN&s(y

(5.a)

(5.b)

  where &x, y( is the geometry of a point inside the domain, ( x0 , y0 ) is coordinate of the

scaling center O, &x&s(, y&s(( represents the corresponding point of &x, y( on the boundary, &x,y) represents the nodal coordinates of the element on the discretized boundary, and N&s( is the shape functions formulated in the circumferential coordinate s. The mapping in the SBFEM, is similar to the mapping exists in formulation of the isoparametric elements in the FEM. The shape functions used in the SBFEM (see Eq. (5)) are a function of only circumferential direction s while in the formulation of the isoparametric elements in FEM, the shape functions are a function of directions both ξ and s . Shape functions for a 3-node quadratic element can be expressed as , &-( = ,. ,/ ,0 ; ,. =

1&12.( /

, ,/ = 1 − -/ , ,0 =

1&14.( /

(6)

The Jacobian matrix and its determinant on the boundary ( ξ = 1) are expressed as J =

78

7: ;

7<

7<

79 6 78

79 = 7: ;

=

x&s( x&s(,>

y&s(  , |J| = x&s(y&s(,> − y&s(x&s(,> y&s(,>

(7)

to transform all the spatial derivatives from the Cartesian coordinate system into the scaled boundary coordinate system ( ξ , s) using the following 7

78 @7A 7: ;

=

. y &s(,> C |B| −x &s(,>

7

−y&s( 79 D @. 7 A x&s( 9 7<

(8)

3.2 Temperature Functions

G&ξ, s( H. Along the The SBFEM looks for an approximate solution for the temperature field ET radial lines passing through the scaling center O and a node on the boundary (ξ = 1), the G &ξ(H are introduced. nodal temperature functions ET (9) E "&I, -(H = , &-(E "&I(H

It means the temperatures of any point ( ξ , s) in a domain is interpolated in the circumferential direction through shape function matrix N&s(, which is the same as that G &ξ(H (Figure 2) denotes the temperatures along the used in FEM. Temperature function ET radial lines which are analytical with respect to the radial coordinate ξ . Note that again, the shape functions N&s( are a function of the variable s only, in contrast to G&ξ, s(H = N&ξ, s(ET GH) where ET G H are the the procedure used in finite-element method (ET temperatures in distinct nodes in the domain.

3.3 Temperature Function Equation in SBFEM

Substituting Eqs. (8) and (9) into Eq. (2), the approximate temperature gradient q&ξ, s( can be presented in terms of the scaled boundary coordinates as 1  &I, -( = − KL. &-(E "&I(H + L/ &-(E "&I(HN ,M I

(10)

where new boundary-dependent functions B.&s( and B / &s( describe the temperature gradient-temperature relationship. B. &s( =

B / &s( =

1 J

 y (s ),s    N&s( − x(s ),s 

1 − y (s )   N&s(,> J  x(s) 

(11a)

(11b)

The method of weighted residual is employed in the direction s into the governing equation to perform the finite-element approximations. The residual and the it’s integral are defined as G = 0; X w&x, y(&LR q + iωρcW T G(dv = 0 R&x, y( = LR q + iωρcW T

(12)

Substituting Equations 10, 9 and 8 into Eq. (12) and making use of integration by parts or Green theory leads to the scaled boundary finite-element equations for temperature field as \ ] I/ E " &I(H,MM + &\ ]  − \.  + \. ^(IE "&I(H,M − &\ /  −  _] I/ ( " &I(   = 0

(13)

\ ]  = X1L. &-( L.&-(|a| b-

(14a)

Eq. (13) is a homogeneous second-order differential equation system with respect to I. The coefficient matrices available in Eq. (13) being functions of only circumferential coordinate s are introduced as \  = X1L &-( L &-(|a| b\ /  = X1L/ &-( L/&-(|a| b.

/



.

_]  = X1,&-(  ,&-(|a|b-

(14b) (14c) (14d)

These coefficient matrices are evaluated element by element and assembled over the discretized boundary to determine the global coefficient matrices. E . and M ]  are positive definite and E /  is symmetric.

3.4 Dynamic-Stiffness Matrix Equation in SBFEM

To present the SBFEM solution of transient heat conduction problems, the dynamic-stiffness matrix K&ξ, ω( is derived in this section. It relates the Fourier transform of the external nodal fluxes q&ξ( to the Fourier transform of the temperature function G&ξ(H (frequency parameter ω is deleted in nodal flux and temperature Eq&ξ( = K&ξ, ω(T function for the sake of simplicity). Using Eqs. (10), (13) and the definition of dynamic-stiffness matrix, the scaled boundary finite-element equation in dynamic-stiffness matrix K&ω( on boundary (ξ = 1) is derived as &K& ( − \. (\ ]2. &f& ( − \.  ( − \ /  + 2 f& (,h −  _]  = 0

(15)

In SBFEM it is generally appropriate to carry out a time-domain analysis of finite domains

with the steady-state stiffness ( K st ) and mass matrices M . It describes the low-frequency behavior of dynamic-stiffness matrix K&ω(. In the following, first the steady-state stiffness matrix and then the mass matrix and transient solution are outlined.

3.5 Steady-State Stiffness Matrix and Mass Matrix

Setting ω = 0 in Eq. (13), leads to the scaled boundary finite-element equations for temperature field in steady state heat conduction as \ ] I/ E " &I(H,MM + &\ ]  − \.  + \. ^(IE "&I(H,M − \ /  "&I(   = 0

(16)

from its analytical solution, the steady-state stiffness matrix to be determined. Note that in comparison to the FEM in which after approximation, a nodal temperature-nodal flux   equation ( q  = [ K st ]  T ) can be obtained, in the SBFEM in Eq. (9), the temperature field is a function of ξ and after finite-element approximation in circumferential direction s, Eq. (16) in terms of radial direction ξ is derived to compute the stiffness matrix and the temperature field at any ξ. Solution of Eq. (16) may be found through a modal combination procedure in the form E "&I(H = .I 2ij k. + / I2il k/  +…

(17)

where the integration constants ci depend on the boundary conditions and are the contribution of each mode to the solution. The temperature of each mode takes the form

G&ξ(H = ξ2 λ ϕ ET

(18)

where the vector k can be interpreted as the modal temperatures at the boundary nodes, while n can be identified as a modal scaling factor for the radial direction [13]. By an inspection and from the derivations from the weighted residual technique for finite element approximation, the heat flux mode is related to the temperature mode as following ;&I( = \ ] IE "&I(H +\.  "&I(  ,M

(19)

Substituting Eq. (18) into Eq. (19) evaluated at the discretized boundary yields G q = &E. R − λ E ](T

(20)

Substituting Eq. (18) into Eq. (16) and making use of Eq. (20), the quadratic equation of Eq. (16) becomes a standard linear eigenvalue problem 

E ] 2. E .R −E / + E. E ]2.E .R

{φ} {φ} −E ] 2. D= λC D . ] 2.  C −E E {qˆ} {qˆ}

(21)

For N degrees of freedom on the discretized boundary, Eq. (21) results in 2N modes. The eigenvectors contain the modal temperatures and the equivalent modal node heat fluxes. For bounded domain, only N modes with non-positive real eigenvalues &λ = diag&λ1 , λ2 , … , λN (( lead to finite temperature at the scaling center. This subset of N temperature modes is denoted by φ = t {φ1} , {φ 2 } , … . , {φ N } v and the corresponding equivalent modal heat fluxes by w. The solution of Eq. (16) becomes |

E "&I(H = x y I 2iz ∅y  = ~I2i €

(22)

K >  = Qφ2.

(23)

y}.

G&ξ( is analytical in the radial direction ξ. On As can be observed, the temperature function T G H = φc from which the integration constants the boundary (ξ = 1), Eq. (22) becomes ET 2. G H. Consequently the steady-state stiffness matrix  K st  c are calculated as c = φ ET with respect to the discretized boundary degrees of freedom is simply obtained as

As described previously, to perform a time-domain analysis of heat conduction in bounded domains, a low-frequency behavior of the SBFE equation including steady-state stiffness matrix and a mass matrix is utilized. The steady-state stiffness matrix was detailed in the previous sub-section. To determine the mass matrix M of a bounded domain with degrees of freedom on the boundary only, two procedures can be followed. Firstly, the standard procedure to determine mass matrix in the finite-element method can be used [22]. The temperature function in Eq. (22) and integration constants c are combined leading to E "&I(H = ~I 2i €~2. E "H

(24)

GH represents the shape function in the ξ direction as The coefficient of matrix ET

N&ξ( = φξ2ƒ €φ2.

(25)

The shape function in the s direction has been presented in , &-( = ,. &-( ,/ &-( ,0 &-(. Substituting Eq. (24) into Eq. (9) leads to

Eq. (6)

as

G&ξ, s(H = N&s(ET G&ξ(H = N&s(~I2i €~2.E "H ET

(26)

N…  = N&s(φ ξ2ƒ €φ2.

(27)

The shape function ,„  in the domain for heat conductions in two dimensions is introduced as E "&I, -(H = ,„ E "H . Comparing with Eq. (26), the shape function N…  is defined as

Mass matrix is presented as M = †N… R ρcW N… dv

(28)

Substituting Eq. (27) into Eq. (28) results in M = &φ2. (R †ξ2ƒ€φ R N&s(R ρcW N&s(φ ξ2ƒ €φ2. dv

(29)

M = &φ2. (R&Xξ2ƒ€φR M ] φξ2ƒ €ξdξ(φ2.

(30)

m]  = φ R M ] φ ; M = φ2. &Xξ2ƒ€m] ξ2ƒ €ξdξ(φ2.

(31)

mŠ € = X] ξ2ƒ €m] ξ2ƒ € ξdξ; M = φ 2. mŠ €φ 2.

(32)

Making use of definition M ]  in Eq. (14d), mass matrix is reformulated as To be able to carry out the integration part in the bracket in Eq. (30) in the ξ-direction analytically, the following abbreviation is introduced R

Each element of the matrix ˆ‰  introduced as the following .

is assessed analytically from .

mŠ‹Œ = † ξ2ƒ €m]‹Œ ξ2ƒŽ € ξdξ = ]

R

m]‹Œ

−λ‹ − λŒ + 2

from which the mass matrix is determined by using Eq. (32) as

(33)

M = &φ2. (RmŠ φ 2.

(34)

Alternatively, the mass matrix can be obtained from the scaled boundary finite-element equation in dynamic stiffness (see [21] for more detailed formulation). To determine the mass matrix M, the dynamic-stiffness matrix K&ω( for low-frequency behavior is expressed as f& ( =  K st  + iω _

(35)

Substituting Eq. (35) into Eq. (15) leads to a constant term independent of iω , a term in iω and higher-order terms in iω which are neglected. The constant term, which is equal to Eq. (15) formulated for ω = 0, vanishes. The coefficient matrix of the remaining term in iω is written as ((-[Kst ] + [E1 ])[E0 ] -1 - [I])[M] + [M]( [E 0 ] -1 (-[Kst ] + [E1 ]T ) - [I]) + [M 0 ] = 0

(36)

This is a linear equation to calculate the mass matrix M. To do so, eigen-vectors obtained from the eigen-value problem used in Eq. (21) is utilized. Some mathematical manipulations (See Ref. [21] for more detailed derivations) using temperature modes of Eq. 21, steady-state stiffness matrix in Eq. (23) and also Eq. (36), lead to an equation for m as & − ny (ˆ + ˆ& − ny ( = ~ _ ] φ

(37)

_ = &~2. ( ˆ~2.

(38)

to be solved by simple back substitution. Once m is calculated from Eq.37, mass matrix M is determined form the following transformation as As M ] is positive definite, M will possess the same property. In this paper, the mass matrix is obtained from the second approach as detailed from Eq. (36) to Eq. (38). 3.6 Transient Analysis using Crank-Nicolson Time Stepping Method

To perform a time-domain analysis of heat conduction, the steady-state stiffness matrix and mass matrix of all sub-domains are assembled (in case of sub-structuring). Substituting Eq. (35) into the definition of dynamic-stiffness matrix followed by inverse Fourier transformation, the nodal heat flux–temperature relationship is determined as a standard time-domain equation as _E H + f 1 T =  

(39)

G€ “ ’  ’ +  ’,’4.  ’4. = f f ’4.

(40)

Standard implicit and explicit numerical time integration algorithms can be employed to solve Eq. (39). Crank-Nicolson method is used in this paper. At time step m + 1, the nodal temperature vector T‘4. is obtained from solving Eq. (40) given the nodal temperatures at previous time step m

with Δt as the time interval and variables G€ f = _ + 0.5 —˜ [ K st ] m+1 ’4.

“’ = _ − 0.5 —˜ [ K ]m f ^ ’,’4. = 0.5—˜’4. + ’ € st

(41.a) (41.b) (41.c)

The initial nodal temperatures are obtained using the steady-state analysis. After calculating T‘4. from Eq. (40), the nodal temperatures on the discretized boundary can be used to calculate the integration constants c from Eq. (22). The interior temperature field is then computed by substituting Eq. (22) into Eq. (9) as G&ξ, s(H = N&s(ET G&ξ(H = N&s(φξ2ƒ €c ET

(42)

3.7 An Algorithm to Implement the Solution Procedures

A step by step algorithm is presented as follow to implement the formulations and the solution procedures adopted for the analysis of heat conduction problems using the SBFEM: Step 1: After choosing suitable scaling center, the boundary is discritized using finiteelements. Having calculated coefficient matrices \ ] ,^ \ .  ,^ \ /  and _]  for each element from Eqs. (14a), (14b), (14c) and (14d), the global coefficient matrices \ ] ,^ \.  ,^ \ /  and _]  for the whole boundary are computed. Step 2: The eigen-value problem of Eq. (21) is constructed. Step 3: The eigen-vectors determined from Step 2 are employed to determine the steady-state heat conduction stiffness matrix f 1  on the boundary using Eq. 23. Step 4: The nodal values of the temperature function and flux on the boundary are obtained from the definition of the stiffness matrix. Step 5: Boundary nodal values determined from Step 4 and the eigen vectors obtained from Step 2 are used to compute the integration constants  from Eq. 22 (note that on the boundary I = 1 ) Step 6: Having computed constant coefficients  , Eq. 22 is again used to determine values of interior temperatures for different values of I starting from 1 on the bounary and ending to 0 in the scaling center. After performing a steady state calculations till Step 6, the following steps are used to carry out a transient analysis: Step 7: Using eigen-values and eigen-vectors determined from Step 2 and coefficient matrix _]  calculated in Step 1, Eq. 37 is employed to find the matrix ˆ  followed by using Eq. 38 to compute mass matrix _ . Alternatively, Eqs. (31) to (34) can be used to obtain mass matrix _ . Step 8: Eq. 39 is constructed and solved for temperatures on the boundary using Eq. 40. Eq. 42 is used to capture the temperatures in the interior domain for different values of I starting from 1 on the bounary and ending to 0 in the scaling center.

4.

NUMERICAL EXAMPLES

In this section five numerical examples are presented to validate the accuracy and performance of the SBFEM in modeling heat conduction problems. Four examples have been selected from the literature and one example is verified using the FEM. The first example models the steady-state heat transfer in a triangular fin attached to a pipe. In the second example, the steady-state heat conduction in an anisotropic hollow eccentric cylinder is addressed. The third example is concerned with the steady-state heat conduction in a thickwall cylinder consisting of two materials. The fourth example deals with thermal analysis of a thick-walled hollow cylinder subjected to transient thermal loading. The last example considers a cracked plate subjected to steady-state heat conduction. A computer program in programming environment MATLAB has been written to implement the SBFE formulations and solution procedures. The results available in the literature and the results of a FE analysis are employed to verify the accuracy and versatility of the SBFEM results. 4.1. Steady state heat conduction in a triangular fin attached to a pipe

For the first example and to investigate the application of the SBFEM in modeling problems with irregular geometries applicable in engineering, a triangular fin attached to a pipe as shown in Fig. 3 is addressed. This example has been chosen from the work by Soleimani et al. [10] analyzed by RBF-DQ and finite-element method. The value of thermal diffusivity ™=

š ›

equals to 8.4 × 10 −5 &ˆ/ ⁄-(. The SBFE mesh discretizing the boundary with its scaling

center is shown in Fig. 4. As shown in this figure, the SBFE mesh consists of 26 three-node line elements and 52 degrees of freedom. The temperature distribution inside the domain obtained from SBFEM, RBF-DQ and FEM (see [10]) are illustrated in Fig. 5. Fig. 6 depicts the temperature variation along line  = 0.5 in Fig. 3. As it can be seen from Figures 5 and 6 there is an excellent agreement between the results obtained from SBFEM, RBF-DQ and FEM. 4.2. Steady state heat conduction in an anisotropic hollow eccentric cylinder

To demonstrate the versatility and advantages of the SBFEM in modeling anisotropic materials, as the second example an anisotropic hollow eccentric cylinder is modeled. The geometry of problem is shown in Figure 7a. This problem has been solved by Zhang et al. [9] using a boundary element method. For the hollow cylinder, uniform and constant temperatures are assumed at the inner and outer surfaces. As it is not possible to locate a scaling center from which whole boundary is visible, the whole computational domain is divided into four sub-domains for each one scaling center is chosen. The SBFE sub-domains with their scaling centers are shown in Figure 7b. The boundary of four sub-domains is discretized using linear finite elements. The SBFE mesh consists of 89 two-node line elements and 85 degrees of freedom. Figs. 8 and 9 show the isotherm distribution of the eccentric cylinder with two different conductivity coefficients by BEM and SBFEM,

respectively. The influence of thermal conductivity tensors yŸ on the temperature distribution inside the domain is clearly demonstrated by comparing Figs. 8 and 9. Excellent agreement between the results obtained from the SBFE and FEM solutions can be observed. 4.3. Steady state heat conduction in a multi-medium thick-wall cylinder

For the third example, SBFEM is applied to solve a multi-medium heat conduction problem. This example is concerned with a thick-wall cylinder consisting of two media. The geometry and boundary conditions of the problem are shown in Fig. 10a. Two materials have the same thickness along the radial direction and the conductivities are 5 and 1 for the inner and outer media, respectively. Due to symmetry only a quarter of the cylinder is analyzed. This problem has been addressed by Gao and Wang [8] using the interface integral boundary element method (IIBEM). The SBFE mesh using 4 sub-domains with their scaling centers are shown in Figure 10b. The SBFE mesh consists of 108 two-node line elements and 105 degrees of freedom. Fig. 11 portrays the contour plot of the temperature inside the domain computed by IIBEM and SBFEM. Excellent agreement is shown. Fig. 12 shows the variation of temperature along the radial direction where d is the radial distance to the inner boundary. For comparison, the FEM results obtained using ANSYS are also provided. From Fig. 12, we can see that the SBFEM results are in a very good agreement with IIBEM and FEM results. 4.4. Isotropic thermal analysis of a thick-walled hollow cylinder subjected to transient thermal loading

To validate the accuracy of the proposed method for transient heat problems, a benchmark problem presented in [23] was carried out. The geometry of the problem is shown in Fig. 13. The values of μ and k are equal to 50.6 Ws⁄cm0 ⁄°C and 42 W⁄cm⁄°C , respectively. The internal transient thermal load is given by H&t( = &1 − e2].¨ ( with convection in the outer surface. The temperature of the surrounding fluid medium is taken to be 0 ℃. The scaled boundary finite-element mesh is displayed in Fig. 14. The SBFE mesh consists of 24 threenode line elements leading to 48 degrees of freedom. Due to symmetry only a quarter of the cylinder is analyzed. This problem has also been solved by Santos et al. [24] using FEM. The results along the radial direction compared with the literature are presented in Fig. 15. Table 1 shows the temperature at r = 0.5 at different time steps found from the SBFEM and 1D and 3D finite-element analyses. Excellent agreement can be observed. 4.5. Isotropic thermal analysis of a cracked plate subjected to steady state thermal loading

To further demonstrate the advantages of the SBFEM, a plate with crack subjected to steadystate heat conduction is modeled. The geometry of the problem is shown in Fig. 16. The thermal conductivity k equals to 40 J/s⁄m⁄°C. The temperature along lines IH and ED are taken to be 200℃ and 20℃, respectively. Lines IA, AC, CD, HG, GO, OF and FE are assumed to have zero normal heat flux. The scaled boundary finite-element mesh is displayed in Fig. 17a. Lines GO and OF are side faces and need not to be discretized. The SBFE mesh consists of 20 three-node line elements leading to 41 degrees of freedom. To have a reference solution a finite element analysis is performed. Three finite element meshes are used. The

finest mesh including 5392 degrees of freedom is shown in Fig. 17b. Fig. 18a and Fig. 18b portray the contour plot of the temperature inside the plate computed by SBFEM and FEM. Excellent agreement is achieved. To highlight the advantages of SBFEM compared with FEM, the unit heat flux along line BOG is computed. Fig. 19 shows that the SBFEM accurately models heat flux singularity. As shown in this figure, despite refining the mesh in FEM, solutions lack convergence. 5.

CONCLUSIONS

The problem of heat conduction for anisotropic and non-homogeneous media is numerically solved using the scaled boundary finite-element method. This technique combines the advantages of both finite-element and boundary element methods and no fundamental solution is required. As only the boundary is discretized fewer elements are needed in the SBFEM than the standard FEM resulting in lower computational costs. Results of this study demonstrate that the SBFEM is an attractive scheme in terms of accuracy and applicability for modeling heat transfer problems with non-homogeneous and anisotropic media. Complex irregular geometries and problems with flux singularity can be easily and accurately analyzed.

References [1] W. Li, B. Yu, X. Wang, P. Wang, S. Sun, A finite volume method for cylindrical heat conduction problems based on local analytical solution, International Journal of Heat and Mass Transfer 55 (2012) 5570–5582. [2] B.L. Wang, J.C. Han, Y.G. Sun, A finite element/finite difference scheme for the nonclassical heat conduction and associated thermal stresses, Finite Elements in Analysis and Design 50 (2012) 201–206. [3] B. Pulvirenti, A. Barletta, E. Zanchini, Finite-difference solution of hyperbolic heat conduction with temperature-dependent properties, Numer. Heat Transfer A 34 (1998) 169– 183. [4] B.L. Wang, Transient one-dimensional heat conduction problems solved by finite element, Int. J. Mech. Sci. 47 (2005) 303–317. [5] S.J. DeSilva, C.L. Chan, Coupled boundary element method and finite difference method for laser drilling, Applied Mathematical Modelling 32 (2008) 2429–2458. [6] S.S. Kutanaei, E. Ghasemi, M. Bayat, Mesh-free modeling of two-dimensional heat conduction between eccentric circular cylinders, International Journal of the Physical Sciences 6(16) (2011) 4044–4052. [7] X.H. Wu, W.Q. Tao, Meshless method based on the local weak-forms for steady-state heat conduction problems, International Journal of Heat and Mass Transfer 51 (2008) 3103– 3112.

[8] X.W. Gao, J. Wang, Interface integral BEM for solving multi-medium heat conduction problems, Engineering Analysis with Boundary Elements 33 (2009) 539–546. [9] Y.M. Zhang, Z.Y. Liu, J.T. Chen, Y. Gu, A novel boundary element approach for solving the anisotropic potential problems, Engineering Analysis with Boundary Elements 35 (2011) 1181–1189. [10] S. Soleimani, M. Jalaal, H. Bararnia, E. Ghasemi, D.D. Ganji, F. Mohammadi, Local RBF-DQ method for two-dimensional transient heat conduction problem, International Communications in Heat and Mass Transfer 37 (2010) 1411–1418. [11] J.P. Wolf, C. Song, Finite-Element Modelling of Unbounded Media, Wiley, Chichester, 1996. [12] C. Song, J.P. Wolf, The scaled boundary finite-element method- alias consistent infinitesimal finite-element cell method-for elastodynamics, Computer Methods in Applied Mechanics and Engineering 147 (1997) 329-355. [13] A.J. Deeks, J.P. Wolf, A virtual work derivation of the scaled boundary finite-element method for elastostatics, Computation Mechanics 28 (2002) 489-504. [14] C. Song, J.P. Wolf, Body loads in the scaled boundary finite-element method, Computational Methods in Applied Mechanics and Engineering 180 (1999) 117-135. [15] M.H. Bazyar, C. Song, Time-harmonic response of non-homogeneous elastic unbounded domains using the scaled boundary finite-element method, Earthquake Engineering and Structural Dynamics 35 (2006) 357-383. [16] M.H. Bazyar, C. Song, Transient analysis of wave propagation in nonhomogeneous elastic unbounded domains by using the scaled boundary finite-element method, Earthquake Engineering and Structural Dynamics 35 (2006) 1787-1806. [17] C. Song, M.H. Bazyar, A boundary condition in Padé series for frequency-domain solution of wave propagation in unbounded domains, International Journal for Numerical Methods in Engineering 69 (2007) 2330-2358. [18] M.H. Bazyar, C. Song, A continued-fraction based high-order transmitting boundary for wave propagation in unbounded domains of arbitrary geometry, International Journal for Numerical Methods in Engineering 74 (2008) 209-237. [19] M.H. Bazyar, A. Graili, Scaled boundary finite-element solution to confined seepage problem, Journal of Computational Methods in Engineering 30(1) (2011) 61-78. [20] M.H. Bazyar, A. Graili, A Practical and Efficient Numerical Scheme for the Analysis of Steady State Unconfined Seepage Flow, International Journal for Numerical and Analytical Methods in Geomechanics, 36(16) (2012) 1793-1812. [21] M.H. Bazyar, A. Talebi, Transient Seepage Analysis in Zoned Anisotropic Soils Based on the Scaled Boundary Finite-Element Method, International Journal for Numerical and Analytical Methods in Geomechanics, 39(1)(2015), 1-22.

[22] J.P. Wolf, The Scaled Boundary Finite-Element Method, Wiley, Chichester, 2003. [23] A.E Segal. Thermoelastic analysis of thick-walled vessels subjected to transient thermal loading. J Pressure Technol, 9 ( 2001) 123–146. [24] H. Santos, Cristovao M. Mota Soares, Carlos A. Mota Soares, J.N. Reddy. A semianalytical finite element model for the analysis of cylindrical shells made of functionally graded materials under thermal shock, Composite Structures 86 (2008) 10–21.

Table 1. Temperature at r = 0.5 for different time steps. Times (s) 1

2

3

4

Theory 3D FEM 1D FEM SBFEM 3D FEM 1D FEM SBFEM 3D FEM 1D FEM SBFEM 3D FEM 1D FEM SBFEM

Temperature (℃)

r = 0.5 0.0899 0.0910 0.0853 0.2120 0.2210 0.2270 0.2928 0.3099 0.3146 0.3426 0.3651 0.3677