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].

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{E}}_0}{\xi ^2}\mathit{\boldsymbol{a}}{{\left( \xi \right)}_{,\xi \xi }} + \left( {{\mathit{\boldsymbol{E}}_0} + {\mathit{\boldsymbol{E}}_1}^{\rm{T}} - } \right.}\\ {\left. {{\mathit{\boldsymbol{E}}_1}} \right)\xi \mathit{\boldsymbol{a}}{{\left( \xi \right)}_{,\xi }} - \left( {{\mathit{\boldsymbol{E}}_2} + {\rm{i}}\omega {\mathit{\boldsymbol{M}}_0}{\xi ^2}} \right)\mathit{\boldsymbol{a}}\left( \xi \right) = 0} \end{array} $

(1)

where E₀, E₁, E₂, and M₀ are 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) can be expressed as^[14]

$ \begin{array}{*{20}{c}} {\int\limits_S {\mathit{\boldsymbol{N}}{{\left( s \right)}^{\rm{T}}}\left( { - {{\bar q}_n}\left( {{\xi _1},s} \right)} \right){\tau ^{{\xi _1}}}{\rm{d}}s} = }\\ {{\mathit{\boldsymbol{E}}_0}{\xi _1}\mathit{\boldsymbol{a}}\left( {{\xi _1}} \right){,_\xi } + {\mathit{\boldsymbol{E}}_1}^{\rm{T}}\mathit{\boldsymbol{a}}\left( {{\xi _1}} \right)} \end{array} $

(2)

where s is the circumferential coordinate, N(s) the shape function matrix, τ^ξ₁ the scaling factors from scaled coordinate to the Cartesian coordinate when ξ=ξ₁=1, and q_n the outer normal heat flux on the boundary.

In SBFEM, the internal normal flow rate Q(ξ) through curved surface S(ξ) can be expressed as ^[14]

$ \mathit{\boldsymbol{Q}}\left( \xi \right) = \xi {\mathit{\boldsymbol{E}}_0}\mathit{\boldsymbol{a}}{\left( \xi \right)_{,\xi }} + {\mathit{\boldsymbol{E}}_1}^{\rm{T}}\mathit{\boldsymbol{a}}\left( \xi \right) $

(3)

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

$ \mathit{\boldsymbol{Q}}\left( \xi \right) = \mathit{\boldsymbol{K}}\left( {\xi ,\omega } \right)\mathit{\boldsymbol{a}}\left( \xi \right) $

(4)

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

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}}\left( {\xi ,\omega } \right)\mathit{\boldsymbol{a}}\left( \xi \right) = \xi {\mathit{\boldsymbol{E}}_0}\mathit{\boldsymbol{a}}{{\left( \xi \right)}_{,\xi }} + }\\ {{\mathit{\boldsymbol{E}}_1}^{\rm{T}}\mathit{\boldsymbol{a}}\left( \xi \right)} \end{array} $

(5)

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

$ \begin{array}{*{20}{c}} {\left( {\mathit{\boldsymbol{K}}\left( {\xi ,\omega } \right) - {\mathit{\boldsymbol{E}}_1}} \right){\mathit{\boldsymbol{E}}_0}^{ - 1}\left( {\mathit{\boldsymbol{K}}\left( {\xi ,\omega } \right) - } \right.}\\ {\left. {{\mathit{\boldsymbol{E}}_1}^{\rm{T}}} \right) - {\mathit{\boldsymbol{E}}_2} + 2\left( {{\rm{i}}\omega {\xi ^2}} \right)\mathit{\boldsymbol{K}}{{\left( {\xi ,\omega } \right)}_{,\left( {{\rm{i}}\omega {\xi ^2}} \right)}} - }\\ {{\rm{i}}\omega {\xi ^2}{\mathit{\boldsymbol{M}}_0} = 0} \end{array} $

(6)

where K(ξ, ω) must be a function of iωξ².By changing the independent variable to x=iωξ², and dynamic stiffness S(x)=K(ξ, ω), the scaled boundary finite-element equation in dynamic stiffness can be rewritten as

$ \begin{array}{*{20}{c}} {\left( {\mathit{\boldsymbol{S}}\left( x \right) - {\mathit{\boldsymbol{E}}_1}} \right){{\left( {{\mathit{\boldsymbol{E}}_0}} \right)}^{ - 1}}\left( {\mathit{\boldsymbol{S}}\left( x \right) - } \right.}\\ {\left. {{\mathit{\boldsymbol{E}}_1}^{\rm{T}}} \right) - {\mathit{\boldsymbol{E}}_2} + 2x{{\left[ {\mathit{\boldsymbol{S}}\left( x \right)} \right]}_{,x}} - x{\mathit{\boldsymbol{M}}_0} = 0} \end{array} $

