盧錦國 梁中剛 吳方良 周軼美

中國艦船研究設(shè)計中心，湖北武漢 430064

水下航行體回轉(zhuǎn)水動力數(shù)值計算研究

盧錦國 梁中剛 吳方良 周軼美

中國艦船研究設(shè)計中心，湖北武漢 430064

準(zhǔn)確預(yù)報水動力對水下航行體的安全操縱設(shè)計而言至關(guān)重要?；谏逃密浖﨔LUENT數(shù)值求解RANS方程，采用定常旋轉(zhuǎn)坐標(biāo)系，運用相對運動理論及運動疊加原理，并對控制方程中向心力源項空間離散誤差進行修正，預(yù)報水下航行體做單平面回轉(zhuǎn)運動的受力及力矩。在Re＝11.7×106條件下，力和力矩預(yù)報精度的偏差在12%以內(nèi)，表明所采用的數(shù)值預(yù)報方法有效、可行，具有較好的工程實用價值。

回轉(zhuǎn)運動；水動力；向心力修正；FLUENT；水下航行體

1 引言

操縱性能預(yù)報是水下航行體設(shè)計的重要環(huán)節(jié)。數(shù)值方法應(yīng)用于操縱性能預(yù)報由來已久，針對不同的模型及不同的航行姿態(tài)，國內(nèi)外眾多研究人員開展了這方面的研究。公開發(fā)表的文獻大多針對直航或偏航狀態(tài)下水動力性能進行研究，對比試驗數(shù)據(jù)來源于直線拖曳水池的試驗結(jié)果。水下航行體旋臂水池試驗通過在一定范圍內(nèi)變換半徑、攻角及舵角來確定旋轉(zhuǎn)導(dǎo)數(shù)、非線性及耦合水動力導(dǎo)數(shù)及力矩導(dǎo)數(shù)。試驗研究表明，在較大半徑旋臂水池試驗條件下得到旋轉(zhuǎn)導(dǎo)數(shù)，通常情況下比相應(yīng)的PMM試驗得到的旋轉(zhuǎn)導(dǎo)數(shù)具有更高的精度［1］。然而，在對旋臂水池試驗的數(shù)值模擬研究方面，由于流動復(fù)雜及試驗數(shù)據(jù)來源少等原因，已經(jīng)開展的工作較為有限，針對這方面的研究主要由美國泰勒水池研究人員發(fā)起。Sung等［2］融合了多塊網(wǎng)格、多重網(wǎng)格、局部加密等技術(shù)求解三維不可壓縮RANS方程，采用標(biāo)準(zhǔn)k－ω及修正的B－L湍流模型，數(shù)值模擬了回轉(zhuǎn)主體在不同偏航角條件下的水下定?；剞D(zhuǎn)運動，作者對影響數(shù)值計算結(jié)果的條件進行了系統(tǒng)研究，求得的垂向力及縱傾力矩與試驗值的誤差在20%以內(nèi)。2002年，Sung 等［3］人采用 standard k－ ω 及 realizable k － ε湍流模型，對ONR Body－1潛艇模型進行回轉(zhuǎn)運動模擬，并與試驗數(shù)據(jù)進行驗證，對比較結(jié)果進行了合理分析。隨著旋臂水池試驗報告的逐步公開［4－6］，模擬水下航行體回轉(zhuǎn)運動的研究將漸漸成為新的研究熱點。

本文基于對旋臂水池試驗的數(shù)值模擬，采用商用軟件FLUENT，首先獲得穩(wěn)定的純旋轉(zhuǎn)運動流場，提出控制方程中向心力源項空間離散誤差修正方法，進而預(yù)報某型水下航行體做單平面回轉(zhuǎn)運動水動力性能，通過與已有試驗數(shù)據(jù)進行比較，對影響數(shù)值計算結(jié)果的條件進行優(yōu)選，最終建立起一套可行的回轉(zhuǎn)水動力數(shù)值預(yù)報方法。

2 數(shù)值計算方法

2.1 物理描述

對于固定于大地的慣性坐標(biāo)系而言，由于水下航行體的定常回轉(zhuǎn)運動存在向心加速度，因此該運動在此坐標(biāo)系下為非定常問題，將其轉(zhuǎn)化為旋轉(zhuǎn)運動坐標(biāo)系下進行求解，可以使問題得到簡化，如圖1所示。

圖1 水下航行體回轉(zhuǎn)運動示意圖Fig.1 Schematic of the submerge d vehicle in turningmotion

