SPH法理論基礎

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

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

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

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

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

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

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

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

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

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

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

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

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

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

15、函數(shù)值。2.2.1函數(shù)的積分表達標準SPH是基于核估計發(fā)展起來的，在SPH方法中，對任意一個函數(shù) f，其積分近似表達式為其中，f(x)為函數(shù)f(x)的近似值，x為位置矢量，Q為包含x的幾分體，W(x-x ,h)為 光滑函數(shù)，其依賴于兩點之間的距離|x-x和光滑長度h。光滑函數(shù)在SPH近似法中其著重要的作用，它決定了函數(shù)表達式的精度和計算效率。光滑函數(shù) W常常選用偶函數(shù)，其應滿足 3個條件。第一個條件為正則化條件，如下所示x,h)dx, = 1Q由于光滑函數(shù)的積分值等于 1,故此條件也稱為歸一化條件。第二個條件是當光滑長度趨向于零時具有狄拉克函數(shù)性質h) - 3(x-xf)這里，cf lt X

16、= x1；0, x第三個條件是緊支性條件爐(X - X：力)=0,當 xxr Kh 時式中，K是與點x處光滑函數(shù)相關的常數(shù)，并確定光滑函數(shù)的有效范圍(非零)。此有效范圍稱作點x處光滑函數(shù)的支持域內的積分。因此，一般來說積分域Q就是支持域。根據(jù)分布積分和散度定理，將第一個式子經(jīng)過變換，對f (x)的導數(shù)的估計為 W) X - J7(X) W(x-xh)dx a由以上方程得知，SPH估計將函數(shù)的空間導數(shù)轉化為光滑函數(shù)的空間導數(shù)，因此，可 以根據(jù)核函數(shù)求取任意場函數(shù)的空間導數(shù)。222函數(shù)的粒子表達SPH方法最終要將函數(shù)積分近似表達式轉化為支持域內所有粒子疊加求和的離散化形式，流體的問題域被離散化為有

17、限數(shù)量的粒子，其中每個粒子具有獨立的質量、密度及其它物理屬性。經(jīng)過離散化，函數(shù)的積分表達式可以寫成如下的粒子近似式5,)*(詠其中，j和i表示粒子編號，mj和p j分別j粒子的質量和密度，N是粒子的總數(shù)，Wij=W(xi-xj , h)=W(|xi-xj | , h)。上式說明，在SPH方法中，對任意一個函數(shù)f，它在某一位置處的函數(shù)值可以通過應用 光滑函數(shù) W對其光滑長度為h的緊支域內所有粒子插值求和的形式來表達【17】。光滑長度h在SPH方法中非常重要，它在計算過程中決定了計算的實用性、有效性和 可靠的適應性。光滑長度h直接影響到計算結果的效率和結果的精度。若h太小，則在以kh為半徑的支持域

18、內沒有足夠的粒子對所求粒子施加影響，這樣會導致計算結果精度較低。如果h太大，太多的粒子對所求粒子產生影響，則可能將粒子所有細節(jié)信息和局部特性忽略，導致粒子趨于雷同，同樣會影響到結果的精度。具體光滑長度應該選擇多大與所研究的問題 相關，沒有一個廣義的光滑長度。2.2.3光滑函數(shù)SPH通過光滑函數(shù)來進行積分表示。光滑函數(shù)非常重要，因為它不僅決定函數(shù)表達式 的形式，定義粒子支持域的大小，而且還決定核近似和粒子近似的一致性和計算效率。常用的光滑函數(shù)有三種， 分別是高斯型光滑函數(shù)， 三次樣條光滑函數(shù)和高次樣條函數(shù)?？蓞⒁妳⒖嘉墨I【3,18,19,】2.3拉格朗日型的Navier-Stokes方程對控制方

19、程的描述有兩種形式：歐拉描述法和拉格朗日描述法。歐拉描述法是對空間點的描述，而拉格朗日描述法是對物質點的描述。流體的流動過程可以用歐拉方程以及適當?shù)臓顟B(tài)方程來表示，這些方程表示了質量、動量和能量守恒。如果用上標a和B來表示坐標方N-S方程組可以表示為:向和求和重復指標，采用時間全導數(shù)形式，連續(xù)性方程:Dp _ dvDt K dxp動量方程：Dv1Dtp dx3能量方程：坐=_蘭遼Dt p在以上方程中，6為總應力張量，有兩部分組成：一部分是壓力 P,另外一部分是粘性 應力T。B =卡護+嚴對于牛頓流體，粘性剪應力與剪應變率S成比例，且比例系數(shù)為粘性系數(shù)卩。其中戶詈+務弓（）嚴如果將壓力與粘性應力

