CFD數值模擬船舶在波浪中的回轉操縱運動

王建華 ，萬德成

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

2高新船舶與深海開發裝備協同創新中心，上海200240

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

0 引 言

船舶操縱運動可以反映出船舶在航行過程中的機動性、回轉特性和航向的糾偏能力。操縱性能的優劣與船舶的航行安全和能耗息息相關，其重要性不言而喻。目前，船舶的操縱性能評估主要通過典型的操縱運動試驗進行，其中對船舶自由回轉操縱運動的數值模擬是評估船舶回轉性能的重要手段。船舶自由回轉操縱運動一般通過執行目標舵角，使船舶在舵提供的回轉力和力矩作用下實現回轉運動，典型的自由回轉操縱運動軌跡及特征參數如圖1所示。一般在船舶設計的初期階段，就需要評估所設計船舶的操縱運動特性，尤其是回轉特性，以指導船舶航行中的操作，保障船舶全壽期的安全性能。

圖1 船舶自由回轉操縱運動軌跡及特征參數［1］Fig.1 Free turning maneuvering motion trajectory and main parameters of ship［1］

一般船舶操縱都是通過后舵操作，實現特定的操縱運動。因此，為了精確地評估船舶操縱運動特性，需要考慮到船—槳—舵的相互耦合作用。目前，應用較為廣泛的船舶操縱性預報方法主要有船模試驗方法和基于CFD的數值模擬方法。船模試驗方法是目前應用最廣泛也是最可靠的方法，尤其是近年來得益于試驗裝置和方法的改進，使得自航船模的操縱試驗成為了可能。該方法可以通過在海洋工程波浪水池中或者天然湖泊中開展特定的操縱運動試驗，來預報船舶操縱性能。但該方法需要較大的試驗水池、精確的螺旋槳和舵控制系統以及用于測量船體六自由度運動的設備。此外，為真實還原船舶操縱運動中的實際流場，需要有與實際環境相似的試驗水池，故船模試驗方法的設備和試驗成本高。而且，鑒于當前的流場測量設備（如粒子圖像測速（PIV）等）普及性及適用性有限，試驗中尚給不出操縱運動過程中船體、螺旋槳和舵周圍精細的流場結構，并無法詳細分析在此過程中船舶的水動力性能變化。

基于CFD的操縱運動數值模擬研究可分為約束船模操縱運動和自航船模操縱運動。前者通過結合操縱運動數學模型，依據數值模擬約束船模操縱運動得到各個水動力導數，進而對典型的船舶操縱性能進行仿真，其被廣泛應用于船舶操縱運動的數值模擬。Simonsen等［2］利用自主開發的求解器CFDShip-Iowa，數值模擬了靜態和動態平面運動機構（PMM）試驗，并且分別采用CFD計算和試驗測量得到的水動力導數，根據分離型操縱性數學模型，對標準集裝箱（KCS）船模在靜水中的回轉操縱和Z形操縱試驗進行模擬。Guo和Zou［3］采用商業軟件 STAR CCM+數值模擬了標準船模ONRT的旋轉試驗、靜態斜拖試驗和純橫搖試驗的操縱運動，數值回歸得到了操縱性水動力導數值，運用四自由度的操縱運動數學模型（MMG）仿真了船舶25°自由回轉和20/20 Z形操縱運動，預報的運動軌跡與試驗值吻合較好，驗證了采用CFD方法數值模擬約束船模操縱運動試驗得到的操縱性導數值的可靠性。通過數值構建船—槳—舵整體耦合運動求解模型，進行自航船舶操縱運動的直接數值模擬，可以更精確地描述船舶操縱運動過程。

