波浪沖擊作用下核電站防浪堤動力響應的數(shù)值模擬

樓云鋒， 楊 勛， 王洪良， 葛鴻輝， 金先龍

(1.上海交通大學 機械與動力工程學院,上海 200240；2.上海交通大學 機械系統(tǒng)與振動國家重點試驗室，上海 200240；3.上海核工程設計研究院,上海 200233)

防浪堤結構廣泛用于沿海港口、沿海核電站等重要設施，日本3.11地震海嘯中釜石港防波堤有效阻擋了海嘯沖擊速度，減輕了大風浪對沿海設施造成的破壞。對防浪堤抗震、抗波浪沖擊及越浪問題已有研究，其中防浪堤抗大風浪沖擊研究由于沖擊過程含較復雜的流-固耦合瞬態(tài)沖擊問題，采用何種有效方法分析該問題成研究重點。

Isaacson等[1]理論研究波浪對水平圓柱的沖擊作用，獲得最大沖擊系數(shù)Cs=3.04～7.79，考慮結構動力響應等因素修正后為4.29。Cuomo等[2]為研究波浪對豎直墻沖擊作用力，通過物理模型試驗總結出豎墻表面波壓力計算經(jīng)驗公式。隨計算方法的改進與計算機性能的提高，數(shù)值模擬可能會發(fā)揮更大作用。Hajivalie等[3]采用RANS模型模擬二維直立式防浪堤受孤波沖擊及越浪過程。Hur等[4]通過VOF模型與多孔體模型相結合數(shù)值模擬研究三維潛堤的波浪力，結果與經(jīng)驗公式吻合。Losada等[5-6]采用COBRAS模型對斜坡堤沖擊越浪過程進行模擬，所得數(shù)據(jù)與實驗一致。而Zhang等[7]運用ALE方法分析船舶碰撞時流固耦合作用，認為ALE方法較好反映液體與固體間作用。Marco等[8]研究箱內(nèi)液面大幅晃動問題，將ALE方法與其它數(shù)值模擬方法進行比較，且與實驗結果吻合較好。通過簡單結構物模試驗研究總結出諸多規(guī)律及理論模型。但對復雜三維模型，理論方法顯然不適用。物模試驗往往既存在試驗繁雜又面臨數(shù)據(jù)檢測困難等問題。因此，選擇合適的流-固耦合計算方法，采用數(shù)值模擬結合物模試驗驗證參數(shù)為主要研究方向。

本文針對波浪沖擊的防浪堤動態(tài)響應進行數(shù)值模擬。通過動態(tài)接觸法[9]建立防浪堤整體結構模型，流體部分基于ALE方法建模[10]，并利用罰函數(shù)法[11]實現(xiàn)流固耦合。通過試驗驗證方法的可靠性，進一步用數(shù)值模擬方法分析波壓力分布、擋浪墻在沖擊過程中承受的最大水平推力及結構動態(tài)響應對波浪沖擊系數(shù)影響。

1 方法與原理

利用多物質ALE方法與罰函數(shù)方法相結合對防浪堤受波浪沖擊現(xiàn)象進行數(shù)值模擬。罰函數(shù)方法通過計算穿透速度與時間計算穿透量，再由穿透量計算耦合力，以保證耦合過程中能量守恒。ALE方法綜合Lagrange方法與Euler方法優(yōu)點，在材料域、空間域外引入?yún)⒖加颍⒃趨⒖季W(wǎng)格上通過算子分裂法求解，可克服Lagrange方法中網(wǎng)格大變形問題及Euler方法移動邊界問題。多物質ALE法采用VOF方法由計算網(wǎng)格單元中體積比函數(shù)F確定物質面，追蹤流體變化。

基于ALE方法的流體連續(xù)性方程、動量方程及網(wǎng)格更新過程為：

(1)

(2)

(3)

式中：χ為參考坐標；vi為流體流動速度；xi，xj為空間坐標；ci，cj為對流速度；ρ為流體密度；bi為流體力；σij為應力張量。

防浪堤結構采用彈性體連續(xù)方程：

(4)

式中：X為Lagrange坐標；ρs為容器密度；fi為體力；u為固體結構位移。

為適應顯式格式求解，采用弱可壓縮流體，引入線性牛頓流體本構方程：

σij=-pδij+2μsij

(5)

式中：p為流體靜水壓力；μ=μ(sij)為動力粘性系數(shù)，牛頓模型中為常數(shù)；sij為應變率張量；δij為克羅內(nèi)克常數(shù)。弱可壓縮流體求解需材料模型與狀態(tài)方程配合使用。

