光滑粒子流體動力學方法SPH課件

1、SPH光滑粒子流體動力學方法(smoothed particle hydrodynamics)報告人：馬天寶2013.4.25無網(wǎng)格方法的主要思想：通過使用一系列任意分布的節(jié)點（或粒子）來求解各種各樣邊界條件的積分方程或偏微分方程組(PDEs)從而得到精確穩(wěn)定的數(shù)值解，這些節(jié)點或粒子之間不需要網(wǎng)格進行連接。 Lucy, Gingold(1977)分別提出了SPH方法，最早用于天體物理現(xiàn)象的模擬，隨后別廣泛地應用于連續(xù)固體力學和流體力學中。 11kg彈丸1418m/s撞靶速度下穿靶過程的數(shù)值模擬1500m/s速度下彈體侵徹混凝土靶變形過程的數(shù)值模擬碎浪與彈性擋墻之間的相互作用近似函數(shù)構造方法偏微

2、分方程的離散形式核估計方法(Kernel Approximation, KA)移動最小二乘法(Moving Least Square, MLS)再生核估計方法(Repuducing Kernel Method, RKM)徑向基函數(shù)方法(Radial Basic Function, RBF)單位分解方法(Patition of Unity, PU)強形式以各種全局或局部加權余量法為統(tǒng)一框架的弱形式兩條主線無網(wǎng)格法(Meshfree Methods)強形式是直接從微分方程及其定解條件出發(fā)，將近似函數(shù)及其導數(shù)的估計形式帶入基本方程、本構方程和初邊值條件中去，聯(lián)立方程進行求解。該方法思路簡單，便于程序

3、編制，應用范圍廣泛，在流體和固體的計算中都有所發(fā)展，適用于計算激波、高速沖擊、爆轟、穿甲等沖擊動力學問題。但此類算法的精度較低，穩(wěn)定性較差，且邊界條件的引入比較困難。弱形式就是從加權余量法或變分原理出發(fā)，把微分方程及其定解條件轉換成弱形式(Weak Form)或Galerkin形式，即用測試函數(shù)(Test function)與控制方程相乘后在全局或部分區(qū)域內積分，并利用高斯散度定理得到不同形式的弱形式，然后進行離散化求解。通過引進新的無網(wǎng)格近似函數(shù)構造方法或采用新的偏微分方程的離散形式，就可以期待開發(fā)出更加高效和精確的無網(wǎng)格方法。 SPH方程的構造常按兩個關鍵步進行。第一步為積分表示法，又稱場

4、函數(shù)近似法；第二步為粒子近似法。光滑粒子流體動力學一種無網(wǎng)格粒子法，湖南大學出版社，G. R. Liu，M. B. Liu著場函數(shù)核近似法（積分表示法）函數(shù)核近似的標準表達式：h是定義光滑函數(shù)W的影響區(qū)域的光滑長度。W被稱為光滑核函數(shù)(smoothing kernel function)或光滑函數(shù)(smoothing function)，簡稱為核(kernel)函數(shù)。粒子近似法 與SPH核近似法相關的連續(xù)積分表示式，可轉化為支持域內所有粒子疊加求和的離散化形式。SPH計算公式光滑函數(shù)最近相鄰粒子搜索法(NNPS)人工粘度邊界處理交界面處理光滑長度的更新SPH方程的求解激波管問題SPH程序結構目

5、 錄1、密度的粒子近似法 由于粒子的分配與光滑長度的變化主要依賴于密度，故在SPH法中密度近似法非常重要。 在SPH法中有兩種方法對密度進行展開，第一種方法是對密度直接用SPH近似法，稱為密度求和法。第二種方法是連續(xù)性密度法，通過應用SPH近似法的概念對連續(xù)性方程進行轉換而得到。SPH計算公式密度求和法：改進方案（正則化）此方法可提高自由邊界處和相同材料粒子密度不連續(xù)交界面處的精度 對于廣義流體問題的模擬，應用修正的密度求和法可得到較好的結果，對于具有強間斷問題的模擬（如爆炸、高速沖擊等），應優(yōu)先選取連續(xù)性密度法。2、動量方程的粒子近似法 將以上兩式相加可得：將動量方程等號右端的梯度項直接應用

