淺海風-浪-流-海床耦合場非定常時空演化規律及評價指標

陳 靜, 柯世堂, 李文杰, 朱庭瑞, 員亦雯, 任賀賀

(南京航空航天大學 土木與機場工程系; 江蘇省機場基礎設施安全工程研究中心, 南京 211106)

隨著我國海洋強國戰略[1]的快速推進,近海岸海洋工程建設需求與日俱增,南海海域經歷復雜的地質演化過程形成不同類型的海床地貌,其三維效應對風-浪-流耦合場演化規律具有顯著影響,導致近岸風、浪、流場的時空分布及統計特性不同于開闊海域,準確預測復雜海洋環境是海岸工程安全保障的重要前提.目前國內外規范[2-3]關于淺海區波浪的計算方法是基于線性波浪理論假定,不能考慮海床面粗糙度對波流運動的影響,工程設計中海洋環境預測結果仍只能過度依賴實測數據校準和關鍵參數取值經驗.隨著計算機技術高速發展,計算流體力學(CFD)數值模擬現在已成為風-浪-流環境耦合場預測的重要手段,但目前模擬采用的邊界造波法忽略了初始場未形成已改變風浪流要素的問題,一定程度限制了淺海地區風-浪-流耦合場的預測精度.

現有關于風-浪-流耦合場研究主要集中在數值方法改進[4-6]、流場演化規律[7-10]、結構水動力模型[11-14]等方面,文獻[9]中基于Fluent平臺二次開發研究發現風-浪-流三者耦合作用下,風會導致波浪傳播過程中波長增大、波幅減小;文獻[6]中結合理論推導將波浪水質點速度化作單位矢量代入風剖面公式之中,改進了風-浪耦合邊界自由液面處風與水質點速度方向的一致性.但上述模擬過程中風-浪-流耦合模型通常假定風速、流速沿高度均勻分布和水深不變,無法考慮淺海地區梯度風、剪切流的垂向速度不均勻和海床地貌不規則對耦合場的影響.針對海床地貌的影響研究,僅有部分學者對淺水波浪的非線性變化[15]或斜坡地形的波浪破碎[16-21]問題進行探索,文獻[19]中結合物理模型實驗結果分析了礁前坡度、坡頂淹沒水深與規則波在島礁地形上的破碎指標之間的相關性;文獻[20]中指出隨海床面坡度的增大,孤立波在傳播過程中受地形因素導致的能量耗散作用為主,淺化作用為輔;文獻[21]中為研究孤立波與斜坡作用,利用新一代CFD軟件STAR CCM+以造波邊界設置多層流體速度方式構建二維內孤立波數值水槽,模擬波浪振幅與設計振幅誤差在0.5%以內.此類研究大都局限于斜坡海床下波浪場的破碎現象,未涉及真實海域典型海床地貌與風-浪-流耦合場的非定常演變物理過程研究,且缺乏基于不同海床對風-浪-流場非定常演變特征的評價標準.

鑒于此,基于STAR-CCM+二次開發考慮了梯度風、剪切流的垂向速度不均勻對耦合場的影響,開展多層三維黏性風-浪-流-海床耦合場仿真模擬,分析4種典型海床地貌下波浪、海流和風場的非定常時空演變規律,對比分析了波高、波長、平均流速、流剖面指數、基本風速及風剖面指數的差異及產生原因,最后結合主成分分析法提出各海床工況下風-浪-流耦合場全生命周期非定常綜合評價指標,可為淺海地區典型海床地貌下的風-浪-流場精準預測提供參考.

1 風-浪-流-海床耦合場模擬方法

1.1 數值模擬方法

1.1.1典型海床概況 我國南海海域包含著紛繁復雜的地貌單元,整體發展趨勢從四周向中部水深逐漸變大,將陸相-海陸過渡相-海相的發展區域地貌單元依次劃分為大陸架、大陸坡、邊緣海盆地[22-23].研究區域位于南海海域的東北位置,屬于大陸架與大陸坡的交界區域,重點以大陸坡典型海床地貌作為研究對象,其中大陸坡的三級地貌共分為海底平原、海底斜坡和海槽,以3種地貌平均深度較低的四級地貌[23]:海底平原的深海水道、海底斜坡的褶皺型“微凹陷”和海槽地帶的“U” 形“微凹陷”作為典型地貌,再以無地貌的平坦海床作為對比工況,以此來研究典型海床地貌下波浪場、風場與海流場演化過程.圖1給出了海底平原、海底斜坡和海槽典型地貌的數值模型,床面高度指海床最大高度差,分別為18、5、10 m;褶皺角度(30°)指褶皺頂部角度,U彎角度(15°)指中央拱起角度.