當(dāng)運動方程在旋轉(zhuǎn)系中求解時，流體的加速度（向心加速度及科氏加速度）將作為源項出現(xiàn)在動量方程中。

2.2 控制方程

在角速度一定的旋轉(zhuǎn)運動坐標(biāo)系下的連續(xù)性方程及不可壓縮RANS方程為：

式中，xj與uj分別為笛卡爾坐標(biāo)系下的坐標(biāo)及速度分量；δij為克羅內(nèi)克符號記法；ν為運動學(xué)粘性系數(shù)；τij為雷諾應(yīng)力張量。

2.3 湍流模型

為了實現(xiàn)對回轉(zhuǎn)運動水流橫向分離復(fù)雜流動的準(zhǔn)確捕捉，結(jié)合文獻［7］給出的水下航行體在偏航狀態(tài)下力和力矩的數(shù)值預(yù)報研究成果，可知SST k－ω湍流模型對帶攻角來流條件下的垂向力和縱傾力矩預(yù)報具有較好的預(yù)報精度。本文采用SST k－ω湍流模型進行數(shù)值模擬研究。

2.4 模型及計算域

2.4.1 模型描述

旋臂試驗?zāi)Ｐ蜑榱骶€型回轉(zhuǎn)體的58系列模型，編號4621，采用6階多項式數(shù)學(xué)表達式描述。體長為4.572m、最大直徑為0.623m、排水體積0.836 2 m3，形心位置位于回轉(zhuǎn)軸線上，距首部2.037 3 m。進行模型試驗時，采用側(cè)裝模型的方法，繞OY軸轉(zhuǎn)動，取角速度Ωy＜0，重心線速度大小為2.572 m／s，以回轉(zhuǎn)體長為特征長度，相應(yīng)的雷諾數(shù)大小為Re＝11.7×106。

圖2 計算模型示意圖Fig.2 Schematic of the comput ation model

2.4.2 計算域選取

計算域的選取需要根據(jù)水下航行體的主尺度和回轉(zhuǎn)半徑來確定［2－3］。選取來流截面距航行體前端約1倍體長，外周邊界距航行體中軸線約1倍體長，出流截面處距航行體尾端約3倍體長，如圖3所示。

2.5 邊界條件

在FLUENT軟件中，復(fù)雜的數(shù)值模擬邊界條件可以通過用戶自定義函數(shù) （UDF）功能編程實現(xiàn)。

圖3 計算域示意圖Fig.3 Schematic of the comput ation domain

水下航行體模型表面設(shè)置為無滑移壁面邊界條件（u＝v＝w＝0），在旋轉(zhuǎn)運動坐標(biāo)系下，將此條件代入式（1），并忽略粘性影響，可以得到模型表面外法線方向應(yīng)滿足壓力梯度為0條件，即

在流場中未放置模型的條件下，由描述流體運動的歐拉方程可知，外周邊界的速度分布在某一固定時刻，同一半徑曲線上任一點的速度方向與曲線在該點的切線方向重合，即外周邊界為由一系列沿同一半徑的流線組成的流管。因此，在得到穩(wěn)定數(shù)值解的條件下，可以得到穩(wěn)定的純旋轉(zhuǎn)運動流場。

2.6 網(wǎng)格劃分及數(shù)值求解

2.6.1 網(wǎng)格劃分

本文使用ICEM－CFD軟件對流域進行多塊結(jié)構(gòu)化網(wǎng)格劃分。對航行體邊界層網(wǎng)格采用O型網(wǎng)格劃分，主要繞流區(qū)域采用C型網(wǎng)格劃分，外部流域采用H型網(wǎng)格劃分，如圖4所示。

圖4 網(wǎng)格劃分示意圖Fig.4 Schematic of the grid partition

邊界層網(wǎng)格考慮各湍流模型及壁面函數(shù)對粘性層的不同要求，參考FLUENT手冊對y＋進行合理取值，y＋＜5，50%以上分布于［0，1］之間；層間增長率均設(shè)置為1.15。沿縱向及周向的網(wǎng)格布置考慮背流面、首部及尾段區(qū)域流動的復(fù)雜性，在這些區(qū)域劃分較細密的網(wǎng)格節(jié)點；同時也在這些區(qū)域采用局部網(wǎng)格加密技術(shù)。

2.6.2 數(shù)值求解

流體運動的控制方程及湍流運動方程采用有限體積法進行離散，壓力項采用標(biāo)準(zhǔn)差分格式，速度、動量及湍流參量采用2階迎風(fēng)差分格式。應(yīng)用Simplec法處理壓力與速度耦合問題，離散方程運用Gauss－Seidel迭代方法求解，并以MutilGrid技術(shù)加速計算收斂。