(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 with x=iωξ² as independent variable. The dynamic stiffness S(x) is usually not calculated. The static-heat conduction matrix K and mass thermal capacity matrix M are used, which follow from the low-frequency expansion of S(x)^[14]

$ \mathit{\boldsymbol{S}}\left( x \right) = \mathit{\boldsymbol{K}} + x\mathit{\boldsymbol{M}} $

(8)

Substituting Eq. (8) into Eq. (7) will obtain the expression of K and M, and then S(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

$ \mathit{\boldsymbol{S}}\left( x \right) = \mathit{\boldsymbol{K}} + x\mathit{\boldsymbol{M}} - {x^2}{\left[ {{\mathit{\boldsymbol{S}}^{\left( 1 \right)}}\left( x \right)} \right]^{ - 1}} $

(9)

where the fraction x²[S⁽¹⁾(x)]^－1 is the residual part of the low frequency expansion. And S⁽¹⁾(x) can be written as

$ {\mathit{\boldsymbol{S}}^{\left( 1 \right)}}\left( x \right) = {\mathit{\boldsymbol{S}}_0}^{\left( 1 \right)} + x{\mathit{\boldsymbol{S}}_1}^{\left( 1 \right)} $

(10)

where S₀⁽¹⁾ and S₁⁽¹⁾ are coefficient matrices of the first-order terms. The second-order term in Eq.(10) is ignored. Substituting Eq.(9) into Eq.(7), and equating terms corresponding to x⁰ and x¹ to zero, respectively, we can obtain the static heat conduction matrix K and mass thermal capacity M. More detail information can be found in Ref ^[14]. Considering term of x² and the zero coefficient matrix, the following equation can be established.

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{S}}_1}^{\left( 1 \right)}\left( x \right){\mathit{\boldsymbol{c}}^{\left( 1 \right)}}{\mathit{\boldsymbol{S}}^{\left( 1 \right)}}\left( x \right) - {\mathit{\boldsymbol{S}}^{\left( 1 \right)}}\left( x \right){{\left[ {{\mathit{\boldsymbol{b}}_0}^{\left( 1 \right)}} \right]}^{\rm{T}}} - }\\ {{\mathit{\boldsymbol{b}}_0}^{\left( 1 \right)}{\mathit{\boldsymbol{S}}^{\left( 1 \right)}}\left( x \right) - x{\mathit{\boldsymbol{S}}^{\left( 1 \right)}}\left( x \right){{\left[ {{\mathit{\boldsymbol{b}}_1}^{\left( 1 \right)}} \right]}^{\rm{T}}} - }\\ {x{\mathit{\boldsymbol{b}}_1}^{\left( 1 \right)}{\mathit{\boldsymbol{S}}^{\left( 1 \right)}}\left( x \right) + 2x{\mathit{\boldsymbol{S}}^{\left( 1 \right)}}{{\left( x \right)}_{,x}} + {x^2}{\mathit{\boldsymbol{a}}^{\left( 1 \right)}}} \end{array} $

(11)

where

$ {\mathit{\boldsymbol{a}}^{\left( 1 \right)}} = {\mathit{\boldsymbol{E}}_0}^{ - 1} $

$ \mathit{\boldsymbol{b}}_0^{\left( 1 \right)} = {\mathit{\boldsymbol{E}}_0}^{ - 1}\left( {\mathit{\boldsymbol{K}} - {\mathit{\boldsymbol{E}}_1}^{\rm{T}}} \right) + 2\mathit{\boldsymbol{I}} $

$ \mathit{\boldsymbol{b}}_1^{\left( 1 \right)} = {\mathit{\boldsymbol{E}}_0}^{ - 1}\mathit{\boldsymbol{M}} $

$ {\mathit{\boldsymbol{c}}^{\left( 1 \right)}} = \mathit{\boldsymbol{M}}{\mathit{\boldsymbol{E}}_0}^{ - 1}\mathit{\boldsymbol{M}} $

Substituting Eq. (10) into Eq. (11) leads to equations for S₀⁽¹⁾ and S₁⁽¹⁾. Setting the terms corresponding to x⁰ equal to zero can obtain the equation for S₀⁽¹⁾, which is

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{b}}_0^{\left( 1 \right)}\mathit{\boldsymbol{S}}_0^{\left( 1 \right)} - \mathit{\boldsymbol{S}}_0^{\left( 1 \right)}{{\left[ {\mathit{\boldsymbol{b}}_0^{\left( 1 \right)}} \right]}^{\rm{T}}} + }\\ {\mathit{\boldsymbol{S}}_0^{\left( 1 \right)}{\mathit{\boldsymbol{c}}^{\left( 1 \right)}}\mathit{\boldsymbol{S}}_0^{\left( 1 \right)} = 0} \end{array} $

(12)

Eq. (12) can be transformed into a Lyapunov equation for [S₀⁽¹⁾]^-1 by pre-multiplying and post-multiplying with [S₀⁽¹⁾]^-1, which is

