周保良， 李志遠(yuǎn)， 黃 丹

（河海大學(xué) 力學(xué)與材料學(xué)院，南京 211100）

引 言

瞬態(tài)熱傳導(dǎo)問題是航空航天、機(jī)械化工、生物工程等領(lǐng)域關(guān)注的重點(diǎn)[1].近年來，相關(guān)領(lǐng)域的研究者們分別采用有限單元法(finite element method, FEM)[2-3]、局部徑向基函數(shù)微分求積法 (radial basis function-based differential quadrature, RBF-DQ)[4]、邊界單元法(boundary element method, BEM)[5]、邊界節(jié)點(diǎn)法(boundary knot method, BKM)[6]、雙重互易邊界元法(double reciprocal boundary element method, DRBEM)[7]、Trefftz 有限元法(Trefftz finite element method, Trefftz FEM)[8]等各種數(shù)值方法求解了材料和結(jié)構(gòu)的瞬態(tài)熱傳導(dǎo)問題.在各種數(shù)值方法中，近場動(dòng)力學(xué)(peridynamics, PD)方法[9-10]用積分方程代替經(jīng)典連續(xù)介質(zhì)力學(xué)中的偏微分方程，避免了基于連續(xù)性假設(shè)建模和求解微分方程時(shí)的傳統(tǒng)方法在不連續(xù)問題中的奇異性.最新的近場動(dòng)力學(xué)微分算子(PDDO)[11-12]是基于PD 非局部思想和PD 函數(shù)的正交性提出的一種非局部算子，通過Taylor 級(jí)數(shù)展開和正交函數(shù)的性質(zhì)建立了局部偏導(dǎo)數(shù)和非局部積分之間的聯(lián)系.PDDO 的優(yōu)勢不僅體現(xiàn)在可以求解復(fù)雜的微分方程，而且能夠在不連續(xù)點(diǎn)或奇異點(diǎn)存在時(shí)計(jì)算光滑或分散數(shù)據(jù)的導(dǎo)數(shù)值.此外，PDDO 不限制對空間和時(shí)間變量求偏導(dǎo)的順序，當(dāng)考慮時(shí)間非局部性或者更廣義的時(shí)空非局部性時(shí)，PDDO 將會(huì)有廣闊的適用性.目前該方法已成功應(yīng)用于液體晃蕩[13]、復(fù)合材料[14-15]、熱流耦合[16-17]、熱力耦合[18]等問題分析.

本文針對二維瞬態(tài)熱傳導(dǎo)問題，基于PDDO 理論進(jìn)行求解.首先將熱傳導(dǎo)方程和邊界條件進(jìn)行Taylor 展開，通過構(gòu)建PD 函數(shù)，利用正交性將微分方程重構(gòu)為積分形式，引入Lagrange 乘子，進(jìn)行變分分析，建立了二維瞬態(tài)熱傳導(dǎo)問題的非局部分析模型，然后通過三個(gè)算例驗(yàn)證了PDDO 計(jì)算二維瞬態(tài)熱傳導(dǎo)問題具有收斂速度快，計(jì)算精度高，適用性好等特點(diǎn).

1 基本理論和模型

1.1 近場動(dòng)力學(xué)微分算子理論

如圖1 所示，在一個(gè)三維空間D中，物質(zhì)點(diǎn)x=(x1,x2,x3)在近場范圍Hx內(nèi)與另一個(gè)物質(zhì)點(diǎn)x′=(x1′,x2′,x3′)相互作用，ξ=x′?x是這兩個(gè)點(diǎn)之間的初始相對位置，δ 表示近場范圍Hx的大小.

圖1 物質(zhì)點(diǎn)在任意形狀近場范圍內(nèi)的相互作用Fig. 1 Interaction of material points within an arbitrary-shape near field

事實(shí)上，PDDO 可考慮在M維空間中對標(biāo)量場函數(shù)f(x′)=f(x+ξ)在點(diǎn)x處進(jìn)行N階Taylor 級(jí)數(shù)展開[11]：

其中R(N,x)是展開余項(xiàng)，ni表示對xi微分的階數(shù).

其中δnipi表示Kronecker 符號(hào).

將PD 函數(shù)構(gòu)造成多項(xiàng)式形式：

由此，可以將式(5)用矩陣的形式表達(dá)為

由式(8)可得待求系數(shù)矩陣a，代入式(4)可得PD 函數(shù)，進(jìn)一步通過式(2)確定局部微分的PD 形式.

1.2 二維瞬態(tài)熱傳導(dǎo)的PDDO 模型

在直角坐標(biāo)系下無熱源的二維瞬態(tài)熱傳導(dǎo)方程由下式確定[19]：

其中T(x,y,t)表示溫度，α表示導(dǎo)溫系數(shù).定解條件如下：

初始條件

