基于CFD流固耦合理論的海上浮式結(jié)構(gòu)物水動(dòng)力性能分析

馬哲，程勇，翟鋼軍

（1.大連理工大學(xué)深海工程研究中心，遼寧大連116024；2.江蘇科技大學(xué)船舶與海洋工程學(xué)院，江蘇鎮(zhèn)江212003）

基于CFD流固耦合理論的海上浮式結(jié)構(gòu)物水動(dòng)力性能分析

馬哲1，程勇2，翟鋼軍1

（1.大連理工大學(xué)深海工程研究中心，遼寧大連116024；2.江蘇科技大學(xué)船舶與海洋工程學(xué)院，江蘇鎮(zhèn)江212003）

文章利用流體動(dòng)力學(xué)控制方程和結(jié)構(gòu)運(yùn)動(dòng)方程的耦合理論，在具有造消波功能的2D數(shù)值水槽中實(shí)現(xiàn)了海上浮式結(jié)構(gòu)物在波浪中運(yùn)動(dòng)過程的數(shù)值模擬。以系泊式浮式方箱防波堤作為工程應(yīng)用實(shí)例，分別采用梯形法和二、三階單步數(shù)值迭代法對(duì)浮體的動(dòng)力特性進(jìn)行數(shù)值分析，并與基于勢(shì)流理論的邊界元法的結(jié)果做對(duì)比，分析后發(fā)現(xiàn)在波浪條件下梯形法和二、三階單步法的計(jì)算精度相當(dāng)，結(jié)果收斂，且與邊界元法結(jié)果吻合較好，滿足要求。該文提出了一種新的全自由度分塊結(jié)構(gòu)移動(dòng)網(wǎng)格技術(shù)，實(shí)現(xiàn)了浮體所有方向的聯(lián)合運(yùn)動(dòng)，而且網(wǎng)格不發(fā)生任何扭曲現(xiàn)象，計(jì)算時(shí)間、網(wǎng)格劃分和精度要求都得到了較好的控制。

流固“全耦合”；梯形法；二、三階單步法；全自由度分塊結(jié)構(gòu)移動(dòng)網(wǎng)格技術(shù)

0 引言

在海洋工程領(lǐng)域，海洋平臺(tái)、船舶和FPSO等海洋結(jié)構(gòu)物在波浪作用下經(jīng)常會(huì)發(fā)生波浪運(yùn)動(dòng)、砰擊和甲板上浪等現(xiàn)象。為了保證海洋結(jié)構(gòu)物在自存和工作工況下的安全性，因此對(duì)各種水動(dòng)力現(xiàn)象進(jìn)行特性分析起著至關(guān)重要的作用。在計(jì)算流體動(dòng)力學(xué)中，研究此問題的方法大致可分為兩類：一類是根據(jù)勢(shì)流理論，首先預(yù)報(bào)出結(jié)構(gòu)在波浪作用下的運(yùn)動(dòng)規(guī)律，然后規(guī)定結(jié)構(gòu)按照此規(guī)律運(yùn)動(dòng)，從而實(shí)現(xiàn)結(jié)構(gòu)的運(yùn)動(dòng)和流場(chǎng)變化，其中對(duì)于結(jié)構(gòu)運(yùn)動(dòng)采用移動(dòng)網(wǎng)格技術(shù)實(shí)現(xiàn)。林兆偉，朱仁傳等學(xué)者[1]在一個(gè)具有造消波功能的二維數(shù)值水槽中提出多層動(dòng)靜網(wǎng)格匹配方法，研究了甲板上浪問題；郭海強(qiáng)[2]采用動(dòng)計(jì)算域和滑移網(wǎng)格匹配法計(jì)算了船體剖面運(yùn)動(dòng)的水動(dòng)力系數(shù)；龔丞[3]等結(jié)合勢(shì)流理論與CFD的優(yōu)勢(shì)，模擬分析了甲板上浪現(xiàn)象及其對(duì)甲板結(jié)構(gòu)的沖擊作用，船體運(yùn)動(dòng)規(guī)律由勢(shì)流理論確定，然后通過FLUENT軟件進(jìn)行船舶的甲板上浪模擬。另一類是根據(jù)流體和浮體的耦合作用，首先算出流體對(duì)浮體的作用力，然后根據(jù)浮體的運(yùn)動(dòng)方程算出浮體的速度和位置，浮體的速度和位置可以改變周圍的流場(chǎng)，即浮體和流體之間是相互耦合的。Hadzic,Henning[4]等學(xué)者基于Reynolds平均化的Navier-Stokes（N-S）方程和剛體運(yùn)動(dòng)方程的耦合迭代求解法，采用網(wǎng)格的變形和滑移技術(shù)預(yù)測(cè)浮體在波流作用下的運(yùn)動(dòng)規(guī)律；Nielsen[5]等基于Navier-Stokes方程，采用VOF自由液面捕捉法研究了FPSO在波浪作用下的甲板上浪問題；Yang和Qiu[6]等基于CIP（Constrained Interpolation Profile）的有限差分法求解N-S方程計(jì)算二維及三維物體的波浪抨擊力，其中流體壓力采用共軛梯度迭代法（conjugate gradient iterative method）求解泊松型方程得到；Kleefsman[7]等基于不可壓縮粘性流體的N-S方程研究了由波浪抨擊力引起的潰壩及物體入水問題，詳細(xì)闡述了在固定笛卡爾坐標(biāo)系下基于有限體積法的N-S方程離散方法，并對(duì)潰壩及楔體入水實(shí)驗(yàn)進(jìn)行了數(shù)值仿真模擬；Shen和Yang[8]等通過實(shí)驗(yàn)及FLUENT軟件研究了聚焦波浪作用下風(fēng)機(jī)基礎(chǔ)的波浪荷載問題。

