SPH法理論基礎(chǔ)參考模板

1、摘要1第一章 緒論21.1 基于網(wǎng)格的方法21.1.1 拉格朗日網(wǎng)格31.1.2 歐拉網(wǎng)格31.1.3 基于網(wǎng)格的數(shù)值方法的局限性31.2 無網(wǎng)格方法41.3光滑粒子動力學(xué)方法4第二章 SPH方法在流體動力學(xué)問題中的應(yīng)用52.1 光滑粒子動力學(xué)原理52.2 光滑粒子動力學(xué)基本方程52.2.1 函數(shù)的積分表達52.2.2 函數(shù)的粒子表達72.2.3光滑函數(shù)72.3 拉格朗日型的Navier-Stokes 方程72.4 Navier-Stokes方程的SPH表達式8第三章 SPH的計算實施方法93.1 粒子的密度近似93.2 核函數(shù)103.3 狀態(tài)方程113.4 人工粘度123.5 邊界處理133

2、.6 時間積分131 / 16摘要計算機數(shù)值仿真逐漸成為解決現(xiàn)代工程和科學(xué)問題的一條重要途徑。數(shù)值仿真能為理論提供測試和檢驗，有助于對復(fù)雜的物理問題加深認識，甚至還能幫助解釋和發(fā)現(xiàn)新現(xiàn)象?；诰W(wǎng)格的數(shù)值方法雖然已經(jīng)有廣泛的應(yīng)用，但是在很多方面仍存在不足之處，比如在計算流體動力學(xué)的大變形、運動物質(zhì)交界面、自由表面等問題時，由于網(wǎng)格產(chǎn)生畸變導(dǎo)致計算誤差過大或無法進行，從而使其在許多問題的應(yīng)用上受到限制。近年來，無網(wǎng)格法倍受關(guān)注，這種方法在許多應(yīng)用中都優(yōu)于傳統(tǒng)的基于網(wǎng)格的有限元法、有限差分法以及有限體積法等數(shù)值方法。本文主要研究新一代無網(wǎng)格方法-光滑粒子流體動力學(xué)方法（SPH）。第一章 緒論1.1

3、基于網(wǎng)格的方法數(shù)值計算方法通常可以分為兩種：基于網(wǎng)格的方法和無網(wǎng)格方法。通常對于物理控制方程的描述有兩種基本方法：歐拉描述法和拉格朗日描述法。歐拉描述法是對空間的描述方法，其典型代表是有限差分法；拉格朗日描述法是對物質(zhì)點的描述方法，其典型代表是有限元法。歐拉描述和拉格朗日描述對應(yīng)著兩種不同的區(qū)域離散化網(wǎng)格：歐拉網(wǎng)格和拉格朗日網(wǎng)格。針對不同類型的問題，這兩種網(wǎng)格在數(shù)值方法中都得到廣泛應(yīng)用。1.1.1 拉格朗日網(wǎng)格拉格朗日網(wǎng)格，基于拉格朗日網(wǎng)格的數(shù)值方法在整個計算過程中網(wǎng)格是固定附著于物質(zhì)上的，網(wǎng)格會隨著物質(zhì)的運動而運動，所以在物質(zhì)點上的所有場變量的整個時間歷程都可以很容易地追蹤，常見的方法如有限

4、元法等?；诶窭嗜站W(wǎng)格方法的優(yōu)點是：由于在相關(guān)的偏微分方程里不存在遷移項，所以程序在方案設(shè)計上會變得相對簡單而且運行較快；由于只需在問題域內(nèi)布置網(wǎng)格，問題域外不需要布置，所以計算效率很高；不規(guī)則或者復(fù)雜的幾何形狀可以用不規(guī)則的網(wǎng)格來處理。由于具有以上這些優(yōu)點，拉格朗日方法得到廣泛的應(yīng)用，并且能成功地求解計算固體力學(xué)問題。然后，基于拉格朗日網(wǎng)格的方法難以應(yīng)用于具有極大網(wǎng)格變形的情況，因為其公式的形式是以網(wǎng)格為基礎(chǔ)的，當(dāng)網(wǎng)格變形太大的時候，公式的精度和求解都會受到很大的影響。另外，由于時間步長是由最小單元尺寸所控制，若網(wǎng)格太小就會影響計算的效率，甚至?xí)?dǎo)致計算失敗。1.1.2 歐拉網(wǎng)格歐拉網(wǎng)格，

