2.2　光滑粒子動力學基本方程

使用SPH方法求解流體動力學問題一般分兩步進行。第一步是積分表示,即使用積分表達式對函數進行近似,函數的積分表達式可以通過對核函數的影響區域積分進行近似。第二步是粒子表示,即通過對最近相鄰粒子的值進行累加求和來近似離散點的函數值^[97]。

2.2.1　函數的積分表達

標準SPH是基于核估計發展起來的,在SPH方法中,對任意一個函數f,其積分近似表達式為

(2?1)

式中,<f(x)>為函數f(x)的近似值,x為位置矢量,Ω為包含x的積分體,W(x-x',h)為光滑函數,其依賴于兩點之間的距離|x-x'|和光滑長度h。

光滑函數在SPH近似法中起著重要的作用,它決定了函數表達式的精度和計算效率。光滑函數W常常選用偶函數,其應滿足3個條件。

第一個條件為正則化條件

(2?2)

由于光滑函數的積分值等于1,故此條件也稱為歸一化條件。

第二個條件是當光滑長度趨向于零時具有狄拉克函數性質

W(x-x',h)=δ(x-x')(2?3)

這里

δ(x-x')=(2?4)

第三個條件是緊支性條件

W(x-x',h)=0,當|x-x'|>κh?時(2?5)

式中,κ是與點x處光滑函數相關的常數,并確定光滑函數的有效范圍(非零)。此有效范圍稱作點x處光滑函數的支持域。應用緊支性條件,在整個問題域內的積分即轉換為在光滑函數的支持域內的積分。因此,一般來說積分域Ω就是支持域。

根據分部積分和散度定理,式(2?1)經過變換,對f(x)的導數的估計為

(2?6)

由上述方程得知,SPH估計將函數的空間導數轉化為光滑函數的空間導數,因此,可以根據核函數求取任意場函數的空間導數。

2.2.2　函數的粒子表達

SPH方法最終要將函數積分近似表達式轉化為支持域內所有粒子疊加求和的離散化形式,流體的問題域被離散為有限數量的粒子,其中每個粒子具有獨立的質量、密度及其他物理屬性。經過離散化,函數的積分表達式(2?1)可以寫成如下的粒子近似式

(2?7)

式中,j和i表示粒子編號,m_j和ρ_j分別表示j粒子的質量和密度,N是粒子的總數,W_ij=W(x_i-x_j,h)=W(|x_i-x_j|,h)。

式(2?7)說明,在SPH方法中,對任意一個函數f,它在某一位置處的函數值可以通過應用光滑函數W對其光滑長度為h的緊支域內所有粒子插值求和的形式來表達^[98],如圖2?1所示。

光滑長度h在SPH方法中非常重要,它在計算過程中決定了計算的實用性、有效性和可靠的適應性。光滑長度h直接影響到計算結果的效率和結果的精度。若h太小,則在以κh為半徑的支持域內沒有足夠的粒子對所求粒子施加影響,這樣會導致計算結果精度較低。如果h太大,太多的粒子對所求粒子產生影響,則可能將粒子所有細節信息和局部特性忽略,導致粒子趨于雷同,同樣會影響到結果的精度。具體光滑長度應該選擇多大與所研究的問題相關,沒有一個廣義的光滑長度。

式(2?6)經過粒子離散后,轉化為

(2?8)

式中,Δ_iW_ij=?=??,r_ij=|x_i-x_j|。由于W存在緊支域,因此方程中的求和實際上只對i粒子影響區域內的粒子進行,而不是在整個計算域上,如圖2?1所示。

2.2.3　光滑函數

無網格方法中的一個關鍵問題就是,如何在一系列不需要預先定義網格來提供節點連接的任意分布的離散點上有效地進行函數近似。函數近似的方法主要有三類:積分表示法、級數表示法和微分表示法。SPH方法通過光滑函數來進行積分表示。光滑函數非常重要,因為它不僅決定函數表達式的形式、定義粒子支持域的大小,而且還決定核近似和粒子近似的一致性、精度和計算效率。許多不同的光滑函數已被應用于SPH方法,許多不同的文獻討論了光滑函數的不同條件和性質。常用的光滑函數有三種,分別是高斯型光滑函數、三次樣條光滑函數和高次樣條函數。

2.2.3.1　高斯型光滑函數

Gingold和Monaghan最早使用了以下高斯型光滑函數來仿真非球形星體^[4],如圖2?2所示。

W(R,h)=α_d(2?9)

式中,α_d在一維空間中為,在二維空間中為,在三維空間中為。

高斯型光滑函數是充分光滑的,即使對于高階導數,也是一種較好的選擇,這是因為它很穩定并且精度很高,特別是對于不規則粒子分布的情況。但是,它并不是真正嚴密的,因為理論上它是不可能為零的,除非R趨向于無窮大。但是由于它數值上趨于零的速度是很快的,所以實際上是緊密的。高斯型光滑函數計算量要求很大,因為光滑函數趨于零時要經歷的距離比較大,特別是對于光滑函數的高階導數。這樣就會產生一個很大的支持區域,區域中用于粒子近似的粒子也會很多,導致運算時間較長。

Morris對高斯型光滑函數進行了修正^[99],修正的高斯光滑函數比原始的高斯型光滑函數更穩定、具有更高的效率,以下是二維的修正高斯光滑函數的表達式。

(2?10)

式中,δ為截斷界限,在二維空間里面一般取δ=3h。

2.2.3.2　三次樣條光滑函數

三次樣條光滑函數又稱B?樣條光滑函數,由Monaghan和Lattanzio在三次樣條函數的基礎上首次使用^[100],如圖2?3所示。

W(R,h)=α_d×(2?11)

在一維空間中α_d=,在二維空間中α_d=,在三維空間中α_d=。到目前為止,三次樣條光滑函數是應用最廣泛的光滑函數,這是因為在狹窄的緊支域中,三次樣條光滑函數與高斯型光滑函數類似。但是,三次樣條光滑函數的二階導數是分段線性函數,相應地,其穩定性就會比那些光滑的核函數要差一些。另外,由于其光滑函數是分段的,所以其使用就會比只有一段的光滑函數要困難一些。

2.2.3.3　高次樣條函數

Morris提出了更加近似于高斯型光滑函數并且更加穩定的高次樣條函數(四次和五次)^[101],其中五次樣條光滑函數如下所示

W(R,h)=α_d×(2?12)

式中,α_d的值在一維空間中為,在二維空間中為,在三維空間中為。