目前，隨著高性能計算機的快速發展以及重疊網格技術的逐步完善，直接數值模擬自航船模的操縱運動已成為現實。Carrica等［4］采用自主開發的水動力學軟件CFDShip-Iowa V4，模擬了不同航速下（Fr=0.25，0.41）船舶自由回轉（35°舵角）和Z形操縱（20/20）特性，在數值計算中，利用動態重疊網格技術處理大幅度的船舶操縱運動，并對波浪工況下的特性進行計算，其數值預報的操縱性參數與試驗值的誤差在10%以內，指出簡化的螺旋槳體積力模型是預報誤差產生的主要原因，該模型忽略了真實情況下螺旋槳旋轉導致的對船體運動的阻力及其受到的側向力影響等因素。Mofidi和 Carrica［5］采用同樣的求解器，但是考慮了真實情況下旋轉的螺旋槳，進行了典型10/10 Z形操縱試驗和修正型的15/1 Z形操縱試驗的數值模擬，數值預報的船體運動及操縱性參數與試驗結果吻合良好，并且對自航操縱運動過程中的詳細流場進行了分析。Broglia等［6］和 Dubbioso等［7］分別進行了單舵及雙舵情況下雙槳推進船舶的自由回轉試驗數值模擬，其中舵和船體運動采用動態重疊網格進行處理，將得到的船體運動軌跡與試驗結果進行對比，并比較了單舵和雙舵情況下自由回轉運動軌跡和回轉降速、漂角及橫搖等時歷曲線，分析了回轉運動全過程中的舵力和船體、附體的側向力變化，指出在雙槳情況下舵會強烈干擾螺旋槳受到的載荷。Shen等［8］基于開源CFD計算平臺OpenFOAM，開發了船舶水動力學求解器 naoe-FOAM-SJTU［9-11］，引入了重疊網格模塊，并擴展到船—槳—舵相互作用下的船舶自航［12］和操縱運動［13-15］模擬計算中，驗證了采用非結構化網格直接進行帶螺旋槳、帶舵船舶的操縱運動數值模擬的可行性。

綜上所述，雖然結合重疊網格技術的CFD方法已經廣泛應用于船舶操縱運動的直接數值模擬，但是大部分研究都是針對靜水工況下的操縱運動，而海上航行的船舶經常處于波浪環境中，因此很有必要開展波浪工況下操縱運動的精確預報，為船舶設計提供更為精確的數據支撐。本文將采用結合重疊網格技術的CFD求解器naoe-FOAM-SJTU，直接數值模擬帶螺旋槳、帶舵船舶在波浪中的自由回轉操縱運動，預報船舶在波浪中的操縱運動特性。通過數值計算船、槳、舵周圍精細的流場，分析船舶在自由回轉過程中的水動力變化和船、槳、舵干擾以及波浪對船舶回轉性能的影響。

1 數值計算方法

1.1 流體計算控制方程

本文計算域流場求解的控制方程為非定常兩相不可壓縮的RANS方程：

式中：▽為求散度；U為速度場；pd=p-ρg·x，為動壓力，其數值等于總壓力值減去靜水壓力；ρ為液體或者氣體的密度；x為空間坐標；t為時間；g為重力加速度向量；μeff=ρ(ν+νt)為有效動力粘性，其中v為運動粘度，vt為渦粘度;fσ為表面張力項。

湍流模型采用SSTk-ω［16］，該模型兼具標準k-ω和k-ε模型的優點，能夠保證壁面處和遠流場求解的精確性和可靠性。自由面求解采用帶有人工可壓縮項的VOF（Volume of Fluid）方法［17］，兩相VOF輸運方程定義為

式中：Ur為用于壓縮界面的速度場；α為兩相流體的體積分數，代表液體部分所占體積的百分比，取值范圍為0～1，0表示氣體，1表示水，介于0到1之間則表征為自由面位置。因此，通過體積分數α便可以將兩相流規化為統一的流體域。

上述RANS方程（式（1）～式（2）），VOF輸運方程（式（3））和湍流方程都采用有限體積法來進行離散。采用OpenFOAM自帶的離散格式進行方程離散，時間項采用隱式Euler格式，對流項采用二階TVD格式，耗散項采用中心差分格式，VOF方程中對流項采用Van Leer格式離散。流體控制方程求解中，速度壓力解耦采用PISO算法［18］。

1.2 自航船舶操縱運動控制模塊

