大方尾三體船型阻力及航態預報

劉 楊，張海華，苗 飛，許凱瑋，張 旭

(中國船舶科學研究中心 上海分部，上海 200011)

0 引 言

在現代船舶中，方尾由于具有增加虛長度、減小阻力、增大排水體積等優點，無論是在單體船舶還是高速三體船舶上都已經得到了廣泛的應用。尤其對于高速三體船，由于大多采用了噴水推進，相比與傳統螺旋槳推進在高速段擁有更好的推進效率，并且利于尾部的布置，也使得方尾的尺寸相對于螺旋槳推進船型更大，在阻力與航態預報中的重要性也更大。因此，在三體船型設計中，能夠準確、快速地對方尾船型的興波阻力、總阻力以及航行姿態進行預報，在方尾船型的設計和優化過程中具有十分重要的作用。

船舶興波阻力及航態的預報方法經過多年發展，主要可分為基于勢流方法和粘流方法。隨著21世紀以來計算機計算效率大幅提高，粘流方法有了顯著發展。但是，粘流方法網格劃分和計算機要求較高。Sherbaz研究表明，與表面壓力和力矩相比，粘性力對航態影響較小，大多在10%以內[16]。在勢流理論的范圍內，對于方尾船型的計算所運用的方法主要有Michell積分法、Havelock源法，Dawson方法和Rapid方法。改進后的Michell積分雖然可以用于方尾船型的計算，但是無法對航態及近場興波進行模擬；Havelock源法運用于高速方尾船舶計算，計算精度較低；而全非線性方法，在方尾船型計算時，自由面迭代容易產生發散。計算近場興波的線性興波理論由Dawson[5]于1977年提出，將速度勢分解為疊模速度勢和擾動速度勢，并對自由面條件進行展開，采用迎風差分格式來實現興波的從前向后傳播；陳京普、Tarafdera等[4,8-10]改進了近場興波方法，采用貼水線自由面網格提高了計算精度及穩定性，使得方法應用更加廣泛，并運用于約束模的興波模擬。Peng等[15]改進自由面船尾后的計算網格，減少了肥大船型尾部流場紊亂及不穩定的問題。由于非線性方法避免了升高面元及自由面迭代等問題，也被廣泛應用于興波計算及船型優化等領域[8-15]。

本文基于勢流理論，針對具有實用價值的大方尾三體船，運用改進后Dawson方法和Rankine源，建立三維的興波阻力及航行姿態的預報方法。針對方尾船型的自身特點，對自由面條件求解的數值計算方法進行了一定的改進，采用自然迭代形成的方尾后尾流區域，避免了虛長度估算所帶來的誤差。由于直接對物體表面進行壓力積分，本方法運用與多體船舶計算時，可以直接反映出片體間的相互作用，同時也可以更為直觀地得到船舶周圍近場的興波情況，利于船型優化工作。本文運用該方法分別對單體、三體方尾船型進行理論預報，并將預報結果與實驗結果進行對比分析，對方法的優缺點進行了具體的研究。

1 數學模型

1.1 基本模型

采用隨船坐標系，坐標系如圖1所示。坐標系原點o位于船舶中心，z軸垂直向上，x軸由船首指向船尾，y軸指向右舷為正，xy平面與靜水面重合。

圖1 坐標系Fig. 1 Coordinate system

假設流體理想，流動定常，無限水深，船舶興波為微幅波。流場速度勢 ψ分解為無自由面影響的重疊模速度勢 φ和考慮自由液面作用的擾動速度勢φ，即

由于重疊模速度勢中包含了來流速度勢，因此假設，φ＜＜φ。流場速度勢 ψ對應的波高為ζ，重疊模速度勢 φ在z=0 處所產生切向速度對應的波高為ξ，擾動速度勢 φ對應的波高為η，且滿足以下關系式：

流場中有

物面上滿足

速度勢 ψ在自由面上分別滿足運動學條件和動力學條件如下：

由式(6)和式(7)可得在z=ζ(x,y)上滿足自由面條件為：

速度勢滿足無窮遠處趨于零的遠方輻射條件。

1.2 自由面邊界條件

考慮到φ關于靜水面對稱，則有：

將式(8)在z=ξ 和z=0 上進行泰勒展開，忽略 φ的非線性項以及高階項。在傳統的流線假設中認為但是并沒有充分的理論來證明擾動速度勢方向與重疊模重疊模相同，因此不再采用傳統的流線假設。可得重疊模流線速度勢表示的自由面條件為：

2 數值計算

2.1 數值處理

計船體網格數為Nb，自由面網格數Nf。在船體網格離散時以四邊形網格為主，船體首部傾斜處采用三角形網格進行處理。自由面采用貼面網格，方尾后網格單獨生成，并進行一定程度的加密。船體處網格長度保持不變，前后以rx的比例遞增，方尾處網格寬度不變，ry的比例遞增。不同的航速所產生的波長不同，即選擇不同的初始網格大小，以得到更加精確的自由面波形以及計算結果。某單體方尾船型自由面網格示意圖如圖2所示。