6、SPH粒子近似法進行變換得：此外，有：動量守恒方程：此對稱方程的優(yōu)點為：可降低粒子不一致問題產生的誤差。3、能量方程的粒子近似法 能量守恒方程：光滑函數(shù)光滑函數(shù)的性質：一、正則化條件由于光滑函數(shù)的積分值等于1，故此條件也稱為歸一化條件。二、當光滑長度趨向于零時具有狄拉克函數(shù)性質三、緊支性條件Monaghan和Lattanzio在三次樣條函數(shù)的基礎上提出了稱為B樣條函數(shù)的光滑函數(shù) ：光滑函數(shù)現(xiàn)有SPH文獻中最為廣泛應用的光滑函數(shù)在一維、二維和三維空間中分別有： ， 和 。三次樣條函數(shù)及其一階導數(shù)四次樣條函數(shù)及其一階導數(shù)五次樣條函數(shù)及其一階導數(shù)最近相鄰粒子搜索法(NNPS)一般將包含在支持域中的粒

7、子稱為相關粒子的最近相鄰粒子(NNP)。尋找最近相鄰粒子的過程通常稱為最近相鄰粒子搜索(NNPS)。在SPH方法中常用的三種NNPS方法為：全配對搜索法(all-pair search)鏈表搜索法(linked-list search algorithm)樹形搜索法(tree search algorithm)全配對粒子搜索法對于給定的粒子i，應用全配對搜索法即是計算粒子i到粒子j的距離。若該距離小于粒子i的支持域的半徑 ，則粒子j為粒子i的支持域內的粒子。 樹形搜索法非常適宜求解具有可變光滑長度的問題。這種算法是通過粒子的位置來構造有序樹。一旦樹形結構構造起來后，便能高效地搜索最近相鄰粒子。

8、 樹形搜索法將最大問題域遞歸分割成一個個卦限，直到每一個卦限內只包含一個粒子為止。樹形結構構造完成后，即可以開始進行最近相鄰粒子搜索。 給定任一粒子i, 并以i為中心，用邊長為 的立方體將粒子包圍起來，然后再檢測粒子i的搜索立方體空間是否與并列層次內的其他節(jié)點所占的空間有重合的地方。若沒有，則中止往下搜索；若有，則繼續(xù)往下一層搜索。直到所搜索到的當前節(jié)點處只有一個粒子為止。接著，檢查此粒子是否在給定粒子的支持域內，若是，則將其記為粒子的相鄰粒子。 迄今為止，與SPH相關的論著中，Monaghan型的人工粘度是最為廣泛使用的人工粘度。它不僅將動能轉換為熱能，提供了沖擊波陣面必不可少的耗散，而且防

9、止了粒子相互接近時的非物理穿透。具體表達式如下： 人工粘度 在以上所有方程中 和 為標準常數(shù)，一般取值在1.0左右。因子 用于防止粒子相互靠近時產生的數(shù)值發(fā)散。式中與 相關的項得到的是體積粘度，而與 相關的項是用于防止在高馬赫數(shù)時粒子的相互穿透。式中所給出的人工粘度被引入到SPH方程的壓力項中。邊界處理 由于在邊界上或鄰近邊界處的粒子存在缺陷，即在積分的時候會被邊界截斷，故而SPH方法不能完全適用于整個區(qū)域。在邊界上或鄰近邊界處的粒子只受到邊界內的粒子的影響作用，而邊界外由于沒有粒子，故而邊界外不對粒子產生影響。這種單邊影響作用會導致求解結果錯誤，因為在固定邊界表面，雖然粒子速度為零，但是其他

10、變量，如密度，則不一定為零。邊界處理 在Liu等的研究中，使用了虛粒子來處理固定邊界條件，提出了兩種類型的虛粒子。第一種類型的粒子（型號I）設置在固定邊界上；第二種類型的粒子（型號II ）分布在邊界的鄰域內。型號II 的虛粒子具有與相應實粒子相同的密度和壓力，但速度相反。型號I的虛粒子被引入到實粒子的核函數(shù)核粒子近似法中。邊界處理 當類型I的虛粒子成為鄰近邊界處的實粒子的相鄰粒子時，則會在沿著兩粒子的中心線處對實粒子產生一個作用力。式中：參數(shù)n1和n2一般取值分別為12和4。D是由具體問題而定的參數(shù)，一般取與速度最大值的平方相等的量級。截止半徑r0在此問題的模擬分析中非常重要，在一般情況下，