圖1 南海東北部海域海床地貌及簡化模型圖Fig.1 Seabed topography and simplified model map in northeastern South China Sea

1.1.2黏性流體流動數學模型 通常假設浮體周圍流域內的流體是不可壓縮性的,不考慮流體的表面張力,可結合連續性方程和N-S方程來表達流體的黏性運動[24]:

(1)

(2)

式中:u為三維流體x、y、z方向的速度矢量;ρ為密度;t為流體運動時間;μ為動力黏度系數與渦動黏度系數之和;p為壓力;g為重力矢量;fs為消波源項.

海床面的運動邊界條件:

(3)

式中:D為海水深度.

1.1.3改進多層風-浪-流耦合模型 線性波浪波面高程方程為

(4)

式中:η為波面高程;H為波高;k為波數;ω為圓頻率.

以五階Stokes波模擬周期性波浪,該波浪波面具有非對稱余弦曲線函數狀態,可直觀觀察到波峰與波谷形態,比線性波浪更接近實際海況.其波面高程方程為

(5)

x方向的速度表達式:

ux=

(6)

z方向的速度表達式:

uz=

(7)

式中:d為水深;c=cosh(kd);λ1=λ,λ2=λ2B22+λ4B24,λ3=λ2B33+λ4B35,λ4=λ2B44,λ5=λ2B55;λ、B系列參數詳見文獻[25].

風場依據《建筑結構荷載規范》[26]采用A類地貌風剖面系數,以指數率梯度風進行定義:

(8)

式中:uw、u10分別為高度z和高度10 m對應的平均風速;α為地面粗糙度指數,A類地貌取值0.12.

自由液面以下的海流流速以1/7定律的指數型分布表示:

(9)

式中:uc、vc分別為海流的流速和平均流速;l為距海床底距離.

實際海域中水流速度在垂向上具有一定的分布差異,從而導致波浪沿不同水深發生變化;同時風在垂向上也呈現梯度分布,導致波面傳播會產生變形現象.為研究近海風-浪-流演化問題,需要同時考慮風、流垂向分布的影響.以多層質點速度耦合方法應用于風浪流共存情況,如圖2所示.圖中:在多層風浪流耦合模型中建立笛卡爾坐標系OXYZ;原點O為海平面位置;Z坐標描述水深;Y坐標描述風-浪-流橫向寬度;風-浪-流沿著X正方向傳播.對計算域進行風波耦合區、波流耦合區、水深一半以上區和水深一半以下區的分層,不同層的風速、流速以風剖面、剪切流的速度進行分解輸入,以此建立新型多層風-浪-流耦合模型.該模型的速度uwj、uci分別為j、i層的風浪耦合速度和浪流耦合速度.風-浪-流在x方向上的速度:

圖2 多層風-浪-流耦合場模型示意圖Fig.2 Diagram of multilayer wind-wave-current coupling field model

Ux=

(10)

1.2 CFD模擬參數設定

圖3給出了風-浪-流-海床三維數值黏性水池模型示意圖.計算域由風-浪-流生成區、阻尼消波區和風-浪-流-海床耦合演化區組成.為保證流體充分發展,計算流體域取400 m×200 m×330 m(流向x×展向y×縱向z),水深按照實際海深取值30 m,尾端設置尺寸為66 m×200 m×330 m的消波區.

圖3 多層風-浪-流-海床三維數值黏性水池模型示意圖Fig.3 3D numerical viscous pool model of multilayer wind-wave-current-seabed