5、相對于拉格朗日網(wǎng)格，歐拉網(wǎng)格剛好相反，它是固定在模擬對象所處的空間上，模擬對象在固定網(wǎng)格單元上運動。因此，在物質(zhì)流過網(wǎng)格時，所有網(wǎng)格節(jié)點和網(wǎng)格單元依然固定在空間上而且不會隨時間的改變而改變。通過模擬質(zhì)量、動量和能量經(jīng)過網(wǎng)格單元邊界的通量，可以計算質(zhì)量、速度和能量等等的分布。在整個計算過程中網(wǎng)格單元的形狀和體積都保持不變。由于歐拉網(wǎng)格在時間和空間上都是固定的，物體的大變形不會引起網(wǎng)格本身的任何變化，所以再以物質(zhì)流動為主體的計算流體力學(xué)領(lǐng)域里，較為廣泛的應(yīng)用歐拉法。但是，歐拉法依然有很多缺點。由于歐拉法不能用固定網(wǎng)格來追蹤物質(zhì)的運動，所以很難分析物質(zhì)上固定點的場變量的變化情況而只能得到空間上固定的

6、歐拉網(wǎng)格的場變量的變化情況；由于歐拉法追蹤的是流過網(wǎng)格單元邊界的質(zhì)量、動量和能量的通量，所以自由表面、變形邊界和運動物質(zhì)交界面的位置就很難精確確定；由于歐拉法需要在整個計算區(qū)域上都覆蓋網(wǎng)格，所以有時為了提高計算效率而使用較粗糙的網(wǎng)格，這樣就會降低離散化區(qū)域的求解精度。1.1.3 基于網(wǎng)格的數(shù)值方法的局限性傳統(tǒng)的基于網(wǎng)格的數(shù)值方法如有限差分法和有限元法已經(jīng)廣泛地應(yīng)用在計算流體力學(xué)和計算固體力學(xué)的各個領(lǐng)域中，是現(xiàn)在進行區(qū)域離散化和數(shù)值離散化模擬的主要方法。但是基于網(wǎng)格的數(shù)值方法在很多方面仍存在不足之處，比如在計算流體動力學(xué)中的大變形、運動物質(zhì)交界面、自由表面等問題時，由于網(wǎng)格產(chǎn)生畸變導(dǎo)致計算誤差過

7、大或無法進行，使其在許多問題的應(yīng)用上受到限制。在基于網(wǎng)格的數(shù)值方法中，數(shù)值模擬的先決條件就是在問題域生成網(wǎng)格。對于基于歐拉網(wǎng)格的方法，在固定的歐拉網(wǎng)格上要精確確定自由表面、變形邊界、運動交界面和不均勻物質(zhì)之間的位置是一項非常困難的工作。并且歐拉方法也不適用于研究如粒子流動這類問題。即需要在固定的物質(zhì)體內(nèi)監(jiān)控材料特性的問題。對于拉格朗日網(wǎng)格法如有限元法，進行模擬前必須要在研究對象上建立網(wǎng)格，這項操作常常占用很大的計算工作量。當(dāng)所研究對象是一系列離散的物質(zhì)點時，同樣不適合使用基于網(wǎng)格的方法，用連續(xù)的基于網(wǎng)格的方法對這些離散系統(tǒng)進行模擬通常不是好的選擇。在計算大變形問題時，比如高速碰撞、金屬加工成型

8、、動態(tài)裂紋擴展、流固耦合和應(yīng)變局部化等，基于網(wǎng)格的方法遇到的困難更大。原因是發(fā)生大變形時，網(wǎng)格畸變過大，使得計算中斷，有限元方法通常采用單元“銷蝕”法或重分網(wǎng)格技術(shù)來克服這種困難。但是單元“銷蝕”法本身缺乏物理依據(jù)，純粹是為了使計算進行下去的一種數(shù)值手段。而網(wǎng)格重分技術(shù)在節(jié)點重新分配物理量時，很難保證系統(tǒng)動量、能量守恒，因而導(dǎo)致計算的精度下降。此外，網(wǎng)格重分技術(shù)不是很容易實現(xiàn)，為了更好地解決大變形問題，必須對網(wǎng)格有新的處理方法或去除網(wǎng)格，所以各種無網(wǎng)格方法相繼被提出來。1.2 無網(wǎng)格方法無網(wǎng)格方法產(chǎn)生于三十多年以前，當(dāng)時發(fā)展很慢，直到近幾年，由于無網(wǎng)格方法的近似函數(shù)不依賴于網(wǎng)格，在分析涉及大變

