SODF-SPH方法及其在熱傳導問題中的應用

熱合買提江·依明 阿合買提江·伊明江 買買提明·艾尼

1.新疆大學數學與系統科學學院，烏魯木齊，8300462.新疆大學機械工程學院，烏魯木齊，830047

SODF-SPH方法及其在熱傳導問題中的應用

熱合買提江·依明1阿合買提江·伊明江1買買提明·艾尼2

1.新疆大學數學與系統科學學院，烏魯木齊，8300462.新疆大學機械工程學院，烏魯木齊，830047

在光滑粒子流體動力學(SPH)原理的基礎上，通過泰勒級數展開提出了計算函數二階偏導數的SODF-SPH方法。選用相等的粒子間距和相同的光滑長度對一維熱傳導問題進行了計算和模擬，與傳統SPH方法進行了對比分析，結果表明新方法的精度高、收斂速度快且穩定性好；對二維和三維各種熱傳導問題進行了計算和模擬，與解析解進行了對比，結果表明新方法得到的計算結果與解析解吻合良好。新方法的計算過程能避免計算核函數導數，致使對核函數的要求降低，可選用更多的核函數且計算量較小，可在工程和數值計算中廣泛應用。

光滑粒子流體動力學方法；SODF-SPH方法；精度；收斂速度; 穩定性

0 引言

光滑粒子流體動力學(smoothed particle hydrodynamics，SPH)方法是一種新的無網格純 Lagrange 粒子方法，由LUCY等[1]和GINGOLD等[2]提出，最初用于解決天體物理學和宇宙學中模擬三維無界空間天體的演化問題，之后他們又提出了適用于SPH方法的人為黏性和人為熱流項[3-4]，并延伸到了計算流體力學的各個領域。LIBERSKY等[5]在SPH方法中引入了材料強度模型，模擬了動態斷裂和破片運動等固體的動態響應，此后SPH也被迅速地應用到固體力學的研究當中。目前，SPH方法已廣泛應用于工程計算和模擬的各個領域[6-8]。

隨著SPH方法應用范圍的擴展和對其研究的不斷深入，SPH方法的計算精度和穩定性等方面的缺陷被人們發現，于是研究者用各種方法構造了不同類型的SPH方法，彌補了經典SPH方法的一些缺點[9-11]。SPH方法在應用于熱傳導問題[12-14]的研究時，二階偏微分方程被分解為兩個一階偏微分方程，而本文在研究SPH方法的基礎上，基于泰勒級數展開構造了直接計算函數二階偏導數的新SODF-SPH方法，并通過計算和模擬各種熱傳導問題，來驗證SODF-SPH方法解決各種熱傳導問題的可行性。

1 經典SPH方法

經典SPH方法基本原理比較簡單，將任意函數f(x)及其梯度在某一點x處的SPH核近似式，用核函數W(x-x′,h)與f(x)乘積的積分來表示。函數f(x)的SPH核近似表達式為

(1)

式中，〈〉表示函數的核積分近似表達式。

(2)

核函數W(x-x′,h)依賴于兩點之間的距離|x-x′|以及核函數影響區域Ω的光滑長度h，因此核函數W(x-x′,h)需滿足若干個基本條件[15]。

在SPH方法中把整個求解域離散成一系列任意分布的有質量和容量的有限多個粒子,如圖1所示(二維情況)。如果核近似式中粒子j處微元體dx′的體積為ΔVj，質量為mj，密度為ρj，位置坐標為xj，記函數f(x)在粒子j處的值為fj=f(xj)，假設粒子i為中心的核函數影響范圍內的粒子總數為M，則函數f(x)在粒子i處的函數及其梯度的粒子和近似式可寫為

(3)

(4)

圖1 SPH核近似計算示意圖Fig.1 SPH kernel approximation

2 SODF-SPH方法的推導

函數f(x) 在x=xi處的泰勒級數展開為

(5)

式中,o(·)表示高階無窮小。

(6)

由核函數對稱性質，式(6)等號右邊第二項積分為零。此時整理式(6)得

