規(guī)則波與流相互作用的數(shù)值模擬與不確定度分析

姚順， 馬寧， 丁俊杰， 顧解忡

(1.上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240；2. 上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海 200240)

在實(shí)際海洋環(huán)境中，波流相互作用是一種非常普遍的現(xiàn)象，流的存在會(huì)使波浪的波高、波長(zhǎng)等運(yùn)動(dòng)學(xué)特性發(fā)生顯著變化。依照波浪理論，疊加規(guī)則波可以得到不規(guī)則波等實(shí)際波浪，可見研究規(guī)則波與流相互作用可為實(shí)際波浪與流相互作用研究打下基礎(chǔ)，對(duì)模擬實(shí)際海域具有重要意義。

研究人員通常結(jié)合數(shù)值模擬和試驗(yàn)的方法對(duì)波浪及波流相互作用問(wèn)題進(jìn)行研究。胡捍紅[1]選擇四階非線性薛定諤方程模擬了過(guò)渡水波在傳播過(guò)程的非線性演變過(guò)程；賈巖等[2]基于波浪模式SWAN和海流模式POM研究了臺(tái)風(fēng)中的波流相互作用；Gao等[3]和Liu等[4]分別利用商業(yè)軟件開發(fā)數(shù)值波浪水池，研究了規(guī)則波與水平圓柱相互作用問(wèn)題。通過(guò)試驗(yàn)?zāi)軌蛑庇^地驗(yàn)證模擬結(jié)果的準(zhǔn)確性，F(xiàn)ernando等[5]、Qi等[6]開展試驗(yàn)研究了波浪與順逆流的非線性相互作用；Soltanpour等[7]和Wei等[8]通過(guò)試驗(yàn)分別研究了波流相互作用下波浪在泥床上的耗散和橋梁基礎(chǔ)的水動(dòng)力響應(yīng)特性；Jiang等[9]利用CFD軟件模擬了不規(guī)則波與淹沒(méi)式防波堤相互作用，并將模擬結(jié)果與實(shí)驗(yàn)進(jìn)行了對(duì)比。然而，對(duì)比模擬結(jié)果與試驗(yàn)結(jié)果驗(yàn)證數(shù)值模型的有效性從本質(zhì)上說(shuō)是定性的。在第22屆國(guó)際船模拖曳水池會(huì)議(ITTC)阻力委員會(huì)上推出了定量評(píng)估數(shù)值模擬結(jié)果的具體方法，即不確定度分析方法[10]，至此相關(guān)學(xué)者開始關(guān)注如何定量評(píng)估計(jì)算流體力學(xué)數(shù)值模擬結(jié)果的可靠性。通常，數(shù)值模型的數(shù)值誤差只能估計(jì)得到，Richardson外推法和網(wǎng)格收斂指數(shù)[11]因此被廣泛應(yīng)用于估計(jì)誤差和數(shù)值不確定性。基于此Wu等[12]對(duì)船模阻力進(jìn)行了不確定度分析，通過(guò)分析結(jié)果發(fā)現(xiàn)有效控制近船模表面第1層網(wǎng)格高度并選擇合適湍流模型可大幅度降低船模阻力數(shù)值試驗(yàn)的不確定度。Zhu等[13]和Deng等[14]分別對(duì)船模橫搖運(yùn)動(dòng)和小水線面雙體船在波浪中的縱向運(yùn)動(dòng)進(jìn)行了數(shù)值模擬與不確定度分析，給出了減小船舶運(yùn)動(dòng)參數(shù)數(shù)值模擬不確定度的有效建議。Silva等[15]、柏君勵(lì)等[16]則對(duì)規(guī)則波和聚焦波的波參數(shù)進(jìn)行了不確定分析。然而，文獻(xiàn)[1-16]沒(méi)有研究流對(duì)波浪模型不確定度的影響，這一問(wèn)題的難點(diǎn)在于開展波流相互作用試驗(yàn)時(shí)很難保證造波精度。隨著丁俊杰等[17]對(duì)上海交通大學(xué)風(fēng)洞循環(huán)水槽消波裝置的改進(jìn)，該水槽能很好的滿足試驗(yàn)要求。本文對(duì)均勻流作用下規(guī)則波的生成與演化進(jìn)行了物理試驗(yàn)，然后基于RANS方程建立數(shù)值波浪水池進(jìn)行波流相互作用下規(guī)則波特性研究，并依照ITTC規(guī)程開展計(jì)算結(jié)果的不確定度分析。給出了規(guī)則波的運(yùn)動(dòng)學(xué)特性隨流速變化的一般規(guī)律，并討論了順流對(duì)波流相互作用數(shù)值模型不確定度的影響情況。

