Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2022.019687

ARTICLE

A Novel Meshfree Analysis of Transient Heat Conduction Problems Using RRKPM

Hongfen Gao1 and Gaofeng Wei2,3,*

1College of Intelligent Engineering, Shandong Management University, Jinan, 250357, China
2School of Mechanical and Automotive Engineering, Qilu University of Technology (Shandong Academy of Sciences), Jinan, 250353, China
3Shandong Institute of Mechanical Design and Research, Jinan, 250353, China
*Corresponding Author: Gaofeng Wei. Email: wgf@qlu.edu.cn
Received: 08 October 2021; Accepted: 18 November 2021

1  Introduction

Many practical engineering structures run in high temperature, such as steam turbines, high-speed diesel engines and nuclear power plants, etc. The temperature field can change the properties of material structure. Therefore, it is an important subject to study the THCP of the structure in the condition of being heated [1,2].

Numerical simulation is an important analysis tool to research the THCP [3]. Finite element method (FEM) is one of the main numerical simulation methods. Many complex and difficult mechanical problems can be solved by the FEM, and valuable results can be derived [4,5]. However, due to the limitation of the correlation conditions between elements in the FEM, it is difficult to deal with the discontinuous problem in practical engineering problems, such as the formation of cracks and their mechanical behavior, the discontinuity in jointed rock mass and the crack propagation with moving boundary [6,7]. In order to improve the restriction of correlation conditions between elements, many novel methods have been proposed in recent years, such as meshfree (or meshless, element-free) method [8,9], numerical manifold method (NMM) [10–12], boundary element method (BEM) [13,14], and numerical method based on least square method [15–17], etc.

The meshfree approximating technique is adopted in the meshfree method, which makes the approximating function free from the constraints of elements and greatly simplifies the analysis and calculation of pretreatment and crack propagation. Meshfree method has attracted attention in mechanics and practical engineering, and been widely used in the study of the THCP. At present, many meshfree methods have been developed, such as smooth particle hydrodymics (Abbreviation: SPH, proposed by Lucyt and Gingold in 1977) [18,19], element-free Galerkin method (Abbreviation: EFGM, proposed by Belytschko in 1994) [20–22], meshfree local Petrov-Galerkin method (Abbreviation: MLPG, proposed by Atluri in 1998) [23–26], reproducing kernel particle method (Abbreviation: RKPM, proposed by Liu in 2005) [27–30], radial basis functions method (Abbreviation: RBF, proposed by Žilinskas in 2010) [31,32], complex variable meshfree manifold method (Abbreviation: CVMMM, proposed by Gao in 2010) [33], the finite point method (Abbreviation: FPM, proposed by Tatari in 2011) [34,35], Hermit-type reproducing kernel particle method (Abbreviation: Hermit-type RKPM, proposed by Ma in 2017) [36–39] and boundary integral equation method (Abbreviation: BIE, proposed by Mantegh in 2010) [40,41], etc.

Because of the advantages of simple form and fast calculation speed, the RKPM is one of the meshfree methods which are widely applied and researched [42–44]. The RKPM is first proposed based on the SPH and the integral reconstruction theory of functions. The RKPM solves the boundary inconsistency, and eliminates the tensile instability of the SPH method. The method has some advantages, such as variable time frequency characteristics and multi-resolution characteristics. Therefore, the RKPM has been widely used in many practical engineering problems, such as large deformation analysis problems, structural dynamics problems, micro-electromechanical system analysis problems, nonlinear problem of hyperelastic rubber materials, high-speed impact problems and so on [45–47].

However, different kernel functions have different effects on the calculating accuracy and computational stability in the analysis of solving the THCP. In order to improve the calculating accuracy and stability of the RKPM, the RBF is introduced into the RKPM, and the meshfree RRKPM is proposed in this paper. Meanwhile, the meshfree RRKPM is applied to the THCP, and the corresponding equations of the meshfree RRKPM for THCP are derived. The numerical results illustrated the effectiveness of the meshfree RRKPM for the THCP.

2  Construction of the Approximating Function of Meshfree RRKPM

