羅天，萬德成

1上海交通大學船舶海洋與建筑工程學院，上海200240

2上海交通大學海洋工程國家重點實驗室，上海200240

3高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海200240

基于CFD的船舶橫搖數(shù)值模擬與粘性效應(yīng)分析

羅天1，2，3，萬德成1，2，3

1上海交通大學船舶海洋與建筑工程學院，上海200240

2上海交通大學海洋工程國家重點實驗室，上海200240

3高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海200240

［目的］船舶橫搖的準確預報對于船舶耐波性、穩(wěn)性以及操縱性的研究具有十分重要的意義，但船舶橫搖運動受流場粘性效應(yīng)的影響很大，計算中存在較多非線性因素，因而并不適于常用于研究船舶運動的傳統(tǒng)勢流理論。為解決這一問題，［方法］采用基于OpenFOAM軟件開發(fā)的naoe-FOAM-SJTU求解器，分別通過歐拉方法和RANS方法對S60船模典型二維橫剖面的強迫橫搖進行數(shù)值模擬，同時，模擬分析DTMB 5512不同三維船模的自由與強迫橫搖運動。［結(jié)果］在成功將橫搖阻尼力矩的不同成分分別計算出來后發(fā)現(xiàn)，其中漩渦阻尼所占的比例最大，摩擦阻尼所占的比例最小，而舭龍骨在一定的橫搖角度范圍內(nèi)則減小了橫搖阻尼力矩。［結(jié)論］該結(jié)果揭示了橫搖參數(shù)對船舶橫搖運動及橫搖阻尼力矩的影響，對準確預報船舶的橫搖運動具有重要意義。

船舶；橫搖運動；阻尼系數(shù)；阻尼力矩；RANS方法；歐拉方法

0 引 言

在現(xiàn)代船舶設(shè)計中，船舶橫搖運動的準確預報十分重要，然而流場的粘性效應(yīng)對船舶橫搖運動影響較大，使得船舶橫搖運動計算中含有很多非線性因素，因此傳統(tǒng)的勢流理論無法準確預報船舶橫搖運動。

船舶橫搖阻尼系數(shù)是決定船舶橫搖預報準確性的關(guān)鍵因素。船舶橫搖阻尼系數(shù)通?？梢苑譃?個部分［1］：摩擦阻尼系數(shù)、波浪阻尼系數(shù)、漩渦阻尼系數(shù)、升力阻尼系數(shù)和附體阻尼系數(shù)。這種分類便于逐個研究外界因素對阻尼系數(shù)不同成分的影響，但在該分類下，5種阻尼系數(shù)之間是假定線性相關(guān)的。對于這5個部分，Ikeda等［2-3］于1977～1978年根據(jù)船模試驗提出了相對應(yīng)的經(jīng)驗公式，但船模試驗的效率不高，經(jīng)驗公式的適用范圍有限，因此這2種方法并不能大規(guī)模地用于船模橫搖研究。

如今，隨著計算機技術(shù)的飛速發(fā)展，計算流體力學［4］（Computational Fluid Dynamics，CFD）被廣泛應(yīng)用于船舶相關(guān)研究。在船舶橫搖運動模擬分析方面，學者們相繼進行了研究。Luca與Stefano［5］利用CFD軟件針對相似橫剖面的單體船、雙體船和小水線面船進行了橫搖運動模擬；Miller等［6］應(yīng)用RANS（Reynolds Averaged Navier Stokes）方法對存在舭龍骨的三維圓柱體進行了預報；Wilson等［7］應(yīng)用CFDSHIP-IOWA求解器預報了船舶在不同船速下的自由橫搖以及強迫橫搖；Yang等［8］針對S60船模的自由與強迫橫搖運動進行了模擬并計算了其阻尼系數(shù)。但以上研究都只得到了總的阻尼系數(shù)，或是僅將摩擦阻尼系數(shù)分離出來進行了分析。

本文將應(yīng)用課題組自行開發(fā)的naoe-FOAMSJTU求解器［9］模擬船模運動與流場的變化。該求解器是基于開源軟件OpenFOAM［10］（Open Field Operation and Manipulation）而開發(fā)，以RANS方程為控制方程，采用有限體積法（Finite Volume Method，F(xiàn)VM）離散控制方程，流體體積法（Volume of Fluid，VOF）技術(shù)捕捉自由液面，用SSTk-ω湍流模型封閉方程，應(yīng)用壓力的隱式算子分割算法（Pressure-Implicit Split-Operator，PISO）［11］耦合求解速度和壓力。而對于船體的運動，則采用動態(tài)網(wǎng)格技術(shù)以及六自由度運動求解器來模擬。naoe-FOAM-SJTU求解器已應(yīng)用于船舶工程、近岸工程、海洋工程等各類復雜流動問題的數(shù)值模擬，包括各類水波的數(shù)值造波與消波、船槳舵匹配自航推進與操縱運動、節(jié)能裝置水動力性能、船舶運動與粘性阻尼、上浪砰擊、浮式平臺運動與氣隙、液化天然氣（LNG）液艙晃蕩、Spar平臺渦激運動等復雜問題，相關(guān)計算結(jié)果也得到了很多試驗的證實。運用該求解器，吳建威等［12］成功預報了螺旋槳敞水性能，Shen等［13］得到了一種非規(guī)則波的生成方法，李鵬飛等［14］基于致動線模型數(shù)值模擬了風力機尾流場，尹崇宏等［15］模擬了實尺度VLCC船的阻力問題，Cao等［16］數(shù)值模擬了多方向海洋波問題，Zha等［17］數(shù)值模擬了雙體船興波阻力問題。Shen等［18］還數(shù)值模擬了船舶自航運動和船—槳—舵的相互作用問題。以上數(shù)值模擬結(jié)果均表明naoe-FOAM-SJTU求解器在解決船舶與海洋工程相關(guān)問題方面是可靠的。