1 數(shù)值模型及驗(yàn)證

1.1 試驗(yàn)安排

本文的試驗(yàn)在上海交通大學(xué)風(fēng)洞循環(huán)水槽中進(jìn)行，該循環(huán)水槽裝置的工作原理如圖1所示。裝置整體長(zhǎng)24.6 m，寬4.5 m，高8.5 m，主要分為工作段、上游整流段、下游整流段、底部回流段以及動(dòng)力段等部分。其中動(dòng)力段能使循環(huán)水槽中的水體沿順時(shí)針?lè)较蜻M(jìn)行循環(huán)流動(dòng)。工作段長(zhǎng)為8.0 m，寬為3.0 m，高1.95 m，工作段的水深根據(jù)循環(huán)水槽工作狀態(tài)的不同可小幅度調(diào)節(jié)，水深范圍為1.60～1.63 m。安裝于工作段前方的搖擺式造波機(jī)，可以生成規(guī)則波、JONSWAP譜不規(guī)則波與聚焦波等多類型波浪。實(shí)際造波時(shí)，根據(jù)試驗(yàn)需要設(shè)置造波機(jī)角度振幅和周期，即可造出目標(biāo)波浪。當(dāng)在靜水中造波時(shí)水槽水深固定為1.63 m，而單獨(dú)造流或同時(shí)造波和造流時(shí)水深固定為1.6 m。循環(huán)水槽能夠?qū)崿F(xiàn)長(zhǎng)時(shí)間穩(wěn)定的造流，穩(wěn)定造流的流速范圍為0.1～3.0 m/s。

試驗(yàn)的目的是獲取無(wú)流及有流情況下規(guī)則波自由液面處波高時(shí)歷曲線，驗(yàn)證波流相互作用數(shù)值模型有效性，并分析流對(duì)規(guī)則波運(yùn)動(dòng)學(xué)特性的影響。試驗(yàn)布置如圖2所示，其中消波裝置距造波端6 m。為滿足試驗(yàn)需要，在水槽工作段共布置5個(gè)浪高儀，1#、2#和3#浪高儀布置在消波裝置前側(cè)，分別距造波端3.2、3.6、4 m。為了減小前方浪高儀對(duì)后方浪高儀的影響，1#、2#和3#浪高儀分別距水槽左端1、2、1.5 m；4#和5#浪高儀平行布置在消波裝置后側(cè)，距造波端6.7 m，且分別距左側(cè)槽壁1、2 m，用于測(cè)試消波裝置的性能。實(shí)驗(yàn)開始時(shí)，造波機(jī)生成的前5～10個(gè)波并不穩(wěn)定，記錄數(shù)據(jù)時(shí)將前10 s時(shí)間段內(nèi)的波浪忽略。所有浪高儀量程為30 cm，測(cè)量頻率為100 Hz，精度為0.5%。

圖2 試驗(yàn)布置示意Fig.2 The arrangement of wave experiments

1.2 數(shù)值模型

本文利用計(jì)算流體力學(xué)商業(yè)軟件STARCCM+建立數(shù)值模型[18]，控制方程為Navier-Stokes方程：

(1)

式中：ui為xi方向的速度分量；ρ為流體的密度；p為壓強(qiáng)；μ為動(dòng)力粘度；gi=g是重力加速度。本文關(guān)注湍流引起的流場(chǎng)在時(shí)間上的平均變化，主要采用雷諾平均的Navier-Stokes(RANS)方程：

(2)

(3)

式中：μt是湍流黏度；τij是Kronecker函數(shù)；k是湍動(dòng)能。可采用Renormalization-Group(RNG)k-ε模型求解湍流黏度和湍動(dòng)能，從而得到雷諾應(yīng)力。

