沈德璋，張 合，常自剛，李豪杰，齊獻山

（1.中國工程物理研究院 電子工程研究所，綿陽 621900；2.南京理工大學 智能彈藥技術國防重點學科實驗室，南京 210094；3.海裝駐重慶地區(qū)軍事代表局，重慶 400042）

0 引 言

水中彈藥的發(fā)射環(huán)境同陸上相比有很大差異，引信環(huán)境敏感及安全保險控制成為水中彈藥引信設計的關鍵技術。采用頭部外置共形渦輪是水下火箭彈引信實現遠距離解除保險的一種途徑［1］。彈丸水下高速飛行時的渦輪轉動情況及流場數據是進行系統參數設計的重要依據。由于模型復雜，且涉及空化等復雜現象，理論分析很難解決此類問題。而試驗研究所需周期長，成本高。就水洞模擬試驗而言，國內目前僅能達到20m／s以內的流速，水下靶道射擊試驗則受測量手段局限，能獲得的數據十分有限。為預測高速發(fā)生空化情況下的渦輪轉速，獲得渦輪葉片附近壓力分布、空化情況、阻力系數等數據，數值仿真是一項經濟且實用的方法。隨著計算機技術的發(fā)展，CFD技術現已成為研究工程流體問題的重要技術手段。

徐長江［2］等提出一種將三維流場數值模擬與電機動力學分析相結合的迭代模擬計算方法，解決了由于引信渦輪電機轉速不確定而導致三維流場數值模擬時初始條件無法確定的問題，其針對的是陸上彈藥引信用側進氣渦輪，模擬流體為空氣的單相流動。關于渦輪機械的仿真研究較多，其中絕大多數為內流問題。王振［3］等應用 S－A、標準k－ε、RNG k－ε、Realizable k－ε以及標準k－ω這五種湍流模型，對渦輪流量傳感器內部流場進行了三維數值模擬，其主要針對工業(yè)用渦輪流量傳感器的工況，不涉及空化情況。黃劍峰［4］等基于歐拉－歐拉方法中均勻多相流假設的混合兩相流體無滑移模型，對某混流式水輪機全流道進行了三維定?？栈牧鲾抵的M。結果表明可以較好地模擬水輪機內真實的有空化發(fā)生的多相流動情況。LIU S H、ZHANG L、WU YL［5－6］等研究了采用RNGk－ε湍流模型及Mixure多相流模型模擬混流式水輪機中氣液兩相湍流流動的方法，認為該方法具有足夠的潛力模擬水輪機空化流。SANSONE E［7］等使用非定常的Navier－Stokes方程模擬了橫流式水輪機空化及非空化流動，分析了氣蝕的動態(tài)失速現象及對渦輪性能的影響，認為空化初期出現了葉片失速現象并放大了渦輪脫落渦結構的尺寸。

本文針對彈藥在水中高速運行過程，開展渦輪高速旋轉且發(fā)生空化時的流場兩相流動數值仿真研究，模擬不同速度及水深情況下渦輪附近流場的變化，為渦輪優(yōu)化設計及安全系統參數設計提供支撐。

1 控制方程及計算模型

1.1 控制方程

不可壓牛頓流體三維Navier－Stokes控制方程的一般形式：

上述方程對層流和湍流均是適用的。但對于湍流如果直接求解，需要采用對計算機內存和計算速度要求很高的直接模擬方法，目前工程上采用的主要方法是對控制方程做時間平均處理，同時補充反應湍流特性的其他方程，如常用的雙方程湍流模型k－ε模型和k－ω 模型。

1.2 旋轉坐標系

在旋轉坐標系下求解質量守恒以及連續(xù)性方程時，動量方程中需要添加流體加速度項。在Fluent中求解旋轉坐標系問題時可使用絕對速度或相對速度。兩者關系為：

式中：v為絕對速度；vr為相對速度；Ω為旋轉角速度；r為旋轉坐標系下位置向量。

旋轉域中質量守恒方程或者連續(xù)性方程可參考下式寫成絕對速度或相對速度形式：

1.3 空化模型

