模擬流體沖擊致結構破壞問題的SPH-PD耦合方法

時浩天，郭 力

(東南大學 土木工程學院，南京 210000)

流固耦合問題具有多場、多物理的特性，廣泛存在于工程應用與生產生活中。在滑坡、泥石流沖擊下游建筑，海嘯沖擊堤壩、海岸結構等自然災害場景，海上平臺受強海浪沖擊作用等工程問題中，往往涉及到自由表面流體與結構的劇烈相互作用。由于流體包含自由表面破碎、浪花飛濺和融合等高度非線性流動過程，同時結構在強沖擊作用下會發生大變形乃至破壞失效等復雜力學行為，使得這類問題的求解極具挑戰性。

在流固耦合問題中，基本方程包括流體和結構各自區域內的控制方程及其邊界條件，以及流固界面上力的平衡條件和幾何相容條件。將這些方程組納入到統一的框架中進行求解非常困難，已有的一些統一求解策略僅適用于簡單的流固耦合問題。相應地，分區求解方法越來越受到人們的關注。分區求解法將流體和結構在各自區域分別求解，通過在流固界面上傳遞數據(壓力、速度、位移等)來實現耦合計算。分區求解法的優勢在于可以最大程度地利用已有的流體和結構求解方法，綜合現有的流體和結構求解模塊，通過求解器之間的數據傳遞實現耦合求解，從而顯著降低了開發成本。

構建流固界面上數據的有效傳遞策略是分區求解法中的關鍵問題。流固界面往往是時變的，為了精確捕捉流固界面的位置，已經發展了任意拉格朗日-歐拉方法[1](arbitrary Lagrangian-Euler，ALE)和歐拉-拉格朗日耦合方法[2](coupled Euler-Lagrangian,CEL)。這些方法主要采用動網格技術，捕捉流固界面的位形，同時保證每個時間步內流體和結構在各自網格區域內正確運算。動網格技術能夠有效處理較小變形下的流固耦合問題，但是在處理任意復雜的流固界面時，追蹤流固界面需要網格進行大量的處理與重構，過程非常繁瑣。浸入邊界法[3](immersed boundary,IB)采用固定網格，避免了動網格技術的缺點，但是在流固界面上的分辨率和精度較低。IB在高雷諾數條件下的應用較為有限，多適用于流體粘性較大的生物醫學問題。

對于分區求解方法，通過組合現有計算流體力學和固體力學中的求解框架，發展了眾多的流固耦合求解方法。Wick[4]將ALE方法和相場法相耦合，用于求解流固耦合問題，同時可以考慮結構的脆性斷裂行為；Zheng 等[5]將有限元法與邊界元法耦合，研究了浸沒在無限可壓縮流體中彈性結構的自由振動問題。上述基于網格的方法可以有效地處理流固耦合問題，但是當流體域包含自由表面、結構域涉及大變形時，對網格的處理會導致高昂的計算成本。此外，當結構域涉及到裂紋的萌生及擴展和材料失效時，會產生幾何上的不連續和奇異性問題，基于網格的方法將存在固有的困難。針對這類問題的求解，無網格法具有獨特的優勢。因此，發展一套具有魯棒性和高精度的無網格類方法，來求解包含自由流體表面、結構大變形和材料斷裂失效的復雜流固耦合問題，無疑具有重要的意義。

光滑粒子流體動力學(smoothed particle hydrodynamics，SPH)是一種具有完全拉格朗日屬性的無網格粒子類方法，廣泛應用于自由表面流、多相流、非牛頓流體、湍流、考慮表面張力的流體流動的模擬以及化學擴散和沉淀分析等[6]。此外，SPH也可以用來模擬固體，比如彈性體、柔性體[7]。沈雁鳴等[8]采用SPH統一算法模擬了自由流體沖擊彈性體問題，證明了SPH具有獨立處理流固耦合問題的能力；施書文等[9]將SPH和光滑點插值法結合，模擬了彈性體入水問題。在泥石流、滑坡等災難仿真以及近海結構受強海浪沖擊作用等場景下，考慮結構的破壞和失效是十分有必要的，但是考慮材料和結構斷裂失效過程的SPH流固耦合模型還很少。

