基于one-fluid模型的波浪-浮體相互作用數值模擬研究

田昀艷，王收軍，陳松貴，王家偉，胡原野

(1.天津理工大學 機電工程國家級實驗教學示范中心，天津 300384；2.交通運輸部天津水運工程科學研究所，天津 300456)

在波浪與結構相互作用的研究中，通常采用物理模型實驗和數值模擬相結合的方法。近年來，隨著高性能計算機技術的發展和數值離散方法的成熟，構建數值波浪水槽并研究水波與結構物的相互作用問題，已成為海岸和海洋工程領域不可缺少的研究方法。

在波浪力學領域，波浪與結構相互作用的建模一直是通過基于勢流理論和粘性流理論的數值波浪水槽方法發展起來的。由于勢流理論是建立在不可壓縮和無粘流動假設的基礎上的，該方法不能模擬粘性效應和破碎波。因此，在實際工程應用中，Navier-Stokes方程被用于模擬波浪與結構物相互作用問題。Park J.C.和Uno Y[1]等人，在粘性三維數值波浪水槽中模擬了規則波、不規則波以及多向隨機波，并討論了在特定波浪作用條件下船舶的水動力特性；Hans Bihs[2]等人采用Level Set法追蹤自由表面，建造三維粘性數值波浪水槽，對波浪在傾斜海灘上傳播及破碎過程、垂直圓柱上的波浪力、孤立波與矩形橋臺之間的相互作用進行了模擬；孫宏月等[3]基于緊致差值曲線方法建立了可壓縮兩相流數值波浪水槽模型，模擬了波浪對結構物的沖擊過程；余陽等[4]基于OpenFOAM平臺，采用重疊網格技術建立了黏性數值波浪水池，驗證了重疊網格在二維和三維兩相流體區域中流體與結構相互作用的準確性和可靠性；張弈澤等[5]采用Flow3D軟件建立數學模型對越浪量進行了模擬研究，對防浪墻頂高程的合理性進行分析，并對不同形式的海堤護面結構進行數值模擬研究；焦方騫等[6]基于OpenFOAM軟件建立了波浪與結構物相互作用的三維數值波浪水槽，模擬了復合筒型基礎在不同波浪條件下的波浪載荷；Patankar[7]以歐拉方式將剛體看作流體，通過兩相流方法對剛體和流體耦合建模問題進行簡化，剛體的慣性完全由流體求解器承擔，但是施加的剛性力是在離散水平上進行的。上述基于NS方程對波浪與結構相互作用的研究大多是基于兩相流的水平集函數，不能用于解決浸沒結構多相多界面問題的演化。與上述大多數可用的浸入邊界方法不同，Yang[8]提出了一種模擬非線性粘性波浪產生全物理過程的one-fluid模型。通過提出一種基于線性最小二乘法求解剛體的運動約束，在連續介質水平上明確表達了浸沒剛體的力,不需要諸如質心、質量和慣性動量之類的信息，只需要知道剛體的密度。這種方法將剛體與流體耦合問題重新表述為一個帶有附加剛體運動約束的三相浸入流動問題。該模型可以用于處理大密度比、多相流問題，已成功應用于水下剛性板的顫振運動[9]以及波能轉換器(WEC)的數值研究[10]等。

鑒于以上背景，本文基于one-fluid方程建立了二維粘性數值波浪水槽，進行了不規則波的模擬并且研究了浮式結構在不規則波作用下的動力響應。第一節介紹了one-fluid方程的數值求解過程以及自由表面的處理方法等；第二節對不規則波進行了數值模擬，并研究了不規則波與浮式結構之間的相互作用；第三節探討了不規則波作用，波高、周期、浮體寬度與波長的相對寬度以及浮體的相對吃水對浮式結構物的垂蕩運動響應的影響；第四節給出了相關的結論。

1 數學模型

圖1 流體與剛體相互作用問題的拉格朗日力學耦合與one-fluid框架的說明

1.1 控制方程

在one-fluid模型中，一個方程里同時考慮了多個不可壓縮相，包括空氣、水和結構。與其它剛體-流體耦合問題處理方法不同之處在于它利用線性最小二乘法，提出了一種用于剛體運動約束的剛體投影算子。該算子在連續介質水平上明確地表達了浸沒剛體的力，不需要諸如質心、質量和慣性動量之類的信息。該模型可以用于處理大密度比、多相流問題。傳統模型中，耦合力是通過轉矩tc和外力fc來實現的，而在one-fluid模型中，耦合是通過界面上的平衡牽引矢量來實現的。圖1說明了流體與剛體相互作用問題傳統解決方法和one-fluid模型之間的差異。

