王 寅，沈火林，王玲玲，計 勇，席芝橙，張雯雯

(1.南昌工程學院 水利與生態(tài)工程學院，江西 南昌 330099；2.河海大學 水利水電學院，江蘇 南京 210098)

在穩(wěn)定分層的流體系統(tǒng)中，外界微小的擾動可在密度分界面(密度躍層)上引起攜巨大能量的內(nèi)波現(xiàn)象[1]。內(nèi)孤立波系波幅大、周期短、所攜能量大的一類非線性內(nèi)波[2]。該類波在傳播時所產(chǎn)生的作用力可對近岸及河口復雜水動力環(huán)境下水下支撐結構物的安全穩(wěn)定性造成嚴重的威脅[3]。

相關學者就柱體在內(nèi)孤立波環(huán)境下的受力問題開展了系列理論研究。蔡樹群等人運用Morison公式估算了內(nèi)孤立波對柱體的作用力，基于該半經(jīng)驗半理論公式求解出小尺度圓柱受力的垂向分布[4-5]。但在采用Morison公式時，形狀阻力系數(shù)CD和慣性力系數(shù)CM雖可通過相關的試驗結果進行選取，但此類試驗都在非內(nèi)孤立波環(huán)境下進行。內(nèi)孤立波所誘發(fā)的分層流環(huán)境下的流場與密度均一流環(huán)境差別較大[6]，兩個受力系數(shù)是否適用于內(nèi)孤立波對柱體的作用問題尚無定論；同樣采用Morsion公式計算內(nèi)孤立波對柱體的作用力，Du[7]得到內(nèi)孤立波作用效果大約與波長為300 m，波高為18 m的表面波所致荷載相等，而Song[8]發(fā)現(xiàn)內(nèi)孤立波作用效果遠弱于表面波，約為9%。可見，即使計算結果可定性計算出柱體所有內(nèi)孤立波的作用力，但在定量上還不夠準確。采用更精細的計算方法獲取內(nèi)孤立波對柱體的作用機制勢在必行。

本研究在構建三維數(shù)值內(nèi)波水槽的基礎上，采用大渦模擬(large-eddy simulation，LES)技術，模擬內(nèi)孤立波的產(chǎn)生和傳播，研究不同內(nèi)孤立波波幅(ηo)對柱體水動力荷載的影響。對比分析各工況下柱周水動力特性，提取相關水動力參數(shù)云圖及矢量圖，深入剖析柱體在內(nèi)孤立波環(huán)境下的受力規(guī)律。

1 數(shù)值方法

1.1 Navier-Stokes方程(N-S方程)

基于連續(xù)性假設，在三維瞬態(tài)流動時，用于描述不可壓縮粘性水體的Navier-Stokes(N-S方程)可表達為

(1)

(2)

式中ρ為密度項；t為時間；xi，xj為笛卡爾坐標系的3個坐標方向；ui，uj為笛卡爾坐標系中的3個流速分量；p為壓力項；μ為動力粘性系數(shù)；fi為i方向的單位體積力。

1.2 標量輸運方程

由于標量(溫度或濃度)而導致的流體分層現(xiàn)象，以密度躍層為界線，上、下層流體之間會進行質(zhì)量交換和對流—擴散，而標量輸運方程則和它們緊密相連，即：

(3)

ρ=C+ρ0，

(4)

式中C為鹽度，kg/m3；k為分子擴散系數(shù)；S為源項；ρ為水體密度，kg/m3；ρo為清水密度，kg/m3。

1.3 大渦模擬(LES)控制方程

大渦模擬技術是對湍流脈動的一種空間平均，經(jīng)濾波函數(shù)將大尺度渦和小尺度渦分離開來[9]。過濾后的動量和質(zhì)量方程可表示為

(5)

(6)

(7)

2 數(shù)學模型

2.1 三維模型的建立、造波方法及工況設置

本研究所建立的三維數(shù)值波浪水槽長(X)、寬(Z)、高(Y)為：12 m×0.6 m×0.8 m，見圖1。所研究的柱體位于槽寬(Z)中心，其底部的中心點坐標為：(x，y，z)=(6，0，0.3)，坐標原點為圖1中的左里下角。柱體直徑D=10 cm。

圖1 三維數(shù)值波浪水槽示意圖

定義該柱體在內(nèi)孤立波環(huán)境下所受到的無量綱水平合力CFn的計算公式為

(8)

式中Fn為數(shù)值模型計算所得的柱體水平合力(N)；g為重力加速度(cm/s2)；A為柱體的迎風面積(cm2)；H為水槽總水深(cm)。

