多物質(zhì)ALE方法在海嘯波沖擊方柱問題中的研究

翟 秋 ，張書劍，王華坤，張稼昊，嚴思寒

(河海大學 港口海岸與近海工程學院，南京210098)

海嘯是一種破壞性極強的自然災害。一些海嘯災害調(diào)查報告[1-2]顯示，海嘯襲擊后受災地區(qū)建筑物的損害程度差異較大，保存完好的建筑物往往位于大型建筑物的正陸側(cè)，損壞嚴重的建筑物通常位于街道的盡頭。以上兩種現(xiàn)象被稱為屏蔽效應(yīng)和聚焦效應(yīng)，它們在濱海城市規(guī)劃、海嘯防災減災等方面有重要啟示作用。

目前大部分物理模型試驗側(cè)重于海嘯波與單個樁柱或墻體構(gòu)件相互作用研究，只有少量工作涉及樁柱群和復雜建筑物的海嘯波水力特性分析。Hayatdavoodi等[3]對比研究了平板與 T梁在海嘯作用下的受力，探討了水深、波高、淹沒深度和抬高高度變化對波浪力的影響，發(fā)現(xiàn)平板上的受力一般是線性分布；楊萬理等[4]研究了門窗及屋面板洞口對低矮房屋海嘯作用力的影響，發(fā)現(xiàn)開洞率越大海嘯水平力越小；荀東亮等[5]研究了整流板、帶切縫的整流板和帶空氣孔的整流板等3種抗海嘯措施，結(jié)果表明整流板和帶切縫的整流板能減小海嘯水平力，帶空氣孔的整流板能有效減小豎向力；陳橙[6]通過模型實驗考察了潰壩波對由平板-樁-斜坡組成的簡易高樁碼頭的沖擊過程，分析了面板壓強隨海嘯波高和斜坡角度等因素的變化規(guī)律。

由于模型實驗成本較高，學者們更多用數(shù)值模擬手段來開展研究工作。景旭斌等[7]采用 ALE方法，分析了漂浮物在海嘯作用下對陸上建筑的作用力，發(fā)現(xiàn)漂浮物的質(zhì)量、撞擊速度對作用力起決定性作用，而被撞建筑的剛度對作用力的影響并不大；楊志瑩等[8]用數(shù)值模擬的方法，分析了海嘯和颶風作用下波高及淹沒系數(shù)對橋梁中T梁、箱梁受力的影響，發(fā)現(xiàn)相比于颶風波，海嘯對橋梁主梁的作用力更大，對橋梁安全性威脅更大。Pringgana[9]采用SPH方法研究了方柱自身旋轉(zhuǎn)角度對屏蔽效應(yīng)及聚焦效應(yīng)的影響，發(fā)現(xiàn)方柱正對海嘯波時受到的波浪力最大；Yang等[10]利用ANSYS Fluent軟件，研究了潰壩波沖擊下前屋對后屋的防護作用，發(fā)現(xiàn)如果前后屋之間的間隙足夠大，則前屋對后屋幾乎沒有防護作用；Wei等[11]用GPUSPH的方法，研究了海嘯沖擊下副公路橋及防波堤對主橋的沖擊減緩作用，發(fā)現(xiàn)防波堤與主橋之間的最佳減緩距離約為當?shù)厮畹?倍或來襲海嘯高度的13倍。

海嘯波與結(jié)構(gòu)物相互作用數(shù)值模擬涉及到流固耦合和自由液面問題，應(yīng)用較多的有SPH方法[9，11]和VOF方法[12-13]等。然而，SPH方法有著計算量大、耗費機時、自由表面計算精度較低、張力不穩(wěn)定等問題[14]；VOF方法根據(jù)體積比函數(shù)F來構(gòu)造和追蹤自由面，在處理F的變化時稍顯繁瑣[15]。相較而言，多物質(zhì)ALE法在流固耦合及自由液面問題上具有計算效率高、穩(wěn)定性好等優(yōu)勢。景旭斌[16]曾采用多物質(zhì)ALE法對孤立波夾雜漂浮物撞擊防波堤的過程進行了數(shù)值模擬，體現(xiàn)出該方法在流體沖擊問題中的適用性，然而其在沿海結(jié)構(gòu)物海嘯防災減災方面的應(yīng)用較少。

本文基于多物質(zhì)ALE方法建立了海嘯波與結(jié)構(gòu)物相互作用數(shù)值模型，研究了海嘯波對不同結(jié)構(gòu)布置形式下方柱的沖擊問題，以探明海嘯災害中的聚焦效應(yīng)及屏蔽效應(yīng)。

