王建華 ，萬德成

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

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

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

0 引 言

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

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

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

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

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

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

1 數(shù)值計算方法

1.1 流體計算控制方程

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

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

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

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

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

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

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

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

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

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

1.3 區(qū)域造波方法

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

圖3 區(qū)域造波圖示Fig.3 Diagram of wave generation zone

2 計算模型和工況

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

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

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

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

圖5 重疊網(wǎng)格布置Fig.5 Overset grid arrangement

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

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

3 波浪中船舶自由回轉(zhuǎn)數(shù)值模擬結(jié)果分析

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

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

圖7 回轉(zhuǎn)圈軌跡對比Fig.7 Comparison of turning circle trajectory

圖8 局部回轉(zhuǎn)運動對比Fig.8 Local comparison of trajectory

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

表2 回轉(zhuǎn)圈軌跡特征參數(shù)對比Table 2 Comparison of main parameters of turning circle trajectory

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 語

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

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