對(duì)于預(yù)先知道浮體運(yùn)動(dòng)形式的水動(dòng)力求解問題，由于運(yùn)動(dòng)規(guī)律根據(jù)勢(shì)流理論分析得出，流體作用在濕表面上的壓力和剪應(yīng)力不會(huì)改變浮體的速度和位置，因此這種方法不是真正意義上的流體和浮體的耦合模擬，而是一種“半耦合”數(shù)值模擬技術(shù)，此方法明顯的缺點(diǎn)是不能模擬出浮體的漂移效果。在數(shù)值造消波的水槽中研究浮體的運(yùn)動(dòng)和受力，最理想的途徑是“全耦合”數(shù)值方法，但由于此方法計(jì)算量大，對(duì)動(dòng)網(wǎng)格的質(zhì)量和數(shù)值計(jì)算精度要求高，因此到目前為止研究甚少。本文采用流體動(dòng)力學(xué)控制方程和浮體運(yùn)動(dòng)方程的耦合技術(shù)，在2D數(shù)值造消波水槽中實(shí)現(xiàn)了對(duì)浮體運(yùn)動(dòng)和受力等水動(dòng)力性能的研究。其中對(duì)浮體的運(yùn)動(dòng)規(guī)律求解分別采用了梯形法和二、三階單步法，并使用VOF法追蹤自由液面來模擬浮體對(duì)流場(chǎng)的耦合作用，浮體所有方向的聯(lián)合運(yùn)動(dòng)通過一種新的全自由度分塊動(dòng)靜結(jié)構(gòu)網(wǎng)格和平動(dòng)、旋轉(zhuǎn)滑移交界面技術(shù)實(shí)現(xiàn)。

海洋工程結(jié)構(gòu)物在流場(chǎng)中的CFD數(shù)值模擬，普遍采用移動(dòng)網(wǎng)格的方法，其實(shí)現(xiàn)形式主要有5種：一是伴隨著物體的運(yùn)動(dòng)，移動(dòng)整個(gè)計(jì)算域的網(wǎng)格，隨之自由面位置不斷變化，因此需要指定外域的邊界條件。此方法只適用于浮體在無限域中的運(yùn)動(dòng)；二是緊貼物體周圍的網(wǎng)格與物體一起做剛體運(yùn)動(dòng)，較遠(yuǎn)部分采用固定網(wǎng)格，這兩個(gè)域之間的調(diào)整區(qū)采用變形網(wǎng)格。此方法的優(yōu)點(diǎn)是易于指定外域邊界條件，而缺點(diǎn)是只能考慮浮體單一方向的平動(dòng)，且在調(diào)整區(qū)的網(wǎng)格變行和扭曲不能過大，即浮體的運(yùn)動(dòng)必須是微小幅值的；三是把方法二中的調(diào)整域改成重構(gòu)網(wǎng)格域，此方法的優(yōu)點(diǎn)是可以應(yīng)用于任意浮體運(yùn)動(dòng)，不足是需采用非結(jié)構(gòu)網(wǎng)格進(jìn)行網(wǎng)格劃分和重構(gòu)，且必須保證相當(dāng)?shù)木龋碌闹貥?gòu)網(wǎng)格還采用原來的求解插值方法，易引入額外的誤差；四是動(dòng)計(jì)算域和滑移交界面的聯(lián)合使用，物體周圍的隨體網(wǎng)格和物體一起做剛體運(yùn)動(dòng)，考慮到浮體的旋轉(zhuǎn)運(yùn)動(dòng)，這個(gè)域的邊界劃分為圓形或者球形，外部網(wǎng)格做壓條運(yùn)動(dòng)。此方法能很好的保證網(wǎng)格質(zhì)量，但它對(duì)網(wǎng)格的劃分要求較高；五是采用重疊網(wǎng)格技術(shù)，浮體周圍的隨體網(wǎng)格與浮體一起做剛體運(yùn)動(dòng)無變形，外部網(wǎng)格采用固定網(wǎng)格，兩個(gè)網(wǎng)格域在浮體附近重疊，這樣浮體可以無限制的運(yùn)動(dòng)，且利于考慮多體運(yùn)動(dòng)，但對(duì)于非結(jié)構(gòu)網(wǎng)格，在重疊區(qū)域的網(wǎng)格耦合和保持守恒特征有一定困難。本文在方法四的基礎(chǔ)上進(jìn)行改進(jìn)，將整個(gè)計(jì)算域由里向外分成4層，所有網(wǎng)格均為結(jié)構(gòu)網(wǎng)格，第一層隨體網(wǎng)格是圓域，可隨浮體一起做任意方向平動(dòng)和轉(zhuǎn)動(dòng)的聯(lián)合運(yùn)動(dòng)，第二層是過渡域，只做平動(dòng)，第三層為矩形網(wǎng)格域，只做平動(dòng)，第四層為靜止網(wǎng)格域（如果外域邊界為動(dòng)邊界條件，則可以有網(wǎng)格運(yùn)動(dòng)，如：造波板邊界條件）。此動(dòng)網(wǎng)格技術(shù)中采用三組交界面，一為圓周邊界的滑移交界面，方便隨體網(wǎng)格做旋轉(zhuǎn)運(yùn)動(dòng)；二為小矩形交界面，用于實(shí)現(xiàn)兩域網(wǎng)格流場(chǎng)信息的過渡；三為動(dòng)網(wǎng)格域與靜網(wǎng)格域的滑移交界面，其中內(nèi)交界面隨動(dòng)網(wǎng)格做平動(dòng)。