溫度邊界條件

絕熱邊界條件

其中L和C分別是二維瞬態(tài)熱傳導(dǎo)方程和定解條件中含PD 函數(shù)的系數(shù)矩陣，向量d是定解條件中離散點(diǎn)上的已知值，向量T表示離散各點(diǎn)處的待求溫度.

在式(18)中引入Lagrange 乘子λ，并利用變分原理的概念[11]，可以得到

將式(19)中第二項(xiàng)展開，得到

2 數(shù) 值 算 例

2.1 誤差與收斂性分析

考慮一塊尺寸為3 m × 3 m 的方板，邊界條件如圖2 所示.其中導(dǎo)溫系數(shù)α=1.25 m2/h，初始溫度為T0=30 ℃，取PD 物質(zhì)點(diǎn)之間離散間距Δx=Δy=0.075 m，Δt=0.075 h.

圖2 方形板及邊界條件Fig. 2 A square plate and its boundary conditions

該問題的解析解[6]為

其中

在該方形板中，取表1 所示的5 個(gè)測點(diǎn)，計(jì)算當(dāng)t=1.2 h 時(shí)測點(diǎn)處溫度，并與解析解和不同數(shù)值方法[5-8]比較.從表中數(shù)據(jù)可見PDDO 與解析解和其他數(shù)值解基本一致，且具有較高精度.

表1 t=1.2 h 時(shí)不同方法的數(shù)值結(jié)果比較Table 1 Comparison of numerical results between different methods at t =1.2 h

圖3 為全局誤差測量方法[20]下，PD 物質(zhì)點(diǎn)間距分別取0. 075, 0.1, 0.15 時(shí)的誤差測量值，誤差ε 計(jì)算公式為

其中，|T(e)|max為解析解求得溫度絕對值的最大值，上標(biāo)e 和c 分別表示精確解和本文模型解，參數(shù)K表示離散后PD 點(diǎn)的總數(shù).圖3 折線中每段直線斜率的平均值表示收斂速度，其值為4.82.從圖中可以看出，隨著物質(zhì)點(diǎn)間距的減小，全局誤差逐漸減小，數(shù)值解逐漸收斂于精確解，說明本文模型具有良好的收斂性.

圖3 瞬態(tài)溫度場的PDDO 誤差測量值Fig. 3 The error measure of the PDDO solution for the transient temperature field

2.2 不規(guī)則邊界板的二維瞬態(tài)熱傳導(dǎo)分析

A、B、C 三種形狀板的邊界條件如圖4 中標(biāo)注所示，其中導(dǎo)溫系數(shù)α=17.95 m2/s，初始溫度為T0=0 ℃,取PD 物質(zhì)點(diǎn)之間離散間距Δx=Δy=0.025 m，Δt=0.025 s.

圖4 A、B、C 板的幾何形狀和邊界條件Fig. 4 Geometric shapes and boundary conditions of plates A, B and C

圖5 為t=0.8 s 時(shí)，A、B、C 三種形狀的板在不同數(shù)值方法下[21]，沿指定路徑的溫度分布曲線.從圖中可以看出PDDO 方法與FEM 和RBF-DQ 方法得出的結(jié)果一致.到達(dá)穩(wěn)態(tài)后，圖6～圖8 分別給出運(yùn)用本文模型與上述方法求解板內(nèi)溫度的對比情況.結(jié)果表明，通過PDDO 計(jì)算二維瞬態(tài)熱傳導(dǎo)問題是一種新穎且有效的方法，即使對于復(fù)雜的不規(guī)則邊界問題也具有較高的精度.

圖5 t=0.8 s 時(shí)，A、B、C 三種形狀的板在不同方法下沿指定路徑的溫度分布Fig. 5 Temperature distributions of plates A, B and C along a specified path for different methods at t=0.8 s

圖6 穩(wěn)態(tài)時(shí)，A 板在不同數(shù)值方法下的溫度分布：（a） FEM[21]；（b） RBF-DQ[21]；（c） PDDOFig. 6 Temperature distributions of plate A for different numerical methods in a steady state: （a） FEM[21]; （b） RBF-DQ[21]; （c） PDDO

圖7 穩(wěn)態(tài)時(shí)，B 板在不同數(shù)值方法下的溫度分布：（a） FEM[21]；（b） RBF-DQ[21]；（c） PDDOFig. 7 Temperature distributions of plate B for different numerical methods in a steady state: （a） FEM[21]; （b） RBF-DQ[21]; （c） PDDO

圖8 穩(wěn)態(tài)時(shí)，C 板在不同數(shù)值方法下的溫度分布：（a） FEM[21]；（b） RBF-DQ[21]；（c） PDDOFig. 8 Temperature distributions of plate C for different numerical methods in a steady state: （a） FEM[21]; （b） RBF-DQ[21]; （c） PDDO

