k－ε雙方程湍流模型對制退機內流場計算的適用性分析＊

張曉東，張培林，傅建平，王 成，楊玉棟

（軍械工程學院一系，河北 石家莊 050003）

火炮會在射擊瞬間產生巨大的后坐沖擊力，制退機是消耗后坐能量、控制平穩后坐的關鍵部件。后坐沖擊時間僅為0.1～0.2s，后坐速度可達到約10m／s［1］。在此工況下，制退機內部液體呈三維非定常湍流流動，并存在運動界面和真空中的射流［2］，流動現象復雜，研究其內部流場規律非常困難。

隨著計算機技術、數值計算技術和湍流理論的快速發展，數值模擬手段已逐漸被應用于制退機內部的復雜流動研究。從制退機局部結構到簡化的二維模型，從層流模型到非定常湍流模型的計算［3－7］都取得了很多有益的結論。但到目前為止，數值模擬結果與實驗結果仍存在較大的誤差，主要是計算的建模誤差所致［8］。究其原因，一是模型過于簡化，二是受計算理論和計算條件所限。

基于Boussinesq假設的k－ε雙方程湍流模型是當前流體機械內部流場計算最常的用湍流模型，主要有標準k－ε 模型［9］，以及改進的 RNGk－ε 模型［10］、Realizable k－ε 模型［11］等。高樂南［4］、陳剛［5］利用標準k－ε模型對制退機簡化模型的定常湍流流動進行了數值模擬。吳利民［6］利用RNGk－ε模型對制退機簡化模型的內部湍流流場進行了數值模擬。以上采用的部分簡化模型如圖1所示。

圖1 2種簡化的制退機模型Fig.1 Two simplified models of recoil brake

本文中基于制退機內部流場三維模型，利用動網格和滑移網格技術，對后坐沖擊過程中制退機內部湍流流場進行數值模擬。為檢驗數值模擬的準確性，分別采用3種不同的k－ε湍流模型進行計算，并將數值模擬結果與實驗結果進行對比，以期為制退機內部流場計算的湍流模型選擇提供參考。

1 基本控制方程與湍流模型

假設制退機內部液體進行不可壓縮三維非定常湍流流動，在直角坐標系下建立矢量形式的流體連續性方程和動量方程分別為

式中：ρ為流體密度，t為時間，u為速度矢量，p 為壓強，xi、xj、xk為坐標分量，Si、Sj、Sk分別為廣義源項分量，μeff為有效粘性系數。μeff＝μ＋μt，μ為分子粘性系數，μt為湍動粘度。

在渦粘模型方法中，不直接處理應力項，而是把湍流應力表示成湍動粘度的函數其中Cμ為計算常數。對以上控制方程進行時間平均，引入Reynolds應力項，就可以基于k－ε湍流模型建立湍動能k與耗散率ε的輸運方程。再聯立μt與湍流模型，進而可對控制方程進行求解。

1.1 標準k－ε模型

標準k－ε模型［9］是目前使用最廣的湍流模型，具有適用范圍廣、計算精度合理的優點，是一種針對高Re數的湍流計算模型。當遇到彎曲壁面流動、強旋流和逆壓梯度較大的問題時，計算精度會降低。當流動不可壓，且不考慮源項時，定義k與ε的輸運方程為

式中：Gk為湍動能k的產生項，各常數C1ε＝1.44，C2ε＝1.92，Cμ＝0.09，σk＝1.0，σε＝1.3。

1.2 RNGk－ε模型

RNGk－ε模型［10］通過重正化群方法，修正了湍動粘度，考慮了平均流動中的旋轉及旋流流動情況，在ε方程中出現了考慮非平衡應變率影響的附加耗散生成項，反映了主流的時均應變率Eij，可以更好地處理具有強曲率影響和壁面約束的湍流分離流動。RNGk－ε模型定義k與ε的輸運方程為

1.3 Realizable k－ε 模型

針對標準模型對時均應變率很大時，可能導致負的正應力的情況，Realizable k－ε模型［11］取消了湍動粘度計算中Cμ的常數取值，引入相關的旋轉和曲率，使其與應變率相關聯。Realizable k－ε模型定義k與ε的輸運方程為

2 建模、網格劃分與數值求解

根據某型火炮制退機的實際結構，建立了三維計算模型，模型參數采用真實尺寸，如圖2所示。

圖2 制退機三維計算模型Fig.2 3－D computation model of recoil brake

2.1 網格劃分

因制退機內部結構復雜，既有規則的圓柱形內腔，又有傾斜的流液孔和變截面的環形縫隙。因此，合理的網格劃分將是提高流場計算準確度的重要前提。本文中采用分塊對接的網格劃分方法，將制退機內部分為若干區域，對形狀簡單規則、會產生大變形的區域，如工作腔、非工作腔、復進節制腔等運用結構網格進行劃分。對形狀復雜、不涉及變形的區域，如制退桿活塞、調速筒等運用適應性更強的非結構網格進行劃分。通過滑移網格界面可實現運動區域與靜止區域間的數據實時傳遞。考慮到制退機內部結構為中心軸對稱，為方便計算，計算模型選為整體結構的1／2，計算網格數為365 842。

2.2 網格運動方法