近來，Silling等[10]提出的近場動力學(Peridynamics,PD)采用非局部積分方程，不連續位移場可以自然地包含在控制方程中，因而在模擬材料中裂紋萌生、擴展方面具有天然的優勢，已經成為研究熱點。目前，PD已經成功用于混凝土板侵徹問題[11]、充液裂縫分析[12]和準脆性材料的沖擊響應和損傷擴展模擬[13]。此外，Gao等[14]利用近場動力學微分算子離散Navier-Stokes方程(簡稱N-S方程)，實現了基于近場動力學統一求解框架的流固耦合數值模型，但是其計算效率還有待于進一步研究。近場動力學的理論基礎不同于無網格法，但是其數值實施策略與無網格法十分類似。Ganzenmüller等[15]詳細討論了在粒子離散條件下PD與SPH的相似性，結果表明，運用SPH或PD模擬基于變形梯度的經典材料模型，采用節點積分時，兩者的數值控制方程是一致的，這意味著在某種情況下PD退化為無網格方法。這為聯合SPH和PD來求解流固耦合問題提供了理論上的依據。

本文將SPH和PD相結合，發展了一套能夠考慮流體自由表面和結構斷裂失效的無網格數值模型。基于動態邊界條件思想，構建了流固界面上的數值傳遞策略，最終達到耦合求解的目的。該模型對于結構斷裂過程中產生的新流固界面無需額外的處理，計算過程簡潔。通過應用人工狀態方程、人工黏性、預測-校正時間步進格式和結構區域的時間凍結技術，改善了計算精度和效率。以潰壩沖擊下游彈性擋板問題為例，驗證了模型在捕捉自由液面和結構破壞方面的有效性，為模擬一定規模的流固耦合問題提供了基礎。

1 物理控制方程

采用分區法求解流固耦合問題時，流體域和固體域采用各自獨立的控制方程，下面分別給出流體域的SPH控制方程和結構域的PD控制方程。

1.1 SPH控制方程

SPH最初用于解決天體運動過程之間的相互作用問題，Monaghan將其引入流體仿真領域，并率先用SPH來模擬自由表面流[16]。SPH采用如下形式來近似場函數f(x)及其在x位置的空間導數?·f(x)

(1)

(2)

式中，Ω為流體區域，W(x-x′,h)為權函數或核函數，h為光滑長度。當定義了核函數的支持域后，支持域半徑一般取光滑長度的k倍。核函數W一般選擇偶函數，并且具有以下三個性質：

歸一化條件

(3)

δ函數性質

(4)

緊支屬性

W(x-x′,h)=0, |x-x′|>kh

(5)

這里采用Wendland[17]提出的五次樣條核函數：

(6)

式中,R=x/h，對于2D問題αd=7/(4πh2)。

上面的〈f(x)〉稱為函數f(x)的一個核近似。通過核近似將原場函數及其空間導數近似為核函數支持域上的積分，由此導出SPH粒子近似：

(7)

(8)

式中,mj為粒子質量，ρj為粒子密度，N為粒子i的鄰域粒子數。粒子近似將空間上連續的積分運算離散為域內有限個粒子的求和運算。

拉格朗日描述下的流體控制方程為非守恒形式的N-S方程，采用指標記法為：

連續性方程

(9)

動量方程(無外力作用)

(10)

能量方程

(11)

利用SPH核近似和粒子近似可以得到SPH離散形式的流體控制方程，簡稱SPH控制方程，其中ρ、m、v、P、e分別表示粒子的密度、質量、速度、壓力、內能。式中，下標a表示與流體粒子a相關的物理量，下標b表示粒子a的支持域內的其他流體粒子b的相關物理量，Wab表示粒子a的核函數在粒子b處的值，vab=va-vb，rab=ra-rb：

連續性方程

(12)

動量方程

(13)

式中，g為重力加速度，大小為9.81 m/s2。P為流體壓力，此處將不可壓縮流體視為弱可壓縮流體，通過引入人工狀態方程來顯式求解壓力項，避免了求解壓力泊松方程。流體粒子的壓力由下式計算

(14)

(15)

Πab為黏性項，此處采用人工黏性，其值由下式確定

(16)

式中，α為第一黏度系數，上劃線和下標ab表示粒子a和b相應物理量的平均值。

(17)

能量方程

(18)

在通過SPH控制方程移動粒子時，需要引入XSPH修正[16]來有效避免粒子堆積，使得粒子分布更加均勻，同時防止粒子之間發生非物理穿透。粒子規則移動由下式定義

