單鐵兵，楊建民，李 欣，肖龍飛

（上海交通大學(xué) 海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海200240）

1 引 言

隨著深海油氣田的不斷開發(fā)，淺海開發(fā)技術(shù)已難以滿足要求，傳統(tǒng)的重力式平臺(tái)和導(dǎo)管架平臺(tái)因其自重和造價(jià)隨水深增加而大幅度增加，已經(jīng)不適應(yīng)深海油氣發(fā)展的要求，取而代之的是浮式生產(chǎn)系統(tǒng)。其中半潛平臺(tái)具有性能優(yōu)良，運(yùn)動(dòng)響應(yīng)較小，工作水深范圍廣，抗風(fēng)浪能力強(qiáng)，甲板面積大及裝載能力強(qiáng)等優(yōu)點(diǎn)，在海洋石油勘探和開采中得到了廣泛的應(yīng)用[1]。隨著石油勘探開采朝深海以及超深海發(fā)展，復(fù)雜的深海環(huán)境條件以及近年來超強(qiáng)臺(tái)風(fēng)等極具破壞力的極端海況的頻繁出現(xiàn)，對平臺(tái)的安全性提出了嚴(yán)峻挑戰(zhàn)。

諸多的安全性評估指標(biāo)中，氣隙預(yù)報(bào)（海洋平臺(tái)下層甲板底部至波面間垂直距離的預(yù)報(bào)）一直是半潛平臺(tái)設(shè)計(jì)過程中一個(gè)關(guān)鍵問題[2]。在波浪與平臺(tái)的相互作用過程中，有一種現(xiàn)象特別引起關(guān)注，即波浪遇平臺(tái)立柱后將沿其表面往上攀升，有時(shí)其高度明顯大于其它區(qū)域的波浪，甚至能爬升至平臺(tái)下甲板，以致嚴(yán)重影響平臺(tái)的氣隙性能，該非線性現(xiàn)象稱為波浪爬升。波浪的爬升過程從理論上可解釋為當(dāng)入射波浪沖擊浮體時(shí)，水波的動(dòng)能轉(zhuǎn)化為勢能同時(shí)波高發(fā)生放大的過程。

導(dǎo)管架平臺(tái)由于相應(yīng)管架直徑小，對入射波擾動(dòng)小，故波浪作用機(jī)理簡單，采用線性理論即可對其周圍的波面分布進(jìn)行預(yù)報(bào)。對于大型的浮式結(jié)構(gòu)物如半潛平臺(tái)而言，因體型大，且下浮體與甲板之間通過若干立柱支撐，當(dāng)波陡較大時(shí)，波浪非線性效應(yīng)較強(qiáng)，柱體周圍波浪爬升現(xiàn)象尤為明顯，波高易發(fā)生放大，并表現(xiàn)出較強(qiáng)的非線性特征，甚至超過平臺(tái)下甲板發(fā)生越浪，此過程中產(chǎn)生的砰擊載荷將對平臺(tái)的安全性構(gòu)成威脅。陡波條件下的波浪爬升問題近年來已成為水動(dòng)力學(xué)研究的熱點(diǎn)問題之一，是平臺(tái)氣隙預(yù)報(bào)的一個(gè)重要方面。

很早以來就有學(xué)者對波浪的繞射和爬升效應(yīng)進(jìn)行了試探性的研究。McCamy等[3]1954年給出了單根圓柱周圍繞射波浪的經(jīng)典線性解。Niedzwecki等[4]1992年針對該圓柱上的波浪爬升問題開展了模型試驗(yàn)，試驗(yàn)結(jié)果表明，波陡（H/λ，其中H為波高，λ為波長）較大時(shí)，采用線性散射理論大大低估了波浪的爬升高度。Huston等[5]對一座四立柱的張力腿平臺(tái)（TLP）在單色規(guī)則波下進(jìn)行了一系列的試驗(yàn)測量，旨在研究波陡、立柱間的間隙等對立柱上的波浪爬坡和甲板下部波浪放大的影響。Lee和Newman等[6]1994年采用基于二階邊界元的方法對平臺(tái)周圍的波浪效應(yīng)進(jìn)行了模擬。盡管二階理論能夠提高波浪爬升預(yù)報(bào)的精確度，但仍不能滿足平臺(tái)詳細(xì)設(shè)計(jì)的要求，且無法捕捉到波浪沖擊立柱時(shí)的高階非線性現(xiàn)象[7-8]。

計(jì)算機(jī)的發(fā)展使平臺(tái)立柱周圍的波浪散射、海浪沿立柱爬升、甚至砰擊平臺(tái)甲板后波浪發(fā)生變形和翻滾等強(qiáng)非線性現(xiàn)象的直接數(shù)值模擬成為可能。與二階理論不同的是，CFD理論基于N-S方程直接求解時(shí)域內(nèi)的流場。近年來有學(xué)者針對一些自由液面大變形的強(qiáng)非線性問題提出了較為精確的追蹤重構(gòu)自由面的數(shù)值方法。其中常見的有Volume-of-Fluid（VOF）、Smooth-Particle-Hydrodynamics（SPH）、Level-Set（LS）以及CIP方法等。Donald[9]采用兩種方法求解重力式平臺(tái)周圍的波面分布，第一種為二階散射程序（WIMIT），基于勢流和小波陡假設(shè),采用截?cái)嗟臄z動(dòng)展開并結(jié)合邊界元理論求解波動(dòng)場，考慮二階非線性特征。第二種為基于VOF理論的CFD方法，該方法直接求解不可壓縮N-S方程，采用的是全非線性自由面邊界條件。通過對比發(fā)現(xiàn)，線性波條件下，二階散射理論可較為準(zhǔn)確地模擬平臺(tái)上的波浪爬升，但當(dāng)波陡較大時(shí)，二階攝動(dòng)理論限制了其預(yù)報(bào)的有效性。相比而言，VOF理論即使在入射波浪波陡較大時(shí)也能夠捕捉到波物相互作用時(shí)的高階非線性效應(yīng)。Roberto等[10]基于雷諾平均的N-S方程，使用僅求解水相的LS方法模擬自由液面，研究粘性效應(yīng)對立柱周圍波面分布的影響。