3 計算結(jié)果與對比分析

3.1 穩(wěn)定的旋轉(zhuǎn)流場

首先驗證上述邊界條件下求解三維不可壓縮RANS方程可以獲得穩(wěn)定的旋轉(zhuǎn)流場。本文采用環(huán)半徑為10m、圓半徑為5m的1／4圓環(huán)計算模型。比較了高度方向Z＝0 m平面內(nèi)，沿半徑方向R＝8、9、10、11、12m 的速度分布， 結(jié)果如表 1 所示。Z＝0截面各半徑處速度分布的不均勻度（（最大值－最小值）／平均值）小于1%，數(shù)值耗散、截斷誤差等可能是造成這一差異的原因。圖5所示為Z＝0截面速度隨半徑的分布曲線。

表1 Z＝0截面各半徑處速度分布的不均勻度分析Tab.1 Nonuniform analysis of velocity distribution on the plane Z＝0 at different radii

圖5 Z＝0截面速度隨半徑變化分布曲線Fig.5 C urves of velocity distribution on the plane Z＝0 at different radi i

由此可見，在本文前述邊界條件下獲得穩(wěn)定旋轉(zhuǎn)流場是完全可行的，符合流體運動規(guī)律。

3.2 向心力源項修正方法

對式（4）分析，當(dāng)矢量點乘積為0時，兩矢量或者相互垂直，或者至少其一為零矢量，由此可得：

對該式兩邊做體積積分，并對左端項應(yīng)用奧高公式轉(zhuǎn)化為面積積分，可以得到壁面壓強分布與向心力源項的關(guān)系。但采用有限體積法建立離散方程時，很重要的步驟是將控制體積界面上的物理量及其導(dǎo)數(shù)通過節(jié)點物理量插值求出。雖然本文采用的離散格式在很大程度上減小了偽擴散現(xiàn)象，但針對本文的研究對象，仍至少存在以下2個問題：動量方程的求解基于直角坐標(biāo)系下的各速度分量，其方向與網(wǎng)格線呈傾斜交叉；建立離散格式時沒有考慮復(fù)雜形式源項的影響，把源項局部線性化的常規(guī)處理方法并不適用。為提高預(yù)報精度，需要對控制方程中的向心力源項做補償修正。

在得到穩(wěn)定旋轉(zhuǎn)流場的基礎(chǔ)上，本文選取與58－4621模型具有相同幾何形狀的流體介質(zhì)體作為計算對象。理論上，在穩(wěn)定的均質(zhì)旋轉(zhuǎn)流場中，這部分流體介質(zhì)體所受的向心力Fc＝ρVrω2由介質(zhì)所處的靜壓場提供。同時，對這部分流體介質(zhì)體形心處的力矩應(yīng)當(dāng)為0。依據(jù)力的相互作用原理，由介質(zhì)體表面壓強積分得到的合力和合力矩應(yīng)與航行體所受向心力及形心力矩一致。

應(yīng)用上述理論方法計算得到的向心力及形心力矩未考慮離散格式所引入的誤差、耗散。為此，首先對整個流域進行數(shù)值求解，在數(shù)值收斂的條件下，選取相同幾何形狀的包絡(luò)邊界作靜壓場積分，可以得到考慮誤差、耗散的向心力及形心力矩。通過與理論值比較，可以對數(shù)值解算引起的誤差做定量估算。運用這一原理，對R＝30.48 m條件下的多個湍流模型做計算，結(jié)果如表2所示。

表2 R＝30.48m條件下不同湍流模型向心力與形心力矩預(yù)報值Tab.2 Predictions of centripetal force and moment in different turbulencemodels at 30.48m turning radi us

可以看出，對于不同的湍流模型，計算結(jié)果與理論值的偏差各不相同，但相對偏差量很小，小于1%。該方法對邊界層網(wǎng)格非常細密的情況，很難實現(xiàn)原有邊界面兩側(cè)網(wǎng)格的均勻過渡，計算收斂困難，因此需要粗化原有的邊界層網(wǎng)格，從而導(dǎo)致修正值在一定程度上失真；但是，該方法從另一個角度驗證了本文選用的計算方法、邊界條件及網(wǎng)格形式等對獲得穩(wěn)定旋轉(zhuǎn)流場是準(zhǔn)確、可行的。