(7)

由SPH方法的粒子和近似表達式原理與式(7)可得在點x=xi計算函數二階偏導數近似值的線性方程：

(8)

xji=xj-xi

fji=fj-fi

由以上推導可知，本方法的計算過程可直接代入核函數而無需求其導數，降低了對核函數的要求，使得可選的核函數范圍更廣；二階偏導數的計算不必分解成兩個一階偏導數，使得計算速度加快。

3 瞬態熱傳導方程

在直角坐標系下，瞬態熱傳導問題的控制方程為

(9)

式中，T為溫度;ρ和c分別為材料密度和材料質量熱容；λα為沿α方向的熱導率；Q為熱源強度。

初始條件：T(x,0)=T0，x∈Ω。邊界條件：T(x,t)=T1，x∈?ΓD(Dirchlet邊界條件)；T(x,t)n+T(x,t)b=0，x∈?ΓN(Neumann邊界條件)。其中，n為區域Ω邊界外單位向量，T(x,t)b為邊界熱源，?ΓD∩?ΓN=?，?ΓD∪?ΓN=?Ω，?Ω為區域邊界。

4 瞬態熱傳導問題的SODF-SPH方法

模擬

4.1 誤差判斷準則

本文通過熱傳導問題的模擬，將新方法與經典SPH方法對比，對新方法的精度、收斂性和穩定性進行了分析。對不同類型的方法進行比較時，本文以問題解析解與數值解的誤差作為判斷準則。誤差為

(10)

式中，M為整個計算區域內的粒子總數；T(xi)為xi處的解析解；TC(xi)為xi處的數值解。

4.2 熱傳導問題的模擬

本文所有實例取ρ=1，c=1，λα=1，核函數光滑長度h=1.1Δd(Δd為初始粒子間距)。

4.2.1 一維熱傳導問題的模擬

實例1 初始和邊界條件為T(x，0)=T0，T(0,t)=T0，T(L,t)=T1，此時該問題的解析解[12]為

模擬時,取L=1，T0=0 ℃，T1=1 ℃。

圖2給出了設立21個粒子，時間步長t=10-4s時，不同時刻溫度曲線的解析解與數值解的對比，結果表明本文方法較經典SPH方法具有更高的精度和穩定性。圖3給出了不同時刻分別設立21、41、81、161個粒子時，本文方法與經典SPH方法誤差的變化，通過對比發現SODF-SPH方法比經典SPH方法具有更好的收斂性。

圖2 不同時刻溫度曲線數值解Fig.2 Numerical solutions at different time step

4.2.2 二維熱傳導問題的模擬

實例2 恒溫邊界條件熱傳導問題的模擬：當T(x,y,0)=T0,T(0,y,t)=T(L,y,t)=T(x,0,t)=T(x,H,t)=T1時，問題的解析解[13]為

模擬中取L=1，H=1，T0=100 ℃，T1=0 ℃，模擬中設立51×51個均勻分布粒子，時間步長dt=10-4s。

(a)SPH方法

(b)SODF-SPH方法圖3 不同時刻不同粒子數的誤差對比Fig.3 Errors at different time step and different particle numbers

圖4為t=0.1 s時解析解與SODF-SPH數值解的云圖，對比發現數值解與解析解幾乎完全一樣。為了進一步驗證數值計算結果的正確性，圖5為不同時刻在直線y=x上解析解與SODF-SPH數值結果的對比圖，由不同時刻的數值結果對比可以看出，計算結果與解析解吻合良好。

式中，αk、βl分別為超非線性方程αktan(αkL)=b和βltan(βlH)=b的非負解。

模擬中取L=1，H=1，T0=100 ℃，b=0.1，模擬中設立51×51個均勻分布粒子，時間步長dt=10-4s。

(a)解析解

(b)SODF-SPH數值解圖4 當t=0.1 s時的溫度分布情況(恒溫邊界問題)Fig.4 Temperature distribution for the constant temperature boundaries at t=0.1 s