注 為了解釋圖中的顏色，讀者可以參考本文的電子網(wǎng)頁版本，后同.

2.3 含裂紋板和含圓孔板的二維瞬態(tài)熱傳導(dǎo)分析

圖9 為含絕熱裂紋和絕熱圓孔的2 m × 2 m 方板幾何模型示意圖.裂紋尺寸為1 m × 0.05 m，圓孔半徑為0.1 m.板的導(dǎo)溫系數(shù)為α=0.1157 m2/d，初始溫度為T0=0 ℃，Γ1和Γ3分別表示T=100 ℃與T= ?100 ℃的溫度邊界條件，Γ2和Γ4表示絕熱邊界條件，測點(diǎn)1 分別位于裂紋尖端處和圓孔的斜下方.

圖9 含裂紋板和含圓孔板模型Fig. 9 The models for plates containing a crack or a hole

該問題中，PD 物質(zhì)點(diǎn)之間離散間距為Δx=Δy=0.05 m，Δt=0.05 d.如圖10 和圖11 所示，將有限元軟件ANSYS 與PDDO 方法模擬的穩(wěn)態(tài)結(jié)果進(jìn)行對比，結(jié)果表明通過兩種方法得到的含裂紋板的溫度分布基本一致.

圖10 穩(wěn)態(tài)時(shí)，含裂紋板在不同數(shù)值方法下的溫度分布：（a） FEM；（b） PDDOFig. 10 Temperature distributions of plates containing cracks for different numerical methods in steady states: （a） FEM; （b） PDDO

圖11 穩(wěn)態(tài)時(shí)，含圓孔板在不同數(shù)值方法下的溫度分布：（a） FEM；（b） PDDOFig. 11 Temperature distributions of plates containing holes for different numerical methods in steady states: （a） FEM; （b） PDDO

圖12 給出了含裂紋板和含圓孔板通過ANSYS 和本文模型計(jì)算的溫度變化曲線.其中圖12(a)和12(c)給出了含裂紋板和含圓孔板在測點(diǎn)1 處溫度隨時(shí)間的變化曲線.圖12(b)給出了達(dá)到穩(wěn)態(tài)時(shí)，溫度沿裂紋下邊界變化的曲線.上述圖中運(yùn)用本文模型與ANSYS 模擬的結(jié)果吻合，證明了本文建立的二維瞬態(tài)熱傳導(dǎo)模型具有較高精度且能夠處理含非連續(xù)結(jié)構(gòu)的熱傳導(dǎo)問題.

圖12 含裂紋板和含圓孔板在不同數(shù)值方法下的溫度變化曲線：（a） 含裂紋板測點(diǎn)處，溫度隨時(shí)間的變化曲線；（b） 穩(wěn)態(tài)時(shí)，溫度沿裂紋的下邊界變化的曲線；（c） 含圓孔板測點(diǎn)處，溫度隨時(shí)間的變化曲線Fig. 12 The temperature variation curves of plates containing a crack or a hole for different numerical methods: （a） temperature change curves with time at the measuring point of the plate containing a crack; （b）temperature curves along the lower boundary of the crack in a steady state; （c） temperature change curves with time at the measuring point of the plate containing a hole

3 結(jié) 論

本文將PDDO 方法應(yīng)用于求解了二維瞬態(tài)熱傳導(dǎo)問題.將二維瞬態(tài)熱傳導(dǎo)方程和邊界條件由其局部微分形式轉(zhuǎn)化為非局部積分形式，引入Lagrange 乘數(shù)法并通過變分分析，求解了二維瞬態(tài)熱傳導(dǎo)問題的溫度場.典型算例表明，基于PDDO 方法構(gòu)建二維瞬態(tài)熱傳導(dǎo)的非局部求解模型，采用均勻布點(diǎn)離散求解域，避免了由于網(wǎng)格依賴性和網(wǎng)格加密而產(chǎn)生的求解效率低的問題，不僅對求解不規(guī)則邊界板的瞬態(tài)熱傳導(dǎo)問題具有較高精度，而且也同樣適用于含裂紋和含圓孔板等不連續(xù)問題.算例結(jié)果顯示，PDDO 方法計(jì)算精度高、適用范圍廣，且具有較好的穩(wěn)定性和收斂性，為求解二維瞬態(tài)熱傳導(dǎo)問題提供了新的思路，同時(shí)也拓寬了PDDO 方法的應(yīng)用范圍.

參考文獻(xiàn)（ References ） ：

［1］王 紅, 李小林. 二維瞬態(tài)熱傳導(dǎo)問題的無單元Galerkin法分析[J]. 應(yīng)用數(shù)學(xué)和力學(xué), 2021, 42（5）: 460-469. （WANG Hong, LI Xiaolin. Analysis of 2D transient heat conduction problems with the element-free Galerkin method[J].Applied Mathematics and Mechanics, 2021, 42（5）: 460-469.（in Chinese））