另一種求取靜壓場的方法是假定水下航行體與流體以同一角速度做回轉(zhuǎn)運動，流體與壁面間不存在相對運動，只受到流體由于回轉(zhuǎn)運動而產(chǎn)生的靜壓場的作用。對 R ＝30.48 m、22.86 m、18.288 m、15.24 m 條件下 SST k － ω 模型進行計算，結(jié)果如表3所示。由表3可見，計算值與理論值的偏差均在5%以內(nèi)。由于該方法不改變原有的網(wǎng)格分布，能較真實地反映出離散格式引入的誤差，更適合作為本文向心力源項修正的方法。

表3 不同回轉(zhuǎn)半徑條件下向心力與形心力矩預(yù)報值Tab.3 Predictions of centripetal force and moment at different turning radi i

3.3 不同半徑下的結(jié)果對比

在前述計算條件優(yōu)選的基礎(chǔ)上，數(shù)值模擬不同半徑下航行體回轉(zhuǎn)運動流場，對所受垂向力與縱傾力矩進行預(yù)報，并與試驗值比較，如圖6、圖7所示。由圖可見，數(shù)值預(yù)報的垂向力和縱傾力矩與試驗值基本吻合，隨著回轉(zhuǎn)半徑的減小、角速度值的增加，垂向力和縱傾力矩的預(yù)報值與試驗值的偏差逐漸增大，且在數(shù)值上與試驗值相比始終偏大，最大偏差不超過12%，但變化趨勢與試驗值相吻合。

圖6 預(yù)報與試驗垂向力系數(shù)隨無因次回轉(zhuǎn)角速度變化分布Fig.6 Effect of nondimensional turning radius，r＇＝（L ／R turn） on the computed and measured coefficients of normal force

圖7 預(yù)報與試驗縱傾力矩系數(shù)隨無因次回轉(zhuǎn)角速度變化分布Fig.7 Effect of nondimensional turning radius，r＇＝（L ／R turn） on the computed and measured coefficients of pitchingmoment

為了進一步分析預(yù)報結(jié)果，同時與分段模型試驗結(jié)果進行比較，將航行體不均等地劃分成10分段，劃分方式與模型試驗相同，見圖2。模擬預(yù)報回轉(zhuǎn)半徑為R＝22.86m條件下，各個分段所受垂向力值、通過比較無因次垂向力沿縱向位置分布如圖8所示。由圖可見，首尾分段預(yù)報值與試驗值偏差較大，導(dǎo)致總的垂向力與力矩存在偏差，尤其對縱傾力矩有較大影響。首部分段的偏差可能由于物理試驗時模型首部裝有激流裝置，而數(shù)值計算則通過人工設(shè)定湍流參數(shù)（經(jīng)驗公式估算，尚無定論）來模擬，兩者可能存在較大差異。進行幾何建模時，考慮結(jié)構(gòu)化網(wǎng)格劃分的需要，對尾端采用小圓弧光順，與模型的尖細尾端有較大差異，直接導(dǎo)致尾部分段值的偏差。同時，由于兩方程模型的固有缺陷，SST k－ω模型對尾部3.5 m附近區(qū)段背水面復(fù)雜分離流動的捕捉能力存在一定的不足，可能導(dǎo)致預(yù)報值偏離試驗值，湍流模型參數(shù)的調(diào)整可能對結(jié)果有所改善。

圖8 R＝22.86m條件下預(yù)報與試驗垂向力系數(shù)沿縱向位置分布曲線Fig.8 Distribution curves of the computed and measured coefficients of normal force along the longitudinal position of the vehicle at 22.86m turning radi us

此外，依據(jù)流體運動規(guī)律，隨著回轉(zhuǎn)半徑變小，航行體首端與來流夾角逐漸增大，尾部去流段流動分離點前移，且由于角速度增大，航行體周圍橫向分離流逐漸變強，分離產(chǎn)生的渦旋運動愈發(fā)強烈，并且可能隨流動下泄，出現(xiàn)嚴(yán)重的非定常分離流動，力和力矩則會出現(xiàn)不穩(wěn)定波動。與文獻［8］中對大攻角直航運動的研究相類似，此時需要考慮對流動進行非定常計算，這也是今后需要進一步研究的內(nèi)容。

4 結(jié) 論

本文采用數(shù)值求解RANS方程及湍流模型模擬水下航行體旋臂水池試驗，在選擇適當(dāng)計算條件的基礎(chǔ)上，可以實現(xiàn)Re＝11.7×106條件下回轉(zhuǎn)水動力的準(zhǔn)確預(yù)報。預(yù)報結(jié)果與試驗值相比，偏差在12%以內(nèi)，具有一定的工程實用價值。本文提出的向心力源項空間離散誤差修正方法，在一定簡化條件下，通過物理對象的轉(zhuǎn)換，較好地實現(xiàn)了水動力計算中對由旋轉(zhuǎn)流動靜壓場產(chǎn)生的向心力修正，避免對控制方程中源項的復(fù)雜數(shù)學(xué)處理，具有一定的新穎性。