9、形的問題中被認為優(yōu)于傳統(tǒng)的基于網(wǎng)格的有限差分法和有限元法，才受到了眾多研究者的高度重視，成為國際計算仿真界的研究熱點之一，并得到了迅速的發(fā)展【1】。無網(wǎng)格方法的主要思想是：通過使用一系列任意分布的節(jié)點（或粒子）來求解具有各種邊界條件的積分方程或偏微分方程組，從而得到精確穩(wěn)定的數(shù)值解，這些節(jié)點或粒子之間不需要網(wǎng)格進行連接。由于無網(wǎng)格法采用基于點的近似，可以徹底或部分地消除網(wǎng)格，不需要網(wǎng)格的初始劃分和重構(gòu)，因此不僅可以保證計算的精度，而且可以減小計算的難度。經(jīng)過幾十年的研究，目前已發(fā)展了十余種無網(wǎng)格方法。最早提出的無網(wǎng)格方法是光滑粒子動力學(xué)方法（Smoothed Particle Hydrodyn

10、amics，簡稱SPH），它是由Lucy、Gingold和Monaghan在1977年分別提出的，并且在天體領(lǐng)域得到成功地應(yīng)用【2,3】。但是由于精度和穩(wěn)定性問題，該方法最初并未得到廣泛的應(yīng)用。80年代，Monaghan等人在該方法的研究與應(yīng)用中作出突出貢獻【4,5】，將光滑粒子動力學(xué)方法應(yīng)用到連續(xù)固體力學(xué)和流體力學(xué)中，他們將SPH方法解釋為核函數(shù)法，模擬了流場中的激波強間斷現(xiàn)象。90年代，Swegle、Dyka和Chen等提出了SPH方法不穩(wěn)定的原因及穩(wěn)定化方法【6-8】。Johnson和Beissel等人也提出了一些用來改善應(yīng)變計算的方法【9】。隨著研究的深入，SPH方法被應(yīng)用于水下爆炸數(shù)

11、值仿真【10】、高速碰撞中材料動態(tài)響應(yīng)數(shù)值仿真等領(lǐng)域【11-13】。近幾年來，我國的學(xué)者也開始關(guān)注SPH計算方法，中科院的張鎖春對SPH方法進行了綜述【14】，國防科學(xué)技術(shù)大學(xué)的鄧方剛研究了SPH在柱坐標(biāo)體系下的應(yīng)用【15】。國防科大的貝新源、岳宗五等將SPH方法用于高速碰撞問題研究【16】。1.3光滑粒子動力學(xué)方法光滑粒子流體動力學(xué)（SPH-Smoothed particle hydrodynamics）方法是近30年發(fā)展起來的一種新的純拉格朗日無網(wǎng)格粒子法，它最初是用于解決三維開放空間的天體物理學(xué)問題【2,3】，現(xiàn)已被廣泛地研究和擴展，并被應(yīng)用于具有材料強度的動態(tài)響應(yīng)問題和具有大變形的流體

12、動力學(xué)問題。SPH方法和它的派生方法是現(xiàn)在粒子法的主要類型。SPH方法通過對鄰近粒子進行加權(quán)平均而得到穩(wěn)定的光滑近似性質(zhì)，該方法應(yīng)用在流體動力學(xué)問題的范圍內(nèi)。由于SPH方法中的自適應(yīng)性、粒子性質(zhì)和拉格朗日性質(zhì)的和諧結(jié)合，使其在工程和科學(xué)領(lǐng)域都得到了實際應(yīng)用。國外對SPH方法研究的比較早，已經(jīng)被應(yīng)用到很多領(lǐng)域中。SPH方法最初在天體研究中應(yīng)用較多，如雙子星和行星碰撞、超新星、星系的形成和瓦解、黑洞和中子星合并、白矮星的單個或多重爆炸甚至是宇宙的進化等問題的研究。第二章 SPH方法在流體動力學(xué)問題中的應(yīng)用2.1 光滑粒子動力學(xué)原理流體動力學(xué)問題的求解主要是解基于密度、速度、能量等變量場的偏微分方程

