吳巧瑞，王 磊，謝永和，白興蘭

(浙江海洋大學(xué)船舶與機(jī)電工程學(xué)院，浙江舟山 316022)

潰壩現(xiàn)象的物理模型是使用擋板在水箱內(nèi)維持一定高度的水柱，在初始時刻迅速撤去擋板，在重力的作用下，水柱會坍塌形成一個自由表面波，并對周圍的障礙物進(jìn)行砰擊作用的一種現(xiàn)象，見圖1。此現(xiàn)象廣泛存在于各個領(lǐng)域，尤其在船舶與海洋結(jié)構(gòu)物設(shè)計制造、海洋水動力學(xué)領(lǐng)域具有廣泛的應(yīng)用背景，經(jīng)常用作甲板上浪、液艙晃蕩、入水沖擊、海洋平臺遭受波浪載荷等工程實(shí)際問題的簡化來研究，其共同點(diǎn)是在砰擊過程中具有流體的大變形和強(qiáng)間斷性。因此，盡管潰壩現(xiàn)象的物理模型非常簡單，卻一直吸引著大量學(xué)者對其進(jìn)行研究，包括實(shí)驗(yàn)研究[1-3]和數(shù)值研究[4-6]，但潰壩問題仍然是大變形自由液面問題處理中的難點(diǎn)之一。

鑒于潰壩現(xiàn)象的特殊性，無網(wǎng)格方法越來越被廣泛地應(yīng)用來處理這類具有大變形自由液面的流動問題，MPS方法就是其中之一。MPS方法最早是由日本東京大學(xué)的KOSHIZUKA S，et al[7]于1995年提出的，相對于其它無網(wǎng)格方法，MPS方法采用預(yù)估-修正的半隱式方法來求解不可壓縮流體的控制方程，從而可以有效地避免離散對流項(xiàng)所引起的數(shù)值發(fā)散問題，另外，粒子具有運(yùn)動性并具有非常便捷的自由液面判別方法，這些都使得MPS方法在處理自由液面問題尤其是大變形問題中有其獨(dú)特的優(yōu)勢。但同時，由于該方法提出不久，MPS方法本身在穩(wěn)定性、精確性方面還存在一定的問題。因此，筆者對原始MPS方法進(jìn)行了改進(jìn)[8]，改進(jìn)后的MPS方法采用弧度判定法和粒子碰撞模型，從而可以有效地避免粒子聚集現(xiàn)象和粒子的誤判現(xiàn)象，大大提高了數(shù)值模擬的穩(wěn)定性和準(zhǔn)確性。本文采用改進(jìn)的MPS方法，通過模擬固定側(cè)壁的潰壩現(xiàn)象驗(yàn)證了該方法在潰壩問題中的有效性，探究了彈性支撐側(cè)壁對潰壩現(xiàn)象中的流動形態(tài)、側(cè)壁受力以及表面波對側(cè)壁的砰擊壓力的影響規(guī)律，旨在為解決實(shí)際工程問題提供一定的借鑒。

圖1 潰壩物理模型圖Fig.1 Physical model of the dam break

1 改進(jìn)的MPS方法數(shù)值模型

對于不可壓縮粘性流體，其控制方程為連續(xù)方程和N-S方程，如下所示：

其中，ρ為流體密度；v為流體質(zhì)點(diǎn)的速度；f為外部的作用力；ρ為流體質(zhì)點(diǎn)的壓力；υ為流體的運(yùn)動學(xué)粘性系數(shù)。

MPS方法是通過預(yù)估-修正求解不可壓縮流體的控制方程，具體實(shí)現(xiàn)方法如下。

1.1 粒子相互作用模型

在MPS方法中，采用可以自由移動的粒子來離散問題域，相鄰之間的相互作用是通過一個核函數(shù)來衡量，本文采用的是Koshizuka和Oka[9]在1996年提出的核函數(shù)，如下：