20、分開寫，我們可以得到如下的能量方程De _ P&B 亠 P rap.p_ It t oDt p dx0 2p在上方式子中，各字母代表的意義如下：V 速度向量E內能P -密度p-壓力t-時間2.4 Navier-Stokes方程的 SPH 表達式根據(jù)光滑函數(shù)的概念和SPH方法的粒子近似原理，流體流動的運動方程可以寫成如下的離散SPH形式:dtdxrI dt /其中i參考粒子。J-粒子i的相鄰粒子。第三章SPH的計算實施方法3.1粒子的密度近似由于在SPH方法中粒子的運動是由相鄰粒子之間的密度差以及壓力差產生的，所以密 度近似在SPH方法中非常重要。在一般的SPH方法中主要有兩種方式對密度進行展開

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

22、必須在計算其他參數(shù)前計算密度，而且還要計算自身的光滑函數(shù)。近似密度的另一種粒子近似法為連續(xù)性密度法，這種近似法通過應用SPH近似法的概念對連續(xù)性方程進行轉換而得。常用的連續(xù)性密度方程為：連續(xù)性密度近似法不需要先計算密度，因此可以提高計算效率， 并適于并行處理算法的編程。具體應用哪一種密度算法要根據(jù)研究的問題而定， 對于廣義流體問題的模擬，應用修正 的密度求和法可得到較好的結果。對于具有強間斷問題的模擬，如爆炸、高速沖擊等， 應優(yōu) 先選取連續(xù)性密度法。3.2核函數(shù)核函數(shù)在SPH近似法中起著重要作用，因為它控制了插值的方式和粒子影響域的效能， 從而決定了函數(shù)表達式的精度和計算效率。由于SPH插值基

