Finite Difference-Peridynamic Differential Operator for Solving Transient Heat Conduction Problems

1  Introduction

Transient heat conduction problems are prevalent in petroleum, chemical, metallurgy, and many other fields. Consequently, effective numerical methods for studying these problems are crucial for practical engineering applications.

Currently, numerical methods for solving transient heat conduction equations are broadly classified into two categories: mesh-based methods, including finite difference method (FDM) [1,2], finite element method (FEM) [3,4], Finite Volume Method [5,6], and Boundary Element Method (BEM) [7]; and meshless methods, such as the generalized finite difference method [8,9], smoothed particle hydrodynamics method (SPH) [10], meshless local Petrov-Galerkin method (MLPG) [11–14], meshless local radial basis function-based differential quadrature (RBF-DQ) [15], peridynamics (PD) [16,17], and peridynamic differential operator (PDDO) [18,19]. Based on time discretization, these methods can be divided into explicit and implicit schemes. In the explicit scheme, unknown quantities are explicitly given in terms of known quantities, offering the advantage of simpler programming but suffering from strict stability conditions. Conversely, in the implicit scheme, unknown quantities cannot be explicitly expressed, making the solution computationally challenging yet allowing for larger time step sizes due to increased stability. Consequently, the implicit scheme is often favored in software.

The PDDO method, a recent advancement based on PD theory [20], is a nonlocal differential operator that bridges local partial derivatives and nonlocal integrals using Taylor series expansions and orthogonal function properties. Capable of solving differential equations and calculating derivatives from smooth functions or scattered data amidst discontinuities or singular points [21], it also features time nonlocality and generalized space-time nonlocality, unrestricted by order of space and time partial derivatives, proving its efficacy in practical applications. For instance, Dorduncu devised a nonlocal stress analysis model for functionally graded sandwich panels using PDDO [22], Gao et al. developed a nonlocal model for fluid flow and heat transfer coupling using PDDO [23], and Li et al. introduced a nonlocal model for steady-state thermoelastic analysis of functionally graded materials with PDDO [24]. Additionally, Li et al. compared the PDDO with the other nonlocal differential operators and proposed some improvements for PDDO [25–27].

In the previous work involving the PDDO method, both time and space derivatives were discretized using the PDDO method [28,29]. Given the relative complexity and computational expense of the PDDO method, it is essential to reduce its complexity. The FDM, being the oldest numerical method, offers the advantage of straightforward implementation. Consequently, it was chosen to discretize the time derivative. Furthermore, considering the capability of the PDDO method to handle complex regions and discontinuous problems in space, the PDDO method was utilized to discretize the spatial derivative.

In this study, the coupling of FDM with the PDDO method (FD-PDDO) has been developed to solve transient heat conduction equations. In order to establish both explicit and implicit methods, the time derivative is approximated using the forward difference formula, backward difference formula, and centered difference formula, respectively. As a result, the forward-in-time and PDDO in space (FT-PDDO) scheme, backward-in-time and PDDO in space (BT-PDDO) scheme, and central-in-time and PDDO in space (CT-PDDO) scheme are developed. The FT-PDDO scheme is explicit, while the BT-PDDO and CT-PDDO schemes are implicit. The stability and convergence of these new schemes are analyzed using the Fourier method and Taylor's theorem, respectively. The developed schemes are applied to solve one-dimensional and two-dimensional transient heat conduction problems, and their accuracy and validity are verified through comparison with the analytical solution.

2  Mathematical Model

The transient heat conduction equation can be expressed as:

