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

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

(1.上海交通大學 海洋工程國家重點實驗室，上海 200240；2. 上海交通大學 船舶海洋與建筑工程學院，上海 200240)

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

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

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

1.1 試驗安排

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

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

圖2 試驗布置示意Fig.2 The arrangement of wave experiments

1.2 數(shù)值模型

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

(1)

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

(2)

(3)

式中：μt是湍流黏度；τij是Kronecker函數(shù)；k是湍動能?？刹捎肦enormalization-Group(RNG)k-ε模型求解湍流黏度和湍動能，從而得到雷諾應力。

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

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

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

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

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

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

(4)

(5)

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

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

1.3 數(shù)值模型驗證

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

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

2 結果及不確定度分析

2.1 結果分析

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

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

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

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

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

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

2.2 不確定度分析

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

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

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

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

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

(6)

(7)

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

(8)

(9)

(10)

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

(11)

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

(12)

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

(13)

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

(14)

接下來在細網(wǎng)格Grid1上進行時間步長收斂性分析，得到3個時間步長對應的各個波高時歷曲線，如圖9所示。

圖9 不同時間步長對應的波高時歷Fig.9 Water surface elevations simulated with different time steps

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

(15)

(16)

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

(17)

(18)

(19)

pT,est同樣取2，時間步長修正因子：

(20)

由時間步長導致的數(shù)值模擬誤差為：

(21)

時間步長導致的數(shù)值模擬不確定度：

(22)

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

(23)

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

(24)

表3 數(shù)值模擬驗證結果Table 3 Numerical results of verification

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

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

(25)

(26)

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

表4 數(shù)值模擬確認結果Table 4 Numerical results of validation

3 結論

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

2)無流和順流工況規(guī)則波波高時歷曲線的比較誤差都小于對應的確認不確定度，數(shù)值模擬結果得到確認，構建的數(shù)值波浪水池具有可靠性。

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

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

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