在該式中，r是相鄰兩個粒子之間的距離，re是核函數(shù)的控制半徑。某個粒子的相鄰粒子即在其核函數(shù)控制半徑內(nèi)的粒子。

粒子數(shù)密度表征的是流場中粒子的分布狀態(tài)，保持粒子數(shù)密度恒定就是保證流體的不可壓縮性。該參數(shù)可以通過粒子在核函數(shù)的控制域內(nèi)與其周圍相鄰粒子的核函數(shù)數(shù)值的疊加表示，即：

式中，N是粒子i的相鄰粒子總數(shù)，ri和rj分別是粒子i和粒子j的位置矢量。

梯度算子和拉普拉斯算子可以分別表示為：

和

1.2 邊界條件

邊界條件可分為兩種類型，即自由表面邊界條件和固壁邊界條件。

對于自由表面邊界條件，首先要進(jìn)行自由表面粒子的判定，原始的判定方法是根據(jù)粒子數(shù)密度進(jìn)行判定，具體判定公式如下[9]：

式中，β是一個小于1的系數(shù)，一般取β=0.97。

上面根據(jù)粒子數(shù)判定自由液面粒子的方法常常導(dǎo)致誤判現(xiàn)象的發(fā)生，為了克服這個問題，改進(jìn)的MPS方法中采用了一種弧度判定方法，這種方法通過對某個粒子的相鄰粒子相對該粒子的弧度進(jìn)行疊加，如果疊加后的弧度大于等于2 π并且在（0，2 π）區(qū)間內(nèi)無間斷，則該粒子判定為流體內(nèi)部粒子，否則為自由表面粒子，這種方法已證明可以有效避免粒子的誤判現(xiàn)象[8]。

對于自由表面粒子，其壓力賦為大氣壓即p=patm=0。

至于固壁條件，首先，固壁用三層邊界粒子來表示，靠近流體的那一層為流體邊界粒子，參與壓力的計算；外層的兩層粒子稱為固體邊界粒子，不參與壓力計算。另外，為了防止內(nèi)部的流體粒子穿透固體邊界，當(dāng)流體粒子靠近邊界時，會收到固壁粒子的一個排斥力，具體計算如下[10]：

其中，r0為初始粒子間距；D=5 gh，h為水柱的初始高度；P1=12，P2=4。

1.3 計算流程

（1）根據(jù)第n個時間步的信息，通過粘性項(xiàng)和外力項(xiàng)顯式預(yù)估粒子的中間速度和中間位置信息，即：

（2）通過式（3）計算粒子的中間粒子數(shù)密度信息n*，然后隱式求解Poisson方程，得到pn+1：

（3）根據(jù)pn+1對中間速度和中間位置信息進(jìn)行修正，得到第n+1時間步的速度和位置信息：

（4）為了避免內(nèi)部流體粒子的聚集現(xiàn)象，引入了一種碰撞模型，這是改進(jìn)MPS方法的另外一個改進(jìn)的措施。即當(dāng)兩個粒子之間的間距滿足d＜βL0時，應(yīng)用的碰撞模型如下[8]：

2 固定側(cè)壁潰壩現(xiàn)象

選取的固定側(cè)壁潰壩模型如圖2所示，水箱的長度為3.22 m，高度為2 m，水柱初始高度0.6 m，位于水柱對面墻壁上的一個點(diǎn)P（3.22 m，0.16 m）作為壓力監(jiān)測點(diǎn)。離散模型的數(shù)據(jù)為：流體粒子7 200個，流體邊界粒子723個，固體邊界粒子2 196個，共計10 116個粒子。

圖2 潰壩模型及離散圖Fig.2 Model of the dam break and the corresponding dicretisizing model

圖3 潰壩過程及壓力分布情況Fig.3 The process of dam break and the pressure distribution