本文將分別采用歐拉方法與RANS方法模擬二維船模橫剖面的強迫橫搖運動，并從總阻尼系數(shù)中分離出波浪阻尼系數(shù)、摩擦阻尼系數(shù)以及漩渦阻尼系數(shù)。對于三維船模，擬通過比較帶舭龍骨的船體與裸船體的橫搖運動，從總阻尼系數(shù)中分離出摩擦阻尼系數(shù)與舭龍骨阻尼系數(shù)。由于本文中所涉及到的橫搖運動并未考慮船模橫蕩運動的影響，因此升力阻尼系數(shù)不在分析范圍內(nèi)。

1 數(shù)值方法

1.1 RANS方法

對于naoe-FOAM-SJTU求解器，以RANS方程為控制方程：

式中：x為網(wǎng)格節(jié)點的位置；t為時間；U為速度場；Ug為網(wǎng)格移動速度；pd=p-ρg·x，為動壓力；ρ為液體或者氣體的密度；g為重力加速度向量；fσ為表面張力項；fs為源項；μeff=ρ(ν-νt)，為有效動力粘度，其中ν為運動粘度，νt為渦黏度。

湍流模型SSTk-ω用來封閉該控制方程，邊界可壓縮的VOF方法用來追蹤和捕捉自由液面的變化，其輸運方程定義如下：

式中：Ur為用于壓縮界面的速度場；α為積分數(shù)，代表網(wǎng)格單元中流體所占的體積比，其值介于0和1之間：

控制方程采用FVM法進行離散，即將計算域離散成一系列小單元，計算流場信息存儲在網(wǎng)格單元中心，最終根據(jù)單元中心的值插值得到單元面的值。對于離散后得到的壓力速度耦合方程，則采用PISO算法進行迭代求解，其主要由分預測—修正—再修正三步組成。具體步驟與算法詳見文獻［9］。

1.2 歐拉方法

本文采用歐拉方法對二維S60船模橫剖面的強迫橫搖進行計算。由于摩擦阻尼系數(shù)和漩渦阻尼系數(shù)均與流場的粘性相關(guān)，同時，本文假定波浪阻尼系數(shù)與流場粘性無關(guān)，因此由歐拉方法最終得到的為波浪阻尼系數(shù)。

歐拉方法的實現(xiàn)同樣是基于naoe-FOAMSJTU求解器?？蓪⑶蠼馄髦姓承皂椀臄?shù)值設(shè)置成0，并改變流場的邊界條件來實現(xiàn)歐拉方法，具體同第1.1節(jié)。

1.3 靜水阻力系數(shù)

為了驗證歐拉方法的可行性，同時驗證波浪與粘性無關(guān)，本文運用歐拉方法以及RANS方法對KCS船模在靜水中的運動進行模擬，并求解靜水阻力系數(shù)。

采用RANS方法可以直接計算得到粘性流場中船舶的總阻力，從而求得阻力系數(shù)

式中：F為船舶總阻力；S為船舶濕表面積；V為船速。

而采用歐拉方法則僅能計算出興波阻力系數(shù)Cw，摩擦阻尼系數(shù)Cf采用1957ITTC摩擦公式：

式中，雷諾數(shù)Re=VL/ν，其中L為船長。最終的阻力系數(shù)C=Cf+Cw。

1.4 橫搖阻尼系數(shù)

本文中，為了與試驗數(shù)據(jù)進行對比驗證，需要計算相關(guān)的橫搖阻尼系數(shù)。強迫橫搖與自由橫搖阻尼系數(shù)的計算方法并不一致，具體如下。

1.4.1 強迫橫搖

對于強迫橫搖，橫搖角度隨時間的變化可以

表示為如下方程：

式中：θ0為初始橫搖角；ω為橫搖頻率。

橫搖阻尼力矩可以表示為

式中，Beq為阻尼系數(shù)。假設(shè)回復力系數(shù)C(θ，t)是關(guān)于θ的多項式，則