計算域邊界條件設置如下:左邊界為風-浪-流速度入口,右邊界設為壓力出口,前后邊界設為對稱邊界,頂部邊界設為滑移壁面,底部邊界設為無滑移壁面,其中為了流體得以更快地形成穩定場,速度入口和壓力出口設置成周期性接觸,形成充分發展性邊界.雷諾平均N-S方程(RANS)非定常計算模型采用更適合于求解雷諾應力分量輸運方程的SSTk-ω湍流模型,使用流體體積方法對自由液面進行追蹤處理.通過自定義場函數接口,將多層風浪流耦合場模型施加在入口邊界上,初始化時停用第一波前設置,可生成對應的風浪流初始場.為防止出口風浪流反射,形成壁面效應,在水池尾部添加動量源阻尼項.根據DNV-RP-C205規范[3]中統計的各海區風浪流長期分布參數,選擇研究區域屬于40號海區,選取海區作業海況作為輸入條件,波高設置為2.1 m,周期為5.3 s,設10 m處風速為11.4 m/s,平均流速為0.68 m/s,通過波浪速度分解x、y和z向,與自定義風剖面、剪切流疊加計算.

在計算域全域沿高度定義波峰以上、波谷以下及波浪發展范圍內的風-浪-流解耦速度(式(10)),初始時刻計算域全域即實現了風-浪-流場解耦,解決現階段數值造波方法存在的初始場形成前,風、浪、流要素已改變問題,同時提高計算效率和精度.圖4給出初始化流場分布圖.

圖4 風-浪-流水池初始化速度場云圖Fig.4 Initial velocity field cloud image of wind-wave-flow tank

圖5給出了風-浪-流-海床耦合數值水池網格劃分示意圖.整體采用混合網格離散形式,以切割型六面體網格加密程度區分液面加密域、液面過渡域和其他區域.其中對液面加密域沿波浪傳播方向、波高方向加密網格,液面過渡區相對于液面加密區兩倍的網格尺寸進行加密.y+值為無量綱壁面距離,關乎近壁面湍流的發展狀態,近壁面區域使用全y+壁面處理函數方法,近壁面第1層網格厚度為0.041 m,y+范圍在150～250內,滿足壁面函數法要求.

圖5 風-浪-流-海床三維數值粘性水池網格劃分示意圖Fig.5 Diagram of wind-wave-stream-seabed three-dimensional numerical viscous pool grid

表1給出了不同水線加密工況下參數對比,其中L為波浪波長.由表可知:隨著網格數量的增加,波高誤差逐漸減小趨勢,波高衰減情況改善.工況4和工況5的波高誤差與規范結果相差較小,且兩者之間也無明顯差異,綜合考慮計算精度與效率的平衡,本文選取網格工況4為數值模擬網格方案,時間步長按照T/(2.4r)取值.其中:T為波浪周期;r為一個波長的網格數量.時間步長為0.022 s.

表1 網格劃分方案工況Tab.1 Working condition of grid partitioning scheme

1.3 有效性驗證

風-浪-流參數模擬精確度與試驗結果有很大的關聯性,為驗證數值水池模擬試驗結果,分別對波浪單獨作用的波浪場及消波效果、對初始時刻風-浪-流耦合作用下的風場及海流場進行分析,通過距入口、出口20 m處監測的波面高程、縱向分布風速及流速的模擬值與理論值進行對比,來驗證試驗模擬的準確性,其中風-浪-流規范值參數如表2所示.

表2 風-浪-流參數Tab.2 Wind-wave-current parameter

圖6(a)為波浪單獨作用時距入口20 m處波浪模擬值與理論值對比曲線,在模擬的60～90 s階段,波浪的模擬值與理論值波形基本吻合,且無衰減的趨勢;圖6(b)為距出口20 m消波區波面高程時程曲線,由圖可見相較于未消波和消波后的理論值,波浪在兩個周期以后已基本消失,驗證了阻尼消波的可行性.

圖6 波高模擬值與理論值對比時程曲線Fig.6 Comparison of time-history of wave height simulation value with theoretical value

圖7給出了風-浪-流耦合作用下初始時刻計算域距入口20 m處沿高度風速和流速模擬結果對比曲線.結果表明入口處風剖面和剪切流剖面與理論解吻合,初始條件下風場和海流場達到穩定.

圖7 距入口20 m處風速及流速剖面模擬值與理論值對比曲線Fig.7 Comparison of simulated and theoretical values of wind speed and current speed at 20 m from inlet

2 耦合場時空演化規律

2.1 波浪場時空演化

