尤天慶，王占瑩，程少華，穆 洲，鮑文春

（北京宇航系統(tǒng)工程研究所，北京，100076）

水下航行體流體阻尼力系數(shù)的cFd計算研究

尤天慶，王占瑩，程少華，穆 洲，鮑文春

（北京宇航系統(tǒng)工程研究所，北京，100076）

水下航行體流體阻尼力關(guān)系到彈道的穩(wěn)定性和操縱性，當(dāng)航行體以較大角速度機動時阻尼力更是不可忽略，因此準(zhǔn)確預(yù)示航行體流體阻尼力具有重要意義。以計算流體動力學(xué)（computational Fluid dynamic，cFd）方法為基礎(chǔ)，結(jié)合動網(wǎng)格和動坐標(biāo)系技術(shù)，進行基于雷諾平均N-S方程的阻尼力計算方法研究，發(fā)展針對水下航行體附加力系數(shù)和俯仰阻尼力矩系數(shù)的數(shù)值計算方法，通過與試驗數(shù)據(jù)對比表明，該方法具有較好精度，在此基礎(chǔ)上分析了空化現(xiàn)象對流體阻尼力系數(shù)的影響。

計算流體動力學(xué)；水下航行體；流體阻尼力；空化

0 引 言

潛射導(dǎo)彈、魚雷、水雷和潛艇等水下航行體在水下運動過程中所受的流體力為彈道或控制系統(tǒng)的主要輸入?yún)?shù)，關(guān)系到航行體彈道穩(wěn)定性、操控性或航程等總體參數(shù)，是各類型航行體設(shè)計的關(guān)鍵參數(shù)。航行體所受流體動力主要取決于流場性質(zhì)、航行體流體動力外形以及相對流體的運動，迄今為止，還難以用理論計算或模型試驗的方法直接得到包含上述所有因素的總流體動力，只能在各種假設(shè)和簡化下，把流體動力分解成若干部分，分別由計算或試驗的手段獲得，然后綜合成總流體動力[1]。通常在航行體設(shè)計過程中，將流體動力分為位置力、阻尼力和慣性力3個部分。航行體以一定速度作平移運動的同時，還作旋轉(zhuǎn)運動，由旋轉(zhuǎn)運動而產(chǎn)生的流體力增量為阻尼力[2]。流體動力中阻尼力（如俯仰阻尼力矩、附加力等）關(guān)系到航行體彈道穩(wěn)定性及操縱性[3，4]，當(dāng)旋轉(zhuǎn)運動角速度較大時阻尼力更是不可忽略，因此阻尼力計算在水動力設(shè)計中有非常重要的作用。

目前，無論是氣動力還是水動力，針對流體阻尼力的計算主要有工程估算、試驗以及計算流體力學(xué)（computational Fluid dynamic，cFd）等方法。其中，工程估算方法計算速度快但精度較差，只能作為初步估算[5]；試驗所得的數(shù)據(jù)較精確，但周期長、費用高；cFd計算方法比勢流理論更具從理論到工程實踐的可擴展性，比試驗方法更容易實現(xiàn)，成本低。

近年來，隨著cFd及計算機技術(shù)發(fā)展，基于N-S方程求解流體阻尼力取得一定的發(fā)展。在氣動力計算方面，分別采用穩(wěn)態(tài)和非穩(wěn)態(tài)求解方法計算飛行器在空氣中飛行的流體阻尼力[6]。其中，穩(wěn)態(tài)求解的方法即通過改造N-S方程，在方程中加入慣性力源項以求解勻速圓周運動飛行器受力情況，得到飛行器的定常拉升運動中阻尼力[7，8]；非穩(wěn)態(tài)求解即計算航行體強迫振蕩運動中的俯仰力矩遲滯曲線，通過辨識得到計算模型的俯仰穩(wěn)定性參數(shù)[9]。不同于氣動力，當(dāng)水下航行體高速運動或通氣時，航行體表面會附著有空泡，空泡發(fā)展過程帶有很強的非定常和非線性特性，這同樣給流體阻尼力的確定造成很大困難。由于水的密度較大，航行體水中運動的流體慣性力（附連水質(zhì)量）不可忽略，計算或試驗的阻尼力系數(shù)會有附連水質(zhì)量的影響。