本文依照上海交通大學(xué)風(fēng)洞循環(huán)水槽構(gòu)建二維數(shù)值波浪水池。數(shù)值水池如圖3所示，主要分為工作段與消波段。由于只考慮波浪沿x方向的傳播，數(shù)值水池的寬設(shè)為0.02 m，小于最小網(wǎng)格y方向的尺寸，采樣點(diǎn)位置與試驗(yàn)中浪高儀位置對(duì)應(yīng)。

圖3 數(shù)值波浪水池Fig.3 Numerical wave flume

網(wǎng)格劃分時(shí)自由液面作為液相與氣相交界面，造波過(guò)程中會(huì)發(fā)生明顯的起伏運(yùn)動(dòng)。通常對(duì)自由液面處的網(wǎng)格進(jìn)行多層加密，從而減少造波誤差。計(jì)算域x、y、z方向?qū)?yīng)的最小網(wǎng)格尺寸分別為0.01、0.1、0.002 m。主要采集x=4 m自由液面處波高信息作為輸出數(shù)據(jù)，采樣頻率為100 Hz。

造波方法采用STARCCM+中造波模塊的源項(xiàng)造波法；消波時(shí)，在尾部2 m添加阻尼項(xiàng)模擬阻尼消波段，并在阻尼消波段前加入一種多層變角度開孔折彎板透水消波裝置[19]一同構(gòu)成消波區(qū)。該消波板能夠保證較高的透水效率的同時(shí)實(shí)現(xiàn)高效的消波。

采用STARCCM+庫(kù)函數(shù)CURRENT設(shè)定水體速度造流。數(shù)值模型空間離散采用二階迎風(fēng)格式，壓力和速度的耦合求解采用SIMPLE算法。自由液面捕捉采用流體體積(VOF)法。入口AC邊界選擇速度入口，出口BD邊界為壓力出口，底部AB邊界為無(wú)滑移固壁邊界而頂部邊界CD為壓力入口邊界，兩側(cè)的邊界選用對(duì)稱邊界。

按照Stokes波浪理論生成規(guī)則波，Stokes波浪理論與線性波類似，波浪運(yùn)動(dòng)也是勢(shì)運(yùn)動(dòng)。對(duì)于弱非線性問(wèn)題，可用攝動(dòng)法進(jìn)行求解，先假設(shè)勢(shì)函數(shù)φ和波面曲線函數(shù)η都是某個(gè)小參數(shù)εs的冪級(jí)數(shù)，則：

(4)

(5)

當(dāng)n=2，5時(shí)，代入自由表面條件泰勒展開式中，即可得到弱非線性的二階和五階Stokes波面方程和勢(shì)函數(shù)。本文數(shù)值模擬的研究對(duì)象為Stokes五階波，在軟件中選擇波浪類型并輸入指定波高和周期，即可得到目標(biāo)波浪。試驗(yàn)及模擬工況如表1所示，其中C為流速，H為波高，L為波長(zhǎng)，T為周期，EXP和NUM分別表示對(duì)該工況進(jìn)行試驗(yàn)或者數(shù)值模擬。每個(gè)工況重復(fù)3組，每組試驗(yàn)或者模擬時(shí)長(zhǎng)為50 s。

表1 規(guī)則波試驗(yàn)與數(shù)值模擬工況

1.3 數(shù)值模型驗(yàn)證

圖4 工況R1數(shù)值模擬結(jié)果的驗(yàn)證Fig.4 Verification of simulation results of case R1

表2 數(shù)值模擬與試驗(yàn)數(shù)據(jù)的平均波高對(duì)比Table 2 Experimental and numerical results of wave height

2 結(jié)果及不確定度分析

2.1 結(jié)果分析

圖5給出了無(wú)流及順流情況下規(guī)則波自由表面分布情況。從圖5可知，C=0 m/s時(shí)規(guī)則波的自由表面比較光滑；而在順流(C=0.3 m/s)作用下，規(guī)則波波浪自由表面變得粗糙。

數(shù)值模擬得到的工況R1～R5波高時(shí)歷曲線如圖6所示。圖6表明，隨著均勻順流流速增加，規(guī)則波波高降低，波峰和波谷相比于C=0 m/s時(shí)更平坦，逆流對(duì)規(guī)則波的影響與順流工況的結(jié)果相反。

圖5 數(shù)值水池中波浪自由表面情況Fig.5 The water surface of waves in the numerical wave flume

圖6 波高時(shí)歷數(shù)值模擬結(jié)果Fig.6 Simulation results of the water surface elevations