圖8給出了計算域中間(x=0 m)4種海床地貌下波浪場波形變化時程曲線.對比分析發現風-浪-流場與平坦地貌、海底斜坡和海槽的耦合作用使得波高大幅度減小、波長變長;平坦地貌、海底斜坡和海槽工況呈現衰減演化過程,海底平原工況呈現波浪破碎過程.平坦地貌、海底斜坡和海槽3種地貌下的衰減演化過程均可分為波面激增、衰減和穩定3個階段,其中衰減階段平坦地貌和海底斜坡工況波浪場處于持續衰減,而海槽地貌由于“淺-深-淺”的地形特征和色散定律,波浪場形成波群性衰減;海底平原工況先經歷外破波和內破波兩個階段的波浪破碎,呈現波高突變和波形衰減的非線性變化,最后在爬波階段受到觸底回流的沖擊,形成最后一次高達6 m的“直立”形波浪破碎現象.

圖8 不同海床地貌下波形高程-時程變化曲線Fig.8 Waveform time-history curves in different seabed topographies

圖9～12給出了演化階段下各海床工況全域空間波形變化示意圖.由圖可知:平坦地貌、海底斜坡和海槽3個工況經歷過激增階段的波面陡立后,都呈現不同程度的衰減現象,其中海槽地形也由于色散定理和水深特征,波浪的前緣速度和后緣速度不同,波面的橫向不對稱性相比其他工況較明顯;海底平原工況外破波階段波浪破碎形態多為激破波與卷破波組合的破碎形式,內破波階段形成出口端波浪折射現象和部分前沿陡立、后坡平坦的段波,與爬坡階段的觸底回流碰撞形成崩破波.

圖9 平坦地貌不同階段三維空間波形演變云圖Fig.9 Three-dimensional spatial waveform variation diagram under flat geomorphic conditions

圖10 海底平原不同階段三維空間波形變化圖Fig.10 Three-dimensional spatial waveform variation diagram under submarine plain working conditions

圖11 海底斜坡不同階段三維空間波形變化圖Fig.11 Three-dimensional spatial waveform variation diagram under submarine slope conditions

圖12 海槽不同階段三維空間波形變化圖Fig.12 Three-dimensional spatial waveform variation diagram under trough working conditions

2.2 海流場時空演化

圖13為各工況的剪切流模擬結果對比,其中每幅下圖對應上圖的最后時刻,β值為流剖面指數,初始為1/7.對比可見:與平坦地貌相反,海底平原、海底斜坡和海槽工況對流剖面要素均存在一定影響,其中平均流速皆呈現隨時間逐漸減小的趨勢;海底平原工況由于觸底回流的原因,演化終期垂向流速基本為負值.分別采用指數型、多項式對一半水深以下、一半水深以上和波流耦合區3個區域進行擬合,并給出如下4個工況擬合后的垂向流速解析式.

教師主要以課堂教學的知識和理論基礎為基礎，與直接獲得實際和實踐能力的生產和科研實踐相分離。因而，通過實踐教育基地的建設，給教師搭建了產學研合作的平臺，不僅擴寬了教師的視野，實現知識的轉化，而且還能提高教師的科研水平；同時，企業也獲得了人才、技術的有力支持，對于提高企業新技術、新產品的開發有了進一步的保障。

圖13 不同海床地貌下剪切流模擬結果對比Fig.13 Comparison of shear flow simulation results in different seabed geomorphologies

平坦地貌:

(11)

海底平原:

(12)

海底斜坡:

(13)

海槽:

(14)

圖14～17給出各工況下速度流線空間分布圖.分析可得:不同于傳統的風-浪-流耦合場波峰下部“V”形、波谷下部“U”形間隔渦區的研究結果,耦合場由于垂向壓差,形成了不同尺度的周期性渦旋;平坦地貌工況地形對渦旋發展無阻礙,形成“多渦積聚”,海底斜坡工況“褶皺”地形的阻礙導致流場形成“多渦共存”的現象;而海底平原工況由于卷破波形成時部分空氣卷入,加劇低空處風和波浪速度,波浪場加劇變化使得全域波形呈現“坡形”,風場與其碰撞形成圖15(b)中風波耦合區的渦旋,爬坡階段渦旋發展過程遇到觸底回流,海流場渦旋體系被破壞;對于海槽工況,海流場發展空間較小,但經過長時間的演化,形成的小渦旋逐漸破壞初始“U+V”渦區形態.

圖14 平坦地貌工況流場時空演化圖Fig.14 Space-time evolution diagram of flow field under flat landform conditions

圖15 海底平原工況流場時空演化圖Fig.15 Space-time evolution diagram of flow field in submarine plain