本文基于通用cFd計算程序，運用動坐標(biāo)系穩(wěn)態(tài)計算求解流體阻尼力，并通過與試驗結(jié)果對比驗證了計算方法的有效性。在此基礎(chǔ)上對數(shù)值計算結(jié)果進行討論，分析了空化現(xiàn)象對流體阻尼力系數(shù)的影響。

1 數(shù)值計算方法

1.1 阻尼力計算方法

航行體所受的流體阻尼力系數(shù)與坐標(biāo)系選取緊密相關(guān)，因此首先介紹本體坐標(biāo)系的定義。本文航行體本體坐標(biāo)系與航行體固連，原點選在航行體質(zhì)心，Ox1軸與航行體縱軸重和指向前方，Oy1軸位于航行體縱對稱面內(nèi)與Ox1軸垂直指向上方，Oz1軸與Ox1y1平面垂直并與Ox1軸和Oy1軸構(gòu)成右手系，如圖1所示。

圖1 本體坐標(biāo)系

無論采用試驗還是數(shù)值計算手段獲得航行體流體阻尼力，其基本思路是構(gòu)造一種運動方式，使試驗測得或數(shù)值計算得到航行體所受的流體力包含所關(guān)注的流體力，并能通過數(shù)據(jù)處理的手段得到所關(guān)注的流體力。這里構(gòu)造的運動方式為勻速圓周運動，數(shù)值模擬實現(xiàn)方案為采用旋轉(zhuǎn)坐標(biāo)系的方法。

令航行體以攻角α繞指定圓心做勻速圓周運動，航行體質(zhì)心與圓周運動圓心距離為R，轉(zhuǎn)動角速度為ωz1，則勻速圓周運動線速度為V=Rωz1。在本體坐標(biāo)系下，航行體所受到的流體力有如下表達式：

式中 q為動壓；S，l分別為參考面積和參考長度；N為法向力；λ為附連水質(zhì)量；M為俯仰力矩。

通過航行體定速圓周運動數(shù)值計算，可求得指定攻角α和無量綱旋轉(zhuǎn)角速度1下的法向力N和俯仰力矩M。式（1）和式（2）等號右側(cè)變量均為已知，由此可求得1和。

由于λ11V ωz1cosα項和λ26V ωz1cosα 項內(nèi)均含有ωz1，因此在彈道計算中通常將含有λ11的λ11Vωz1cosα項記入內(nèi)，將含有λ26的λ26Vωz1cosα 項記入1內(nèi)。1和的計算公式為

數(shù)值計算控制方程為基于雷諾平均N-S方程的質(zhì)量守恒方程、動量守恒方程和kε-湍流模型方程。該計算流場建立在旋轉(zhuǎn)坐標(biāo)系上，在此坐標(biāo)系下，網(wǎng)格在計算中保持靜止。其中：

質(zhì)量守恒方程為

動量守恒方程為式中rU為相對旋轉(zhuǎn)坐標(biāo)系的速度矢量；ω為旋轉(zhuǎn)坐標(biāo)系旋轉(zhuǎn)角速度矢量；r為流場質(zhì)點在旋轉(zhuǎn)坐標(biāo)系中的位置矢量；p為壓力；t為時間；effμ為有效粘度。

本文采用標(biāo)準(zhǔn)kε-湍流模型來建立流場中時均統(tǒng)計量和脈動量之間的關(guān)系。kε-模型是目前比較通用的一種湍流模型，在計算上易于處理，魯棒性也較強。

本文數(shù)值離散方法采用有限體積法，湍流輸運方程采用一階差分格式，時間離散為一階隱式。壓力速度耦合方式采用SIMpLe算法，壓力及動量方程離散采用二階迎風(fēng)格式。

1.2 空化模型