當θ=0時，可以推導得到

式中：M0為船舶橫搖角度為0°時船舶所受到的阻尼力矩；B為船寬。

由以上方程可以看出，對于某一種確定的工況，阻尼系數(shù)與橫搖阻尼力矩相關(guān)，因此，橫搖阻尼力矩相應(yīng)地也可以分為5個部分。為直觀分析，本文將優(yōu)先從橫搖阻尼力矩曲線以及最大力矩來進行分析。

本文采用naoe-FOAM-SJTU求解器求得的總阻尼力矩可以分為剪切和法向2個部分。剪切部分是由船體周圍的剪切力所產(chǎn)生，即摩擦力產(chǎn)生的阻尼力矩，相對應(yīng)的為摩擦阻尼系數(shù)。而法向部分則是由船體周圍的不平衡壓力所造成，對應(yīng)于波浪阻尼系數(shù)和漩渦阻尼系數(shù)，其中假定波浪阻尼系數(shù)與粘性無關(guān)，是采用歐拉方法計算得到。舭龍骨的存在增加了船體摩擦，同時也改變了船體周圍的渦量分布和渦量大小，因此對于剪切與法向2個部分均有影響。4個部分（漩渦、波浪、摩擦、附體）的橫搖阻尼力矩的關(guān)系如下：

式中：Mp和Mf為采用RANS方法求得的阻尼力矩法向部分與剪切部分，其中Mp是由漩渦、波浪與舭龍骨共同影響產(chǎn)生，Mf是由摩擦部分產(chǎn)生；Mpw為采用歐拉方法求得的阻尼力矩法向部分，是由波浪影響產(chǎn)生的；Mpe為由漩渦影響產(chǎn)生的橫搖力矩；Mpb為由舭龍骨單獨影響產(chǎn)生的橫搖力矩的法向部分；Mff為由裸船體的摩擦產(chǎn)生的橫搖剪切力矩；Mfb為由舭龍骨單獨的摩擦產(chǎn)生的橫搖剪切力矩。

對于自由橫搖運動，除了需要通過橫搖角度的時歷曲線來估算橫搖衰減周期外，還需通過該曲線來計算船舶自由橫搖的阻尼系數(shù) μ??［19］。圖1所示為典型的船舶自由衰減橫搖曲線。

根據(jù)船舶耐波性理論［20］，由該曲線可以得到相鄰兩振幅的差Δφ=φk-φk+1以及均值φm=(φk+φk+1)/2，以 φm為橫軸、Δφ為縱軸作消滅曲線如圖2所示。

則船舶自由橫搖阻尼系數(shù)可通過下式求得：

2 計算模型和網(wǎng)格

本文一共使用了3種船型：KCS船型、S60船型以及DTMB 5512船型。計算網(wǎng)格由OpenFOAM自帶的網(wǎng)格生成工具SnappyHexMesh生成。為了提高計算精度，捕捉自由面興波以及船體附近的渦量變化情況，對自由面和船體壁面附近處的網(wǎng)格進行了適當加密。除此外，為了實現(xiàn)船模的橫搖運動，還在計算中引入了動態(tài)網(wǎng)格（Dynamic mesh）技術(shù)。

2.1 KCS船型

在通過靜水阻力驗證歐拉方法的可行性時，為便于與試驗數(shù)據(jù)［21］進行對比，采用了與試驗相同的船型——KCS船型，其船模參數(shù)如表1所示。

KCS船的船體視圖如圖3所示。

由于船體的對稱性，僅在船模中縱剖面的右舷側(cè)劃分了網(wǎng)格，計算域尺寸為：-1.0L≤x≤4.0L，0≤y≤1.5L，-1.0L≤z≤1.0L。網(wǎng)格總量為108×104，船體艏、艉網(wǎng)格以及背景網(wǎng)格如圖4所示。

（1）通過對國內(nèi)外現(xiàn)有軟件的應(yīng)用分析調(diào)研，確定基于BIM技術(shù)的綠色建筑規(guī)劃設(shè)計軟件架構(gòu)和數(shù)據(jù)結(jié)構(gòu)，建立基于BIM技術(shù)的三維建筑模型

2.2 S60船型

對于本文中討論的典型的二維船型橫剖面，為便于與試驗數(shù)據(jù)［22］進行對比，采用了與試驗相同的船型及剖面——S60船型的典型橫剖面，船型主要參數(shù)如表2所示。

對于該二維橫剖面的橫搖模擬，橫搖軸為x軸，橫搖的軸心O為船舶橫剖面中心線與船模橫剖面靜止時水線的交點。計算域尺寸為：-2.5B≤y≤2.5B，-2B≤z≤2B。另外，還對該橫搖算例進行了網(wǎng)格收斂性分析，網(wǎng)格量分別為4×104，7×104和14×104，計算網(wǎng)格如圖5所示。

2.3 DTMB 5512船型

