周 芃，李 為，陳韜穎

中國艦船研究設(shè)計中心，上海201108

0 引 言

隨著計算機(jī)性能的逐漸提高，用計算流體動力學(xué)（CFD）軟件模擬船體粘性興波流場已成為可能。與試驗(yàn)驗(yàn)證相比，CFD 數(shù)值模擬技術(shù)具有信息量大、成本低、易并行化及能快速響應(yīng)等特點(diǎn)，但其預(yù)報精度一直是船舶CFD 研究者及船舶設(shè)計人員高度關(guān)注的問題。

船舶在靜水中航行時，其表面壓力分布會隨之改變，導(dǎo)致船體表面受力不平衡，故而導(dǎo)致航態(tài)變化。表面壓力主要是船舶在航行中受到的流體對其表面動壓力及粘性力的影響，表面壓力會導(dǎo)致垂向力和縱傾力矩發(fā)生變化，從而引起升沉與縱傾［1］?？梢姶皩?shí)際航行的姿態(tài)與理論計算時的設(shè)計吃水正浮狀態(tài)是有區(qū)別的，尤其是對于高速艦船，其升沉和縱傾與低速狀態(tài)相比更為明顯。因此，在進(jìn)行理論計算時就應(yīng)該考慮航態(tài)的變化。由于船舶形狀左右對稱，因此，只考慮在縱向上的升沉及垂向上的縱傾。

目前，在計及航態(tài)影響的船舶阻力預(yù)報研究方面，國內(nèi)外學(xué)者已進(jìn)行了一些研究。在國外，Dawson［2］早 在1979年就采用Rankine 面元法，并考慮升沉和縱傾的影響研究了興波問題，但這些方法都是建立在弱非線性基礎(chǔ)上，依靠設(shè)計浮態(tài)下的船舶水動力系數(shù)對姿態(tài)進(jìn)行預(yù)報，而非基于船體水動力平衡狀態(tài)下的計算；Yang 等［3］采用非結(jié)構(gòu)化網(wǎng)格和有限元法，對S60 船型計及升沉與縱傾的興波流場進(jìn)行了計算，但未考慮流體粘性的影響，計算結(jié)果與試驗(yàn)結(jié)果尚有一定的差距。在國內(nèi)，張志榮等［4］對水面艦船繞流場進(jìn)行了數(shù)值模擬，并將阻力計算結(jié)果和興波波形與試驗(yàn)結(jié)果進(jìn)行對比，驗(yàn)證了CFD 在水面艦船阻力方面的可行性；曹洪建［5］對航態(tài)變化更大的滑行艇繞流場進(jìn)行了數(shù)值模擬，其航態(tài)變化是通過試驗(yàn)獲得，而非數(shù)值計算得到。近年來，部分學(xué)者將船舶航態(tài)作為求解的一部分對水面艦船的繞流場進(jìn)行數(shù)值計算。董文才等［6］對中高速深V 船型進(jìn)行了計及升沉和縱傾的船舶阻力預(yù)報，其方法是首先計算正浮狀態(tài)的船體升力與縱搖力矩，根據(jù)該升力和縱搖力矩與靜浮力和縱傾矩的平衡來調(diào)整船舶航態(tài)，然后根據(jù)調(diào)整后的航態(tài)再次建模并進(jìn)行數(shù)值計算，如此重復(fù)多次得到穩(wěn)定、平衡的船舶航態(tài)，進(jìn)而得到船舶的阻力，其計算阻力值與試驗(yàn)值的誤差均在5%以內(nèi)。倪崇本等［7］采用動網(wǎng)格技術(shù)對高速多體船進(jìn)行了計及航態(tài)的阻力預(yù)報，該方法無需多次建模和重復(fù)求解，計算更加簡便。

