三維數(shù)值波流水槽的構(gòu)建及驗(yàn)證

崔栗銘

(西南交通大學(xué)土木工程學(xué)院， 四川成都 610031)

在跨海橋梁的建設(shè)過(guò)程中，往往面臨著復(fù)雜的海洋環(huán)境，將會(huì)遇到風(fēng)、浪、流的嚴(yán)峻考驗(yàn)。傳統(tǒng)的波浪及水流對(duì)結(jié)構(gòu)物的作用已有大量學(xué)者展開(kāi)研究并得到了計(jì)算方法。而波流聯(lián)合作用的機(jī)理十分復(fù)雜，并不能簡(jiǎn)單地將兩者進(jìn)行迭加而得到。計(jì)算機(jī)數(shù)值模擬的方法由于其不需實(shí)驗(yàn)費(fèi)用，且可視化程度較高，因此很多學(xué)者采用數(shù)值方法對(duì)流體作用進(jìn)行研究。其中T·Kawamura[2]運(yùn)用大渦模擬研究了穿透液面的圓柱繞流；Jungsoo Suh[3]研究了氣液交界面對(duì)圓柱繞流漩渦脫離的影響；康啊真[4]建立了三維波浪與結(jié)構(gòu)物相互作用的數(shù)學(xué)模型；馮誼武[5]通過(guò)CFX研究了圍堰所受波浪荷載。目前關(guān)于波流場(chǎng)的數(shù)值模擬，主要方法有利用邊界條件造波造流法、質(zhì)量源法。李勇[6]基于RANS方程，使用邊界造波和邊界入流相迭加的方法實(shí)現(xiàn)數(shù)值波流場(chǎng)的模擬；康啊真[7]基于LES建立了三維波流與結(jié)構(gòu)物相互作用的數(shù)值模型；Li[8]建立了二維數(shù)值水槽來(lái)計(jì)算波流作用下淹沒(méi)圓柱的水動(dòng)力系數(shù)。

本文使用數(shù)值計(jì)算軟件Fluent，基于Reynolds平均的Navier-Stocks方程，以及速度入口造波造流的方法模擬波流聯(lián)合作用，并使用動(dòng)量源項(xiàng)對(duì)下游波流場(chǎng)進(jìn)行消波，通過(guò)考察波浪傳播的完整性和規(guī)則性，以及沿水深的時(shí)均速度分布驗(yàn)證了數(shù)值水槽的合理性。

1 數(shù)值計(jì)算理論

現(xiàn)有的流體數(shù)值計(jì)算理論基礎(chǔ)主要包括兩類，一是基于勢(shì)流理論，二是基于粘性流體理論。其中基于勢(shì)流理論的邊界元法被廣泛應(yīng)用于波浪與浮體的相互作用問(wèn)題，但是勢(shì)流方法忽略了流體粘性，而本文需要模擬波流聯(lián)合的結(jié)構(gòu)物作用，因此采用基于粘流理論的計(jì)算軟件Fluent，并基于雷諾平均的NS方程(RANS)建立波流聯(lián)合作用數(shù)值水槽。其控制方程主要如下。

雷諾平均方程的質(zhì)量守恒方程以及動(dòng)量方程分別如下：

(1)

(2)

對(duì)于湍流較強(qiáng)的數(shù)值模型，在Fluent中常采用k-ω或k-ε模型，傳統(tǒng)的兩方程模型k-ω模型更適用于存在逆壓梯度的邊界層問(wèn)題，而k-ε模型能更好地模擬遠(yuǎn)離壁面處的充分發(fā)展的湍流流動(dòng)。k-ωSST湍流模型結(jié)合了k-ε和k-ω湍流模型兩者的優(yōu)勢(shì)，在近壁面處保留原始的k-ω，在遠(yuǎn)離壁面處采用k-ε模型。k-ωSST湍流模型的控制方程如下：

(3)

=αS2-βω2-(F1-1)CDkω

(4)

(5)

(6)

(7)

湍流粘度定義為：

(8)

S為應(yīng)變速率的不變量：

G=2vTSijSij

