images
Computer Modeling in Engineering & Sciences
images

DOI: 10.32604/cmes.2024.050003

ARTICLE

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

Chunlei Ruan1,2,*, Cengceng Dong1, Zeyue Zhang1, Boyu Chen1 and Zhijun Liu1

1School of Mathematics & Statistics, Henan University of Science & Technology, Luoyang, 471023, China
2Longmen Laboratory, Luoyang, 471023, China
*Corresponding Author: Chunlei Ruan. Email: ruanchunlei@haust.edu.cn
Received: 24 January 2024; Accepted: 03 April 2024; Published: 08 July 2024

Abstract: Transient heat conduction problems widely exist in engineering. In previous work on the peridynamic differential operator (PDDO) method for solving such problems, both time and spatial derivatives were discretized using the PDDO method, resulting in increased complexity and programming difficulty. In this work, the forward difference formula, the backward difference formula, and the centered difference formula are used to discretize the time derivative, while the PDDO method is used to discretize the spatial derivative. Three new schemes for solving transient heat conduction equations have been developed, namely, the forward-in-time and PDDO in space (FT-PDDO) scheme, the backward-in-time and PDDO in space (BT-PDDO) scheme, and the central-in-time and PDDO in space (CT-PDDO) scheme. The stability and convergence of these schemes are analyzed using the Fourier method and Taylor’s theorem. Results show that the FT-PDDO scheme is conditionally stable, whereas the BT-PDDO and CT-PDDO schemes are unconditionally stable. The stability conditions for the FT-PDDO scheme are less stringent than those of the explicit finite element method and explicit finite difference method. The convergence rate in space for these three methods is two. These constructed schemes are applied to solve one-dimensional and two-dimensional transient heat conduction problems. The accuracy and validity of the schemes are verified by comparison with analytical solutions.

Keywords: Peridynamic differential operator; finite difference method; stability; transient heat conduction problem

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)

The stability conditions of the FT-PDDO scheme in two-dimensional shown in Eqs. (34) and (36) are equal to the one-dimensional shown in Eqs. (20) and (23), respectively. This means that the stability conditions of the FT-PDDO scheme are independent of dimensionality. The values of stability conditions of the FT-PDDO scheme are less strict than those of explicit finite element method of λ=DΔt/Δx2≤2/π2 and explicit finite difference method of λ=DΔt/Δx2≤1/4.

3.6.2 Stability and Convergence Analysis of the BT-PDDO Method

For the BT-PDDO scheme in Eq. (30), in the case of polynomial order N=2 and integer parameter m=2, the detailed form is:

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

Furthermore, the amplification factor is:

κ =1/[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)]](39)

The stability requirement is |κ|≤1 and it can be seen that this always holds for the amplification factor in Eq. (39). 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 E. 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. (35) and (37) in the case of polynomial order N=2 and integer parameter m=2 and m=3, respectively.

3.6.3 Stability and Convergence Analysis of the CT-PDDO Method

For the CT-PDDO scheme in Eq. (31), in the case of polynomial order N=2 and integer parameter m=2, the detailed form is:

Tik+1 =Tik+λ2[0.003Ti−2,j−2k+1+0.0345Ti−1,j−2k+1+0.0703Ti,j−2k+1+0.0345Ti+1,j−2k+1+0.003Ti+2,j−2k+1 +0.0345Ti−2,j−1k+1+0.1738Ti−1,j−1k+1+0.0021Ti,j−1k+1+0.1738Ti+1,j−1k+1+0.0345Ti+2,j−1k+1 +0.0703Ti−2,jk+1+0.0021Ti−1,jk+1−1.2728Ti,jk+1+0.0021Ti+1,jk+1+0.0703Ti+2,jk+1 +0.0345Ti−2,j+1k+1+0.1738Ti−1,j+1k+1+0.0021Ti,j+1k+1+0.1738Ti+1,j+1k+1+0.0345Ti+2,j+1k+1 +0.003Ti−2,j+2k+1+0.0345Ti−1,j+2k+1+0.0703Ti,j+2k+1+0.0345Ti+1,j+2k+1+0.003Ti+2,j+2k+1] +λ2[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](40)

Moreover, the amplification factor is as follows:

κ ={1+λ2[λ2−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)λ2]}/{1−λ2[λ2−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)λ2]}(41)

The stability requirement is |κ|≤1 and it can be observed that this always holds for the amplification factor in Eq. (41). 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 F. 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. (35) and (37) in the case of polynomial order N=2 and integer parameter m=2 and m=3, respectively.

