波浪作用下錨鏈系泊浮式防波堤動力響應(yīng)的數(shù)值模型研究

金哲飛，張金鳳, 2，張慶河，蔡 威，季超群

(1. 天津大學 水利工程安全與仿真國家重點實驗室，天津 300072； 2. 中國地震局地震工程綜合模擬與城鄉(xiāng)抗震韌性重點實驗室(天津大學)，天津 300350；3. 中交第三航務(wù)工程勘察設(shè)計院有限公司，上海 200032)

波浪在傳播過程中能量分布不均，其中98%的波能分布在距離水面三倍波高的水深范圍內(nèi)[1]。相比于傳統(tǒng)固定式防波堤，浮式防波堤能夠充分利用波能分布特性，在滿足工程消浪要求的同時減少對工程水域水沙交換的影響。由于浮式防波堤具備施工快捷、拆卸方便、對環(huán)境適應(yīng)性強且影響小等優(yōu)點，在近岸工程中得到了較為廣泛的應(yīng)用[2]。

浮式防波堤的動力響應(yīng)研究主要圍繞消浪、透浪性能展開。在理論研究方面，Bouwmeester和Van[3]與Drimer等[4]分別基于線性波理論和二維勢流理論對不同波況下單方箱浮式防波堤的透射系數(shù)計算進行了理論推導(dǎo)，Roul等[5]和何超勇等[6]在物理試驗基礎(chǔ)上對浮式防波堤透射系數(shù)公式進行了不同程度修正。物理試驗方面，Uzaki等[7]和Sannasiraj等[8]分別通過改進浮堤形式、錨泊方式提升浮式防波堤消浪性能。侯勇等[9-10]通過改變浮堤相對寬度、相對吃水深度、錨鏈拖地長度、系泊角度等物理試驗研究，分析總結(jié)了入射波周期、波高、水深、錨鏈系泊角度、拖地長度對浮式防波堤消浪性能的影響。

浮式防波堤運動過程受到波浪條件、浮堤結(jié)構(gòu)形式、系泊方式等多種因素影響，經(jīng)驗公式難以準確計算其水動力特性，而物理模型試驗操作復(fù)雜，成本較高。因此諸多學者提出不同數(shù)值算法對浮式防波堤水動力特性進行模擬。在波浪與浮式防波堤相互作用的數(shù)值模擬研究中，浮堤運動邊界的流固耦合處理方法及錨鏈模型的選取至關(guān)重要。其中浮堤運動邊界根據(jù)邊界條件處理方式的不同，主要分為動網(wǎng)格變形、浸入邊界、重疊網(wǎng)格法等。動網(wǎng)格方法適用于小幅度浮堤位移模擬，Ji等[11]使用分塊移動網(wǎng)格建立二維雙圓柱浮式防波動力響應(yīng)數(shù)值模型，并利用物理模型試驗驗證模型在計算浮堤動力響應(yīng)、錨鏈張力變化方面的準確性。王夢賢等[12]采用浸入邊界法處理波浪與起伏板防波堤的耦合作用問題，開展規(guī)則波與起伏水平板防波堤相互作用的數(shù)值模擬研究。浸入邊界法在三維問題中標志點的追蹤較為復(fù)雜，在三維模型中使用較為困難。重疊網(wǎng)格法最早由Benek等[13]提出，具備網(wǎng)格生成質(zhì)量高、適應(yīng)剛體大幅度運動的特點，廣泛應(yīng)用于復(fù)雜結(jié)構(gòu)物與流體相對運動的流固耦合問題求解。錨鏈張力模型建立主要基于動力學分析與靜力學分析兩方面展開。靜力學分析中錨鏈張力的數(shù)值求解主要有分段外推法等。當忽略錨鏈拉伸變形及流場作用力，靜力平衡方程退化為懸鏈線方程。繆泉明等[14]基于懸鏈線方程對于浮體動力響應(yīng)采用六自由度運動方程求解，研究極限海況下浮體的動力響應(yīng)及錨鏈受力，錨鏈張力峰值與實測結(jié)果存在較大差距。相比于靜力學分析，動力學分析考慮錨鏈運動過程中所受的動力荷載，在錨鏈劇烈運動狀態(tài)下纜繩錨鏈張力求解精度更高，錨鏈張力求解數(shù)值方法主要有集中質(zhì)量法、節(jié)點單元法等。Vanden[15]采用集中質(zhì)量法對錨鏈在動態(tài)響應(yīng)下的最大張力進行了細致研究。Huang[16]將模型拓展至三維，更真實地模擬實際工況，提升計算精度。