圖16 海底斜坡工況流場時空演化圖Fig.16 Temporal and spatial evolution diagram of flow field under submarine slope conditions

圖17 海槽工況流場時空演化圖Fig.17 Temporal and spatial evolution diagram of flow field in trough conditions

2.3 風場時空演化

圖18為不同工況下風剖面擬合結果對比,其中每幅下圖對應上圖的最后時刻.對比發現4種海床地貌引起的波浪場變化均會造成風剖面指數的增大和基本風速的減小,且參考表1,發現海床床面高度與風剖面指數成正比關系,與A類地貌取值0.12[26]相比,平坦地貌、海底平原、海底斜坡和海槽的演化終期風剖面指數分別為0.164、0.197、0.173和 0.176;海底平原工況由于外破波階段波高全域性突變至7～8 m,海表面的粗糙度大幅度提高,此時風場也在波面抬升中損失大量能量,導致圖18(b)中在t=7 s附近基本風速降到極小值5 m/s,風剖面指數提高到極大值0.4,在后面的演化階段,風場恢復至穩定,但其基本風速的穩定值小于其他3個工況.

圖18 不同海床地貌下風剖面模擬結果對比Fig.18 Comparison of simulation results of wind profile in different seabed topographies

3 耦合場非定常效應評價指標

(15)

結合主成分分析方法,根據風-浪-流演化過程中每0.02 s模擬出的6個非定常環境因子作為樣本,提取最大貢獻主成分.圖19分別為4個海床工況的非定常環境因子在主成分1×主成分2 上的投影.圖中:圓點的大小與演化時間成正比;箭頭終點的坐標值代表各非定常因子對主成分的貢獻值.表3為不同海床地貌下主成分最大貢獻數量及累計貢獻率.其中平坦地貌、海底斜坡及海槽工況前3個主成分的累計貢獻率分別達到91.7%、85.6%和94.6%;海底平原工況所需4個主成分,累計貢獻率達到90.0%,已滿足主成分分析法降維精度要求[27].

表3 不同工況主成分貢獻數量及累計貢獻率Tab.3 Number and cumulative contribution rate of principal components under different working conditions

圖19 非定常環境因子和樣本數據對主成分貢獻數值圖Fig.19 Unsteady environmental factors and sample data contribute to numerical diagram of principal components

據此可建立風-浪-流全生命周期非定常綜合評價函數:

F=

(16)

式中:X1～X6依次為波高、波長、平均流速、流剖面指數、基本風速與風剖面指數的非定常環境因子.基于此取演化終期的非定常環境因子數值,建立各海床地貌工況下風-浪-流耦合非定常評價指標參考值.得到平坦地貌、海底平原、海底斜坡和海槽的非定常評價指標分別為0.268、4.612、0.672和0.926.

4 結論

系統研究了典型海床地貌下風-浪-流耦合過程非定常演化規律及評價指標,尤其在多層風-浪-流耦合模型、風-浪-流場階段性發展形態演化和非定常綜合評價指標方面取得了原始創新,具體如下:

(1) 新型多層風-浪-流耦合模型能精確模擬出垂向風速、流速不均勻對風-浪-流耦合場造成的影響.以風-浪-流解耦模型運用于CFD初始條件中,有效提高計算效率和精度.

(2) 各海況的風-浪-流-海床演化過程:波浪場演變分為三階段,其中平坦地貌、海底斜坡和海槽工況分為波面激增、衰減和穩定階段;海底平原工況共分為外破波、內破波和爬坡階段.海流場隨演變過程在空間上呈多段式分布,一半水深以下的指數率函數、一半水深以上及波流耦合區的多次項函數.

(3) 海流速度場由于垂向壓差形成渦旋,各工況根據海床地貌不同形成多渦積聚或多渦共存的現象;風-浪-流-海床耦合演化對風剖面指數產生放大作用,海床床面高度與風剖面系數成正比關系.

(4) 建立了基于主成分分析法的耦合場全生命周期非定常評價指標, 并提出演化終期各典型海床下非定常評價指標數值:平坦地貌、海底平原、海底斜坡和海槽分別為0.268、4.612、0.672和0.926.非定常評價指標可精準預測淺海地區典型海床地貌下風、浪、流時空分布變異性特征, 可為我國海岸工程安全建設中海洋環境預測提供參考.