Increment-Dimensional Scaled Boundary Finite Element Method for Solving Transient Heat Conduction Problem

, , ,

1. College of Aerospace Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016,P. R. China; 2. College of Aeronautic Science and Engineering, Beihang University, Beijing 100191, P. R.China; 3. The 610 Research Institute of Chinese Aeronautical Establishment, Xiangyang 441000, P. R. China

Abstract: An increment-dimensional scaled boundary finite element method (ID-SBFEM) is developed to solve the transient temperature field. To improve the accuracy of SBFEM, the effect of high frequency factor on dynamic stiffness is considered, and the first-order continued fraction technique is used. After the derivation, the SBFE equations are obtained, and the dimensions of thermal conduction, the thermal capacity matrix and the vector of the right side term in the equations are doubled. An example is presented to illustrate the feasibility and good accuracy of the proposed method.

Key words: heat conduction; scaled boundary finite element method (SBFEM); temperature field; accuracy

0 Introduction

The heat conduction problem is frequently encountered in the field of aerospace, that can be investigated by analytical or numerical methods. The analytical methods are limited to the problems of simple geometry and material properties. In many cases of engineering practice, heat conduction problems have to be simulated by numerical methods, such as the finite element method (FEM)[1,2], the finite difference method (FDM)[3,4], the finite volume method (FVM)[5,6]and the boundary element method (BEM)[7,8]et al. However, FEM, FDM and FVM need to discretize the whole computational domain, which will lead to more discretization workload. In contrast to FEM and FDM, BEM only requires a boundary mesh and significantly reduces the computing times. Nevertheless, the traditional BEM needs to find fundamental solutions. The scaled boundary finite element method (SBFEM)[9-13]has both the advantages of FEM and BEM. Several studies on heat conduction (or diffusion) problems have been reported. Refs. [9-11] presented SBFEMs for steady-state heat conduction analysis. Birk and Song[12]proposed a temporally local method for the numerical solution of transient diffusion problems in unbounded domains. In this method, SBFEM and a novel solution procedure for fractional differential equations were combined, but only the solution at the boundary of SBFE can be obtained by this method, and the inner temperature filed of SBFE was not discussed. Bazyar and Talebi[13]used SBFEM to solve non-homogeneous anisotropic heat conduction problems. In their studies, the heat conduction matrix determined from the eigenvalue problem and the mass capacity matrix determined from the low frequency behavior formed a system of the first-order ordinary differential equations, from which the temperature at SBFE boundary were solved by using the time finite difference schemes. Meanwhile, the inner temperature in SBFE was obtained by using the steady state formula based on the boundary temperatures, so the effect of the thermal capacity on the inner temperature was not taken into account. These disadvantages undoubtedly lead to some calculation errors and even wrong results, especially when the large size of SBFE was divided. To overcome the disadvantages of previous work, Li and Ren[14]presented an extending SBFEM, in which the effect of heat capacity on the inner temperature field of SBFE was taken into account, and meanwhile, the algorithms were developed by combining SBFEM and precise integration method (PIM) or finite difference scheme. The sub-domain method was used to improve the computational accuracy. However, the division of sub-domain also increased the discretization workload. Summarizing the advantages of SBFEM in heat conduction problem, we can found that: Firstly, like BEM, it only discretizes the boundaries and has no fundamental solutions, and the solution is analytical in radial direction for steady state problems; secondly, for the transient state problem, if we are only interested in the boundary temperature values, we can only solve it by SBFEM. Whereas, if we are also interested in the inner temperature filed of the problem domain, we can transform the equation for describing SBFE interior temperature field to an initial value problem. This problem can be solved by the finite difference scheme, in which the mesh can be generated automatically, without discretization workload. Therefore, SBFEM is a very attractive algorithm for the engineering and technical personnel.

The purpose of this paper is to seek a method for solving transient heat conduction problem by using SBFEM, which not only does not increase the discretization workload, but also can improve the calculation accuracy. An increment-dimensional scaled boundary finite element method (ID-SBFEM) is proposed to solve the transient heat conduction problem.

1 SBFEM Equations of Transient Heat Conduction Problem

In Ref. [14], the SBFEM equations have been derived by transforming the governing differential equations into the scaled boundary coordinates, in which a weighted residual method and Green theorem are applied. Full details of the formulations can be found in Ref. [14]. Here some key equations of SBFE for the analysis of transient heat conduction problem without sources in frequency domain are given directly[14].

E0ξ2a(ξ),ξξ+(E0+E1T-

E1)ξa(ξ),ξ-(E2+iωM0ξ2)a(ξ)=0

(1)