13、組。除了一些簡單問題可以得到偏微分方程組的解析解以外，大部分的流體動力學(xué)問題只能尋求其數(shù)值解。光滑粒子動力學(xué)的原理如下：1） 使用SPH方法求解流體動力學(xué)問題時，問題域用一系列任意分布的粒子來表示，粒子之間沒有任何連接，因此SPH方法是無網(wǎng)格性質(zhì)的，通常為了求解方便，初始粒子都均勻排列；2） 場函數(shù)用積分表示法來近似，在SPH方法中稱為核近似法；3） 然后應(yīng)用支持域內(nèi)的相鄰粒子對應(yīng)的值疊加求和取代場函數(shù)的積分表達式來對場函數(shù)進行粒子近似，由于在每一個時間步內(nèi)都要進行粒子近似，支持域內(nèi)的有效粒子為當(dāng)前時刻支持域內(nèi)的粒子，因此SPH方法具有自適應(yīng)性；4） 將粒子近似法應(yīng)用于所有偏微分方程組的場函數(shù)

14、相關(guān)項中，將偏微分方程組進行離散，可知SPH是一種純拉格朗日方法；5） 粒子被附上質(zhì)量后，則意味著這些粒子是真實的具有材料特性的粒子；最后應(yīng)用顯式積分法得到所有粒子的場變量隨時間的變化值。從以上分析可以看到，SPH方法是一種純拉格朗日的具有無網(wǎng)格、自適應(yīng)屬性的流體動力學(xué)求解方法。2.2 光滑粒子動力學(xué)基本方程使用SPH方法求解流體動力學(xué)問題一般分兩步進行。第一步是積分表示，即使用積分表達式對函數(shù)進行近似，函數(shù)的積分表達式可以通過對核函數(shù)的影響區(qū)域積分進行近似。第二步是粒子表示，即通過對最近相鄰粒子的值進行累加求和來近似離散點的函數(shù)值。2.2.1 函數(shù)的積分表達標(biāo)準SPH是基于核估計發(fā)展起來的，

15、在SPH方法中，對任意一個函數(shù)f，其積分近似表達式為其中，<f(x)>為函數(shù)f(x)的近似值，x為位置矢量，為包含x的幾分體，W(x-x,h)為光滑函數(shù)，其依賴于兩點之間的距離|x-x和光滑長度h。光滑函數(shù)在SPH近似法中其著重要的作用，它決定了函數(shù)表達式的精度和計算效率。光滑函數(shù)W常常選用偶函數(shù)，其應(yīng)滿足3個條件。第一個條件為正則化條件，如下所示由于光滑函數(shù)的積分值等于1，故此條件也稱為歸一化條件。第二個條件是當(dāng)光滑長度趨向于零時具有狄拉克函數(shù)性質(zhì)這里，第三個條件是緊支性條件式中，K是與點x處光滑函數(shù)相關(guān)的常數(shù)，并確定光滑函數(shù)的有效范圍（非零）。此有效范圍稱作點x處光滑函數(shù)的支持

16、域內(nèi)的積分。因此，一般來說積分域就是支持域。根據(jù)分布積分和散度定理，將第一個式子經(jīng)過變換，對f（x）的導(dǎo)數(shù)的估計為由以上方程得知，SPH估計將函數(shù)的空間導(dǎo)數(shù)轉(zhuǎn)化為光滑函數(shù)的空間導(dǎo)數(shù)，因此，可以根據(jù)核函數(shù)求取任意場函數(shù)的空間導(dǎo)數(shù)。2.2.2 函數(shù)的粒子表達SPH方法最終要將函數(shù)積分近似表達式轉(zhuǎn)化為支持域內(nèi)所有粒子疊加求和的離散化形式，流體的問題域被離散化為有限數(shù)量的粒子，其中每個粒子具有獨立的質(zhì)量、密度及其它物理屬性。經(jīng)過離散化，函數(shù)的積分表達式可以寫成如下的粒子近似式其中，j和i表示粒子編號，mj和j分別j粒子的質(zhì)量和密度，N是粒子的總數(shù)，Wij=W(xi-xj，h)=W(|xi-xj |，h

17、)。上式說明，在SPH方法中，對任意一個函數(shù)f，它在某一位置處的函數(shù)值可以通過應(yīng)用光滑函數(shù)W對其光滑長度為h的緊支域內(nèi)所有粒子插值求和的形式來表達【17】。光滑長度h在SPH方法中非常重要，它在計算過程中決定了計算的實用性、有效性和可靠的適應(yīng)性。光滑長度h直接影響到計算結(jié)果的效率和結(jié)果的精度。若h太小，則在以kh為半徑的支持域內(nèi)沒有足夠的粒子對所求粒子施加影響，這樣會導(dǎo)致計算結(jié)果精度較低。如果h太大，太多的粒子對所求粒子產(chǎn)生影響，則可能將粒子所有細節(jié)信息和局部特性忽略，導(dǎo)致粒子趨于雷同，同樣會影響到結(jié)果的精度。具體光滑長度應(yīng)該選擇多大與所研究的問題相關(guān)，沒有一個廣義的光滑長度。2.2.3光滑函