流固耦合需在液體、結構界面滿足幾何相容條件及力的平衡條件，即：

(6)

式中：vi,vs為流固耦合界面水流動速度、固體材料速度；Ff,Fs為流、固體結構作用在流固耦合界面的力，可由罰函數(shù)方法計算獲得：

(7)

2 模型

2.1 幾何設計圖與物理試驗模型

防浪堤抗波浪沖擊試驗斷面設計圖見圖1。防浪堤主要由素混泥土擋浪墻、堆石、護面結構堆砌而成。物理試驗模型見圖2，物模模型與實物間Froude數(shù)、Strouhal數(shù)相同?？紤]堤身高度、波浪要素、水深及水槽尺寸等因素，采用模型比尺M=25。

圖1 防浪堤設計圖

圖2 物模模型

2.2 有限元模型

基于核電站防浪堤設計圖紙，1∶1于物模模型建立波浪-防浪堤耦合系統(tǒng)數(shù)值模型。為準確描述防浪堤受波浪沖擊后的動力響應，防浪堤模型對實際試驗模型基本不簡化，各組件間通過動態(tài)接觸現(xiàn)實連接關系；結構Lagrange網(wǎng)格既滿足動態(tài)接觸所需尺寸匹配，又滿足與流體ALE網(wǎng)格耦合所需尺寸匹配。圖3為防浪堤結構模型，前后擋墻波壓力記錄點高程分別為A1∶0.44 m，A2∶0.4 m，A3∶0.36 m，A4∶0.32 m，B1∶0.498 m，B2∶0.472 m，B3:0.448 m，B4∶0.418 m。

圖3 防浪堤結構有限元模型

波浪沖擊耦合系統(tǒng)有限元模型及沖擊區(qū)域局部放大有限元模型分別見圖4、圖5。據(jù)對波浪傳播至防浪堤前距離估計建立10 m×1.2 m×1.7 m數(shù)值水槽模型，一端用速度入口，另一端模擬玻璃圍墻用法向約束，流場其余側面設無滑移邊界，防浪堤結構底部全約束，側面與后端面設置滑移邊界。考慮波浪傳播性質及計算精度，空氣與水體交界處網(wǎng)格劃分密集，流固耦合區(qū)單元劃分密集，其它區(qū)域稀疏，既能保證計算精度，亦能縮短計算時間。整個波浪沖擊系統(tǒng)有限元節(jié)點數(shù)1 187 140，單元數(shù)1 090 272。

圖4 整體數(shù)值模型

圖5 沖擊耦合區(qū)數(shù)值模型

2.3 材料模型

2.3.1 水體

本文水體用空材料模型描述應力偏量與應變偏量關系，描述弱體積變形與壓力關系狀態(tài)方程用Gruneisen方程：

其中：E0為單位體積內(nèi)能；C為us-up曲線截距(us為沖擊波波速，up為質點速度)；S1，S2，S3為us-up曲線斜率系數(shù)；γ0為Gruneisen系數(shù)；α為對γ0一階體積修正。

2.3.2 空氣

空氣采用空材料模型，狀態(tài)方程用壓力隨單位體積內(nèi)能線性變化，為體積變化多項式的Linear_Polynomial方程：

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

2.3.3 防浪堤結構

素混泥土采用線彈性模型(Elastic)。各組件材料見表1。

表1 防浪堤各組塊主要材料參數(shù)

3 試驗驗證及結果討論

本文用表2工況進行物理模型試驗及數(shù)值模擬：①通過物理模型試驗驗證本文所用建模方法與耦合方法的有效性；②通過數(shù)值模擬對波浪沖擊規(guī)律進行討論。用上海超級計算中心“蜂鳥”計算平臺及顯式動力學軟件LS-DYNA完成波浪對防浪堤沖擊動態(tài)響應計算。

3.1 數(shù)值模擬與試驗對比

建模時數(shù)值、物理模型在尺寸、結構布置及邊界條件已達一致。數(shù)值模擬與物理試驗均用規(guī)則波沖擊防浪堤結構，為保證數(shù)值模擬與物模試驗初始條件一致；波要素均用設計工況；通過波浪在堤上波浪形態(tài)與記錄各自預先設置在前、后擋浪墻8個波壓力記錄點的波浪沖擊壓力對數(shù)值模擬與物模試驗進行比較。對波浪上堤后傳播進化規(guī)律進行比較，包括平臺波浪變形、沖擊前側擋浪墻、波浪直接越過后擋浪墻及后擋浪墻二次越浪仿真數(shù)據(jù)與實驗室圖像對比，可定性認為數(shù)值模擬與物模試驗一致，見圖6。