(9)

F2為混合函數(shù)：

(10)

αi=α1F1+α2(1-F1)

(11)

βi=β1F1+β2(1-F1)

(12)

其中系數(shù)αk1=0.85034,αk2=1,αω1=0.5,αω2=0.85616,γ1=0.5532,γ2=0.4403,β1=0.075,β2=0.0828,β*=0.09。

2 計(jì)算方法及初始條件

Fluent采用有限體積法對(duì)控制方程進(jìn)行離散，本文采用基于壓力求解的瞬態(tài)計(jì)算方法，瞬態(tài)格式采用一階隱式，湍動(dòng)能和紊動(dòng)能耗散率采用一階迎風(fēng)，并采用PISO算法進(jìn)行迭代計(jì)算。對(duì)于自由液面的追蹤，采用流體體積方法(VOF)，計(jì)算初始條件將水體定義體積分?jǐn)?shù)為1，水面以上部分體積分?jǐn)?shù)定義為0，純流及波流工況下，定義整個(gè)水體部分初始流速Uo。對(duì)于純波工況，除了湍動(dòng)能和紊動(dòng)能耗散率定義為小值，其余均為零。

3 邊界條件的設(shè)置

為了建立水槽模擬波流聯(lián)合作用，本文采用速度入口方法進(jìn)行造流與造波，即在水槽的上游邊界通過(guò)用戶自定義函數(shù)輸入波流聯(lián)合作用下的速度函數(shù)來(lái)造波流，在下游通過(guò)明渠出口邊界條件使水流自由流出，其中各個(gè)邊界的設(shè)置及物理意義分別如下。

兩側(cè)邊界：采用對(duì)稱(Symmetry)邊界條件，沿該邊界法向，速度為零，所有物理量梯度為零。

頂部邊界：采用壓力出口(Pressure)邊界條件，滿足動(dòng)力學(xué)邊界條件，表面壓力為大氣壓力Pa。

底部邊界：由于隨著時(shí)間變化，底部邊界的紊動(dòng)會(huì)從邊壁向上部擴(kuò)散， 形成充分的紊流，因此底部邊條件設(shè)置為不可滑移墻邊界(Wall)，壁面粗糙度為常數(shù)0.5。

上游邊界：對(duì)于波流工況，采用速度入口邊界(Velocity-inlet)輸入波流的速度分布函數(shù)，設(shè)置x，y方向速度為波流作用下的水平速度ux和豎向速度uy，z方向的速度設(shè)置為0，同時(shí)由于入口邊界的液面會(huì)隨波面的變化而變化，因此還需要自定義液面體積分?jǐn)?shù)函數(shù)Wavesurface。波流作用的速度ux，uy，以及液面體積分?jǐn)?shù)函數(shù)Wavesurface公式如下所示：

(13)

(14)

wavesurface=H/2·sin(k·x-σ·t)

(15)

其中：Uo為純流流速；H為波高；k為波數(shù)；d為水深；σ為波頻；z為沿水深方向高度；C為波速。

下游邊界：將下游出口邊界設(shè)置為壓力出口(Pressure-outlet)，并打開(kāi)明渠流設(shè)置(Openchannel)，邊界的壓力判別方法采用自由液面(FreeSurfaceLevel)判別，出口邊界的體積分?jǐn)?shù)采用相鄰單元的體積分?jǐn)?shù)，邊界的靜壓及總壓計(jì)算方法如下：

(16)

(17)

源項(xiàng)消波：為了避免下游波浪反射影響，因此設(shè)置消波條件，采用源項(xiàng)消波，通過(guò)在消波段定義動(dòng)量源x_source和y_source來(lái)進(jìn)行消波，其中：

con=C_R(c,t)·(x[0]-x1)/(x2-x1)·θ1

(18)

純波：x_source=-C_U(c,t)·con

(19)

波流：x_source=(C_U(c,t)-Uo)·(-con)

(20)

y_source=-C_V(c,t)·con

(21)