本文以連續(xù)性方程和動(dòng)量方程為控制方程，按照與深水試驗(yàn)池相同的搖板造波方式，采用VOF方法捕捉自由液面，建立了三維數(shù)值波浪水槽，為了防止水池末端的波浪發(fā)生反射而對工作區(qū)域造成影響，在動(dòng)量方程中添加源項(xiàng)進(jìn)行消波。采用六面體結(jié)構(gòu)化網(wǎng)格對數(shù)值模型進(jìn)行離散，保證了網(wǎng)格的質(zhì)量，提高了計(jì)算的精度。文中還對多種網(wǎng)格劃分形式進(jìn)行了詳細(xì)研究，同時(shí)與試驗(yàn)測量值進(jìn)行對比，確定了最優(yōu)的網(wǎng)格劃分方案，保證了數(shù)值方法的可行性。實(shí)際的模型試驗(yàn)中，因浪高儀測量設(shè)備數(shù)量有限，模型試驗(yàn)前確定平臺(tái)哪些位置安裝浪高儀對研究平臺(tái)周圍的波浪分布更具價(jià)值就顯得十分重要。此外，入射波陡較大時(shí)，波浪的非線性特征較強(qiáng)，與平臺(tái)之間的水動(dòng)力作用加劇，將導(dǎo)致平臺(tái)立柱周圍的波浪爬升效應(yīng)變得更為復(fù)雜，甚至產(chǎn)生砰擊、越浪等對平臺(tái)結(jié)構(gòu)和安全性構(gòu)成一定威脅的強(qiáng)非線性現(xiàn)象，相比平緩的波浪，研究波陡較大的波浪更有實(shí)際意義。介此，在該波浪水槽中，對陡波條件下半潛平臺(tái)周圍的波浪分布進(jìn)行數(shù)值模擬。計(jì)算時(shí)，在平臺(tái)周圍設(shè)定多個(gè)波浪監(jiān)測點(diǎn)，用以詳細(xì)預(yù)報(bào)不同位置處波浪的繞射和爬升效應(yīng)，確定波浪非線性較強(qiáng)的區(qū)域，為下一步模型試驗(yàn)中確定浪高儀的安裝位置提供重要依據(jù)。此外，文中還從物理現(xiàn)象上模擬并分析了波浪沿立柱的繞射和爬升過程，為今后進(jìn)一步研究波物之間的非線性相互作用，如波浪參數(shù)對立柱周圍波浪爬升的影響、極限波浪下半潛平臺(tái)的強(qiáng)非線性砰擊和越浪、半潛平臺(tái)的氣隙響應(yīng)等提供參考。

2 數(shù)值波浪

2.1 控制方程

半潛平臺(tái)波浪繞射和爬升的數(shù)值模擬是在一個(gè)類似于海洋深水試驗(yàn)池的數(shù)值水池中進(jìn)行的。數(shù)值模擬包括搖板造波的動(dòng)態(tài)區(qū)域，過渡區(qū)域，工作區(qū)域以及消波區(qū)域。

設(shè)定流體不可壓，則流場的連續(xù)性方程為：

相應(yīng)的動(dòng)量方程：

式中：x1，x2，x3分別為x，y，z方向的坐標(biāo)；u1，u2，u3分別為x，y，z方向的速度；f1=0，f2=0，f3=-g；S1，S2，S3分別為x，y，z方向上的附加源項(xiàng)；μ為動(dòng)力學(xué)粘性系數(shù)；ρ為流體的密度。

自由液面采用VOF（Volume of Fluid）方法進(jìn)行求解。該方法的基本原理是通過網(wǎng)格單元中的各流體和網(wǎng)格體積比函數(shù)aq來求解自由面的位置。aq=1為該網(wǎng)格單元全部由指定相流體所占據(jù)；aq=0說明該網(wǎng)格單元內(nèi)無指定相流體；0＜aq＜1說明該網(wǎng)格單元內(nèi)有部分指定相流體。氣液兩相流的自由面波動(dòng)的求解方程可寫為：

其中：q=1表示空氣相，q=2表示水相。

2.2 搖板造波

數(shù)值造波和消波的技術(shù)是數(shù)值波浪水池模擬成功的關(guān)鍵，一直是國內(nèi)外學(xué)者關(guān)注的問題。造波的方式有很多種，Boo等[11]通過在入口邊界給定速度或速度勢，然后基于時(shí)間步進(jìn)方式造出二階和三階stokes波。Brorsen和Larsen等[12]提出了一種適合積分方程法的域內(nèi)源造波方法，即在計(jì)算區(qū)域內(nèi)沿水深方向設(shè)定一造波源，在源項(xiàng)兩邊同時(shí)產(chǎn)生方向相反的兩列波。周勤俊和王本龍等[13]從忽略粘性的歐拉方程出發(fā)，以N-S方程為控制方程，將源項(xiàng)添加到動(dòng)量方程進(jìn)行造波和消波。董志等[14]采用仿物理的推板造波方式，通過設(shè)定造波板邊界按簡諧運(yùn)動(dòng)實(shí)現(xiàn)規(guī)則波的模擬，該方法造波原理簡單，物理造波機(jī)理論能較容易地應(yīng)用到數(shù)值造波中。