$ {\left[ {{\mathit{\boldsymbol{S}}_0}^{\left( 1 \right)}} \right]^{ - 1}}\mathit{\boldsymbol{b}}_0^{\left( 1 \right)} + {\left[ {\mathit{\boldsymbol{b}}_0^{\left( 1 \right)}} \right]^{\rm{T}}}{\left[ {{\mathit{\boldsymbol{S}}_0}^{\left( 1 \right)}} \right]^{ - 1}} = {\mathit{\boldsymbol{c}}^{\left( 1 \right)}} $

(13)

After solving [S₀⁽¹⁾]^-1 from the Lyapunov Eq. (13), S₀⁽¹⁾ can be calculated.

Equating terms corresponding to x¹ as zero yields

$ \begin{array}{*{20}{c}} {\left( { - \mathit{\boldsymbol{b}}_0^{\left( 1 \right)} + \mathit{\boldsymbol{S}}_0^{\left( 1 \right)}{\mathit{\boldsymbol{c}}^{\left( 1 \right)}}} \right)\mathit{\boldsymbol{S}}_1^{\left( 1 \right)} + }\\ {\mathit{\boldsymbol{S}}_1^{\left( 1 \right)}\left( { - {{\left[ {\mathit{\boldsymbol{b}}_0^{\left( 1 \right)}} \right]}^{\rm{T}}} + {\mathit{\boldsymbol{c}}^{\left( 1 \right)}}\mathit{\boldsymbol{S}}_0^{\left( 1 \right)}} \right) + 2\mathit{\boldsymbol{S}}_1^{\left( 1 \right)} = }\\ {\mathit{\boldsymbol{b}}_1^{\left( 1 \right)}\mathit{\boldsymbol{S}}_0^{\left( 1 \right)} + \mathit{\boldsymbol{S}}_0^{\left( 1 \right)}{{\left[ {\mathit{\boldsymbol{b}}_1^{\left( 1 \right)}} \right]}^{\rm{T}}}} \end{array} $

(14)

From the Lyapunov Eq. (14), S₁⁽¹⁾ can be obtained.

After determining S₀⁽¹⁾ and S₁⁽¹⁾, we can establish an increment-dimensional SBFE equation. Using Eqs. (4) and (9), the relationship between the nodal flux Q(ξ) and nodal temperature vector a(ξ) can be written as

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}}\left( \xi \right) = \mathit{\boldsymbol{Ka}}\left( \xi \right) + x\mathit{\boldsymbol{Ma}}\left( \xi \right) - }\\ {x{\mathit{\boldsymbol{a}}^{\left( 1 \right)}}\left( \xi \right)} \end{array} $

(15)

Here, the auxiliary variable a⁽¹⁾(ξ) is defined in

$ x\mathit{\boldsymbol{a}}\left( \xi \right) = {\mathit{\boldsymbol{S}}^{\left( 1 \right)}}\left( x \right){\mathit{\boldsymbol{a}}^{\left( 1 \right)}}\left( \xi \right) $

(16)

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

$ {\xi ^2}{\mathit{\boldsymbol{M}}_{\rm{h}}}\mathit{\boldsymbol{\dot y}}\left( \xi \right) + {\mathit{\boldsymbol{K}}_{\rm{h}}}\mathit{\boldsymbol{y}}\left( \xi \right) = \mathit{\boldsymbol{F}}\left( \xi \right) $

(17)

where

$ {\mathit{\boldsymbol{K}}_{\rm{h}}} = {\rm{diag}}\left( {\mathit{\boldsymbol{K}},{S_0}^{\left( 1 \right)}} \right) $

(18a)

$ {\mathit{\boldsymbol{M}}_{\rm{h}}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{M}}&{ - \mathit{\boldsymbol{I}}}\\ { - \mathit{\boldsymbol{I}}}&{{S_1}^{\left( 1 \right)}} \end{array}} \right] $

(18b)

$ \mathit{\boldsymbol{y}}\left( \xi \right) = {\left[ {\mathit{\boldsymbol{a}}\left( \xi \right);{\mathit{\boldsymbol{a}}^{\left( 1 \right)}}\left( \xi \right)} \right]^{\rm{T}}} $

(18c)

$ \mathit{\boldsymbol{F}}\left( \xi \right) = {\left[ {\mathit{\boldsymbol{Q}}\left( \xi \right);\mathit{\boldsymbol{0}}} \right]^{\rm{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

For comparison, we chose the same example as in Ref.[14].As shown in Fig. 1, the left and bottom boundaries are insulated, while the top and right boundaries are located in a convection environment. We assume that a=b=1 m, $\alpha = \frac{k}{{\rho {c_p}}}$=1 m²/s, $\frac{h}{k}$=0.5 m^－1, T_∞=0 ℃, and the initial temperature T₀ is set to 100 ℃.

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.

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. 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.

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 at t=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.

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.