在空化流模擬時，采用Schnerr－Sauer［8］空化模型。Schnerr－Sauer空化模型氣相運輸控制方程：

不同于完整空化模型，Schnerr和Sauer推導出了精確的從液體到空泡的凈相變率方程：

以上兩式源于方程：

其中，R為質量遷移率，?B為氣泡半徑，v為氣相，a為氣相體積比，ρv為氣相密度，Re和Rc分別為與氣泡生長和潰滅相關的質量遷移源項。

2 網格劃分及邊界條件

2.1 幾何模型及網格劃分

整個渦輪機構零部件較多，主要包括：小端蓋、渦輪、磁體、霍爾器件、連接件、引信體、密封圈、傳感器座、緊固圈、軸承等。本文主要分析外部流場的情況，如完全按照某型火箭彈引信原型樣機建模，無疑十分復雜，且將對網格劃分工作帶來極大的困難。因此本文將在計算時對模型做如下簡化：

（1）忽略渦輪內部軸承、小磁體、緊固圈等結構，只考慮渦輪外部流場；

（2）原型樣機中渦輪與引信體之間的間隙忽略不計；

（3）不考慮渦輪加工工藝中產生的倒角；

（4）忽略渦輪端面用于安裝而加工的淺槽。

以渦輪入口葉頂直徑φ為參考量，渦輪部分的主要幾何參數如表1所示。

表1 渦輪主要幾何參數Table 1 The main geometric parameters of turbine

表中Φhi為入口輪轂直徑；Φhe為出口輪轂直徑；Φte為出口葉頂直徑；Tb為葉片厚度；Lh為輪轂長度；Ld為葉片導程；N為葉片數。

為模擬渦輪旋轉，需采用三維數值仿真方法。由于本文涉及的是外部繞流問題，為保證計算結果精度，將渦輪直徑φ作為特征尺寸，取整個計算流場長度L＝30φ，直徑D＝20φ。

圖1 全流場計算區(qū)域Fig.1 Flow field computational domain

整個流場主要可分為如圖1幾個區(qū)域，a、d為渦輪前后區(qū)，b、c為渦輪區(qū)，其中b為旋轉區(qū)域。整個計算流場體積較大，為控制網格數量，幾何模型采用周期性邊界，根據渦輪葉片數，取全流場的1／8為計算區(qū)域。

計算區(qū)域以45°扇形截面沿一定路徑進行切割產生，切割路徑為葉片前后端連線在扇形角平分線垂直平面內的投影方向。

結構剖分方法如圖2所示，圖中扇形面SOAB沿上述直線路徑平移至面SO′A′B′，該直線投影 OO′與∠AOB平分線OX垂直。

圖2 結構剖分方法示意圖Fig.2 Schematic diagram of the last division method

平移后邊線投影之間的距離A′C＝B′C′，而OA′＝OB′，即有：sin∠AOA′＝sin∠BOB′，∠AOA＝∠BOB′。從而：∠AOB＝∠A′OB′，易得平移后所得圓弧，滿足周期性條件。

為劃分結構性網格，須保證對邊節(jié)點數一致的原則。本文涉及的扇形體，對三角形部分Y型剖分，外部圓環(huán)部分進行O型剖分，將渦輪部分劃分為一個單獨的連續(xù)體，以便后續(xù)使用旋轉參考系計算流場旋轉的情況。本文對幾何模型全部按結構性六面體網格劃分，以保證計算精度。為控制計算量，在保證網格質量的同時須控制網格總數。經大量調整，完成網格劃分，網格總數約90萬。渦輪區(qū)域部分網格劃分最為復雜，劃分完成的渦輪附近部分細節(jié)如圖3所示。

圖3 渦輪附近部分細節(jié)Fig.3 Details of turbine nearby area

2.2 邊界條件及求解算法

采用壓力基求解器進行穩(wěn)態(tài)下的定常計算，采用速度入口、壓力出口邊界，對動區(qū)采用旋轉坐標系模擬渦輪轉動情況。入口速度設置為彈速，出口壓力定義為彈丸所處水深位置的壓力，假設整個出口壓力為恒定值。由于計算前渦輪轉速未知，旋轉區(qū)域轉速初始條件確定相對復雜一點。