圖5 不同時刻溫度變化情況(y=x)Fig.5 Temperature distribution at different time step (y=x)

圖6為t=0.05 s時解析解與SODF-SPH數值解的云圖，結果顯示絕熱邊界和對流邊界熱傳導問題的溫度等值線分布情況基本一致。此時本文提出方法得到的最高和最低溫度分別是Tmax=99.993 205 ℃和Tmin=95.000 652 ℃，而相應的解析解分別為Tmax=99.994 653 ℃和Tmin=95.113 279 ℃，兩者也非常接近。以上結果表明用本文提出方法模擬混合邊界條件瞬態熱傳導問題的可靠性。

(a)解析解

(b)SODF-SPH數值解圖6 當t=0.05 s時的溫度分布情況(對流邊界問題)Fig.6 Temperature distribution for the convection boundaries at t=0.05 s

4.2.3 三維熱傳導問題的模擬

實例4 初始條件為T(x,y,z,0)=sin(πx)+sin(πy)+sin(πz)；邊界條件為T(0,y,z,t)=T(1,y,z,t)=(sin(πy)+sin(πz))exp(-π2t)，T(x,0,z,t)=T(x,1,z,t)=(sin(πx)+sin(πz))exp(-π2t)，T(x,y,0,t)=T(x,y,1,t)=(sin(πx)+sin(πy))·exp(-π2t)時，該問題的解析解[16]為

T(x,y,z,t)=(sin(πx)+sin(πy)+
sin(πz))exp(-π2t)

模擬中設立26×26×26個均勻分布粒子，時間步長dt=10-6s。

(a)解析解

(b)SODF-SPH數值解圖7 當t=0.001 s時的溫度分布情況(三維問題)Fig.7 Temperature distribution for the 3D case att=0.001 s

圖7給出了當t=0.001 s時的解析解與SODF-SPH數值解，此時本文方法得到的最高溫度是Tmax=2.9647 ℃，而相應的解析解為Tmax=2.9516 ℃。因此，對三維熱傳導問題用SODF-SPH方法所得結果與解析解基本一致。

5 結論

(1)在設立相同的粒子間距、核函數和光滑長度的情況下計算函數二階偏導數，本文方法比經典SPH方法誤差小、穩定性好且收斂速度快。

(2) 用SODF-SPH方法對一維、二維和三維瞬態熱傳導問題在不同的邊界條件下進行計算和模擬，數值解與解析解吻合良好，驗證了該方法模擬瞬態熱傳導問題的可行性。

(3)用SODF-SPH方法進行計算，對核函數導數的要求較低，使得可選的核函數更多，且二階偏導數的計算速度更快。這一方法具有良好的計算精度和求解能力，有望在復雜的熱傳導等工程問題中得到應用。

[1]LUCYLB,ANumericalApproachtotheTestingofFissionHypothesis[J].AstrophysicalJournal, 1977, 82(12):1013-1020.

[2]GINGOLDRA,MONAGHANJJ.SmoothedParticleHydrodynamics:TheoryandApplicationtoNon-sphericalStars[J].MonthlyNoticesoftheRoyalAstronomicalSociety, 1977, 181(3):375-389.

[3]GINGOLDRA,MONAGHANJJ.KernelEstimatesasaBasisforGeneralParticleMethodsinHydrodynamics[J].JournalofComputationalPhysics, 1982, 46(3):429-453.

[4]MONAGHANJJ.SimulatingFreeSurfaceFlowswithSPH[J].JournalofComputationalPhysics, 1994, 110(2):399-406.

[5]LIBERSKYLD,PETSCHEKAG.SmoothParticleHydrodynamicswithStrengthofMaterials[J].LectureNotesinPhysics, 2006, 395:248-257.

[6]STOWED,KUPSCHELLAR,PANH,etal.InvestigationofS-SPHforHypervelocityImpactCalculations[J].ProcediaEngineering, 2015, 103:585-592.

[7]NAIRP,TOMARG.AnImprovedFreeSurfaceModelingforIncompressibleSPH[J].Computers&Fluids, 2014, 102:304-314.

