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

（大連理工大學(xué)海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室，大連116023）

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

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

（大連理工大學(xué)海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室，大連116023）

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

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

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

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ī)則波作用下開孔沉箱反射率和波浪總水平力的計(jì)算結(jié)果，為研究波浪與開孔沉箱相互作用探索了一種新的數(shù)值模式。

1　數(shù)值計(jì)算方法

1.1控制方程及其離散

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

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

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

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

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

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

至此，方程組封閉?？刹捎肎alerkin加權(quán)余量法進(jìn)行離散求解。

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

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

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

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

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

2.1線性波浪模擬

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

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

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

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

式中：us為海綿吸收層內(nèi)單元節(jié)點(diǎn)的速度；μ（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，計(jì)有8種波要素組合。波浪傳播穩(wěn)定后，各工況的線性規(guī)則波浪與理論解吻合良好。限于篇幅，本文以波高0.10 m、周期1.0 s工況為例給出了水槽工作區(qū)末端位置的波高歷時(shí)曲線，如圖3所示，圖中實(shí)線表示理論值，圓圈表示本文計(jì)算值。

圖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個(gè)孔，開孔率為23.8%，消浪室寬度為0.3 m，消浪室內(nèi)的水深為0.2 m。

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

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

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

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

3.1反射率的比較

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

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

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

3.2水平力的比較

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

圖6、圖7分別給出了波高為0.10 m和0.12 m工況下（a）波峰總水平力比的比較和（b）波谷總水平力比的比較。其中橫坐標(biāo)為波長，縱坐標(biāo)為力比；實(shí)線表示經(jīng)驗(yàn)公式（17）或經(jīng)驗(yàn)公式（18）的結(jié)果，圓點(diǎn)表示本文的數(shù)值計(jì)算結(jié)果，三角形點(diǎn)表示陳雪峰［2］的數(shù)值計(jì)算結(jié)果。由圖可知，本文計(jì)算的波峰（谷）總水平力與陳雪峰數(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　工程實(shí)例驗(yàn)證

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

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

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

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

5　結(jié)語

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

參考文獻(xiàn)：

［1］陳雪峰，李玉成，孫大鵬，等.波浪與開孔沉箱作用的實(shí)驗(yàn)研究［J］.中國海洋平臺(tái)，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］張玉彬，孫大鵬，夏志盛，等.明基床上開孔沉箱水平力極值對(duì)應(yīng)的垂直力分析［J］.中國海洋平臺(tái)，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)兵，等.開孔沉箱消浪性能試驗(yàn)研究［J］.水運(yùn)工程，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é)會(huì).第十二屆中國海岸工程學(xué)術(shù)討論會(huì)論文集.北京:海洋出版社，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)水動(dòng)力數(shù)值計(jì)算［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í).港工建筑物的防浪設(shè)計(jì)（中譯本）［M］.北京：海洋出版社，1983.

全球最大波浪實(shí)驗(yàn)水槽首次實(shí)驗(yàn)成功

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

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

中圖分類號(hào)：TV 139.2；O 242.1

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1005-8443（2014）05-0497-06

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

基金項(xiàng)目：國家自然科學(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