本研究所制造出的內(nèi)孤立波波形為下凹型，數(shù)值模型中分層水體基本參數(shù)設定為：水槽總水深H=80 cm，其中上層水體密度ρ1=0.998(g·cm-3)，水層厚度h1=20 cm；下層水體密度ρ2=1.017(g·cm-3)，水層厚度h2=60 cm，數(shù)模整體布局如圖2所示。采用重力塌陷法制造內(nèi)孤立波[10]。為制造出下凹型內(nèi)孤立波，需把水槽在縱向分為造波區(qū)和傳播計算區(qū)兩部分，并使造波區(qū)密度躍層位置低于傳播計算區(qū)密度躍層位置，如圖3所示。定義造波區(qū)與傳播計算區(qū)密度躍層高度差ho稱為階躍高度，造波區(qū)長度Lo為階躍長度。具體計算參數(shù)及工況設置見表1，其中ηo為內(nèi)孤立波波幅，ηo/H為相對波幅。

圖2 數(shù)值模型布局

圖3 重力塌陷法制造下凹型內(nèi)孤立波示意圖

表1 數(shù)值模擬工況及相關參數(shù)

2.2 數(shù)值方法和邊界條件

采用大渦模擬(LES)技術實現(xiàn)三維數(shù)值波浪水槽中下凹型內(nèi)孤立波的形成和傳播；為保證數(shù)值計算域的質(zhì)量守恒，研究采用SIMPLE算法耦合N-S方程中的流速—壓力項，以此得到壓力場；利用二階中心差分格式對空間進行離散，時間的離散采用隱式格式。造波區(qū)左邊壁、槽底、柱體表面采用無滑移固壁邊界，水槽側壁采用滑移固壁邊界。水槽右邊壁采用Sommerfeld輻射型邊界[11]，以此可防止內(nèi)孤立波再反射，保證數(shù)值計算的精度；由于本研究中表面波波幅很小，影響甚微，水面采用“剛蓋”假定忽略表面波的作用[12]。

2.3 模型驗證

網(wǎng)格無關性分析及物理模型試驗對數(shù)模的驗證參照筆者之前發(fā)表的文章[13-14]，本文所采用的模型設計參數(shù)與兩篇文獻[13-14]一致。

3 數(shù)值結果及分析

3.1 各工況下CFn歷時曲線對比圖分析

圖4給出了不同波幅ηo/H工況下，柱體所受內(nèi)孤立波無量綱水平合力CFn隨時間的變化曲線。當內(nèi)孤立波傳經(jīng)柱體時，柱體所受的阻力與波幅成正比，均經(jīng)歷了先減后增、再減再增、最后趨于0的過程。由此判斷，內(nèi)孤立波會對柱體產(chǎn)生剪切及振蕩作用。此外，各工況下柱體的受力趨勢相似，且ηo/H越大，柱體達到其受力峰值時間越早。說明波幅越大，波的傳播速度越快，柱體所受阻力曲線的峰值越大，其衰減也越快。另外，隨著內(nèi)孤立波的傳播，其波動程度逐漸減小。這是由于以密度躍層為界的兩類不同密度水體的相互作用效果遠強于氣水交界面處氣體與水體的相互作用，更強的粘性耗散導致內(nèi)孤立波波幅更快的衰減。

柱體所受水平力變化可分解為3個過程：(1)內(nèi)孤立波抵達柱前：隨著內(nèi)孤立波逐漸接近柱體，柱體所受水平力逐漸增大。當內(nèi)孤立波與柱體接觸時，水平力達到其第一個峰值，即圖中的CFn-max；(2)在孤立波穿越柱體時：上、下兩層流體流速反向，柱體隨之承受一個方向為右的反作用力，逐漸抵消柱體左向水平力。隨后水平力變?yōu)?，繼而達到水平力谷值。(3)內(nèi)孤立波穿越柱體后：左向作用力持續(xù)增大，逐漸抵消右向作用力。水平力逐漸變大，在內(nèi)孤立波離開柱體時水平力達到其第2個峰值。最后由于尾波的存在，引發(fā)水平力的震蕩作用。不同波幅工況下，水平力均出現(xiàn)峰值及谷值。

3.2 各工況下壓強云圖對比分析

圖5給出了槽寬中心(Z=30 cm)截面的壓強云圖。各工況下(E1～E6)壓強分布圖均相同。數(shù)值水槽水深為80 cm，水面相對壓強為0 Pa，槽底相對壓強為7 840 Pa，符合實際情況。

圖4 各工況下CFn歷時曲線對比圖

圖5 橫向(Z=30 cm)截面壓強云圖

3.3 各工況下密度云圖對比分析

圖4中CFn-max表示內(nèi)孤立波無量綱水平合力峰值，本文數(shù)模數(shù)據(jù)提取及云圖繪制的特征時刻皆為柱體受力達到其峰值CFn-max所對應的最不利時刻，即CFn=CFn-max時刻。表3給出了各工況之間躍階高度ho之比hoi/hoi+1與波幅ηo之比ηoi/ηoi+1的對應關系。圖6(a)～(f)給出了不同波幅工況下柱體CFn=CFn-max時刻槽寬中心(Z=30 cm)截面的密度云圖。由圖6可知，隨著躍階高度ho的增大，波幅ηo不斷增加。