基于OpenFOAM中olaFlow求解器，采用基于sugger++[17-18]重疊網(wǎng)格方法模擬箱式浮式防波堤動力響應(yīng)，并在剛體六自由度運動求解過程中植入笛卡爾坐標系下MOODY錨鏈求解模塊[19]進行耦合求解，錨鏈張力分析基于動力學平衡方程展開，采用間斷有限元數(shù)值方法進行計算。建立鏈式結(jié)構(gòu)浮體運動的二維數(shù)值模型，并通過試驗對比驗證模型計算的準確性。

1 波浪作用下浮式防波堤運動數(shù)學模型

1.1 基本控制方程

采用雷諾時均的Navier-Stokes方程(RANS方程)作為基本控制方程，其連續(xù)性方程及動量方程分別為：

(1)

(2)

式中：U為流體速度；prgh為流體的修正壓強(動水壓強)，prgh=p-ρgh，p為總壓強，ρgh為靜水壓強；ρ為流體密度；μeff為考慮分子動力黏性和紊流作用的有效黏性系數(shù)，μeff=μl+μt，μl為層流黏性系數(shù)，μt為紊動黏性系數(shù)，由紊流模型聯(lián)合求解；C為流體表面張力系數(shù)，一般取C=0.07；κ為界面處曲率；α為VOF法引入的體積分數(shù)。

在求解RANS方程組時使用Menter[20]提出的SST k-ω兩方程紊流模型閉合，自由液面采用Hirt和Nichols[21]提出的VOF(volume of fuid)法進行追蹤捕捉。采用PIMPLE算法對壓力泊松方程進行迭代求解。

1.2 數(shù)值波浪水槽

利用OpenFOAM中olaFlow求解器建立數(shù)值波浪水槽，采用速度入口造波，出流邊界設(shè)置阻尼消波區(qū)，在該計算區(qū)域內(nèi)動量方程在右端加入式(3)所示阻尼項求解[22]。

Sf=ρUθ(x)

(3)

(4)

式中：Sf為阻尼源項；θ(x)為阻尼系數(shù)函數(shù)；x為網(wǎng)格單元坐標值；x0為阻尼消波區(qū)起始位置坐標值，Lf為阻尼消波區(qū)長度；c為經(jīng)驗系數(shù)，取c=20。

波浪水槽底部邊界為無滑移邊界條件，頂部邊界為壓力邊界。入流邊界設(shè)置為主動吸收式速度入口造波邊界，出流邊界為消波邊界[23]。側(cè)邊界為空邊界。為減少計算網(wǎng)格，在水槽y方向利用OpenFOAM中extrudeMesh進行網(wǎng)格拉伸處理成單層網(wǎng)格。二維浮堤運動共分為3個自由度，分別為垂蕩(Heave，沿z軸上下運動)，縱蕩(Surge，沿x軸前后運動)，縱搖(Pitch，繞y軸轉(zhuǎn)動)。浮堤運動計算采用六自由度運動方程求解，對其余三個自由度的運動進行限制，數(shù)值求解方法為Newmark算法。

1.3 重疊網(wǎng)格方法

為合理捕捉浮堤的運動過程，將浮式防波堤作為剛體處理，并采用重疊網(wǎng)格法進行模擬，如圖1所示。該方法由劃分固體的無變形貼體網(wǎng)格(虛線部分)、劃分流場的背景網(wǎng)格(實線部分)兩部分聯(lián)合求解完成，網(wǎng)格中插值節(jié)點流場信息計算方程：

(5)

式中：φL為邊界(邊緣)節(jié)點的某一流場信息；αk為其對應(yīng)的第k個插值貢獻節(jié)點的插值系數(shù)；φk為對應(yīng)第k個插值貢獻節(jié)點的某一流場信息。