當渦輪在穩(wěn)定工作狀態(tài)下時，旋轉的角加速度等于0，渦輪力矩平衡方程為：式中Md為水流作用在渦輪葉片上引起的旋轉驅動力矩；Mf10和Mf11分別為無負荷時的軸承摩擦力矩和由負荷引起的摩擦力矩。Mh為輪轂端面的粘性摩擦阻力矩；Ms輪轂外表面粘性摩擦阻力矩；Mi為輪轂內表面粘性摩擦阻力矩；Mt為葉片頂端粘性摩擦阻力矩。

其中軸承摩擦力矩采用Palmgren提出軸承摩擦力矩計算公式進行計算［9］，輪轂內表面粘性摩擦阻力矩Mi根據趙學端［10］建立的輪轂內表面粘性摩擦阻力矩公式計算，渦輪驅動力矩、渦輪及葉片表面摩擦力矩均由數值仿真獲得。

渦輪轉速邊界由以下方法確定：（1）根據經驗設置初值R0，進行三維數值仿真；（2）根據仿真結果計算驅動力矩及各阻力矩之和，并確定下一次計算的轉速值R1；（3）重復步驟（1）～（2），直至上式達到平衡，渦輪葉片力矩之和小于ε，此處取ε＝1×10－6。

湍流模型采用Realizable k－ε模型，近壁處理采用標準壁面函數。多相模擬采用Mmixture多相流模型。壓力－速度耦合方程求解算法采用半隱式連接壓力方程方法，即SIMPLE算法，該方法是默認算法，穩(wěn)健性好。單元中心的變量梯度基于單元體的最小二乘法插值（Least Squares Cell Based方法），主要用于多面體網格，與基于節(jié)點的格林－高斯格式具有相同的精度和格式。壓力插值算法采用PRESTO！格式，該格式主要用于高旋流，壓力急劇變化流（如多孔介質、風扇模型等），或劇烈彎曲的區(qū)域。對流項的插值采用QUICK格式，此格式適用于四邊形／六面體以及混合網格，對旋轉流動有用，在均勻網格上能達到三階精度。

3 計算結果及分析

本文模擬了水深10m、30m及50m條件下，彈速5～100m／s情況下的渦輪高速旋轉時的兩相流動流場數據。

3.1 空化情況分析

由于外形根據原水下平頭彈丸的設計，渦輪表面及葉片部分為彈丸頭部最易發(fā)生空化的部位。在模擬水深位置為10m情況下，彈速約30m／s時在引信渦輪表面初生空化。不同速度下渦輪部分的空化程度如圖4所示。

從圖中可以看出，渦輪表面處的空化尺度隨彈速增大而增大。彈速30m／s時渦輪表面氣相最大體積分數約為34.4%，彈速大于50m／s時渦輪表面部分區(qū)域完全空化，氣相體積分數最大處接近100%。彈速80m／s時渦輪表面大部分區(qū)域發(fā)生空化，葉片背面的空化程度比葉片正面更為嚴重。彈速100m／s的情況下則整個渦輪表面均存在不同程度的空化現象。空化的發(fā)生將對渦輪的驅動力矩產生一定影響，進而影響渦輪轉動速度。渦輪表面氣相體積分數Vp在各區(qū)間內的單元數En及所占比例Ep如表2所示。

3.2 壓力分布及阻力系數

不同速度下渦輪輪轂表面軸向的表面壓力分布曲線如圖5所示。

圖4 不同彈速下的空化情況Fig.4 The cavitation under different speeds

表2 渦輪表面氣相體積分數Vp在各區(qū)間內的單元數En及所占比例EpTable 2 Unit number En and proportion Ep of turbine surface gas phase volume fraction Vp in each interval

圖5 渦輪輪轂表面壓力軸向分布Fig.5 Axial pressure distribution of turbine hub surface