動網(wǎng)格的一個關(guān)鍵技術(shù)是邊界運(yùn)動更新時，網(wǎng)格也隨之更新，并且不能出現(xiàn)負(fù)體積網(wǎng)格。但由于船舶表面形狀較復(fù)雜，目前的動網(wǎng)格技術(shù)還不能滿足計算要求。本文的研究是基于靜態(tài)平衡法，但在該方法的基礎(chǔ)上作出假設(shè)進(jìn)行了改進(jìn)，無需多次建模，只需進(jìn)行迭代計算即可，大大降低了工作量。通過數(shù)值計算，得到了傅氏數(shù)在0.176～0.441 范圍內(nèi)的船模阻力，并在傅氏數(shù)大于0.35 時計及了航態(tài)迭代計算，其與船模試驗(yàn)值的誤差均在5%以內(nèi)。另外，通過分析不同傅氏數(shù)下的自由表面波形圖和波形等值線圖，分析了其艏艉興波隨航速的變化情況，符合試驗(yàn)驗(yàn)證結(jié)果。

1 流場數(shù)值模擬

1.1 控制方程與湍流模型

船舶繞流場屬高度復(fù)雜的三維流動，其主要原因是船首尾形狀復(fù)雜、曲率變化大等，這一復(fù)雜的特性主要體現(xiàn)在船的艉部伴流中。在數(shù)值計算的準(zhǔn)確程度方面，湍流模型的選取會直接影響其精度［8-9］。本文將采用商業(yè)軟件CFX 中RNG k-ε湍流模型進(jìn)行模擬，湍流模型及控制方程［10］如下。

在直角坐標(biāo)系下，不可壓縮牛頓流體連續(xù)性方程與RANS 方程為

式中：ρ 為密度；μ 為流體粘性系數(shù)；p 為平均壓力；Fi為外力項(xiàng)；ui為時均速度；u'i 為脈動速度；與脈動速度相關(guān)的項(xiàng)稱為雷諾應(yīng)力。

文中采用的湍流模型是RNG k-ε模型，這種模型是標(biāo)準(zhǔn)k-ε模型的改進(jìn)形式，它修正了湍動粘度，考慮了平均流動中的旋轉(zhuǎn)及旋流流動情況，增加了反映主流的時均應(yīng)變率。該模型是一種適合船舶流場計算的湍流模型，其方程形式如下。

湍流動能方程與耗散率方程［11］：

1.2 浮態(tài)調(diào)整方程

將船舶看成是一個剛體，其在靜水航行穩(wěn)定時，必須滿足力和力矩平衡。

式中：F 為外力；M 為外力對重心處的力矩。以上式（5）和式（6）稱為浮態(tài)平衡方程。

本文采取的是一種簡化形式的迭代計算。首先，假設(shè)船舶在升沉?xí)r水線面面積保持設(shè)計水線面面積不變；其次，浮態(tài)調(diào)整運(yùn)動十分緩慢。

基于以上兩點(diǎn)假設(shè)，升沉和縱傾的公式可由以下等式給出：

式中：Fz為船體所受壓力在z 方向（垂向）的分量；m 為船體的質(zhì)量；Δd 為船體的升沉值（下沉為正）；S 為設(shè)計水線面面積；M 為壓力對船體重心處的力矩在y 方向（橫向）的分量；Δ 為船體排水量；為船體的縱穩(wěn)性高；θ 為船體縱傾角（以艉傾為正）。根據(jù)計算算出的垂向力和力矩驗(yàn)證平衡方程，如果不滿足平衡方程（5）和方程（6），便帶入式（7）和式（8），算出浮態(tài)調(diào)整量，改變船舶的浮態(tài)，重新進(jìn)行計算。

1.3 計算模型與網(wǎng)格劃分

1.3.1 計算模型

本文以某高速艦船為研究對象，CFD 計算采用1∶19.756 8 的縮小模型，模型建立至設(shè)計水線上的一定高度處。模型主要參數(shù)如表1 所示。

表1 船模主要參數(shù)Tab.1 Principle parameters of ship model