重疊網(wǎng)格中的貼體網(wǎng)格與背景網(wǎng)格分別建立在兩套笛卡爾坐標系下，通過域連接信息(domain connectivity information，簡稱DCI)在背景網(wǎng)格與貼體網(wǎng)格間建立信息傳遞。計算時首先在貼體網(wǎng)格中尋找不參與計算的節(jié)點單元(洞單元，如圖1黑點)，并分別在背景網(wǎng)格中尋找邊緣節(jié)點(與洞單元相鄰的計算節(jié)點)，在貼體網(wǎng)格中尋找邊界節(jié)點(貼體網(wǎng)格中未被定義物理邊界條件的計算節(jié)點)；在邊界節(jié)點、邊緣節(jié)點尋找合適的插值貢獻節(jié)點并利用其位置關(guān)系求解插值系數(shù)；最后通過式(5)在兩套網(wǎng)格中交換流場數(shù)據(jù)，并根據(jù)上一時間步浮堤運動調(diào)整貼體網(wǎng)格區(qū)域，重新建立DCI連接。

圖1 背景網(wǎng)格與貼體網(wǎng)格組成的重疊網(wǎng)格系統(tǒng)

1.4 錨鏈求解模型

對于錨鏈受力計算部分采用節(jié)點單元法求解受力，在六自由度求解中耦合Palm開發(fā)的MOODY模塊[19]。錨鏈的受力計算考慮懸空段所受流體的各項作用力及拖地段受海床面的作用力，忽略錨鏈彎曲剛度，并假定錨鏈質(zhì)量不隨時間發(fā)生變化。錨鏈在笛卡爾坐標系中的運動方程為：

(6)

錨鏈在實際運動過程中應(yīng)力—應(yīng)變關(guān)系較為復(fù)雜。在受到波浪力沖擊時錨鏈易產(chǎn)生較大張力，致使錨鏈單元應(yīng)力—應(yīng)變產(chǎn)生非線性關(guān)系。計算中假設(shè)錨鏈應(yīng)力—應(yīng)變關(guān)系遵循胡克定律[24]，始終處于彈性變形區(qū)[19]：

T=EA0ε

(7)

f=fa+fb+fc+fd

(8)

式中：fa為附加質(zhì)量力，包括入射波產(chǎn)生的Froude-Krylov力；fb為浮力與重力的合力；fc為接觸力；fd為黏性拖曳力。

圖2 錨鏈錨鏈節(jié)點受力示意

附加質(zhì)量對錨鏈的動力學影響可分為：流體加速度引起的慣性力fa+和流體隨錨鏈加速變化所引起的慣性力fa-兩部分組成：

fa=fa+-fa-

(9)

根據(jù)Morison方程[25]，附加質(zhì)量力分別從法向、切向兩個方向計算：

(10)

(11)

式中：Ac為錨鏈截面面積；ρ為流體密度；U為流體速度；v為錨鏈速度；下標tan為速度的切向分量；CMn，CMt分別為法向、切向附加質(zhì)量力系數(shù)。

浮力項的計算表達式為：

(12)

式中：ρc為錨鏈密度；γ0為錨鏈單位長度質(zhì)量。

錨鏈尾端固定，其中拖地段與地面的接觸力基于雙線性彈簧和阻尼器進行計算，接觸力項計算表達式：

(13)

(14)

(15)

黏性拖拽力項根據(jù)Morison方程計算所得：

(16)

式中：CDt和CDn分別為切向和法向的拖拽力系數(shù)；v*為流體與錨鏈的相對速度；下標n為法向速度。

錨鏈求解采用間斷有限元方法，對式(6)展開并在錨鏈域內(nèi)進行空間離散。采用三階Runge-Kutta對離散后的式(6)進行時間步遞進，求出錨鏈單元節(jié)點的加速度、速度、位置信息，并更新單元節(jié)點張力。

1.5 浮堤運動求解過程