1 數(shù)值模型的構(gòu)建

流體運動及自由液面的處理采用多物質(zhì)ALE方法，而在固體域離散時使用拉格朗日描述，固體運動由線彈性結(jié)構(gòu)動力方程或剛體運動方程來控制，流體與固體之間的耦合采用罰函數(shù)法，基于上述算法構(gòu)建海嘯波與結(jié)構(gòu)物相互作用數(shù)值模型。

1.1 ALE描述流體運動控制方程

流體運動Navier-Stokes方程和連續(xù)方程的ALE形式可表述如下

(1)

(2)

式中：ξ為ALE坐標，xi和xj為空間坐標，t為時間，ρ為流體密度，vi為流體速度，ci是流體質(zhì)點相對于網(wǎng)格點的對流速度，σij為流體應(yīng)力張量，可用應(yīng)力偏量和壓力來表示

(3)

1.2 ALE多物質(zhì)材料模型

在多物質(zhì)ALE算法中，不同的材料能在計算區(qū)域的網(wǎng)格中自由流動，一個網(wǎng)格可以存在兩種及以上的材料，在計算過程中跟蹤每種材料的邊界，并在相應(yīng)的單元中進行物質(zhì)交換及輸送，在處理多相流與結(jié)構(gòu)之間相互作用的問題上表現(xiàn)較好。本文涉及水體及空氣兩種流體材料，水的狀態(tài)方程采用 Gruneisen方程

(4)

式中：p為壓力；E為單位體積內(nèi)能；C為us-up曲線截距；μ為體積變化率；S1，S2和S3為us-up曲線斜率系數(shù)；γ0為Gruneisen系數(shù)；α為對γ0的一階體積修正系數(shù)。

空氣的狀態(tài)方程則采用線性多項式狀態(tài)方程

p=C0+C1μ+C2μ2+C3μ3+(C4+C5μ+C6μ2)E0

(5)

式中：E0為初始比內(nèi)能，C0～C6為自定義常數(shù)。相關(guān)參數(shù)均可參考景旭斌[7]。

1.3 流固耦合處理方法

采用罰函數(shù)法處理流體與固體之間的耦合接觸。流體與結(jié)構(gòu)體相互接觸時，設(shè)流體為主面，結(jié)構(gòu)體為從面。每一個時步先檢查各從節(jié)點是否穿透主表面，沒有穿透則不對該從節(jié)點做任何處理；如果穿透，在該從節(jié)點與主表面間、主節(jié)點與從表面間引入大小與穿透深度和接觸剛度成正比的截面接觸力，其物理意義相當于在其中放置一系列限制穿透的法向彈簧。

2 數(shù)值模型的驗證

2.1 潰壩波沖擊彈性板算例

根據(jù)Rafiee和Thiagarajan[17]提供的算例建立數(shù)值水槽，如圖1所示。整個水槽由蓄水體、空氣域、壁面三部分構(gòu)成，在空氣域中放置一彈性板，其底端固定在床面。彈性板采用拉格朗日單元，材料為線彈性模型；蓄水體及空氣域采用多物質(zhì)ALE單元；水槽壁面按滑移邊界處理；對整個系統(tǒng)施加重力。

具體參數(shù)如下：蓄水體長L取14.6 cm，高2L，整個水槽長4L；彈性板與水體相距L，厚b取1.2 cm，高20/3b，密度ρ為2 500 kg/m3，楊氏模量E為106N/m2，重力加速度g取9.8 m/s2。當水體右側(cè)閘門被迅速抽離之后，生成潰壩波沖擊彈性板。

圖1 潰壩波沖擊彈性板數(shù)值模型 圖2 潰壩波沖擊彈性板的SPH與多物質(zhì)ALE法結(jié)果對比Fig.1 Numerical model of elastic plate impacted by dam break Fig.2 Comparison of SPH and multi-material ALE simulation for breaking dam on elastic plate

圖3 彈性板自由端位移歷時曲線比較Fig.3 Comparison between time histories of displacement of the free end of the elastic plate

圖2給出了潰壩波自由液面的歷時變化。圖左為文獻[17]SPH結(jié)果，右為本文結(jié)果，可見兩者的潰壩波液面形態(tài)吻合很好。將彈性板自由端的位移歷時曲線與文獻[17-20]結(jié)果進行比較，見圖3。可以看出，不同方法得到的位移在0.25 s達到最大值，之后振動相位互有偏差，但大體趨勢基本一致，即在振動幅值逐漸衰減波動。與各文獻數(shù)據(jù)相比，本文結(jié)果處于合理的范圍內(nèi)，證明了本文數(shù)值模型的可靠性。

