邢飛，孫大鵬，夏志盛，吳浩

（大連理工大學(xué)海岸和近海工程國家重點實驗室，大連116023）

波浪與開孔沉箱相互作用的CLEAR?VOF數(shù)值模擬

邢飛，孫大鵬，夏志盛，吳浩

（大連理工大學(xué)海岸和近海工程國家重點實驗室，大連116023）

摘要：文章利用三步有限元方法離散Navier?Stokes方程，并借助CLEAR?VOF方法、大渦模擬方法追蹤水體自由表面和模擬湍流效應(yīng)，建立了二維數(shù)值波浪水槽。籍此研究了規(guī)則波作用下開孔沉箱的反射率和波浪總水平力，并與物模試驗結(jié)果和其他學(xué)者的數(shù)模結(jié)果以及工程實例進行了對比，吻合良好，為波浪與開孔沉箱相互作用問題的研究探索了一種新的數(shù)值模式。

關(guān)鍵詞：三步有限元；CLEAR?VOF方法；開孔沉箱；反射率；總水平力比

自1961年Jarlan提出開孔沉箱以來，國內(nèi)外大量采用了這種形式的防波堤來減小波浪的反射及其對建筑物的作用力。對于波浪與開孔沉箱相互作用的研究，目前主要有物模試驗和數(shù)值模擬兩種途徑，陳雪峰、李玉成等［1-2］利用物理模型試驗研究了單層和雙層開孔沉箱在規(guī)則波和不規(guī)則波作用下的反射率、水平力及點壓力分布；張玉彬、孫大鵬等［3］利用物理模型試驗研究了明基床上開孔沉箱在規(guī)則波作用下的受力情況；施曉迪等［4］利用物理模型試驗研究了開孔率和開孔范圍對開孔沉箱消浪性能的影響；朱大同［5］得出了波浪與開孔沉箱相互作用的反射率的嚴(yán)格封閉解析公式；陳雪峰［2］采用改進的VOF方法結(jié)合k?ε模型建立二維數(shù)值波浪水槽，完成了線性規(guī)則波和不規(guī)則波與單層開孔沉箱相互作用的數(shù)值模擬。目前，對于波浪與開孔沉箱相互作用的物模試驗成果較為豐富，但是數(shù)值模擬研究成果相對較少，其中，時域內(nèi)準(zhǔn)確地完成對自由水面的追蹤，是建立波浪與開孔沉箱相互作用數(shù)值模式過程中的難點問題之一。

Ashgriz等［6］于2004年提出了CLEAR?VOF（A combination of a Computational Lagrangian?Eulerian Advec?tion Remap and the Volume of the Fluid method）自由表面重構(gòu)技術(shù)，其精度較高且適用于任意多邊形網(wǎng)格。本文借助CLEAR?VOF方法追蹤水體自由表面，建立了波浪與開孔沉箱相互作用的三步有限元數(shù)值模式。并分析了規(guī)則波作用下開孔沉箱反射率和波浪總水平力的計算結(jié)果，為研究波浪與開孔沉箱相互作用探索了一種新的數(shù)值模式。

1　數(shù)值計算方法

1.1控制方程及其離散

基本控制方程為二維不可壓縮粘性流體的連續(xù)性方程和Naiver?Stokes方程

由于直接進行湍流模擬在現(xiàn)有的計算條件下難以實現(xiàn)，本文采用大渦模擬的方式模擬湍流效應(yīng)，通過濾波的手段將實際流場分解為大尺度和小尺度的渦流場，對其中大尺度的流場物理量進行直接數(shù)值模擬，小尺度的物理量通過模型化的手段加以反應(yīng)。經(jīng)濾波處理并引入廣義Boussinesq渦粘假設(shè)對控制方程進行封閉之后得到

式中：有效粘性系數(shù)νe為運動粘性系數(shù)ν與湍流粘性系數(shù)vt之和。