錨鏈力的求解通過修改OpenFOAM六自由度剛體運動求解順序，植入基于節(jié)點單元法的錨鏈求解模塊MOODY，并在六自由度計算過程中完成對錨鏈力的求解。植入MOODY模塊后重疊網(wǎng)格剛體運動計算具體流程：

1) 流場初始化，建立重疊網(wǎng)格域連接信息(DCI)；

2) 根據(jù)DCI，進行網(wǎng)格插值計算；

3) 依據(jù)全場流場信息，調(diào)用MOODY模塊求解錨鏈受力，并進行六自由度計算求解浮堤運動方程，對浮堤運動、錨鏈力耦合迭代求解，最后更新浮體周圍流場信息，更新浮堤、錨鏈信息；

4) 更新貼體網(wǎng)格信息；

5) 根據(jù)更新的網(wǎng)格、流場信息，重新建立DCI；

6) 根據(jù)邊界條件、上一時間步全場信息，PIMPLE迭代循環(huán)求解當前間步全場流場信息；

7) 進行新一輪迭代循環(huán)，直至殘差符合要求進入新一時間步進行循環(huán)求解，或達到計算時長停止計算。

2 數(shù)值模擬結(jié)果

2.1 模型設(shè)置

數(shù)值模擬計算模型參照單方箱錨鏈式浮式防波堤水動力學特性試驗，物理模型試驗布置及各項參數(shù)如圖3和表1所示,其余尺寸詳見侯勇物理模型試驗[9-10]。浮堤長(W)、高(H)分別為0.3 m、0.18 m，質(zhì)量為17.82 kg。錨鏈物理參數(shù)、波高測點布置均與侯勇試驗設(shè)置相同，錨鏈初始時刻拖地長度均為0。為避免造波邊界對浮堤的影響，同時減少造波過程中的二次反射影響，數(shù)值波浪水槽造波邊界與浮堤左側(cè)之間距離設(shè)為4倍入射波長，出流邊界附近設(shè)置2倍波長阻尼消波區(qū)，阻尼消波區(qū)距浮堤右側(cè)2倍波長，其設(shè)置如圖4所示。背景網(wǎng)格x方向長度取150分之波長，z方向取20分之波高。波浪水槽在y方向經(jīng)過extrudeMesh拉伸成寬度為0.44 m的單層網(wǎng)格。

圖3 物理模型試驗示意

表1 模型及錨鏈各參數(shù)設(shè)置

圖4 數(shù)值波浪水槽斷面示意

設(shè)置8組算例對二維錨鏈系泊浮式防波堤動力響應(yīng)數(shù)值型進行驗證，算例設(shè)置如表2所示。首先對數(shù)值波浪水槽進行空水槽造波校驗。在未放置防波堤的數(shù)值波浪水槽中，提取預(yù)留浮式防波堤位置處(距速度入口邊界4倍波長處)波高歷時變化作為對比，計算時長為30倍波周期，圖5以第二組波況為例顯示了模擬結(jié)果與二階Stokes波理論波面。由圖5可知，波高變化在第20T以后漸趨穩(wěn)定。表2列出了第20T～30T時刻波高均值與理論波高的對比情況，模擬結(jié)果與理論誤差不大于2%，可以認為浮堤附近波高值達到算例設(shè)置波高值。

數(shù)值計算采用6塊Intel Xeon E7 8870 v3處理器，每塊處理器擁有18個核心，每組算例計算時長為10.23小時,CPU時長為10.21小時。8組對比算例中，選取相對寬度W/L、波高h，作為模擬工況的控制變量。選取各組第20T～30T周期浮堤透射系數(shù)(Ct)、縱蕩(Surge)、垂蕩(Heave)、縱搖(Pitch)、迎浪測浮堤系泊處錨鏈張力(F1)、背浪側(cè)浮堤系播處錨鏈張力(F2)平均值與侯勇試驗值進行對比。

圖5 未放置防波堤的波浪水槽中x=4L位置處波高歷時與波浪理論值對比(組Ⅱ)

表2 波浪參數(shù)及數(shù)值模擬結(jié)果與理論解的比較

2.2 網(wǎng)格收斂性分析