2.2 潰壩波沖擊立柱算例

利用Gomez-Gesteira等[21]和Crespo等[22]提供的模型實驗數(shù)據(jù)進一步驗證本文算法，建立的水槽模型如圖4所示。數(shù)值水槽構(gòu)建方法同上一小節(jié)，閘門左側(cè)蓄水體高30 cm，右側(cè)床面以上設(shè)置高1 cm的淺水層，在空氣域中放置一剛性方柱，高75 cm，底端全約束。

潰壩波沖擊方柱的自由液面演變?nèi)鐖D5所示。圖6比較了柱體所受的潰壩波沖擊力。可以看出，本文數(shù)值結(jié)果與實驗數(shù)據(jù)接近，柱體所受沖擊力計算值為34 N，實驗結(jié)果為32.5 N。但在數(shù)值模擬中潰壩波沖擊柱體時刻比實驗中提前約0.07 s。這是因為本數(shù)值模型將水槽底面作滑移邊界處理，而實驗中存在底床摩擦，雖在閘門右側(cè)區(qū)域設(shè)置了淺薄水層以減輕水槽摩阻的影響，但無法將其完全消除，實驗中摩擦的存在消耗了部分潰壩水體動能，從而導致潰壩波沖擊立柱時刻稍晚于數(shù)值模擬，且峰值沖擊力也略有降低。總體上來看，基于多物質(zhì)ALE方法構(gòu)建的數(shù)值模型能較為準確地模擬結(jié)構(gòu)物受到的潰壩波沖擊力。

圖4 潰壩波沖擊剛體柱數(shù)值模型 圖5 潰壩波沖擊剛體柱的各時刻狀態(tài)Fig. 4 Numerical model of rigid cylinder impacted by breaking dam Fig.5 Time evolution of dam break wave impacting a rigid cylinder

圖6 潰壩波沖擊力數(shù)值模擬與實驗數(shù)據(jù)對比Fig.6 Comparison between numerical and experimental results of tsunami impact force

3 海嘯波沖擊不同排列下的多柱體研究

3.1 雙柱順排布置

采用多物質(zhì)ALE法對海嘯波沖擊雙柱順排布置下的剛性方柱進行探究，數(shù)值模型如圖7所示。以潰壩波來模擬海嘯波，前柱(方柱1)距潰壩水體保持0.5 m不變。后柱(方柱2)與后水槽壁保持0.56 m的距離不變，后柱與前柱之間的距離L分別取0.2 m、0.5 m、0.8 m、1 m四種工況，即無量綱間距比L/d(d為方柱邊長)分別為1.67、4.17、6.67、8.33。工況1～4及工況5(單柱)如圖8所示。

圖7 雙柱順排布置數(shù)值模型Fig.7 Numerical model for two cylinders in tandem arrangement

方柱1受到的海嘯波沖擊力如圖8所示，五種工況下的沖擊力幾乎相同，最大值均為24 N左右。各工況下的海嘯波沖擊力歷時曲線差別不大，這表明方柱1受后柱影響很小。

方柱2和單柱受到的波浪沖擊力如圖9所示，四種工況下方柱2受到的波浪沖擊力最大值分別為10.5 N、21.8 N、24 N、24 N；單柱受到的最大沖擊力為25 N。工況1下的最大沖擊力比其它工況下的減小了50%以上，當L/d由1.67增加到4.17時，最大沖擊力由單柱的42%上升到單柱的87%；當L/d繼續(xù)增加到6.67時，最大沖擊力與單柱的最大沖擊力25 N幾乎相同，即屏蔽效應(yīng)在L/d為6.67時基本消失。

圖8 方柱1受到的沖擊力 圖9 方柱2受到的沖擊力Fig.8 Impact force on cylinder 1 Fig.9 Impact force on cylinder 2

不同工況下方柱2受力最大時刻水體自由表面形態(tài)如圖10所示。在工況1及工況2中，因受到方柱1的屏蔽作用，水體繞過方柱1向后傳播，前后柱之間有一明顯的低洼水域，沖擊在方柱2上的水體動量比方柱1要小。在工況4中，沖擊在方柱2上的水體動量與單柱相比幾乎相同，前后柱之間觀察不到明顯的低洼區(qū)域。

10-a 單柱10-b 工況110-c 工況210-d 工況4圖10 不同工況下方柱2受力最大時刻水體自由表面形態(tài)Fig.10 Free surface of water corresponding to the maximum force on cylinder 2