The approximating function uh(z) of Meshfree RRKPM can be written as a combination of the RKPM constructed by n nodes in the local problem domain and the RBFs constructed by m terms.

uh(z)=∑j=1nSjn(z)cj+∑i=1mSim(z)ai,z=(x,y)T(1)

where cj and ai represent the undetermined coefficients, n represents the number of the local influence domain, m represents the number of RBFs, Sim and Sjn represent the RBFs and the reproducing kernel function (RKF), respectively.

The RBF Sim is the function of the distance ri from calculating point z to the node zi, ri=||z−zi||

Sim(z)=(1−riδ)6(3+18riδ+35ri2δ2)(2)

with δ denoting the scaling parameter.

The RKF Sjn can be expressed as

Sjn(z)=C(z,zj)w(z−zj)Sjn(zj)ΔVj(3)

where Sjn(zj) is the parameter of node zj, ΔVj is the area or volume of domain of influence, w is kernel function.

C(z,zj)=bT(z)p(z−zj)(4)

The coefficient matrix b(z) can be given by

b(z)=[b1(z),b2(z),⋯,bn(z)]T(5)

The polynomial basis function p(z−zj) is expressed as

pT(z−zj)=[1,x−xj,y−yj,(x−xj)2,(x−xj)(y−yj),(y−yj)2,⋯](6)

The Eq. (1) can be rewritten as the following form:

uh(z)=∑j=1nSjn(z)cj+∑i=1mSim(z)ai=∑I=1MSI(z)aI(z)=S(z)a(z)(M=m+n)(7)

in which SI(z) is basis function, aI(z) is corresponding coefficient, given by

S(z)=[S1m(z),⋯,Smm(z),S1n(z),⋯,Snn(z)](8)

a(z)=[a1,⋯,am,c1,⋯,cn]T(9)

The approximating function uh(z) can be locally approximated in the neighborhood of calculating point z

uh(z,z¯)=∑I=1MSI(z¯)aI(z)=S(z¯)a(z)(10)

in which z¯ is point of the neighborhood in calculating point z.

The weighted least squares method is used to obtain approximating functions uh(z) accurately in this paper. The weighted least squares function J is defined as

J=∑K=1Nw(z−zK)[uh(z,zK)−u(zK)]2=∑K=1Nw(z−zK)[∑I=1MSI(zK)aI(z)−u(zK)]2(11)

where w(z−zK) is weighted function in the domain of influence, zK(K=1,2,⋯,N) are the nodes in the domain of influence.

The Eq. (11) can be rewritten as the following form:

J=(Sa−u)TW(z)(Sa−u)(12)

where

uT=(u1,u2,⋯,uN)(13)

S=[S1(z1)S2(z1)⋯SM(z1)S1(z2)S2(z2)⋯SM(z2)⋮⋮⋱⋮S1(zN)S2(zN)⋯SM(zN)](14)

W(x)=[w(z−z1)0⋯00w(z−z2)⋯0⋮⋮⋱⋮00⋯w(z−zN)](15)

Let J take the minimum, that is

∂J∂a=0(16)

The following form can be obtained

A(z)a(z)=B(z)u(17)

where

A(z)=STW(z)S(18)

B(z)=STW(z)(19)

The Eq. (17) can be given as

a(z)=A−1(z)B(z)u(20)

Substituting Eq. (20) into Eq. (7), the approximating function uh(z) is obtained

uh(z)=Φ(z)u=∑K=1NΦK(z)uK(21)

in which shaped function Φ(z) is expressed as

Φ(z)=(Φ1(z),Φ2(z),⋯,ΦN(z))=ST(z)A−1(z)B(z)(22)

3  Governing Equation of the THCP for Meshfree RRKPM

3.1 Fundamental Equations for the THCP

From the theory of transient heat conduction, the differential equation of THCP in orthotropic plane can be expressed as

∂∂x(kx∂T(z,t)∂x)+∂∂y(ky∂T(z,t)∂y)+qv=ρcp∂T(z,t)∂t(23)