渦輪輪轂表面壓力在彈速10m／s時變化很小，且均大于0，故不會發(fā)生空化。彈速在30m／s以上時，最低壓力均為－9.78×104Pa，此即為液態(tài)水常溫下的飽和蒸汽壓?？梢钥闯觯S著速度的增大，輪轂表面為最低壓力值的軸向距離也在增大，這一趨勢正好反映了空化區(qū)域尺度的變化，與圖4所表達的情況一致。最大壓力出現在渦輪后端，隨速度增大而增大。

整個彈丸頭部軸向阻力系數如表3所示。表中Cpx、Cvx、Ctx分別為頭部壓力系數、頭部粘性阻力系數及頭部總阻力系數。

根據表3數據繪制彈丸頭部阻力系數曲線如圖6所示。

表3 頭部軸向阻力系數數據Table 3 Head axial drag coefficients

圖6 彈丸頭部阻力系數曲線Fig.6 Projectile head drag coefficient curve

該彈丸頭部裝有渦輪，在水下以10m／s速度運動時頭部總阻力系數為0.5353，100m／s時的總阻力系數為0.0358，可以看出總阻力系數隨速度的增大而減小。10m／s時粘性阻力在總阻力中所占比例為1.8%，隨速度增大，粘性阻力的影響逐漸明顯，100m／s速度下粘性阻力占總阻力的18.7%。

3.3 渦輪轉速計算

通過迭代方法獲得了不同水深（10m、30m及50m）及彈速條件下的渦輪轉速，部分數據如表4所示。

表4 不同水深及速度下的渦輪轉速（r／min）Table 4 Turbine rotation speed at different water depths and speeds（unit：r／min）

在10m、30m及50m水深條件下，繪制出渦輪轉速與彈速的關系曲線如圖7所示。

圖7 不同水深條件下彈速與渦輪轉速關系曲線Fig.7 The relationship curve of projectile speed and turbine rotation speed at different water depths

根據上述曲線可以看出，靜壓力對渦輪轉速存在一定影響。在彈速為100m／s時，50m水深條件下的渦輪轉速比10m水深時小約5.96%。

在彈速范圍內，渦輪轉速與彈速接近為線性關系。對水深10m條件下的彈速與渦輪轉速關系曲線采用最小二乘法進行線性擬合，得：y＝230.8x－1454，非線性度γ＝0.998。

4 結 論

由于對高速下渦輪表面的空化情況、彈丸頭部阻力系數以及渦輪轉速規(guī)律沒有定量研究結論可供參考，一直以來水中彈藥引信渦輪設計及系統參數選取均以經驗設計為主。本文通過對某水下火箭彈引信渦輪高速旋轉兩相流動的三維流場數值仿真，得到以下結論：

（1）引信渦輪表面空化程度隨彈速的增加而嚴重，空化尺度與彈速成正比例關系，彈速30m／s左右時空化初生；

（2）水中彈藥的粘性阻力比空氣中大很多，彈速100m／s下頭部粘性阻力占總阻力的18.7%；

（3）高速下渦輪轉速與彈速仍具有較好的線性關系，10m水深條件下的R－v曲線非線性度γ＝0.998；

（4）不同水深條件對渦輪轉動有一定影響，文中R－v曲線對于安全系統參數設計、控制解保距離具有一定參考意義。

［1］SHEN D Z，ZHANG H，LI H J.A delay arming technical scheme for underwater rocket fuze based on external conformal turbine［C］.The 2ndInternational Conference on Computer Application and System Modeling，Xiamen，China.2011

［2］XU C J，ZHANG H.Numerical simulation 3－D of flow－field of a turbine alternator with side intake ducts for the fuze coupled［J］.Journal of Projectiles，Rockets，Missiles and Guidance，2006，26（1）：110－113.（in Chinese）徐長江，張合.引信用側進氣渦輪電機三維流場數值模擬［J］.彈箭與制導學報，2006，26（1）：110－113.