11、r0的取值與粒子的初始間距的大小相近。 型號II的虛粒子沒有固定的參數(shù)。它們是在每個計算步中由對應的實粒子對稱產生的?？梢詰眯吞朓I的虛粒子來處理固定邊界和自由表面。 通過數(shù)值算例的測試驗證了應用虛粒子來處理邊界的可靠性和有效性。其不僅僅提高了SPH近似法在邊界區(qū)域處的精度，而且防止了粒子非物理穿透邊界。交界面處理 由于SPH方法具有拉格朗日性質和粒子性質，在整個演變過程中，來自不同介質的相互接觸的粒子可能會隨著運動而分離，甚至有可能不再相鄰。交界面處理 如果求和僅僅是使用相同材料的粒子，則離得很近但是材料不同的兩粒子就不再是相鄰粒子，因此，在交界面附近的區(qū)域，模擬就會遇到如一般邊界的粒子缺

12、陷問題。相反，如果求和使用不同材料的粒子，那么離得很近但是材料不同的兩粒子就可以認為是相鄰粒子，因此能夠減少在SPH核近似和粒子近似中的邊界缺陷問題。然而，這種處理將會導致一些非物理穿透或摻雜的問題。在大部分的情況下非物理穿透都不是很嚴重，僅僅有一或兩層不同材料的粒子發(fā)生相互穿透。對于具有高強度載荷作用的問題，穿透對精度來說可能是致命的，有時甚至會造成SPH程序運行的崩潰。 Liu等人提出了一種粒子對粒子的交界面算法，無論是對于高速沖擊問題還是水下爆炸問題這種算法都是非常有效的。在此算法中，物質交界面通過核近似求和來處理，即在求解守恒方程時考慮了不同材料粒子之間的相互作用。然而，這樣處理還不是

13、很完善，因為在水下爆炸過程中，高壓爆炸氣體和周圍的水介質之間的激烈相互作用通常會導致交接面附近的非物理穿透和摻雜。因此在交界面附近的不同材料的粒子，當它們彼此之間趨于穿透時就對其應用一個力的特殊處理方法。如果以下條件滿足就認為穿透發(fā)生： 所施加的懲罰力是成對地作用在兩個相互接近的粒子上，力的方向是沿著粒子的中心線方向。式中的參數(shù) 和n1、n2 分別取105，6和4。 懲罰力與接觸粒子求和的結合應用，雖然在交接面附近仍然存在著數(shù)值振蕩，但是卻能夠很好地預防在水下爆炸模擬中出現(xiàn)的非物理穿透問題。光滑長度的更新 在SPH方法中，光滑長度非常重要，直接影響到求解的計算效率和精度。在下一個時間步，光滑長

14、度變?yōu)椋?Benz(1989)提出了一種對光滑長度h進行動態(tài)變換的方法，即在連續(xù)性方程中將光滑函數(shù)對時間求導： 為了模擬具有極大密度不均勻的問題，Liu等推導出了以下優(yōu)化和松弛過程，此過程的目的是使每個粒子與數(shù)目近似為常量的最近相鄰粒子相互作用，由兩步組成：預測步和修正步。預測步：修正步：一旦得到 ，則可以確定當前的相鄰粒子數(shù) 。若 與 大致相近，則其是可取的。若誤差為百分之幾，則要在1.0附近對松弛因子 進行調整，以得到新的 和 。重復此過程，直到兩者大致相近為止。其中 為松弛因子，初始取值為1.0，然后在隨后的優(yōu)化和松弛過程步驟中在1.0附近進行調節(jié)。優(yōu)化和松弛過程 SPH方程的求解 上述的SPH方程是一系列與時間相關的聯(lián)立的常微分方程，可以使用標準的數(shù)值分析算法進行求解，如：leapfrog(LF)、Predictor-corr