［2］Z IENIUK E, SAWICKI D. Modification of the classical boundary integral equation for two dimensional transient heat conduction with internal heat source, with the use of NURBS for boundary modeling[J].Journal of Heat Transfer, 2017, 139（8）: 81-95.

［3］B URLAYENKO V N, ALTENBACH H, SADOWSKI T, et al. Modelling functionally graded materials in heat transfer and thermal stress analysis by means of graded finite elements[J].Applied Mathematical Modelling,2017, 45（5）: 422-438.

［4］W U X H, TAO W Q. Meshless method based on the local weak-forms for steady state heat conduction problems[J].International Journal of Heat and Mass Transfer, 2008, 51（11/12）: 3103-3112.

［5］B REBBIA C A, TELLES J C F, WROBEL L C.Boundary Element Techniques:Theory and Applications in Engineering[M]. Berlin: Springer Verlag, 1984.

［6］師 晉紅, 傅卓佳, 陳文. 邊界節(jié)點(diǎn)法計(jì)算二維瞬態(tài)熱傳導(dǎo)問題[J]. 應(yīng)用數(shù)學(xué)和力學(xué), 2014, 35（2）: 111-120. （SHI Jinhong, FU Zhuojia, CHEN Wen. Boundary knot method for 2D transient heat conduction problems[J].Applied Mathematics and Mechanics, 2014, 35（2）: 111-120.（in Chinese））

［7］P ARTRIDGE P W, BREBBIA C A, WROBEL L C.The Dual Reciprocity Boundary Element Method[M].Southampton: Computational Mechanics Publications, 1992.

［8］J IROUSEK J, QIN Q H. Application of hybrid-Trefftz element approach to transient heat conduction analysis[J].Computers & Structures, 1996, 58（1）: 195-201.

［9］S ILLING S A, EPTON M, WECKNER O, et al. Peridynamic states and constitutive modeling[J].Journal of Elasticity, 2007, 88（2）: 151-184.

［10］黃 丹, 章青, 喬丕忠, 等. 近場動(dòng)力學(xué)方法及其應(yīng)用[J]. 力學(xué)進(jìn)展, 2010, 40（4）: 448-459. （HUANG Dan, ZHANG Qing, QIAO Pizhong, et al. A review on peridynamics(PD)method and its applications[J].Advances in Mechanics, 2010, 40（4）: 448-459.（in Chinese））

［11］M ADENCI E, BARYT A, DORDUNCU M.Peridynamic Differential Operator for Numerical Analysis[M].Switzerland: Springer Nature Switzerland AG, 2019.

［12］M ADENCI E, BARUT A, FUTCH M. Peridynamic differential operator and its applications[J].Computer Methods in Applied Mechanics and Engineering, 2016, 304: 408-451.

［13］B AZAZZADEH S, SHOJAEI A, ZACCARIOTTO M, et al. Application of the peridynamic differential operator to the solution of sloshing problems in tanks[J].Engineering Computations, 2018, 36（1）: 45-83.

［14］D ORDUNCU M. Stress analysis of laminated composite beams using refined zigzag theory and peridynamic differential operator[J].Composite Structures, 2019, 218: 193-203.

［15］D ORDUNCU M. Stress analysis of sandwich plates with functionally graded cores using peridynamic differential operator and refined zigzag theory[J].Thin-Walled Structures, 2020, 146: 106468.

［16］G AO Y, OTERKUS S. Nonlocal modeling for fluid flow coupled with heat transfer by using peridynamic differential operator[J].Engineering Analysis With Boundary Elements, 2019, 105: 104-121.

［17］G AO Y, OTERKUS S. Nonlocal numerical simulation of low Reynolds number laminar fluid motion by using peridynamic differential operator[J].Ocean Engineering, 2019, 179: 135-158.

［18］L I Z Y, HUANG D, XU Y P, et al. Nonlocal steady-state thermoelastic analysis of functionally graded materials by using peridynamic differential operator[J].Applied Mathematical Modelling, 2021, 93: 294-313.

［19］C UI M, XU B B, FENG W Z, et al. A radial integration boundary element method for solving transient heat conduction problems with heat sources and variable thermal conductivity[J].Numerical Heat Transfer, 2018, 73:1-18.

［20］S OLEIMANI S, JALAAL M, BARARNIA H, et al. Local RBF-DQ method for two dimensional transient heat conduction problems[J].International Communications in Heat and Mass Transfer, 2010, 37（9）: 1411-1418.

［21］M UKHERJEE Y X, MUKHERJEE S. On boundary conditions in the element-free Galerkin method[J].Computational Mechanics, 1997, 19（4）: 264-270.