三步有限元方法比傳統(tǒng)的有限元方法精度更高，數(shù)值實現(xiàn)過程相對簡單且高效。對Taylor時間展開取三階精度，將每個Δt時間步長分為Δt/3，Δt/2和Δt三步來完成，其基本思想如下

本文計算中，三步有限元展開N?S方程的形式為

在Δt/3步，利用式（8）由上一時刻求出的速度un和壓力pn求出該時刻的速度un+1/3；在Δt/2步，由式（9）易求得該時刻的速度un+1/2；而在Δt步，由于式（10）中存在未知量pn+1，因此引入壓力泊松方程［7-8］

至此，方程組封閉。可采用Galerkin加權(quán)余量法進行離散求解。

1.2CLEAR?VOF追蹤流體自由表面

與傳統(tǒng)的基于有限差分方法提出的VOF方法不同，CLEAR?VOF方法不需要通過求解輸運方程來更新F函數(shù)，而是通過歐拉網(wǎng)格節(jié)點拉格朗日運動和計算幾何手段實現(xiàn)。在非結(jié)構(gòu)化歐拉網(wǎng)格計算區(qū)域中，當(dāng)F=0時，單元內(nèi)無水，不存在流體多邊形；當(dāng)F=1時，流體多邊形與網(wǎng)格重合，流體多邊形的頂點即網(wǎng)格單元節(jié)點；當(dāng)0＜F＜1時，流體多邊形的頂點包括網(wǎng)格單元節(jié)點和自由表面界線與網(wǎng)格邊線的交點。

以圖1為例，圖中實線為計算區(qū)域的不規(guī)則歐拉網(wǎng)格，對于某時刻t，各單元節(jié)點的速度已知，網(wǎng)格邊線上任意點的速度可由插值方法求出，單元a-b-c-d-e為自由表面單元（其中e-d為自由表面界線）；根據(jù)拉格朗日思想計算流體多邊形的運動，則t+Δt時刻，該單元運動至a1-b1-c1-d1-e1，初步實現(xiàn)了流體向周邊網(wǎng)格的輸運；通過復(fù)雜的計算幾何方法得出t+Δt時刻下流體多邊形留在母單元內(nèi)的面積和進入母單元相鄰單元的面積，對所有網(wǎng)格單元進行計算，得到新時刻的F函數(shù)、自由表面界線和各流體多邊形的頂點速度，并將此作為下一時刻的已知條件循環(huán)計算，直到計算結(jié)束。

圖1CLEAR?VOF方法中的流體多邊形示例Fig.1Sketch of fluid polygon of CLEAR?VOF

2　數(shù)值波浪水槽的建立

2.1線性波浪模擬

利用上述數(shù)值計算方法，如圖2所示建立二維數(shù)值波浪水槽，水槽前端為造波區(qū)，為得到理想的線性規(guī)則波，采用Dong C M等［9］提出的無反射造波理論，理論波面η和造波板運動速度U分別為

圖2　波浪數(shù)值水槽模擬計算域Fig.2Computational domain for numerical flume

式中：H為理論波高；k為波數(shù)；ω為頻率；d為水槽水深；ηRe為造波板前實際波面高程；ηTh為造波板前理論波面高程。

水槽末端安放由Larsen J等［10］提出的海綿吸收層消波邊界，在海綿吸收層內(nèi)按下式對速度進行衰減

式中：us為海綿吸收層內(nèi)單元節(jié)點的速度；μ（x）為衰減函數(shù)；α為常系數(shù)，本文取1.1；Ls為網(wǎng)格單元中心到海綿吸收層末端邊界的距離；λ為海綿吸收層的厚度，本文取最大波長的1.5倍，為3.6 m。