1 數(shù)學(xué)模型

1.1 流體運(yùn)動(dòng)控制方程

本文采用VOF法捕捉自由液面，根據(jù)網(wǎng)格單元內(nèi)流體所占體積分?jǐn)?shù)來追蹤自由液面的變化過程。這種方法的優(yōu)點(diǎn)在于可以跟蹤翻轉(zhuǎn)、吞并、破碎、上浪、沖擊等復(fù)雜的非線性自由面現(xiàn)象。

在本文的研究中，計(jì)算域包含空氣和水兩種介質(zhì)，流體為不可壓縮的，忽略粘性，則流場(chǎng)控制方程為：

式中：u,y分別為x,y方向的速度；ρ為流體密度；p為流體壓強(qiáng)；g為重力加速度；aq為第q相流體單元內(nèi)所占的體積分?jǐn)?shù)，q=1，2分別代表空氣和水。

1.2 配有錨鏈系泊的浮式防波堤計(jì)算模型

為了研究浮體在波浪中的運(yùn)動(dòng)和受力，本文對(duì)錨鏈系泊的浮式防波堤進(jìn)行數(shù)值模擬。假設(shè)錨鏈對(duì)稱布置，初始為自然伸長(zhǎng)狀態(tài)，導(dǎo)纜孔位于防波堤兩角點(diǎn)，且只承受拉力，應(yīng)力應(yīng)變服從胡克定律，錨鏈剛度466.0 N/m，忽略其質(zhì)量、阻尼的影響和浮體轉(zhuǎn)動(dòng)引起的轉(zhuǎn)動(dòng)慣量變化。規(guī)定所有變量垂直向上和水平向右為正，彎矩逆時(shí)針為正，以浮式方箱防波堤為例進(jìn)行受力分析。在波浪作用下，防波堤受到重力、流體作用力和錨鏈力，對(duì)防波堤建立運(yùn)動(dòng)方程如下：

式中：x，y分別為防波堤x和y方向上的位移；θ為防波堤的轉(zhuǎn)角；fx，fy分別為防波堤所受流體作用力在x和y方向上的分量，通過聯(lián)立（1）、（2）式和（3）式求出的流體壓力沿浮體濕表面積分獲得；I為防波堤繞重心的轉(zhuǎn)動(dòng)慣量，Mfg為流體作用力對(duì)防波堤重心的力矩，根據(jù)濕表面每個(gè)單元x和y方向上的流體作用力所產(chǎn)生的力矩沿濕表面積分計(jì)算出；Fx錨鏈、Fy錨鏈分別為x和y方向上的錨鏈力，表示為：

Mx和My分別為x和y方向上的錨鏈力對(duì)防波堤重心的力矩，當(dāng)x≥0，y≥0時(shí)：

這里Y為導(dǎo)纜孔y方向坐標(biāo)；Y_cg為重心縱坐標(biāo)；k為錨鏈剛度；X_left，Y_left，X_right，Y_right，X_cg，Y_cg分別為左右角點(diǎn)及重心的全局坐標(biāo)。由以上分析得，聯(lián)立（1）～（7）式可求解浮式防波堤所受的流體作用力、運(yùn)動(dòng)速度、位置及錨鏈力。

2 數(shù)值計(jì)算方法

2.1 造消波技術(shù)

根據(jù)推板造波理論[9]，沖程為X0，圓頻率為ω，平衡位置在原點(diǎn)的造波機(jī)做簡(jiǎn)諧運(yùn)動(dòng)的水平速度為：

在水深為d的2D水槽中，產(chǎn)生的波面高程為：

式中：η（ x,t）為波面的瞬時(shí)高程；k為規(guī)則波的波數(shù)，由色散方程ω2=gktanh（ kd）得出；ki（i=1,2…,n）由色散方程ω2=-gktan（ kd）得出。等號(hào)右邊的第1項(xiàng)為造波板運(yùn)動(dòng)產(chǎn)生的傳遞波，第2項(xiàng)為造波板運(yùn)動(dòng)產(chǎn)生的局部振蕩，隨距造波板距離增大呈指數(shù)衰減，達(dá)到一定距離后可忽略不計(jì)，因此需設(shè)置前端緩沖區(qū)來過濾掉非傳播模態(tài)，得到所需的規(guī)則波波面。