［3］WANG Z，ZHANG T，ZHENG D D.Comparison of different turbulence models for numerical simulation of internal flow field of turbine flow sensor［J］.Journal of Tianjin University，2007，40（12）：1447－1451.（in Chinese）王振，張濤，鄭丹丹.渦輪流量傳感器內部流場數值模擬中湍流模型比較［J］.天津大學學報，2007，40（12）：1447－1451.

［4］HUANG J F，ZHANG L X，YAO J，et al.Numerical simulation of 3D cavitation turbulent flow in francis hydro－turbine passage on mixture modeling technology［J］.Proceedings of the CSEE，2011，31（32）：115－121.（in Chinese）黃劍峰，張立翔，姚激，等.混流式水輪機三維空化湍流場混合數值模擬［J］.中國電機工程學報，2011，31（32）：115－121.

［5］LIU S H，ZHANG L，NISHI M，et al.Cavitating turbulent flow simulation in a francis turbine based on mixture model［J］.Journal of Fluids Engineering，2009，131（5）.

［6］WU Y L，LIU S H，DOU H S，et al.Simulations of unsteady cavitating turbulent flow in a Francis turbine using the RANS method and the improved mixture model of two－phase flows［J］.Engineering with Computers，2011，27（3）：235－250.

［7］SANSONE E，PELLONE C，MAITRE T.Modeling the unsteady cavitating flow in a cross－flow water turbine［J］.Journal of Fluids Engineering，2010，132（7）.

［8］SCHNERR G H，SAUER J.Physical and numerical modeling of unsteady cavitation dynamics［C］.Fourth International Conference on Multiphase Flow，New Orleans.USA，2001.

［9］LIU Z J，HE S Q，LIU H.Antifriction bearing application［M］.Beijing：Mechanic Industry Press，2007，3：650－654.（in Chinese）劉澤九，賀士荃，劉暉.滾動軸承應用［M］.北京：機械工業(yè)出版社，2007，3：650－654.

［10］ZHAO X D.Turbine flowmeter mathematical model and optimal design［J］.Academic Journal of Shanghai Mechanic College.1985，（2）：1－15.（in Chinese）趙學端.渦輪流量計數學模型與優(yōu)化設計［J］.上海機械學院學報，1985，（2）：1－15.

［11］XIA Z F，LUO S，YANG Y.Numerical simulations of propeller slipstream flows using actuator disk theory［J］.ACTA Aerodynamica Sinica，2012，30（2）：219－222.（in Chinese）夏貞鋒，羅淞，楊永.基于激勵盤理論的螺旋槳滑流數值模擬研究［J］.空氣動力學學報，2012，30（2）：219－222.

［12］ZHU M H，YE Z Y，JIN L.Numberical simulation of separation equilibrium of flow around 2D wing configuration at high angle of attack［J］.ACTA Aerodynamica Sinica，2012，30（4）：477－482.（in Chinese）祝明紅，葉正寅，金玲.二元翼型大迎角繞流的平衡態(tài)（解）的數值研究［J］.空氣動力學學報，2012，30（4）：477－482.

［13］WANG Z，ZHANG T，SUN L J.Numerical simulation on turbine flowmeter for measurement of gas－liquid two－phase flow and experiment verification［J］.Journal of Tianjin University，2008，41（11）：1303－1308.（in Chinese）王振，張濤，孫立軍.渦輪流量傳感器測量氣液兩相流的數值仿真與實驗驗證［J］.天津大學學報，2008，41（11）：1303－1308.

［14］WANG Z，ZHANG T，XU Y.Simulation and experiment on tangential type turbine flow sensor for measurement of low flow rate［J］.Journal of Tianjin University，2007，40（9）：1048－1053.（in Chinese）王振，張濤，徐英.測量小流量的切向式渦輪流量傳感器的仿真與實驗［J］.天津大學學報，2007，40（9）：1048－1053.

［15］WU Y D，OUYANG H，XU K，et al.Numerical investigation on the mechanism of micro tab to the airfoil［J］.ACTA Aerodynamica Sinica，2012，30（6）：786－791.（in Chinese）吳亞東，歐陽華，許坤，等.基于數值方法研究帶小插片的翼型的流動機理［J］.空氣動力學學報，2012，30（6）：786－791.