實(shí)際的深水試驗(yàn)池大多采用搖板式造波方法生成波浪，即通過機(jī)械驅(qū)動(dòng)使搖板繞固定軸擺動(dòng)，強(qiáng)迫水質(zhì)點(diǎn)運(yùn)動(dòng)，達(dá)到模擬相應(yīng)波動(dòng)的目的。波浪波高通過搖板的擺動(dòng)幅度控制，其擺動(dòng)頻率則決定波浪的周期大小。實(shí)驗(yàn)室對規(guī)則波的模擬中，搖板驅(qū)動(dòng)信號(hào)通過線性波浪理論計(jì)算得出。根據(jù)波浪理論，若搖板按固定角度做簡諧運(yùn)動(dòng)，在離開搖板一段距離后波浪將呈現(xiàn)二維規(guī)則波特性，即波浪周期與搖板簡諧運(yùn)動(dòng)相同，波高與擺幅有關(guān)。本文基于CFD理論，采用與上海交通大學(xué)深水試驗(yàn)池相同的搖板造波方式生成波浪。

對于規(guī)則波浪，造波板自由液面處的運(yùn)動(dòng)方程可寫為：

其中：S0為造波板自由液面處的運(yùn)動(dòng)幅值，k為波數(shù)。

通過勢流理論和造波板的運(yùn)動(dòng)方程，可推導(dǎo)出造波板運(yùn)動(dòng)幅值和波高H之間的轉(zhuǎn)換關(guān)系[15-16]，如（5）式。已知規(guī)則波浪的相關(guān)參數(shù)，為了確定造波板的運(yùn)動(dòng)規(guī)律，需要知道搖板造波方式下的水力傳遞函數(shù)T（ω,d），該值同樣可由勢流理論得出，如（6）式所示。

其中：d為水深，k為波數(shù)，ω為波浪圓頻率。而k可由下面的色散關(guān)系求得。

造波板的運(yùn)動(dòng)方程可寫為：

根據(jù)線性造波機(jī)理論，在水深為d的水池中，波面方程為：

本文利用動(dòng)網(wǎng)格模塊中的鋪層（Laying）功能中的網(wǎng)格更新技術(shù)，指定造波板為動(dòng)邊界，并繞某一固定點(diǎn)轉(zhuǎn)動(dòng)來實(shí)現(xiàn)數(shù)值造波。按照深水試驗(yàn)池中所采用的線性造波機(jī)原理，推導(dǎo)出搖板擺動(dòng)的角度信號(hào)，如（10）式所示。其中R為初始時(shí)刻造波板的水下長度。

2.3 源項(xiàng)消波

對于數(shù)值波浪水池，消波功能的實(shí)現(xiàn)有很多。如主動(dòng)式消波、阻尼層消波、Sommerfeld輻射邊界消波、動(dòng)量源項(xiàng)消波和多孔介質(zhì)消波等。主動(dòng)式消波[17]是指在造波的同時(shí)，通過實(shí)時(shí)的修正造波板的運(yùn)動(dòng)，在生成目標(biāo)波浪的同時(shí)產(chǎn)生額外的波動(dòng)抵消反射波浪。而被動(dòng)式消波方法如動(dòng)量源項(xiàng)消波，通過在動(dòng)量方程中添加源項(xiàng)以消除反射波浪。海綿層阻尼消波[17]主要是通過在特定的計(jì)算區(qū)域設(shè)定1-2倍波長的人工阻尼區(qū)，對波浪垂向速度進(jìn)行強(qiáng)迫衰減，同時(shí)開邊界處滿足Sommerfeld輻射條件。多孔介質(zhì)消波方法[19]是一種仿物理消波法，通過在N-S方程中增加動(dòng)量衰減的源項(xiàng)達(dá)到消波的目的。

模擬計(jì)算過程中，采用王本龍等[13]提出的消波技術(shù)，在消波區(qū)域內(nèi)添加源項(xiàng)來消除水池尾部池壁的波浪反射，確保工作區(qū)域波浪不受干擾。消波段內(nèi)波動(dòng)場可寫為：

i=1方向，已加源項(xiàng)的動(dòng)量方程按照忽略粘性的歐拉方程出發(fā)可離散為：

未添加源項(xiàng)的動(dòng)量方程可離散成：

將（11）式代入方程聯(lián)立求解可得源項(xiàng)S1的表達(dá)式：

同理i=3方向的源項(xiàng)S3的表達(dá)式：

其中下標(biāo)C代表計(jì)算值；上標(biāo)N+1、N分別代表N+1、N時(shí)刻的值；λ為與空間位置有關(guān)的光滑過渡加權(quán)函數(shù)，避免速度的劇烈變化影響流場的求解精度。λ=λ（x）在數(shù)值水池消波段的表達(dá)式為：

式中：xb為消波區(qū)的起始位置，xe為消波區(qū)的末端，因此可知：

2.4 幾何模型的建立

本文采用ICEMCFD軟件建立了波浪模擬所需的幾何模型，其示意圖如圖1所示，相應(yīng)的三維立體模型見圖8和圖9。設(shè)定坐標(biāo)原點(diǎn)位于平臺(tái)重心在初始自由液面的投影處。x軸正向與波浪傳播方向一致，y軸正方向與平臺(tái)橫向方向一致，z軸沿平臺(tái)的吃水方向。該數(shù)值波浪水池長24 m，寬5 m，水深8 m，靜水面以上的高度為1 m，造波板位于X軸-10 m位置處，各區(qū)域長分別為搖板及前端過渡區(qū)5 m，工作區(qū)10 m，尾部過渡區(qū)5 m，尾部消波區(qū)4 m，在水池試驗(yàn)區(qū)域的中央位置設(shè)置一虛擬浪高儀來獲取波浪的時(shí)歷。

