無網(wǎng)格對稱粒子法中兩類熱邊界條件的處理

肖毅華，張浩鋒，平學(xué)成

（華東交通大學(xué)機(jī)電工程學(xué)院，江西南昌330013）

無網(wǎng)格對稱粒子法中兩類熱邊界條件的處理

肖毅華，張浩鋒，平學(xué)成

（華東交通大學(xué)機(jī)電工程學(xué)院，江西南昌330013）

第二、三類邊界條件對無網(wǎng)格對稱粒子法求解瞬態(tài)熱傳導(dǎo)問題的計(jì)算精度和穩(wěn)定性有重要影響。提出一種有效的方法來處理這兩類邊界條件。該方法采用統(tǒng)一表達(dá)式描述兩類邊界條件，并將邊界粒子分為光滑域內(nèi)的邊界粒子和非光滑域的邊界粒子。對于前者，在溫度導(dǎo)數(shù)的對稱粒子近似中引入邊界條件方程，精確地施加邊界條件；對于后者，根據(jù)其相鄰粒子的溫度導(dǎo)數(shù)，采用具有一階連續(xù)性的粒子近似，計(jì)算出其溫度導(dǎo)數(shù)。應(yīng)用該方法計(jì)算了一個(gè)數(shù)值算例。計(jì)算結(jié)果與解析解和有限元法的計(jì)算結(jié)果符合良好，且計(jì)算誤差隨時(shí)間呈下降趨勢，說明該方法具有良好的精度和穩(wěn)定性。

瞬態(tài)熱傳導(dǎo)；邊界處理；無網(wǎng)格對稱粒子法

熱傳導(dǎo)問題［1-2］普遍存在于機(jī)械、材料、土木等研究領(lǐng)域。由于實(shí)際熱傳導(dǎo)問題的復(fù)雜性，人們往往需要借助數(shù)值方法進(jìn)行研究。常用的數(shù)值方法是基于網(wǎng)格的方法，如有限元法［1］、有限體積法［2］。無網(wǎng)格法發(fā)展較晚，目前在熱傳導(dǎo)問題中也得到了一定的應(yīng)用［3-8］。相對于基于網(wǎng)格的方法，無網(wǎng)格法顯示出一些優(yōu)勢。它不涉及網(wǎng)格，可以避免網(wǎng)格畸變，且無需復(fù)雜和耗時(shí)的網(wǎng)格劃分等步驟。但是，多數(shù)無網(wǎng)格法需要復(fù)雜的數(shù)值積分和形函數(shù)計(jì)算，計(jì)算量大、效率低。無網(wǎng)格法中的粒子法，如光滑粒子法（SPH），以帶有質(zhì)量、密度等物理量的粒子模擬物體，并用配點(diǎn)法離散控制微分方程，無需數(shù)值積分，具有概念簡單、程序?qū)崿F(xiàn)容易、計(jì)算效率較高等優(yōu)點(diǎn)。近年來，一些研究者［9-10］探討了此類方法在熱傳導(dǎo)問題中的應(yīng)用。由于無網(wǎng)格粒子法是配點(diǎn)型無網(wǎng)格法，該類方法對邊界特別是導(dǎo)數(shù)邊界的處理很敏感。因此，邊界條件的處理是關(guān)于此類方法研究的一個(gè)熱點(diǎn)［11］。在熱傳導(dǎo)問題中，第二、三類邊界條件都是導(dǎo)數(shù)邊界條件，其處理對于保證無網(wǎng)格粒子法的精度和穩(wěn)定性具有重要意義。

1 無網(wǎng)格對稱粒子法的基本方程

在求解瞬態(tài)熱傳導(dǎo)問題時(shí)，無網(wǎng)格對稱粒子法采用一組粒子離散問題域。每個(gè)粒子攜帶密度ρ，質(zhì)量m和溫度T等物理量。物理量的導(dǎo)數(shù)通過對稱粒子近似進(jìn)行計(jì)算。為了得到某一物理量f的一階導(dǎo)數(shù)的對稱粒子近似，對于任意粒子i，在其鄰域內(nèi)定義如下誤差項(xiàng)

式中：ξ=d h；d為粒子間的距離；h為影響域半徑。令式（1）定義的誤差最小化并求解得到的線性方程組，即給出的對稱粒子近似式為

式中

