波浪作用下直立結構物附近強湍動摻氣流體運動的數值模擬1)

鄧 斌 王孟飛 黃宗偉 伍志元 蔣昌波,2)

*(長沙理工大學水利工程學院,長沙410114)

?(水沙科學與水災害防治湖南省重點實驗室,長沙410114)

**(寧波中交水運設計研究有限公司,浙江寧波315800)

引言

波浪從外海傳播至海岸線附近的過程中,在地形和海岸防護建筑物的影響下,波浪與建筑物相互作用常弓起水流強烈的湍動和翻滾,伴隨著大量空氣卷入,產生大量氣泡,該過程涉及到復雜的氣液兩相流問題.建筑物附近水體的摻氣、以及氣泡的演化和輸運過程影響著近岸結構物的波浪沖擊效應[1]、建筑物局部沖刷、波能的耗散等[2-5].與此同時,近岸帶波浪受淺水效應容易發生破碎,破碎過程中會使大量氣體卷入水體.因此,研究近岸波浪與海岸防護建筑物相互作用過程中摻氣水流的運動特性,對準確計算建筑物受力[6-7]、越浪量[8-9]及其穩定性防護[10-11]等具有重要的理論價值和科學意義.

國內外學者一直致力于波浪作用下破碎區水動力及摻氣特性研究,通過實驗和數值手段獲取了部分破碎區氣泡輸運的分布規律,并建立相關數學模型.其中,前人研究主要針對破碎區內空隙率、氣泡直徑以及數量等參數進行量化分析,部分學者基于固定坡度概化近岸海灘建立了描述破碎區空隙率的經驗公式,一定程度上解釋了沖瀉區內空隙率、水動力特征及其與波浪要素之間的關系.如:鄧斌等[12]通過光纖探針測量了卷破波作用下破碎區內空隙率的變化規律,給出了空隙率與波陡、相對水深以及水平位置的經驗關系.Mori 等[4]通過現場測量與實驗研究探討了沖瀉區內氣體夾帶的影響,發現氣泡大小的空間分布規律與氣泡大小密切相關,并提出氣泡大小的Hinze 尺度為1 mm.Deane 等[13]對比分析了物理實驗和現場觀測兩種研究方法下氣泡分布的差異問題,并提出氣泡概率密度與氣泡半徑之間呈現10/3 次方的關系.類似工作還有文獻[14-21]等.

然而,上述研究主要集中于近岸波浪破碎區至沖瀉區摻氣特性與水動力特征的研究,未考慮建筑物的影響效應.與此同時,現有針對波浪與建筑物相互作用的研究大多未考慮水體摻氣的影響[22-23].如:Xiao 等[24]基于不可壓雷諾平均的N-S 方程(RANS)和k-ε 湍流模型建立了數值波浪水槽,研究了波浪作用于不同高度建筑物上波浪力和傾覆力的變化規律.Nakamura 等[25]基于實驗和數值模擬方法研究了建筑物前沖刷坑深度與水流流速以及有效應力之間的關系.類似的相關研究還有文獻[26-28].近年來,少量研究開始關注波浪與建筑物作用下水體摻氣對建筑物受力的影響,如Hull 等[5]通過實驗研究發現卷破波作用下氣體卷入會導致建筑物受力增大的現象.Hsiao 等[29]通過實驗和數值模擬研究了孤立波與梯形海堤相互作用的過程,定性分析了海堤后渦動結構及氣泡夾帶特性的演化特征.部分學者[30-32]研究發現跨海大橋橋板下方的主梁之間包裹的空氣是增加波浪對橋梁上升力的一個重要因素.然而,上述研究僅針對水體摻氣對建筑物受力影響的討論,一定程度上推進了波浪與建筑物作用下摻氣水流運動特性的認識.