表2 試驗工況

對前后擋浪墻表面波壓力進行比較。在物模前、后擋浪墻波壓力記錄點相同位置，數(shù)值仿真所得在規(guī)則波沖擊作用下前側擋浪墻表面波浪壓力曲線，減去靜水壓力后為波浪沖擊壓力，見圖7。

試驗數(shù)據(jù)見表3，數(shù)值仿真所得到A3，A4兩點靜水壓力為380 Pa，785 Pa，與理論水壓公式所得靜水壓一致；數(shù)值仿真所得前后擋浪墻8個波壓力記錄點最大波浪沖擊壓力較物模試驗略小，原因為數(shù)值計算中流體壓力由耦合單元計算獲得，而波壓力記錄單元尺寸與試驗中布置的波壓力傳感器大小有出入，導致計算數(shù)值誤差。本文通過控制耦合位置單元尺寸將誤差縮小在較小范圍，保證數(shù)值計算結果定量與試驗結果一致。

圖6 仿真數(shù)據(jù)、實驗室圖像間波浪演化定性比較

圖7 波壓力曲線

表3 波壓力峰值

3.2 數(shù)值模擬結果分析及討論

由波浪波高Hm、波周期Tm及防浪堤擋浪墻結構剛度要素對防浪堤抗沖擊動態(tài)響應進行數(shù)值模擬，獲得前、后擋浪墻表面波壓力分布規(guī)律、擋浪墻結構在波浪沖擊過程所受最大水平推力及波浪沖擊系數(shù)隨結構剛度變化情況。

3.2.1 擋浪墻表面波壓力分布

引入無量綱物理量z/d，Pmax/ρgHm，其中z為前擋浪墻記錄點相對前擋浪墻曲面底端高程，d為靜水面至前擋浪墻曲面底端水深。Hm為后擋浪墻記錄點至堤頂高程，Pmax為最大波壓力。

圖8 前后擋浪墻表面最大波壓力分布

設計工況前、后擋浪墻表面最大波壓力分布曲線見圖8。結合試驗數(shù)據(jù)知，前擋浪墻波壓力在靜水位處最大，靜水位以上呈線性減小，靜水位以下呈拋物線減??；由數(shù)值模擬、物理模型試驗的波浪在堤頂演化過程知，波浪在后擋浪墻形成二次越浪后大量水體沖擊后擋浪墻造成后墻底部波壓力最大，并沿高程方向減小。前、后擋浪墻表面波壓力隨波高、波周期變化見圖9。由圖9看出，前擋浪墻最大波壓力出現(xiàn)在靜水處，后擋浪墻最大波壓力位置在后墻底部，隨波浪強度的減小，波浪越過前擋浪墻后未直接沖擊后擋浪墻而在堤頂流動后沖擊后擋浪墻底部，后墻最大波壓力逐漸減小且波壓力分布集中于底部。

圖9 前、后擋浪墻表面最大波壓力隨波浪要素變化

3.2.2 擋浪墻承載能力

設計工況前后擋浪墻承所受最大一階主應力變化見圖10。由圖10看出，前、后擋浪墻先后達到最大值，且后擋浪墻出現(xiàn)兩次峰值，與物模試驗沖擊先后順序及后擋浪墻二次越浪現(xiàn)象一致。其中后擋浪墻應力最大值幾乎為前擋浪墻的兩倍。

沖擊中單位寬度擋浪墻水平推力計算式為：

(8)

式中：Pk為擋浪墻迎浪面波壓力記錄點值；ΔZ=0.04 m為擋浪墻迎浪面波壓力記錄點間距離；n為波壓力記錄點數(shù)，計算前擋浪墻水平推力時，含A1、A2兩點，計算后擋浪墻推力時，含B1、B2、B3、B4四點。

通過模型比尺等效后獲得現(xiàn)實工況的前后擋浪墻最大水平推力隨波浪周期、波浪高度變化見圖11。前、后擋浪墻水平推力均隨波浪強度下降大幅下降，且因后擋浪墻迎波面積大于前擋浪墻，故所受水平推力高于前擋浪墻。

3.2.3 波浪沖擊系數(shù)

據(jù)Cuomo[2]對波浪沖擊與準靜載的研究，引入無量綱物理量Fh/ρgHmd，Z/d為沖擊與準靜載,見圖12。

式中:d為靜水面至前擋浪墻曲面底端水深。