利用式（3）對瞬態(tài)熱傳導(dǎo)問題的控制方程離散，得到無網(wǎng)格對稱粒子法的基本方程為

式中：t為時(shí)間，c為比熱，λx、λy和λz為3個(gè)坐標(biāo)方向的導(dǎo)熱系數(shù)，為熱源強(qiáng)度。本文考慮各個(gè)方向的導(dǎo)熱系數(shù)相同，即λx=λy=λz=λ。

2 邊界條件處理

熱傳導(dǎo)問題涉及3類邊界條件：第一類為給定邊界上的溫度值，第二類為給定邊界上的熱流密度分布，第三類為給定邊界上物體與周圍介質(zhì)間的對流換熱系數(shù)及介質(zhì)溫度。3類邊界條件分別表述如下

式中：Ts為物體邊界面上的溫度，為給定的邊界面上的溫度分布函數(shù)；n為物體邊界面的外法向量，為給定的邊界面上的熱流密度分布函數(shù)；α為對流換熱系數(shù)，Ta為周圍介質(zhì)的溫度。

對于第一類邊界條件，無網(wǎng)格對稱粒子法可以直接施加，即在每一時(shí)間步內(nèi)將要滿足此種邊界條件的粒子賦予指定的溫度值。第二類和第三類邊界條件都包含導(dǎo)數(shù)，其處理相對比較困難。將這2種邊界條件寫為如下統(tǒng)一形式

式中：cosθx，cosθy和cosθz為邊界面外法向量的方向余弦。

式（8）需要計(jì)算邊界面外法向量的方向余弦，而外法向量在不光滑的邊界位置（如圖1所示）無法確定。因此，為了避免計(jì)算困難，本文將滿足第二類和第三類邊界條件的粒子分為兩類：光滑域內(nèi)的邊界粒子為I類邊界粒子，非光滑域的邊界粒子為II類邊界粒子。

對于I類邊界粒子，根據(jù)邊界條件方程式（8）易得

式中：A=-cosθyi/cosθxi，B=-cosθzi/cosθxi，C=κ/cosθxi。利用上式可將兩粒子處的溫度差近似表示為

圖1 邊界粒子的分類Fig.1 Classification of boundary particles

再將上式代入誤差定義式（1）并令其最小化，求解得到

式中。根據(jù)上述過程求得的邊界粒子的溫度導(dǎo)數(shù)是精確滿足邊界條件的。值得指出的是，為了保證數(shù)值計(jì)算的穩(wěn)定性，避免式（10）出現(xiàn)分母為0的情況，應(yīng)選擇將方向余弦絕對值最大的方向的溫度導(dǎo)數(shù)用其余兩個(gè)方向的溫度導(dǎo)數(shù)表示出來。

對于II類邊界粒子，先根據(jù)式（3）和式（10）、（12）分別計(jì)算完內(nèi)部粒子和I類邊界粒子的溫度導(dǎo)數(shù)，再根據(jù)已計(jì)算出的粒子溫度導(dǎo)數(shù)，通過具有一階連續(xù)性的粒子近似計(jì)算得到其溫度導(dǎo)數(shù)，計(jì)算表達(dá)式為

式中

3 數(shù)值算例

下面通過一個(gè)數(shù)值算例來驗(yàn)證上節(jié)提出的邊界處理算法的有效性。算例考慮邊長為a=2的立方形物體在兩種不同邊界條件下的溫度場變化。計(jì)算時(shí)采用無量綱化參數(shù)，物體的相關(guān)熱物理參數(shù)為密度ρ=1；比熱c=1；導(dǎo)熱系數(shù)λx=λy=λz=1；物體的初始溫度為1；熱源強(qiáng)度=0。求解區(qū)域的原點(diǎn)設(shè)在其中心處，采用均勻分布的粒子對其離散；粒子間距為0.05，粒子總數(shù)為68 921；粒子的影響域半徑為1.5倍粒子間距。模擬總時(shí)間設(shè)為0.2，計(jì)算時(shí)間步長取為0.001。