(19)

ε為修正系數，取ε=0.5。

1.2 PD控制方程

根據牛頓第二定律，近場動力學的控制方程可以表達為如下形式

b(x,t)

(20)

式中，H為空間內物質點x的近場范圍，u為物質點的位移矢量場，b為一個預定義的體力密度場，ρ為物質密度，f是PD模型中的點對力函數，它包含了材料常數，類似于傳統理論中的本構關系，此處的f反映了物質點x′對物質點x的作用效果。

為了能夠模擬裂紋在變形體中自然演化的過程，Silling[10]在2005年提出了經典微觀彈脆性(prototype microelastic brittle，PMB)本構模型，適用于各向同性材料。記ξ=x′-x，η=u′-u，則PMB模型中點對力函數f(η,ξ)可表達為如下形式

(21)

式中,c是微彈模量。通過對比PD理論中的應變能密度表達式和經典連續介質力學中的應變能密度表達式可得

(22)

(17)中S是物質點對的伸長率，其定義式為

(23)

鍵型近場動力學中，材料的損傷可以通過物質點之間“鍵”的斷裂來描述。當兩個物質點之間的伸長率S達到臨界值Scr時，“鍵”即斷裂，物質點之間的相互作用消失。該過程是不可逆的，臨界伸長率的具體數值可以通過材料宏觀的抗壓強度Sc0、抗拉強度St0來定義

(24)

將上述斷鍵準則引入點對力函數的表達式，可得更新后的點對力函數

(25)

當鍵的破壞累計成一個面的時候，也就形成了宏觀的裂紋。每個點定義了局部損傷φ(x,t)

(26)

(27)

通過統計一點近場范圍內的斷鍵數目來計算φ(x,t)，從而模擬裂紋的萌生和擴展。鍵型近場動力學求解的數值過程不依賴網格，通過離散節點攜帶的信息，求解域被離散為有限的物質點，邊界條件可以直接施加到材料的邊界物質點層上，其施加邊界條件的具體形式如下：

u(x,t)=U0,x∈Rc

b(x,t)=p(x,t),x∈Rc

(28)

2 基于動態邊界條件的耦合策略

2.1 動態邊界條件

根據動態邊界條件的思想，圖1(a)中的剛體粒子與流體粒子滿足同樣的SPH控制方程。對于靠近邊界的流體粒子a1，當其支持域內出現剛體粒子時，剛體粒子對其產生的影響與支持域內其他流體粒子對其產生的影響按照同一套控制方程來計算，即按照式(12)更新密度，依據式(13)更新加速度。這時剛體粒子被看作具備流體粒子的特性，這類粒子稱為動態邊界粒子。與流體粒子不同，剛體粒子的位置不隨時間改變，速度始終保持為零，從而實現固定邊界的效果。由于剛體邊界粒子被當做流體粒子參與計算，不需要做額外的判別，動態邊界條件可以在程序中高效統一實現。

(a) 流體粒子受動態邊界粒子作用 (b) 彈性體粒子受流體粒子的反作用圖1 SPH-PD耦合策略Fig.1 SPH-PD coupling strategy

2.2 流固界面耦合策略

考慮到動態邊界條件的便捷性，在SPH-PD模型中同時將其應用到流固界面的處理上。在圖1(a)中，對于靠近流體-彈性體界面的流體粒子a2，將彈性體粒子也作為動態邊界粒子處理，此時彈性體粒子對流體粒子的影響與剛體粒子對流體粒子的影響完全一樣。與剛體粒子不同的是，彈性體粒子的位置和速度隨時間變化，其變化量由彈性體區域內的PD控制方程所確定。在彈性體粒子作為動態邊界粒子對流體粒子產生作用力時，流體粒子對彈性體粒子有大小相等、方向相反的反作用力，如圖1(b)所示。為簡明起見，這里只給出了一個流體粒子a2對彈性體的反作用力，實際彈性體受到的反作用力應該是全部受其影響的流體粒子所產生反作用力的總和。計算該作用力，并在一個時間步內作為載荷施加到彈性體上，參與PD控制方程的運算，從而實現流固界面上的數據交換。