where T(z,t) represents transient temperature, t denotes transient heat transfer time, kx and ky represent thermal conductivities of material in plane principal axes, ρ represents density of material, cp represents constant pressure specific heat and qv is internal heat source intensity.

The material is assumed to be isotropic with kx=ky=k, Eq. (23) can be simplified as

k(∂2T(z,t)∂x2+∂2T(z,t)∂y2)+qv=ρcp∂T(z,t)∂t(24)

or

∇2T(z,t)+qv/k=(1/αT)∂T(z,t)/∂t(25)

in which αT=k/ρcp represents thermal diffusivity, ∇2 denotes Laplace operator.

∇2=∂2/∂x2+∂2/∂y2(26)

Assuming no heat source, the Eq. (25) is rewritten as the Fourier equation:

∇2T=(1/αT)∂T(z,t)/∂t(27)

In order to obtain the unique solution of the THCP, boundary conditions and initial conditions must be applied. There are three kinds of boundary conditions as follows:

1.    The first kind of boundary condition is that the temperature on the boundary is known, and the formula is

T(z,t)|Γ1=T¯(z,t)(28)

T(x,y,t)|Γ1=f(x,y,t)(29)

(2)   The second kind of boundary condition is that the heat flux density on the boundary is known. Because the direction of the heat flux is the exterior normal direction of the boundary, the formula is as following:

−(∂T(z,t)/∂n)|Γ2=h¯(30)

−(∂T(x,y,t)/∂n)|Γ2=h(x,y,t)(31)

(3)   The third kind of boundary condition is that the convection or radiant heat transfer on the boundary is known. For convection heat transfer conditions

−k(∂T(z,t)/∂n)|Γ3=g(T(z,t)−Tq(z,t))|Γ3(32)

For radiant heat transfer conditions, it can be written as

−k(∂T(z,t)∂n)|Γ3=εfσ0(T4(z,t)−Tr4(z,t))|Γ3(33)

with ε representing blackness coefficient, f denoting shaped factor, σ0 being Stefan-Bolzman constant and Tr(z,t) representing temperature of radiant source.

The initial condition is the known value of the temperature at the beginning of the heat transfer process, and the formula is

T(z,t)|t=0=0(34)

or

T(z,t)|t=0=T0(x,y)(35)

From the heat conduction equation and boundary conditions, it can be seen that there is only one partial differential equation (PDE) and only one temperature as an unknown variable, therefore, the THCP is actually solving the PDE.

3.2 Integral Weak Form of the THCP

In a certain instantaneous state, T(z,t) and ∂T(z,t)∂t can be considered as deterministic functions of plane coordinates. The THCP can be transformed to the elliptic equation of boundary value problem, and the formula is

Π=∫Ω[T(z,t)(ρc∂T(z,t)∂t−qv)]dΩ+∫Ω[12k1(∂T(z,t)∂x)2+12k2(∂T(z,t)∂y)2]dΩ+∫Γ2T(z,t)⋅q¯dΓ+∫Γ3h(T2(z,t)2−T(z,t)⋅Tα(z,t))dΓ(36)

The field function, which makes the variational of the function Π equal zero, is the solution which satisfies the governing differential Eq. (23) and boundary condition of this problem.

Taking the first kind of boundary problem as an example, the temperature function T(z,t) must satisfy the essential boundary condition (29), which is T(z,t)−T¯(z,t)=0 on the boundary Γ1, so the function Π is conditional function. Introducing the penalty function method into the essential boundary conditions (28)–(32), and another modified function Π∗ can be constructed as

Π∗=Π+12∫Γ1(T(z,t)−T¯(z,t))T⋅β⋅(T(z,t)−T¯(z,t))dΓ(37)

with β representing penalty factor of generally 103{∼}105 in the THCP. After introducing essential boundary conditions, the conditional stationary value problem of original function Π transforms into the unconditional stationary value problem of modified function Π∗. The first variation of the stationary condition for modified function Π∗ equals zero.

δΠ∗=δΠ+∫Γ1δ(T(z,t)−T¯(z,t))T⋅β⋅(T(z,t)−T¯(z,t))dΓ=0(38)