C_R(c,t)為水槽消波段任意水質(zhì)點(diǎn)的密度；Uo為波流工況中均勻流流速；C_U(c,t)；C_V(c,t)為水槽消波段任意水質(zhì)點(diǎn)的x，y方向速度；源項(xiàng)消波阻尼參數(shù)θ1=2～4。

4 水槽品質(zhì)檢驗(yàn)

為了對(duì)數(shù)值水槽進(jìn)行驗(yàn)證，因此建立數(shù)值模型與其他學(xué)者所做數(shù)值水槽及實(shí)驗(yàn)進(jìn)行對(duì)比，本文數(shù)值驗(yàn)證采用與李勇文章中相同的工況進(jìn)行驗(yàn)證，計(jì)算域?yàn)閤yz=36m×0.8m×0.5m。模型邊界條件及網(wǎng)格劃分方式如圖1所示，分別在液面高度和靠近水槽底部處進(jìn)行網(wǎng)格加密，而在速度入口一個(gè)波長(zhǎng)范圍也同樣進(jìn)行網(wǎng)格加密，這樣做的目的分別是為了更好的捕捉液面，更好地模擬水槽底部的剪切作用使得沿水深方向的速度分布更加接近于實(shí)際，以及防止在造波邊界處出現(xiàn)波浪破碎。波浪工況采用李勇文章中的波流工況，同時(shí)也是Van Rijn[9]實(shí)驗(yàn)中的工況，工況波浪參數(shù)如表1所示。

圖1 模型邊界

表1 Van Rijn實(shí)驗(yàn)條件

考察數(shù)值水槽的品質(zhì)是否優(yōu)良,需要考察水槽波浪變化是否規(guī)則,波浪傳播是否完整,因此需要知道水槽瞬時(shí)的三維波面變化情況。根據(jù)Mehaute波浪理論選擇的判斷方法,判定工況中的波浪參數(shù)非線性較強(qiáng),且水深較淺,因此采用橢圓余弦波浪理論較為合適。從圖2中也可以看到,三維數(shù)值波面?zhèn)鞑ネ暾乙?guī)則,呈現(xiàn)波峰高窄而波谷緩長(zhǎng)的形式,說(shuō)明波流聯(lián)合作用下的波浪形態(tài)模擬較好。

圖2 三維瞬時(shí)波面

再對(duì)水槽波面變化情況進(jìn)行考察，由于波流水槽中需要滿足水槽下游反射波浪對(duì)上游無(wú)影響，且水槽內(nèi)流量守恒，因此需要在尾部進(jìn)行消波，并且需要波面無(wú)抬升。由圖3波流水槽波面變化圖可以看出，隨著波浪進(jìn)入消波區(qū)域，波高迅速衰減至靜水面位置，且無(wú)反射波的出現(xiàn)。同時(shí)在水槽中四個(gè)波浪的波高不隨位置的改變而衰減，說(shuō)明沒(méi)有發(fā)生數(shù)值耗散，在整個(gè)水槽中波浪都保持非常完整的形態(tài)。同時(shí)也可以看出，四個(gè)波高沒(méi)有產(chǎn)生依次的增高，說(shuō)明整個(gè)水槽內(nèi)流量保持守恒的狀態(tài)，在下游邊界，水流可以自由地流出。

圖3 波流水槽波面變化

再對(duì)數(shù)值水槽中的速度分布進(jìn)行考察，由于建立的波流水槽在消波前的區(qū)域?yàn)椴髀?lián)合作用，而在消波區(qū)為了消除波浪的反射等，因此需要通過(guò)源項(xiàng)消波使下游區(qū)域到出口邊界逐漸由波流聯(lián)合形態(tài)轉(zhuǎn)變?yōu)榧兞餍螒B(tài)。因此由圖4水槽速度云圖可以看出，在消波區(qū)域之前，由于波流作用，使得在波浪位置處，速度變化較為劇烈，能量主要聚集在四個(gè)波浪位置處，而進(jìn)入消波區(qū)域后,隨著波浪的消減，波面變得平靜。而在靠近下游邊界的地方，整個(gè)水槽深度方向速度變得一致，說(shuō)明由于動(dòng)量源項(xiàng)的作用，使得水體沿深度方向變成了速度一致的均勻流。