船舶操縱一般都會有大幅度的船舶運動，傳統的變形網格在模擬物體大幅度運動時網格質量會下降，影響求解精度；而重疊網格技術允許多個相互獨立的網格之間產生無約束的相對位移，在計算過程中能夠保證網格不發生變形，從而保證計算過程中網格的質量，因此非常適用于帶槳、帶舵船舶操縱運動問題的數值求解。采用動態重疊網格技術離散的船、槳、舵多級物體運動模型如圖2所示。螺旋槳和舵根據不同船體運動形式，可以按照指定的控制參數（如螺旋槳轉速、最大轉舵角度等），繞著旋轉軸進行自身的旋轉運動，船體則在槳、舵自身運動以及船體受力情況下在自由面環境下做六自由度的運動。

圖2 船—槳—舵多級物體運動示意圖Fig.2 Diagram of motions in ship-propeller-rudder system

依托重疊網格方法以及多級物體運動模塊，可以很方便地實現自航船舶的操縱運動控制，即通過對舵角的控制，實現特定船舶操縱運動的數值模擬。35°滿舵向右舷進行自由回轉操縱運動的舵角控制方程為

式中：δ(t)為舵角；k為轉舵速率；tp表征進行回舵的時刻，從而結束回轉運動，即何時回到初始零度舵角。在起始時刻，按照轉舵速率k進行轉舵直到滿舵狀態，之后維持此舵角完成回轉操縱運動，根據模擬需求進行回舵操作，結束回轉運動。

1.3 區域造波方法

區域造波方法與速度入口邊界造波方法直接的區別就是前者不僅需要邊界造波，同時還需要在特定的區域范圍內對流場進行改造。具體實現方式是通過采用松弛區域，保證外部邊界處沒有波浪反射，同時還能夠確保計算域內部的波浪反射不會對造波邊界產生干擾，這也是邊界造波方法所不具備的特點。本文采用開源造波工具包waves2Foam［19］在移動計算域中進行波浪場生成。采用環形造波區，如圖3所示，進行回轉操縱運動過程中的波浪生成，環形區域中通過松弛方式即可實現造波，同時也能完成消波。該造波區在計算中可以跟隨計算域進行移動，因此可以保證波浪在船舶進行360°回轉運動過程中可以傳播到整個計算域中。圖中，L為船長。

圖3 區域造波圖示Fig.3 Diagram of wave generation zone

2 計算模型和工況

本文計算船型采用全附體雙槳、雙舵的ONRT船模，該船模是被廣泛應用于CFD驗證的標準船型，被列為Tokyo2015 CFD研討會上的自航模問題的標準船型。對于該船型，有非常豐富的操縱試驗數據，從而可以驗證當前數值預報手段的可靠性。船體的幾何模型如圖4所示。船體的主尺度如表1所示。

圖4 ONRT船幾何模型Fig.4 Geometry model of ONR Tumblehome ship

表1 ONRT船體模型主尺度Table 1 Main particulars of ONR Tumblehome ship model

數值計算中采用重疊網格方法進行船、槳、舵網格的直接劃分，重疊網格的布置如圖5所示。計算域共分為6部分，即背景網格、船體周圍網格、2套螺旋槳網格和2套舵的網格，劃分完成的網格如圖6所示，計算網格總數量為711萬。

圖5 重疊網格布置Fig.5 Overset grid arrangement

圖6 槳和舵周圍網格分布Fig.6 Local grid distribution around twin propellers and rudders

本文進行了在35°舵角下，船舶在波浪中的自由回轉操縱運動的直接數值模擬，船舶初始航速為1.11 m/s，對應于Fr=0.2，數值計算中螺旋槳的轉速設置為對應于這個航速下的模型自航點值，為 529.14 r/min［13］。入射波浪根據 IIHR 的試驗［20］進行設置，入射波浪的波長λ等于船長LWL，波陡H λ為 0.02。

3 波浪中船舶自由回轉數值模擬結果分析