為保證計(jì)算的穩(wěn)定性，造波板的運(yùn)動(dòng)規(guī)律采用如下形式：

根據(jù)王本龍和宋帥[10-11]等的阻尼海綿層消波理論進(jìn)行數(shù)值消波，即在數(shù)值水槽的消波區(qū)設(shè)置海綿層來吸收穿過工作區(qū)的能量，使速度場(chǎng)和壓力場(chǎng)逐漸衰減，水槽邊界處幾乎無波浪反射，以消除反射波對(duì)工作區(qū)的影響，從而保證數(shù)值水槽模擬的準(zhǔn)確性。消波區(qū)的速度和壓強(qiáng)分別為：

式中：下標(biāo)C表示在進(jìn)入海綿層之前的波動(dòng)場(chǎng)；下標(biāo)M表示修正后的波動(dòng)場(chǎng)；λ=λ（x）為與空間位置有關(guān)的從0～1之間的光滑過渡加權(quán)函數(shù)。本文?。?/p>

式中：x1為消波區(qū)左端點(diǎn)的橫坐標(biāo)，xl為海綿層長(zhǎng)度。

2.2 梯形法和二、三階單步數(shù)值計(jì)算方法

一階常微分方程初值問題的梯形數(shù)值解法和二、三階單步法分別具有二階和四階精度，本文分別采用這兩種方法進(jìn)行時(shí)間推進(jìn)，并作比較分析。

設(shè)一階常微分方程初值問題為

梯形法的求解公式為：

二、三階單步求解公式同時(shí)包含二階和三階兩個(gè)單步公式，比單獨(dú)使用二階或三階單步法具有更高的準(zhǔn)確性。該方法的計(jì)算分為三個(gè)階段，需求四個(gè)斜率ki（i=1,2,3,4）。k1的初始值可由t0，u0計(jì)算得到，以后每積分一步，將k4作為下一時(shí)間步k1。該算法的關(guān)鍵步驟如下：

其中：u0表示初始時(shí)刻a的值；tn=tn-1+h；un，un-1分別表示此刻和積分前一步時(shí)刻的函數(shù)值；h表示時(shí)間步長(zhǎng)。在時(shí)間步長(zhǎng)h中，將前一時(shí)間步的x、y方向上的位移x0，y0，角點(diǎn)和重心位置作為錨鏈力和彎矩計(jì)算的初值帶入等式右邊，算出這一時(shí)刻的位移x1，y1和坐標(biāo)，設(shè)定目標(biāo)精度為ε，則當(dāng)且時(shí)認(rèn)為收斂，停止計(jì)算，否則將此刻位移x1，y1和坐標(biāo)賦給x0，y0和初始坐標(biāo)代入等式右邊重復(fù)上述步驟，直到滿足精度要求為止。

在計(jì)算中，將浮式防波堤的運(yùn)動(dòng)方程（5）、（6）和（7）改寫成具有初值問題的一階常微分方程組，即：

式中：u，v，ω1分別為浮式防波堤的水平，垂直方向的線速度和角速度，數(shù)值迭代計(jì)算時(shí)，初值都設(shè)置為0。

2.3 全自由度分塊結(jié)構(gòu)移動(dòng)網(wǎng)格技術(shù)

目前運(yùn)用移動(dòng)網(wǎng)格方法實(shí)現(xiàn)流場(chǎng)中運(yùn)動(dòng)浮體數(shù)值模擬的CFD技術(shù)中，多數(shù)學(xué)者采取的是將縱蕩、垂蕩和縱搖（2D數(shù)值模擬）分別處理或?qū)⒖v蕩、垂蕩兩者之一與縱搖一起聯(lián)合處理。網(wǎng)格都可劃分為結(jié)構(gòu)網(wǎng)格，計(jì)算精度高，但卻無法考慮所有自由度運(yùn)動(dòng)的相互耦合影響。有人提出在隨體圓形區(qū)域網(wǎng)格的外圍設(shè)置非結(jié)構(gòu)三角形網(wǎng)格，通過剛體運(yùn)動(dòng)和網(wǎng)格重構(gòu)的方法模擬浮體上述三種方式的聯(lián)合運(yùn)動(dòng)，但網(wǎng)格再生耗費(fèi)時(shí)間長(zhǎng)，分辨率不易控制，而且再生網(wǎng)格局部可能產(chǎn)生過大的扭曲度或網(wǎng)格形狀發(fā)生突變，引起收斂困難和精度降低。