從圖6也發(fā)現(xiàn)當(dāng)逆流流速達(dá)到-0.3 m/s時(shí)，規(guī)則波傳播變得緩慢，落后于無(wú)流及順流工況；逆流流速過(guò)大時(shí)(C=-0.6 m/s)，規(guī)則波波高幾乎為0，這是由于當(dāng)逆流流速過(guò)大時(shí)，規(guī)則波在傳播時(shí)會(huì)在某個(gè)位置被阻隔。圖7給出了C=-0.6 m/s時(shí)波浪被阻隔的情況，被阻隔后規(guī)則波波面趨于平穩(wěn)。

圖7 波浪阻隔現(xiàn)象Fig.7 The wave blocking

2.2 不確定度分析

之前對(duì)數(shù)值模型的驗(yàn)證從本質(zhì)上來(lái)說(shuō)是定性的，無(wú)法定量評(píng)估數(shù)值模擬的誤差。本節(jié)擬采用ITTC推薦的不確定度分析方法，定量評(píng)估數(shù)值模擬的誤差。ITTC推薦的CFD數(shù)值模型的不確定度分析包括驗(yàn)證和確認(rèn)2個(gè)部分，其中評(píng)估數(shù)值不確定度的過(guò)程稱為驗(yàn)證，評(píng)估模型不確定度的過(guò)程稱為確認(rèn)，具體過(guò)程在ITTC規(guī)程中有詳細(xì)介紹[6]。

評(píng)估數(shù)值不確定度時(shí)，需要對(duì)數(shù)值模型的參數(shù)進(jìn)行收斂性分析。通常考慮的參數(shù)包括迭代次數(shù)、網(wǎng)格尺寸、時(shí)間步長(zhǎng)和其他參數(shù)，本文規(guī)則波與流相互作用屬于非穩(wěn)態(tài)問(wèn)題，網(wǎng)格尺寸和時(shí)間步長(zhǎng)引起的誤差比迭代次數(shù)等其他參數(shù)引起的誤差要大，因此本節(jié)的不確定度分析將網(wǎng)格尺寸和時(shí)間步長(zhǎng)作為主要因素加以考慮。

首先對(duì)工況R1中的網(wǎng)格參數(shù)進(jìn)行收斂性研究，時(shí)間步長(zhǎng)為0.001 s。用3套網(wǎng)格分別進(jìn)行數(shù)值模擬后，得到3組波高時(shí)歷曲線，如圖8所示。

圖8 不同網(wǎng)格對(duì)應(yīng)的波高時(shí)歷曲線Fig.8 Water surface elevations simulated with different grids

相鄰網(wǎng)格參數(shù)波高時(shí)歷之差的平均L2范數(shù)為：

(6)

(7)

其中ηG1i、ηG2i、ηG3i分別表示在細(xì)網(wǎng)格、中等網(wǎng)格和粗網(wǎng)格上模擬得到的波高時(shí)間序列中的第i個(gè)波高值，全局收斂因子為：

(8)

(9)

(10)

網(wǎng)格收斂因子：

(11)

式中pG,est是極限階數(shù)的估計(jì)值，一般取2。由網(wǎng)格尺寸導(dǎo)致的數(shù)值模擬誤差為：

(12)

當(dāng)CG遠(yuǎn)小于或遠(yuǎn)大于1時(shí)，網(wǎng)格尺寸引起的數(shù)值模擬不確定度為：

(13)

當(dāng)|1-CG|≥0.25，網(wǎng)格尺寸引起的數(shù)值模擬不確定度為：

(14)

接下來(lái)在細(xì)網(wǎng)格Grid1上進(jìn)行時(shí)間步長(zhǎng)收斂性分析，得到3個(gè)時(shí)間步長(zhǎng)對(duì)應(yīng)的各個(gè)波高時(shí)歷曲線，如圖9所示。

圖9 不同時(shí)間步長(zhǎng)對(duì)應(yīng)的波高時(shí)歷Fig.9 Water surface elevations simulated with different time steps

相鄰時(shí)間步長(zhǎng)參數(shù)的波高時(shí)歷之差的平均L2范數(shù)為：

(15)

(16)