Substituting Eq. (36) into Eq. (38), the integral weak form of the THCP is

∫ΩδT⋅ρc⋅∂T∂tdΩ+∫Ωδ(LT)Tk^(LT)dΩ−∫ΩδT⋅qvdΩ−∫Γ2δT⋅h¯dΓ−∫Γ3δT⋅g(Tα−T)dΓ+∫Γ1δT⋅β⋅TdΓ−∫Γ1δT⋅β⋅T¯dΓ=0(39)

where

L(⋅)=(∂∂x∂∂y)(⋅)(40)

k^=[k100k2](41)

3.3 The Meshfree RRKPM for the THCP

The THCP is a function in space domain Ω and in time domain t, and these two domains are not coupled. Therefore, the meshfree RRKPM and finite difference method (FDM) can be used to solve the problem, that is, the THCP is solved by the meshfree RRKPM in space domain and by the FDM in time domain. Firstly, the domain Ω is discretized into a finite number of nodes, and then the temperature of any point in the domain at any time t is approximated by the node temperature TI(t) in its influence domain.

T(t)=TI(zI,t)(42)

It should be noted that T(zI,t) of any field point z in the domain is a scalar at any time, so the temperature can be given as

T(z,t)=∑J=1NΦJ(z)TJ(z,t)=Φ(z)⋅T(t)(43)

in which Φ(z) represents a shaped function vector, which is just a function in the space domain.

T=(T1(t),T2(t),⋯,TN(t))T(44)

and

∂T(z,t)∂t=∂∂t∑J=1NΦJ(z)⋅TJ(t)=∑J=1NΦJ(z)⋅∂TJ(t)∂t=Φ(z)⋅T˙(t)(45)

T˙=(∂T1(t)∂t,∂T2(t)∂t,⋯,∂TN(t)∂t)T(46)

LT(z,t)=∑J=1N(∂∂x∂∂y)⋅[ΦJ(z)⋅TJ(t)]=∑J=1NBJ(z)⋅TJ(t)=B(z)⋅T(t)(47)

with

B(z)=(B1(z),B2(z),⋯,BN(z))(48)

BJ(z)=(∂∂x∂∂y)⋅ΦJ(z)=[ΦJ,x(z)ΦJ,y(z)](49)

Substituting Eqs. (43), (45) and (47) into Eq. (39), the following form can be obtained:

∫Ωδ[Φ(z)T]T⋅ρc⋅[Φ(z)T˙]dΩ+∫Ωδ(B(z)T)Tk^(B(z)T)dΩ−∫Ωδ[Φ(z)T]T⋅qvdΩ−∫Γ2δ[Φ(z)T]T⋅h¯dΓ−∫Γ3δ[Φ(z)T]T⋅g⋅TαdΓ+∫Γ3δ[Φ(z)T]T⋅g⋅[Φ(z)T]dΓ+∫Γ1δ[Φ(z)T]T⋅β⋅[Φ(z)T]dΓ−∫Γ1δ[Φ(z)T]T⋅β⋅T¯dΓ=0(50)

In order to solve the discrete system equations, the integral Eq. (50) is discussed separately below.

The first term of Eq. (50) is

∫Ωδ[Φ(z)T]T⋅ρc⋅[Φ(z)T˙]dΩ=δTT[∫ΩΦT(z)⋅ρc⋅Φ(z)dΩ]⋅T˙=δTT⋅C⋅T˙(51)

where C represents heat capacity matrix, and can be expressed as

C=[C11C12⋯C1NtC21C22⋯C2Nt⋮⋮⋱⋮CNt1CNt2⋯CNtNt](52)

CIJ=∫ΩΦIT(z)⋅ρc⋅ΦJ(z)dΩ(53)

The second term of Eq. (50) is

∫Ωδ(B(z)T)Tk^(B(z)T)dΩ=δTT[∫ΩBT(z)⋅k^⋅B(z)dΩ]⋅T=δTT⋅K⋅T(54)

where K represents heat conduction matrix

