劉李為，余嘉威，馮大奎*，張志國，陳美霞

1 華中科技大學 船舶與海洋工程學院，湖北 武漢 430074

2 華中科技大學 船舶和海洋水動力湖北省重點實驗室，湖北 武漢 430074

0 引 言

水面艦船在隨浪或艉斜浪中航行時，由于遭遇周期的增大，會使船舯位于波浪波峰處的時間變長，致使其穩(wěn)性力臂顯著降低，進而導致穩(wěn)性失效模式的發(fā)生。這一失穩(wěn)模式被定義為純穩(wěn)性喪失，是船舶在波浪中典型的3 種傾覆模式（參數(shù)橫搖、純穩(wěn)性喪失和騎浪橫甩）之一[1]，涉及非線性橫搖運動、非線性復原力矩、非線性波浪力矩及其與船體的相互作用等強非線性動力學問題[2]。因此，由純穩(wěn)性喪失引發(fā)的船舶傾覆現(xiàn)象成為國際海事組織（IMO）針對船舶第二代完整穩(wěn)性的研究問題之一。

在IMO 船舶與設(shè)計建造分委會（SDC）第四次會議中，完整船舶典型非線性失穩(wěn)模式的直接穩(wěn)性評估方法成為關(guān)注的重點。直接穩(wěn)性評估方法的建立能夠以更加精確的方式考慮所有相關(guān)因素，可以有效實現(xiàn)不同船舶的安全隱患評估，這在船舶設(shè)計階段是較為重要的。目前，國內(nèi)外有部分學者針對純穩(wěn)性喪失的預報方法以及運動特性進行了較為深入的研究。朱軍等[3]針對船舶在波浪中的大傾角穩(wěn)性問題，計算得到了規(guī)則波浪中船舶的穩(wěn)性變化規(guī)律，結(jié)果顯示船舶在正橫浪與斜浪條件下穩(wěn)性大幅度喪失的主要原因是波面相對船體橫剖面的不對稱性。魯江等[4]構(gòu)建了一種縱蕩?垂蕩?橫搖?縱搖四自由度的耦合運動模型，其運用該模型對隨浪規(guī)則波中的船舶穩(wěn)性喪失進行直接評估，分析了縱蕩運動對純穩(wěn)性喪失預報的影響。Umeda 等[5]基于一種縱蕩?橫蕩?橫搖?艏搖四自由度的耦合運動模型，對艉斜浪中船舶的純穩(wěn)性喪失以及傾覆現(xiàn)象進行了模擬。Lu 等[6]針對艉斜浪中船舶的純穩(wěn)性喪失問題分別進行了模型試驗與六自由度數(shù)學模型預報研究，并進一步闡明了艉斜浪中船舶艏搖與橫搖之間耦合作用的重要性。

以上研究主要是以基于勢流理論的數(shù)學模型為預報方法，通過模型試驗，將船舶操縱運動水動力導數(shù)、線性以及非線性橫搖阻尼的測量作為參數(shù)輸入來對船舶在波浪中的失穩(wěn)運動進行直接評估。由IMO 對船舶典型失穩(wěn)模式直接評估方法的建議與要求可知，穩(wěn)性直接評估方法需要能對波浪中船舶穩(wěn)性變化、橫搖周期特性、主參數(shù)共振特性、騎浪平衡、波浪中的轉(zhuǎn)艏特性、多自由度耦合運動以及艉斜浪中的直航特性等水動力性能實現(xiàn)可靠的模擬?；谌珪r域黏性流的數(shù)值模擬（CFD）方法正好能夠滿足上述要求，且無需依據(jù)試驗測量或是經(jīng)驗公式估計的數(shù)據(jù)作為參數(shù)輸入，因此可以考慮將其應用至船舶非線性失穩(wěn)模式的直接評估中。