圖3給出了改進(jìn)MPS方法模擬的潰壩過程中的六個典型時刻的流動現(xiàn)象及對應(yīng)的壓力分布。從圖中可以看出，當(dāng)潰壩現(xiàn)象發(fā)生后，水柱在重力的作用下發(fā)生坍塌形成了一個自由表面波或稱為潰壩波，該表面波在慣性的作用下向右傳播，然后與右側(cè)壁面發(fā)生砰擊現(xiàn)象，之后出現(xiàn)了自由表面的翻卷、破碎和飛濺現(xiàn)象，之后接著反向向左傳播，能量逐漸減弱直至達(dá)到新的靜平衡狀態(tài)。

圖4（a）給出的是潰壩波前端位置隨時間的變化曲線，圖4（b）給出了P點(diǎn)的壓力隨時間的變化情況，這兩幅圖均與對應(yīng)的實(shí)驗(yàn)值[11]進(jìn)行了比較。從圖中可以看出，潰壩波大約在T=2.5附近到達(dá)右面的側(cè)壁并緊接著發(fā)生了砰擊作用，整個過程共發(fā)生了兩次比較明顯的砰擊現(xiàn)象。另外，改進(jìn)的MPS方法的計算結(jié)果較為平滑，與實(shí)驗(yàn)值吻合更好。但是，改進(jìn)的MPS方法的計算結(jié)果，無論是潰壩波前端位置的變化還是P點(diǎn)的壓力時歷值，相對于實(shí)驗(yàn)值存在一定的延時，這種情況可能是由于數(shù)值模擬中運(yùn)動學(xué)粘性系數(shù)的取值較大造成的。

圖4 潰壩波前端位置和P點(diǎn)壓力的時歷曲線Fig.4 Time history of the leading edge of collapsed water column wave and the pressure at the point of P

3 彈性支撐側(cè)壁潰壩現(xiàn)象

實(shí)際工程中，當(dāng)發(fā)生潰壩時，由于潰壩波的巨大砰擊力的作用，會導(dǎo)致一些建筑物的移動或旋轉(zhuǎn)。因此，本文選取如圖5所示的具有彈性支撐側(cè)壁的潰壩模型，水柱右面的墻壁為彈性支撐，在潰壩波的砰擊力作用下將圍繞彈性支撐點(diǎn)做旋轉(zhuǎn)運(yùn)動，其轉(zhuǎn)動方程為：

式中，I為慣性矩；M為轉(zhuǎn)動體的質(zhì)量；H為轉(zhuǎn)動體質(zhì)心到旋轉(zhuǎn)點(diǎn)的距離；K為彈性支撐的剛度系數(shù)；θ為彈性支撐側(cè)壁的轉(zhuǎn)角；Mf為流體壓力作用在彈性支撐側(cè)壁上的合力矩，可按照下式計算：

其中，pi為彈性支撐側(cè)壁粒子i處的壓力；N為彈性支撐側(cè)壁離散的粒子總數(shù)；l為粒子i距彈性支撐點(diǎn)的距離；dl可取為粒子初始分布間距。

具體尺寸選取如圖2所示的潰壩現(xiàn)象的物理模型，只是右側(cè)固壁變成了彈性支撐側(cè)壁。在整個數(shù)值模擬中運(yùn)動性粘性系數(shù)v取0.001 m/s2，計算時間步長dt為0.000 5 s。圖6給出了具有彈性支撐側(cè)壁的潰壩現(xiàn)象模擬結(jié)果，包括6個在整個潰壩現(xiàn)象中比較典型時刻點(diǎn)的圖像，分布對應(yīng)T=3.859 9、6.147 3、7.576 9、9.649 8、11.293 8和21.444 0。圖7顯示了在整個潰壩過程中彈性支撐側(cè)壁的旋轉(zhuǎn)角度θ隨時間的變化情況，圖8中給出了P點(diǎn)壓力隨時間的變化規(guī)律，圖9比較了右側(cè)側(cè)壁在固定和彈性支撐兩種情況下的所受流體合力矩隨時間的變化情況。

圖5 彈性支撐側(cè)壁的潰壩模型示意圖Fig.5 Illustration of dam break with elastic-supported rigid wall