波浪中的自由回轉操縱運動數值模擬從最終穩定的自航數值計算開始，然后開始放開船舶的六自由度運動，與試驗的一致，均在入射波浪的波峰到達船艏時進行操舵，舵按照自由回轉操縱運動進行控制。所有數值計算均在上海交通大學船海計算水動力學研究中心高性能計算集群進行，采用40個進程并行計算，計算時間步長為Δt=0.000 5 s，對應于每個時間步螺旋槳轉過1.5°。完成波浪中的自由回轉操縱運動共耗時1 206 h，對應于155 000個時間步。

圖7所示為數值預報的波浪中船舶自由回轉得到的運動軌跡以及與試驗值［17］的對比。從圖中可以看出，當前的數值預報結果與試驗結果吻合較好，但是數值預報的回轉圈會比試驗的結果更大，這主要是由于數值計算中為了保證重疊網格間足夠的插值單元，而對舵的幾何模型進行了修正，減小了有效的舵面積，因此使得舵效減小。此外，從圖中還可以看出，在船舶航向角改變90°和270°時，回轉曲線會產生明顯的波動現象。圖8所示為對應的局部放大對比圖，在對應的2個時間段，船舶的運動軌跡會產生明顯波動，并且試驗和CFD預報結果都顯示有這一現象。數值計算得到的波動幅值明顯小于試驗中的波動值，這也說明了CFD模擬中船舶的回轉性較試驗稍差，這也解釋了圖7中數值預報的船舶回轉圈更大的原因。

圖7 回轉圈軌跡對比Fig.7 Comparison of turning circle trajectory

圖8 局部回轉運動對比Fig.8 Local comparison of trajectory

數值預報的船舶回轉圈特征參數及其與試驗值［20］的對比如表2所示。

表2 回轉圈軌跡特征參數對比Table 2 Comparison of main parameters of turning circle trajectory

仿真過程中，為了保證對比的可靠性，對時間尺度進行了調整，使得CFD模擬和物理試驗滿足在同一個時刻執行轉舵操作。從表中的對比結果可以看出，所有特征參數與試驗值的誤差均在10%以內，當前的數值計算可以較高的精度預報出波浪中自由回轉船舶的操縱運動特性。

圖9所示為數值預報船舶在波浪中自由回轉過程中六自由度運動的時歷曲線。從圖中可以看出，船舶的垂蕩、縱搖和橫搖運動會產生較為明顯的波頻振蕩特性（圖9（a），（b）（c））。另外，由于回轉運動過程中船舶遭遇的浪向角也在時刻變化，因此在高頻的波頻運動下還存在由于回轉操縱運動導致的低頻波動。整個回轉運動過程中，最大的縱搖幅值可達2.5°，橫搖運動的幅度為-4.4°～8°。此外，從橫搖運動的時歷曲線（圖9（c））可以看出，由于波浪導致的橫搖運動幅值較由于初始操舵導致的橫搖運動更大。而3個平面的運動，縱蕩、橫蕩和艏搖運動則展現較小的波頻運動特性。從艏搖運動的曲線（圖9（d））上可以看出較小的波動，這也可能導致了圖8中展現的平面運動軌跡中的局部波動。

圖9 波浪中自由回轉操縱中船舶六自由度運動時歷曲線Fig.9 Time history curves of ship motions for turning circle maneuvering in waves

圖10所示為船舶在波浪中自由回轉運動過程中的航速以及艏搖速率的變化曲線。從圖中可以看出，波浪中的船舶在回轉運動過程中會出現明顯的回轉降速現象，并且最大降速可達40%。初始的航速降低是由于轉舵導致的，進入回轉運動以后會維持在平均30%的降速范圍內。而對于艏搖速率來說，初始的明顯速率變化是由于受到轉舵的影響，而后期的波動則是由于船舶遭遇變化的浪向導致，最大的艏搖速率可達12.2（°）/s。

圖10 波浪中自由回轉操縱中船舶航速和首搖速率時歷曲線Fig.10 Time history curves of ship speed and yaw rate for turning circle maneuvering in waves