在背景網(wǎng)格固定情況下，對組Ⅱ算例選取5種不同加密方式對貼體網(wǎng)格進行網(wǎng)格無關(guān)性分析，其中貼體網(wǎng)格最外層網(wǎng)格大小均為0.04 m×0.04 m，每級加密都對上一級加密中最靠近浮堤的四層網(wǎng)格尺寸進行減半。選取第20T～30T時刻垂蕩、縱蕩、縱搖峰值平均作為對比。從表3可以看出，當對貼體網(wǎng)格進行2級漸變加密時(圖6)，浮堤縱蕩、垂蕩、縱搖均值計算結(jié)果與3、4級加密計算結(jié)果相差較近達到收斂，且計算耗時較短。因此選用2級加密方案(圖6)對貼體網(wǎng)格進行局部加密。

表3 貼體網(wǎng)格收斂性研究

圖6 貼體網(wǎng)格加密方案

2.3 浮式防波堤透射系數(shù)

計算結(jié)果分別提取圖3中四個測點水面高度數(shù)據(jù)，并采用Goda兩點法[26]消除堤前波浪反射影響，并計算透射系數(shù)，各組透射系數(shù)如圖7所示,其中組Ⅳ兩者相差最小，僅為1.2%；組Ⅴ相差最大，為10.72%。圖7中各組透射系數(shù)均比侯勇物模試驗值偏小，最大偏差發(fā)生在組Ⅴ，達10.72%。造成這種誤差的原因可能與數(shù)值波浪水槽浮堤兩側(cè)無波浪透過，而物理模型試驗浮堤兩側(cè)存在一定透水現(xiàn)象有關(guān)。

圖7 各組透射系數(shù)(Ct)計算值與試驗值對比

2.4 浮式防波堤動力響應(yīng)

數(shù)值模擬試驗中浮堤運動共三個自由度，分別為縱蕩、垂蕩和縱搖。圖8以組Ⅱ為例顯示了縱蕩、垂蕩和縱搖模擬計算結(jié)果。取第20T～30T周期浮塊運動幅值平均值，并對三個自由度的振幅分別做無因次化處理，其運動特性通過響應(yīng)振幅算子RAO加以描述：

(10)

式中：縱蕩、垂蕩和縱搖響應(yīng)振幅算子RAO中的響應(yīng)幅值分別取浮堤20T～30T時刻縱蕩、垂蕩和縱搖運動幅值均值，相應(yīng)的響應(yīng)算子單位分別為cm/cm，cm/cm，°/cm。

圖8 組Ⅱ模擬歷時

圖9顯示了組Ⅱ第28T各時刻浮堤運動響應(yīng)。當波谷到達浮堤位置(t=0T)，浮堤垂蕩達到谷值。隨波浪傳播垂蕩、縱搖逐漸恢復(fù)至平衡位置，隨后縱蕩達到谷值，縱搖達到峰值(t=5/20T)。當波峰到達浮堤位置(t=10/20T)，浮堤縱蕩、垂蕩均達到峰值，縱搖接近平衡位置。波浪在浮堤迎浪側(cè)形成反射造成浮堤縱搖出現(xiàn)次峰。隨后浮堤垂蕩再次回歸平衡位置，縱搖到達谷值(t=15/20T)。當浮堤垂蕩再次經(jīng)過平衡位置(t=15/20T)，縱搖逐漸回歸平衡位置，垂蕩逐漸到達谷值(t=20/20T)。由圖8可知浮堤的縱蕩和垂蕩基本與波面同相位，縱搖與波面相位相差約T/4。

圖9 組Ⅱ第28T(35.84～37.12 s)浮堤運動響應(yīng)

圖10顯示了不同工況下浮堤運動響應(yīng)振幅算子與侯勇試驗值的對比情況，兩者偏差最大者為組Ⅶ垂蕩響應(yīng)振幅算子，為10.44%；組Ⅳ縱蕩響應(yīng)振幅算子模擬與試驗值偏差最小，為0.6%?？傮w而言，模擬結(jié)果與試驗結(jié)果吻合較好。

圖10 響應(yīng)振幅算子模擬值與試驗值對比

2.5 錨鏈張力