{∂T(x,t)∂t=D∇2T(x,t)+Q(x,t)x∈Ω,T(x,0)=g(x)x∈Ω,T(x,t)=T¯(x,t)x∈Γ(1)

where T is the temperature; t is the time; x is the spatial variable; Q(x,t) is the heat source; T¯(x,t) is the temperature on the boundary Γ of the computational zone Ω; g(x) is the initial temperature; D=kpρc where ρ is the density; c is the specific heat capacity; kp is the thermal conductivity coefficient.

3  Numerical Method

3.1 Basic Theory of the FDM

The FDM is a mesh-based algorithm, and its basic idea is to transform the derivative into numerical differentiation. Taking the time derivative dTdt as an example and using the uniform grid shown in Fig. 1, the forward difference formula can be obtained from [1] as follows:

dTdt|t=tk≈Tk+1−TkΔt(2)

The backward difference formula in [1] is as follows:

dTdt|t=tk+1≈Tk+1−TkΔt(3)

The centered difference formula in [1] is as follows:

dTdt|t=tk+1/2≈Tk+1−TkΔt(4)

where Tk represents the numerical solution of T at the time point t=tk. The truncation error of the forward difference formula and the backward difference formula is O(Δt), while the truncation error of the centered difference formula is O((Δt)2). The truncation errors can be obtained using Taylor's theorem.

Figure 1: Uniform grid diagram

3.2 The PDDO Method

The PD theory is a nonlocal theory proposed by Silling et al. [20]. Point x interacts with Point x ′ within an interaction domain Hx as shown in Fig. 2. The relative position vector between these points is defined as ξ=x ′−x. Each point has its interaction domain (family). The interaction can be specified as δ=mΔx with m being an integer and Δx representing the grid spacing between the points. The interaction domain of points may have different sizes and shapes. The degree of nonlocal interaction between the points is specified by the weight function ω(ξ).

Figure 2: PD interaction domains for the discretized points x and x ′

Madenci et al. [29,30] proposed the PDDO method based on the PD theory. For the M-dimensional scalar function, the N-th order Taylor expansions of f(x ′)=f(x+ξ) are expressed as Eq. (5).

f(x+ξ)=∑n1=0N∑n2=0N−n1⋯∑nM=0N−n1−⋯nM−11n1!n2!⋯nM!ξ1n1ξ2n2⋯ξMnM∂n1+n2+⋯+nMf(x)∂x1n1∂x2n2⋯∂xMnM+R(N,x)(5)

where R(N,x) representing the remainder of the M-dimensional approximation.

Multiplying each term in Eq. (5) by PD functions gNp1p2⋯pM(ξ) and integrating over the domain of interaction Hx forms Eq. (6).

∂p1+p2+⋯+pMf(x)∂x1p1∂x2p2⋯∂xMpM=∫Hxf(x+ξ)gNp1p2⋯pM(ξ)dV(6)

where pi =0,1,⋯,N (i=1,2,⋯,M) represents the partial derivative with respect to variable xi and p1+p2+⋯+pM≤N. The function gNp1p2⋯pM(ξ) is a PD function constructed based on the orthogonality principle in the domain of interaction Hx of the point x, which must possess the orthogonality property of

1n1!n2!⋯nM!∫Hxξ1n1ξ2n2⋯ξMnMgNp1p2⋯pM(ξ)dV=δn1p1δn2p2⋯δnMpM(7)

where ni=0,⋯,N with n1+n2+⋯+nM≤N; δnipi represents the Kronecker symbol. The PD function gNp1p2⋯pM(ξ) can be constructed as [29,30]:

gNp1p2⋯pM(ξ)=∑q1=0N∑q2=0N−q1⋯∑q M=0N−q1−⋯qM−1aq1q2⋯qMp1p2⋯pMωq1q2⋯qM(|ξ|)ξ1q1ξ2q2⋯ξMqM(8)

where ωq1q2⋯qM(|ξ|) is the weight function related to term ξ1q1ξ2q2⋯ξMqM.

The unknown coefficient aq1q2⋯qMp1p2⋯pM in Eq. (8) can be obtained from the following equations:

∑q1=0N∑q2=0N−q1⋯∑q M=0N−q1−⋯qM−1∑q1+q2+⋯+qM=0NA(n1n2⋯nM)(q1q2⋯qM)aq1q2⋯qMp1p2⋯pM=bq1q2⋯qMp1p2⋯pM(9)

where qi=0,⋯,N with i=1,2,⋯,M and q1+q2+⋯+qM≤N. The coefficient matrix is constructed as follows:

A(n1n2⋯nM)(q1q2⋯qM)=∫Hxωq1q2⋯qM(|ξ|)ξ1n1+q1ξ2n2+q2⋯ξMnM+qMdV(10)

The right-hand term is as follows:

bq1q2⋯qMp1p2⋯pM=n1!n2!⋯nM!δn1p1δn2p2⋯δnMpM(11)

3.3 The FD-PDDO Method for Solving One-Dimensional Transient Heat Conduction Equations

Using the FDM to discrete the time derivative and the PDDO to discrete the spatial derivative, the new scheme of FD-PDDO is obtained for solving the one-dimensional transient heat conduction equation.

The FT-PDDO scheme is obtained by using the forward difference formula in time and the PDDO method in space as follows:

Tik+1 =Tik+DΔt∫Hx iTk(x+ξ)gN2(ξ)dξ =Tik+DΔt∑j=1N(i)Tk(xj)gN2(xj−xi)Δxj(12)

The BT-PDDO scheme is obtained by using a backward difference formula in time and the PDDO method in space, namely

Tik+1 =Tik+DΔt∫HxiTk+1(x+ξ)gN2(ξ)dξ =Tik+DΔt∑j=1N(i)Tk+1(xj)gN2(xj−xi)Δxj(13)

The CT-PDDO scheme is obtained by using the centered difference formula in time and the PDDO method in space, which is:

Tik+1 =Tik+DΔt∫HxiTk+1(x+ξ)+Tk(x+ξ)2gN2(ξ)dξ =Tik+DΔt2∑j=1N(i)Tk+1(xj)gN2(xj−xi)Δxj+DΔt2∑j=1N(i)Tk(xj)gN2(xj−xi)Δxj(14)

3.4 Stability and Convergence Analysis of the FD-PDDO Method for Solving One-Dimensional Transient Heat Conduction Equations

The weight function is taken as ω(ξ)=e−(2ξ/δ)2, where the interaction domain is δ=mΔx with m the integer parameter and Δx the uniform spatial step size. This study assumes the polynomial order as N=2, integer parameter as m=2 and m=3 (the range for m is suggested as N≤m≤N+2 [29,30]). Fig. 3 shows the PD interaction domains for the discretized Point x and x ′ in the one-dimensional case, with the interaction domain of δ=mΔx(m=2).

Figure 3: PD interaction domains for the discretized points x and x ′ in one-dimensional case

3.4.1 Stability and Convergence Analysis of the FT-PDDO Scheme

In the FT-PDDO scheme of Eq. (12), in the case of polynomial order N=2 and integer parameter m=2, the right-hand term ∑j=1N(i)Tk(xj)gN2(xj−xi)Δxj can be expanded as follows:

∑j=1N(i)Tk(xj)g22(xj−xi)Δxj =Ti−2k[ω(−2Δx)∑q=02aq2(−2Δx)q]Δx+Ti−1k[ω(−Δx)∑q=02aq2(−Δx)q]Δx +Tik[ω(0)∑q=02aq2(0)q]Δx +Ti+1k[ω(Δx)∑q=02aq2(Δx)q]Δx+Ti+2k[ω(2Δx)∑q=02aq2(2Δx)q]Δx(15)

where a02,a12,a22 can be obtained from the following equations:

(A00A01A02A10A11A12A20A21A22)(a02a12a22)=(b02b12b22)(16)

where A00=1.7724Δx; A01=A10=A12=A21=0; A02=A20=A11=0.8823Δx3; A22=1.3219Δx5; b02=b12=0,b22=2. The coefficients for the PD function are obtained as given in Eq. (17).

{a02=−1.12791Δx3a12=0a22=2.26591Δx5(17)

Therefore, the FT-PDDO scheme is obtained, as given in Eq. (18).

Tik+1=Tik+λ[0.1453Ti−2k+0.4186Ti−1k−1.1279Tik+0.4186Ti+1k+0.1453Ti+2k](18)

where λ=DΔt/Δx2.

The Fourier method is now utilized to show the stability [1]. Let Tik=ωkerxiI (where I=−1 and r is the positive constant). Substituting it into Eq. (18) yields the amplification factor as shown in Eq. (19).

κ=ωk+1/ωk=1+λ[−1.1279+0.2907cos⁡(2rΔx)+0.8372cos⁡(rΔx)](19)

The stability requirement is |κ|≤1. Thus, the stability condition is shown in Eq. (20).

λ≤1.1(20)

Taylor's theorem is now utilized to show the truncation error [1]. The truncation error caused by the PDDO method in space is shown only here. The term can be obtained from Eq. (18) as follows:

1Δx2[0.1453Ti−2k+0.4186Ti−1k−1.1279Tik+0.4186Ti+1k+0.1453Ti+2k](21)

Eq. (21) is the approximation of Txx(xi,tk). By removing the time index, the function of T(xi−2), T(xi−1), T(xi+1), T(xi+2) can be expanded by Taylor’s theorem. Take T(xi+1) for example, it can be expanded as follows:

T(xi+1)=T(xi+Δx)=T(xi)+ΔxTx(xi)+(Δx)22Txx(xi)+(Δx)33!Txxx(xi)+(Δx)44!Txxxx(xi,tk)+⋯

After calculation, the truncation error is as follows:

τ=14!⋅5.4868Txxxx|x=ηi(Δx)2≈0.2286Txxxx|x=ηi(Δx)2(22)

where ηi is a point between xi−2 and xi+2.

Similarly, for the case of polynomial order N=2 and integer parameter m=3, the detailed derivation process of the stability analysis of the FT-PDDO scheme is presented in Appendix A. The stability condition is as follows:

λ≤2.8(23)

The truncation error is as follows:

τ≈0.5098Txxxx|x=ηi(Δx)2(24)

The stability conditions of the FT-PDDO scheme in Eqs. (20) and (23) are less strict than that of explicit FEM and explicit FDM, which the former one is λ=DΔt/Δx2≤2/π2 and the latter one is λ=DΔt/Δx2≤1/2, respectively.

3.4.2 Stability and Convergence Analysis of the BT-PDDO Scheme

For the BT-PDDO scheme in Eq. (13), when the polynomial order is taken as N=2 and the integer parameter is taken as m=2, it has the following form:

Tik+1=Tik+DΔt/Δx2[0.1453Ti−2k+1+0.4186Ti−1k+1−1.1279Tik+1+0.4186Ti+1k+1+0.1453Ti+2k+1](25)

and the amplification factor is:

κ=11−λ[−1.1279+0.2907cos⁡(2rΔx)+0.8372cos⁡(rΔx)](26)

The stability requirement is |κ|≤1 and it is not hard to show that this always holds for the amplification factor in Eq. (26). Therefore, in this case, the BT-PDDO scheme is unconditionally stable.

For the case of polynomial order N=2 and integer parameter m=3, the detailed derivation process of the stability analysis of the BT-PDDO scheme is presented in Appendix B. In this case, the BT-PDDO scheme is unconditionally stable.

The truncation error for the BT-PDDO scheme is the same as that of the FT-PDDO scheme, which is Eqs. (22) and (24) in the case of polynomial order N=2 and integer parameter m=2 and m=3, respectively.

3.4.3 Stability Analysis of the CT-PDDO Scheme

For the CT-PDDO scheme in Eq. (14), when the polynomial order is taken as N=2 and the integer parameter is taken as m=2, it has the following form:

Tik+1 =Tik+λ2[0.1453Ti−2k+1+0.4186Ti−1k+1−1.1279Tik+1+0.4186Ti+1k+1+0.1453Ti+2k+1] +λ2[0.1453Ti−2k+0.4186Ti−1k−1.1279Tik+0.4186Ti+1k+0.1453Ti+2k](27)

The amplification factor is as follows:

κ=1+λ2[−1.1279+0.2907cos⁡(2rΔx)+0.8372cos⁡(rΔx)]1−λ2[−1.1279+0.2907cos⁡(2rΔx)+0.8372cos⁡(rΔx)](28)

The stability requirement is |κ|≤1 and it is not hard to show that this always holds for the amplification factor in Eq. (28). Therefore, in this case, the CT-PDDO scheme is unconditionally stable.

For the case of polynomial order N=2 and integer parameter m=3, the detailed derivation process of the stability analysis of the CT-PDDO scheme is presented in Appendix C. In this case, the CT-PDDO scheme is unconditionally stable.

The truncation error for the CT-PDDO scheme is the same as that of the FT-PDDO scheme, which is Eqs. (22) and (24) in the case of polynomial order N=2 and integer parameter m=2 and m=3, respectively.

3.5 FD-PDDO Method for Solving Two-Dimensional Transient Heat Conduction Equations

This section focuses on the schemes for solving two-dimensional transient heat conduction equations. The FT-PDDO scheme is obtained using the forward difference formula in time and the PDDO method in space as follows:

Tik+1=Tik+DΔt[∑j=1N(i)Tk(x1(j))gN20(ξ1(j,i),ξ2(j,i))ΔAj+∑j=1N(i)Tk(x2(j))gN02(ξ1(j,i),ξ2(j,i))ΔAj](29)

where ξ1(j,i)=x1(j)−x1(i), ξ2(j,i)=x2(j)−x2(i), ΔAj=Δx1Δx2.

The BT-PDDO scheme is obtained using a backward difference formula in time and the PDDO method in space, which is:

Tik+1=Tik+DΔt[∑j=1N(i)Tk+1(x1(j))gN20(ξ1(j,i),ξ2(j,i))ΔAj+∑j=1N(i)Tk+1(x2(j))gN02(ξ1(j,i),ξ2(j,i))ΔAj](30)

The CT-PDDO scheme is obtained by using the centered difference formula in time and the PDDO method in space, namely

Tik+1 =Tik+DΔt2[∑j=1N(i)Tk+1(x1(j))gN20(ξ1(j,i),ξ2(j,i))ΔAj+∑j=1N(i)Tk+1(x2(j))gN02(ξ1(j,i),ξ2(j,i))ΔAj] +DΔt2[∑j=1N(i)Tk(x1(j))gN20(ξ1(j,i),ξ2(j,i))ΔAj+∑j=1N(i)Tk(x2(j))gN02(ξ1(j,i),ξ2(j,i))ΔAj](31)

3.6 Stability and Convergence Analysis of the FD-PDDO Method for Solving Two-Dimensional Transient Heat Conduction Equations

The weight function is chosen as ω(ξ1,ξ2)=e−(2ξ1+2ξ2)2/δ2, where the interaction domain is δ=mΔx with m the integer parameter and Δx the spatial step size. Herein, the polynomial order N=2, integer parameter m=2 and m=3 are considered. Fig. 4 shows the PD interaction domains for the discretized Point X(i) with the interaction domain of δ=mΔx (m=2).

Figure 4: PD interaction domains for the discretized point X(i)

3.6.1 Stability and Convergence Analysis of the FT-PDDO Method

For the FT-PDDO scheme in Eq. (29), in the case of polynomial order N=2 and integer parameter m=2, the right-hand terms ∑j=1N(i)Tk(x1(j))gN20(ξ1(j,i),ξ2(j,i))ΔAj and ∑j=1N(i)Tk(x2(j))gN02(ξ1(j,i),ξ2(j,i))ΔAj can be expanded as follows:

 ∑j=1N(i)Tk(x1(j))gN20(ξ1(j,i),ξ2(j,i))ΔAj=Δx1Δx2∑q=−22Ti−2,j+qkg220(−2Δx1,qΔx2) +Δx1Δx2∑q=−22Ti−1,j+qkg220(−Δx1,qΔx2)+Δx1Δx2∑q=−22Ti,j+qkg220(0,qΔx2) +Δx1Δx2∑q=−22Ti+1,j+qkg220(Δx1,qΔx2)+Δx1Δx2∑q=−22Ti+2,j+qkg220(2Δx1,qΔx2)

and

 ∑j=1N(i)Tk(x2(j))gN02(ξ1(j,i),ξ2(j,i))ΔAj=Δx1Δx2∑q=−22Ti+q,j−2kg202(qΔx1,−2Δx2) +Δx1Δx2∑q=−22Ti+q,j−1kg202(qΔx1,−Δx2)+Δx1Δx2∑q=−22Ti+q,jkg202(qΔx1,0) +Δx1Δx2∑q=−22Ti+q,j+1kg202(qΔx1,Δx2)+Δx1Δx2∑q=−22Ti+q,j+2kg202(qΔx1,2Δx2)

where

g220(ξ1,ξ2) =a0020ω(ξ1,ξ2)+a1020ω(ξ1,ξ2)ξ1+a0120ω(ξ1,ξ2)ξ2 +a2020ω(ξ1,ξ2)ξ12+a0220ω(ξ1,ξ2)ξ22+a1120ω(ξ1,ξ2)ξ1ξ2

g202(ξ1,ξ2) =a0002ω(ξ1,ξ2)+a1002ω(ξ1,ξ2)ξ1+a0102ω(ξ1,ξ2)ξ2 +a2002ω(ξ1,ξ2)ξ12+a0202ω(ξ1,ξ2)ξ22+a1102ω(ξ1,ξ2)ξ1ξ2

The unknown coefficient can be obtained from the following equation:

Aa=b,

with A=∫Hxω[1ξ1ξ2ξ12ξ22ξ1ξ2ξ1ξ12ξ1ξ2ξ13ξ1ξ22ξ12ξ2ξ2ξ1ξ2ξ22ξ12ξ2ξ23ξ1ξ22ξ12ξ13ξ12ξ2ξ14ξ12ξ22ξ13ξ2ξ22ξ1ξ22ξ23ξ12ξ22ξ24ξ1ξ23ξ1ξ2ξ12ξ2ξ1ξ22ξ13ξ2ξ1ξ23ξ12ξ22]dV ′, a=[a0020a0002a1020a1002a0120a0102a2020a2002a0220a0202a1120a1102], and b=[000000200200].

The focus herein is solely on a uniform grid with Δx=Δx1=Δx2. Hence, the non-zero elements in matrix A can be obtained as A0000=3.1414Δx2, A0020=A0002=A1010=A0101=A2000=A0200=1.5638Δx4, A2020=A0202=2.3429Δx6, and A2002=A0220=A1111=0.7784Δx6. Thus, the coefficients for the PD function are as follows:

{a0020=−0.63641Δx4a1020=0a0120=0a2020=1.27841Δx6a0220=0a1120=0 and {a0002=−0.63641Δx4a1002=0a0102=0a2002=0a0202=1.27841Δx6a1102=0.

The FT-PDDO scheme is obtained as shown in Eq. (32).

Ti,jk+1 =Ti,jk+λ[0.003Ti−2,j−2k+0.0345Ti−1,j−2k+0.0703Ti,j−2k+0.0345Ti+1,j−2k+0.003Ti+2,j−2k +0.0345Ti−2,j−1k+0.1738Ti−1,j−1k+0.0021Ti,j−1k+0.1738Ti+1,j−1k+0.0345Ti+2,j−1k +0.0703Ti−2,jk+0.0021Ti−1,jk−1.2728Ti,jk+0.0021Ti+1,jk+0.0703Ti+2,jk +0.0345Ti−2,j+1k+0.1738Ti−1,j+1k+0.0021Ti,j+1k+0.1738Ti+1,j+1k+0.0345Ti+2,j+1k +0.003Ti−2,j+2k+0.0345Ti−1,j+2k+0.0703Ti,j+2k+0.0345Ti+1,j+2k+0.003Ti+2,j+2k](32)

where λ=DΔt/Δx2.

The Fourier method is now utilized to show the stability [1]. Let Ti,jk=ωker1x1iIer2x2jI (where I=−1 and r is the positive constant). Substituting it into Eq. (32) yields the amplification factor as shown in Eq. (33).

κ=ωk+1ωk =1−1.2728λ+λ[0.012cos⁡(2r1Δx)cos⁡(2r2Δx)+0.138cos⁡(r1Δx)cos⁡(2r2Δx) +0.138cos⁡(2r1Δx)cos⁡(r2Δx)+0.6952cos⁡(r1Δx)cos⁡(r2Δx) +0.1406cos⁡(2r2Δx)+0.1406cos⁡(2r1Δx)+0.0042cos⁡(r2Δx)+0.0042cos⁡(r1Δx)](33)

The stability requirement is |κ|≤1. Thus, the stability condition is shown in Eq. (34).

λ≤1.1(34)

The Taylor's theorem is used to show the truncation error. It is shown as follows:

τ=(Δx)2[0.228625Tx1x1x1x1+0.228625Tx2x2x2x2+0.4978Tx1x1x2x2]|x1=ηi,x2=ξj(35)

where (ηi,ξj) is the point in the area enclosed by (x1,i−2,x2,j−2), (x1,i+2,x2,j−2), (x1,i−2,x2,j+2) and (x1,i+2,x2,j+2).

Similarly, for the case of polynomial order N=2 and integer parameter m=3, the detailed derivation process of the stability analysis of the FT-PDDO scheme is presented in Appendix D. The stability condition is as follows:

λ≤2.8(36)

The truncation error is:

τ=(Δx)2[0.509975Tx1x1x1x1+0.509975Tx2x2x2x2+1.1165Tx1x1x2x2]|x1=ηi,x2=ξj(37)