[8]CHENZ,ZONGZ,LIUMB,etal.AnSPHModelforMultiphaseFlowswithComplexInterfacesandLargeDensityDifferences[J].JournalofComputationalPhysics, 2015, 283(C):169-188.

[9]CHENJK,BERAUNJE,CARNEYTC.ACorrectiveSmoothedParticleMethodforBoundaryValueProblemsinHeatConduction[J].InternationalJournalforNumericalMethodsinEngineering, 1999, 46(2):231-252.

[10]ZHANGGM,BATRARC.ModifiedSmoothedParticleHydrodynamicsMethodandItsApplicationtoTransientProblems[J].ComputationalMechanics, 2004, 34(2):137-146.

[11]ZHANGGM,BATRARC.SymmetricSmoothedParticleHydrodynamics(SSPH)MethodandItsApplicationtoElasticProblems[J].ComputationalMechanics, 2009, 43(3):321-340.

[12]JEONGJH,JHONMS,HALOWJS,etal.SmoothedParticleHydrodynamics:ApplicationstoHeatConduction[J].ComputerPhysicsCommunications, 2003, 153(1):71-84.

[13]ZHANGX,ZHANGP,ZHANGL.AnImprovedMeshlessMethodwithAlmostInterpolationPropertyforIsotropicHeatConductionProblems[J].EngineeringAnalysiswithBoundaryElements, 2013, 37(5):850-859.

[14] 蔣濤, 歐陽潔, 栗雪娟,等. 瞬態熱傳導問題的一階對稱SPH方法模擬[J]. 物理學報, 2011, 60(9):49-58.JIANGTao,OUYANGJie,LIXuejuan,etal,TheFirstOrderSymmetricSPHMethodforTransientHeatConductionProblems[J].ActaPhysicaSinica, 2011, 60(9):49-58.

[15]LIUGR,LIUMB.SmoothedParticleHydrodynamics:aMeshlessParticleMethod[M].Singapore:WorldScientificPublishing, 2003:150-156.

[16]CHANGCW,LIUCS.ANewAlgorithmforDirectandBackwardProblemsofHeatConductionEquation[J].InternationalJournalofHeat&MassTransfer, 2010, 53(23/24):5552-5569.

(編輯 王旻玥)

SODF-SPH Method and Its Applications in Heat Conduction Problems

Rahmatjan IMIN1Ahmatjan IMINJAN1Mamtimin GHENI2

1.College of Mathematics and System Science, Xinjiang University,Urumqi,830046 2.School of Mechanical Engineering, Xinjiang University,Urumqi,830047

Based on the principles of SPH method, a SODF-SPH method to compute second order derivatives was constructed through Taylor series expansion. Equal particle distance and same smoothed length were chosen to compute and simulate 1D heat conduction problems. After compared with conventional SPH method, results show that the accuracy, convergence rate and stability of the new method are better. To compute and simulate 2D and 3D heat conduction problems and compared with analytical solutions, the numerical results are in accord with analytical solutions. The derivatives may avoid calculating differential coefficients of kernel function. Therefore, the demands of kernel function are reduced, more kernel functions may be used and the calculation amounts are decreased. The new method may be widely used in engineering fields and numerical calculations.

smoothed particle hydrodynamics (SPH) method; second order derivative free SPH (SODF-SPH) method; accuracy; convergence rate; stability

2016-06-16

國家自然科學基金資助項目( 51565054, 51075346)；新疆大學博士畢業生科研啟動基金資助項目(BS150210)

TH212；TH213.3

10.3969/j.issn.1004-132X.2017.04.007

熱合買提江·依明，男，1974年生。新疆大學數學與系統科學學院副教授、博士。主要研究方向為機械優化設計、工程數值計算及其模擬。發表論文20余篇。E-mail:rehmatjanim@xju.edu.cn。阿合買提江·伊明江，男，1969年生。新疆大學數學與系統科學學院副教授。買買提明·艾尼，男，1958年生。新疆大學機械工程學院教授、博士研究生導師。