一種改進(jìn)的kgf-sph方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及一種改進(jìn)的KGF-SPH方法，屬于計(jì)算流體力學(xué)技術(shù)領(lǐng)域。
【背景技術(shù)】
[0002] SPH (Smoothed Particle Hydrodynamics)是一種拉格朗日型無(wú)網(wǎng)格粒子方法，由 Lucy，Gingold與Monaghan于1977年提出并用于求解三維開放空間的天體物理學(xué)問(wèn)題，現(xiàn) 在已經(jīng)被廣泛用于多相流、重力流、熱傳導(dǎo)等流體力學(xué)問(wèn)題的數(shù)值模擬中。在SPH方法中， 用一系列粒子來(lái)代替連續(xù)介質(zhì)，這些粒子可以在計(jì)算域內(nèi)自由移動(dòng)，因此SPH方法可以很 好的處理傳統(tǒng)有網(wǎng)格方法難以解決的問(wèn)題，如具有移動(dòng)邊界、自由表面或邊界極大變形的 流動(dòng)問(wèn)題。雖然SPH方法已經(jīng)在流體力學(xué)和固體力學(xué)領(lǐng)域有了非常廣泛的應(yīng)用，但是SPH 方法在對(duì)流問(wèn)題的應(yīng)用仍較少。
[0003] 對(duì)流現(xiàn)象無(wú)處不在，如日常生活中的熱水沸騰，自然界中的海洋環(huán)流，工程應(yīng)用 中的核反應(yīng)堆冷卻等，因此研究對(duì)流問(wèn)題具有非常的意義。但是傳統(tǒng)SPH方法存在粒子 不一致性、應(yīng)力不穩(wěn)定性等缺點(diǎn)，無(wú)法很好的模擬對(duì)流問(wèn)題。在SPH方法的基礎(chǔ)上，Huang 等人在 2015 年提出了 一種 KGF-SPH 方法（kernel gradient free-smoothed particle hydrodynamics)。該方法在核近似和粒子近似過(guò)程中只用到了核函數(shù)本身，不包含核函數(shù) 的導(dǎo)數(shù)，因此無(wú)論核函數(shù)是否可導(dǎo)都可用于該方法，是一種無(wú)核梯度的SPH方法。KGF-SPH 方法降低了無(wú)網(wǎng)格方法對(duì)于核函數(shù)的要求，具有較高的精度。但是KGF-SPH方法中二階導(dǎo) 數(shù)通過(guò)對(duì)一階導(dǎo)數(shù)再求導(dǎo)來(lái)求解，這樣雖然提高了計(jì)算精度，但是存在以下缺點(diǎn)：
[0004] (1)計(jì)算量較大。因?yàn)槎A導(dǎo)數(shù)的求解需要在對(duì)一階導(dǎo)數(shù)求導(dǎo)，增加了一倍的計(jì)算 量。
[0005] (2)不利于處理邊界條件。因?yàn)樵谇蠼膺吔绺浇黧w粒子的物理量時(shí)，需要先求解 出邊界上粒子物理量的導(dǎo)數(shù)值，然而為了精確的求出邊界上粒子物理量的導(dǎo)數(shù)值，還需要 在邊界外布置更多的虛粒子。

【發(fā)明內(nèi)容】