首先考慮在立方體表面作用第三類邊界條件，即立方體的外表面與周圍介質(zhì)對流換熱，設(shè)對流換熱系數(shù)為α=1，周圍介質(zhì)溫度為Ta=0。在此邊界條件下，可以獲得該問題的解析解［12］，從而便于檢驗(yàn)算法的正確性和精度。圖2給出了3個(gè)特征點(diǎn)A、B、C的溫度變化歷程，這3個(gè)點(diǎn)的坐標(biāo)分別為（0，0，0）、（0.5，0.5，0.5）和（1.0，1.0，1.0）。可見，無網(wǎng)格對稱粒子法計(jì)算的特征點(diǎn)溫度與解析解在整個(gè)模擬時(shí)間段內(nèi)都能良好地吻合。圖3為時(shí)間t=0.2時(shí)立方體上表面（z=1.0）的溫度等值線圖。無網(wǎng)格對稱粒子法計(jì)算的溫度分布與解析解基本一致。

圖2 第三類邊界條件下3個(gè)特征點(diǎn)的溫度變化Fig.2 Temperature histories of three representative points for boundary conditions of the third kind

圖3 t=0.2時(shí)上表面（z=1.0）的溫度分布對比Fig.3 Comparison of temperature distributions of the top surface（z=1.0）att=0.2

圖4給出了無網(wǎng)格對稱粒子法的計(jì)算結(jié)果與解析解的總體相對誤差Re1和單個(gè)點(diǎn)最大相對誤差Re2隨時(shí)間的變化。兩誤差的定義分別為

式中：Tp(xi)和Tf(xi)為分別為根據(jù)粒子法和解析解計(jì)算的粒子點(diǎn)xi處的溫度值，n為粒子數(shù)。由圖4可見，總體相對誤差始終小于1%，且隨模擬時(shí)間的增加而減少；單個(gè)點(diǎn)的最大相對誤差在開始時(shí)為8.4%，隨后也呈遞減趨勢。以上結(jié)果說明了本文方法對于第三類邊界條件具有良好的精度和穩(wěn)定性。

圖4 第三類邊界條件下的總體相對誤差和單個(gè)節(jié)點(diǎn)最大相對誤差Fig.4 Relative error of the whole region and themaximum relative error at single points for boundary conditions of the third kind

下面考慮在立方體表面作用第二類邊界條件的情況。設(shè)立方體各表面上都有均勻的熱量輸入，熱流密度的大小為=1。在此邊界條件下，難以得到解析解。為了驗(yàn)證邊界處理算法的有效性，建立單元尺寸與粒子間距相同的有限元模型，以有限元法的求解結(jié)果作為參考依據(jù)，將無網(wǎng)格對稱粒子法的結(jié)果與其對比。圖5給出了與前述位置相同的3個(gè)特征點(diǎn)的溫度變化歷程。兩種方法的計(jì)算結(jié)果在整個(gè)模擬時(shí)間段內(nèi)也都符合良好。圖6為兩種方法計(jì)算結(jié)果的總體相對誤差和單個(gè)節(jié)點(diǎn)最大相對誤差的變化曲線。兩種誤差均隨模擬時(shí)間的增加而減少，總體相對誤差的最大值僅為0.6%，單個(gè)節(jié)點(diǎn)最大相對誤差的最大值約為6%。以上結(jié)果可以說明本文方法對于第二類邊界條件也是可行和準(zhǔn)確的。

圖5 第二類邊界條件下三個(gè)特征點(diǎn)的溫度變化Fig.5 Temperature histories of three representative points for boundary conditions of the second kind

圖6 第二類邊界條件下的總體相對誤差和單個(gè)節(jié)點(diǎn)最大相對誤差Fig.6 Relative error of the whole region and themaximum relative error at single points for boundary conditions of the second kind

4 結(jié)論

研究了無網(wǎng)格對稱粒子法求解瞬態(tài)熱傳導(dǎo)問題時(shí)第二類和第三類邊界條件的處理，提出了一種統(tǒng)一的方法施加這兩類邊界條件。數(shù)值計(jì)算的結(jié)果表明該方法是十分有效的，能使無網(wǎng)格對稱粒子法準(zhǔn)確、穩(wěn)定地求解具有第二類和第三類邊界條件的瞬態(tài)熱傳導(dǎo)問題。同時(shí)，該方法無需添加輔助節(jié)點(diǎn)或在邊界附近進(jìn)行節(jié)點(diǎn)加密等，實(shí)施過程較簡單，可以為無網(wǎng)格對稱粒子法模擬復(fù)雜熱邊界條件提供了一種有效的途徑，促進(jìn)其在解決實(shí)際工程問題中的應(yīng)用。