whereE0,E1,E2, andM0are coefficient matrices dependent on the geometry and materials of the element,but they are independent on the normalized radial coordinateξ.a(ξ) is nodal temperature vector, andωis the frequency.

The heat flow rate on the boundary (ξ=ξ1=1) can be expressed as[14]

E0ξ1a(ξ1),ξ+E1Ta(ξ1)

(2)

In SBFEM, the internal normal flow rateQ(ξ) through curved surfaceS(ξ)can be expressed as[14]

Q(ξ)=ξE0a(ξ),ξ+E1Ta(ξ)

(3)

The definition of the dynamic-stiffness matrixK(ξ,ω) atξcan be introduced as

Q(ξ)=K(ξ,ω)a(ξ)

(4)

Substituting Eq.(4) into Eq. (3) gives

K(ξ,ω)a(ξ)=ξE0a(ξ),ξ+

E1Ta(ξ)

(5)

From Eqs. (1) and (5), the equation forK(ξ,ω) can be derived. After using the dimensionless analysis method,K(ξ,ω) can be obtained as[14]

(K(ξ,ω)-E1)E0-1(K(ξ,ω)-

E1T)-E2+2(iωξ2)K(ξ,ω),(iωξ2)-

iωξ2M0=0

(6)

whereK(ξ,ω) must be a function of iωξ2.By changing the independent variable tox=iωξ2, and dynamic stiffnessS(x)=K(ξ,ω)，the scaled boundary finite-element equation in dynamic stiffness can be rewritten as

(S(x)-E1)(E0)-1(S(x)-

E1T)-E2+2x[S(x)],x-xM0=0

(7)

Eq. (7) represents a SBFE equation in dynamic stiffness formulated in the frequency domain for heat conduction. It is a nonlinear first-order ordinary differential equation withx=iωξ2as independent variable. The dynamic stiffnessS(x) is usually not calculated. The static-heat conduction matrixKand mass thermal capacity matrixMare used, which follow from the low-frequency expansion ofS(x)[14]

S(x)=K+xM

(8)

Substituting Eq. (8) into Eq. (7) will obtain the expression ofKandM, and thenS(x) can be calculated. Finally, substituting Eq. (8) into Eq.(4) and utilizing the inverse Fourier’s transform, we can obtain the SBFE equation in the time domain, which can be solved. Here, we call this method as conventional method, which is developed in our previous work[14].

To improve the computational accuracy, the low-frequency expansion is not enough, and the high-frequency expansion needs to be considered. The continued-fraction approach is an effective method for the high-frequency expansion of the dynamic stiffness. For example, Ch Song et al[12,14]used a continued-fraction approach to solve the diffusion problem in an unbounded domain and the problem of structural dynamics. Here, the same idea is used to solve the transient heat conduction problem in the bounded domain. However, the high-order continued fraction needs to increase computational costs, and the ill-conditioned matrix may occur. Considering both the computational cost and accuracy factors, in this paper, a first-order continued fraction solution is used for the matrix equation as

(9)

S(1)(x)=S0(1)+xS1(1)

(10)

(11)

where

a(1)=E0-1

c(1)=ME0-1M

(12)

(13)

Equating terms corresponding tox1as zero yields

b1(1)S0(1)+S0(1)[b1(1)]T

(14)

Q(ξ)=Ka(ξ)+xMa(ξ)-

xa(1)(ξ)

(15)

Here, the auxiliary variablea(1)(ξ) is defined in

xa(ξ)=S(1)(x)a(1)(ξ)

(16)

After utilizing the inverse Fourier transform, Eqs. (15) and (16) are combined into a matrix expression for the first order continued fraction as

(17)

where

Kh=diag(K,S0(1))

(18a)

(18b)

(18c)

F(ξ)=[Q(ξ)；0]T

(18d)

From Eq.(17), we can see that the first-order continued fraction method presented in this paper is applied not only to the boundary, but also to the inner domain. There is the obvious difference between previous studies[12,15]and this work. From Eqs.(18a)—(18d), we also can see that the dimensions of matrix and vector are doubled, comparing to conventional low frequency expansion of the dynamic stiffness in Ref.[14]. Therefore, the method presented in this paper is called increment-dimensional SBFEM (ID-SBFEM). To solve the first-order SBFE Eq.(17), a two-step method was developed in our previous work[14]. The first step is to solve the time domain equation of heat conduction problem at SBFE boundaries (i.e.,ξ=1) by PIM, and the second step is to transform the equation for describing SBFE interior temperature field to an initial value problem, which is solved by a finite difference scheme. Herein, the same numerical scheme is used to obtain the solution of Eq.(17).

2 Numerical Examples

Fig.1 2-D transient heat conduction problem and boundary conditions

