規則波與流相互作用的數值模擬與不確定度分析

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

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

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

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

1 數值模型及驗證

1.1 試驗安排

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

試驗的目的是獲取無流及有流情況下規則波自由液面處波高時歷曲線，驗證波流相互作用數值模型有效性，并分析流對規則波運動學特性的影響。試驗布置如圖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個波并不穩定，記錄數據時將前10 s時間段內的波浪忽略。所有浪高儀量程為30 cm，測量頻率為100 Hz，精度為0.5%。

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

1.2 數值模型

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

(1)

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

(2)

(3)

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

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

圖3 數值波浪水池Fig.3 Numerical wave flume

網格劃分時自由液面作為液相與氣相交界面，造波過程中會發生明顯的起伏運動。通常對自由液面處的網格進行多層加密，從而減少造波誤差。計算域x、y、z方向對應的最小網格尺寸分別為0.01、0.1、0.002 m。主要采集x=4 m自由液面處波高信息作為輸出數據，采樣頻率為100 Hz。

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

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

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

(4)

(5)

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

表1 規則波試驗與數值模擬工況

1.3 數值模型驗證

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

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

2 結果及不確定度分析

2.1 結果分析

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

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

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

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

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

圖7 波浪阻隔現象Fig.7 The wave blocking

2.2 不確定度分析

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

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

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

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

相鄰網格參數波高時歷之差的平均L2范數為：

(6)

(7)

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

(8)

(9)

(10)

網格收斂因子：

(11)

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

(12)

當CG遠小于或遠大于1時，網格尺寸引起的數值模擬不確定度為：

(13)

當|1-CG|≥0.25，網格尺寸引起的數值模擬不確定度為：

(14)

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

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

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

(15)

(16)

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

(17)

(18)

(19)

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

(20)

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

(21)

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

(22)

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

(23)

數值模擬誤差估計值為：

(24)

表3 數值模擬驗證結果Table 3 Numerical results of verification

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

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

(25)

(26)

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

表4 數值模擬確認結果Table 4 Numerical results of validation

3 結論

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

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

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

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

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