隨著高性能計算技術(shù)的飛速發(fā)展，直接CFD預報方法已能應用到艦船復雜水動力性能研究中。目前，CFD 方法已被用于船舶的回轉(zhuǎn)運動[7]、停船操縱[8]、騎浪運動[9]以及參數(shù)橫搖現(xiàn)象[10]的直接模擬，說明直接模擬純穩(wěn)性喪失這類集船舶操縱性與耐波性于一體的強非線性問題已成為可能。因此，本文將基于CFD 方法對帶螺旋槳、帶舵的自航船舶在艉斜浪下保持航向操縱運動進行直接模擬，預報典型失穩(wěn)模式下船舶的六自由度耦合運動，并對船舶的失穩(wěn)運動及傾覆特性進行直接評估。

1 數(shù)值方法

1.1 流場求解方法

使用課題組自主開發(fā)的艦船水動力CFD 求解器HUST-Ship 進行數(shù)值仿真計算，通過對非定常RANS 方程的求解來獲取流場特征，然后采用有限差分法進行數(shù)值離散，并基于結(jié)構(gòu)網(wǎng)格模型在每個網(wǎng)格點上對差分方程進行求解?？刂品匠倘缦拢?/p>

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

CFD 求解器在對流場壓力速度特性進行求解的同時，也會耦合剛體六自由度運動方程來進行力和力矩以及運動的實時預報。考慮到船舶在波浪中的大幅度運動、船后螺旋槳旋轉(zhuǎn)以及艉舵操縱運動特性，在CFD 求解器中應用了重疊網(wǎng)格技術(shù)（HUST-Overset），用以實現(xiàn)船體網(wǎng)格模型與計算域之間的大幅相對運動，以及自航船舶船、槳、舵多級耦合運動的數(shù)值仿真。另一方面，重疊網(wǎng)格模塊的應用允許在前處理網(wǎng)格生成過程中單獨對不同的對象（計算域、船體、各附體以及推進器）進行網(wǎng)格劃分，并通過重疊區(qū)域的構(gòu)建實現(xiàn)網(wǎng)格之間的組裝和數(shù)值計算模型的生成，這大大降低了全附體船舶結(jié)構(gòu)化網(wǎng)格的生成難度，提高了結(jié)構(gòu)化網(wǎng)格的質(zhì)量。

重疊網(wǎng)格模塊的重點是生成不同網(wǎng)格塊之間的物面重疊以及體重疊區(qū)域。在HUST-Overset中，首先通過洞映射方法進行洞點的處理以及附近插值點的識別，包括洞邊界點與外邊界點；然后，基于交替數(shù)字樹算法（ADT）為插值點搜尋相對應的貢獻點，其中插值點與貢獻點之間通過三線性插值方法進行流場數(shù)據(jù)的傳遞；最后，基于插值點與貢獻點之間的體積匹配及點距離進行重疊區(qū)域的優(yōu)化，在保證網(wǎng)格單元體積匹配且至少2 層網(wǎng)格互相重疊的前提下盡可能縮小重疊區(qū)域的面積。

1.3 操縱控制模塊

基于動態(tài)重疊網(wǎng)格方法，自航船舶在波浪中進行操縱運動過程中的船、槳、舵多級耦合運動得以實現(xiàn)。在CFD 求解器中集成了操縱運動控制模塊，可以實現(xiàn)自航船舶在波浪條件下的航向保持操縱運動，然后通過使用PID 控制技術(shù)對螺旋槳轉(zhuǎn)速進行控制來使船舶達到目標航速，并以船舶艏向角為輸入對艉舵操縱角度進行反饋控制，從而保證船舶不偏離目標航向。基于PID 的控制形式如下：

式中：P，I，D分別為比例項、積分項和微分項的控制常數(shù)；e為實測值與目標值之間的差值。對于自航船舶在波浪條件下的航向保持運動，在計算過程中可通過不斷調(diào)整舵角角度來使船舶在航行過程中保持目標航向。

1.4 移動數(shù)值試驗水池