圖1 數(shù)值波浪水池的示意圖Fig.1 Schematic show of numerical wave tank

圖2 電容式浪高儀測量系統(tǒng)Fig.2 Measurement system of capacitive wave probes

2.5 波浪時(shí)歷的試驗(yàn)測量

為了驗(yàn)證數(shù)值波浪水池的精確度，在上海交通大學(xué)海洋深水池內(nèi)對波浪時(shí)歷進(jìn)行了實(shí)際測量。用于該試驗(yàn)研究的測量儀器主要為電容式精密浪高儀，該儀器由探頭、弓形支架、電纜以及信號(hào)輸出系統(tǒng)等組成，如圖2所示，其通常安裝在水池拖車下方，通過測量因水位引起的電容變化來獲取波浪時(shí)歷。試驗(yàn)前，首先需對浪高儀進(jìn)行標(biāo)定，以確定電壓轉(zhuǎn)換系數(shù)。待水體充分平靜后，通過數(shù)據(jù)采集系統(tǒng)對浪高儀進(jìn)行采零處理，正式采集得到的波浪時(shí)歷均減去該采零值。隨后，在控制程序中輸入波浪相關(guān)參數(shù)并轉(zhuǎn)化為搖板信號(hào)來驅(qū)動(dòng)造波板實(shí)現(xiàn)物理造波。待波形穩(wěn)定之后進(jìn)行數(shù)據(jù)采集，該次試驗(yàn)的采樣頻率選定為40 Hz，采樣時(shí)間3分鐘，采樣點(diǎn)數(shù)大約為7 200點(diǎn)。

2.6 網(wǎng)格尺度的測試

CFD計(jì)算中，網(wǎng)格的質(zhì)量、尺度和布置方式等將直接影響數(shù)值計(jì)算結(jié)果的精度。網(wǎng)格質(zhì)量較差勢必影響數(shù)值結(jié)果的收斂性，甚至可能造成解的發(fā)散。六面體結(jié)構(gòu)網(wǎng)格具有存儲(chǔ)空間小、數(shù)據(jù)結(jié)構(gòu)簡單、迭代運(yùn)算速度快、網(wǎng)格質(zhì)量高的優(yōu)點(diǎn)，且對于VOF模型來說，排列整齊，形狀規(guī)則的六面體網(wǎng)格不僅可以提高求解精度，同時(shí)可以捕捉到更為光順的自由液面。因此，本文全部采用H-H型六面體結(jié)構(gòu)化網(wǎng)格對數(shù)值波浪水槽進(jìn)行離散，以確保網(wǎng)格質(zhì)量良好且方便布置網(wǎng)格節(jié)點(diǎn)。此外，過密的網(wǎng)格將大大增加迭代計(jì)算所需的時(shí)間，過于稀疏的網(wǎng)格不僅無法捕捉到流場的相關(guān)細(xì)節(jié)，而且將引起嚴(yán)重的數(shù)值耗散，以致計(jì)算結(jié)果失真。一般而言，在進(jìn)行網(wǎng)格的生成時(shí)需遵循的原則是：在不影響計(jì)算精度的情況下，盡量減少網(wǎng)格數(shù)量，提高計(jì)算的速度。因此本文在網(wǎng)格劃分時(shí)對一些流場梯度變化較大的區(qū)域進(jìn)行了加密，而對物理量不敏感的區(qū)域采用了粗網(wǎng)格，如自由液面附近采用較細(xì)的網(wǎng)格，遠(yuǎn)離這些區(qū)域網(wǎng)格間距按指數(shù)分布逐漸拉大，以起到疏密有致的效果。但即便如此，計(jì)算精度也不能完全得到保證，因此必須同試驗(yàn)測量值進(jìn)行對比，才能確定最終的網(wǎng)格劃分方式。

結(jié)合試驗(yàn)室的造波能力，本文選取了陡度較大的規(guī)則波浪作為研究的對象，即波陡（H/λ，其中λ為入射波波長）為1/15，周期T為1.0 s，波高H為0.104 m，采用三種密度由疏到密的網(wǎng)格劃分方式對該波浪進(jìn)行了數(shù)值模擬?？紤]到AB段網(wǎng)格對數(shù)值造波是否能夠成功有重要的影響，為了減少波浪傳播過程中的數(shù)值耗散，對其等距離布置較細(xì)密的網(wǎng)格，并按一個(gè)波長內(nèi)100個(gè)、50個(gè)和25個(gè)節(jié)點(diǎn)三種方式來檢驗(yàn)網(wǎng)格的不同間距對計(jì)算結(jié)果帶來的影響。此外，考慮到自由液面附近，相應(yīng)的物理量梯度變化大，在EF段同樣采用等距離的網(wǎng)格布置方式。且三種方案中，單個(gè)波高內(nèi)的網(wǎng)格數(shù)分別為40、20和10個(gè)。BC段為模型區(qū)域，該段的網(wǎng)格參數(shù)無論對計(jì)算迭代的收斂時(shí)間和精度，還是對模型周圍的波動(dòng)場細(xì)節(jié)的模擬等都起著至關(guān)重要的作用，因此該處的網(wǎng)格應(yīng)更為細(xì)密。其余段的網(wǎng)格劃分均按指數(shù)分布的形式，尺度由梯度變化大的區(qū)域往梯度小的區(qū)域逐漸拉大，具體如表1所示。

表1 網(wǎng)格尺度的測試Tab.1 Grid size test