圖11所示為船舶在波浪中進行自由回轉操縱運動過程中螺旋槳推力和扭矩的變化曲線。從圖中可以看出，螺旋槳的推力和扭矩呈現出明顯的波頻振動特性，這主要是由于船舶運動過程中導致螺旋槳的進流產生變化，進而使得推進性能產生波動。從局部放大圖中可以看到更為高頻的振蕩現象，這是由于真實旋轉螺旋槳葉片切割流場導致。

圖11 波浪中自由回轉操縱中船舶推進性能時歷曲線Fig.11 Time history curves of propulsion coefficients for turning circle maneuvering in waves

圖12所示為轉舵時刻舵所受到的水動力載荷的變化曲線。從圖中可以看出，執行轉舵操作之前，作用在兩側舵上的阻力基本一致，并且側向力對稱，而執行完操舵以后，舵阻力增加明顯，而側向力變成同向，產生較大的側向合力，而側向力的合力也使得船舶產生回轉運動。

圖12 轉舵過程中舵受到的水動力時歷曲線Fig.12 Time history curves of rudder forces during rudder execution

圖13所示為在轉舵過程中的槳舵周圍的渦量場變化。從圖可以看出，初始時刻，舵角為0°時，槳舵周圍的渦量分布基本為對稱形式，而隨著舵角的增加，舵對前面螺旋槳的瀉渦會產生明顯的干擾，由于舵向左舷轉動，因此左舷舵會對螺旋槳的槳轂渦產生干擾，而右舷舵則會影響到右舷槳的葉梢渦。這種現象的區別也解釋了圖11和圖12中兩側螺旋槳和舵水動力的區別。而舵周圍會發生明顯的流動分離現象，但是現在采用的是RANS方法，無法精確地捕捉這種情況下的周圍流動，因此會對舵力的計算產生誤差，這也是導致目前計算中的回轉圈變大的原因之一。

圖13 轉舵過程中槳舵周圍的渦量場Fig.13 Snapshots of vortical field around twin propellers and rudders during rudder execution

圖14所示為船舶在波浪中自由回轉過程中，4個典型時刻的自由面波形變化，分別對應于0°、120°、240°和 360°航向角的時刻。從圖中可以看出，在沒有轉向時，船舶周圍的波浪環境基本對稱，但是在回轉角度達到360°時，船艏和船艉處均能看出由于轉動導致的兩側波面的差別；而從120°和240°航向角時的自由面可以看出，兩側波面存在明顯的高度差別，這也導致了船體兩側的壓力分布不均；從圖14（d）同樣可以看出船艏會抬出水面，這證明了在該波浪情況下船舶會產生大幅度的六自由度運動。

圖14 波浪中自由回轉過程中船舶周圍自由面Fig.14 Snapshots of wave elevation around ship hull for turning circle maneuvering in waves

4 結 語

本文采用結合重疊網格技術的CFD求解器naoe-FOAM-SJTU，對船—槳—舵相互作用下的波浪中船舶自由回轉操縱運動進行了直接數值模擬。數值預報的波浪中船舶回轉運動的回轉圈特征參數（如縱距、橫距、戰術直徑、回轉直徑等）與已有試驗結果吻合較好，誤差均在10%以內，驗證了當前求解器對船—槳—舵相互作用下的波浪中船舶自由回轉操縱運動預報的適用性和可靠性。此外，根據計算結果顯示，船舶的垂蕩、縱搖和橫搖運動展現出明顯的波頻運動響應，而縱蕩、橫蕩和艏搖3個平面運動的波頻振動特征不明顯。在波浪中船舶進行自由回轉時的最大船舶失速可達40%。同時，給出了整個操縱運動過程中的推進性能和舵力的變化。并且通過詳細的流場信息，如不同時刻自由面變化和槳、舵周圍渦量場變化等，分析了波浪中回轉操縱運動下水動力變化的原因。

由于當前數值模擬采用時均的RANS方法進行流場求解，因此對于槳、舵周圍大分離流動現象捕捉的精度較差，這也導致了目前的數值預報存在一定的誤差，將來的工作將主要開展基于更為精確的分離渦模擬方法進行該問題的求解，以給出更精細的流場模擬，獲取更高精度的數值預報結果。