one-fluid模型中，一個方程中同時考慮了多個不可壓縮相，包括空氣、水和結構。在整個區域Ω=Ωa∪Ωw∪Ωr(在沒有表面張力的情況下)三相系統的線性動量守恒方程和連續性方程用一個方程表示如下

(1)

式中：u是速度矢量(u，v)，ρ=ρrHr+ρwHw+ρaHa，其中下標a、w、r分別表示空氣、水和剛體的相，H是不同相的Heaviside分布。式(1)中的體積力f定義在每個相位上表示如下

(2)

其中應變率張量D(u)=(u+(u)T)/2，μ表示動力粘度系數。P(u)是剛性速度投影算子。

1.2 自由液面捕捉

one-fluid模型的表述依賴于連續體不同相之間界面的正確識別，可以通過在界面上引入合適的跳轉條件進行。在采用浸沒邊界方法的情況下，可以通過界面的平滑正則化避免明確使用跳轉條件。水、空氣和結構這三個不同的區域由三個Heaviside函數表示。通過三個Heaviside函數來描述界面的具體步驟如下：

忽略剛體，與兩相流問題類似，引入了水平集函數來識別空氣和水之間的界面。帶符號的距離函數φ給出了最接近界面的距離Γ，并且通過符號的變化來區分這兩個相位，因此滿足下式

Γ={x∈R2|φ(x,t)=0}

(3)

(4)

式中：phase1表示水；Γ表示自由液面；phase2表示空氣。當界面Γ在外部生成的速度場u下移動時，將獲得水平集函數的對流方程如下

(5)

①引入剛性Heaviside函數Hr(x,t)=Hr(X,0)來識別剛體區域。其中：X表示拉格朗日坐標；歐拉坐標x=Ψ(X，t)；Ψ表示從X到x的映射；

②“剛體”、“水”、“空氣”所占據的區域定義如下

Ωrigid={x∈R2|Hr(x,t)=1}
Ωwater={x∈R2|Hφ(x,t)-Hr(x,t)=1}

(6)

Ωair={x∈R2|1-Hφ(x,t)-Hr(x,t)=1}

使用兩個Heaviside函數表示三相的圖示如圖2所示。

圖2 使用兩個Heaviside函數表示三相的圖示

1.3 Heaviside函數的求解

①平滑流體Heaviside函數：

Heaviside函數的演化遵循Ha(x,t)=H(φ(x,t))，其近似值可以表示如下

(7)

其中∈代表拖尾帶寬的參數。因此，在涂抹界面上，Heaviside函數的近似值從零變為一，從而描述了從一個相到另一個相的平滑過渡。

②平滑剛體Heaviside函數：

在初始時間t=0時，引入卷積程序來規定任意的Heaviside函數。

步驟1：給定笛卡爾網格，將值1賦給描述剛體的多邊形，該多邊形是非光滑離散的Heaviside函數Hr。

步驟3：將非平滑的Heaviside函數與光滑的Dirac Delta函數卷積得到平滑的Heaviside函數

(8)

圖3 邊界和松弛區示意圖

1.4 邊界條件

多相流數值水槽模型采用速度邊界入口法和松弛消波法來實現造波和消波功能，邊界和松弛區如圖3所示。模型左側為速度入口邊界條件，通過自定義函數給定入口邊界的速度和波面變化；計算域底部為滑移壁面邊界條件；水槽右端為封閉邊界條件。圖中松弛區Ⅰ的功能是協助波浪的產生以及吸收從結構物反射回來的波浪；松弛區Ⅱ的功能是消除出口邊界的反射波。波浪在造波區和消波區的速度、波面以及壓力的計算公式見文獻[11],不規則波水質點的波面方程、垂直速度與水平速度方程見文獻[12]。

1.5 多相流求解器求解過程

(1)通過卷積初始化流體速度u0、水平集方程φ和Heaviside函數H。

(2)開始非定常計算的時間步循環。

1)在不考慮粘性的情況下，通過one-fluid動量方程來計算中間速度場u*，ρ(u*-un/Δt)=RHSn，其中RHSn=-ρ(un)un+ρg，上標n表示時間步長。

2)計算區域Ωr的剛性力fn：

①線性最小二乘法：通過最小化‖P(u*)-u*‖2來找到剛性速度場P(u*)；

②將剛體力更新為：fn=ρ[p(u*)-un/Δt]-RHSn。

