光滑粒子流體動(dòng)力學(xué)方法SPH課件

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

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

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

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

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

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

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

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

9、止了粒子相互接近時(shí)的非物理穿透。具體表達(dá)式如下： 人工粘度 在以上所有方程中 和 為標(biāo)準(zhǔn)常數(shù)，一般取值在1.0左右。因子 用于防止粒子相互靠近時(shí)產(chǎn)生的數(shù)值發(fā)散。式中與 相關(guān)的項(xiàng)得到的是體積粘度，而與 相關(guān)的項(xiàng)是用于防止在高馬赫數(shù)時(shí)粒子的相互穿透。式中所給出的人工粘度被引入到SPH方程的壓力項(xiàng)中。邊界處理 由于在邊界上或鄰近邊界處的粒子存在缺陷，即在積分的時(shí)候會(huì)被邊界截?cái)?，故而SPH方法不能完全適用于整個(gè)區(qū)域。在邊界上或鄰近邊界處的粒子只受到邊界內(nèi)的粒子的影響作用，而邊界外由于沒有粒子，故而邊界外不對(duì)粒子產(chǎn)生影響。這種單邊影響作用會(huì)導(dǎo)致求解結(jié)果錯(cuò)誤，因?yàn)樵诠潭ㄟ吔绫砻妫m然粒子速度為零，但是其他

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

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

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

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

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