本文采取全自由度分塊結(jié)構(gòu)移動(dòng)網(wǎng)格技術(shù)，實(shí)現(xiàn)2D數(shù)值水槽中運(yùn)動(dòng)浮體縱蕩、垂蕩和縱搖所有方向上的聯(lián)合運(yùn)動(dòng)，數(shù)值模擬時(shí)，浮式防波堤置于水槽工作區(qū)中央，重心在靜水面以下0.4 m處，如圖1所示。工作區(qū)的動(dòng)靜網(wǎng)格細(xì)劃分情況如圖2所示。定義圍繞浮式防波堤的圓域3為剛體網(wǎng)格，既有平動(dòng)又有轉(zhuǎn)動(dòng)，網(wǎng)格尺度與扭曲度均不變，且與物體的相對(duì)位置不變；結(jié)構(gòu)平動(dòng)網(wǎng)格域2只有平動(dòng)，網(wǎng)格尺度與扭曲度不變；結(jié)構(gòu)平動(dòng)網(wǎng)格域1采用動(dòng)態(tài)網(wǎng)格層變法實(shí)現(xiàn)網(wǎng)格的平動(dòng)變形，且保持扭曲度不變。這里三個(gè)移動(dòng)網(wǎng)格域均采用結(jié)構(gòu)網(wǎng)格，其中域1采用矩形結(jié)構(gòu)網(wǎng)格。在設(shè)置邊界條件時(shí)，域1與工作靜止區(qū)，域1與2以及圓域3與域2分別設(shè)置一對(duì)interface，以保證運(yùn)動(dòng)時(shí)各個(gè)子區(qū)域互不干擾，同時(shí)相鄰區(qū)域的流場(chǎng)信息通過插值方式實(shí)現(xiàn)相互傳遞。此外，定義域1的左右邊界做垂向運(yùn)動(dòng)，同域1的垂向速度，而上下邊界做水平運(yùn)動(dòng)，同域1的水平速度，來實(shí)現(xiàn)防波堤縱蕩、垂蕩和縱搖三種方式的聯(lián)合運(yùn)動(dòng)。

圖1 2D數(shù)值水槽示意圖Fig.1 The sketch of 2D numerical tank

圖2 結(jié)構(gòu)移動(dòng)網(wǎng)格的分塊劃分Fig.2 Meshing of multiply structured grids

3 數(shù)值試驗(yàn)與結(jié)果分析

按以下要求劃分計(jì)算域網(wǎng)格以保證數(shù)值模擬的精度：波浪傳播方向上的網(wǎng)格數(shù)量必須充足，以避免數(shù)值阻尼引起入射波浪衰減；造波區(qū)、前端緩沖區(qū)和工作區(qū)水平網(wǎng)格尺寸約為波長(zhǎng)的1%，尾端消波區(qū)可適當(dāng)增大，這樣既可減少計(jì)算時(shí)間，又可配合阻尼海綿達(dá)到更好的消波目的；為精確捕捉自由面的變化過程，垂直方向的網(wǎng)格在自由面附近約5%波高，距離越遠(yuǎn)網(wǎng)格尺寸越大；精細(xì)劃分浮體附近的網(wǎng)格，可更好地模擬物體和流體之間的相互耦合作用。

3.1 數(shù)值造消波的驗(yàn)證

為驗(yàn)證數(shù)值水槽的造消波效果，本文按比例尺1: 10模擬波高為0.4 m，周期2.21 s（波長(zhǎng)7.62 m）的2D行進(jìn)波。根據(jù)Ursell和Dean[12]的線性造波機(jī)理論，計(jì)算得出沖程X0=0.205。x=2 m、22 m、62 m處的波面時(shí)間歷程，如圖3所示。由圖中可測(cè)得三處波高分別為0.45 m，0.39 m，0.001 4 m，周期約2.20 s。在x=2 m處，造波區(qū)內(nèi)的波浪包含傳播模態(tài)和非傳播模態(tài)，尚未完全衰減，故波高偏大；在x=22 m處，工作區(qū)內(nèi)的傳播模態(tài)已基本衰減，因此模擬波高近似理論值；在x=62 m處，消波區(qū)內(nèi)的波面經(jīng)過阻尼海綿已基本消減，因此近似靜水面，消浪效果良好。

3.2 全自由度分塊結(jié)構(gòu)移動(dòng)網(wǎng)格技術(shù)的驗(yàn)證

本文通過基于流固“全耦合”原理的全自由度分塊結(jié)構(gòu)移動(dòng)網(wǎng)格技術(shù)分別采用梯形法和二、三階單步法進(jìn)行比較和驗(yàn)證，這里以一浮式方箱防波堤為例，模型尺寸和環(huán)境要素如前所述。浮式方箱防波堤在靜水中浮力等于和小于重力時(shí)，垂蕩速度和流體垂向作用力的時(shí)間歷程曲線如圖4～5所示。其中圖5所示的防波堤質(zhì)量為1.787e3 kg。由圖4（a）可見，采用梯形法和二、三階單步法計(jì)算得到的浮體垂蕩速度隨時(shí)間變化趨勢(shì)相同，但二、三階單步法的曲線較梯形法更加光滑均勻，二者隨著時(shí)間的推進(jìn)，都逐漸衰減趨于零，這是由于初始階段流體和空氣存在相對(duì)壓強(qiáng)使得浮力略小于重力，浮體做微小幅值的垂蕩運(yùn)動(dòng)。由圖4（b）可見，由于初始階段流體和空氣存在密度差，使得流體垂向作用力為13 654.70 N，與重力差值為13 704.57-13 654.70=49.87 N，因此防波堤所受的流體垂向作用力大致以浮力13 654.70為中心線做上下振蕩，兩種方法對(duì)應(yīng)的曲線變化趨勢(shì)基本相同，但與梯形法相比，二三階單步法的振動(dòng)幅值略微偏小，可見二、三階單步法的數(shù)值模擬更加精確。由圖5（a）可見，當(dāng)浮力小于重力時(shí)，浮式方箱防波堤做垂向往復(fù)衰減運(yùn)動(dòng)，兩種方法對(duì)應(yīng)的浮體垂蕩速度幾乎完全重合。由圖5（b）中可見，初始階段由于浮力小于重力，浮式方箱防波堤垂直向下運(yùn)動(dòng)，隨著浮體吃水增加，浮力也隨之增加，運(yùn)動(dòng)由變加速到變減速再到變加速，如此循環(huán)下去，最終流體垂向作用力大致以重力17 530.47 N為中心線做衰減簡(jiǎn)諧運(yùn)動(dòng)。圖中兩條曲線趨勢(shì)大致相同，二、三階單步數(shù)值迭代方法光滑度略高于梯形法。從圖4和圖5的浮式防波堤水動(dòng)力性能分析，可以看出本文的流固全“耦合”方法和全自由度分塊結(jié)構(gòu)移動(dòng)網(wǎng)格技術(shù)算法穩(wěn)定，收斂良好，精度滿足要求，與工程實(shí)際情況相符，且二、三階單步數(shù)值迭代法的數(shù)值結(jié)果要更加均勻和光滑。