為研究防浪堤結構響應對波浪沖擊影響，本文設三種材料：剛體、C20混凝土(E=2.0E10)、C60混凝土(E=3.6E10)。波浪沖擊系數(shù)隨擋浪墻材料剛度及波浪條件變化見圖13。由圖13 看出，波浪沖擊系數(shù)隨結構剛度的減小而增大，隨波浪強度的增大而增大；且小剛度結構，沖擊系數(shù)隨水文條件變化更明顯。

圖10 前、后擋浪墻動態(tài)響應

圖12 沖擊與準靜載

圖13 不同結構剛度下沖擊系數(shù)比較

4 結 論

在對某沿海核電站防浪堤進行物理模型試驗基礎上，采用基于多物質ALE流固耦合方法與三維有限元模型對防浪堤沖擊響應進行研究。結合試驗數(shù)據(jù)，結論如下：

(1) 前擋浪墻波浪沖擊壓力最大處始終為靜水面位置，向上呈線性減小，向下呈拋物線減小；后擋浪墻波壓力最大位置在后墻底部，隨波浪強度的減小，后墻最大波壓力減小且較大波壓力集中于后墻底部；

(2) 通過數(shù)值計算獲得不同水文條件前、后擋浪墻承受的在最大波浪沖擊水平力及不同結構剛度的擋浪墻波浪沖擊系數(shù)。設計工況的后墻結構應力最大值約為前墻的兩倍；后擋浪墻始終較前擋浪墻承受的水平力大；結構動態(tài)響應對波浪沖擊系數(shù)有放大效應，且小剛度結構，沖擊系數(shù)隨水文條件變化更明顯。所得數(shù)據(jù)可為防浪堤結構安全設計及控制越浪設計提供依據(jù)。

參 考 文 獻

[1]Isaacson M,Prasad S. Wave slamming on a horizontal circular cylinder[C]. Proc. Civil Engineering in the Oceans V,1992.

[2]Cuomo G, Allsop W, Bruce T. Breaking wave loads at vertical seawalls and breakwaters[J]. Coastal Engineering, 2010,57(4):424-439.

[3]Hajivalie F, Yeganeh-Bakhtiary A. Numerical simulation of the interaction of a broken wave and a vertical breakwater[J]. International Journal of Civil Engineering, 2011,9(1):71-79.

[4]Hur D S, Mizutani N. Numerical estimation of the wave forces acting on a three dimensional body on submerged breakwater[J]. Coastal Engineering, 2003, 47(3): 329-345.

[5]Losada I J, Lara J L, Guanche R. Numerical analysis of wave overtopping of rubble mound breakwaters[J]. Coastal Engineering, 2008, 55(1): 47-62.

[6]Hsiao S C, Lin T C. Tsunami-like solitary waves impinging and overtopping an impermeable seawall: experiment and RANS modeling[J]. Coastal Engineering, 2010, 57(1):1-18.

[7]Zhang A, Suzuki k. A comparative study of numerical simulations for fluid-structure interaction of liquid-filled tank during ship collision[J]. Ocean Engineering, 2007, 34(5-6): 645-652.

[8]Marco A, Castelletti L M L. Fluid-structure interaction of water filled tanks during the impact with the ground[J]. International Journal of Impact Engineering, 2005, 31(3): 235-254.

[9]郭毅之，金先龍，丁峻宏．沉管隧道地震響應分析中的三維接觸模型與算法研究[J]．應用力學學報，2006，23(1):48-52．

GUO Yi-zhi, JIN Xian-long, DING Jun-hong. Three dimensional contact model and algorithm for earthquake response of immersed tunnel[J]. Chinese Journal of Applied Mechanics, 2006, 23 (1):48-52.

[10]張 雄,陸明萬,王建軍. 任意拉格朗日-歐拉描述法研究進展[J].計算力學學報,1997,14(1):91- 102.

ZHANG Xiong, LU Ming-wan, WANG Jian-jun. Research progress in arbitrary Lagrangian-Eulerian method[J]. Chinese Journal of Computatlonal Mechanics,1997, 14(1):91-102.

[11]張曉敏，金先龍，陳向東. 風荷載作用下大型柔性儲液結構動態(tài)響應仿真方法研究[J].振動與沖擊,2009, 28(5):115-118．

ZHANG Xiao-min, JIN Xian-long, CHEN Xiang-dong. Numerical simulation of dynamic characteristics of flexible container under wind load[J]. Journal of Vibration and Shock,2009, 28(5): 115-118．