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

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

0 引 言

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

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

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

1 數(shù)學(xué)模型

1.1 基本模型

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

圖1 坐標(biāo)系Fig. 1 Coordinate system

假設(shè)流體理想，流動(dòng)定常，無(wú)限水深，船舶興波為微幅波。流場(chǎng)速度勢(shì) ψ分解為無(wú)自由面影響的重疊模速度勢(shì) φ和考慮自由液面作用的擾動(dòng)速度勢(shì)φ，即

由于重疊模速度勢(shì)中包含了來(lái)流速度勢(shì)，因此假設(shè)，φ＜＜φ。流場(chǎng)速度勢(shì) ψ對(duì)應(yīng)的波高為ζ，重疊模速度勢(shì) φ在z=0 處所產(chǎn)生切向速度對(duì)應(yīng)的波高為ξ，擾動(dòng)速度勢(shì) φ對(duì)應(yīng)的波高為η，且滿足以下關(guān)系式：

流場(chǎng)中有

物面上滿足

速度勢(shì) ψ在自由面上分別滿足運(yùn)動(dòng)學(xué)條件和動(dòng)力學(xué)條件如下：

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

速度勢(shì)滿足無(wú)窮遠(yuǎn)處趨于零的遠(yuǎn)方輻射條件。

1.2 自由面邊界條件

考慮到φ關(guān)于靜水面對(duì)稱，則有：

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

2 數(shù)值計(jì)算

2.1 數(shù)值處理

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

p

離散單元布置源匯后，物面控制點(diǎn) 重疊模流場(chǎng)速度勢(shì)勢(shì)表示為：

圖2 自由面網(wǎng)格Fig. 2 Free surface mesh

其中：Sb為物面；σB0為物面控制點(diǎn)處源強(qiáng)；rpq為任意源點(diǎn)q到控制點(diǎn)p的距離。

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

擾動(dòng)速度勢(shì)表示為：

為了保證自由液面不發(fā)生震蕩以得到更加光順、貼合實(shí)際的自由面波形，進(jìn)而使計(jì)算結(jié)果更加精確，在計(jì)算時(shí)對(duì)自由面單元的控制點(diǎn)進(jìn)行升高，升高高度與波長(zhǎng)和各個(gè)單元各自的長(zhǎng)度成一定比例。同時(shí)，為了更好滿足遠(yuǎn)方輻射條件，使得向前差分不受單元對(duì)控制中心切向誘導(dǎo)速度為零對(duì)數(shù)值計(jì)算造成的影響，將控制點(diǎn)統(tǒng)一前移一定距離，移動(dòng)距離為單元長(zhǎng)度的0.1倍。

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

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

其中：ηdry為方尾中心與兩側(cè)未浸水長(zhǎng)度系數(shù)；Bt，Dt為方尾的寬度與吃水；Ft為方尾傅汝德數(shù)；RNT為方尾雷諾數(shù)。參照文獻(xiàn)數(shù)據(jù)對(duì)系數(shù)C1，C2，C3，C4進(jìn)行取值，求得方尾中心與兩側(cè)的吃水后，對(duì)方位進(jìn)行靜水力修正。

2.3 興波阻力及航態(tài)計(jì)算

船體表面控制點(diǎn)壓力系數(shù)為：

其中，Zi為控制點(diǎn)到靜水面間距離。

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

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

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

僅對(duì)船體水線以下近水面處局部網(wǎng)格進(jìn)行實(shí)時(shí)劃分，而大部分網(wǎng)格始終保持不變，這樣既減小了網(wǎng)格重新離散帶來(lái)的誤差，也節(jié)省了計(jì)算時(shí)間。通過(guò)反復(fù)迭代直到滿足收斂條件（Fz＜ε，Ny＜ε）即可得到考慮船體姿態(tài)變化影響的興波阻力，進(jìn)行船體姿態(tài)預(yù)報(bào)。

3 算例及結(jié)果分析

3.1 數(shù)值計(jì)算及驗(yàn)證

以具有較大方尾的某三體船為例，對(duì)該船型的興波阻力、航行姿態(tài)以及方尾吃水進(jìn)行理論預(yù)報(bào)并與相應(yīng)的試驗(yàn)結(jié)果對(duì)比。在總阻力計(jì)算時(shí)，參照文獻(xiàn)[12]取粘壓阻力系數(shù)為0.1。

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

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