當計算彈性體粒子對流體粒子的影響時，依據動態邊界條件的思想，處于邊界的彈性體粒子被視為流體粒子。因此，在此過程中彈性體粒子的質量、密度、體積均取SPH粒子的值，但其位移和速度不按照SPH控制方程計算；當任一時間步內流體粒子的運動計算完成后，再開始計算彈性體粒子在流體粒子反作用下的運動。這時，彈性體粒子的參數應取實際的固相材料參數，通過求解PD控制方程，即可計算得到結構在流體反作用力下的變形和位移。

3 數值方法

3.1 流體域預測-校正時間步進格式

流體域由SPH離散形式的N-S方程組控制。在第n+1個時間步內，流體在第n步的密度、速度、能量的變化率由式(12)、(13)、(18)計算得到：

(29)

根據第n步的物理量來預測第n+1/2步的物理量

(30)

由預測得到的第n+1/2步的物理量來校正第n+1/2步的物理量

(31)

最終計算得到n+1步的物理量

(32)

3.2 流體域基于CFL條件的變時間步長法

為了防止粒子間的非物理穿透，流體域內的時間步長不能任意選取，必須滿足時間步長CFL條件(courant-friedrichs-levy)。進行時間積分時，采用滿足CFL條件的最大時間步長。每一步的時間步長Δt由下式確定：

Δt=0.3·min(Δtf,Δtcv)

(33)

3.3 結構域時間凍結技術

結構域由PD積分方程控制，在流體域和結構域選擇近似粒子間距的情況下，依據CFL條件確定的結構域最大時間步長ΔtPD往往遠小于流體域的最大時間步長ΔtSPH。進行時間積分時，應該取較小的時間步長ΔtPD作為耦合算法整體的時間步長。由于流體域的計算量占總體計算量的較大部分，而且耗時比結構域更長，因此對于流體域，選擇時間步長ΔtPD進行運算，會增加很多不必要的迭代次數，降低了算法整體的效率。為了提高計算效率，這里在流體域和結構域分別采用不同的時間步長進行計算。具體的時間步長選擇和計算策略為：

在第n個時間步內，若流體域和結構域根據CFL條件計算得到的最大時間步長分別為ΔtSPH和ΔtPD，則流體域按照ΔtSPH迭代一次，結構域按照ΔtPD迭代M次。

(34)

這樣，相對于結構域，流體域內的時間就仿佛被凍結了一樣，因此稱為時間凍結技術。在流體域時間凍結期間，流體粒子對彈性體粒子的反作用保持不變。結構在該載荷作用下，通過求解PD控制方程，迭代M步來計算變形和位移。值得注意的是，在這M次迭代過程中，若任一彈性體粒子的位移累計超過光滑長度h，則結構域的迭代終止，流體域解除時間凍結，進行下一步迭代，重新計算對結構的反作用。

4 數值算例

4.1 模型參數

計算參數的選取對計算結果往往會產生一定影響，這里的算例采用統一模型參數，詳見表1。

表1 SPH-PD模型統一參數設定Tab.1 Common Parameters of SPH-PD model

計算中流體域和結構域采用相同粒子間距，流體域時間步長由式(33)確定；結構域時間步長在初始時刻通過CFL條件確定，迭代中不動態改變，通過迭代次數M的變化來間接調整結構域迭代的總次數。

4.2 彈性擋板在水流作用下的變形模擬

為了驗證所提出的SPH-PD模型的有效性，首先模擬流體流動和結構變形的相互作用。Antoci等[18]試驗研究了彈性擋板控制下的排水問題，并用SPH-SPH統一法進行了數值模擬。采用所提出的SPH-PD模型對該問題進行模擬，并與實驗結果數據和SPH-SPH統一法模擬結果進行對比。模型的初始條件設置如圖2所示，高為H1的水柱被限制在寬為W1的水槽內，水槽右側和下方為剛性邊界，左側有一高為L1的彈性擋板，彈性擋板的上端固定，下端自由。模型的幾何參數，水的密度ρF，彈性板的密度ρS、彈性模量E、泊松比μ見表2。

圖2 彈性擋板下排水問題初始條件設置Fig.2 Initial condition settings of drainage problem with elastic plate

表2 彈性擋板下排水問題參數設定Tab.2 Parameters of drainage problem with elastic plate