p

離散單元布置源匯后，物面控制點 重疊模流場速度勢勢表示為：

圖2 自由面網格Fig. 2 Free surface mesh

其中：Sb為物面；σB0為物面控制點處源強；rpq為任意源點q到控制點p的距離。

對式(12)離散，考慮物面條件條件可得：

擾動速度勢表示為：

為了保證自由液面不發生震蕩以得到更加光順、貼合實際的自由面波形，進而使計算結果更加精確，在計算時對自由面單元的控制點進行升高，升高高度與波長和各個單元各自的長度成一定比例。同時，為了更好滿足遠方輻射條件，使得向前差分不受單元對控制中心切向誘導速度為零對數值計算造成的影響，將控制點統一前移一定距離，移動距離為單元長度的0.1倍。

2.2 方尾后自由面條件處理及靜水力計算

由于航速較低的情況下，方尾并未完全出水，因此采用理論計算方法進行數值模擬時，必須對方尾在低速時未出水部分所受到的靜水力進行估算。通過經驗公式，對方尾船型在不同航速下方尾中心及兩側的吃水情況進行估算，并對方尾網格重新加密劃分，積分得到方尾所受到的靜水力，進而提高阻力計算的精度，具體公式如下：

其中：ηdry為方尾中心與兩側未浸水長度系數；Bt，Dt為方尾的寬度與吃水；Ft為方尾傅汝德數；RNT為方尾雷諾數。參照文獻數據對系數C1，C2，C3，C4進行取值，求得方尾中心與兩側的吃水后，對方位進行靜水力修正。

2.3 興波阻力及航態計算

船體表面控制點壓力系數為：

其中，Zi為控制點到靜水面間距離。

船體受到的興波阻力Rw、興波阻力系數Cw、升力Fz和縱傾力矩Ny表示為：

其中：Si表示船體的濕表面積；xi，zi為控制點坐標；xg，zg為船舶的重心坐標；nxi，nzi為單元法線x，z方向分量；Rhx，Rhx，Nhy為傅汝德數Fn=0時由于船體表面的網格離散劃分使船體受到的靜水阻力、升力、力矩修正值。

給定時間步長Δt，可以得到船舶下一時刻的吃水St+Δt以及縱傾值Tt+Δt為：

僅對船體水線以下近水面處局部網格進行實時劃分，而大部分網格始終保持不變，這樣既減小了網格重新離散帶來的誤差，也節省了計算時間。通過反復迭代直到滿足收斂條件（Fz＜ε，Ny＜ε）即可得到考慮船體姿態變化影響的興波阻力，進行船體姿態預報。

3 算例及結果分析

3.1 數值計算及驗證

以具有較大方尾的某三體船為例，對該船型的興波阻力、航行姿態以及方尾吃水進行理論預報并與相應的試驗結果對比。在總阻力計算時，參照文獻[12]取粘壓阻力系數為0.1。

計算過程中，片體間距L，B，D為主體水線長、水線寬和吃水，B/L=0.08，D/L=0.04；L1，B1，D1為側體水線長、水線寬和吃水，B1/L1=0.05，D1/L1=0.05。BT，DT分別為主體方尾水線寬和吃水，Bt/Dt=3.65；Bt1，Dt1分別為側體方尾水線寬和吃水，Bt1/Dt1=3.55。在計算過程中，保持片體與主體尾部平齊，即a/L=0.0，片體間距p/L=0.1，三體船片布局示意圖如圖3所示。

圖3 三體船片體布局示意圖Fig. 3 Layout diagram of trimaran hull

主體縱向網格60個，垂向網格10個；側體縱向網格30個，垂向網格5個。自由表面網格采用隨船長傅汝德數變化而變化，單個波長保持25個網格，船首、尾及橫向網格保持1.05的擴張比例，動態保證計算域船首往前1L，船尾往后1.5L，橫向1L。船體網格如圖4所示，在Fn=0.3時自由面網格如圖5所示。

圖4 船體表面網格示意圖Fig. 4 Schematic diagram of hull surface grid

圖5 自由面網格示意圖(Fn = 0.3)Fig. 5 Schematic diagram of free surface grid (Fn = 0.3)

將該船型興波阻力計算結果與試驗結果進行對比，升沉、縱傾、阻力計算結果如圖6～圖8所示。

圖6 升沉計算結果與試驗值對比Fig. 6 Comparison of heave calculation results and test values

圖7 縱傾計算結果與試驗值對比Fig. 7 Comparison of trim calculation results and test values