在船舶大尺度操縱運動（例如，自由自航運動、回轉(zhuǎn)運動和Z 形操縱運動）的數(shù)值模擬過程中，如果類似于物理試驗水池一樣采用靜止的數(shù)值水池進行仿真，會大大增加數(shù)值計算的網(wǎng)格量，從而降低計算速率。因此，可在CFD 求解器中采用隨船一起運動的移動數(shù)值試驗水池模型，以便節(jié)省數(shù)值水池網(wǎng)格量，提高數(shù)值計算的效率。

在移動數(shù)值試驗水池模型的實現(xiàn)過程中，將船舶主體網(wǎng)格和數(shù)值水池網(wǎng)格分別予以處理：通過計算得到船舶受力和力矩，結(jié)合船舶六自由度運動方程的求解得到船舶在大地坐標系下的姿態(tài)和位置變化，提取船舶的縱蕩和橫蕩速度以及位移信息，并賦值于數(shù)值試驗水池上，從而保證在數(shù)值模擬過程中數(shù)值水池與船舶主體網(wǎng)格之間的縱向和橫向相對位移足夠小，具體實現(xiàn)流程圖如圖1 所示。在此過程中，可將速度入口造波方式與移動數(shù)值水池模型進行耦合處理，并將數(shù)值試驗水池的速度入口邊界與入射波浪的浪向設(shè)置為始終保持垂直，用以保證數(shù)值試驗水池在移動過程中水池內(nèi)的數(shù)值波浪模擬精度不受影響。圖2直觀地示出了船舶在不同類型操縱運動過程中數(shù)值波浪水池的移動方式。

圖1 移動數(shù)值波浪水池模型與船舶運動耦合求解流程圖Fig. 1 Flowchart of the coupling solution between the moving numerical wave tank and ship motions

圖2 船舶不同操縱運動過程中數(shù)值波浪水池的移動方式Fig. 2 The movement mode of the moving numerical wave tank during different ship maneuvering motions

2 計算對象與工況

選取全附體雙槳、雙舵的ONRT 船模為研究對象，該船模為被廣泛應用于數(shù)值仿真算法驗證的標準模型。文獻[6]以O(shè)NRT 自航船模為對象對其進行了純穩(wěn)性喪失模型試驗，相關(guān)的試驗數(shù)據(jù)可用于驗證本文CFD 方法的可靠性。模型的縮尺比為40.526，其幾何形狀如圖3 所示，主尺度如表1 所示。根據(jù)文獻[6]中的自航模試驗工況，對自航船舶在30°艉斜浪規(guī)則波中的航向保持操縱運動進行數(shù)值模擬，其中規(guī)則波波長λ =1.25LWL，波高H= 0.05λ，波速C= 2.722 m/s，波浪頻率f= 0.573 Hz。船模航速為1.831 m/s，對應的弗勞德數(shù)Fr= 0.30。在數(shù)值模擬過程中考慮船舶的六自由度運動，并通過航向角反饋控制器不斷操控艉舵舵角來使船舶保持航向航行。

表1 ONRT 模型主尺度參數(shù)Table 1 Main parameters of ONRT model

圖3 ONRT 模型幾何圖Fig. 3 Geometry of ONRT model

在數(shù)值模擬的前處理過程中，分別對計算域、船體以及各附體（呆木、艉舵、舭龍骨、軸系及螺旋槳）進行結(jié)構(gòu)化網(wǎng)格劃分，計算域網(wǎng)格與數(shù)值模擬邊界條件的設(shè)置如圖4 所示。全附體模型網(wǎng)格如圖5 所示。運用重疊程序HUST-Overset 實現(xiàn)網(wǎng)格組裝與重疊區(qū)域的生成，所建立的數(shù)值計算網(wǎng)格模型如圖6 所示。在進行數(shù)值模擬之前，對螺旋槳網(wǎng)格進行網(wǎng)格收斂性分析，具體的結(jié)果參見文獻[14]，從中可知中等數(shù)量網(wǎng)格與密網(wǎng)格方案下的螺旋槳敞水特性預報結(jié)果差別較小，且與試驗結(jié)果間的誤差均小于2.5%，而稀疏網(wǎng)格方案下螺旋槳敞水特性預報結(jié)果與試驗結(jié)果間的誤差則為6%。由于全附體模型網(wǎng)格收斂性分析結(jié)果與螺旋槳網(wǎng)格收斂性分析結(jié)果的現(xiàn)象一致，因此選用中等數(shù)量網(wǎng)格方案進行數(shù)值模擬。根據(jù)網(wǎng)格收斂性結(jié)果選取了最終適用的計算網(wǎng)格模型，網(wǎng)格節(jié)點總數(shù)量為978 萬。