18、數(shù)SPH通過光滑函數(shù)來進行積分表示。光滑函數(shù)非常重要，因為它不僅決定函數(shù)表達式的形式，定義粒子支持域的大小，而且還決定核近似和粒子近似的一致性和計算效率。常用的光滑函數(shù)有三種，分別是高斯型光滑函數(shù)，三次樣條光滑函數(shù)和高次樣條函數(shù)?？蓞⒁妳⒖嘉墨I【3,18,19,】2.3 拉格朗日型的Navier-Stokes 方程對控制方程的描述有兩種形式：歐拉描述法和拉格朗日描述法。歐拉描述法是對空間點的描述，而拉格朗日描述法是對物質(zhì)點的描述。流體的流動過程可以用歐拉方程以及適當(dāng)?shù)臓顟B(tài)方程來表示，這些方程表示了質(zhì)量、動量和能量守恒。如果用上標(biāo)和來表示坐標(biāo)方向和求和重復(fù)指標(biāo)，采用時間全導(dǎo)數(shù)形式，N-S方程組可

19、以表示為：連續(xù)性方程：動量方程：能量方程：在以上方程中，為總應(yīng)力張量，有兩部分組成：一部分是壓力p，另外一部分是粘性應(yīng)力。對于牛頓流體，粘性剪應(yīng)力與剪應(yīng)變率成比例，且比例系數(shù)為粘性系數(shù)。其中如果將壓力與粘性應(yīng)力分開寫，我們可以得到如下的能量方程在上方式子中，各字母代表的意義如下：V -速度向量E-內(nèi)能-密度p-壓力t-時間2.4 Navier-Stokes方程的SPH表達式根據(jù)光滑函數(shù)的概念和SPH方法的粒子近似原理，流體流動的運動方程可以寫成如下的離散SPH形式：其中 i-參考粒子。 J-粒子i的相鄰粒子。第三章 SPH的計算實施方法3.1 粒子的密度近似由于在SPH方法中粒子的運動是由相鄰

20、粒子之間的密度差以及壓力差產(chǎn)生的，所以密度近似在SPH方法中非常重要。在一般的SPH方法中主要有兩種方式對密度進行展開：一種是密度求和法，另一種是連續(xù)密度法。在SPH方法中粒子的運動是由相鄰粒子之間的密度差以及壓力差產(chǎn)生的，所以密度近似在SPH方法中非常重要。在一般的SPH方法中主要由兩種方式對密度進行展開：一種是密度求和法，另一種是連續(xù)密度法。在SPH方法的實際應(yīng)用中，密度求和法似乎使用得更為廣泛，主要是因為密度求和法體現(xiàn)了SPH近似法的本質(zhì)?，F(xiàn)在已經(jīng)有不少改進此種算法精度的修正方案。其中一種有效地改進是將鄰近粒子的SPH光滑函數(shù)累加式，即上式的等號右端正則化。上式可提高自有邊界處和密度不連

21、續(xù)材料交界面處的精度，此式也非常適合模擬非連續(xù)廣義流體的流動問題。密度求和法的缺陷是所需要的計算量大，因為必須在計算其他參數(shù)前計算密度，而且還要計算自身的光滑函數(shù)。近似密度的另一種粒子近似法為連續(xù)性密度法，這種近似法通過應(yīng)用SPH近似法的概念對連續(xù)性方程進行轉(zhuǎn)換而得。常用的連續(xù)性密度方程為：連續(xù)性密度近似法不需要先計算密度，因此可以提高計算效率，并適于并行處理算法的編程。具體應(yīng)用哪一種密度算法要根據(jù)研究的問題而定，對于廣義流體問題的模擬，應(yīng)用修正的密度求和法可得到較好的結(jié)果。對于具有強間斷問題的模擬，如爆炸、高速沖擊等，應(yīng)優(yōu)先選取連續(xù)性密度法。3.2 核函數(shù)核函數(shù)在SPH近似法中起著重要作用，