對于本文中討論的典型三維船型，為便于與試驗數(shù)據(jù)進行對比，采用了與試驗相同的船型——DTMB 5512船型，船型主要參數(shù)如表3所示。

該船體的裸船模型以及帶舭龍骨的船模如圖6所示。

對于橫搖運動，該船型的橫搖中心與重心重合，橫搖軸為中縱剖面上通過重心且平行于水線面的直線。計算域尺寸為：-1.0L≤x≤3.0L，-2.0L≤y≤2.0L，-1.0L≤z≤0.5L。為了驗證網(wǎng)格的可靠性，該算例也做了網(wǎng)格收斂性驗證，網(wǎng)格量分別為97.8×104，232×104和646×104，其中，232×104的計算網(wǎng)格如圖7所示。

3 驗證與計算結(jié)果

3.1 S60典型二維橫剖面橫搖計算

本文結(jié)合歐拉以及RANS這2種方法對S60二維橫剖面的強迫橫搖運動進行了模擬計算，最終將橫搖阻尼力矩分為了波浪、摩擦和漩渦3個部分。

3.1.1 歐拉方法可行性驗證

為了驗證歐拉方法的可行性以及其對船體周圍波浪預報的準確性，本文使用naoe-FOAMSJTU求解器，采用RANS和歐拉2種方法對不同傅汝德數(shù)Fr下的船舶靜水阻力進行了計算并與試驗數(shù)據(jù)進行了對比，計算結(jié)果如表4所示。

從表4的結(jié)果可以看出，由該求解器應(yīng)用歐拉方法得到的結(jié)果與試驗結(jié)果相近，即歐拉方法計算得到的興波阻尼部分是可信的，可以由該方法計算橫搖阻尼力矩中的波浪阻尼力矩。

當Fr=0.27時，由2種方法計算得到的波形比較如圖8所示。圖中，上半部分為歐拉計算結(jié)果，下半部分為RANS計算結(jié)果。

由上述阻力系數(shù)比較可以看出，采用歐拉方法和RANS方法所得結(jié)果與試驗結(jié)果的誤差在允許的范圍之內(nèi)，說明歐拉方法對于非粘性流場的計算結(jié)果是可靠的。從采用歐拉方法與RANS方法模擬得到的穩(wěn)定波形圖上來看，兩者在船體附近的波形圖基本相同，因此可以推斷流場的粘性效應(yīng)對于船體附近的興波并沒有太大的影響，繼而可以推論，在橫搖模擬中，采用歐拉方法計算波浪阻尼部分的做法是可行的。

3.1.2 RANS方法可行性驗證

為了與Ikeda等［2-3］的試驗數(shù)據(jù)進行對比，本文采用RANS方法在相同的橫搖工況下對S60典型二維橫剖面進行了強迫橫搖計算，橫搖工況如表5所示。

首先選用0.17 rad進行收斂性驗證，結(jié)果如表6所示。

由表6中數(shù)據(jù)可知，網(wǎng)格2與網(wǎng)格3的誤差均在允許范圍內(nèi)，但是考慮到效率問題，在后面的計算中，將選擇網(wǎng)格2進行計算。

為了進一步驗證網(wǎng)格2計算結(jié)果的可靠性，使用網(wǎng)格2對另外3種橫搖角度下的強迫橫搖進行了計算模擬，計算結(jié)果如表7所示。

由表7中的數(shù)值結(jié)果可以看出，誤差均在允許范圍內(nèi)，可以證實采用RANS方法以及網(wǎng)格2計算該二維船模橫剖面強迫橫搖運動所得結(jié)果是可以信賴的。

3.1.3 S60典型二維橫剖面強迫橫搖計算結(jié)果

為了單獨分析二維船模橫剖面強迫運動時橫搖阻尼系數(shù)的不同成分，本文運用歐拉方法以及RANS方法對于同一種橫搖工況下的船舶橫搖進行了模擬，具體如表8所示。

在該工況下，二維船模橫剖面強迫橫搖的力矩時歷曲線如圖9所示。

圖9中，黑色曲線為采用歐拉方法算得的船舶所受力矩的垂向部分 Mpw，紅色曲線為采用RANS方法算得的船舶所受摩擦力矩Mff，藍色曲線為采用RANS方法求得的力矩垂向部分與Mpw之差，即Mpe。

粘性流場中，船舶強迫橫搖時船體橫剖面周圍流場的渦量分布如圖10所示。

圖10為第（n-1）個周期內(nèi)4張典型時刻的渦量分布圖，其中T表示強迫橫搖周期，t表示該圖所對應(yīng)的時刻。船舶強迫橫搖時粘性流場與非粘性流場中的渦量分布圖如圖11所示。