從圖4可以看出，網(wǎng)格C的計(jì)算結(jié)果與試驗(yàn)值偏差嚴(yán)重，且隨著模擬時(shí)間的增加，波形因數(shù)值的耗散逐漸衰減，因此不符合高品質(zhì)波浪的要求。而網(wǎng)格A和網(wǎng)格B所得到的結(jié)果均與試驗(yàn)值較為接近，說明網(wǎng)格B的劃分方案可以使模擬結(jié)果達(dá)到要求。相比網(wǎng)格B，劃分最為精細(xì)的網(wǎng)格A并沒有帶來更好的精度，且由于網(wǎng)格數(shù)量巨大，將大大地消耗數(shù)值運(yùn)算的時(shí)間和資源。因此，從計(jì)算的精確度和時(shí)間消耗兩方面權(quán)衡考慮，最終我們選取網(wǎng)格B作為數(shù)值波浪水池網(wǎng)格劃分的基礎(chǔ)。

圖5為網(wǎng)格B方案下，數(shù)值波浪水池中央位置處未經(jīng)擾動(dòng)的波形時(shí)歷曲線，造波板周圍的波形分布如圖7。通過與線性Airy波理論值對比發(fā)現(xiàn)，計(jì)算得到的波峰高而尖，波谷較為平坦，波浪出現(xiàn)了非線性的特征，已經(jīng)不再滿足線性波理論。

圖3 數(shù)值波浪水池網(wǎng)格劃分方案的示意圖Fig.3 Schematic show of mesh strategy of numerical wave tank

圖4 不同網(wǎng)格劃分方案的波浪時(shí)歷分布Fig.4 Time series of wave elevation model in different mesh scheme

圖5 試驗(yàn)值與線性波理論值對比Fig.5 Comparison of wave elevation between test and linear wave theory

圖6 計(jì)算值與線性波理論值的波浪譜密度對比Fig.6 Comparison of wave spectrum between calculation and linear theory

圖7 網(wǎng)格B方案下?lián)u板附近的瞬時(shí)波形圖Fig.7 Instantaneous wave shape around flap wavemaker in MESH-B

采用FFT理論可將一段時(shí)歷的波列轉(zhuǎn)化為基于波頻的譜密度分布形式，從而得到波浪的各階諧頻成份以及相應(yīng)的貢獻(xiàn)[22]。從圖6中可以看出，airy波為線性規(guī)則波，其僅包含一階諧頻成份，由于入射波周期T為1.0 s，故譜峰值出現(xiàn)在頻率ω約為6.28 rad/s處。此外，我們對網(wǎng)格B劃分方案計(jì)算得到的時(shí)歷曲線進(jìn)行頻譜分析發(fā)現(xiàn)，譜密度曲線中出現(xiàn)兩個(gè)峰值：第一個(gè)峰值出現(xiàn)在頻率ω約6.28 rad/s位置，即主頻位置，表示一階諧頻成份（First-harmonic component）；第二個(gè)譜峰對應(yīng)的頻率約為12.56 rad/s，其值的大小恰好為主頻的2倍，對應(yīng)的譜峰值稱為二階諧頻成份（Second-harmonic component）。故認(rèn)為該入射波具有較明顯的二階非線性特征。盡管波浪的一階線性成份所占的比例較大，其高階諧頻成份卻是影響波浪非線性爬升效應(yīng)的重要因素，在以后的研究中將對此重點(diǎn)考察。

3 波浪繞射和爬升的數(shù)值計(jì)算

3.1 模型建立

本文在上述理論建立的數(shù)值波浪水池內(nèi)，以中海油3 000 m深水半潛式鉆井平臺(tái)為研究對象，對其周圍的波浪繞射和波浪爬升現(xiàn)象進(jìn)行了數(shù)值模擬。文章重點(diǎn)研究波浪散射效應(yīng)（包括繞射、反射等效應(yīng)）對平臺(tái)周圍波浪爬升的影響，因此設(shè)定平臺(tái)為剛性固定，以消除運(yùn)動(dòng)引起的輻射效應(yīng)的干擾。

表2 中海油半潛平臺(tái)主要參數(shù)Tab.2 Main pariticulars of CNOOC’semisubmersible platform

首先建立了計(jì)算所需的幾何模型，如圖8和9所示。該半潛平臺(tái)包含水下浮箱、四根直立柱、橫撐和平臺(tái)甲板等，模型縮尺比為1:60，平臺(tái)實(shí)型和模型的具體尺度見表2。半潛平臺(tái)置于數(shù)值波浪水池工作區(qū)域的中央位置，平臺(tái)吃水為34.5 cm，初始?xì)庀叮o水面至平臺(tái)下甲板的垂直距離）為15.6 cm。

為了測量平臺(tái)周圍波面分布情況，在平臺(tái)前立柱周圍設(shè)置5排浪高儀，各排浪高儀間的角度為45°，后立柱周圍設(shè)有1排浪高儀。為了揭示波浪沿立柱爬升的過程，每排均設(shè)置四根浪高儀，其總長（R）為立柱的等效半徑?？紤]到波浪越靠近立柱，波浪爬升效應(yīng)越復(fù)雜，故設(shè)置浪高儀沿徑向的布置間距（r）由外到內(nèi)逐漸變小，分別為0.4r/R，0.4r/R，0.193 75r/R和0.006 25r/R，且單排內(nèi)浪高儀的編號(hào)依次為iA，iB，iC和iD（i為浪高儀的排號(hào)），如圖10所示。

3.2 邊界條件

圖8 半潛平臺(tái)三維數(shù)值模型Fig.8 3D numerical model of semi-submersible

圖9 計(jì)算區(qū)域的幾何模型Fig.9 Geometrical model of calculated domains