圖3 不同位置的波面時(shí)間歷程Fig.3 Time series of wave profile at x=2 m,22 m,62 m

圖4 浮式方箱防波堤在靜水中浮力等于重力時(shí)的垂蕩速度和流體作用力Fig.4 Heaving velocity of the floating breakwater and fluid force when buoyancy force is equal to force of gravity in still water

圖5 浮式方箱防波堤在靜水中浮力小于重力時(shí)的垂蕩速度和流體作用力Fig.5 Heaving velocity of the floating breakwater and fluid force when buoyancy force is less than force of gravity in still water

3.3 系泊式浮式方箱防波堤在波浪條件下的水動(dòng)力特性分析

為了研究系泊式海上建筑物在波浪作用下的動(dòng)力過程，這里對(duì)近年來海洋中日益增多錨鏈系泊的浮式方箱防波堤運(yùn)動(dòng)過程進(jìn)行數(shù)值模擬，關(guān)于模型參數(shù)和海況條件如前文所述。圖6分別描述了規(guī)則波中，系泊式方箱防波堤的運(yùn)動(dòng)速度，所受錨鏈力（矩）和流體作用力（矩）（均以無因次量表示）。其中為進(jìn)一步驗(yàn)證基于“全耦合”方法所采用的全自由度分塊結(jié)構(gòu)移動(dòng)網(wǎng)格技術(shù)對(duì)實(shí)現(xiàn)結(jié)構(gòu)所有方向上運(yùn)動(dòng)的準(zhǔn)確性，本文將數(shù)值計(jì)算所得的防波堤速度和錨鏈力分別與基于勢(shì)流理論的邊界元法所得結(jié)果進(jìn)行比較。由圖6（a）可見，本文數(shù)值方法與邊界元法所得結(jié)果吻合較好，且在三個(gè)方向上梯形法與二、三階單步法計(jì)算出的浮式防波堤的速度基本相同。其中水平方向和垂直方向上的線速度均只呈現(xiàn)出波頻約2.86 rad/s的振蕩特性，而繞重心的角速度的時(shí)程曲線則呈現(xiàn)出波頻約為2.86 rad/s和低頻約為0.348 rad/s相結(jié)合的振蕩特性，可見縱蕩和垂蕩運(yùn)動(dòng)主要受行進(jìn)波的影響，而縱搖運(yùn)動(dòng)與行進(jìn)波和浮體固有特性兩者都有關(guān)。由圖6（b）可見，對(duì)錨鏈力的計(jì)算，本文“全耦合”方法與邊界元法所得結(jié)果除了峰值和谷值略有差別，其他都吻合較好。隨波浪穩(wěn)定傳播，防波堤在水平方向、垂直方向和縱搖方向受到的錨鏈力（矩）都以波頻約2.86 rad/s做振蕩運(yùn)動(dòng)，且梯形法和二、三階單步法結(jié)果基本相同。由圖6（c）可見，防波堤水平方向上和繞重心的流體作用力（矩）大致分別以-0.05和0為平衡點(diǎn)做有規(guī)律的往復(fù)波頻振動(dòng)，且波谷的絕對(duì)值要大于波峰，其與波浪的傳播方向所引起的漂移運(yùn)動(dòng)有關(guān)；防波堤垂直方向上所受的流體作用力大致以1.94（約為14 438 N）為中心線做往復(fù)波頻振動(dòng)，此中心線14 438 N略大于防波堤靜水平衡時(shí)所受的浮力13 704 N，這是因?yàn)殄^鏈線對(duì)防波堤的垂直向下的作用使得防波堤下沉，浮力也隨之增加，由圖可知梯形法與二、三階單步法的曲線趨勢(shì)基本一致。

圖6 波浪中防波堤的運(yùn)動(dòng)速度、錨鏈力和流體作用力的時(shí)間歷程Fig.6 Velocity component,mooring forces and fluid force of floating breakwater in wave