數(shù)值水槽工作區(qū)長7.7 m，高0.6 m，水深d為0.4 m，波高H取0.10 m，0.12 m，周期T取0.86 s，1.0 s，1.2 s，1.4 s，計有8種波要素組合。波浪傳播穩(wěn)定后，各工況的線性規(guī)則波浪與理論解吻合良好。限于篇幅，本文以波高0.10 m、周期1.0 s工況為例給出了水槽工作區(qū)末端位置的波高歷時曲線，如圖3所示，圖中實線表示理論值，圓圈表示本文計算值。

圖4　二維數(shù)值波浪水槽中開孔沉箱的布置圖Fig.4Sketch of perforated caisson in two?dimensional numerical flume

2.2開孔沉箱結(jié)構(gòu)物的設(shè)置

水槽工作區(qū)末端安放開孔沉箱如圖4-a所示，開孔沉箱的開孔板尺寸如圖4-b所示。寬68 cm，高40 cm的開孔面上布置大小相同的高3.8 cm，寬14.2 cm的12個孔，開孔率為23.8%，消浪室寬度為0.3 m，消浪室內(nèi)的水深為0.2 m。

數(shù)值模擬過程中采集距開孔板2.5 m、2.9 m和3.5 m處的波面時間過程信息用以計算分析反射率。

在數(shù)值計算的過程中，將開孔沉箱的后實體直墻和開孔板的固壁部分定義為自由滑移邊界（全反射邊界），即邊界法向速度設(shè)置為u=0。

3　數(shù)值計算結(jié)果及分析

開孔沉箱的消浪效果常用反射率和水平力來加以評價。為驗證本文數(shù)值模式研究波浪與開孔沉箱相互作用的適用性和準(zhǔn)確性，將本文計算結(jié)果與物理模型試驗結(jié)果（物理模型試驗在大連理工大學(xué)海岸和近海工程國家重點實驗室的波流水槽進行，數(shù)值試驗和物模試驗的結(jié)構(gòu)物尺寸、試驗波要素等保持一致）［2］、陳雪峰［2］改進的VOF方法結(jié)合k?ε模型的數(shù)值計算結(jié)果進行對比驗證。

3.1反射率的比較

反射率按照合田良實［11］提出的兩點法來計算。陳雪峰等［2］通過分析物理模型試驗數(shù)據(jù)，得到了反射率Kr與入射波長L、水深d、消浪室寬度B和開孔率α之間的經(jīng)驗關(guān)系

圖5-a、圖5-b分別給出了波高為0.10 m和0.12 m工況下反射率的比較。其中橫坐標(biāo)為波長，縱坐標(biāo)為反射率；實線表示經(jīng)驗公式（16）的結(jié)果，圓點表示本文的數(shù)值計算結(jié)果，三角形點表示陳雪峰的數(shù)值計算結(jié)果。由圖5可知，本文數(shù)模的反射率結(jié)果比陳雪峰數(shù)模結(jié)果略小，但與式（16）吻合較好。

圖5　反射率的比較Fig.5Comparison of reflection coefficient

3.2水平力的比較

以F1/F0作為開孔沉箱總水平力的分析指標(biāo)，其中開孔沉箱所受水平總力F1為開孔板的迎浪面、背浪面及后實體墻的合力，實體沉箱所受水平總力F0為其迎浪面所受水平力。陳雪峰等［2］通過分析物理模型試驗數(shù)據(jù)，得到了波峰（谷）總水平力比F1/F0與入射波長L、波高H、水深d、消浪室寬度B和開孔率α之間的經(jīng)驗關(guān)系

圖6、圖7分別給出了波高為0.10 m和0.12 m工況下（a）波峰總水平力比的比較和（b）波谷總水平力比的比較。其中橫坐標(biāo)為波長，縱坐標(biāo)為力比；實線表示經(jīng)驗公式（17）或經(jīng)驗公式（18）的結(jié)果，圓點表示本文的數(shù)值計算結(jié)果，三角形點表示陳雪峰［2］的數(shù)值計算結(jié)果。由圖可知，本文計算的波峰（谷）總水平力與陳雪峰數(shù)模結(jié)果以及式（17）和式（18）吻合較好且規(guī)律性明顯。