由圖9可以看出，在強迫橫搖阻尼力矩中，漩渦影響部分相較于其他部分占比要大得多，因為在船體橫搖過程中周圍會產(chǎn)生很多漩渦，造成船體周圍壓力的不均衡，相比于摩擦以及由波浪引起的力矩，這部分產(chǎn)生的力矩影響更大。同時，這種結(jié)果與Ikeda等［2-3］的試驗結(jié)果和Korpus等［23］的數(shù)值計算結(jié)果均很相似。另外值得注意的是，波浪部分產(chǎn)生的力矩與其他2個部分力矩的方向是相反的，其原因分析如圖12所示。

由圖12可知，波浪對船舶形成的力矩使船趨向于傾覆方向，因此與另外2個部分力矩的方向是相反的，但因其所占總力矩的比例較小，故總的橫搖力矩還是與漩渦產(chǎn)生的橫搖力矩方向相同。

另外，從力矩時歷曲線圖（圖9）中可以看出，漩渦阻尼部分和摩擦阻尼部分達到峰值與零點是有時間差的，摩擦阻尼部分出現(xiàn)峰值的時間比漩渦阻尼部分要遲。船模所受到的漩渦阻尼部分由二維船模橫剖面橫搖時流場的具體狀態(tài)所決定，由上面的數(shù)據(jù)可以看出，峰值出現(xiàn)在橫搖角度最大的時刻，對應(yīng)于計算結(jié)果的渦量圖，發(fā)現(xiàn)此時船模橫剖面周圍的漩渦最多，尤其是舭部，造成了船舶周圍壓力的極度不平衡，因此，力矩的漩渦阻尼部分達到峰值。而摩擦阻尼部分是由船模橫剖面臨近的液體與船模表面的相對速度所決定，當船模橫剖面處于橫搖角度最大時，船模表面的速度為0，而臨近的流體由于慣性仍有速度，之后，船模向反方向運動，使得相對速度繼續(xù)增大，即力矩的摩擦阻尼部分會繼續(xù)增大，因此摩擦阻尼部分峰值出現(xiàn)的時間會比漩渦阻尼部分晚。

分析船舶橫搖渦量圖（圖10）可知，漩渦產(chǎn)生于二維船模橫剖面的壁上，在船模橫搖角度從最大開始逐漸變小時，由于相對速度在船模周圍開始形成一層很薄的、方向一定的漩渦，因此隨著船舶逐漸左傾（或右傾），在船模橫搖角度到另一個方向達到最大時，該方向上的漩渦已全部產(chǎn)生并附在船體表面形成較厚的一層。之后，船模反方向轉(zhuǎn)動，立刻會在船模周圍形成反方向的漩渦，中和一部分外層不同方向的漩渦，剩下沒有被中和的部分則沿著船體底部和船側(cè)方向分別排向船底與自由液面，在船模周圍形成4處漩渦中心。這4處漩渦的方向由最初船體轉(zhuǎn)動的方向所決定，當船體初始轉(zhuǎn)動方向變動時，漩渦的方向?qū)⒆優(yōu)榉聪颉?/p>

根據(jù)歐拉方法和RANS方法計算所得的無粘性與粘性流場中的壓力與渦量分布圖（圖11）可以明顯看出，當無粘性時，流場中并未產(chǎn)生明顯的渦量，因此可以證實阻尼系數(shù)的漩渦部分僅由粘性影響產(chǎn)生，即漩渦阻尼系數(shù)是根據(jù)法向力矩的粘性部分而求得。

3.2 DTMB 5512三維船模橫搖計算

本文采用RANS方法對DTMB 5512三維船型的橫搖運動進行了模擬計算，最終將摩擦阻尼部分以及舭龍骨阻尼部分從總阻尼中分離了出來。

3.2.1 三維船模自由橫搖計算驗證

對于DTMB 5512三維船模，僅有在穩(wěn)定航速下的自由橫搖試驗結(jié)果，因此本節(jié)采用同樣的計算工況來模擬該船模的自由橫搖運動以驗證RANS方法和網(wǎng)格的可行性。計算工況具體如表9所示。收斂性驗證結(jié)果如表10所示。

從表10中的數(shù)據(jù)可知，網(wǎng)格2與網(wǎng)格3的誤差均在允許范圍內(nèi)，但是考慮到效率問題，在后面的計算中將選擇網(wǎng)格2進行計算。

為了進一步驗證使用該方法以及該網(wǎng)格所得結(jié)果的可靠性，使用網(wǎng)格2，根據(jù)試驗數(shù)據(jù)采用了2種來流速度進行計算驗證，2種工況具體如表11所示。

計算結(jié)果如圖13和表12所示。由上面的結(jié)果可以看出，本文計算與試驗間的誤差均在6%以內(nèi)，可以證明采用RANS方法以及網(wǎng)格2計算該三維船模橫剖面強迫橫搖運動結(jié)果可信。此外，從圖13中可以發(fā)現(xiàn)，數(shù)值模擬結(jié)果與試驗結(jié)果橫搖周期的相對誤差一定，但絕對誤差值會隨時間累積，在自由橫搖的橫搖角度時歷曲線上，就表現(xiàn)為數(shù)值模擬結(jié)果與試驗結(jié)果間的誤差是隨時間逐步增大的。