如圖9所示，水池頂部設(shè)置為壓力入口邊界條件，半潛平臺(tái)和造波板為無滑移壁面。因數(shù)值計(jì)算中，所選的波浪為迎浪狀態(tài)（與平臺(tái)艏部方向夾角為180°），且平臺(tái)沿中縱剖面對稱，因此水池內(nèi)流場亦關(guān)于中縱剖面對稱，為了節(jié)省網(wǎng)格的數(shù)量，設(shè)定水池中縱剖面為對稱面，其余均默認(rèn)為壁面。波動(dòng)場采用層流模式（Laminar）求解，通過VOF方法追蹤自由面；壓力—速度耦合項(xiàng)采用SIMPLEC算法進(jìn)行迭代計(jì)算，動(dòng)量方程中的對流項(xiàng)和擴(kuò)散項(xiàng)均采用二階迎風(fēng)格式；VOF方程采用幾何重構(gòu)法[20]（Geo-Reconstruct）求解。在數(shù)值模擬時(shí)，設(shè)定左邊界上部的搖板為運(yùn)動(dòng)邊界，使搖板繞其末端轉(zhuǎn)動(dòng)，采用鋪層技術(shù)控制動(dòng)態(tài)造波區(qū)域的網(wǎng)格變形、潰滅和再生來產(chǎn)生波浪。

圖10 虛擬浪高儀的布置位置Fig.10 Layout of virtual wave probes

3.3 網(wǎng)格劃分

圖11 計(jì)算區(qū)域網(wǎng)格分布Fig.11 Meshes of calculated domains

圖12 造波板和自由液面附近網(wǎng)格劃分Fig.12 Meshes around wavemaker and free surface

圖13 平臺(tái)周圍的網(wǎng)格劃分Fig.13 Meshes around semi-submersible

圖14 立柱周圍的網(wǎng)格劃分Fig.14 Meshes around columns

CFD計(jì)算中，網(wǎng)格離散化的質(zhì)量將對數(shù)值計(jì)算結(jié)果的精度產(chǎn)生重要影響，本文采用分塊耦合網(wǎng)格劃分技術(shù)，將整個(gè)流場劃分成若干區(qū)域，各流域均采用六面體結(jié)構(gòu)化網(wǎng)格進(jìn)行離散。多塊結(jié)構(gòu)化網(wǎng)格是單塊結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格的一個(gè)折中方案，它簡單而又能處理較復(fù)雜的形狀。按照上面闡述的方案生成流場網(wǎng)格，部分區(qū)域相應(yīng)的網(wǎng)格分布如圖11-14所示，整個(gè)模型網(wǎng)格單元在360萬左右。

3.4 計(jì)算結(jié)果與分析

圖15為Mavrakos等[21]的模型試驗(yàn)中，規(guī)則波浪條件下平臺(tái)立柱周圍的波浪爬升情形。本文也得到了類似的波面分布，如圖16所示，且波浪沿立柱往上爬升過程中形成的“水舌”也極為相似，基于VOF方法求解得到的自由液面分布也較光順，因此文中采用的數(shù)值計(jì)算方法能夠較好地模擬這一非線性特征。

圖15 模型試驗(yàn)中波浪爬升情形Fig.15 Wave run-up in model tests

圖16 數(shù)值模擬中立柱周圍的波浪爬升Fig.16 Wave run-up in numerical simulation

圖17為前端立柱、后端立柱迎浪面中央位置，波浪沿柱體爬升時(shí)最大波面分布情況。從圖中可以看出，由于正面迎浪，兩處的波浪爬升都較大，且隨著波浪向前傳播，最大波面升高逐漸增大，在靠近立柱壁面時(shí)波浪迅速爬升至最高位置，這與Morris-Thomas[22]以及Nielsen等[23]研究結(jié)論相一致。

圖18給出的是前后立柱壁面附近的波浪爬升時(shí)歷分布。從圖中可以看出，相比平臺(tái)前立柱壁面附近的3D位置，6D處（后立柱的壁面附近）波浪時(shí)歷曲線顯示出更強(qiáng)的非線性特征，且波浪爬升幅值更大，部分波浪已經(jīng)爬升至下層甲板，此時(shí)上端水體仍具有一定能量，將對甲板局部產(chǎn)生載荷沖擊。同時(shí)結(jié)合圖24（b）可知，文中采用的計(jì)算方法可以捕捉到波浪砰擊后立柱表面時(shí)所產(chǎn)生的水體分離和破碎等強(qiáng)非線性現(xiàn)象。由于該類砰擊載荷可能造成平臺(tái)甲板的結(jié)構(gòu)破壞，因此已成為國際上研究的熱點(diǎn)問題之一，今后將對該強(qiáng)非線性現(xiàn)象展開詳細(xì)分析。

圖17 最大波面升高隨r/R的變化曲線Fig.17 Maximum wave elevation along wave probe row

圖18 立柱壁面附近的波浪爬升時(shí)歷Fig.18 Time history of wave elevation near column wall

圖19顯示的是虛擬浪高儀組2（虛擬浪高儀組3順時(shí)針方向45°位置）和4處（虛擬浪高儀組3逆時(shí)針方向45°位置）的最大波面分布情況。盡管波浪爬升幅度比柱體的迎浪面中央位置小，但仍較為明顯，并且表現(xiàn)出較強(qiáng)的非線性特性。虛擬浪高儀組2處最大的波面升高規(guī)律與虛擬浪高儀組3與6處基本一致，越靠近立柱，波浪爬坡現(xiàn)象越明顯。然而，虛擬浪高儀組4處，最大的波浪爬升并未發(fā)生在立柱壁面附近，而是離立柱較遠(yuǎn)的4B位置，造成這種現(xiàn)象的原因可能是由于波浪遇前立柱后發(fā)生繞射和反射，部分反射的水體與入射波在遠(yuǎn)離柱體位置發(fā)生疊加所致，該現(xiàn)象同時(shí)由Mavrakos等[21]在模型試驗(yàn)中得以證實(shí)。