圖6　波峰（谷）總水平力比的比較（h=0.10 m）Fig.6Comparison of ratio of the total horizontal forces（h=0.10 m）

圖7　波峰（谷）總水平力比的比較（h=0.12 m）Fig.7Comparison of ratio of the total horizontal forces（h=0.12 m）

4　工程實例驗證

大連港大窯灣港區(qū)#11～#16泊位的規(guī)劃設(shè)計中，采用圖8所示的開孔沉箱。靜水水深14.44 m，沉箱迎浪面規(guī)則分布16個方形開孔，沉箱內(nèi)部以碎石緊密填充。本文以該工程為例，選取未發(fā)生越浪的工況進行比較，由于原始資料數(shù)據(jù)有限，本文僅給出實際尺度開孔沉箱結(jié)構(gòu)波浪反射系數(shù)的數(shù)值驗證。

表1給出了不同工況規(guī)則波作用下實際尺度開孔沉箱反射率的本文計算值和設(shè)計值的比較。由表1可知，本文計算結(jié)果與實際工程反射率吻合良好，說明本文所建立的數(shù)值模式可適用于開孔沉箱結(jié)構(gòu)的工程設(shè)計應(yīng)用。

圖8　開孔沉箱防波堤工程實例Fig.8Engineering project of perforated caisson

表1　實際尺度開孔沉箱反射率的比較Tab.1Comparison of reflection coefficient of engineering project

5　結(jié)語

本文基于Navier?Stokes方程，借助大渦模擬方法模擬湍流效應(yīng)，并成功地實現(xiàn)CLEAR?VOF追蹤水體自由表面，建立了波浪與開孔沉箱相互作用的三步有限元數(shù)值模式。進而在數(shù)值波浪水槽中對開孔沉箱前的波浪反射率和沉箱所受波浪水平力進行了數(shù)值模擬，本文數(shù)值計算結(jié)果與物理模型試驗結(jié)果［2］吻合較好且規(guī)律性明顯，并且進行了實際尺度開孔沉箱防波堤的數(shù)值驗證，表明本文開發(fā)的波浪數(shù)值模式對于開孔沉箱的波浪反射率和波浪力數(shù)值模擬具有良好的適用性和準(zhǔn)確性，為該問題的研究探索了一種新的數(shù)值模式。

參考文獻：

［1］陳雪峰，李玉成，孫大鵬，等.波浪與開孔沉箱作用的實驗研究［J］.中國海洋平臺，2001，16（5）:1-6. CHEN X F，LI Y C，SUN D P，et al.An experimental study of wave acting on perforated caisson［J］.China Offshore Platform，2001，16（5）:1-6.

［2］陳雪峰.波浪與開孔沉箱的相互作用［D］.大連：大連理工大學(xué)，2003.

［3］張玉彬，孫大鵬，夏志盛，等.明基床上開孔沉箱水平力極值對應(yīng)的垂直力分析［J］.中國海洋平臺，2012（5）:12. ZHANG Y B，SUN D P，XIA Z S，et al.Analysis on total vertical forces of perforated caisson sitting on rubble mound foundation while total horizontal forces reach its peak［J］.China Offshore Platform，2012（5）:12.

［4］施曉迪，琚烈紅，馮衛(wèi)兵，等.開孔沉箱消浪性能試驗研究［J］.水運工程，2011（3）:4. SHI X D，JU L H，F(xiàn)ENG W B，et al.Experimental research on wave attenuation effectiveness of perforated caission［J］.Port& waterway engineering，2011（3）:4.

［5］朱大同.開孔沉箱防波堤的波浪反射系數(shù)［C］//中國海洋工程學(xué)會.第十二屆中國海岸工程學(xué)術(shù)討論會論文集.北京:海洋出版社，2005.