SST k－ω湍流模型對回轉(zhuǎn)運動復(fù)雜流動具有良好的表征能力，但對首尾段強烈分離流動的捕捉能力仍不足。對較大無因次回轉(zhuǎn)角速度條件下，流動非定常特性的描繪，還需要具有更優(yōu)良性能的湍流模型來實現(xiàn)。

將本文數(shù)值預(yù)報方法進一步推廣應(yīng)用，可以對操縱運動方程中與定?；剞D(zhuǎn)運動相關(guān)的角速度導(dǎo)數(shù)、舵角導(dǎo)數(shù)及系列耦合水動力導(dǎo)數(shù)進行預(yù)報，最終建立水下航行體操縱性水動力的數(shù)值回轉(zhuǎn)水池。

［1］ FELDMAN J.Straight－line and rotating arm captive－model experiments to investigate the stability and control characteristics of submarine and other submerged vehicles［R］.Taylor Naval Ship Research and Development Center，1989.

［2］ SUNG C，F(xiàn)U T，GRIFFIN M，et al.Validation of incompressible flow computation of forces and moments on axisymmetric bodies undergoing constant radius turning［C］//Twenty－First Symposium on Naval Hydrodynamics，1997.

［3］ SUNG C H，JIANG M Y，RHEE B，et al.Validation of the flow around a turning submarine ［C］//The Twenty－Fourth Symposium on Naval Hydrodynamics， Fukuoka， Japan：2002.

［4］ FU T C，ATSAVAPRANEE P，HESSD E.Measurements of the c ross－f low w ake of a t urning s ubmarine m odel（ONR Body－1）［C］//25th ONR Symposium on Naval Hydrodynamics，F(xiàn)ukuoka， Japan，2002.

［5］ ATSAVAPRANEE P，F(xiàn)ORLINIT，F(xiàn)UREY D.Experimental m easurements for CFD v alidation of the f low about a s ubmarine m odel （ONR Body－1）［C］//25th ONR Symposium on Naval Hydrodynamics，St John’s，Canada，2004.

［6］ ETEBARI A， ATSAVAPRANEE P， et al.Experimental measurements on a SUBOFFmodel in a turningmaneuver［C］//27th Symposium on Naval Hydrodynamics，Seoul，Korea，2008.

［7］ 柏鐵朝，梁中剛，周軼美，等.潛艇操縱性水動力數(shù)值計算中湍流模式的比較與運用［J］.中國艦船研究，2010，5（2）： 22－28.

BAIT C，LIANG Z G，ZHOU Y M，et al.Comparison and a pplication of t urbulence m odes in s ubmarine m aneuvering h ydrodynamic f orces c omputation ［J］.Chinese Journal of Ship Research，2010，5（2）： 22－28.

［8］ BRIDGES D H.A detailed study of the flow field of a submarine at large angle drifts ［R］.US：Office of Naval Research， 2001.

Numerical Calculation on Hydrodynam ic Performance of the Submerged Vehicle in Turning Motion

Lu Jin－guo Liang Zhong－gang Wu Fang－liang Zhou Yi－mei
China Ship Developmentand Design Center，Wuhan 430064，China

A num erical procedure based on solving the incompressible RANSequations in a steady rotating coordinate system was developed for the prediction of forces and moments on submerged vehicle in turningmotion with the commercial code FLUENT.New method considered numerical dissipation was developed for the correction of centripetal force，and several key elements influenced the solution were also taken into account.In the case of Re＝11.7×106， themaximum relative deviation between the calculation results and measured data is less than 12%.The good agreement proves that the method is practical for numerical prediction of the hydrodynamic performance on the submerged vehicle in turningmotion.

turnin gmotion； hydrodynamics； centripetal force correction； FLUENT； submerged vehicle

U661.3

A

1673－3185（2011）06－08－05

10.3969/j.issn.1673-3185.2011.06.002

2011-06-28

盧錦國（1986－），男，碩士研究生。研究方向：潛艇操縱性水動力數(shù)值計算。E-mail：ljgfatherland＠126.com

梁中剛（1966－），男，研究員，碩士生導(dǎo)師。研究方向：船舶總體設(shè)計。E-mail：lzg701papers＠126.com