此外，通過波浪時(shí)歷對比發(fā)現(xiàn)，盡管4B與2D兩位置處最大波面升高相接近，但在4B處，由于立柱內(nèi)側(cè)區(qū)域受到的干擾因素多，使得波谷位置單個(gè)周期內(nèi)出現(xiàn)了明顯的二次波峰，具有比2D處更強(qiáng)的高階諧頻非線性特征。

當(dāng)波浪繞過虛擬浪高儀組2和4處繼續(xù)往后傳播時(shí)，其爬升幅值迅速減小。且一個(gè)立柱半徑R范圍內(nèi)，最大波面升高基本無變化，但其非線性特征卻越明顯,如圖21、22所示。

圖19 最大波面升高隨r/R的變化曲線Fig.19 Maximum wave elevation along wave probe row

圖20 立柱壁面附近的波浪爬升時(shí)歷Fig.20 Time history of wave elevation near column wall

圖21 最大波面升高隨r/R的變化曲線Fig.21 Maximum wave elevation along wave probe row

圖22 立柱壁面附近的波浪爬升時(shí)歷Fig.22 Time history of wave elevation near column wall

對比沿平臺(tái)內(nèi)側(cè)布置的浪高儀組3、4、5發(fā)現(xiàn)，3處的波浪爬升最為顯著，隨著波浪繞立柱繼續(xù)向前傳播，波高逐漸減小，波浪爬升效應(yīng)明顯減弱，這與Nielsen等[23]在對半潛平臺(tái)的單立柱展開模型試驗(yàn)時(shí)所得出的結(jié)論相同。值得注意的是，3處波浪爬升幅度雖相對較大，但波形較為規(guī)整，波谷的形狀較為對稱，隨著波浪繞立柱往前傳播至4處，波谷位置出現(xiàn)了二次波峰，在波浪傳至5處時(shí)，二次波峰幅值不斷的增大，同時(shí)主波峰（一個(gè)波浪周期內(nèi)第一次出現(xiàn)的波峰）的形狀更尖瘦，該過程將使波浪非線性效應(yīng)增強(qiáng)，與平臺(tái)間的相互作用變得更為復(fù)雜。

為了更為清晰地理解波物之間的相互作用，文中對沿平臺(tái)立柱的繞射和爬升現(xiàn)象進(jìn)行了模擬。圖23為平臺(tái)前立柱周圍波浪的爬升過程。該情景可描述為：波浪遇立柱時(shí)，部分水流繞過立柱繼續(xù)朝前傳播，而另一部分水體因立柱的阻礙開始沿立柱爬升，此時(shí)波浪的能量很大。隨著波浪繼續(xù)往上爬升，其能量不斷耗散，興起的波峰不斷增高，波峰周圍的速度逐漸減小。待波峰增至一定位置后停止升高，在重力的作用下開始回落，但此時(shí)緊貼立柱表面的水體卻繼續(xù)往上爬升，如圖（c）所示，直至其速度接近0時(shí)，波浪開始下滑，此時(shí)可以觀察到，回落的水體沿立柱繞射至側(cè)面時(shí)將興起波浪，且其波峰周圍的速度顯著增大，如圖（d）所示。

圖24顯示的是后立柱迎浪面波浪的爬升現(xiàn)象。從圖（a）可以看出，后立柱周圍的波浪已上升到平臺(tái)下甲板并與之發(fā)生非線性砰擊。砰擊結(jié)束后，下甲板上的部分水流由于重力作用直接脫落，剩余部分沿立柱回落，同時(shí)有小部分水體與尾部水舌發(fā)生分離，如圖（b），這種現(xiàn)象可能是因水流下滑過程中，兩部分水體速度不同引起的。相比前立柱，該規(guī)則波浪下，后立柱周圍的波浪爬升情況更為嚴(yán)重，非線性效應(yīng)更強(qiáng)，原因可能因后立柱位于前立柱的尾流中，其周圍的波面更易受擾動(dòng)影響。

圖23 前立柱迎浪面波浪的爬升過程Fig.23 Wave run-up process on front side of forward columns

圖24 后立柱迎浪面波浪爬升現(xiàn)象Fig.24 Wave run-up process on front side of aft columns

4 結(jié) 論

文中基于N-S方程和VOF方法建立了三維數(shù)值波浪水槽，運(yùn)用類似于深水試驗(yàn)池的搖板造波方式產(chǎn)生波浪。采用H-H型六面體結(jié)構(gòu)化網(wǎng)格對含平臺(tái)模型的數(shù)值波浪水槽進(jìn)行離散，保證了網(wǎng)格的質(zhì)量，提高了計(jì)算的精度。此外，文中還對多種網(wǎng)格劃分形式進(jìn)行了詳細(xì)研究，同時(shí)與試驗(yàn)測量值進(jìn)行對比，確定了最優(yōu)的網(wǎng)格劃分方案，并保證了數(shù)值方法的可行性。隨后，在該水池中對陡波條件下半潛平臺(tái)周圍的波浪繞射和爬升效應(yīng)進(jìn)行了數(shù)值模擬。結(jié)果表明，本文采用的計(jì)算方法能夠?qū)ζ脚_(tái)周圍的波面分布給出詳細(xì)的數(shù)值模擬，同時(shí)可預(yù)報(bào)出波浪沿平臺(tái)立柱的繞射以及從波浪開始爬升至回落的完整過程，為下一步有關(guān)平臺(tái)氣隙的模型試驗(yàn)前確定浪高儀的安裝位置提供了重要參考。計(jì)算獲得的波浪爬升情形和模型試驗(yàn)相似，且文中得到的重要結(jié)論與相關(guān)文獻(xiàn)一致，進(jìn)一步證明該方法模擬波浪繞射和爬升等現(xiàn)象是可行的。通過計(jì)算和分析可以得到以下結(jié)論：