3.2.2 DTMB 5512三維船模強迫橫搖計算結(jié)果

為了單獨分析三維船模強迫運動時橫搖阻尼系數(shù)的不同成分，采用RANS方法對同一工況下帶舭龍骨的船模與裸船模進行了模擬計算，工況具體如表13所示。

在該工況下，強迫橫搖的力矩時歷曲線如圖14所示。

圖14中，黑色線為帶舭龍骨船模所受的總力矩M，綠色曲線為裸船模所受力矩總力矩的垂向部分 Mpw+Mpe，藍色曲線為裸船模所受的摩擦力矩 Mff，紅色曲線為受舭龍骨影響的力矩的垂向部分Mpb，粉色曲線為受舭龍骨影響的力矩的切向部分Mfb。

粘性流場中船舶強迫橫搖時船體橫剖面周圍流場的渦量分布如圖15所示。

由圖15中結(jié)果可以看出，在三維船模的強迫橫搖中，與二維類似，仍然是Mpw+Mpe占總橫搖阻尼力矩的大部分，而裸船模的摩擦力矩和帶舭龍骨船模的摩擦力矩所占的比例則非常小。此外值得注意的是，對于帶舭龍骨的船模，計算得到的橫搖總阻尼力矩減小了，其具體原因如圖16所示。

由舭龍骨造成的阻尼力矩變化的原因來源于2個方面：一是舭龍骨增加了船模表面積從而增加了摩擦阻尼；二是舭龍骨在船舶強迫橫搖中對流場中渦量的數(shù)值以及分布有一定的影響從而改變了船舶受到的橫搖阻尼力矩。其中，摩擦阻尼占據(jù)很小的一部分，在此階段的分析中可以忽略不計。第2個因素首先考慮的是由舭龍骨附近產(chǎn)生的壓力不平衡所造成的阻尼力矩。如圖16所示，在船模轉(zhuǎn)到最高點往平衡位置轉(zhuǎn)的時候，總阻尼力矩是有利于船模轉(zhuǎn)回平衡位置的，但此時由舭龍骨受到的流場的力造成的橫搖阻尼力矩則是阻止船模轉(zhuǎn)回平衡位置的。其次，由圖15所示的渦量分布圖可以看出，在有舭龍骨的情況下，在舭龍骨的范圍內(nèi)船體舭部形成了2個漩渦中心，改變了船模周圍流場渦量的數(shù)值與分布。在這2種因素的綜合影響下，流場對船模造成的 Mp在一定的橫搖角度范圍內(nèi)減小了，其中也包括Mp的最大值。

與二維船模橫剖面強迫橫搖類似，三維船模的摩擦阻尼力矩時歷曲線與阻尼力矩垂向部分的時歷曲線具有一定的相位差。Mp與船體周圍的渦量有關(guān)，在船模橫搖達到最大角度時，船模附近的渦量分布最不均勻，渦量數(shù)值最大，因此Mp也達到最大值。當船模從橫剖面處于橫搖角度最大處向反方向運動時，臨近流體與船模表面的相對速度會繼續(xù)增大，即力矩的剪切力部分會繼續(xù)增大，因此兩者的時歷曲線會有一定的相位差。

由圖15所示的渦量圖還可看出，船模附近漩渦以及漩渦中心形成的原理也與二維船模橫剖面類似。漩渦在船模橫搖時從船模表面開始形成，隨著船模的運動，部分與新形成的漩渦相互中和消散，部分受船模運動的影響而向外移動。由于算例中的強迫橫搖角度偏小，裸船模最終僅形成了2個漩渦中心，而帶舭龍骨的船模在舭龍骨附近也會形成渦中心。圖15還展示了船模橫搖時自由表面的變化，可以看出，帶舭龍骨船模與裸船模形成的行波基本相似。

4 結(jié) 語

本文采用RANS方法對S60船模二維橫剖面的橫搖進行模擬計算，驗證了網(wǎng)格以及計算方法的可行性。在此基礎(chǔ)上，結(jié)合歐拉方法下同一模型與工況的計算結(jié)果，將船舶阻尼系數(shù)分為了摩擦、波浪和漩渦3個部分進行分析，其中摩擦阻尼部分在總阻尼里占的比重最少，而漩渦阻尼部分占的比重最多，兩者的相關(guān)力矩曲線具有一定的相位差。另外，波浪阻尼力矩的方向與摩擦和漩渦這2種阻尼力矩的方向是相反的。