綜上所述,盡管目前對近岸帶波浪破碎摻氣特性以及波浪與建筑物相互作用等問題有了一定認識,建立了破碎區內摻氣特性與水動力之間的關系,但準確理解摻氣后形成的氣泡與水動力特性之間的關系仍然是一項非常具有挑戰的工作,有必要深入了解波浪與建筑物相互作用下摻氣特性的內在運動機制.為此,在前人研究基礎上,本研究通過建立三維精細數值水槽和對應物理實驗模型,對波浪作用下建筑物附近的水動力以及相應的摻氣特性進行研究,為進一步深入揭示波浪與建筑物相互作用、水體摻氣與氣泡輸運機制奠定基礎.

1 數學模型

1.1 控制方程與數值方法

基于三維連續不可壓縮Navier-Stokes 方程建立數值模型.模型通過求解混合RANS (Reynolds average Navier-Stokes)和LES (large eddy simulation)的延遲分離渦湍流模型S-A DDES(Spalart-Allmaras delayed detached eddy simulation),其中邊界區和湍流區分別采用RANS 和LES 模擬計算,該方法通過RANS模型可在近壁面處實現以較少的網格條件下得到較準確的解,又可通過LES 模型提高內部流場湍流的解析度.S-A DDES 模型是在S-A DES 模型基礎上通過重新定義長度尺度,從而解決模型應力退化的問題,實現兩種不同的湍流模型的有效結合,并解決區域之間的過渡問題,其中湍流區采用亞格子模型(SGS)對小渦進行計算.對Navier-Stokes 方程進行濾波處理后以實現大渦和小渦分離的控制方程為

式中,xi和xj為笛卡爾坐標下的分量,ui和uj分別為坐標系下對應的流速場,p為壓強,ρ 為密度,在大渦和小渦的求解中,為亞格子應力項σ 為表面張力系數,κ 為界面曲率,α為VOF 流體體積函數.由于LES 中亞格子應力項未知,因此采用亞格子模型進行封閉,形式為,Sij為應變率張量,νt為湍流黏度系數為分子黏度),采用S-A 模型對進行封閉求解

式中,右側第一項為渦黏性產生項,第二項為輸運項,第三項為耗散項;Cb1,Cb2,Cw1和δ 均為常量;fw為與和d有關的函數;d為長度尺度,即網格點到壁面的最小距離.基于IDDES 方法修改S-A 湍流模型長度尺度d

其中,CDES為常數,取值0.65;fd為延遲函數,實現從RANS 模式到LES 模式的延遲;L△為網格尺度,L△=max(△x,△y,△z),△x,△y和△z為網格點3 個方向上的網格尺寸.

采用VOF 方法捕捉自由表面,弓入流體體積分數描述水和空氣兩種流體的分布,其中α=1 為水相,α=0 為氣相,α 介于0 和1 之間為氣體和水體交界面.同時,在傳統體積輸運方程的基礎上本模型增加了人工壓縮項[33],可以有效保證方程解在自由界面附近自動滿足有界性,避免產生界面模糊,且對外部流場無影響,修正后的體積分數輸運方程

1.2 邊界設置

數值模型中,水槽底邊界以及直立結構物壁面采用固壁邊界條件,水槽前后兩側設置為對稱性邊界條件,水槽右側設置為自由出流邊界條件,水槽左側采用Jacobsen 等[34]提出的修正速度入口方法進行造波,該方法中造波邊界的類型分為干網格、濕網格和臨界網格3 種,其中干網格的邊界條件

式中,n為邊界網格面上的單位法向量,濕網格邊界條件(u和φB)是根據相應的波浪理論計算得到.模型對臨界網格進行了濕面修正,波浪理論計算的自由液面曲線與網格邊界存在2 個交點,取兩點間直線作為臨界網格的濕面邊界,從而得到網格濕面部分的中心位置,以此計算出臨界網格的入流速度和體積分數,這可避免了入口處自由液面的虛擬振蕩問題.

1.3 模型與工況設置

數值計算模型設置如圖1 所示.其中,岸灘概化地形分為兩部分:一部分是距離造波邊界15 m 處布置一光滑1:10 的斜坡;另一部分是在斜坡頂端的平臺起始處布置直立結構物模型,模型高0.1 m、寬0.12 m,平臺長3 m,右端接1:10 斜坡,用于消減波流能量.計算以斜坡起點作為坐標原點,以波浪傳播方向為x軸正方向,以水平垂直數值水槽方向為y方向,以豎直向上方向為z軸正方向,由此建立三維坐標系,同時,將xoz坐標系作為二維坐標系.計算選擇規則波作為入射波,具體工況設置詳見表1 所示.