模型的建立在三維建模軟件UG NX7.5 中進(jìn)行。將各橫剖線、水線聯(lián)合形成船艏、船舯和船艉曲面，其三維模型如圖1 所示。

圖1 三維模型Fig.1 The 3D model of computation object

1.3.2 計算域及邊界條件

由于計算模型是沿中縱剖面完全對稱，故為便于計算，僅取一半模型進(jìn)行幾何建模和網(wǎng)格劃分，另一半做鏡像處理。為了使計算結(jié)果盡量不受邊界條件的影響，將外流域場設(shè)為長方體，具體設(shè)置如下：

1）水域：將船艏向上游延伸至1 倍船長處設(shè)為水域速度入口；船艉向下游延伸至3 倍船長處設(shè)為水域壓力出口；距中縱剖面5 倍半寬值處為水域外邊界，設(shè)為開放式邊界；設(shè)計水線面向下延伸5 倍半寬處為水域底部，設(shè)為滑移壁面。其速度大小為來流速度，方向指向船艉。

2）空氣域：在設(shè)計水線面上增加一個深度為艦船吃水高度1.5 倍的空氣層。邊界條件與水域完全類似。

3）船表面：設(shè)為不可滑移壁面邊界。

1.3.3 網(wǎng)格劃分

計算網(wǎng)格的質(zhì)量會直接影響到計算的可行性、收斂性以及其計算精度。本文采用的是計算精度更高、收斂性更好的全船結(jié)構(gòu)化網(wǎng)格，總數(shù)約156 萬。計算網(wǎng)格在ICEM 中生成，采用O-H 多塊結(jié)構(gòu)化網(wǎng)格，即縱向?yàn)镠 型網(wǎng)格，橫向?yàn)镺 型網(wǎng)格。為了更精確地計算摩擦阻力，在船體表面向外生成一層O-grid 劃分邊界層網(wǎng)格，使y+控制在30 左右。在網(wǎng)格劃分過程中，保證靠近船體表面的網(wǎng)格密，遠(yuǎn)離船體表面的逐漸變疏；在船體表面，保證艏、艉部的網(wǎng)格與中部和頂部相比更密。這樣，就能保證在控制網(wǎng)格數(shù)量的同時達(dá)到理想的計算精度。圖2 所示為計算域網(wǎng)格劃分情況，圖3～圖5 所示為局部網(wǎng)格劃分情況。

圖2 計算域網(wǎng)格Fig.2 Grid system of computational region

圖3 艏部網(wǎng)格Fig.3 The grid of foreship

圖4 艉部網(wǎng)格Fig.4 The grid of stern

圖5 船艏的對稱面及邊界層網(wǎng)格Fig.5 The boundary grid of symmetry at foreship

1.4 計算結(jié)果及分析

1.4.1 阻力計算結(jié)果及分析

1）不考慮航態(tài)變化。

按照上述網(wǎng)格劃分方法及邊界條件設(shè)置，將網(wǎng)格導(dǎo)入計算流體動力學(xué)軟件CFX 中進(jìn)行氣、液兩相流數(shù)值計算，湍流模型采用RNG k-ε模型。船模試驗(yàn)結(jié)果與數(shù)值結(jié)果按照二因次法進(jìn)行阻力換算，采用柏—許公式計算相當(dāng)平板的摩擦阻力系數(shù)（Cf=0.455/（logRe）2.58，其中Re 為雷諾數(shù)）。

取ρ = 998 kg/m3，ρ海水= 1 025.91 kg/m3，ν =0.892 92×10-6m-2·s-1，ν海水= 1.188 31×10-6m-2·s-1，計算得到的試驗(yàn)值各速度下的阻力如表2 所示。

表2 不考慮航態(tài)下阻力計算值與試驗(yàn)值的比較Tab.2 The comparison of resistances from computation results and experiment data in fix pose