4  Numerical Examples

The global error ε is defined as follows:

ε=1|Te|max1K∑i=1K(Tie−Ti)2(42)

where |Te|max denotes the maximum of the absolute value of the temperature in the analytical solution; Tie and Ti are the analytical solution and numerical solution at the node xi, respectively. The parameter, K, represents the total number of PD points in the discretization.

4.1 One-Dimensional Transient Heat Conduction Problem

The first example is the one-dimensional transient heat conduction problem with Dirichlet boundary condition. The equation, initial condition, and boundary condition are shown as follows:

{∂T∂t=∂2T∂x20<x<π,0<t≤10T(x,0)=sin⁡xT(0,t)=0T(π,t)=0(43)

The analytical solution for this problem is T(x,t)=e−tsin⁡(x). In the FD-PDDO scheme, the weight function is taken as ω(ξ)=e−(2ξ/δ)2 and the interaction domain is chosen as δ=mΔx (m=2 or m=3); the initial condition can be discretized as: Ti0=sin⁡(xi),i=1,2,⋯,K; the boundary condition can be discretized as: ∑j=1N(i)Tk(xj)g20(xj−xi)Δx=0, for xi=0+Δx/2 and xi=π−Δx/2.

Tables 1 and 2 show the comparison of global error and the rate of convergence when taking interaction domain constant m=2 and m=3 in FD-PDDO schemes with polynomial order N=2, respectively. It is evident that all the schemes have a convergent rate of around 2, which is consistent with the theoretical analysis. The convergent solution can be obtained when the time step size satisfies the stability condition in the FT-PDDO scheme. In addition, with the increase in the interaction domain constant m, the stability region of the FT-PDDO scheme is enlarged, and the convergent solution can be obtained even if the time step size is large. The maximum time step size of FT-PDDO has increased five times and two times, respectively, with interaction domain constant m=2 compared with explicit FEM (λ=DΔt/Δx2≤2/π2) and explicit FDM (λ=DΔt/Δx2≤1/2), respectively. Meanwhile, the maximum time step size can be magnified thirteen times and five times in FT-PDDO with interaction domain constant m=3 compared with explicit FEM (λ=DΔt/Δx2≤2/π2) and explicit FDM (λ=DΔt/Δx2≤1/2), respectively.

4.2 Two-Dimensional Transient Heat Conduction Problem

The second example is the two-dimensional transient heat conduction problem with Dirichlet boundary condition. The equation, initial condition, and boundary condition are shown as follows:

{∂T∂t=D(∂2T∂x12+∂2T∂x22)0<x1<Lx1,0<x2<Lx2,0<t≤tendT(x1,x2,0)=30T(0,x2,t)=T(Lx1,x2,t)=T(x1,0,t)=T(x1,Lx2,t)=0(44)

The analytical solution is [3] T(x1,x2,t)=∑i=1∞∑j=1∞Aijsin⁡iπx1Lx1sin⁡jπx2Lx2exp⁡[−(kpi2π2Lx12+kpj2π2Lx22)t], where Aij=120ijπ2[(−1)i−1][(−1)j−1]. Herein, the parameters are taken as: Lx1=Lx2=3m, tend=1.2s, kp=1.25W/(m⋅K), ρc=1J/(m3⋅K) which means D=kpρc=1.25m2/s. In the FD-PDDO scheme, the weighting function is taken as ω(ξ1,ξ2)=e−(2ξ1+2ξ2)2/δ2 and the interaction domain is considered as δ=mΔx (m=2 or m=3). The boundary condition can be discretized as follows: ∑j=1N(i)Tk(x1(j))g200(ξ1(j,i),ξ2(j,i))ΔAj=0 for x1(i)=0+Δx/2, x1(i)=Lx1−Δx/2, and ∑j=1N(i)Tk(x2(j))g200(ξ1(j,i),ξ2(j,i))ΔAj=0 for x2(i)=0+Δx/2, x2(i)=Lx2−Δx/2.

Tables 3 and 4 show the comparison of global error and the rate of convergence when taking interaction domain constant m=2 and m=3 in FD-PDDO schemes with polynomial order N=2, respectively. It can be seen that all the schemes have convergent rates around 2 which is consistent with the theoretical analysis. The FT-PDDO scheme gives the convergent solution even if the time step is large. The maximum time step size can be magnified five times and four times in FT-PDDO with interaction domain constant m=2 compared with explicit FEM (λ=DΔt/Δx2≤2/π2) and explicit FDM (λ=DΔt/Δx2≤1/4), respectively. Meanwhile, the maximum time step size can be magnified ten times and eight times in FT-PDDO with interaction domain constant m=3 compared with explicit FEM (λ=DΔt/Δx2≤2/π2) and explicit FDM (λ=DΔt/Δx2≤1/4), respectively.

4.3 Two-Dimensional Transient Heat Conduction Problem with Both Dirichlet and Neumann Boundary Condition

The third example is the two-dimensional transient heat conduction problem with both Dirichlet and Neumann boundary conditions. The equation, initial condition, and boundary condition are shown as follows:

{∂T∂t=D(∂2T∂x12+∂2T∂x22)0<x1<Lx1,0<x2<Lx2,0<t≤tendT(x1,x2,0)=30Tx(0,x2,t)=0,T(Lx1,x2,t)=T(x2,0,t)=T(x1,Lx2,t)=0(45)

The analytical solution is [3] T(x1,x2,t)=∑i=1∞∑j=1∞Bijcos⁡(2i−1)πx12Lx1sin⁡jπx2Lx2exp⁡[−(k(2i−1)2π24Lx12+kj2π2Lx22)t], where Bij=240(2i−1)jπ2(−1)i+2[(−1)j−1]. Herein, the parameters are taken as follows: Lx1=Lx2=3m, tend=1.2s, kp=1.25W/(m⋅K), ρc=1J/(m3⋅K) which means D=kpρc=1.25m2/s. In the FD-PDDO scheme, the boundary condition can be discretized as: ∑j=1N(i)Tk(x1(j))g210(ξ1(j,i),ξ2(j,i))ΔAj=0 for x1(i)=0+Δx/2, ∑j=1N(i)Tk(x1(j))g200(ξ1(j,i),ξ2(j,i))ΔAj=0 for x1(i)=Lx1−Δx/2, and ∑j=1N(i)Tk(x2(j))g200(ξ1(j,i),ξ2(j,i))ΔAj=0 for x2(i)=0+Δx/2, x2(i)=Lx2−Δx/2.

Tables 5 and 6 show the comparison of global error and the rate of convergence when taking interaction domain constant m=2 and m=3 in FD-PDDO schemes with polynomial order N=2, respectively. As shown in Tables 5 and 6, the rate of convergence of the FT-PDDO scheme is around 2, which is consistent with the theoretical analysis. The FT-PDDO scheme gives the convergent solution even if the time step is large.

5  Conclusions

In this study, the FD-PDDO schemes for solving one-dimensional and two-dimensional transient heat conduction equations are constructed. These schemes utilize the finite difference method to discretize the time derivative and the PDDO method to discretize the spatial derivative. The FD-PDDO schemes, which include the FT-PDDO scheme, the BT-PDDO scheme, and the CT-PDDO scheme, are developed. The stability and convergence of these schemes are analyzed using the Fourier method and Taylor's theorem, respectively. The performance of the schemes in solving transient heat conduction equations is investigated, and the results are compared to those of the analytical solutions. The conclusions are as follows:

(1) The FT-PDDO scheme is conditionally stable, with the stability condition λ=DΔt/Δx2≤1.1 in the case of polynomial order N=2 and interaction domain constant m=2, and λ=DΔt/Δx2≤2.8 in the case of polynomial order N=2 and interaction domain constant m=3 in both one-dimensional and two-dimensional cases. Compared to the explicit FEM and FDM, the FT-PDDO method has less strict stability conditions. In addition, both the BT-PDDO scheme and the CT-PDDO scheme are unconditionally stable.

(2) The FD-PDDO schemes, including the FT-PDDO scheme, the BT-PDDO scheme, and the CT-PDDO scheme, have a convergence rate of 2 in space when the polynomial order N=2.

This study introduces three new schemes, namely the FT-PDDO scheme, the BT-PDDO scheme, and the CT-PDDO scheme, for solving one-dimensional and two-dimensional transient heat conduction equations. Numerical examples demonstrate their effectiveness. The algorithm's approach can also be extended to solve more complex differential equations. Furthermore, given that the PDDO method can handle complex geometries [17,28] and discontinuity problems [31], it is anticipated that our method will find wider practical applications.

Acknowledgement: The authors acknowledge the School of Mathematics and Statistics, Henan University of Science and Technology, for allowing the use of their high-performance computing. The authors acknowledge the reviewer for the valuable revision suggestions.

Funding Statement: This work was financially supported by the Key Science and Technology Project of Longmen Laboratory (No. LMYLKT-001), Innovation and Entrepreneurship Training Program for College Students of Henan Province (No. 202310464050).

Author Contributions: The authors confirm contribution to the paper as follows: study conception and design: Chunlei Ruan; data collection: Cengceng Dong, Zeyue Zhang, Boyu Chen; analysis and interpretation of results: Chunlei Ruan, Zhijun Liu; draft manuscript preparation: Chunlei Ruan. Cengceng Dong. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available upon reasonable request from the corresponding author.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

1. Holmes MH. Introduction to numerical methods in differential equations. New York: Springer; 2007.

2. Özişik MN, Orlande HR, Colaço MJ, Cotta RM. Finite difference methods in heat transfer. Boca Raton: CRC Press; 2017.

3. Bruch Jr JC, Zyvoloski G. Transient two-dimensional heat conduction problems solved by the finite element method. Int J Numer Meth Eng. 1974;8(3):481–94. doi:10.1002/nme.v8:3.

4. Wilson EL, Nickell RE. Application of the finite element method to heat conduction analysis. Nucl Eng Des. 1966;4(3):276–86. doi:10.1016/0029-5493(66)90051-3.

5. Minkowycz W. Advances in numerical heat transfer. Boca Raton: CRC Press; 1996.

6. Chai JC, Lee H, Patankar SV. Finite-volume method for radiation heat transfer. J Thermophys Heat Tran. 1994;8(3):419–25. doi:10.2514/3.559.

7. Brebbia CA, Telles JC, Wrobel LC. Boundary element techniques: theory and applications in engineering. New York: Springer-Verlag; 2012.

8. Sun W, Ma H, Qu W. A hybrid numerical method for non-linear transient heat conduction problems with temperature-dependent thermal conductivity. Appl Math Lett. 2024;148:108868. doi:10.1016/j.aml.2023.108868.

9. Sun W, Qu W, Gu Y, Li PW. An arbitrary order numerical framework for transient heat conduction problems. Int J Heat Mass Tran. 2024;218:124798. doi:10.1016/j.ijheatmasstransfer.2023.124798.

10. Liu GR, Liu MB. Smoothed particle hydrodynamics: a meshfree particle method. Singapore: World Scientific Publishing; 2003.

11. Sladek J, Sladek V, Hellmich C, Eberhardsteiner J. Heat conduction analysis of 3-D axisymmetric and anisotropic FGM bodies by meshless local Petrov-Galerkin method. Comput Mech. 2007;39:323–33.

12. Wu XH, Tao WQ. Meshless method based on the local weak-forms for steady-state heat conduction problems. Int J Heat Mass Tran. 2008;51(11–12):3103–12.

13. Zhang XH, Xiang H. Variational multiscale element free Galerkin method for convection-diffusion-reaction equation with small diffusion. Eng Anal Bound Elem. 2014;46:85–92. doi:10.1016/j.enganabound.2014.05.010.

14. Zhang XH, Zhang P. Meshless modeling of natural convection problems in non-rectangular cavity using the variational multiscale element free Galerkin method. Eng Anal Bound Elem. 2015;61:287–300. doi:10.1016/j.enganabound.2015.08.005.

15. Soleimani S, Jalaal M, Bararnia H, Ghasemi E, Ganji DD, Mohammadi F. Local RBF-DQ method for two-dimensional transient heat conduction problems. Int Commun Heat Mass. 2010;37(9):1411–8. doi:10.1016/j.icheatmasstransfer.2010.06.033.

16. Zhou LM, Zhu ZD, Que XC. Simulation of non-fourier heat conduction in discontinuous heterogeneous materials based on the peridynamic method. Therm Sci. 2023;27:917–31. doi:10.2298/TSCI220803157Z.

17. Wen ZX, Hou C, Zhao MY, Wan XP. A peridynamic model for non-Fourier heat transfer in orthotropic plate with uninsulated cracks. Appl Math Model. 2023;115:706–23. doi:10.1016/j.apm.2022.11.010.

18. Guo YT, Huang XH, Jin YL, Feng MS, Zheng Z, Su GS. Discussion on the form of construction function in the peridynamic differential operator based on relative function. Eng Anal Bound Elem. 2023;151:136–63. doi:10.1016/j.enganabound.2023.02.042.

19. Zhou BL, Li ZY, Xu YP, Huang D. Analysis of nonlinear heat conduction problems with temperature-dependent conductivity using peridynamic differential operator. Int J Appl Mech. 2022;14(5):2250047. doi:10.1142/S1758825122500478.

20. Silling SA, Epton M, Weckner O, Xu J, Askari E. Peridynamic states and constitutive modeling. J Elasticity. 2007;88:151–84. doi:10.1007/s10659-007-9125-1.

21. Jiang L, Zhang XH, Zheng BJ, Peng H, Gao XW. A reduced-order peridynamic differential operator for unsteady convection-diffusion problems. Eng Anal Bound Elem. 2024;161:1–10. doi:10.1016/j.enganabound.2024.01.010.

22. Dorduncu M. Stress analysis of sandwich plates with functionally graded cores using peridynamic differential operator and refined zigzag theory. Thin Wall Struct. 2020;146:106468. doi:10.1016/j.tws.2019.106468.

23. Gao Y, Oterkus S. Nonlocal modeling for fluid flow coupled with heat transfer by using peridynamic differential operator. Eng Anal Bound Elem. 2019;105:104–21. doi:10.1016/j.enganabound.2019.04.007.

24. Li ZY, Huang D, Xu YP, Yan KH. Nonlocal steady-state thermoelastic analysis of functionally graded materials by using peridynamic differential operator. Appl Math Model. 2021;93:294–313. doi:10.1016/j.apm.2020.12.004.

25. Bergel GL, Li SF. The total and updated lagrangian formulations of state-based peridynamics. Comput Mech. 2016;58:351–70. doi:10.1007/s00466-016-1297-8.

26. Kan XY, Yan JL, Li SF, Zhang AM. On differences and comparisons of peridynamic differential operators and nonlocal differential operators. Comput Mech. 2021;68:1349–67. doi:10.1007/s00466-021-02072-8.

27. Yu HC, Li SF. On approximation theory of nonlocal differential operators. Int J Numer Meth Eng. 2021;122(23):6984–7012. doi:10.1002/nme.v122.23.

28. Zhou BL, Li ZY, Huang D. Peridynamic differential operator analysis of two-dimensional transient heat conduction. Appl Math Mech. 2022;43(6):660–8.

29. Madenci E, Barut A, Dorduncu M. Peridynamic differential operator for numerical analysis. Berlin: Springer International Publishing; 2019.

30. Madenci E, Barut A, Futch M. Peridynamic differential operator and its applications. Comput Method Appl Mech Eng. 2016;304:408–51. doi:10.1016/j.cma.2016.02.028.

31. Ruan CL, Dong CC, Liang KF, Liu ZJ, Bao XR. Euler’s first-order explicit method-peridynamic differential operator for solving population balance equations of the crystallization process. Comput Model Eng Sci. 2024;138:3033–49. doi:10.32604/cmes.2023.030607.

Appendix A. Stability analysis of the FT-PDDO method for solving one-dimensional transient heat conduction equations

For the case of polynomial order N=2 and integer parameter m=3, the FT-PDDO scheme is:

Tik+1 =Tik+λ[0.0455Ti−3k+0.1535Ti−2k−0.0233Ti−1k−0.3513Tik−0.0233Ti+1k+0.1535Ti+2k +0.0455Ti+3k]

and the amplification factor is:

κ=ωk+1/ωk=1+λ[−0.3513+0.0909cos⁡(3rΔx)+0.3070cos⁡(2rΔx)−0.0467cos⁡(rΔx)]

The stability requirement is |κ|≤1. Thus, the stability condition can be obtained as follows:

λ≤2.8

Appendix B. Stability analysis of the BT-PDDO method for solving one-dimensional transient heat conduction equations

For the case of polynomial order N=2 and integer parameter m=3, the BT-PDDO scheme is as follows:

Tik+1 =Tik+λ[0.0455Ti−3k+1+0.1535Ti−2k+1−0.0233Ti−1k+1−0.3513Tik+1−0.0233Ti+1k+1+0.1535Ti+2k+1 +0.0455Ti+3k+1]

Additionally, the amplification factor is:

κ=11−λ[−0.3513+0.0909cos⁡(3rΔx)+0.3070cos⁡(2rΔx)−0.0467cos⁡(rΔx)]

The stability requirement always holds and is |κ|≤1. Therefore, in this case, the BT-PDDO scheme is stable.

Appendix C. Stability analysis of the CT-PDDO method for solving one-dimensional transient heat conduction equations

For the case of polynomial order N=2 and integer parameter m=3, the CT-PDDO scheme is:

Tik+1 =Tik+λ2[0.0455Ti−3k+0.1535Ti−2k−0.0233Ti−1k−0.3513Tik−0.0233Ti+1k+0.1535Ti+2k+0.0455Ti+3k] +λ2[0.0455Ti−3k+1+0.1535Ti−2k+1−0.0233Ti−1k+1−0.3513Tik+1−0.0233Ti+1k+1+0.1535Ti+2k+1+0.0455Ti+3k+1]

Moreover, the amplification factor is as follows:

κ=1+λ2[−0.3513+0.0909cos⁡(3rΔx)+0.3070cos⁡(2rΔx)−0.0467cos⁡(rΔx)]1−λ2[−0.3513+0.0909cos⁡(3rΔx)+0.3070cos⁡(2rΔx)−0.0467cos⁡(rΔx)]

The stability requirement is |κ|≤1. It can be seen that it always holds for the amplification factor in the above equation. Therefore, in this case, the CT-PDDO scheme is stable.

Appendix D. Stability analysis of the FT-PDDO method for solving two-dimensional transient heat conduction equations

For the case of polynomial order N=2 and integer parameter m=3, the FT-PDDO scheme is as follows:

Ti,jk+1 =Ti,jk+λ[0.0006Ti−3,j−3k+0.004Ti−2,j−3k+0.0108Ti−1,j−3k+0.0147Ti,j−3k+0.0108Ti+1,j−3k+0.004Ti+2,j−3k+0.0006Ti+3,j−3k +0.004Ti−3,j−2k+0.0195Ti−2,j−2k+0.0356Ti−1,j−2k+0.0354Ti,j−2k+0.0356Ti+1,j−2k+0.0195Ti+2,j−2k+0.004Ti+3,j−2k +0.0108Ti−3,j−1k+0.0356Ti−2,j−1k−0.0113Ti−1,j−1k−0.0936Ti,j−1k−0.0113Ti+1,j−1k+0.0356Ti+2,j−1k+0.0108Ti+3,j−1k +0.0147Ti−3,jk+0.0354Ti−2,jk−0.0936Ti−1,jk−0.2644Ti,jk−0.0936Ti+1,jk+0.0354Ti+2,jk+0.0147Ti+3,jk +0.0108Ti−3,j+1k+0.0356Ti−2,j+1k−0.0113Ti−1,j+1k−0.0936Ti,j+1k−0.0113Ti+1,j+1k+0.0356Ti+2,j+1k+0.0108Ti+3,j+1k +0.004Ti−3,j+2k+0.0195Ti−2,j+2k+0.0356Ti−1,j+2k+0.0354Ti,j+2k+0.0356Ti+1,j+2k+0.0195Ti+2,j+2k+0.004Ti+3,j+2k +0.0006Ti−3,j+3k+0.004Ti−2,j+3k+0.0108Ti−1,j+3k+0.0147Ti,j+3k+0.0108Ti+1,j+3k+0.004Ti+2,j+3k+0.0006Ti+3,j+3k]

Moreover, the amplification factor is:

κ=ωk+1ωk =1+λ[0.0024cos(3r1Δx)cos(3r2Δx)+0.016cos(2r1Δx)cos(3r2Δx)+0.0432cos(r1Δx)cos(3r2Δx) +0.016cos(3r1Δx)cos(2r2Δx)+0.078cos(2r1Δx)cos(2r2Δx)+0.1424cos(r1Δx)cos(2r2Δx) +0.0432cos(3r1Δx)cos(r2Δx)+0.1424cos(2r1Δx)cos(r2Δx)−0.0452cos(r1Δx)cos(r2Δx) +0.0294cos(3r2Δx)+0.0708cos(2r2Δx)−0.1872cos(r2Δx) +0.0294cos(3r1Δx)+0.0708cos(2r1Δx)−0.1872cos(r1Δx)−0.2644]

The stability requirement is |κ|≤1. Thus, the stability condition is given by:

λ≤2.8

Appendix E. Stability analysis of the BT-PDDO method for solving two-dimensional transient heat conduction equations

For the case of polynomial order N=2 and integer parameter m=3, the BT-PDDO scheme is:

Ti,jk+1 =Ti,jk+λ[0.0006Ti−3,j−3k+1+0.004Ti−2,j−3k+1+0.0108Ti−1,j−3k+1+0.0147Ti,j−3k+1+0.0108Ti+1,j−3k+1+0.004Ti+2,j−3k+1+0.0006Ti+3,j−3k+1 +0.004Ti−3,j−2k+1+0.0195Ti−2,j−2k+1+0.0356Ti−1,j−2k+1+0.0354Ti,j−2k+1+0.0356Ti+1,j−2k+1+0.0195Ti+2,j−2k+1+0.004Ti+3,j−2k+1 +0.0108Ti−3,j−1k+1+0.0356Ti−2,j−1k+1−0.0113Ti−1,j−1k+1−0.0936Ti,j−1k+1−0.0113Ti+1,j−1k+1+0.0356Ti+2,j−1k+1+0.0108Ti+3,j−1k+1 +0.0147Ti−3,jk+1+0.0354Ti−2,jk+1−0.0936Ti−1,jk+1−0.2644Ti,jk+1−0.0936Ti+1,jk+1+0.0354Ti+2,jk+1+0.0147Ti+3,jk+1 +0.0108Ti−3,j+1k+1+0.0356Ti−2,j+1k+1−0.0113Ti−1,j+1k+1−0.0936Ti,j+1k+1−0.0113Ti+1,j+1k+1+0.0356Ti+2,j+1k+1+0.0108Ti+3,j+1k+1 +0.004Ti−3,j+2k+1+0.0195Ti−2,j+2k+1+0.0356Ti−1,j+2k+1+0.0354Ti,j+2k+1+0.0356Ti+1,j+2k+1+0.0195Ti+2,j+2k+1+0.004Ti+3,j+2k+1 +0.0006Ti−3,j+3k+1+0.004Ti−2,j+3k+1+0.0108Ti−1,j+3k+1+0.0147Ti,j+3k+1+0.0108Ti+1,j+3k+1+0.004Ti+2,j+3k+1+0.0006Ti+3,j+3k+1]

The amplification factor is as follows:

κ =1/{1−λ[0.0024cos⁡(3r1Δx)cos⁡(3r2Δx)+0.016cos⁡(2r1Δx)cos⁡(3r2Δx)+0.0432cos⁡(r1Δx)cos⁡(3r2Δx) +0.016cos⁡(3r1Δx)cos⁡(2r2Δx)+0.078cos⁡(2r1Δx)cos⁡(2r2Δx)+0.1424cos⁡(r1Δx)cos⁡(2r2Δx) +0.0432cos⁡(3r1Δx)cos⁡(r2Δx)+0.1424cos⁡(2r1Δx)cos⁡(r2Δx)−0.0452cos⁡(r1Δx)cos⁡(r2Δx) +0.0294cos⁡(3r2Δx)+0.0708cos⁡(2r2Δx)−0.1872cos⁡(r2Δx) +0.0294cos⁡(3r1Δx)+0.0708cos⁡(2r1Δx)−0.1872cos⁡(r1Δx)−0.2644]}

The stability requirement is |κ|≤1. It can be observed that it always holds. Therefore, in this case, the BT-PDDO scheme is stable.

Appendix F. Stability analysis of the CT-PDDO method for solving two-dimensional transient heat conduction equations

For the case of polynomial order N=2 and integer parameter m=3, the CT-PDDO scheme is:

Tik+1 =Tik+λ2[0.0006Ti−3,j−3k+0.004Ti−2,j−3k+0.0108Ti−1,j−3k+0.0147Ti,j−3k+0.0108Ti+1,j−3k+0.004Ti+2,j−3k+0.0006Ti+3,j−3k +0.004Ti−3,j−2k+0.0195Ti−2,j−2k+0.0356Ti−1,j−2k+0.0354Ti,j−2k+0.0356Ti+1,j−2k+0.0195Ti+2,j−2k+0.004Ti+3,j−2k +0.0108Ti−3,j−1k+0.0356Ti−2,j−1k−0.0113Ti−1,j−1k−0.0936Ti,j−1k−0.0113Ti+1,j−1k+0.0356Ti+2,j−1k+0.0108Ti+3,j−1k +0.0147Ti−3,jk+0.0354Ti−2,jk−0.0936Ti−1,jk−0.2644Ti,jk−0.0936Ti+1,jk+0.0354Ti+2,jk+0.0147Ti+3,jk +0.0108Ti−3,j+1k+0.0356Ti−2,j+1k−0.0113Ti−1,j+1k−0.0936Ti,j+1k−0.0113Ti+1,j+1k+0.0356Ti+2,j+1k+0.0108Ti+3,j+1k +0.004Ti−3,j+2k+0.0195Ti−2,j+2k+0.0356Ti−1,j+2k+0.0354Ti,j+2k+0.0356Ti+1,j+2k+0.0195Ti+2,j+2k+0.004Ti+3,j+2k +0.0006Ti−3,j+3k+0.004Ti−2,j+3k+0.0108Ti−1,j+3k+0.0147Ti,j+3k+0.0108Ti+1,j+3k+0.004Ti+2,j+3k+0.0006Ti+3,j+3k] +λ2[0.0006Ti−3,j−3k+1+0.004Ti−2,j−3k+1+0.0108Ti−1,j−3k+1+0.0147Ti,j−3k+1+0.0108Ti+1,j−3k+1+0.004Ti+2,j−3k+1+0.0006Ti+3,j−3k+1 +0.004Ti−3,j−2k+1+0.0195Ti−2,j−2k+1+0.0356Ti−1,j−2k+1+0.0354Ti,j−2k+1+0.0356Ti+1,j−2k+1+0.0195Ti+2,j−2k+1+0.004Ti+3,j−2k+1 +0.0108Ti−3,j−1k+1+0.0356Ti−2,j−1k+1−0.0113Ti−1,j−1k+1−0.0936Ti,j−1k+1−0.0113Ti+1,j−1k+1+0.0356Ti+2,j−1k+1+0.0108Ti+3,j−1k+1 +0.0147Ti−3,jk+1+0.0354Ti−2,jk+1−0.0936Ti−1,jk+1−0.2644Ti,jk+1−0.0936Ti+1,jk+1+0.0354Ti+2,jk+1+0.0147Ti+3,jk+1 +0.0108Ti−3,j+1k+1+0.0356Ti−2,j+1k+1−0.0113Ti−1,j+1k+1−0.0936Ti,j+1k+1−0.0113Ti+1,j+1k+1+0.0356Ti+2,j+1k+1+0.0108Ti+3,j+1k+1 +0.004Ti−3,j+2k+1+0.0195Ti−2,j+2k+1+0.0356Ti−1,j+2k+1+0.0354Ti,j+2k+1+0.0356Ti+1,j+2k+1+0.0195Ti+2,j+2k+1+0.004Ti+3,j+2k+1 +0.0006Ti−3,j+3k+1+0.004Ti−2,j+3k+1+0.0108Ti−1,j+3k+1+0.0147Ti,j+3k+1+0.0108Ti+1,j+3k+1+0.004Ti+2,j+3k+1+0.0006Ti+3,j+3k+1]

Moreover, the amplification factor is as follows:

κ ={1+λ2[0.0024λ2cos⁡(3r1Δx)cos⁡(3r2Δx)+0.016cos⁡(2r1Δx)cos⁡(3r2Δx)+0.0432cos⁡(r1Δx)cos⁡(3r2Δx)+0.016cos⁡(3r1Δx)cos⁡(2r2Δx)+0.078cos⁡(2r1Δx)cos⁡(2r2Δx)+0.1424cos⁡(r1Δx)cos⁡(2r2Δx)+0.0432cos⁡(3r1Δx)cos⁡(r2Δx)+0.1424cos⁡(2r1Δx)cos⁡(r2Δx)+(−0.0452)cos⁡(r1Δx)cos⁡(r2Δx)+0.0294cos⁡(3r2Δx)+0.0708cos⁡(2r2Δx)+(−0.1872)cos⁡(r2Δx)+0.0294cos⁡(3r1Δx)+0.0708cos⁡(2r1Δx)+(−0.1872)cos⁡(r1Δx)+(−0.2644)λ2]}/{1−λ2[0.0024λ2cos⁡(3r1Δx)cos⁡(3r2Δx)+0.016cos⁡(2r1Δx)cos⁡(3r2Δx)+0.0432cos⁡(r1Δx)cos⁡(3r2Δx)+0.016cos⁡(3r1Δx)cos⁡(2r2Δx)+0.078cos⁡(2r1Δx)cos⁡(2r2Δx)+0.1424cos⁡(r1Δx)cos⁡(2r2Δx)+0.0432cos⁡(3r1Δx)cos⁡(r2Δx)+0.1424cos⁡(2r1Δx)cos⁡(r2Δx)+(−0.0452)cos⁡(r1Δx)cos⁡(r2Δx)+0.0294cos⁡(3r2Δx)+0.0708cos⁡(2r2Δx)+(−0.1872)cos⁡(r2Δx)+0.0294cos⁡(3r1Δx)+0.0708cos⁡(2r1Δx)+(−0.1872)cos⁡(r1Δx)+(−0.2644)λ2]}

The stability requirement is |κ|≤1. It can be noticed that it always holds for the above equation. Therefore, in this case, the CT-PDDO scheme is stable.

images
This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.