模擬結(jié)果取錨鏈-浮堤連接端節(jié)點張力值作與試驗結(jié)果進行對比。圖11為組Ⅱ錨鏈張力變化歷時曲線。浮體動力響應(yīng)呈穩(wěn)定周期性變化后(20T)，迎浪側(cè)錨鏈張力(F1)第20T～30T周期平均峰值為40.12 N，背浪側(cè)錨鏈張力(F2)平均值為29.85 N。迎浪側(cè)與背浪側(cè)錨鏈張力峰值交替出現(xiàn)，且迎浪側(cè)錨鏈張力峰值更大。

圖11 組Ⅱ錨鏈張力歷時曲線

圖12為浮堤錨固處錨鏈張力振幅算子與試驗值對比(侯勇試驗值無組Ⅳ，Ⅶ，Ⅷ)，誤差最大為21.8%(組Ⅰ背浪側(cè))，最小為2.8%(組Ⅱ迎浪測)。錨鏈背浪側(cè)張力模擬和試驗值的偏差相對較大，其原因可能與浮堤運動響應(yīng)過程中背浪側(cè)錨鏈松弛時間更長(圖9)，拖地段錨鏈與海床的接觸力模擬誤差較大，且錨鏈在張緊瞬時受沖擊峰值存在一定震蕩現(xiàn)象相關(guān)。模擬中當錨鏈拖地段運動速度達到0.1 m/s時到達靜摩擦最大值并近似認為進入滑動摩擦階段，滑動摩擦力不變。實際錨鏈運動時拖地段地面接觸力變化比模擬設(shè)置條件更加復(fù)雜。此外，模擬中對于錨鏈的應(yīng)力—應(yīng)變本構(gòu)模型采用了線性化假設(shè)，在錨鏈承受劇烈沖擊荷載時，應(yīng)力—應(yīng)變關(guān)系與實際情況存在一定偏差。

圖12 錨鏈張力計算值與試驗對比

3 結(jié) 語

基于OpenFOAM中olaFlow求解器，利用重疊網(wǎng)格在計算過程中適應(yīng)剛體大幅運動的優(yōu)勢，采用重疊網(wǎng)格方法模擬浮堤動力響應(yīng)，并在剛體運動求解過程中植入笛卡爾坐標系下MOODY錨鏈求解模塊進行耦合求解，建立錨鏈系泊浮式防波堤動力響應(yīng)的二維數(shù)值模型。浮式防波堤錨鏈尾端在海床中固定，因此在數(shù)值模擬中可簡化成MOODY錨鏈模型進行錨鏈張力的求解計算。數(shù)值模型通過侯勇物理模型試驗進行對比，并分別驗證透射系數(shù)、浮堤運動響應(yīng)、錨鏈張力峰值模擬結(jié)果。

在透射系數(shù)對比中，由于數(shù)值波浪水槽在y方向?qū)W(wǎng)格進行了拉伸處理，相比物模試驗不存在入射波從浮堤兩側(cè)透過的情況，因此透射系數(shù)模擬值整體偏小(誤差為1.2%～10.72%)；在浮堤運動響應(yīng)方面，垂蕩、縱蕩、縱搖運動響應(yīng)模擬誤差在0.6%～10.44%；由于錨鏈拖地段與海床接觸力基于雙線性彈簧和阻尼器進行求解計算，實際錨鏈受力更加復(fù)雜，因此背浪側(cè)錨鏈張力峰值模擬誤差部分達到21.8%(組Ⅰ背浪側(cè)錨鏈張力)。此外，模擬中錨鏈應(yīng)力-應(yīng)變模型采用了線性化假設(shè)，在浮堤動力響應(yīng)劇烈、錨鏈受到瞬時沖擊荷載較大時，應(yīng)力-應(yīng)變關(guān)系存在一定誤差。

通過數(shù)值模型模擬結(jié)果，該模型在計算波浪作用下錨鏈系泊浮式防波堤動力響應(yīng)中具有比較高的模擬精度，可用于錨鏈材料應(yīng)力-應(yīng)變模型改進，浮式防波堤三維數(shù)學模型開發(fā)、浮堤消浪性能評價與結(jié)構(gòu)優(yōu)化計算。