根據(jù)上述計算結(jié)果可知，對于高速艦船的阻力，在傅氏數(shù)低于0.35 左右時，由于航態(tài)變化較小，對阻力計算影響的誤差不大，因此誤差在可接受范圍內(nèi)。但隨著航速的增加，計算得出的阻力誤差在傅氏數(shù)大于0.4 時甚至?xí)^10%。這是因?yàn)殡S著航速的增加，船體的升沉和縱傾變化越來越明顯，而計算時默認(rèn)的設(shè)計吃水正浮狀態(tài)與實(shí)際的航行姿態(tài)差距較大，從而導(dǎo)致船體興波和船體表面的壓力分布與真實(shí)狀態(tài)不相符合，因此計算達(dá)不到工程應(yīng)用的精度。

2）在傅氏數(shù)大于0.35 時考慮航態(tài)變化。

根據(jù)正浮狀態(tài)各速度下的理論計算結(jié)果，可分別求解船體表面力和力矩的積分，并將其帶入力和力矩平衡方程進(jìn)行驗(yàn)證，若不滿足，便將其帶入式（7）、式（8）進(jìn)行升沉和縱傾值的調(diào)節(jié)，然后再改變航態(tài)重新進(jìn)行計算。

表3 所示為浮態(tài)調(diào)整計算的基本參數(shù)。以基線與中橫剖面的交點(diǎn)為原點(diǎn)，坐標(biāo)系定義如下：x軸為中縱剖面與基平面的交線，以指向船首為正；y 軸為中橫剖面與基平面的交線方向，以指向右舷為正；z 軸按右手法則確定。

表3 浮態(tài)調(diào)整計算基本參數(shù)Tab.3 Basic parameters of pose adjustment

靜態(tài)平衡法的人工實(shí)現(xiàn)：根據(jù)正浮狀態(tài)下的理論計算結(jié)果，分別求解船體表面壓力沿z 方向（垂向力）的分量，以及壓力在重心處的力矩沿y軸的分量。計算得到調(diào)整的升沉值和縱傾值后，在ICEM 中進(jìn)行模型與塊的平移和旋轉(zhuǎn)，重新進(jìn)行計算。此方法無需重新建模和網(wǎng)格劃分，大大減少了工作量。由此，便可得到新的垂向力及壓力在重心處的力矩沿y 軸的分量。重新調(diào)整航態(tài)計算，直至垂向力與重力之差小于重力的3%，便可認(rèn)為壓力與重力平衡了，無需繼續(xù)調(diào)整。一般經(jīng)過1～2 次迭代即可達(dá)到平衡。調(diào)整航態(tài)后的阻力計算結(jié)果與試驗(yàn)值的對比如表4 和圖6、圖7所示。

圖6 航態(tài)變化前后與試驗(yàn)值剩余阻力比較Fig.6 The comparison of residual resistances

表4 航態(tài)調(diào)整后阻力計算值與試驗(yàn)值比較Tab.4 The comparison of resistances from computation results and experiment data in free pose

圖7 航態(tài)變化前后與試驗(yàn)值總阻力比較Fig.7 The comparison of total resistances

從表4、圖6 和圖7 可以看出，計及升沉和縱傾后，計算值與試驗(yàn)值的誤差均在5%以內(nèi)，吻合度較好。同時也可以看出，航態(tài)的變化使摩擦阻力更接近試驗(yàn)換算值，其變化對摩擦阻力影響較小，而對興波阻力影響則較大；低速時航態(tài)的變化對總阻力影響較小，而高速時對總阻力的影響則較大。從阻力的對比可以看出，與試驗(yàn)換算出的相當(dāng)平板摩擦阻力相比，計算出的摩擦阻力稍偏小，這可能是由船體曲率變化所致，同時也說明邊界層O 網(wǎng)格的劃分還有待改進(jìn)。但摩擦阻力計算與試驗(yàn)的誤差較小，說明此計算方法對于摩擦阻力的計算具有一定的精度。總阻力在大于18 kn以后，與試驗(yàn)值相比均偏小，且隨著航速的增加差距逐漸增大，但精度基本控制在5%以內(nèi)，對于工程應(yīng)用能起到一定的指導(dǎo)作用。