圖6 具有彈性支撐側(cè)壁的潰壩現(xiàn)象Fig.6 Phenomenon of dam break with elastic-supported side wall

結(jié)合圖6-9可以看出，當(dāng)潰壩現(xiàn)象發(fā)生之后，水柱坍塌并形成一個潰壩波，在T=3.8附近到達(dá)右側(cè)的墻壁，右側(cè)側(cè)壁所受流體的合力矩從0開始增長，在該力矩的作用下，彈性支撐側(cè)壁開始逐漸向右側(cè)旋轉(zhuǎn)；緊接著潰壩波與右側(cè)側(cè)壁發(fā)生了比較激烈的抨擊現(xiàn)象，右側(cè)側(cè)壁所受流體的合力矩在T=8左右達(dá)到最大值，與此同時其轉(zhuǎn)角也達(dá)到最大值0.25左右；此后，潰壩波在左右兩個墻壁之間反復(fù)的傳播直至靜止?fàn)顟B(tài)，右側(cè)側(cè)壁的旋轉(zhuǎn)角度也在這個往復(fù)震蕩過程中逐漸減小至0。

圖7 轉(zhuǎn)角θ隨時間的變化曲線Fig.7 The time history of the rotation angle

圖8 P點(diǎn)的壓力隨時間的變化曲線Fig.8 The time history of the pressure at the point of P

另外，結(jié)合前節(jié)固定側(cè)壁的潰壩現(xiàn)象，彈性支撐側(cè)壁的潰壩模型主要有以下特征：首先，具有彈性支撐側(cè)壁的潰壩現(xiàn)象從潰壩發(fā)生到重新恢復(fù)平衡狀態(tài)所消耗的時間要大于固定側(cè)壁所需的時間，前者是在T大于20以后才達(dá)到靜平衡狀態(tài)，后者是在T=11左右。與固定側(cè)壁的潰壩模型相比，彈性支撐側(cè)壁的潰壩模型的整個潰壩過程和側(cè)壁上流體壓力合力矩的峰值的出現(xiàn)均具有一定的延遲性。其次，這兩種潰壩模型右側(cè)側(cè)壁的壓力時歷曲線有明顯的不同。雖然這兩種潰壩模型的最大抨擊壓力峰值均在1.2左右，但固定側(cè)壁模型中有兩次比較明顯的砰擊現(xiàn)象，而對于彈性支撐側(cè)壁的潰壩模型，當(dāng)潰壩波到達(dá)右側(cè)側(cè)壁發(fā)生砰擊現(xiàn)象時，右側(cè)墻壁在砰擊力矩的作用下向右側(cè)旋轉(zhuǎn)，彈性支撐起到了明顯的延遲用，從而沒有明顯的第一次砰擊現(xiàn)象。

圖9 右側(cè)側(cè)壁的流體合力矩隨時間的變化曲線Fig.9 The fluid force momentum on the right-side wall

4 小結(jié)

本文基于改進(jìn)的MPS方法，采用Fortran語言進(jìn)行自主編程，分別模擬了固定側(cè)壁潰壩現(xiàn)象和彈性支撐側(cè)壁潰壩現(xiàn)象。研究證明，改進(jìn)的MPS方法可以成功地模擬潰壩現(xiàn)象中自由液面的翻卷、破碎和飛濺現(xiàn)象，驗(yàn)證了該方法在處理這種大變形自由液面流動問題中的優(yōu)勢。另一方面，通過這兩種潰壩模型的對比發(fā)現(xiàn)，雖然兩種模型的最大砰擊壓力相近，但固定側(cè)壁的潰壩模型有兩次明顯的砰擊現(xiàn)象，彈性支撐側(cè)壁的潰壩現(xiàn)象只有一次明顯的砰擊現(xiàn)象。其次，由于彈性支撐側(cè)壁的緩沖作用，整個潰壩過程的進(jìn)行出現(xiàn)了顯著的延遲效應(yīng)。