時浩天，郭 力

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

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

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

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

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

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

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

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

1 物理控制方程

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

1.1 SPH控制方程

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

(1)

(2)

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

歸一化條件

(3)

δ函數(shù)性質(zhì)

(4)

緊支屬性

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

(5)

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

(6)

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

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

(7)

(8)

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

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

連續(xù)性方程

(9)

動量方程(無外力作用)

(10)

能量方程

(11)

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

連續(xù)性方程

(12)

動量方程

(13)

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

(14)

(15)

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

(16)

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

(17)

能量方程

(18)

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

(19)

ε為修正系數(shù)，取ε=0.5。

1.2 PD控制方程

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

b(x,t)

(20)

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

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

(21)

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

(22)

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

(23)

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

(24)

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

(25)

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

(26)

(27)

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

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

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

(28)

2 基于動態(tài)邊界條件的耦合策略

2.1 動態(tài)邊界條件

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

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

2.2 流固界面耦合策略

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

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

3 數(shù)值方法

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

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

(29)

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

(30)

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

(31)

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

(32)

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

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

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

(33)

3.3 結(jié)構(gòu)域時間凍結(jié)技術(shù)

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

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

(34)

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

4 數(shù)值算例

4.1 模型參數(shù)

計算參數(shù)的選取對計算結(jié)果往往會產(chǎn)生一定影響，這里的算例采用統(tǒng)一模型參數(shù)，詳見表1。

表1 SPH-PD模型統(tǒng)一參數(shù)設(shè)定Tab.1 Common Parameters of SPH-PD model

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

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

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

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

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

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

(a) 0.04 s

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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