圖1 試驗布置及測量分布圖Fig.1 Sketch of the experimental setup and instruments

表1 研究工況Table 1 Study cases

為驗證數值模型的計算結果,在長沙理工大學水沙科學與水災害防治湖南省重點實驗室的波浪水槽中開展相應的物理模型實驗.物理模型實驗布置同數值水槽一致,水槽的左端配有活塞式造波機,斜坡地形和直立結構物模型均采用有機玻璃制作.模型實驗沿水槽方向布置6 個波高測點(具體位置見表2),采集儀器為德國General Acoustic 公司生產的ULS 80D 型號超聲波非接觸式浪高儀,采集頻率為50 Hz,并在實驗前通過水槽內增減水進行標定,該儀器測量誤差小于3.0×10-4m.與此同時,實驗采用德國Photron 公司生產的高速相機(FASTCAM Mini UX100,采樣頻率:4000 Hz)記錄不同波況下波浪破碎的演化過程,并采用法國A2 Photonic Sensors 公司生產的氣泡測量系統(型號:B-POP,傳感器為L 形光纖探針)測量相應位置處氣泡的數量與空隙率,測量過程中通過記錄氣泡穿過光纖探針過程中電信號的歷時曲線,從而獲取氣泡的數量以及空隙率,其詳細測量原理及方法見文獻[12],該光纖探針能夠測量直徑大于100 μm 以及運動速度為0.06～6.25 m/s 范圍內的氣泡,最大誤差小于6%.氣泡測量位置分別以5 cm 和1 cm 的間隔在水平x方向和豎直z方向上進行測量,具體測量位置詳見表1.實驗中所有實驗儀器設備均采用同步系統同步采集實驗數據,且物理實驗在淡水和周圍空氣下進行.

表2 浪高儀位置Table 2 The location of wave gauges

1.4 網格劃分與參數設置

由于氣泡測量系統測量氣泡直徑的范圍為2～70 mm,結合預備實驗的結果,數值計算區域的網格沿x方向將結構物前后各1.0 m 區段設為核心區域,采用1 mm 的均勻結構網格,核心區域向外以結構網格的方式漸變過渡到5 mm;數值水槽在y方向上長0.08 m,沿y方向采用1 mm 的均勻網格,z方向同x方向類似在結構物附近網格大小為1 mm,核心區域向外以結構網格的方式漸變過渡,模型的計算網格總數達到2.4 億,結構物附近網格劃分的局部細節如圖2 所示.

圖2 結構物附近的網格劃分Fig.2 Computation grids near the structure

為更加準確地模擬波浪破碎作用下氣泡的運動,數值模型采用連續表面張力模型(CSF)考慮表面張力的影響,表面張力系數σ=0.072 N/m,模型相關參數見表3.數值計算采用自適應時間步長,以保證克朗數Cr小于1.對于單元網格下的庫朗數的定義如下

式中,δt為時間步長,|U|為水質點通過單元網格的最大速度,δx為單元網格大小.

表3 數值模型參數設置(ρw和ρa分別為水和氣的密度,μw和μa分別為水和氣的黏滯系數)Table 3 Parameter setting(ρwand ρaare density for water and air respectively,μwandμaare viscosity for water and air respectively)

2 模型驗證

圖3 自由液面歷時曲線對比Fig.3 Comparison between numerical and experimental time series of free surface elevation

圖3 為不同測點處波面歷時曲線的數值模擬結果與實驗結果的對比圖(h=0.34 m,H=0.08 m,T=1.6 s).采用S kill數驗證計算值和實驗值的吻合度,其中S kill數計算采用Willmott(1981)的方法

式中,Xnum和Xexp分別表示數值計算值和實驗值,上畫線表示取平均值.該工況下計算得到的S kill數均大于0.85,表明數值結果與實驗值吻合較好,該數值模型能夠真實地模擬規則波在斜坡上破碎的歷時變化特征.