本文還采用RANS方法對DTMB 5512三維船模進行了自由衰減的橫搖運動模擬，驗證了網(wǎng)格以及計算方法的可行性。在此基礎(chǔ)上，對該三維裸船模與帶舭龍骨船體進行了強迫橫搖運動的模擬計算，將摩擦阻尼部分與舭龍骨阻尼部分從總阻尼中分離出來進行了分析，與二維結(jié)果類似，摩擦阻尼部分占比非常小，波浪與漩渦阻尼的影響較大。而舭龍骨則在一定的橫搖角度范圍內(nèi)減小了總阻尼力矩，即阻尼力矩的舭龍骨影響部分曲線與總阻尼力矩曲線具有一定的相位差。

在以后的研究中，還將分析舭龍骨在船模自由橫搖運動中的影響，并將歐拉方法應(yīng)用于三維船模強迫橫搖中，以分離得到波浪阻尼系數(shù)。同時，船模橫搖相關(guān)參數(shù)的變化對于船模橫搖阻尼系數(shù)各個部分的影響也有待進一步的研究。

［1］ CHAKRABARTI S.Empirical calculation of roll damp?ing for ships and barges［J］.Ocean Engineering，2001，28（7）：915-932.

［2］ IKEDA Y，HIMENO Y，TANAKA N.On eddy making component of roll damping force on naked hull［J］. Journal of the Society of Naval Architects of Japan，1977，142：54-64.

［3］ IKEDA Y，HIMENO Y，TANAKA N.Prediction meth?od for ship roll damping：00405［R］.Osaka：Depart?ment of Naval Architecture，University of Osaka，1978.

［4］ ALESSANDRINI B，DELHOMMEAU G.A fully cou?pled Navier-Stokes solver for calculation of turbulent incompressible free surface flow past a ship hull［J］.In?ternational Journal for Numerical Methods in Fluids，1999，29（2）：125-142.

［5］ LUCA B，STEFANO B.Influence of viscosity on radia?tion forces：a comparison between monohull，catamaran and SWATH［C］//Proceedings of the Twenty-third In?ternational Offshore and Polar Engineering.Anchorage， Alaska，USA：International Society of Offshore and Po?lar Engineers，2013.

［6］ MILLER R，GORSKI J，F(xiàn)RY D.Viscous roll predic?tions of a circular cylinder with bilge keels［C］//Pro?ceedings of the 24th Symposium on Naval Hydrodynam?ics.Fukuoka，Japan：NRNAS，2002.

［7］ WILSON R V，PATERSON E G，STERN F.Unsteady RANS CFD method for naval combatants in waves［C］// Proceedings of the 22nd ONR Symposium on Naval Hy?drodynamics.Washington，DC：［s.n.］，2006.

［8］ YANG C L，ZHU R C，MIAO G P，et al.Numerical sim?ulation of rolling for 3-D ship with forward speed and nonlinear damping analysis［J］.Journal of Hydrody?namics，Series B，2013，25（1）：148-155.

［9］ SHEN Z R，CAO H J，YE H X，et al.Manual of CFD solver naoe-FOAM-SJTU：2012SR118110［P］.Shang?hai：Shanghai Jiaotong University，2012.

［10］ JASAK H，JEMCOV A，ZELJKO T，et al.Open?FOAM：a C++library for complex physics simulations［C］//International Workshop on Coupled Methods in Numerical Dynamics.Dubrovnik.Croatia：IUC，2007：1-20.

［11］ ISSA R I.Solution of the implicitly discretised fluid flow equations by operator-splitting［J］.Journal of Computational Physics，1986，62（1）：40-65.

［12］ 吳建威，尹崇宏，萬德成.基于三種方法的螺旋槳敞水性能數(shù)值預報［J］.水動力學研究與進展，2016，31（2）：177-187. WU J W，YIN C H，WAN D C.Numerical prediction of the propeller open-water performance based on three numerical methods［J］.Chinese Journal of Hy?drodynamics，2016，31（2）：177-187（in Chinese）.

［13］ SHEN Z R，WAN D C.An irregular wave generating approach based on naoe-FOAM-SJTU solver［J］.Chi?na Ocean Engineering，2016，30（2）：177-192.

［14］ 李鵬飛，萬德成，劉建成.基于致動線模型的風力機尾流場數(shù)值模擬［J］.水動力學研究與進展，2016，31（2）：127-134.

LI P F，WAN D C，LIU J C.Numerical simulations of wake flows of wind turbine based on actuator line model［J］.Chinese Journal of Hydrodynamics，2016，31（2）：127-134（in Chinese）.

［15］ 尹崇宏，吳建威，萬德成.基于IDDES方法的模型尺度和實尺度VLCC阻力預報與流場分析［J］.水動力學研究與進展，2016，31（3）：259-268.

YIN C H，WU J W，WAN D C.Model-and full-scale VLCC resistance prediction and flow field analysis based on IDDES method［J］.Chinese Journal of Hy?drodynamics，2016，31（3）：259-268（in Chinese）.

［16］ CAO H J，WAN D C.Development of multidirectional nonlinear numerical wave tank by naoe-FOAM-SJTU solver［J］.International Journal of Ocean System En?gineering，2014，4（1）：49-56.