圖4 水槽速度云圖

最后再對(duì)水槽中不同位置處沿水深方向的速度分布進(jìn)行比較，將x=4 m、5.5 m、8 m、11 m、15 m、20 m位置處的速度進(jìn)行時(shí)間平均，然后與李勇建立的水槽速度分布以及其他學(xué)者實(shí)驗(yàn)中的數(shù)據(jù)進(jìn)行對(duì)比。由圖5可以看出，在中上游位置處，流速的分布比較穩(wěn)定，流速變化趨勢(shì)基本一致。在水深0.2 m以上位置，本文建立水槽速度分布介于純流以及波流實(shí)驗(yàn)之間，在水深0.05～0.2 m之間，本文與李勇所構(gòu)建水槽速度分布及波流實(shí)驗(yàn)數(shù)據(jù)較為一致。而在0.05 m位置以下，本文速度分布與李勇及純流實(shí)驗(yàn)數(shù)據(jù)較為一致；在靠近下游x=15 m、20 m位置處，可以看出流速相對(duì)中上游位置速度分布發(fā)生較大改變。在水深0.05 m以上位置，流速逐漸減小，而在底部位置，速度分布變化相對(duì)較小，整體更趨向于均勻分布，說(shuō)明在靠近于下游位置處，由于消波區(qū)域的影響，使得流速分布更加趨向于均勻流。

總體而言，本文所建立數(shù)值水槽在能夠保證波流聯(lián)合作用下波浪要素的完整性，并不發(fā)生數(shù)值耗散，同時(shí)水槽中流量守恒，并沒(méi)有出現(xiàn)波面的抬升。同時(shí)下游消波段能夠完整消波，使得消波區(qū)域的流體形態(tài)由波流轉(zhuǎn)變?yōu)榫鶆蛄?。整個(gè)水槽沿水深的速度來(lái)看，在中上游位置處，水槽速度分布介于純流及波流實(shí)驗(yàn)數(shù)據(jù)之間，在靠近下游消波段位置處，流速分布受其影響，更加趨近于均勻流，只要將結(jié)構(gòu)物放置在遠(yuǎn)離下游的位置處，可以忽略下游對(duì)上游的影響，因此本文所建立水槽能夠較好模擬波流聯(lián)合作用。

(a)x=4m

(b)x=5.5m

(c)x=8m

(d)x=11m

(e)x=15m

(f)x=20m圖5 x=4m、5.5m、8m、11m、15m、 20m處沿水深方向速度分布對(duì)比

5 結(jié)論

本文采用邊界造波和造流方法所得到的三維數(shù)值水槽模型，經(jīng)過(guò)液面變化形態(tài)考察和速度分布定量對(duì)比兩方面，驗(yàn)證了數(shù)值水槽消波的效果以及流場(chǎng)中速度分布與實(shí)際情況的吻合程度。從瞬時(shí)三維波面變化圖以及二維液面變化圖得出，水槽造波完整且穩(wěn)定，并且水槽內(nèi)部流量守恒，液面無(wú)抬升，消波效果迅速，在消波段內(nèi)，流場(chǎng)形態(tài)由波流聯(lián)合作用場(chǎng)迅速消減為均勻流場(chǎng)。從時(shí)均速度分布對(duì)比圖中可以看出，在中上游位置，流速分布介于純流實(shí)驗(yàn)及波流實(shí)驗(yàn)數(shù)據(jù)之間。總體而言水槽流速分布與實(shí)際情況吻合較好，靠近消波段位置的流速分布受消波效果影響，沿水深上部流速降低。因此在研究波流聯(lián)合作用時(shí)，應(yīng)考慮中上游位置的流場(chǎng)，且結(jié)構(gòu)物應(yīng)遠(yuǎn)離下游消波段。綜上所述，本文所建三維數(shù)值水槽對(duì)于波流聯(lián)合作用下的流場(chǎng)形態(tài)模擬較好，波浪傳播穩(wěn)定，參數(shù)完整，且無(wú)反射波的影響。