23、于積分插值理論，用核函數(shù)來近似狄拉克S函數(shù)，因此核函數(shù)應該具有某些特性，如非負性，緊支性，歸一性，衰減性和狄拉克S函數(shù)條件。在實際應用中，通常使用兩種形式的核函數(shù)，一種是如下所示的高斯核函數(shù)。1 1公式中的a d在一維空間中為，在一維空間中為，在二維空間中為另一種是三次樣條核函數(shù)，到目前為止，三次樣條函數(shù)是最廣泛應用的核函數(shù)。l-R2+R3 07?10(/?/) =勺 x *(2-幻1R2在一維空間中15，在二維空間中，在本文的三維空間中,3.3狀態(tài)方程在求解可壓縮流體問題的標準SPH方法中，粒子的運動是由于壓力梯度的作用而產生的，而粒子的壓力是通過狀態(tài)方程由粒子自身的密度和內能來計算的。然而

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

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

26、使用比真實值要小很多的聲速。所以對聲速大小的要求為：一方面，必須足夠大，以至于人工可壓縮流體的特性與真實流體充分接近；另一方面，應足夠小，可使時間步的增量在容許范圍內。3.4人工粘度對于某些數(shù)值問題， 誤差跳躍很大，會產生數(shù)值振蕩，為了提高算法的穩(wěn)定性，并且防止粒子間相互接近時的非物理穿透，Monaghan和Gingold在SPH方法的動量方程中引入了人工粘度【21】，它不是物理的粘度，只是數(shù)值的修正。0,式中,5 - 2（Ci + 勺）Pij =寺（Q +0）%弓他+勺）二旳一卩廠在式中，an, 3 n為常數(shù)，取值一般為 1.0左右【22】。因子p =0.1hij用于防止粒子 相互靠近時產生

27、的數(shù)值發(fā)散。c和v分別表示聲速和粒子的速度矢量。將人工粘度引入到動量方程中，最終動量方程的SPH粒子近似式可寫為如下形式：dt3.5邊界處理在SPH方法的應用中，邊界條件的處理既是該方法的優(yōu)點，也是目前的薄弱環(huán)節(jié)，因 為SPH方法最初是應用在沒有邊界條件限制的領域內，后來逐漸擴展到一些具有邊界的工 業(yè)應用中。我們使用虛粒子來處理固定邊界條件，通過設置邊界虛粒子可以方便地模擬固體邊界，同時也要防止粒子的非物理穿透。根據(jù)不同的數(shù)值仿真條件，我們在SPH數(shù)值模擬中一共使用了兩種類型的虛粒子第一種類型的虛粒子設置在固定邊界上，與Monaghan所使用的相似【4】。二種類型的虛粒子分布在邊界的外部，與L

28、ibersky所使用的相似 【23】。第二種類型的虛粒子按以下的方式構造， 即給定一個實粒子i，則在邊界外與實粒子對稱處分布一個虛粒子，這些虛粒子具有與相對 應實粒子相同的壓力和密度，但速度方向相反。第二種類型的虛粒子通常在邊界條件不斷變 化的場合下使用【24】。3.6時間積分可以使用一些標準的計算方法如跳蛙法 （LF）、預估校正法、龍格一庫塔法（RK）對普通微 分方程進行數(shù)值積分來求解每一個粒子的物理變化。我們使用跳蛙法進行時間積分，跳蛙法的優(yōu)點是計算時所需要的存儲量低， 而且在每一次計算中只需要進行一次優(yōu)化估值。粒子的密度、速度、內能和位移按照下面的公式遞推 ：t = t +A/pit +

29、Af/2)= pj( t 山/2)+ M - Dpt)Vi(t +Af/2)= vf( t ZV/2)+ A t Dvt(t)Ut(t +Af/2)= Ui( t A/2)+ A t - Dui(t)xi(t +40= xi(t)+ Ar vi(t +M/2)應該注意的是，跳蛙法是條件穩(wěn)定的， 穩(wěn)定的條件為 CFL（Courant Friedrichs 一 Levy） 條件，此條件一般會導致時間步長與光滑長度互成比例。在本文中，時間步長按如下方式取值：A/ = min(f勺 /肉W 片 + q + t2(anc- + 0口 |S vj)式中，E為Courant系數(shù)，一般取值在 0.3左右【4】

30、?！?】 宋康祖 緊支函數(shù)無網(wǎng)格方法研究【D】。清華大學博士學位論文。2000.【2】 Lucy L B.Numerical approach to testing the fission hypothesisj,Astronomical Journal.1977(82):1013-1024【3 】Gi ngold R A,M on agha n J J.Smoothed Particle Hydrod yn amics:Theory and Applicatio n toNon-sphericalstarsC.M on thlyNotices of the RoyalAstr ono mic

31、alSociety.1977(181)375-389.【4 】 Mon aghan J J.A n in troduction to SPHj.Comput.Phys.Commu n.1988,48:89-96.【5 】 Mo naghan 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 impactcomputationsJ.Int J Num Meth Eng,19

32、96,39:2725-2741.【7】 Swegle J W,Hicks D L, Attaways S W. Smoothed particle hydrodynamics stability an alysisJ .J Comput Phys,1995,116:123-134.【8 】 Dyka C T. Addressi ng ten sion in stability in SPH methodsR.tech nical Report Nrl/MR/6384,NRL,1994.【9】Johnson G R,stryk R A. Beissel S R.SPH for high velo

33、city impact computationsJ.Comput.Methods Appl.Mech.E ngrg.1996,139:347-373.【10 】Swegle J W, Attaway S W. On the feasibilility of using smoothed particle hydrody namics for underwater explosion calculationsj.Comput Mtch,1995,17:151-168.【11】 Joh nson G R,stryk R A,Beissel S R.SPH for high velocity imp

34、act computatio nsJ.ComputerMethods in Applied Mecha nics andEn gi neeri ng,1996,139:347-373【12 】Libersky L D,petschek A G,et al. High strain lagrangian hydrodynamics:A three-dime nsi onal SPH code for dyn amic material resp on seJ.J Comput Phys,1993.109:67-75【13 】Johnson G R,Beissel S R.Normalized smoothing functions for SPH impact computatio nsJ.l ntern ati onalJour nalfor Numerical Methods inEn gi neeri ng,1996,39:2725-2741.【14】