[1]倪昀,黃志超,熊國良,等.基于ANSYS汽車后橋殼體焊接溫度場有限元分析[J].華東交通大學(xué)學(xué)報(bào),2006,23(2):115-118.

[2]劉仍通,潘陽.自然對流影響結(jié)冰的數(shù)值模擬以及實(shí)驗(yàn)研究[J].華東交通大學(xué)學(xué)報(bào),2010,27(5):22-27.

[3]王冰冰,高欣,段慶林.穩(wěn)態(tài)熱傳導(dǎo)問題的二階一致無網(wǎng)格法[J].應(yīng)用數(shù)學(xué)與力學(xué),2013,34(7):750-755.

[4]戴艷俊,吳學(xué)紅,陶文銓.三維不規(guī)則區(qū)域熱傳導(dǎo)問題無網(wǎng)格方法的數(shù)值模擬[J].工程熱物理學(xué)報(bào),2011,32(7):1173-1177.

[5]吳學(xué)紅,李增耀,申勝平,等.不規(guī)則區(qū)域熱傳導(dǎo)問題無網(wǎng)格Ptrov-Galerkin方法的數(shù)值模擬[J].工程熱物理學(xué)報(bào),2009,30(8): 1350-1352.

[6]SHIBAHARA M,ATLURI S N.Themeshless local Petrov-Galerkinmethod for the analysis of heat conduction due to amoving heat source,in welding[J].International Journal of Thermal Sciences,2011,50(6):984-992.

[7]CHENG R J,Ge H X.Meshless analysis of three-dimensional steady-state heat conduction problems[J].Chinese Physics B, 2010,19(9):090201-090201.

[8]蔣濤,歐陽潔,栗雪娟,等.瞬態(tài)熱傳導(dǎo)問題的一階對稱SPH方法模擬[J].物理學(xué)報(bào),2011,60(9):090206-090206.

[9]章文捷,王小惠,江順亮.瞬態(tài)熱傳導(dǎo)方程的SPH法[J].化學(xué)工程與裝備,2010(10):18-23.

[10]金文佳,章文捷,王小惠,等.熱傳導(dǎo)問題的SPH方法及其邊界處理[J].化學(xué)工程與裝備,2011(10):37-41.

[11]張宏偉,李美香,李衛(wèi)國.關(guān)于配點(diǎn)型無網(wǎng)格法邊界條件處理技術(shù)[J].大連理工大學(xué)學(xué)報(bào),2010,50(4):614-618.

[12]傅里葉.熱的解析理論[M].北京:北京大學(xué)出版社,2008.

Treatment of Two Kinds of Thermal Boundary Conditions in Meshless Symmetric Particle Method

Xiao Yihua,Zhang Haofeng,Ping Xuecheng
(School of Mechatronics Engineering,East China Jiaotong University,Nanchang 330013,China)

Boundary conditions of the second and the third kind have important effects on the accuracy and stabili?ty ofmeshless symmetric particlemethod for the solution of transient heat conduction.An effectivemethod is pre?sented in this paper to treat boundary conditions of both kinds.Two kinds of boundary conditions are expressed in a uniform form,and boundary particles are divided into two categories,the particles in smooth boundary regions and particles in non-smooth boundary regions.For the former,boundary conditions are enforced exactly by consid?ering the equations of boundary conditions as constraints in the symmetric particle approximation of derivatives of temperature.For the later,their derivatives of temperature are approximated with a particle approximation of firstorder consistency based on the derivatives of temperature of neighboring particles.The numerical results calculat?ed by a numerical example are in good agreement with analytical solutions and results obtained from the finite ele?mentmethod,which proves that the proposedmethod is accurate and stable.

transient heat conduction;boundary treatment;meshless symmetric particlemethod

O242

A

2014-05-25

國家自然科學(xué)基金項(xiàng)目（11302077）；江西省“井岡之星”青年科學(xué)家培養(yǎng)計(jì)劃項(xiàng)目（20112BCB23013）；江西省教育廳科技項(xiàng)目（GJJ14398）

肖毅華（1984—），男，講師，博士，主要研究方向?yàn)闊o網(wǎng)格法及其應(yīng)用。

1005-0523（2014）04-0065-06