（1）立柱正面迎浪位置，波浪爬升較大，且越靠近立柱，波峰越尖瘦，波浪爬升越為嚴(yán)重。因此在下一步的模型試驗(yàn)中，需在該位置附近安裝浪高儀進(jìn)行進(jìn)一步研究。

（2）一般而言，越靠近立柱波浪爬升效應(yīng)越明顯，但也有例外，如偏離前立柱中縱剖面45°的內(nèi)側(cè)位置附近（浪高儀組4附近），可能的原因是由于內(nèi)側(cè)位置發(fā)生繞射和反射，部分反射的水體與入射波在遠(yuǎn)離柱體位置發(fā)生疊加所致。

（3）平臺(tái)立柱的內(nèi)側(cè)，波面受干擾的因素多，波浪爬升現(xiàn)象比相應(yīng)的外側(cè)更為明顯，爬升幅度更大，成因復(fù)雜，應(yīng)加以詳細(xì)研究。

（4）由于波浪沿前立柱的繞射、遇后立柱的反射以及波浪疊加等效應(yīng)的相互作用，平臺(tái)前立柱的背面附近以及前后立柱之間區(qū)域的波面分布將變得較為復(fù)雜，波浪的非線性效應(yīng)增強(qiáng)，應(yīng)給予足夠重視。

[1]王世圣,謝 彬,曾恒一,馮 瑋,李曉平,張海濱.3 000米深水半潛式鉆井平臺(tái)運(yùn)動(dòng)性能研究[J].中國海上油氣,2007,19(4):277-280.

[2]Stansberg C T,Baarholm R,Kristiansen T.Extreme wave amplification and impact loads on offshore structures[C].Offshore Technology Conference,OTC 17487,2004.

[3]McCamy R,Fuchs R.Wave forces on piles:A diffraction theory[R].Tech.Memo No.69,Beach Erosion Board,US Army Corps of Engineers,Washington DC,USA,1954.

[4]Niedzwecki J N,Duggal A S.Wave run-up and forces on cylinders in regular and random waves[J].J Waterways,Port,Coastal,Ocean Eng,1992,118:615-634.

[5]Niedzwecki J M,Huston J R.Wave interaction with tension leg platforms[J].Ocean Engineering,1992,19(1):21-37.

[6]Lee C H,Newman J N.Second-order wave effects on offshore structures[C].Proceedings of the Seventh International Conference on Behaviour of Offshore structures,1994,2.

[7]Donald G Danmeier,Robert K M Seah,Timothy Finnigan,Dominique Roddier.Validation of wave run-up calculation methods for a gravity based structure[C].Proceedings of the 27th OMAE,2008.

[8]Teigen P,Trulsen K.Numerical investigation of non-linear wave effects around multiple cylinders[C].11th International Offshore and Polar Engineering Conference,2001.

[9]Danmeier G D,Robert K M S,Finnigan T,et al.Timothy Finnigan,Dominique Roddier.Validation of wave run-up calculation methods for a gravity based structure[C].Proceedings of the 27th OMAE,2008.

[10]Roberto M,Andrea D M,Riccardo B.Numerical simulation of the flow around an array of free-surface piercing cylinders in waves[C].Proceedings of the 26th OMAE,2007.

[11]Boo S Y,Kim C H.Simultation of fully nonlinear irregular waves in a 3-D numerical wave tank[C].Proceedings of the 17th ISOPE,1994.

[12]Brorsen M,Larsen.Source generation of nonlinear gravity waves with the boundary integral equation method[J].Coastal Engineering,1987,11(2):93-113.

[13]周勤俊,王本龍,蘭雅梅,劉 樺.海堤越浪的數(shù)值模擬[J].力學(xué)季刊,2005,26(4):629-633.

[14]董 志,詹杰民.基于VOF方法的數(shù)值波浪水槽以及造波、消波方法研究[J].水動(dòng)力研究與進(jìn)展A輯,2009,24(1):15-21.

[15]孫昭晨.推搖混合式造波機(jī)理論曲線[J].港口工程,1988(4):30-32.

[16]陶建華.水波的數(shù)值模擬[M].天津:天津大學(xué)出版社,2004.

[17]李雪臨,任 冰,王永學(xué).VOF方法中主動(dòng)吸收式無反射數(shù)值造波研究[C].第二十屆全國水動(dòng)力學(xué)研討會(huì)文集,2006.

[18]韓 朋,任 冰,李雪臨,王永學(xué).基于VOF方法的不規(guī)則波數(shù)值波浪水槽的阻尼消波研究[J].水道港口,2009,30(1).

[19]詹杰民,董 志.黏性數(shù)值波浪水槽的多孔介質(zhì)消波方法[C].第十四屆中國海洋（岸）工程學(xué)術(shù)討論會(huì)論文集（上冊）,2009.

[20]Youngs D L.Time-dependent multi-material flow with large fluid distortion[C]//Proceedings of Numerical Methods for Fluid Dynamics.New York,1982.

[21]Mavrakos S A,Chatjigeorgiou I K,Grigoropoulos G,Maron A.Scale experiment of motions and wave run-up on a TLP model,subjected to monochromatic waves[C].Proceedings of the 23rd OMAE,2004.

[22]Morris-Thomas M T,Thiagarajan K P.The run-up on a cylinder in progressive surface gravity waves:Harmonic components[J].Applied Ocean Research,2004,26(3-4):98-113.

[23]Nielsen F G.Comparative study on airgap under floating platforms and run-up along platform columns[J].Marine Structures,2003,16(2):97-134.