將相關的模擬結果與實驗結果進行對比，選取了6個典型時刻的結果，如圖3所示：第一列為實驗結果，第二列為基于本文提出的SPH-PD模型模擬的結果，第三列為Antoci等采用SPH-SPH統一方法模擬的結果。起初彈性擋板在水壓作用下發生變形，彈性板自由端向左上方偏移，導致下方出現空隙，從而形成流體排出的通道。在0.16 s以前彈性板的變形在水壓作用下不斷增大，但隨著槽內流體的排出，水位逐漸降低，水壓漸漸變小，彈性板的變形開始逐漸減小。通過對比可以看出：基于SPH-PD模型的模擬結果很好地再現了實驗中水位的變化和彈性板的變形過程。

(a) 0.04 s

彈性板自由端在水平方向(x方向)和垂直方向(y方向)的位移變化曲線見圖4，相應的水位變化過程見圖5。可以看出模擬結果與實驗數據在變化趨勢上吻合良好，因此所提出的SPH-PD模型是有效的。

圖4 彈性板自由端水平、垂直方向位移Fig.4 Horizontal and vertical displacement of free end of elastic plate

圖5 彈性板右側水位Fig.5 Water level on the right side of elastic plate

4.3 水壩坍塌沖擊作用下彈性擋板的破壞模擬

為了驗證SPH-PD模型模擬結構損傷斷裂的能力，模擬潰壩沖擊下游彈性擋板的破壞過程。在二維平面上有一正方形剛性封閉容器如圖6所示。容器邊長a=4 m，初始狀態下，水柱位于容器內左側1/4處，水柱的寬度和高度分別為，W=1 m和H=3 m。一彈性擋板位于方形容器中部，底端固定，彈性擋板的厚度和高度分別為，l=0.1 m和h=1 m，與水柱的距離為L=1 m。擋板采用線彈性本構模型，考慮損傷，密度ρ=1 878 kg/m3，彈性模量E=20 MPa，泊松比μ=0.333 3，當鍵的伸長率超過拉壓臨界伸長率St0和Sc0時，鍵即斷裂，取St0=0.05和Sc0=-0.09。為了增加流動過程的復雜程度，同時在下游設置一個寬W2=1.9 m、高H2=0.7 m流體區域。水柱在重力作用下自由坍塌沖擊中部彈性擋板，利用SPH-PD模型計算流體和結構的運動變形過程，計算總時長為3 s。

圖6 潰壩沖擊下游彈性擋板模型初始條件設置Fig.6 Initial condition settings of elastic plate with water downstream impacted by dam break

選取了潰壩沖擊過程中幾個典型時刻的粒子位形及壓力分布情況，如圖7所示。

圖7 彈性擋板在潰壩沖擊作用下的斷裂過程及t=1.2 s時結構區域損傷場分布Fig.7 Fracture process of elastic plate under dam break impact and damage distribution of structure at t=1.2 s

可以看出，左側水柱被釋放之后撞擊彈性擋板，同時形成高達3 m的浪尖。彈性擋板在水流沖擊作用下先發生大變形，變形過程中由于局部鍵的伸長率超出臨界值而發生斷裂。1.2 s時可以清楚地看到結構區域發生斷裂，斷裂截面上粒子的損傷量為50%。擋板斷裂之后，坍塌的水柱與下游水域發生融合，并攜帶斷裂部分的擋板共同運動。SPH-PD模型可以模擬出流體沖擊致結構破壞的典型過程和細節，是有效可靠的。

5 結 論

本文將SPH和PD相結合，提出了一套模擬流固耦合問題的無網格數值方法。所提出的SPH-PD模型采用分區求解策略：SPH用于模擬自由表面流，PD用于模擬結構的變形、斷裂行為，通過流固界面上基于動態邊界條件的思想實現耦合運算。采用預測-校正時間步進格式，提高了時間積分精度。利用時間凍結技術，改善了流體域內的計算效率，同時采用人工黏性和XSPH修正，改善了粒子分布的均勻性。利用提出的SPH-PD模型，模擬了二維潰壩沖擊下游彈性擋板問題，可以高效捕捉到流體自由表面破碎后出現的浪花翻滾、飛濺、融合等復雜現象，同時實現了擋板中裂紋的萌生、擴展直到斷裂失效的全過程模擬，從而驗證了模型的有效性和魯棒性。

SPH-PD模型具有無網格粒子類方法的典型優勢，可以拓展到更多包含自由表面流、結構損傷失效的流固耦合問題，可以為災難仿真、影視特效等領域提供了一種新的模擬方法。