航行體在水下高速航行時，周圍流體與航行體相對運動，當(dāng)航行體表面壓力下降至飽和蒸汽壓力以下并維持一定時間后，空泡將產(chǎn)生。流場中是否存在空泡，依賴于流動條件和航行體的幾何形狀。在空化現(xiàn)象中，通常用無量綱參數(shù)σ來表征自然空化狀態(tài)。

應(yīng)力張量τ與應(yīng)變率關(guān)系為

式中 p∞為流場中未受擾動處的壓力；vp為水的飽和蒸汽壓；V為航行體運動速度。

數(shù)值計算采用基于輸運方程的空化模型?？栈紤]為水與水蒸汽之間的質(zhì)量轉(zhuǎn)換，多相流中的輸運方程為

式中 Γib為b相進入i相的單位體積質(zhì)量源（本文中b相為水，i相為水蒸汽），Γib=m˙ibAib，其中，b為b到i相單位相間交界面質(zhì)量流率；Aib為相間界面面積。

根據(jù)空泡動力學(xué)可知，氣泡壁面的運動方程為

式中BR為氣泡半徑；τ為表面張力；p為流體參考壓力；lρ為液體密度。忽略式（10）中表面張力項以及高次項，化簡后得：

氣體體積含量為γb=VBNB=4πNB/3，其中NB，VB分別為單位體積氣泡數(shù)和氣泡體積。

水和水蒸汽間單位體積總質(zhì)量傳遞率為

式中 F為實驗獲得的經(jīng)驗系數(shù)；gρ為水蒸汽密度。

2 數(shù)值計算結(jié)果討論

2.1 數(shù)值計算結(jié)果的有效性

為驗證阻尼力計算方法的準(zhǔn)確程度，計算一航行體在不同無量綱旋轉(zhuǎn)角速度ωz1條件下的法向力N及俯仰力矩M，圖2為航行體周圍流場速度云圖。通過對計算數(shù)據(jù)處理得到附加力系數(shù)Cωz1和俯仰阻尼力矩

。數(shù)值計算結(jié)果與實驗數(shù)據(jù)無量綱化對比結(jié)果如表1所示。

圖2 α=2?1=0.46航行體周圍流場速度云圖

表1 和m數(shù)值計算結(jié)果與試驗對比

表1 和m數(shù)值計算結(jié)果與試驗對比

系數(shù) 攻角/(°) 數(shù)值計算 試驗結(jié)果 誤差0 0.660 0.628 5% Cωz1y12 0.655 0.618 6% 6 0.527 0.500 5% 0 -0.376 -0.385 2% mωz1z12 -0.394 -0.385 2% 6 -0.380 -0.394 4%

表1中的數(shù)值計算結(jié)果表明，所采用的流體阻尼力數(shù)值計算方法具有較高的精度。

以往研究結(jié)果表明圓柱航行體長徑比，對航行體阻尼力有較大的影響，并且其影響近似為線性[1]。為驗證數(shù)值計算結(jié)果規(guī)律的正確性，在保持質(zhì)心相對位置cg/xL和航行體頭部外形不變的情況下，計算了不同長徑比/L d航行體的附加力系數(shù)和俯仰阻尼力矩系數(shù)（見圖3、圖4）。數(shù)值計算結(jié)果表明隨著長徑比的增大，附加力系數(shù)和俯仰阻尼力矩系數(shù)均逐漸增大。其中附加力系數(shù)增加幅度較小，俯仰阻尼力矩系數(shù)增幅較大，長徑比/L d在7.1～14.2范圍內(nèi)逐漸增加時，附加力系數(shù)在0.905～0.951內(nèi)緩慢遞增，增幅約為5%，俯仰阻尼力矩系數(shù)絕對值在0.453～0.573內(nèi)近似程線性逐漸遞增的趨勢，增幅約為26%。

圖3 不同長徑比的附加力系數(shù)

圖4 不同長徑比的俯仰阻尼力矩系數(shù)

2.2 空化現(xiàn)象的影響