［17］ ZHA R S，YE H X，SHEN Z R，et al.Numerical com?putations of resistance of high speed catamaran in calm water［J］.Journal of Hydrodynamics，2014，26（6）：930-938.

［18］ SHEN Z R，WAN D C，CARRICA P M.Dynamic over?set grids in OpenFOAM with application to KCS self-propulsion and maneuvering［J］.Ocean Engi?neering，2015，108：287-306.

［19］ 楊波，石愛國，吳明.基于計算流體力學理論的船舶橫搖阻尼系數(shù)計算［J］.中國航海，2012，35（3）：76-80.

YANG B，SHI A G，WU M.Calculation of ship roll damping coefficient based on CFD［J］.Navigation of China，2012，35（3）：76-80（in Chinese）.

［20］ 李積德.船舶耐波性［M］.哈爾濱：哈爾濱工程大學出版社，2007.

LI J D.Ship seakeeping［M］.Harbin：Harbin Engi?neering University Press，2007（in Chinese）.

［21］ OLIVIERI A，PISTANI F，AVANZINI A，et al.Tow?ing tank experiments of resistance，sinkageand trim，boundary layer，wake，and free surface flow around a naval combatant INSEAN 2340 model：IIHR-TR-421［R］.Iowa：The University of Iowa，2001.

［22］ 黃昊，郭海強，朱仁傳，等.粘性流中船舶橫搖阻尼計算［J］.船舶力學，2008，12（4）：568-573.

HUANG H，GUO H Q，ZHU R C，et al.Computations of ship roll damping in viscous flow［J］.Journal of Ship Mechanics，2008，12（4）：568-573（in Chinese）.

［23］ KORPUS R A，F(xiàn)ALZARANO J M.Presiction of vis?cous ship roll damping by unsteady Navier-Stokes techniques［J］.Journal of Offshore Mechanics and Arctic Engineering，1997，119：108-113.

Numerical analysis of viscous effect on ship rolling motions based on CFD

LUO Tian1，2，3，WAN Decheng1，2，3
1 School of Naval Architecture，Ocean and Civil Engineering，Shanghai Jiao Tong University，Shanghai 200240，China
2 State Key Laboratory of Ocean Engineering，Shanghai Jiao Tong University，Shanghai 200240，China
3 Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration，Shanghai 200240，China

During the ship design procedure,the analysis of ship rolling motions is of great significance because the rolling motions have extraordinary effects on the sea-keeping,maneuverability and stability of a ship.It is difficult to simulate rolling motions due to the effect of viscosity,which causes many nonlinear components in computation.As such,the potential theory used for other ship motions cannot be used for rolling motions.This paper simulates the rolling motions of the DTMB 5512 ship model and the ship transverse section of the S60 ship model with a naoe-FOAM-SJTU solver using the Reynolds Averaged Navier Stokes（RANS）method based on the OpenFOAM.The results of rolling motions are compared with the experimental data,which confirms the reliability of the meshes and results.For the ship transverse section of the S60 ship model,the damping coefficient is divided into three parts with the Euler and RANS methods:friction,vorticity and wave parts.For the DTMB 5512 ship model,the damping coefficient is also respectively analyzed,including the friction,vorticity,wave and bilge keel parts.The results in this paper show that the vorticity part accounts for the greatest proportion,while the friction part accounts for the least, and the bilge keels reduces the damping moment to a certain extent which shows the effect of rolling parameters on rolling motions and moments.

ship；rolling motions；damping coefficient；damping moment；RANS method；Euler method

U661.3

A

10.3969/j.issn.1673-3185.2017.02.001

http://kns.cnki.net/kcms/detail/42.1755.TJ.20170313.1618.036.html

羅天，萬德成.基于CFD的船舶橫搖數(shù)值模擬與粘性效應(yīng)分析［J］.中國艦船研究，2017，12（2）：1-11，48. LUO T，WAN D C.Numerical analysis of viscous effect on ship rolling motions based on CFD［J］.Chinese Journal of Ship Research，2017，12（2）：1-11，48.

2016-09-01 < class="emphasis_bold"> 網(wǎng)絡(luò)出版時間：

時間：2017-3-13 16:18

國家自然科學基金資助項目（51379125，51490675，11432009，51579145）；長江學者獎勵計劃（T2014099）；

上海高校東方學者特聘教授崗位跟蹤計劃（2013022）；工信部數(shù)值水池創(chuàng)新專項VIV/VIM項目（2016-23/09）

羅天，女，1991年生，碩士生。研究方向：粘性流場中船舶的橫搖運動分析。

E-mail：sky360652272@163.com

萬德成（通信作者），男，1967年生，博士，教授。研究方向：船舶與海洋工程水動力學及計算流體力學的基礎(chǔ)理論及其應(yīng)用。E-mail：dcwan@sjtu.edu.cn

期刊網(wǎng)址：www.ship-research.com