[0006] 本發(fā)明的目的是為了提高數(shù)值模擬的精度和穩(wěn)定性，減小計(jì)算量，提出了一種改 進(jìn)的KGF-SPH方法。
[0007] 本發(fā)明是通過(guò)以下技術(shù)方案實(shí)現(xiàn)的：
[0008] -種改進(jìn)的KGF-SPH方法，具體實(shí)現(xiàn)步驟如下：
[0009] 步驟1、在所研究問(wèn)題的計(jì)算域中布置粒子；
[0010] 步驟2、根據(jù)所研究問(wèn)題的初始條件和邊界條件對(duì)粒子的物理屬性進(jìn)行初始化；
[0011] 步驟3、利用下述空間離散格式對(duì)所研究問(wèn)題的控制方程進(jìn)行近似，得到函數(shù)的時(shí) 間導(dǎo)數(shù)#的粒子近似式，其近似原理如下： dt
[0012] 對(duì)任意函數(shù)f (r])在點(diǎn)r =巧處進(jìn)行泰勒展開，并且只保留到一階導(dǎo)數(shù)，可得：
[0013] f (r j) = f (r；) +Vf (r；) · Cr jTi) +〇 (h2) (I)
[0014] 忽略式⑴中的高階項(xiàng)o (h2)，在式⑴兩邊同時(shí)乘以核函數(shù)W(R，h)和（rfrj W (R, h)并在計(jì)算域上進(jìn)行積分得：
[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] 聯(lián)立式⑵和式（3)，并求解方程組，可得：
[0019] 上式對(duì)應(yīng)的粒子近似式為：
[0021] 式中，t= f(r ;)，fj= f(r .)，（ ▽ :〇; = ▽ f(r J，N表示在i粒子的支持域內(nèi)的 j粒子（即i粒子的最近相鄰粒子）總量，P ,分別表示j粒子的質(zhì)量和密度。
[0022] 在一維空間下，式（5)的具體表達(dá)式為：
[0024] 其中，Xjl= Xj-X1Jlix為粒子i在X方向上的一階導(dǎo)數(shù)。修正矩陣厶立一維空間 下的表達(dá)式為：
[0026] 在一維空間下，對(duì)函數(shù)f (Xj)在點(diǎn)Xj= X 1進(jìn)行二階泰勒展開并略去高階項(xiàng)可得：
[0028] 式中，fliXX為粒子i在X方向上的二階導(dǎo)數(shù)。式⑶兩邊同時(shí)乘以核函數(shù)W并移 項(xiàng)可得：
[0030] 式（9)兩邊同時(shí)積分可得：
[0032] 由核函數(shù)的對(duì)稱特性可知，上式右邊第一項(xiàng)為0,因此式（10)可化簡(jiǎn)為：
[0034] 對(duì)（11)式進(jìn)行移項(xiàng)可以推導(dǎo)出一維空間下拉普拉斯算子的離散格式為：
[0036] 式（12)對(duì)應(yīng)的粒子近似式為：
[0038] 在二維空間下，式（5)的具體表達(dá)式為：
[0040] 其中，Xji= X fXi，Yji= yAjP f。分別為粒子i在X和y方向上的一階導(dǎo) 數(shù)。修正矩陣Ai在二維空間下的表達(dá)式為：
[0042] 在二維空間下，對(duì)函數(shù)f(Xj，yj)在點(diǎn)（Xl，yi)進(jìn)行二階泰勒展開并略去高階項(xiàng)可 得：
[0044] 式中，fiiXy為粒子i的混合導(dǎo)數(shù)，f iiX)^P f iiyy分別為粒子i在X和y方向上的二階 導(dǎo)數(shù)。式（16)兩邊同時(shí)乘以核函數(shù)W并移項(xiàng)可得：
[0046] 式（17)兩邊同時(shí)積分得：
[0048] 由核函數(shù)的對(duì)稱特性可知，式（18)右邊第一項(xiàng)和第二項(xiàng)為0,因此上式可化簡(jiǎn)為：
[0050] 當(dāng)粒子分布均勻且規(guī)則時(shí)將
在極坐標(biāo)下進(jìn)行求解可得：
[0055] 式（19)移項(xiàng)并利用式（22)對(duì)其進(jìn)行化簡(jiǎn)，可以推導(dǎo)出二維空間下拉普拉斯算子 的離散格式為：
[0057] 式（23)對(duì)應(yīng)的粒子近似式為：
[0059] 為了減弱粒子分布不均勻的影響，采用下式代替式（24):
[0061] 同理可推導(dǎo)出三維空間下任意函數(shù)及其一階導(dǎo)數(shù)和拉普拉斯算子的近似方案：

[0064]其中，Xji= X .j-Xi，yji= y .jii，Zji= Z fZi，f iiZ分別為粒子 i 在 x、y 和 z方向上的一階導(dǎo)數(shù)，fiiXX、4"和f iiZZ分別為粒子i在x、y和z方向上的二階導(dǎo)數(shù)，x k、yk 和zk分別為粒子k在x、y和z方向上的坐標(biāo)，k = i或j。修正矩陣A ""在三維空間下的表 達(dá)式為：
[0066] 式（6)、式（7)和式（13)給出了一維空間下任意函數(shù)及其一階導(dǎo)數(shù)和拉普拉斯算 子的近似方案，式（14)、式（15)和式（25)給出了二維空間下任意函數(shù)及其一階導(dǎo)數(shù)和拉普 拉斯算子的近似方案，式（26)、式（27)和式（28)給出了三維空間下任意函數(shù)及其一階導(dǎo)數(shù) 和拉普拉斯算子的近似方案。該方案具有較高的數(shù)值精度和穩(wěn)定性，便于邊界條件的處理。 [0067] 步驟4、根據(jù)所研究問(wèn)題的控制方程進(jìn)行時(shí)間迭代，迭代一定步數(shù)后，得到模擬結(jié) 果，其中迭代公式為：
[0069] 式中，dt為迭代的時(shí)間步長(zhǎng)。
[0070] 有益效果
[0071] 本發(fā)明對(duì)比已有技術(shù)具有如下優(yōu)點(diǎn)：
[0072] (1)改進(jìn)的KGF-SPH方法計(jì)算量較小；
[0073] (2)改進(jìn)的KGF-SPH方法提高了數(shù)值模擬的精度和穩(wěn)定性；
[0074] (3)改進(jìn)的KGF-SPH方法便于邊界條件的處理。
【附圖說(shuō)明】
[0075] 圖1為本發(fā)明實(shí)施例1初始粒子分布示意圖；
[0076] 圖2為本發(fā)明實(shí)施例1基于（a) KGF-SPH方法、（b)改進(jìn)的KGF-SPH方法以及（c) 解析解得到的一維熱傳導(dǎo)問(wèn)題溫度隨時(shí)間變化的分布圖；
[0077] 圖3為本發(fā)明實(shí)施例1基于KGF-SPH方法和改進(jìn)的KGF-SPH方法模擬一維熱傳導(dǎo) 問(wèn)題得到的t = 0. 2s時(shí)的溫度分布曲線與解析解的對(duì)比；
[0078] 圖4為本發(fā)明實(shí)施例2初始粒子分布示意圖；
[0079] 圖5為本發(fā)明實(shí)施例2基于（a) SPH方法、（b)改進(jìn)的KGF-SPH方法和（c) FVM方 法模擬封閉方腔自然對(duì)流得到的溫度分布云圖對(duì)比；
[0080] 圖6為本發(fā)明實(shí)施例2基于（a) SPH方法、（b)改進(jìn)的KGF-SPH方法和（C) FVM方 法模擬封閉方腔自然對(duì)流得到的水平方向速度分布云圖對(duì)比；
[0081] 圖7為本發(fā)明實(shí)施例2基于（a) SPH方法、（b)改