從圖6和圖7可以看出，理論方法對于方尾三體船舶航行姿態的模擬有著較高的精度，變化趨勢與試驗值基本一致，幅值也較為接近，但是都略小于試驗結果，這可能是由于忽略了邊界層粘性切向力的影響。文獻[13]也認為隨著航速的增大，粘性作用對于航行姿態的模擬存在10%左右的影響。由圖8可以看出，理論計算方法引入了方尾處理后，對于傅汝德數小于0.35的方尾未完全出水的模擬具有一定精度。而隨著航速增大，在傅汝德數大于0.35方尾完全出水后，計算精度也相對更高，誤差可以控制在2%以內。

3.2 不同布局計算結果

在三體船型設計與優化過程中，如何對不同片體布局方案的阻力及航行姿態進行快速而準確的預報起著十分關鍵的作用。為了驗證本文方法在片體布局優化過程中的計算精度，對上述三體船型進行片體間距和縱向位置的變化，并將計算結果與試驗結果進行對比分析。計算工況如表1所示。

方案1計算結果與上述對比結果相同，將方案2、方案3計算結果與試驗值對比，結果如圖9～圖11所示。

由圖9可以看出，加入方尾假設后的理論方法計算所得的總阻力在低速段（Fn＜0.35）與試驗值相符良好，而在方尾完全出水的較高傅汝德數階段與試驗值接近，最大誤差發生在傅汝德數0.63左右，約為6%，這主要是由于傅汝德數高于0.6后將產生一定的噴濺及破波會對阻力造成影響。

表1 片體布局方案Tab. 1 Different layout options

圖9 阻力計算結果與試驗值對比Fig. 9 Comparison of resistance calculation results and test values

由圖10和圖11的縱傾、升沉結果可以看出，理論方法對于不同片體布局下航行姿態的預報結果與試驗值在大小和趨勢上相符良好，同樣由于忽略了粘性的影響，計算結果都略小于試驗值。圖11的升沉結果在傅汝德數0.48左右的峰值處，理論計算結果與試驗結果誤差最大，在其余航速范圍內誤差相對較小，這可能時由于在該該航速下片體與主體間興波干擾較強對于升沉產生較大的影響。從圖10和圖11不同片體布局結果的對比可以看出，片體的前移使得縱傾值、升沉值總體是增大的，縱傾值的上升趨勢和峰值位置也發生了一定的變化，而理論與試驗結果所反映出來的趨勢變化時相同的，因此證明理論方法可以很好反映出片體布局變化所造成的深沉和縱傾結果的趨勢變化。

圖12～圖15為傅汝德數0.46時，不同片體布局下的近場興波云圖。

圖10 縱傾計算結果與試驗值對比Fig. 10 Comparison of trim calculation results and test values

圖11 升沉計算結果與試驗值對比Fig. 11 Comparison of heave calculation results and test values

圖12 方案1波形云圖(Fn=0.46)Fig. 12 Scheme 1 waveform nephogram (Fn=0.46)

圖13 方案2波形云圖(Fn=0.46)Fig. 13 Scheme 2 waveform nephogram (Fn=0.46)

圖14 方案3波形云圖(Fn=0.46)Fig. 14 Scheme 3 waveform nephogram (Fn=0.46)

圖15 方案4波形云圖(Fn=0.46)Fig. 15 Scheme 4 waveform nephogram (Fn=0.46)

由圖12～圖15可以看出，與方案1、方案2相比，當側體位于船中時興舶高度有著明顯的增大，而隨著側體寬度增加，方案4所產生的興波又相對方案3有所增加，這與理論計算結果和試驗結果是完全相符的。由方案1與方案2的波形對比可以看出，理論方法很好捕捉到了由于片體橫向位置變化所造成的波形變化，圖13尾部橫向散波數相對于圖12有明顯的增加，這對于片體布局設計階段是十分有意義的。同時，理論方法中對于差分條件、控制點處理等數值方法，可以得到光順的興波波形并很好滿足遠方輻射條件。由于不再使用傳統的方尾虛長度模擬，通過方尾后自由面條件的處理使得方尾后氣穴、興波波形可以自然形成，更加貼近實際情況。

4 結 語

本文基于近場興波理論，探討具有較大方尾的三體船型的阻力及航態預報方法，給出三體船近場興波及阻力、航態計算的算法、網格劃分等數值處理方法。通過與試驗值對比，驗證了本方法的計算精度；通過對不同側體布局船型的算例進行數值模擬，探討和分析了本方法對不同布局三體船型的普適性。可以得出以下結論：

1）本方法對方尾后自由面條件的處理，可以增加理論方法在較小傅汝德數方尾未完全出水情況下對于方尾三體船型阻力及航態數值模擬的精度；

2）理論方法在較大傅汝德數范圍內，對于不同片體布局的方尾三體船型的阻力及航態預報有較高的計算精度，理論預報結果與試驗結果整體相符良好；

3）理論方法計算所得的三體船型周圍流場及近場興波很好反映出了片體布局變化對于興波造成的變化，這對于三體船型設計階段具有實際意義；

4）由于片體布局變化對于粘壓阻力也產生一定的影響，因此在后續工作中需要將這一變化代入總阻力預報中，進而提升計算精度。