3)計算粘性項：fn=μ·2μD(un)。

4)考慮粘性，通過one-fluid動量方程來計算中間速度場u**：ρ(u**-un/Δt)=RHSn+fn。

5)用Neumann邊界條件計算壓力以滿足不可壓縮約束條件：·(pn/ρ)=-·u**/Δt。

6)應用速度校正：un+1=u**-Δt·pn/ρ。

7)計算水平集對流方程并重新初始化。

8)基于高階插值的剛性Heaviside函數對流化。

9)更新下一時間步的速度un+1和密度ρn+1。

2 數值水槽模型驗證

為了驗證數值水槽模型的有效性，本文對不規則波進行了模擬，并研究了浮式結構在不規則波作用下的動力響應。

表1 不規則波模擬目標參數

2.1 不規則波浪的數值模擬

為了驗證one-fluid模型的準確性，應用合田改進后的JONSWAP譜對在交通運輸部天津水運工程科學研究所大比尺波浪水槽實驗室實測的不規則波數據進行了分析和數值再現。選擇了三組不同工況下浪高儀所記錄的波面數據進行比較分析，模擬工況如表1所示。數值波浪水槽計算模型的設置與交通運輸部天津水運工程科學研究所的大比尺波浪水槽一致，如圖4所示。模型尺寸為450 m×12 m，消波段長60 m。計算網格尺寸為3 800×100，時間步長取0.02 s。在數值水槽距離造波邊界150 m處設置虛擬浪高儀G1來監測波面變化。

圖4 數值波浪水槽示意圖

圖5給出了按JONSWAP譜進行模擬得到的x=150 m處的不規則波面時歷曲線與實測不規則波面時歷曲線。對波浪時歷進行統計分析可得有效波高及有效周期的計算值見表2，所得結果與實驗值吻合較好。還可以得出one-fluid模型對于波高為Hs=1.55 m和Hs=1.8 m的模擬結果要好于波高為Hs=1 m的，隨著波高和周期的增大，數值計算所得結果與實驗值的誤差越小。

5-a 工況a

表2 計算值和實驗值的比較

表3 計算所得波能譜和實驗所得波能譜的比較

圖6給出了波面經過傅里葉變換之后得到的頻譜與目標靶譜的對比圖。從圖中可以看出，得到的頻譜與目標靶譜基本吻合。表3給出了由實驗所得的波能譜與one-fluid模型計算所得波能譜的值，可以看出one-fluid模型模擬的譜峰頻率與實驗所得的譜峰頻率兩者吻合較好，表明本文模擬的不規則波具有一定的準確性。

6-a 工況a

2.2 預測不規則波與結構物相互作用

為了驗證所建立的數學模型預測結構與不規則波相互作用的準確性，本文數值模擬了William等人[13]的案例。利用one-fluid模型的“剛體”流固耦合特性對浮體結構進行了建模。

本驗證計算所選取的數值水槽幾何尺寸、波浪參數、浮式結構尺寸以及距水槽的初始位置等均參照William等人的數學模型設置。模擬過程中的波面數據來源于大西洋海能測試站點浮標所記錄的三個不同時段的數據，時間分別為2010年12月17日10:30、2010年10月10日4:30、2011年2月9日17:00。空氣、水和浮式方箱的密度分別取1.0 kg/m3、1 000 kg/m3和750 kg/m3。空氣和水的粘度分別為17.9×10-6Pa·s和1.0×10-3Pa·s。本文僅針對結構物的垂蕩運動進行研究分析。

圖7顯示了不規則波與結構物之間在多個時間步長上的相互作用，包括波浪剖面和結構物的動態響應。

7-a t=120 s時

為了評估模型的準確性，將數值波浪水槽的垂蕩運動動態響應與結構水動力分析得出的解析解以及William等人數值分析所得到的結果進行了比較，結果如圖8所示。可以看出：本文數值模擬得到的結果與William等人數模的結果在頻率和幅值方面吻合度都很高。而由one-fluid模型得到的結果與結構物水動力分析得到的結果在頻率方面匹配的很好，但兩者的響應幅度所得結果有所不同。這是由于數值模型從穩定狀態開始，因此在模擬初始部分與水動力計算軟件AQWA得到的結果存在很大差異。這種差異可能是由于浮式結構物的位置變化和傾斜情況的影響，從而導致整個運行過程中的潤濕面發生變化所造成的。在采用one-fluid耦合模型時，考慮了這一因素對浮體水動力荷載的影響，而由結構水動力分析得到的解析解中忽略了這一影響。此外，結構物本身壁面或波浪水槽的數值模型中存在粘性非線性影響，會導致阻尼力增加，然而AQWA是基于勢流理論的水動力計算軟件，在計算過程中不考慮流體粘性，結果有一定誤差，這可能導致one-fluid模型與解析解之間的差異。