為真實模擬后坐沖擊過程中制退機內部液體流動，本文利用動態層變網格技術，通過非結構網格區域整體運動與結構網格區域內網格的增減，實現了在給定速度下的制退機內流場運動。結構網格更新采用動態層變法［12］，該方法根據與運動物面相鄰的網格層高度來決定增加或減少網格層數，動網格模型指定一個理想層高，鄰近動邊界的網格單元根據層高度來分裂出新單元層或與鄰近層合并。以復進節制腔區域為例，后坐過程中處于被拉伸狀態，則動邊界相鄰的單元層被擴展至hmin，hmin為單元層最小高度，滿足hmin＞ （1 ＋αs）hi，αs為分裂因子，hi為理想單元層高。若條件滿足，該層將根據指定條件（恒定高度或恒定比例）進行分裂，即在層j將分裂為2層：1層高為hi，1層高為（hmin–hi）。反之，被壓縮的單元層將與鄰近的單元層合并成1個新層。

2.3 數值求解方法

利用有限體積法將控制方程離散為計算模型網格單元上的代數方程，擴散項、壓力項采用二階中心差分格式離散，動量項采用一階迎風格式離散，湍動能、湍流耗散率采用二階迎風格式離散，近壁區域采用增強型壁面函數方法處理，壓力速度耦合方程采用SIMPLEC算法求解。

3 計算結果與分析

采用標準k－ε模型、RNGk－ε模型、Realizable k－ε模型分別對制退機內部流場進行數值計算，得到工作腔、復進節制腔的壓力數據，并與實驗測試的壓力結果進行了對比，如圖3所示。

從圖3可以看出，后坐過程中制退機內的工作腔和復進節制腔壓力的計算結果均較好地反映出實際壓力的變化規律，采用標準k－ε模型的計算結果與實驗結果具有更好的總體符合度，其最高壓力點也比其余2種模型的計算結果更接近實驗值，而RNGk－ε模型和Realizable k－ε模型的計算結果基本一致，但與實驗結果的誤差大于標準k－ε模型的，并且相對誤差隨壓力的增減而相應地增減。

圖3 壓力的計算值與實驗值的比較Fig.3 Pressures calculated by different k－εmodels compared with experimental results

后坐過程中，在0～13.5ms是彈丸膛內運動期，即后坐加速運動期，3種模型計算結果與實驗值符合較好；在13.5～116.2ms的后效期及慣性后坐期，即后坐減速運動期，標準k－ε模型的計算結果最接近實驗值，其他2種模型的計算誤差明顯增大；在116.2～172.0ms的后坐結束期，速度變化較為平穩，3種模型的計算值與實驗值基本吻合，誤差都較小。表1給出了工作腔、復進節制腔壓力計算結果與實驗結果的誤差。經對比分析可知，采用標準k－ε模型的數值計算結果與實驗結果吻合最好。

表1 不同湍流模型下的制退機內部壓力計算誤差Table1 Pressure computation errors by different turbulence models of recoil brake

4 結 論

（1）利用動網格與滑移網格技術，對后坐沖擊過程中，分別采用標準k－ε模型、RNGk－ε模型和Realizable k－ε模型對真實結構的制退機內部三維湍流運動流場進行了數值模擬，結果與實驗曲線的變化趨勢基本一致，較好地反映了制退機內部壓力的變化規律。

（2）對比分析工作腔、復進節制腔壓力計算結果與實驗結果，得出在彈丸膛內運動期與后坐結束期，3種模型的壓力計算值均與實驗值較好吻合，誤差較小，而在后效期與慣性后坐期，標準k－ε模型計算值最接近實驗值，誤差不超過5%，小于其他2種模型。

［1］二Ο二研究所.火炮射擊參數的常規測試［M］.北京：第五機械工業部，1980：147－149.

［2］高樹滋，陳運生，張月林，等.火炮反后坐裝置設計［M］.北京：兵器工業出版社，1995：54－63.

［3］Chen C J，Neshat H N，Li P.The Finite Analytic Method［M］.Iowa Institute of Hydraulic Research，The University of Iowa，1983.

［4］高樂南.火炮駐退機湍流流場的預測分析［D］.南京：華東工學院，1990.

［5］陳剛.駐退機內不可壓湍流流動研究［D］.南京：華東工學院，1991.

［6］吳利民.駐退機湍流場檢測技術研究及其數值模擬［D］.南京：南京理工大學，1999.

［7］鄭建國.火炮制退機流場的數值模擬［J］.力學與實踐，2001，23（2）：30－32.ZHENG Jian－guo.Numerical simulation of flow field in a gun recoil mechanism［J］.Mechanics in Engineering，2001，23（2）：30－32.

［8］張楠，沈泓萃，姚惠之.阻力和流場的CFD不確定度分析探討［J］.船舶力學，2008，12（2）：211－224.ZHANG Nan，SHEN Hong－cui，YAO Hui－zhi.Uncertainty analysis in CFD for resistance and flow field［J］.Journal of Ship Mechanics，2008，12（2）：211－224.

［9］Launder B E，Spalding D B.Lectures in mathematical models of turbulence［M］.London：Academic Press，1972.

［10］Yakhot V，Orzag S A.Renormalization group analysis of turbulence：Basic theory［J］.Journal of scientific Computing，1986，1（1）：3－11.

［11］Shih T H，Liou W W，Shabbir A，et al.A new k－εeddy viscosity model for high reynolds number turbulent flows－model development and validation［J］.Computers ＆ Fluids，1995，24（3）：227－238.

［12］張志榮，冉景煜，張力，等.內燃機缸內CFD瞬態分析中動態網格劃分技術［J］.重慶大學學報：自然科學版，2005，28（11）：97－100.ZHANG Zhi－rong，RAN Jing－yu，ZHANG Li，et al.Techniques of dynamic mesh in the CFD transient analysis of an engine［J］.Journal of Chongqing University：Natural Science Edition，2005，28（11）：97－100.