K=[K11K12⋯K1NtK21K22⋯K2Nt⋮⋮⋱⋮KNt1KNt2⋯KNtNt](55)

KIJ=∫ΩBI(z)⋅k~IJ⋅BJ(z)dΩ(56)

The third term of Eq. (50) is

∫Ωδ[Φ(z)T]T⋅qvdΩ=δTT∫ΩΦT(z)⋅qvdΩ=δTT⋅F(1)(57)

F(1)=(f1(1),f2(1),⋯,fNt(1))T(58)

fI(1)=∫ΩΦIT(z)⋅qvdΩ(59)

The fourth term of Eq. (50) is

∫Γ2δ[Φ(z)T]T⋅h¯dΓ=δTT∫Γ2Φ~T(z)⋅h¯dΓ=δTT⋅F(2)(60)

F(2)=(f1(2),f2(2),⋯,fNt(2))T(61)

fI(2)=∫Γ2ΦIT(z)⋅h¯dΓ(62)

The fifth term of Eq. (50) is

∫Γ3δ[Φ(z)T]T⋅g⋅TβdΓ=δTT∫Γ3ΦT(z)⋅g⋅TβdΓ=δTT⋅F(3)(63)

F(3)=(f1(3),f2(3),⋯,fNt(3))T(64)

fI(3)=∫Γ3ΦIT(z)⋅g⋅TβdΓ(65)

where F(1), F(2) and F(3) represent thermal load vectors known as heat source, given heat flow and heat exchange, respectively.

The sixth term of Eq. (50) is

∫Γ3δ[Φ(z)T]T⋅g⋅[Φ(z)T]dΓ=δTT∫Γ3ΦT(z)⋅g⋅Φ(z)dΓ⋅T=δTT⋅G⋅T(66)

where

G=[G11G12⋯G1NtG21G22⋯G2Nt⋮⋮⋱⋮GNt1GNt2⋯GNtNt](67)

GIJ=∫Γ3ΦI(z)⋅g⋅ΦJ(z)dΓ(68)

The seventh term of Eq. (50) is

∫Γ1δ[Φ(z)T]T⋅β⋅[Φ(z)T]dΓ=δTT∫Γ1ΦT(z)⋅β⋅Φ(z)dΓ⋅T=δTT⋅Kβ⋅T(69)

where

Kβ=[K11βK12β⋯K1NtβK21βK22β⋯K2Ntβ⋮⋮⋱⋮KNt1βKNt2β⋯KNtNtβ](70)

KIJβ=∫Γ1ΦI(z)⋅β⋅ΦJ(z)dΓ(71)

The eighth term of Eq. (50) is

∫Γ1δ[Φ(z)T]T⋅β⋅T¯dΓ=δTT∫Γ1ΦT(z)⋅β⋅T¯dΓ=δTT⋅Fβ(72)

where

Fβ=(f1β,f2β,⋯,fNtβ)T(73)

fIβ=∫Γ1ΦIT(z)⋅β⋅T¯dΓ(74)

Substituting Eqs. (51), (54), (57), (60), (63), (66), (69) and (72) into Eq. (50), and the following form can be given:

δTT(CT˙+KT+GT+KβT−F(1)−F(2)−F(3)−Fβ)=0(75)

From the arbitrariness of δTT, the final ordinary differential equations (ODE) can be given as

CT˙+K^T−F^=0(76)

where

K^=K+H+Kβ(77)

F^=F(1)+F(2)+F(3)+Fβ(78)

The PDE problems of the THCP have been discretized into initial value problems of the ODE with nodal temperature T(t) in space domain Ω.

The above is the meshfree RRKPM for the THCP.

4  Time Integral Scheme

Eq. (76) is the linear ODE with time t being independent variable. In order to discretize the time domain of the ODE, the traditional two-point difference method is used in this paper.

In space domain Ω, the temperature vector T is a function of time t, and can be divided into several elements. In any element, T(z,t) is given as

T(z,t)≈T¯(z,t)=∑Ni(z)Ti(t)(79)

where Ti(t)=T(ti) represents nodal temperature vector at time ti. The interpolating function Ni(z) takes as same form for each component of the vector T(z, t).

