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

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

（北京宇航系統工程研究所，北京，100076）

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

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

（北京宇航系統工程研究所，北京，100076）

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

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

0 引 言

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

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

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

本文基于通用cFd計算程序，運用動坐標系穩態計算求解流體阻尼力，并通過與試驗結果對比驗證了計算方法的有效性。在此基礎上對數值計算結果進行討論，分析了空化現象對流體阻尼力系數的影響。

1 數值計算方法

1.1 阻尼力計算方法

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

圖1 本體坐標系

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

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

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

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

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

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

質量守恒方程為

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

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

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

1.2 空化模型

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

應力張量τ與應變率關系為

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

數值計算采用基于輸運方程的空化模型。空化考慮為水與水蒸汽之間的質量轉換，多相流中的輸運方程為

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

根據空泡動力學可知，氣泡壁面的運動方程為

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

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

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

式中 F為實驗獲得的經驗系數；gρ為水蒸汽密度。

2 數值計算結果討論

2.1 數值計算結果的有效性

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

。數值計算結果與實驗數據無量綱化對比結果如表1所示。

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

表1 和m數值計算結果與試驗對比

表1 和m數值計算結果與試驗對比

系數 攻角/(°) 數值計算 試驗結果 誤差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中的數值計算結果表明，所采用的流體阻尼力數值計算方法具有較高的精度。

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

圖3 不同長徑比的附加力系數

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

2.2 空化現象的影響

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

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

圖5 0.24σ=頭部空泡形態及周圍流場壓力云圖

圖6 不同空化條件附加力系數分布

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

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

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

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

圖8 不同空化條件的附加力系數

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

3 結 論

本文基于cFd計算，利用動坐標系技術，發展水下航行體流體阻尼力的數值計算方法，在此基礎上，研究了空化現象的影響，得到如下結論：

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

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

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

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

[3] 蔣勝矩, 劉玉琴, 黨明利. 基于定常N-S方程的飛行器滾轉阻尼力矩系數導數計算方法[J]. 彈箭與制導學報, 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的帶控制舵導彈的動導數數計算[J]. 航空動力學報, 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

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

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