由表5 可以看出，艉傾對應(yīng)的是一個最低速度，在這之前是艏傾，之后，艉傾便隨著速度的增大而增大。升沉在一定速度范圍內(nèi)也是隨速度的增大而增大。

1.4.2 壓力及波形計算結(jié)果分析

圖8 給出了船體表面壓力分布云圖在固定姿態(tài)和計及航態(tài)變化的比較（左圖為固定航態(tài)，右圖計及了升沉和縱傾）。從中可看出：隨著在航行過程中出現(xiàn)艉傾，艏部球鼻艏處的吃水變淺，最大壓力值降低，艏部至船舯部之間的區(qū)域壓力均增大，船體艉部處的壓力也稍增大，但增加值不及艏部明顯。在考慮了升沉與縱傾的影響后，船體表面壓力分布的準(zhǔn)確性能更好地反映船體姿態(tài)變化的影響，從而使阻力預(yù)報值更準(zhǔn)確。

表5 各速度下的升沉、縱傾值Tab.5 The sinkage and trim results at different speeds

圖8 船體壓力對比云圖Fig.8 The comparison of hull pressure contours

圖9 所示分別為實(shí)船航速為24，27，30 kn 時的船模試驗(yàn)波面效果圖。由圖中可以看出波面效果逼真，峰谷清晰，艏肩波及艉部雞尾流均能捕捉到。

圖9 波面效果圖Fig.9 The effect drawings of wave

圖10 所示為實(shí)船航速分別為24，27，30 kn 時的波高對比云圖。分析船模試驗(yàn)波面效果圖和波高對比云圖（左圖為計及升沉和縱傾，右圖為固定航態(tài)）可知，計及航態(tài)后，艏、艉興波高度增加，這也與計及航態(tài)后算得的興波阻力比不計航態(tài)算得的興波阻力更大的規(guī)律相符。艏部波高隨著航速的增加逐漸升高，并且艏肩波有逐漸向船舯后移的現(xiàn)象。從船艉興波波高對比圖可以看出，船體在低速時不會產(chǎn)生雞尾流，只有在高速的情況下才會產(chǎn)生，并且在一定速度范圍內(nèi)，隨著速度的升高雞尾流的長度逐漸變長。

圖10 波高對比云圖Fig.10 The comparison of wave height contours

2 結(jié) 論

通過對計及航態(tài)的艦船阻力及自由面興波的預(yù)報，可以得出以下結(jié)論：

1）在低速狀態(tài)下，航行狀態(tài)的變化對船舶水動力性能影響較??；而當(dāng)航速增加到一定程度后，其阻力計算誤差會達(dá)到10%以上。在計及航態(tài)的情況下，阻力計算誤差降低到了5%以內(nèi)，由此可知，航態(tài)對船舶水動力性能的影響不可忽略，此時，計及航態(tài)是必要的。通過數(shù)值計算的摩擦阻力、總阻力結(jié)果與船模試驗(yàn)值換算結(jié)果的比較，表明其符合程度令人滿意，驗(yàn)證了本文數(shù)值計算方法的正確性和有效性。

2）本文的方法基于靜態(tài)平衡法，在其方法的實(shí)現(xiàn)上進(jìn)行了改進(jìn)。每個航速下，改變航態(tài)后通常能節(jié)省1～2 次塊和網(wǎng)格劃分過程工作量，在保證計算精度的同時，大大縮短了工作時間，對實(shí)際工程應(yīng)用具有積極的指導(dǎo)作用。

3）航態(tài)的變化對摩擦阻力影響較小，對興波阻力影響較大。同時，在低速時，航態(tài)的變化對總阻力影響較小，在高速時影響較大。

4）在一定速度范圍內(nèi)，隨著航速的增加，升沉值逐漸增大。艏傾和艉傾存在一個速度分界點(diǎn)，在此分界點(diǎn)后，隨著速度的增大，艉傾逐漸變大。對于本文的研究，此分界點(diǎn)的傅氏數(shù)約為0.37。