4 結(jié)語

本文采用耦合求解流體動(dòng)力學(xué)控制方程和浮式結(jié)構(gòu)運(yùn)動(dòng)方程的方法，在具有數(shù)值造消波功能的2D數(shù)值水槽中，實(shí)現(xiàn)了海上浮式結(jié)構(gòu)物在波浪中運(yùn)動(dòng)過程的數(shù)值模擬。以配有錨鏈系泊系統(tǒng)的浮式方箱防波堤作為算例，研究了浮式防波堤的水動(dòng)力性能。在數(shù)值計(jì)算中，采用梯形法和二、三階單步數(shù)值迭代法進(jìn)行浮體運(yùn)動(dòng)參數(shù)的時(shí)間推進(jìn)，并且采用全自由度分塊結(jié)構(gòu)移動(dòng)網(wǎng)格技術(shù)，實(shí)現(xiàn)了浮體所有方向的聯(lián)合運(yùn)動(dòng)，從而保證了在計(jì)算過程中網(wǎng)格不發(fā)生任何扭曲現(xiàn)象，計(jì)算時(shí)間、網(wǎng)格劃分和精度要求都得到了較好的控制。主要得出以下結(jié)論：

（1）對(duì)靜水中自由漂浮的浮式方箱防波堤，當(dāng)浮力分別等于和小于重力時(shí)，用梯形法和二、三階單步數(shù)值迭代法計(jì)算得浮體水動(dòng)力特性變化趨勢(shì)基本一致，結(jié)果收斂，但二、三階單步法得出的結(jié)果較梯形法更加平順、光滑，從而驗(yàn)證了流固耦合方法和全自由度分塊結(jié)構(gòu)移動(dòng)網(wǎng)格技術(shù)的可行性。

（2）建立2D數(shù)值波浪水槽，根據(jù)推板造波理論在水槽中模擬規(guī)則波浪，經(jīng)驗(yàn)證工作區(qū)的波浪條件滿足理論值的精度要求，且消波區(qū)的數(shù)值海綿層消浪效果較為理想。

（3）在已建的波浪水槽工作區(qū)中央放入系泊式浮式方箱防波堤模型，采用本文數(shù)值方法對(duì)其水動(dòng)力特性進(jìn)行分析，并將結(jié)果與基于勢(shì)流理論的邊界元法進(jìn)行比較。防波堤的縱蕩和垂蕩運(yùn)動(dòng)主要為波頻運(yùn)動(dòng)，而縱搖運(yùn)動(dòng)為波頻運(yùn)動(dòng)和與防波堤固有低頻特性相關(guān)的低頻運(yùn)動(dòng)相結(jié)合的運(yùn)動(dòng)方式，在波浪中所受的錨鏈力（矩）的時(shí)程曲線均表現(xiàn)為波頻振蕩趨勢(shì)，且結(jié)果與邊界元分析所得吻合較好。

（4）系泊式浮式方箱防波堤在波浪作用條件下，水平方向上和繞重心的流體作用力（矩）大致分別以-0.05和0為平衡點(diǎn)做有規(guī)律的往復(fù)波頻振動(dòng)，且波谷的絕對(duì)值要大于波峰，其與波浪在傳播方向所引起的漂移運(yùn)動(dòng)有關(guān)；垂直方向上的流體作用力的中心線值比靜水平衡條件下浮體所受的浮力略大，原因是錨鏈對(duì)防波堤的垂向作用使得防波堤下沉，浮力也隨之增加。

（5）二、三階單步求解公式在數(shù)值原理上比單個(gè)使用二階或三階單步法具有更高的準(zhǔn)確性，這一點(diǎn)也在靜水條件下浮體的運(yùn)動(dòng)中得到驗(yàn)證，但在波浪中它與梯形法的計(jì)算結(jié)果基本一致，數(shù)值收斂情況和精度也基本相同，因此防波堤的動(dòng)力特性運(yùn)用梯形法計(jì)算即可。

以上計(jì)算結(jié)果表明，本文采用的流固耦合方法和移動(dòng)網(wǎng)格技術(shù)可以用來分析一系列海上浮式結(jié)構(gòu)物在波浪條件下的水動(dòng)力問題，但在數(shù)值模擬系泊系統(tǒng)時(shí)做了一些假定，并且未考慮錨鏈和浮體的耦合作用以及錨鏈對(duì)流場(chǎng)的影響。因此，如何更加精確、更加全面地計(jì)算結(jié)構(gòu)在波浪中的動(dòng)力問題還有待于進(jìn)一步研究。

[1]林兆偉,朱仁傳,繆國(guó)平.甲板上浪問題的二維數(shù)值模擬[J].船舶力學(xué),2009,13(1):1-8. Lin Z W,Zhu R C,Miao G P.2-D Numerical simulation for green water on oscillating ships[J].Journal of Ship Mechanics,2009,13(1):1-8.

[2]郭海強(qiáng).船舶性能數(shù)值水池的研究及其若干應(yīng)用[D].上海:上海交通大學(xué)船舶海洋與建筑工程學(xué)院,2009. Guo H Q.The development and some applications of ship performance numerical tank[D].Master thesis,Shanghai Jiao Tong University,School of Naval Architecture,Ocean and Civil Engineering,Shanghai,2009.