表3 ho之比與ηo之比對應關系

由此可知hoi/hoi+1越大，ηoi/ηoi+1越大。越大的躍階高度ho可激發(fā)更大波幅的內(nèi)孤立波，且波幅增幅隨躍階高度的增大呈拋物線分布。

圖6 不同波幅工況下橫向(Z=30cm)截面密度云圖(CFn=CFn-max時刻)

3.4 各工況下速度云圖對比分析

圖7給出了不同波幅工況下柱體CFn=CFn-max時刻Z=30 cm截面的速度云圖。如圖7所示，以密度躍層為界，上下水層流速反向。上層速度為正，下層速度為負，上層水體的流向與內(nèi)孤立波的波向一致，且上層水體速度絕對值要大于下層。此外，對比圖7(a)～(f)可以看出，隨ho不斷增大，上下水層流體質(zhì)點速度不斷增大，其中當ho=6.5 cm(工況E1)時，上層最大流體質(zhì)點速度為0.024 9 m/s，下層最小流體質(zhì)點速度為-0.011 85 m/s；ho=17.15 cm(工況E6)時，上層最大流體質(zhì)點速度為0.058 6 m/s，下層最小流體質(zhì)點速度為-0.037 6 m/s。表4給出了各工況波幅之比ηoi/ηoi+1與上、下水層質(zhì)點最小流速Vmin之比、最大流速Vmax之比的對應關系。

圖7 不同波幅工況下橫向(Z=0.3 m)截面速度云圖(CFn=CFn-max時刻)

表4 ηo之比與V對應關系

由此可知，流體速度與波幅成正比；hoi/hoi+1越大，流體Vmin和Vmax越大，且兩幅值的增幅也逐漸變大。說明波幅越大的內(nèi)孤立波能激發(fā)出更劇烈的流場。

3.5 各工況下流線圖及速度矢量圖對比分析

圖8(a)～(f)及圖9(a)～(f)分別給出了不同波幅工況下柱體CFn=CFn-max時刻Z=30cm截面的流線圖及速度矢量圖對比。如圖8所示，密度躍層上、下層水體流動反向，流動轉(zhuǎn)向發(fā)生在躍層附近；而在繞流作用下，無論是上層水體或下層水體，柱前、柱后均會產(chǎn)生漩渦，從而導致柱前、柱后流向一致；截面速度矢量總體上呈現(xiàn)出一個順時針方向的旋渦，旋渦的存在加劇水體紊亂，產(chǎn)生能量損失。柱前下層流體質(zhì)點速度較小，柱后層流體質(zhì)點速度較大。此外，隨著波幅η的不斷增大，流體各部分質(zhì)點速度矢量均不斷增大。

同時，距離密度躍層越近，內(nèi)孤立波波動影響越大。在流場中的流體質(zhì)點的運動速度越小，即在垂向位置上，隨著與密度分層面距離的增大，水平速度值不斷減小。綜合圖8～9的流線圖及速度矢量圖可知，密度躍層的上、下水層流速反向，這也是在內(nèi)孤立波環(huán)境中可誘導出強烈剪切流的根本原因所在，剪切作用致使水下結構物剪切疲勞甚至發(fā)生破壞。

圖8 不同波幅工況下橫向(Z=0.3 m)截面流線圖(CFn=CFn-max時刻)

圖9 不同波幅工況下橫向(Z=0.3 m)截面速度矢量圖(CFn=CFn-max時刻)

4 結束語

借助三維數(shù)值波浪水槽，通過繪制密度云圖、速度云圖、流線圖及速度矢量圖深入研究了不同內(nèi)孤立波波幅(ηo)對圓柱作用力的影響。相比過往研究，本文給出了一定的量化結果，細致對比了不同波幅工況下的流場特性，得到如下結論：

(1)兩工況躍階高度ho之比hoi/hoi+1越大，波幅ηo之比ηoi/ηoi+1越大。說明躍階高度是波幅大小的主要控因；

(2)對比不同波幅下柱體的受力歷時曲線，ho越大，ηo越大，波的傳播速度越快，柱體所受阻力的峰值也越大，表明波幅與柱體所受水平作用力成正比。同時大波幅工況下，柱體受力的衰減幅度也越大。

(3)從提取的流線圖及速度矢量圖可知，以密度躍層為界，上層水體質(zhì)點速度為正，下層水體質(zhì)點速度為負，導致內(nèi)孤立波具有巨大的剪切力，易對水下結構物造成破壞。隨著波幅ηo的不斷增大，流動區(qū)域的水動力越劇烈，上、下水層流體質(zhì)點速度不斷增大。