8-a 2010年12月17日10:30

3 數值水槽模型應用

入射波的周期、波高不同，波長也不同，對浮體的受力也不盡相同。結構物在與波浪相互接觸時，不僅會受到波浪的沖擊載荷，而且由于波浪的形態變化，結構物浮心位置也會隨之發生變化，因此結構物在波浪中的運動響應是由多種因素共同作用產生的。

表4 波高、周期的取值

本節探討了在不規則波作用下，波高、周期、相對寬度(浮箱寬度與波長的比值W/L)以及浮箱相對吃水對浮式方箱垂蕩動態響應的影響。二維數值水槽模型長28 m，高1 m，其中水深0.6 m,消波段長5 m。水平方向網格為896個單元，垂向網格為32個單元，時間步長取0.015 s。數值模擬波要素的選取見表4～表6。

表5 相對寬度的取值

表6 浮箱相對吃水的取值

從圖9中可以看出，浮體的垂蕩值隨著波高的增大而增大。在不同的波高下，浮體的垂蕩曲線呈現出一定的周期性。當波高Hs=0.07 m時，垂蕩峰值0.048 m；當波高Hs=0.08 m時，垂蕩的峰值為0.056 m；當波高Hs=0.1 m時，垂蕩的峰值為0.075 m。通過對比三種工況可以看出，隨著入射波高的增加，浮式方箱垂蕩運動響應也增大。

從圖10中可以看出，在不同的周期下，浮體的垂蕩曲線也呈現出一定的周期性。運動大都在30 s后保持穩定，周期Ts=0.91 s、1.1 s、1.28 s對應的垂蕩峰值分別為0.086 m、0.075 4 m、0.06 m。從圖中可以觀察到，入射波周期Ts=0.91 s的情況下，浮體的垂蕩峰值最大，然后隨著入射波周期變長，運動響應峰值減小。

從圖11中可以看出，浮箱的寬度與波長之比W/L=0.163 9、0.294、0.361 4所對應的垂蕩峰值分別為0.078 m、0.058 4 m、0.056 m。可以看出相對寬度W/L=0.163 9時浮體垂蕩動態響應最大，垂蕩運動響應隨著相對寬度的增加逐漸減小。

從圖12中可以看出，浮體的吃水深度h=0.1 m、0.12 m、0.14 m時，對應的垂蕩的峰值分別為0.117 m、0.095 m、0.082 m。可以看出吃水深度h=0.1 m時浮體垂蕩運動響應最大，隨著浮體吃水深度的增大，垂蕩運動響應峰值減小。

圖9 不同波高下,高a=0.18 m,寬b=0.3 m的浮式方箱垂蕩運動響應曲線

圖10 不同周期下,高a=0.18 m,寬b=0.3 m的浮式方箱垂蕩運動響應曲線

圖11 不同相對寬度下,高a=0.18 m,寬b=0.3 m的浮箱垂蕩運動響應曲線

圖12 不同吃水深度下，高a=0.2 m，寬b=0.3 m的浮式方箱垂蕩運動響應曲線

4 結論

(1)本文使用one-fluid數學模型建立了二維粘性數值波浪水槽。通過對濾波后的實測波與one-fluid模型輸出值在時域和頻域的比較，發現它們在頻率方面吻合的非常好。

(2)將數值波浪水槽模型得到的結構垂蕩運動動力響應與結構水動力分析得到的解析解以及william的數值研究結果進行了比較，結果表明三種結果吻合度較高，特別是在頻率方面吻合度很高。

(3)研究了不同波高、周期以及浮式結構物的大小對不規則波作用下結構物垂蕩動態響應的影響。結果表明隨著波高的增加，浮式方箱垂蕩運動響應也增大；隨著入射波周期變長，運動響應峰值減小；隨著浮體相對寬度和相對吃水增大，垂蕩運動響應峰值減小。

(4)與現有方法不同的是，one-fluid模型是從剛體的分布拉格朗日乘子法出發，剛體被建模為流體相，并通過類似于空氣和水的廣義方程求解，剛體的空間速度是通過線性最小二乘投影來計算的。結果表明，浸入式one-fluid公式是適合于模擬波浪-浮體相互作用問題的，而且不需要明確地引入運動學和動力學邊界條件。