[3]龔丞,朱仁傳,繆國(guó)平,范菊.基于CFD的高速船甲板上浪載荷的工程計(jì)算方法[J].船舶力學(xué),2014,18(5):524-531. Gong C,Zhu R C,Miao G P,Fan J.A numerical method of the simulation of green water on the deck of a vessel[J].Journal of Ship Mechanics,2014,18(5):524-531.

[4]Hadzic I,Hennig J,Peric M,Xing-Kaeding Y.Computation of flow-induced motion of floating bodies[J].Applied Mathematical Modelling,2005,29:1196-1210.

[5]Nielsen K B,Stefan M.Numerical prediction of green water incidents[J].Ocean Engineering,2004,31(3-4):363-399.

[6]Yang Q Y,Qiu W.Numerical simulation of water impact for 2D and 3D bodies[J].Ocean Engineering,2012,43:82-89.

[7]Kleefsman K M T,Fekken G,Veldman A E P,et al.A Volume-of-Fluid based simulation method for wave impact problems[J].Journal of Computational Physics,2005,206(1):363-393.

[8]沈玉稿,楊建民,李欣,等.風(fēng)機(jī)基礎(chǔ)所受波浪抨擊力的數(shù)值模擬和實(shí)驗(yàn)研究（英文）[J].船舶力學(xué),2013,17(9): 1009-1020. Shen Y G,Yang J M,Li X,Peng T.Numerical and experimental investigations of the impact force on offshore wind turbine pile[J].Journal of Ship Mechanics,2013,17(9):1009-1020.

[9]王文華,王言英.圓柱在波浪中入水的數(shù)值模擬[J].上海交通大學(xué)學(xué)報(bào),2010,44(10):1393-1399. Wang W H,Wang Y Y.Numerical study on cylinder entering water in wave[J].Journal of Shanghai Jiaotong University, 2010,44(10):1393-1399.

[10]王本龍,劉樺.一種適用于非均勻地形的高階Boussinesq水波模型[J].應(yīng)用數(shù)學(xué)和力學(xué),2005,26(6):714-722. Wang B L,Liu H.Higher order boussinesq-type equations for water waves on uneven bottom[J].Applied Mathematics and Mechanics,2005,26(6):714-722.

[11]宋帥,尤云祥,魏崗.孤立波與直墻式多孔介質(zhì)結(jié)構(gòu)相互作用數(shù)值分析[J].海洋工程,2007,25(4):7-14. Song S,You Y X,Wei G.The interaction of the solitary wave with a vertically walled porous structure[J].Ocean Engineering,2007,25(4):7-14.

[12]Ursell F,Dean R G,Yu Ys.Forced small-amplitude water waves:a comparison of theory and experiment[J].Fluid Mechanics,1960,7(Pt1):33-52.

[13]FLUENT INC.Fluent 6.2 Manual[M].Fluent Inc.,2004.

[14]王福軍.計(jì)算流體動(dòng)力學(xué)分析[M].北京:清華大學(xué),2004. Wang F J.Computational Fluid Dynamics Analysis[M].Beijing:Tsinghua University,2004.

[15]施吉林,張宏偉,金光日.計(jì)算機(jī)科學(xué)計(jì)算[M].北京:高等教育出版社,2005. Shi J L,Zhang H W,Jin G R.Computer Science Computing[M].Beijing:Higher Education Press,2005.

Hydrodynamic performance of the floating body based on CFD fluid-body interaction theory

MA Zhe1,CHENG Yong2,ZHAI Gang-jun1
(1.Deepwater Engineering Research Center,Dalian University of Technology,Dalian 116024,China;2.Ship and Ocean Engineering School,Jiangsu University of Science and Technology,Zhenjiang 212003,China)

Using fluid-body interaction mechanics,the motion proceeds of a floating body was numerically simulated in regular wave,which was implemented in a wave making and absorbing numerical flume.As the application case of a box-type floating breakwater with mooring systems,dynamic characteristics of the floating body among the trapezoid,2,3order single step method and BEM based on potential flow theory were compared.The results indicate that the trapezoid and 2,3order single step method which were in good agreement with the BEM results has prospective accuracy and convergence in regular wave.This paper proposed an new technique of dynamic mesh with full 2D-3 degrees of freedom multiply structured grids, which realized the joint movement of all directions,and made sure that grids had not any distortion.The computation time,precision and meshing remained fairly well-contained.

fluid-body perfectly coupling;trapezoid;2,3order single step method; dynamic mesh with full freedom multiply structured grids

TV131.2

A

10.3969/j.issn.1007-7294.2017.08.003

1007-7294（2017）08-0950-10

2016-12-22

國(guó)家自然科學(xué)基金應(yīng)急管理項(xiàng)目（51651902）；國(guó)家自然科學(xué)基金青年基金項(xiàng)目（51609109）；中央高?；究蒲袠I(yè)務(wù)費(fèi)（DUT16QY22）；江蘇省自然科學(xué)青年基金（BK20160556）

馬哲（1983-），男，博士，E-mail:deep_mzh@dlut.edu.cn；程勇（1986-），男，博士。