22、因為它控制了插值的方式和粒子影響域的效能，從而決定了函數(shù)表達式的精度和計算效率。由于SPH插值基于積分插值理論，用核函數(shù)來近似狄拉克函數(shù)，因此核函數(shù)應(yīng)該具有某些特性，如非負性，緊支性，歸一性，衰減性和狄拉克函數(shù)條件。在實際應(yīng)用中，通常使用兩種形式的核函數(shù)，一種是如下所示的高斯核函數(shù)。公式中的d在一維空間中為，在二維空間中為，在三維空間中為。另一種是三次樣條核函數(shù)，到目前為止，三次樣條函數(shù)是最廣泛應(yīng)用的核函數(shù)。在一維空間中，在二維空間中，在本文的三維空間中，。3.3 狀態(tài)方程在求解可壓縮流體問題的標(biāo)準SPH方法中，粒子的運動是由于壓力梯度的作用而產(chǎn)生的，而粒子的壓力是通過狀態(tài)方程由粒子自身的密度

23、和內(nèi)能來計算的。然而，在不可壓縮流體問題中，流體實際的狀態(tài)方程限制了時間步長的大小，即時間步長不能太小。有效地計算動量方程中的壓力項是仿真計算不可壓縮流體的一個主要任務(wù)。事實上，理論不可壓縮流體實際上是可壓縮的，所以可以引進人工壓縮率。人工壓縮率的引進就是把所有理論不可壓縮流都考慮為實際上是可壓縮的。因此，可用仿真可壓縮狀態(tài)方程去仿真不可壓縮流。不同的物質(zhì)應(yīng)該具有不同的狀態(tài)方程。在本文中應(yīng)用較多的是如下所示的狀態(tài)方程:式中：是常數(shù)，在許多情況下取=7；0是參照密度；B是由具體問題而定的參數(shù)，用于限制密度的最大改變量。在許多情況下，一般用B作為初始壓力。上式中減去1能消除自由表面流動的邊緣效應(yīng)。

24、由上式可以看出，密度的微小變化都會導(dǎo)致壓力值產(chǎn)生較大的波動。另外一個可選擇的人工狀態(tài)方程為式中：c為聲速。1997年，Morris使用此人工狀態(tài)方程通過SPH方法數(shù)值仿真了低雷諾數(shù)不可壓縮流動【20】。Zhu等使用此人工狀態(tài)方程通過SPH方法用于微孔尺度的數(shù)值模型，用于數(shù)值仿真穿越多孔介質(zhì)的流動。在人工壓縮率技術(shù)中，聲速是一個需要慎重考慮的因素。如果使用實際的聲速，則意味著用理想不可壓縮的人工流體來近似真實的流體。相對密度差為：式中：Vb和M分別是流體的整體速度和馬赫數(shù)。由于真實的聲速相當(dāng)大，所以相應(yīng)得到的馬赫數(shù)非常小，相對密度差幾乎可以忽略。因此為了用人工壓縮流體近似真實流體，應(yīng)該使用比真實

25、值要小很多的聲速。所以對聲速大小的要求為：一方面，必須足夠大，以至于人工可壓縮流體的特性與真實流體充分接近；另一方面，應(yīng)足夠小，可使時間步的增量在容許范圍內(nèi)。3.4 人工粘度對于某些數(shù)值問題，誤差跳躍很大，會產(chǎn)生數(shù)值振蕩，為了提高算法的穩(wěn)定性，并且防止粒子間相互接近時的非物理穿透，Monaghan和Gingold在SPH方法的動量方程中引入了人工粘度【21】，它不是物理的粘度，只是數(shù)值的修正。式中，在式中，為常數(shù)，取值一般為1.0左右【22】。因子=0.1hij用于防止粒子相互靠近時產(chǎn)生的數(shù)值發(fā)散。c和v分別表示聲速和粒子的速度矢量。將人工粘度引入到動量方程中，最終動量方程的SPH粒子近似式可