［6］Ashgriz N，Barbat T，Wang G.A computational Lagrangian?Eulerian advection remap for free surface flows［J］.International Jour?nal for Numerical Methods in Fluids，2004，44（1）:1-32.

［7］呂林.海洋工程中小尺度物體的相關(guān)水動力數(shù)值計算［D］.大連：大連理工大學(xué)，2006.

［8］孫英偉.直立堤前的破碎波流場特征及波浪力研究［D］.大連：大連理工大學(xué)，2012.

［9］Dong C M，Huang C J.Generation and propagation of water waves in a two?dimensional numerical viscous wave flume［J］.Journal of waterway，port，coastal，and ocean engineering，2004，130（3）:143-153.

［10］Larsen J，Dancy H.Open boundaries in short wave simulations?a new approach［J］.Coastal Engineering，1983，7（3）:285-297.

［11］合田良實.港工建筑物的防浪設(shè)計（中譯本）［M］.北京：海洋出版社，1983.

全球最大波浪實驗水槽首次實驗成功

本刊從交通運輸部天津水運工程科學(xué)研究院獲悉，2014年7月29日，現(xiàn)今全球最大規(guī)模的大比尺波浪水槽成功進行了首次實驗——新型雙箱浮式防波堤結(jié)構(gòu)穩(wěn)定和消浪效果原型實驗。該大型水動力實驗中心大比尺波浪水槽長450 m、寬5 m、深8～12 m，能產(chǎn)生3.5 m的波浪和20 m3/s的水流，能進行1:5到1:1的大比尺模型試驗，是目前世界上尺度最大、造波能力最強的波浪試驗水槽。大型水動力實驗中心將建設(shè)成為具有國際領(lǐng)先水平的中國水運工程科研創(chuàng)新平臺。大比尺波浪水槽通過接近原型的實驗研究，克服比尺效應(yīng)，突破水運工程建設(shè)中涉及建筑物安全、生態(tài)環(huán)保、防災(zāi)減災(zāi)的水流、波浪、泥沙、地基、工程結(jié)構(gòu)等基礎(chǔ)理論技術(shù)制約，形成強大的自主創(chuàng)新能力。該次實驗按照1∶1比尺進行了波高1.5 m的驗證性實驗和波高2.25 m的破壞性實驗。（殷缶，梅深）

Biography：XING Fei（1990-），male，master student.

中圖分類號：TV 139.2；O 242.1

文獻標(biāo)識碼：A

文章編號：1005-8443（2014）05-0497-06

收稿日期：2013-12-27；修回日期：2014-02-27

基金項目：國家自然科學(xué)基金（51279027）；國家自然科學(xué)基金（51221961）

作者簡介：邢飛（1990-），男，山東省濰坊市人，碩士研究生，主要從事波浪與結(jié)構(gòu)物相互作用的研究。

Numerical simulation of wave interaction with perforated caisson breakwaters by CLEAR?VOF method

XING Fei，SUN Da?peng，XIA Zhi?sheng，WU Hao
（State Key Laboratory of Coastal and Offshore Engineering，Dalian University of Technology，Dalian 116023，China）

Abstract:A two?dimensional numerical flume based on the Navier?Stokes equations which were discretized with a three?step finite element method was developed to simulate water wave problems.In this numerical model，large eddy simulation method was used to describe the turbulent effect and computational Lagrangian?Eulerian ad?vection remap and the volume of the fluid method（CLEAR?VOF）was used to advect free surfaces.The reflection coefficient of perforated caisson and the ratio of the total horizontal force acting on the perforated caisson to them acting on the solid caisson had been studied under linear regular waves by using this numerical method.Through the comparison with physical model test，numerical calculation of other scholars and engineering project，it can be known that the numerical method may be useful for studying wave interaction with perforated caisson breakwaters.

Key words:three?step finite element method;CLEAR?VOF;perforated caisson;reflection coefficient;total horizontal force