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

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

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

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

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

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

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

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

0 引 言

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

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

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

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

1 數(shù)值計(jì)算方法

1.1 阻尼力計(jì)算方法

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

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

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

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

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

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

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

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

質(zhì)量守恒方程為

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

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

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

1.2 空化模型

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

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

式中 p∞為流場(chǎng)中未受擾動(dòng)處的壓力；vp為水的飽和蒸汽壓；V為航行體運(yùn)動(dòng)速度。

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

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

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

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

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

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

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

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

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

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

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

圖2 α=2?1=0.46航行體周圍流場(chǎng)速度云圖

表1 和m數(shù)值計(jì)算結(jié)果與試驗(yàn)對(duì)比

表1 和m數(shù)值計(jì)算結(jié)果與試驗(yàn)對(duì)比

系數(shù) 攻角/(°) 數(shù)值計(jì)算 試驗(yàn)結(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ù)值計(jì)算結(jié)果表明，所采用的流體阻尼力數(shù)值計(jì)算方法具有較高的精度。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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

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

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

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

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

[3] 蔣勝矩, 劉玉琴, 黨明利. 基于定常N-S方程的飛行器滾轉(zhuǎn)阻尼力矩系數(shù)導(dǎo)數(shù)計(jì)算方法[J]. 彈箭與制導(dǎo)學(xué)報(bào), 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òng)導(dǎo)數(shù)數(shù)計(jì)算[J]. 航空動(dòng)力學(xué)報(bào), 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

水動(dòng)力學(xué)教育部重點(diǎn)試驗(yàn)室資助項(xiàng)目

尤天慶（1984-），男，工程師，主要研究方向?yàn)橛?jì)算流體力學(xué)