主體縱向網(wǎng)格60個(gè)，垂向網(wǎng)格10個(gè)；側(cè)體縱向網(wǎng)格30個(gè)，垂向網(wǎng)格5個(gè)。自由表面網(wǎng)格采用隨船長(zhǎng)傅汝德數(shù)變化而變化，單個(gè)波長(zhǎng)保持25個(gè)網(wǎng)格，船首、尾及橫向網(wǎng)格保持1.05的擴(kuò)張比例，動(dòng)態(tài)保證計(jì)算域船首往前1L，船尾往后1.5L，橫向1L。船體網(wǎng)格如圖4所示，在Fn=0.3時(shí)自由面網(wǎng)格如圖5所示。

圖4 船體表面網(wǎng)格示意圖Fig. 4 Schematic diagram of hull surface grid

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

將該船型興波阻力計(jì)算結(jié)果與試驗(yàn)結(jié)果進(jìn)行對(duì)比，升沉、縱傾、阻力計(jì)算結(jié)果如圖6～圖8所示。

圖6 升沉計(jì)算結(jié)果與試驗(yàn)值對(duì)比Fig. 6 Comparison of heave calculation results and test values

圖7 縱傾計(jì)算結(jié)果與試驗(yàn)值對(duì)比Fig. 7 Comparison of trim calculation results and test values

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

3.2 不同布局計(jì)算結(jié)果

在三體船型設(shè)計(jì)與優(yōu)化過(guò)程中，如何對(duì)不同片體布局方案的阻力及航行姿態(tài)進(jìn)行快速而準(zhǔn)確的預(yù)報(bào)起著十分關(guān)鍵的作用。為了驗(yàn)證本文方法在片體布局優(yōu)化過(guò)程中的計(jì)算精度，對(duì)上述三體船型進(jìn)行片體間距和縱向位置的變化，并將計(jì)算結(jié)果與試驗(yàn)結(jié)果進(jìn)行對(duì)比分析。計(jì)算工況如表1所示。

方案1計(jì)算結(jié)果與上述對(duì)比結(jié)果相同，將方案2、方案3計(jì)算結(jié)果與試驗(yàn)值對(duì)比，結(jié)果如圖9～圖11所示。

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

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

圖9 阻力計(jì)算結(jié)果與試驗(yàn)值對(duì)比Fig. 9 Comparison of resistance calculation results and test values

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

圖12～圖15為傅汝德數(shù)0.46時(shí)，不同片體布局下的近場(chǎng)興波云圖。

圖10 縱傾計(jì)算結(jié)果與試驗(yàn)值對(duì)比Fig. 10 Comparison of trim calculation results and test values

圖11 升沉計(jì)算結(jié)果與試驗(yàn)值對(duì)比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相比，當(dāng)側(cè)體位于船中時(shí)興舶高度有著明顯的增大，而隨著側(cè)體寬度增加，方案4所產(chǎn)生的興波又相對(duì)方案3有所增加，這與理論計(jì)算結(jié)果和試驗(yàn)結(jié)果是完全相符的。由方案1與方案2的波形對(duì)比可以看出，理論方法很好捕捉到了由于片體橫向位置變化所造成的波形變化，圖13尾部橫向散波數(shù)相對(duì)于圖12有明顯的增加，這對(duì)于片體布局設(shè)計(jì)階段是十分有意義的。同時(shí)，理論方法中對(duì)于差分條件、控制點(diǎn)處理等數(shù)值方法，可以得到光順的興波波形并很好滿足遠(yuǎn)方輻射條件。由于不再使用傳統(tǒng)的方尾虛長(zhǎng)度模擬，通過(guò)方尾后自由面條件的處理使得方尾后氣穴、興波波形可以自然形成，更加貼近實(shí)際情況。

4 結(jié) 語(yǔ)

本文基于近場(chǎng)興波理論，探討具有較大方尾的三體船型的阻力及航態(tài)預(yù)報(bào)方法，給出三體船近場(chǎng)興波及阻力、航態(tài)計(jì)算的算法、網(wǎng)格劃分等數(shù)值處理方法。通過(guò)與試驗(yàn)值對(duì)比，驗(yàn)證了本方法的計(jì)算精度；通過(guò)對(duì)不同側(cè)體布局船型的算例進(jìn)行數(shù)值模擬，探討和分析了本方法對(duì)不同布局三體船型的普適性。可以得出以下結(jié)論：

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

2）理論方法在較大傅汝德數(shù)范圍內(nèi)，對(duì)于不同片體布局的方尾三體船型的阻力及航態(tài)預(yù)報(bào)有較高的計(jì)算精度，理論預(yù)報(bào)結(jié)果與試驗(yàn)結(jié)果整體相符良好；

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

4）由于片體布局變化對(duì)于粘壓阻力也產(chǎn)生一定的影響，因此在后續(xù)工作中需要將這一變化代入總阻力預(yù)報(bào)中，進(jìn)而提升計(jì)算精度。