The SBFE divisions are shown in Fig.2. As seen from Fig.2, there are two SBFE divisions selected to calculate. The dividing whole domain method is called DWD, as shown in Fig.2(a), and the dividing sub-domain method is called DSD, as shown in Fig.2(b). During calculating the interior temperature fields, 8 equal parts along the directionξfor both the DWD and DSD methods are divided. In order to easily compare with the computational results, 10 representative points are selected as marked in Figs.2(a), (b). In computation, the time step size is chosen as 0.01 s.

Fig.2 Geometry mesh model of 2-D transient heat conduction problem

To compare the results of conventional low frequency expansion method with those of the dynamic stiffness and increment-dimensional method conveniently, the conventional method (CM) with DWD or DSD is denoted as CM-DWD or CM-DSD, whereas the increment-dimensional method with DWD is denoted as IDM-DWD.

Fig.3 shows the temperature distribution at different time obtained by using IDM-DWD. It is shown that the isothermal lines at different time are perpendicular to the left and bottom boundaries, which means that no heat flow crosses the left and bottom boundaries. While the temperature gradient exists in normal direction at the top and right boundary, meaning that the heat flow crosses the top and right boundaries. The temperature distribution results simulated by the present IDM-DWD method are reasonable according to the heat transfer theory.

Fig.3 Temperature distribution obtained by IDM-DWD at different time

Fig.4 shows a comparison of theoretical solutions and results from the present IDM-DWD. From Fig.4, we can see that the temperature decreases with time variation due to the convection heat exchange at the top and right boundaries. And the proposed IDM-DWD has good consistency with the theoretical solution, which shows that the method has high precision at both interior and boundary points.

Fig.4 Comparison of theoretical solutions and results from IDM-DWD

In order to compare the accuracy of the proposed IDM-DWD and other numerical methods, we calculated the theoretical solution (TS) of the problem and numerical solution from CM-DWD, CM-DSD and IDM-DWD for different points at various time, which are listed in Table 1. And the computational error is shown in Table 2.It is shown that the results by the proposed IDM-DWD method are more close to those by the theoretical solution. The maximum relative error of the IDM-DWD method is 0.282%, while the larger error occurs by other methods, especially the maximum relative error by CM-DWD method is 2.208%, and it is about 0.376% by CM-DSD method. It is shown that the effect by dividing sub-domain method and increment-dimensional method to improve the SBFEM convergence is obvious. The reason lies in that the dynamic stiffness matrix in conventional method is dealt with only a low-frequency approximation of dynamic property, which causes the effect of high-frequency not to be simulated in the case of large domain size. When the whole domain is subdivided into several sub-domains, the size of the element decreases, then the high-frequency responses can be modeled. Whereas, IDM-DWD can improve the accuracy by introducing high-frequency terms into dynamic stiffness directly.

To investigate the effect of time step size on the calculation results and its errors, we use the IDM-DWD and select time step sizes τ=0.1 s, and τ=0.01 s to simulate the temperatures, respectively. The simulated temperatures at interior points and errors att=0.5 s are listed in Table 3.

From Table 3, it can be seen that decreasing the time step size can decrease obviously the calculating errors. The maximum relative error with τ=0.1 s is 0.259%, while that with τ=0.01 s is 0.039%. The reason is that the influence of high-frequency in small time step is more obvious.

Table 1 Result comparison of theoretical and numerical methods ℃

Table 2 Relative error comparison of the results obtained by numerical methods %

Table3Effectoftimestepsizeoncalculationresultsantitserrors

Methods anderrorTime stepsizeP(x,y)1 (1/4,1/4)2(1/4, 1/2)3(1/4, 3/4)4(1/2, 3/4)5 (3/4,3/4)TS /℃72.71469.82165.07062.48258.229ID-DWD /℃(relative error) /%τ=0.1 s72.812(0.135)69.975(0.221)65.177(0.164)62.644(0.259)58.344(0.197)τ=0.01 s72.718(0.006)69.848(0.039)65.073(0.005)62.505(0.037)58.233(0.007)

3 Conclusions

An increment-dimensional SBFEM (ID-SBFEM) is developed to solve the transient heat conduction problems. The first-order continued fraction technique is used to consider the effect of high-frequency, and the dimension of the SBFE equation is doubled. From the analysis of an example, we can conclude that in contrast to CM-DSD, IDM-DWD not only does not require dividing the sub-domain mesh, but also has a higher accuracy. When the smaller computational time step is taken, the accuracy of the solution will be higher.

Acknowledgements

This work was supported by the Innovation Training Project for Students in NUAA (No.2016C-X0010-129) and the Key Laboratory of Aircraft Environment Control and Life Support (NUAA), Ministry of Industry and Information Technology.