無網格對稱粒子法中兩類熱邊界條件的處理

肖毅華，張浩鋒，平學成

（華東交通大學機電工程學院，江西南昌330013）

無網格對稱粒子法中兩類熱邊界條件的處理

肖毅華，張浩鋒，平學成

（華東交通大學機電工程學院，江西南昌330013）

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

瞬態熱傳導；邊界處理；無網格對稱粒子法

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

1 無網格對稱粒子法的基本方程

在求解瞬態熱傳導問題時，無網格對稱粒子法采用一組粒子離散問題域。每個粒子攜帶密度ρ，質量m和溫度T等物理量。物理量的導數通過對稱粒子近似進行計算。為了得到某一物理量f的一階導數的對稱粒子近似，對于任意粒子i，在其鄰域內定義如下誤差項

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

式中

利用式（3）對瞬態熱傳導問題的控制方程離散，得到無網格對稱粒子法的基本方程為

式中：t為時間，c為比熱，λx、λy和λz為3個坐標方向的導熱系數，為熱源強度。本文考慮各個方向的導熱系數相同，即λx=λy=λz=λ。

2 邊界條件處理

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

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

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

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

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

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

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

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

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

式中。根據上述過程求得的邊界粒子的溫度導數是精確滿足邊界條件的。值得指出的是，為了保證數值計算的穩定性，避免式（10）出現分母為0的情況，應選擇將方向余弦絕對值最大的方向的溫度導數用其余兩個方向的溫度導數表示出來。

對于II類邊界粒子，先根據式（3）和式（10）、（12）分別計算完內部粒子和I類邊界粒子的溫度導數，再根據已計算出的粒子溫度導數，通過具有一階連續性的粒子近似計算得到其溫度導數，計算表達式為

式中

3 數值算例

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

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

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

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

圖4給出了無網格對稱粒子法的計算結果與解析解的總體相對誤差Re1和單個點最大相對誤差Re2隨時間的變化。兩誤差的定義分別為

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

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

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

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

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

4 結論

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

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

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

[3]王冰冰,高欣,段慶林.穩態熱傳導問題的二階一致無網格法[J].應用數學與力學,2013,34(7):750-755.

[4]戴艷俊,吳學紅,陶文銓.三維不規則區域熱傳導問題無網格方法的數值模擬[J].工程熱物理學報,2011,32(7):1173-1177.

[5]吳學紅,李增耀,申勝平,等.不規則區域熱傳導問題無網格Ptrov-Galerkin方法的數值模擬[J].工程熱物理學報,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]蔣濤,歐陽潔,栗雪娟,等.瞬態熱傳導問題的一階對稱SPH方法模擬[J].物理學報,2011,60(9):090206-090206.

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

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

[11]張宏偉,李美香,李衛國.關于配點型無網格法邊界條件處理技術[J].大連理工大學學報,2010,50(4):614-618.

[12]傅里葉.熱的解析理論[M].北京:北京大學出版社,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

國家自然科學基金項目（11302077）；江西省“井岡之星”青年科學家培養計劃項目（20112BCB23013）；江西省教育廳科技項目（GJJ14398）

肖毅華（1984—），男，講師，博士，主要研究方向為無網格法及其應用。

1005-0523（2014）04-0065-06