一種改進的kgf-sph方法
【技術領域】
[0001] 本發明涉及一種改進的KGF-SPH方法，屬于計算流體力學技術領域。
【背景技術】
[0002] SPH (Smoothed Particle Hydrodynamics)是一種拉格朗日型無網格粒子方法，由 Lucy，Gingold與Monaghan于1977年提出并用于求解三維開放空間的天體物理學問題，現 在已經被廣泛用于多相流、重力流、熱傳導等流體力學問題的數值模擬中。在SPH方法中， 用一系列粒子來代替連續介質，這些粒子可以在計算域內自由移動，因此SPH方法可以很 好的處理傳統有網格方法難以解決的問題，如具有移動邊界、自由表面或邊界極大變形的 流動問題。雖然SPH方法已經在流體力學和固體力學領域有了非常廣泛的應用，但是SPH 方法在對流問題的應用仍較少。
[0003] 對流現象無處不在，如日常生活中的熱水沸騰，自然界中的海洋環流，工程應用 中的核反應堆冷卻等，因此研究對流問題具有非常的意義。但是傳統SPH方法存在粒子 不一致性、應力不穩定性等缺點，無法很好的模擬對流問題。在SPH方法的基礎上，Huang 等人在 2015 年提出了 一種 KGF-SPH 方法（kernel gradient free-smoothed particle hydrodynamics)。該方法在核近似和粒子近似過程中只用到了核函數本身，不包含核函數 的導數，因此無論核函數是否可導都可用于該方法，是一種無核梯度的SPH方法。KGF-SPH 方法降低了無網格方法對于核函數的要求，具有較高的精度。但是KGF-SPH方法中二階導 數通過對一階導數再求導來求解，這樣雖然提高了計算精度，但是存在以下缺點：
[0004] (1)計算量較大。因為二階導數的求解需要在對一階導數求導，增加了一倍的計算 量。
[0005] (2)不利于處理邊界條件。因為在求解邊界附近流體粒子的物理量時，需要先求解 出邊界上粒子物理量的導數值，然而為了精確的求出邊界上粒子物理量的導數值，還需要 在邊界外布置更多的虛粒子。
【發明內容】
[0006] 本發明的目的是為了提高數值模擬的精度和穩定性，減小計算量，提出了一種改 進的KGF-SPH方法。
[0007] 本發明是通過以下技術方案實現的：
[0008] -種改進的KGF-SPH方法，具體實現步驟如下：
[0009] 步驟1、在所研究問題的計算域中布置粒子；
[0010] 步驟2、根據所研究問題的初始條件和邊界條件對粒子的物理屬性進行初始化；
[0011] 步驟3、利用下述空間離散格式對所研究問題的控制方程進行近似，得到函數的時 間導數#的粒子近似式，其近似原理如下： dt
[0012] 對任意函數f (r])在點r =巧處進行泰勒展開，并且只保留到一階導數，可得：
[0013] f (r j) = f (r；) +Vf (r；) · Cr jTi) +〇 (h2) (I)
[0014] 忽略式⑴中的高階項o (h2)，在式⑴兩邊同時乘以核函數W(R，h)和（rfrj W (R, h)并在計算域上進行積分得：
[0015] / f(rj)ffdV = f (r;) f WdV+V f (r;) · f (r^ri) WdV (2)
[0016] f f (Tj) (T-Ti)WdV = f (r J f (rj-r^ffdV+V f (r;) · f (r ^ri) 2WdV (3)
[0017] 聯立式⑵和式（3)，并求解方程組，可得：
[0019] 上式對應的粒子近似式為：
[0021] 式中，t= f(r ;)，fj= f(r .)，（ ▽ :〇; = ▽ f(r J，N表示在i粒子的支持域內的 j粒子（即i粒子的最近相鄰粒子）總量，P ,分別表示j粒子的質量和密度。
[0022] 在一維空間下，式（5)的具體表達式為：
[0024] 其中，Xjl= Xj-X1Jlix為粒子i在X方向上的一階導數。修正矩陣厶立一維空間 下的表達式為：
[0026] 在一維空間下，對函數f (Xj)在點Xj= X 1進行二階泰勒展開并略去高階項可得：
[0028] 式中，fliXX為粒子i在X方向上的二階導數。式⑶兩邊同時乘以核函數W并移 項可得：
[0030] 式（9)兩邊同時積分可得：
[0032] 由核函數的對稱特性可知，上式右邊第一項為0,因此式（10)可化簡為：
[0034] 對（11)式進行移項可以推導出一維空間下拉普拉斯算子的離散格式為：
[0036] 式（12)對應的粒子近似式為：
[0038] 在二維空間下，式（5)的具體表達式為：
[0040] 其中，Xji= X fXi，Yji= yAjP f。分別為粒子i在X和y方向上的一階導 數。修正矩陣Ai在二維空間下的表達式為：
[0042] 在二維空間下，對函數f(Xj，yj)在點（Xl，yi)進行二階泰勒展開并略去高階項可 得：
[0044] 式中，fiiXy為粒子i的混合導數，f iiX)^P f iiyy分別為粒子i在X和y方向上的二階 導數。式（16)兩邊同時乘以核函數W并移項可得：
[0046] 式（17)兩邊同時積分得：
[0048] 由核函數的對稱特性可知，式（18)右邊第一項和第二項為0,因此上式可化簡為：
[0050] 當粒子分布均勻且規則時將
在極坐標下進行求解可得：
[0055] 式（19)移項并利用式（22)對其進行化簡，可以推導出二維空間下拉普拉斯算子 的離散格式為：
[0057] 式（23)對應的粒子近似式為：
[0059] 為了減弱粒子分布不均勻的影響，采用下式代替式（24):
[0061] 同理可推導出三維空間下任意函數及其一階導數和拉普拉斯算子的近似方案：
[0064]其中，Xji= X .j-Xi，yji= y .jii，Zji= Z fZi，f iiZ分別為粒子 i 在 x、y 和 z方向上的一階導數，fiiXX、4"和f iiZZ分別為粒子i在x、y和z方向上的二階導數，x k、yk 和zk分別為粒子k在x、y和z方向上的坐標，k = i或j。修正矩陣A ""在三維空間下的表 達式為：
[0066] 式（6)、式（7)和式（13)給出了一維空間下任意函數及其一階導數和拉普拉斯算 子的近似方案，式（14)、式（15)和式（25)給出了二維空間下任意函數及其一階導數和拉普 拉斯算子的近似方案，式（26)、式（27)和式（28)給出了三維空間下任意函數及其一階導數 和拉普拉斯算子的近似方案。該方案具有較高的數值精度和穩定性，便于邊界條件的處理。 [0067] 步驟4、根據所研究問題的控制方程進行時間迭代，迭代一定步數后，得到模擬結 果，其中迭代公式為：
[0069] 式中，dt為迭代的時間步長。
[0070] 有益效果
[0071] 本發明對比已有技術具有如下優點：
[0072] (1)改進的KGF-SPH方法計算量較小；
[0073] (2)改進的KGF-SPH方法提高了數值模擬的精度和穩定性；
[0074] (3)改進的KGF-SPH方法便于邊界條件的處理。
【附圖說明】
[0075] 圖1為本發明實施例1初始粒子分布示意圖；
[0076] 圖2為本發明實施例1基于（a) KGF-SPH方法、（b)改進的KGF-SPH方法以及（c) 解析解得到的一維熱傳導問題溫度隨時間變化的分布圖；
[0077] 圖3為本發明實施例1基于KGF-SPH方法和改進的KGF-SPH方法模擬一維熱傳導 問題得到的t = 0. 2s時的溫度分布曲線與解析解的對比；
[0078] 圖4為本發明實施例2初始粒子分布示意圖；
[0079] 圖5為本發明實施例2基于（a) SPH方法、（b)改進的KGF-SPH方法和（c) FVM方 法模擬封閉方腔自然對流得到的溫度分布云圖對比；
[0080] 圖6為本發明實施例2基于（a) SPH方法、（b)改進的KGF-SPH方法和（C) FVM方 法模擬封閉方腔自然對流得到的水平方向速度分布云圖對比；
[0081] 圖7為本發明實施例2基于（a) SPH方法、（b)改