26、寫為如下形式：3.5 邊界處理在SPH方法的應(yīng)用中，邊界條件的處理既是該方法的優(yōu)點，也是目前的薄弱環(huán)節(jié)，因為SPH方法最初是應(yīng)用在沒有邊界條件限制的領(lǐng)域內(nèi)，后來逐漸擴展到一些具有邊界的工業(yè)應(yīng)用中。我們使用虛粒子來處理固定邊界條件，通過設(shè)置邊界虛粒子可以方便地模擬固體邊界，同時也要防止粒子的非物理穿透。根據(jù)不同的數(shù)值仿真條件，我們在SPH數(shù)值模擬中一共使用了兩種類型的虛粒子第一種類型的虛粒子設(shè)置在固定邊界上，與Monaghan所使用的相似【4】。二種類型的虛粒子分布在邊界的外部，與Libersky所使用的相似【23】。第二種類型的虛粒子按以下的方式構(gòu)造，即給定一個實粒子i，則在邊界外與實粒子對稱

27、處分布一個虛粒子，這些虛粒子具有與相對應(yīng)實粒子相同的壓力和密度，但速度方向相反。第二種類型的虛粒子通常在邊界條件不斷變化的場合下使用【24】。3.6 時間積分可以使用一些標(biāo)準的計算方法如跳蛙法(LF)、預(yù)估校正法、龍格一庫塔法(RK)對普通微分方程進行數(shù)值積分來求解每一個粒子的物理變化。我們使用跳蛙法進行時間積分，跳蛙法的優(yōu)點是計算時所需要的存儲量低，而且在每一次計算中只需要進行一次優(yōu)化估值。粒子的密度、速度、內(nèi)能和位移按照下面的公式遞推: 應(yīng)該注意的是，跳蛙法是條件穩(wěn)定的，穩(wěn)定的條件為CFL(Courant一Friedrichs一Levy)條件，此條件一般會導(dǎo)致時間步長與光滑長度互成比例。在

28、本文中，時間步長按如下方式取值:式中，為Courant系數(shù)，一般取值在0.3左右【4】。【1】 宋康祖 緊支函數(shù)無網(wǎng)格方法研究【D】。清華大學(xué)博士學(xué)位論文。2000.【2】 Lucy L B.Numerical approach to testing the fission hypothesisj,Astronomical Journal.1977(82):1013-1024【3】 Gingold R A,Monaghan J J.Smoothed Particle Hydrodynamics:Theory and Application to Non-spherical starsC.Mon

29、thly Notices of the Royal Astronomical Society.1977(181)375-389.【4】 Monaghan J J.An introduction to SPHj.Comput.Phys.Commun.1988,48:89-96.【5】 Monaghan J J.Why particle methods workj. SIAM j.Sci.Stat.Comput.1982,3(4);422-433【6】 Johnson G R,Beissel S R.Normalized smoothing functions for SPH impact com

30、putationsJ.Int J Num Meth Eng,1996,39:2725-2741.【7】 Swegle J W,Hicks D L, Attaways S W. Smoothed particle hydrodynamics stability analysisJ.J Comput Phys,1995,116:123-134.【8】 Dyka C T. Addressing tension instability in SPH methodsR.technical Report Nrl/MR/6384,NRL,1994.【9】 Johnson G R,stryk R A. Bei

31、ssel S R.SPH for high velocity impact computationsJ. Comput.Methods Appl.Mech.Engrg.1996,139:347-373.【10】 Swegle J W, Attaway S W. On the feasibilility of using smoothed particle hydrodynamics for underwater explosion calculationsj.Comput Mtch,1995,17:151-168.【11】 Johnson G R,stryk R A,Beissel S R.S

32、PH for high velocity impact computationsJ.Computer Methods in Applied Mechanics and Engineering,1996,139:347-373【12】 Libersky L D,petschek A G,et al. High strain lagrangian hydrodynamics:A three-dimensional SPH code for dynamic material responseJ.J Comput Phys,1993.109:67-75【13】 Johnson G R,Beissel S R.Normalized smoothing functions for SPH impact computationsJ.International Journal for Numerical Methods in Engineering,1996,39:2725-2741.【14】 張鎖春.光滑質(zhì)點流體動力學(xué)（SPH）方法（綜述）J.計算物理，1996,13（4）：385-397.【1