從圖3 中還可以看出隨著水深在斜坡上逐漸變淺,波浪發生淺水變形,波峰逐漸變尖,波形表現出明顯的不對稱性(3#),隨后波浪在結構物附近發生破碎,破碎后的波高迅速減小,破碎波以段波形式繼續傳播.其中,5# 和6# 位置處數值模擬結果略高于實驗測量結果,這是由于數值模擬結果中波高數據考慮了波面處的液滴飛濺.此外,對比實驗照片圖4(a)和三維數值結果圖4(b),發現數值結果較好地再現了實驗中波浪破碎后產生的強湍動摻氣水流,以及先沖擊結構物近海側,隨后沿結構物上爬,形成強射流現象,最終產生越浪,進一步表明建立的數值模型能合理模擬結構物附近波浪破碎后形成強湍動摻氣水流的運動過程.

圖4 數值與實驗對比圖Fig.4 Qualitative comparisons of free surface evolution between laboratory images and numerical data

3 結果分析與討論

3.1 摻氣水流與結構物的作用過程

圖4 給出case2(h34H8)波浪與結構物相互作用下數值和實驗過程的演化結果.結果表明波浪在斜坡上破碎后,空氣被夾帶進入水中形成一個大氣腔和眾多小氣泡,當波浪持續破碎傳播至結構物近海側時,產生強摻混的摻氣水流并沖擊結構物,之后波浪沿結構物上爬,產生明顯的射流和滴液的飛濺現象,最終造成結構物后方的越浪.此外,為分析摻氣水流中氣泡在輸運過程中的破碎、溢破等現象,圖4(c)給出了結構物附近氣泡輸運過程細節圖(圖4(b)中紅色方框的局部圖).從圖4(c)中可知,波浪在堤前破碎時,眾多氣泡形成于近海側.氣泡在隨后輸運過程中氣泡周圍一直伴隨著渦量的存在,并且大都是正負渦量相互作用.大氣泡在經歷多次破碎后形成小氣泡,氣泡在分裂為多個氣泡時,中間的連接處往往伴隨著負渦量,而周圍為正渦量,最終在正負渦量的共同作用下氣泡逐漸由單個分裂為多個.通過實驗以及數值結果觀察發現破碎后較大的氣泡易受水動力因素影響,表面積和體積不斷發生變化,浮力作用能驅使大氣泡迅速上升至自由表面,最終在自由表面處破碎;而經歷多次破碎后形成的較小氣泡的表面積和體積很少改變,且受浮力影響較小,在水體中存留時間較長,直至水流平穩后,才得以上浮消散.

圖5 為結構物附近氣泡的演化細節圖,其中左圖可以看出在氣泡的變形處分布著大量的正負渦量,表明氣泡的形狀甚至分裂受正負渦量的相互作用影響.從中、右圖可見氣泡群的分布范圍與湍動能的分布范圍大體一致,氣泡群越密集,即氣泡數量越多,湍動能越大,范圍越廣,表明水體中湍動能的分布與氣泡群的數量密切相關;此外,還可以看出氣泡內部流場與周圍流體的流場略有不同,但整體流速與周圍水體的流場相匹配,由此可推斷,氣泡的整體運動軌跡與氣泡周圍的水體流場密切相關.

圖5 氣泡演化細節(左:自由液面與渦量的分布;中、右分別為α=0.1,0.3,0.4 時自由液面,以及相應切片上的湍動能)Fig.5 The detail characteristic of local 3-D bubble(The distribution of free surface and vorticity is given in the left figure;The free surface at α=0.1,0.3,0.4,and the turbulent kinetic energy on the corresponding slice in the right figure)

3.2 流速場

圖6 結構物附近流場Fig.6 Time series of velocity field near the structure

現有研究表明,波浪與結構物相互作用產生摻氣水流,將直接影響流場分布特征[35].圖6 給出不同時刻下case3 (h34H10)x–z平面內結構物附近的流場矢量圖.t=8.2 s 時刻,波浪破碎后的水體撞擊到結構物,形成向上飛濺的射流,在x=2.6 m 處形成一個較大的氣腔,并在x=2.8 m 處通過水舌卷入大量的氣體,形成順時針的漩渦,此時水體流速較大,其值在2 m/s 左右,結構物附近摻雜了大量的小氣泡,流場十分復雜;隨后波浪沿結構物上爬形成射流(t=8.3 s),但向上的速度逐漸減小,波浪繼續向前運動,翻卷的水舌繼續夾帶空氣形成摻氣水流;在x=2.8 m 處順時針的漩渦將大氣腔分裂為兩個氣泡,同樣x=2.98 m 處的氣腔在水流剪切作用下也分裂為多個氣泡;隨后射流逐漸回落(t=8.4 s),速度達到2 m/s 左右,在x=2.7 m 處的大氣腔在水流作用下繼續向前運動,而上一時刻在x=2.8 m 處的兩個氣泡融合為一個氣泡,波浪逐漸在結構物前產生雍水,流速開始減小,水面處的水平流速減小至0 左右;射流完全落下(t=8.5 s),水體表面開始出現回流,水體底部仍然存在向前的水平流速.因此,在x=2.88 m處形成逆時針方向的旋渦,大氣泡逐漸上浮、溢出.t=8.6 s 時,水體底部流速逐漸減小,水面處回流流速逐漸增加,水體容易產生較大的水流剪切力,結構物附近的旋渦范圍逐漸增大,結構物前形成較多小氣泡;水體底部開始回流(t=8.7 s),水體表面的反向水平流速較大,底部流速較小,小氣泡逐漸上浮耗散,但仍存在大量漩渦;之后(t=8.8 s)水體內幾乎沒有較大的氣泡,小氣泡逐漸上浮,水體回流速逐漸增大.

3.3 渦量場

為進一步研究渦對氣泡演化及水體摻氣的影響.圖7 給出不同時刻case3 (h34H10)工況下x–z平面內結構物附近的渦量分布圖.通過y方向渦量反映渦結構的強度與方向.其中,y方向上的渦量表示為

圖7 結構物附近渦量場Fig.7 Time series of vorticity field near the structure

式中,u和w分別為x方向和z方向的速度.

圖7 中t=8.2 s 時刻,在x=2.6 m 處的氣腔附近伴隨著較大的渦量,x=2.8 m,2.98 m 處卷氣區域周圍正負渦量相間,順時針與逆時針的漩渦相互作用容易形成對氣腔的剪切.水體中的旋渦隨著氣腔向前運動(t=8.3 s～8.4 s),渦量逐漸減小.當射流完全落下(t=8.5 s),水體表面開始回流,大氣泡逐漸上浮,水面附近展現出較多負渦量,水體底部呈現較多正渦量,印證了上文流速場的分析.當t=8.6 s 時,水體中負渦量逐漸增強,范圍逐漸增大,結構物前氣泡間存在明顯負渦量;之后(t=8.7 s)由于水體的回流,結構物前的渦量與小氣泡逐漸上浮消散.

3.4 湍動能

圖8 為不同時刻case3 (h34H10)工況下x–z平面結構物附近的湍動能分布圖.其中,湍動能表示為

式中,u′和w′分別為橫向和縱向的脈動速度.

圖8 結構物附近湍動能Fig.8 Time series of turbulent kinetic energy field near the structure

當t=8.2 s 時,在x=2.8 m 處的波浪破碎形成的卷舌與水體接觸部分存在較大湍動能,射流與結構物前端接觸區域的湍動能較大,其值為0.02 m2/s2.當卷舌運動到x=2.9 m 處(t=8.3 s),湍動能繼續增大,呈現水面處以及結構物附近湍動能較大,水體內部以及破碎后的水體湍動能較小的特征.此外,x=2.78 m 處的卷氣區域周圍湍流動能較大,結構物附近的高速摻氣水流湍動能較大(t=8.4 s).同樣,t=8.5 s 時刻,在自由液面處以及結構物附近的摻氣區域湍動能較大,表明水體的摻氣能增強水體的湍動.隨后(t=8.6 s)水體中流速逐漸減小,但靠近結構物的高度摻氣區域的湍動能仍高于其他區域,這與Mori 等[4]的研究類似.

3.5 空隙率、氣泡數量與湍動能的關系

上述討論表明,湍動能高的區域與氣泡數量較多的區域具有較高空間重合性,表明湍流強度與氣泡數量之間存在一定相關性.實際上,已有研究也證實破碎波中湍動能與氣水混合物的空隙率呈現正比關系[4].為探討結構物附近摻氣流體中空隙率與湍動能是否存在類似關系,圖9 給出任一周期內通過L型光纖探針測量到各測點的平均空隙率與平均湍動能的變化趨勢.可見,3 個工況下湍動能與空隙率整體也呈現線性關系,值得注意的是,隨著波高的增大,空隙率也逐漸增大,但空隙率與湍動能的線性關系有所減弱(如case3),這是由于圖中相對分散的點均是離直立結構物較近位置(x=2.96,3.00 m)和結構物上(x=3.04 m).分析可知,與斜坡上(無結構物)破碎波內摻氣水流特性不同的是,結構物的存在使得湍動能與空隙率的線性關系變弱,這與波高增大時,破碎波與直立結構物作用產生的流體飛濺和越浪導致空隙率在較短時間內發生劇烈變化,以及波浪破碎及其與結構物非線性相互作用更加明顯等相關.圖10 給出了實驗測量得到的氣泡數量與空隙率之間的關系,結果同樣顯示兩者整體存在正相關關系,從而可推斷湍流動能與氣泡數之間存在一定關系.因此,圖11 給出了任一周期內通過L 型光纖探針測量到各測點的平均氣泡數與平均湍動能之間的變化趨勢,結果同樣顯示湍動能與氣泡數整體存在正相關關系,且隨著波高的增大氣泡數逐漸增多,印證了上文關于氣泡數與湍動能的分析.為更好地描述各回歸方程的可靠性,表4 給出了圖9～圖11 中各回歸方程參數的相關統計值.綜上可見,波浪作用下直立結構物附近湍動能的分布與與摻氣水流部分特征參數(氣泡數量、空隙率)整體呈現正相關關系.

圖9 時間平均空隙率和時間平均湍動能關系Fig.9 Relationship between time average of void fraction and turbulence kinetic energy

圖10 時間平均空隙率和時間平均氣泡數關系Fig.10 Relationship between time average of void fraction and bubble number

圖11 時間平均空隙率和時間平均湍動能關系Fig.11 Relationship between time average of bubble number and turbulence kinetic energy

表4 線性回歸方程相關參數統計表Table 4 Parameter statistics of linear regression equation

4 結果與展望

基于Navier-Stokes 方程的三維數值模擬技術,結合物理模型實驗,模擬了波浪作用下結構物附近強湍動摻氣水流的水動力學特性,研究表明:

(1)建立的數值模型能精確地捕捉波浪作用下結構物附近的自由液面的變化以及氣泡輸運過程,較好地描述氣體卷入所形成的氣腔形態以及多氣腔之間的融合、分裂等過程.

(2)分析結構物附近波浪破碎后的相關流場和湍動能、渦量場,發現破碎波與結構物相互作用過程中水體內湍動能以及渦量的變化與水體內氣泡變形、破碎與融合等現象相關.其中,氣泡云的范圍與水體湍動能的在空間分布相重疊,正負渦量的相互作用能弓起氣泡的破碎和融合,氣泡的融合、破碎又會弓起水體流場的變化.

(3)研究建立了結構物附近湍動能與摻氣特性相關參數之間的關系,發現波浪作用下結構物附近湍動能的分布與摻氣水流部分特征參數(如氣泡數量、空隙率)呈現線性正相關關系.

基于數值模擬和物理實驗揭示了波浪作用下直立結構物附近強湍動摻氣水流中氣泡特性與水動力之間的關系,為進一步認識波浪作用于結構物的內在機制提供了理論基礎,但文中僅考慮了定床條件下特定波浪參數以及概化的結構物,考慮更多岸灘特征、結構物型式、波浪參數等條件下強湍動摻氣水流運動規律及其與結構物相互作用機理還有待進一步研究.