5）艏部波高隨著航速的增加逐漸升高，并且艏肩波有逐漸向船舯后移的現(xiàn)象。船體在低速時不會產(chǎn)生雞尾流，只有在高速的情況下才會產(chǎn)生雞尾流，并且在一定速度范圍內(nèi)，隨著速度的升高，雞尾流的長度逐漸變長。對于本文的研究，當(dāng)傅氏數(shù)達(dá)到0.36 左右時會產(chǎn)生雞尾流。

［1］姚朝幫. 計及航態(tài)影響的船舶阻力計算方法及其應(yīng)用研究［D］.武漢：海軍工程大學(xué)，2010.YAO Chaobang. Study on methods to calculate resis?tance of a ship taking the effect of dynamic sinkage and trim［D］. Wuhan：Naval University of Engineering，2010.

［2］DAWSON C. Calculations with the XYZ free surface program for five ship models［C］//Proceeding of Work?shop on Ship Wave-Resistance Computations，1979：232-255.

［3］YANG C，L?HNER R. Calculation of ship sinkage and trim using a finite element and unstructured grids［J］. International Journal of Computational Fluid Dy?namics，2002，16（3）：217-227.

［4］ZHANG Zhirong，ZHAO Feng，LI Baiqi. Numerical calculation of viscous free surface flow about ship hull［J］.Journal of Ship Mechanics，2002，6（6）：11-18.

［5］曹洪建. 基于FLUENT 的滑行艇阻力計算研究［D］.哈爾濱：哈爾濱工程大學(xué)，2008.CAO Hongjian. The computation and research on resis?tance of planing craft based on the software FLUENT［D］.Harbin：Harbin Engineering University，2008.

［6］董文才，姚朝幫.中高速深V 型船阻力預(yù)報方法及尾板減阻機(jī)理［J］. 哈爾濱工程大學(xué)學(xué)報，2011，32（7）：848-852.DONG Wencai，YAO Chaobang. Study on resistance prediction method and resistance reduction mechanism of medium & high speed deep-Vee ship by stern flap［J］. Journal of Harbin Engineering University，2011，32（7）：848-852.

［7］倪崇本，朱仁傳，繆國平，等. 計及航行姿態(tài)變化的高速多體船阻力預(yù)報［J］. 水動力學(xué)研究與進(jìn)展（A輯），2011，26（1）：101-107.NI Chongben，ZHU Renchuan，MIAO Guoping，et al.The resistance prediction for high speed multi-hull ves?sels with consideration of hull gesture variation during voyage［J］. Chinese Journal of Hydrodynamics（Ser.A），2011，26（1）：101-107.

［8］張志榮，李百齊，趙峰.船舶粘性流動計算中湍流模式應(yīng)用的比較［J］. 水動力學(xué)研究與進(jìn)展（A 輯）：2004，19（5）：637-642.ZHANG Zhirong，LI Baiqi，ZHAO Feng. The compari?son of turbulence models applied to viscous ship flow computation［J］. Chinese Journal of Hydrodynamics（Ser.A），2004，19（5）：637-642.

［9］徐杰，俞汲，黃少鋒，等.散貨船阻力預(yù)報中湍流模型的應(yīng)用［J］.中國造船，2011，52（增刊1）：53-58.XU Jie，YU Ji，HUAGN Shaofeng，et al. Application of Turbulence models to prediction of powering performa?ce for ship［J］. Shipbuilding of China，2011，52（s1）：53-58.

［10］張師帥.計算流體動力學(xué)及其應(yīng)用—CFD 軟件的原理與應(yīng)用［M］.武漢：華中科技大學(xué)出版社，2010.

［11］王福軍.計算流體動力學(xué)分析—CFD 軟件原理與應(yīng)用［M］.北京：清華大學(xué)出版社，2004.