其中ηT1i、ηT2i、ηT3i分別表示時(shí)間步長(zhǎng)0.001 s、0.001 41 s、0.002 s時(shí)數(shù)值模擬得到的波高時(shí)間序列中的第i個(gè)波高值,故全局收斂因子：

(17)

(18)

(19)

pT,est同樣取2，時(shí)間步長(zhǎng)修正因子：

(20)

由時(shí)間步長(zhǎng)導(dǎo)致的數(shù)值模擬誤差為：

(21)

時(shí)間步長(zhǎng)導(dǎo)致的數(shù)值模擬不確定度：

(22)

綜上，在x=4 m自由液面處t∈[15,25]s時(shí)間段內(nèi)波高時(shí)歷的數(shù)值模擬不確定度為：

(23)

數(shù)值模擬誤差估計(jì)值為：

(24)

表3 數(shù)值模擬驗(yàn)證結(jié)果Table 3 Numerical results of verification

表3給出，R1和R22個(gè)工況數(shù)值模擬的網(wǎng)格不確定度都大于時(shí)間步長(zhǎng)不確定度，可知數(shù)值模擬無(wú)流及順流情況下的規(guī)則波運(yùn)動(dòng)時(shí)，數(shù)值模擬結(jié)果對(duì)網(wǎng)格尺寸依賴程度大于時(shí)間步長(zhǎng)，實(shí)際模擬時(shí)計(jì)算資源有限情況下應(yīng)盡量滿足網(wǎng)格精度要求。對(duì)比R1和R2工況數(shù)據(jù)可知，順流工況R2數(shù)值模擬的網(wǎng)格不確定度是無(wú)流工況R1的3.46倍，時(shí)間步長(zhǎng)不確定度是無(wú)流時(shí)的1.60倍，可見相比于無(wú)流工況R1，順流提高了規(guī)則波時(shí)歷曲線數(shù)值模擬結(jié)果對(duì)網(wǎng)格尺寸和時(shí)間步長(zhǎng)的依賴程度。并且，順流時(shí)規(guī)則波的平均波高小于無(wú)流工況，考慮順流與規(guī)則波相互作用問(wèn)題時(shí)，可適當(dāng)減少自由液面加密區(qū)的高度，但需要注意選擇合適的網(wǎng)格尺寸和時(shí)間步長(zhǎng)，減小數(shù)值模擬帶來(lái)的誤差。

接下來(lái)根據(jù)ITTC規(guī)程，對(duì)數(shù)值模擬結(jié)果進(jìn)行確認(rèn)，工況R1中x=4 m自由液面處模擬的波高時(shí)歷曲線確認(rèn)不確定度UV為：

(25)

(26)

式中ηEXP1,x=4 m和ηNUM1,x=4 m分別表示工況R1在x=4 m自由液面處波高時(shí)歷曲線的試驗(yàn)和模擬值。

表4 數(shù)值模擬確認(rèn)結(jié)果Table 4 Numerical results of validation

3 結(jié)論

1)均勻順流使規(guī)則波波高降低，傳播速度加快，波峰和波谷與無(wú)流工況相比更平坦；逆流對(duì)規(guī)則波的影響與順流的影響相反，逆流流速過(guò)大時(shí)，規(guī)則波會(huì)被阻隔。

2)無(wú)流和順流工況規(guī)則波波高時(shí)歷曲線的比較誤差都小于對(duì)應(yīng)的確認(rèn)不確定度，數(shù)值模擬結(jié)果得到確認(rèn)，構(gòu)建的數(shù)值波浪水池具有可靠性。

3)無(wú)流及順流情況下規(guī)則波的數(shù)值模擬結(jié)果對(duì)網(wǎng)格尺寸依賴程度大于時(shí)間步長(zhǎng)，實(shí)際模擬時(shí)計(jì)算資源有限情況下應(yīng)盡量滿足網(wǎng)格精度要求。

4)相比無(wú)流工況，順流加劇了規(guī)則波波高時(shí)歷曲線計(jì)算結(jié)果對(duì)網(wǎng)格尺寸和時(shí)間步長(zhǎng)的依賴程度。

目前在循環(huán)水槽中進(jìn)行逆流與規(guī)則波相互作用的試驗(yàn)時(shí)很難控制造波精度，下一步將從改善造波板水力傳遞函數(shù)等方面提高逆流中的造波質(zhì)量，開展逆流與規(guī)則波相互作用的不確定度分析。