圖4 計算域網(wǎng)格Fig. 4 Grid of the computational domain

圖5 全附體模型網(wǎng)格Fig. 5 Grid of the fully appended ship model

圖6 數(shù)值計算網(wǎng)格模型Fig. 6 Numerical grid model

3 數(shù)值模擬結(jié)果分析

基于CFD 求解器HUST-Ship，耦合動態(tài)重疊網(wǎng)格方法，對ONRT 模型在30°艉斜浪規(guī)則波（λ = 1.25LWL，H= 0.05λ）中的典型失穩(wěn)模式純穩(wěn)性喪失進行數(shù)值模擬。在數(shù)值模擬過程中放開了船舶的六自由度運動，并通過航向角反饋控制器不斷操控艉舵舵角來使船舶保持航向航行。在文獻[15]中，IIHR 水池運用PID 控制技術(shù)中的比例控制器（P= 1，I=D= 0）對ONRT 模型在不同浪向下的航向保持進行了試驗，并取得了良好的效果。因此，本文進行數(shù)值模擬的控制設(shè)定將與該試驗中的控制方法一致。

由文獻[6] 中的自航模試驗結(jié)果可知，模型在Fr達到0.30 時發(fā)生了強烈的穩(wěn)性喪失現(xiàn)象，并導致船舶出現(xiàn)了傾覆。因此，在數(shù)值模擬過程中將參考該工況，固定螺旋槳轉(zhuǎn)速以使模型航速與試驗保持一致。數(shù)值計算工作在國家超級計算廣州中心高性能計算機集群上進行，共采用120 個計算核心進行了并行計算，所花費的總計算時間約208 h，對應的CPU 計算時間為2.5 萬計算機核時。圖7 所示為純穩(wěn)性喪失過程中船后周圍的渦結(jié)構(gòu)。模型在波浪中自由自航時航速的時歷曲線如圖8 所示（圖中，U為船舶航速）。由圖可知，在該工況下船舶的瞬時航速接近于波速在船長方向的分量，這就意味著船舶與波浪之間的相對速度較小，即船舶在航行過程中，在波峰處的停留時間非常長，導致船舶在艉斜浪中的穩(wěn)性長時間降低。由于缺乏足夠的復原力，在持續(xù)的外部橫傾力矩作用下，船舶會發(fā)生純穩(wěn)性喪失現(xiàn)象，從而導致大幅度的非線性橫搖運動。

圖7 純穩(wěn)性喪失發(fā)生時船后周圍渦結(jié)構(gòu)Fig. 7 Vorticity structure around the ship stern when the pure loss of stability occurs

圖8 船舶航速時歷曲線Fig. 8 Time histories of the ship speed

圖9 所示為自航船舶在艉斜浪規(guī)則波中的縱蕩、橫蕩以及垂蕩運動時歷曲線（圖中，x，y，z分別代表縱蕩、橫蕩與垂蕩），船舶在螺旋槳的推進作用下始終保持向前航行。圖10 所示為自航船舶在艉斜浪規(guī)則波中的橫搖、縱搖以及艏搖運動時歷曲線（圖中， ?，ψ，θ 分別為橫搖角、艏搖角和縱搖角）。由圖10 可知，船舶在持續(xù)的穩(wěn)性喪失情況下發(fā)生了大幅度的橫搖運動，最終，由于橫搖角過大而導致了傾覆現(xiàn)象的發(fā)生。與此同時，船舶的艏搖角也顯著增大，并與船舶橫搖運動發(fā)生了耦合作用。當船舶橫搖角達90°時，其艏搖角達到了16°，這說明船舶在穩(wěn)性喪失狀態(tài)下依靠操舵已無法有效控制船舶航向，從而導致了較為明顯的橫甩現(xiàn)象的發(fā)生。最終，純穩(wěn)性喪失、大幅度橫搖、劇烈的橫甩以及外部載荷的耦合作用導致了自航船舶的傾覆。

