模擬流體沖擊致結(jié)構(gòu)破壞問(wèn)題的SPH-PD耦合方法

時(shí)浩天，郭 力

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

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

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

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

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

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

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

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

1 物理控制方程

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

1.1 SPH控制方程

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

(1)

(2)

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

歸一化條件

(3)

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

(4)

緊支屬性

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

(5)

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

(6)

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

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

(7)

(8)

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

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

連續(xù)性方程

(9)

動(dòng)量方程(無(wú)外力作用)

(10)

能量方程

(11)

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

連續(xù)性方程

(12)

動(dòng)量方程

(13)

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

(14)

(15)

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

(16)

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

(17)

能量方程

(18)

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

(19)

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

1.2 PD控制方程

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

b(x,t)

(20)

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

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

(21)

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

(22)

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

(23)

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

(24)

將上述斷鍵準(zhǔn)則引入點(diǎn)對(duì)力函數(shù)的表達(dá)式，可得更新后的點(diǎn)對(duì)力函數(shù)

(25)

當(dāng)鍵的破壞累計(jì)成一個(gè)面的時(shí)候，也就形成了宏觀的裂紋。每個(gè)點(diǎn)定義了局部損傷φ(x,t)

(26)

(27)

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

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

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

(28)

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

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

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

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

2.2 流固界面耦合策略

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

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

3 數(shù)值方法

3.1 流體域預(yù)測(cè)-校正時(shí)間步進(jìn)格式

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

(29)

根據(jù)第n步的物理量來(lái)預(yù)測(cè)第n+1/2步的物理量

(30)

由預(yù)測(cè)得到的第n+1/2步的物理量來(lái)校正第n+1/2步的物理量

(31)

最終計(jì)算得到n+1步的物理量

(32)

3.2 流體域基于CFL條件的變時(shí)間步長(zhǎng)法

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

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

(33)

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

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

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

(34)

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

4 數(shù)值算例

4.1 模型參數(shù)

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

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

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

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

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

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

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

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

(a) 0.04 s

彈性板自由端在水平方向(x方向)和垂直方向(y方向)的位移變化曲線見(jiàn)圖4，相應(yīng)的水位變化過(guò)程見(jiàn)圖5?？梢钥闯瞿M結(jié)果與實(shí)驗(yàn)數(shù)據(jù)在變化趨勢(shì)上吻合良好，因此所提出的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 水壩坍塌沖擊作用下彈性擋板的破壞模擬

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

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

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

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

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

5 結(jié) 論

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

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