可見，在雙柱順排布置的形式下，兩柱的間距越小，前面方柱1的屏蔽作用越明顯，方柱2受到的最大沖擊力明顯減小；當兩柱的間距超過一定范圍后，將不再產(chǎn)生屏蔽效應(yīng)。這種現(xiàn)象是因為海嘯波沖擊方柱1后受到阻礙，在其后方一定區(qū)域內(nèi)形成遮蔽區(qū)域，此區(qū)域內(nèi)波速、波高都顯著減小；在遠離方柱1一定距離之后波浪能量有所恢復，此時前后兩柱受到的最大沖擊力幾乎相同。

3.2 三柱錯排布置

對海嘯波沖擊三柱錯排布置下的方柱結(jié)構(gòu)物進行探究，數(shù)值模型如圖11所示。工況1只設(shè)置1根方柱，工況2、工況3、工況4均設(shè)置3根方柱，方柱的底面中心點連接起來為一個正三角形，邊長L分別為0.5 m、0.7 m、0.9 m，即L/d分別為2.5、3.5、4.5。對各工況下方柱1和方柱3受到的海嘯波沖擊力進行分析，由于方柱1和方柱2對稱布置，此處不再給出方柱2的受力結(jié)果。

圖11 三柱錯排布置數(shù)值模型Fig.11 Numerical model for three cylinders in staggered arrangement

工況1～4下方柱1受到的沖擊力如圖12所示。可見，各工況峰值沖擊力基本相同，均在106 N左右，各沖擊力歷時曲線也非常接近。工況2～4的結(jié)果顯示，隨著L/d的增加，方柱1的受力變化較小，這說明方柱1受到的最大沖擊力受周圍柱體干擾很小。

四種工況下方柱3受到的沖擊力如圖13所示，沖擊力最大值分別為106 N、112 N、120 N、106 N。工況2和工況3中方柱3的沖擊力最大值比單柱分別增大了5.6%和13.2%，工況4中方柱3的沖擊力最大值與單柱相近。基本規(guī)律是：隨著L/d的增加，方柱3的受力先增大再減小，最后和單柱接近。當L/d在一定范圍內(nèi)時，方柱1、2的存在會使方柱3受到的峰值沖擊力增大。這是因為水流在方柱1、2之間形成很強的間隙流，導致沖擊方柱3的水體動量增加。當L/d超過一定值時，間隙流逐漸消失，方柱3的受力接近單柱。

圖12 方柱1受到的沖擊力 圖13 方柱3受到的沖擊力Fig.12 Impact force on cylinder 1 Fig.13 Impact force on cylinder 3

不同工況下方柱3受力最大時刻潰壩波態(tài)如圖14所示。在工況2和工況3下，能觀察到方柱1、2之間較強的間隙流，但在工況4下方柱1、2間隙較大，水流的聚集效果已不明顯。

14-a 工況114-b 工況214-c 工況314-d 工況4圖14 不同工況下方柱3受力最大時刻水體自由表面形態(tài)Fig.14 Free surface of water corresponding to the maximum force on cylinder 3

4 結(jié)論

本文基于多物質(zhì)ALE法，綜合考慮了空氣、水體與結(jié)構(gòu)的耦合作用，對海嘯波沖擊彈性板及剛體柱兩個經(jīng)典算例進行了驗證，證明了該方法不僅能夠很好地模擬自由液面，而且能準確地計算水流沖擊力和結(jié)構(gòu)動力響應(yīng)。在此基礎(chǔ)上，應(yīng)用該方法探究了海嘯波對雙柱順排布置及三柱錯排布置形式下方柱的沖擊力規(guī)律。結(jié)果表明，前柱受力與結(jié)構(gòu)布置形式無關(guān)，當L/d較小時(L/d=1.67)，雙柱順排布置形式下后柱受到的波浪沖擊力相比單柱能降低50%以上，隨著L/d增大，后柱受力逐漸增大，且最終接近于單柱受到的沖擊力(L/d≥6.67)，此時意味著屏蔽效應(yīng)基本消失；三柱錯排布置形式下后柱受到的波浪沖擊力相比單柱會增大，且隨著L/d的增加呈現(xiàn)先增大再減小的趨勢，其受力相比單柱最多可增長10%以上(L/d=3.5)，此時聚焦效應(yīng)顯著。因此，濱海城市規(guī)劃及海工建筑物設(shè)計均應(yīng)考慮結(jié)構(gòu)物的布置形式，盡量利用順排布置下前方建筑物的遮擋效應(yīng)，避免錯排布置下聚焦水流的沖擊，以使海嘯災害減到最小。