When the ODE only contains the first-order derivative of time t, the interpolate function is a linear polynomial, and the two-point first-order interpolation can be used.

For the time interval Δt, T(z, t) can be obtained by interpolation of node values Tn(t) and Tn+1(t) in an interval

T(z,t)=Nn(z)Tn(t)+Nn+1(z)Tn+1(t)(80)

The first-order derivative of T(z, t) is

T˙(z,t)=N˙n(z)Tn(t)+N˙n+1(z)Tn+1(t)(81)

The interpolate function and the first-order derivative can be expressed by the local variable λ

{λ=tΔt(0≤λ≤1)Nn=1−λ,N˙n=−1ΔtNn+1=λ,N˙n+1=1Δt(82)

Using approximating interpolation of Eqs. (80) and (81), the Eq. (76) inevitably produces residual in a time interval Δt. A weighted residual expression is derived as

∫01w[C(N˙nTn+N˙n+1Tn+1)+K^(NnTn+Nn+1Tn+1)−F^]dλ=0(83)

Substituting Eq. (82) into Eq. (83), the residual relation of two time intervals can be given as

(K^∫01wλdλ+C∫01w1Δtdλ)Tn+1+(K^∫01w(1−λ)dλ−C∫01w1Δtdλ)Tn−∫01wF^dλ=0(84)

Eq. (84) can be seen as a general form applicable to any weighted function.

(CΔt+K^ζ)Tn+1+[−CΔt+K^(1−ζ)]Tn=F¯(85)

F¯ is supposed to use the same interpolation as the unknown temperature function T(z, t).

F¯=F^n+1ζ+F^n(1−ζ)(86)

Substituting Eq. (86) into Eq. (85), the following form can be expressed

(CΔt+ζK^n+1)Tn+1=(CΔt−(1−ζ)K^n)Tn+ζF^n+1+(1−ζ)F^n(87)

here

ζ=∫01wλdλ∫01wdλ(88)

The above is the time difference scheme for the THCP.

5  Numerical Examples

5.1 Transient Temperature Field of the THCP in Rectangular Domain

The governing equation of the THCP in the rectangular domain is

∂T∂t−∂2T∂x2−∂2T∂y2=(1+t2)T+(2π2−t2−2)e−tsin⁡(πx)cos⁡(πy),(x,y)∈(0,1)(89)

According to the boundary conditions, written as

T(x,y,t)|x=0=T(x,y,t)|x=1=0(90)

T(x,y,t)|y=0=T(x,y,t)|y=1=e−tsin⁡(πx)(91)

and the initial condition, given by

T(x,y,t)|t=0=sin⁡(πx)cos⁡(πy)(92)

The analytical solution of the THCP can be obtained as

T(x,y,t)=e−tsin⁡(πx)cos⁡(πy)(93)

As shown in Fig. 1, 11 × 11 nodes are uniformly distributed in the rectangular THCP domain Ω, and time interval Δt=0.001s. The penalty factor is taken as α=1.0×108 and the scaling parameter is taken as δ=2.0. The regular quadrilateral background mesh is applied to the governing equation of THCP, and 4×4 Gauss integral scheme is used.

In order to discuss the influence of kernel functions on the calculation accuracy and stability, the kernel function is taken as the following two forms:

w1(ri)={1−6ri2+8ri3−3ri4ri≤10ri>1(94)

w2(ri)={2/3−4ri2+4ri3ri≤1/24/3−4ri+4ri2−4ri3/31/2<ri≤10ri>1(95)

Figure 1: Distribution of nodes in rectangular domain for the THCP

In order to illustrate the validity of the proposed method, the temperature in x = 0.5 and at t = 0.1 s is calculated by the analytical solution, RKPM and RRKPM, respectively. Fig. 2 gives the comparison of the temperature between two methods using different kernel functions which the kernel function of Eq. (94) is defined as the kernel function 1 and the kernel function of (95) is defined as the kernel function 2. The relative error is defined as

erelativeerror=|Tanalytical−Tnumerical||Tanalytical|(96)