航行體表面空泡發(fā)生在初生空化臨界狀態(tài)之后，水流從物面脫流，在壁面附近形成相對穩(wěn)定的氣水交接面。但空泡發(fā)展具有一定的非定常特性，即在來流穩(wěn)定條件下，空泡形態(tài)也會存在一定的波動，同時航行體的受力也存在一定的波動。在此忽略這種空泡脈動的影響，采用均質(zhì)多相流模型結(jié)合空化模型，分析了不同空化數(shù)條件下的附加力系數(shù)和俯仰阻尼力矩系數(shù)。

在數(shù)值計算中采用減小環(huán)境壓力的方法，使航行體自然空化數(shù)逐漸降低，自然空泡覆蓋的區(qū)域逐漸增大。在空泡覆蓋航行體表面的情況下，空泡覆蓋區(qū)域內(nèi)航行體表面壓力為恒定，空泡末端與航行體壁面接觸時，水流損失徑向速度形成壓力峰值，如圖5所示。

圖5 0.24σ=頭部空泡形態(tài)及周圍流場壓力云圖

圖6 不同空化條件附加力系數(shù)分布

圖7 不同空化條件的俯仰阻尼力矩系數(shù)分布

針對空泡覆蓋區(qū)域，由于空泡內(nèi)部壓力恒定，空泡覆蓋區(qū)域內(nèi)不存在分布法向力，且隨著旋轉(zhuǎn)角速度的增加法向力不會發(fā)生變化，因此在空泡覆蓋區(qū)域，分布附加力系數(shù)和俯仰阻尼力矩系數(shù)均為零，在圖6和圖7中表現(xiàn)為分布阻尼力系數(shù)曲線的平直段。

針對回射壓力影響區(qū)域，航行體圓周運動過程中，頭部空泡覆蓋區(qū)域存在局部攻角，使空泡末端閉合處迎、背水面空泡回射壓力存在明顯不同，迎水面壓力較背水面壓力峰值在位置上略靠前，但空間上分布范圍較窄。因此在空泡末端，沿航行體軸線方向分布法向力在空間上先后出現(xiàn)負(fù)正兩個峰值。隨著旋轉(zhuǎn)角速度變化，空泡區(qū)的局部攻角不斷發(fā)生變化，空泡迎背水面回射壓力的位置和強度隨之而發(fā)生變化，即空泡末端的負(fù)正法向力峰隨著旋轉(zhuǎn)角速度變化而大幅度變化，因此在此區(qū)域的分布附加力系數(shù)和俯仰阻尼力矩系數(shù)也先后存在負(fù)正兩個峰值。這種阻尼力系數(shù)分布情況，在圖6和圖7中體現(xiàn)為曲線先下降而后緩慢上升，最終與全濕流曲線變化規(guī)律一致。

隨著空泡數(shù)降低，空泡覆蓋的長度逐漸增大，表現(xiàn)為圖6和圖7中，阻尼力系數(shù)曲線平直段逐漸增大。不同空泡數(shù)條件下，空泡閉合區(qū)位置發(fā)生移動，相應(yīng)的回射區(qū)局部攻角發(fā)生變化，從而導(dǎo)致分布流體阻尼力系數(shù)正負(fù)峰值存在差異。以上兩方面因素共同作用，使得航行體整體的附加力系數(shù)和俯仰阻尼力矩系數(shù)隨著空泡數(shù)降低而增大，如圖8和圖9所示。

圖8 不同空化條件的附加力系數(shù)

圖9 不同空化條件的俯仰阻尼力矩系數(shù)

3 結(jié) 論

本文基于cFd計算，利用動坐標(biāo)系技術(shù)，發(fā)展水下航行體流體阻尼力的數(shù)值計算方法，在此基礎(chǔ)上，研究了空化現(xiàn)象的影響，得到如下結(jié)論：

a）基于旋轉(zhuǎn)坐標(biāo)下雷諾平均N-S方程，求解水下航行體流體阻尼力系數(shù)，通過與試驗數(shù)據(jù)以及以往規(guī)律性認(rèn)識對比驗證計算方法的準(zhǔn)確性；本文數(shù)值計算方法能考慮流動非線性、渦以及粘性等因素，同時適用于任何幾何外形的航行體，其對于工程問題更具有適用性和可擴展性。