圖9 船舶縱蕩、橫蕩與垂蕩運動的時歷曲線Fig. 9 Time histories of the ship surge, sway, and heave motion

圖10 船舶橫搖、縱搖與艏搖運動的時歷曲線Fig. 10 Time histories of the ship roll, pitch, and yaw motion

為了進行數(shù)值模擬可靠性的驗證，將船舶橫搖、縱搖以及艏搖角時域信號的數(shù)值模擬結(jié)果（CFD）與文獻[6]中的模型試驗數(shù)據(jù)（EFD）進行了對比，結(jié)果如圖10 所示。由圖可知，CFD 與EFD結(jié)果之間吻合較好，且這2 種方法對船舶運動隨時間變化趨勢的預報非常一致，說明采用CFD 方法模擬船舶失穩(wěn)運動與傾覆現(xiàn)象準確性較高。

在CFD 數(shù)值模擬過程中，不同時刻下船舶姿態(tài)的實時變化及其周圍自由面的波形如圖11 所示，圖中，Te為波浪遭遇周期。由圖可以看出，艉斜浪的不對稱性會造成船舶左右舷所受的外部載荷不同，從而導致船舶發(fā)生不對稱的橫搖運動，船艏處的渦結(jié)構(gòu)同樣也表現(xiàn)出了不對稱性。同時，由圖還可以發(fā)現(xiàn)，采用CFD 方法能準確模擬出船舶因穩(wěn)性喪失而導致實時橫搖直至最終傾覆的過程，以及在橫搖過程中耦合發(fā)生的橫甩現(xiàn)象。這表明運用CFD 方法能夠直接模擬船舶的失穩(wěn)運動與傾覆現(xiàn)象，且模擬結(jié)果已較為接近真實場景，有利于船舶典型失穩(wěn)模式直接評估工作的開展。

圖11 不同時刻下船舶姿態(tài)的實時變化及其周圍自由面波形Fig. 11 The real-time change of ship attitude and the free surface wave around the ship at different time instances

4 結(jié) 語

本文基于黏性流理論建立了CFD 數(shù)值模擬方法，然后結(jié)合動態(tài)重疊網(wǎng)格技術(shù)與操縱運動反饋控制模塊，對帶螺旋槳和帶舵的自航船舶ONRT在艉斜浪下保持航向操縱運動進行了直接模擬，預報了純穩(wěn)性喪失失穩(wěn)模式下船舶的六自由度運動，并對船舶的傾覆特性進行了直接評估。結(jié)果表明，運用CFD 方法模擬船舶失穩(wěn)運動與傾覆現(xiàn)象可靠性和準確性較高。

此外，CFD 數(shù)值模擬結(jié)果表明，船舶在持續(xù)的穩(wěn)性喪失情況下會發(fā)生大幅度的橫搖運動，最終會因橫搖角過大而導致傾覆現(xiàn)象的發(fā)生。同時，船舶的艏搖角也將顯著增大，說明船舶在穩(wěn)性喪失狀態(tài)下依靠操舵已無法有效控制船舶航向，從而發(fā)生了較為明顯的橫甩現(xiàn)象，這也是船舶發(fā)生傾覆的關(guān)鍵因素之一。

未來，將借助直接CFD 數(shù)值模擬方法對不同航速下自航船舶在艉斜浪中的失穩(wěn)運動特性進行分析與比較，并對船舶由純穩(wěn)性喪失引發(fā)的大幅度非線性失穩(wěn)運動的發(fā)生機理進行研究。