b）針對空化現(xiàn)象的影響，在空泡內(nèi)部，壓力恒定并且隨著旋轉(zhuǎn)角速度增加分布法向力保持不變，分布流體阻尼力系數(shù)為零；在回射壓力影響區(qū)域，分布阻尼力系數(shù)沿航行體軸線先下降而后緩慢上升；隨著空泡數(shù)減小，空泡覆蓋長度增加，航行體整體的附加力系數(shù)增大，俯仰阻尼力矩系數(shù)絕對值減小。

[1] 榮建徳. 水下運載器性能的分析與設(shè)計[M]. 北京: 國防工業(yè)出版社, 2008.

[2] 張宇文. 魚雷總體設(shè)計原理與方法[M]. 西安: 西北工業(yè)大學(xué)出版社, 1998.

[3] 蔣勝矩, 劉玉琴, 黨明利. 基于定常N-S方程的飛行器滾轉(zhuǎn)阻尼力矩系數(shù)導(dǎo)數(shù)計算方法[J]. 彈箭與制導(dǎo)學(xué)報, 2008, 28(1): 180-182.

[4] Sen d. a study on sensitivity of maneuverability performance on the hydrodynamic coefficients for submerged bodies[J]. Journal of Ship Research, 2000, 44(3), 186-196.

[5] danberg J e, Weinacht p. approximate computation of pitch-damping coefficients[c]. Monterey: aIaa atmospheric Flight Mechanics conference and exhibit, 2002.

[6] park S H, Kwon J H. comparisons of steady and unsteady methods for pitch-damping predictions[c]. Orlando: 21stapplied aerodynamics conference, 2003.

[7] Weinacht p, Sturek W B, Schiff L B. Navier-stokes predictions of pitch-damping for axisymmetric shell using steady coning motion[R]. ad-a285866, 1994

[8] Weinacht p. Navier-stokes predictions of the individual components of the pitch-damping coefficient sum[R]. aRL-TR-3169, 2004.

[9] 陶洋, 等. 基于cFd的帶控制舵導(dǎo)彈的動導(dǎo)數(shù)數(shù)計算[J]. 航空動力學(xué)報, 2010, 25(1): 102-105.

[10] Knapp R T, daily J W, Hammitt F G. cavitation[M]. New York: Mcgraw-Hill Book company, 1970.

[11] Kunz R F, Boger d a, chyczewski T S, et al. Multi-phase cFd analysis of natural and ventilated cavitation about submerged bodies[c]. San Francisco: Third aSMe/JSMe Joint Fluids engineering conference, 1999.

[12] Knapp R T, daily J W, Hammitt F G. cavitation[M]. New York: Mcgraw-Hill Book company, 1970.

CFd Research on Underwater Vehicle Hydrodynamic damping Force Coefficient

You Tian-qing, Wang Zhan-ying, cheng Shao-hua, Mu Zhou, Bao Wen-chun
(Beijing Institute of aerospace System engineering, Beijing, 100076)

Underwater vehicle hydrodynamic damping force is crucial to trajectory stability and maneuverability, especially for the vehicle with high angle velocity, thus the accurate prediction of hydrodynamic damping force is important. Based on the cFd method, Reynolds averaged N-S equations based hydrodynamic damping force calculation research is conducted. combining with dynamic mesh and rotating reference frames, pitch damping Force coefficient calculation methods are introduced. compared with experiment result, the validity and accuracy is proofed. Furthermore, the effect of cavitation is discussed.

computational Fluid dynamics; Underwater vehicle; Hydrodynamic damping force; cavitation

O351.2

a

1004-7182(2016)01-0066-05

10.7654/j.issn.1004-7182.20160115

2015-04-23；

2015-06-02

水動力學(xué)教育部重點試驗室資助項目

尤天慶（1984-），男，工程師，主要研究方向為計算流體力學(xué)