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

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

火炮會在射擊瞬間產(chǎn)生巨大的后坐沖擊力，制退機(jī)是消耗后坐能量、控制平穩(wěn)后坐的關(guān)鍵部件。后坐沖擊時(shí)間僅為0.1～0.2s，后坐速度可達(dá)到約10m／s［1］。在此工況下，制退機(jī)內(nèi)部液體呈三維非定常湍流流動(dòng)，并存在運(yùn)動(dòng)界面和真空中的射流［2］，流動(dòng)現(xiàn)象復(fù)雜，研究其內(nèi)部流場規(guī)律非常困難。

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

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

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

本文中基于制退機(jī)內(nèi)部流場三維模型，利用動(dòng)網(wǎng)格和滑移網(wǎng)格技術(shù)，對后坐沖擊過程中制退機(jī)內(nèi)部湍流流場進(jìn)行數(shù)值模擬。為檢驗(yàn)數(shù)值模擬的準(zhǔn)確性，分別采用3種不同的k－ε湍流模型進(jìn)行計(jì)算，并將數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行對比，以期為制退機(jī)內(nèi)部流場計(jì)算的湍流模型選擇提供參考。

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

假設(shè)制退機(jī)內(nèi)部液體進(jìn)行不可壓縮三維非定常湍流流動(dòng)，在直角坐標(biāo)系下建立矢量形式的流體連續(xù)性方程和動(dòng)量方程分別為

式中：ρ為流體密度，t為時(shí)間，u為速度矢量，p 為壓強(qiáng)，xi、xj、xk為坐標(biāo)分量，Si、Sj、Sk分別為廣義源項(xiàng)分量，μeff為有效粘性系數(shù)。μeff＝μ＋μt，μ為分子粘性系數(shù)，μt為湍動(dòng)粘度。

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

1.1 標(biāo)準(zhǔn)k－ε模型

標(biāo)準(zhǔn)k－ε模型［9］是目前使用最廣的湍流模型，具有適用范圍廣、計(jì)算精度合理的優(yōu)點(diǎn)，是一種針對高Re數(shù)的湍流計(jì)算模型。當(dāng)遇到彎曲壁面流動(dòng)、強(qiáng)旋流和逆壓梯度較大的問題時(shí)，計(jì)算精度會降低。當(dāng)流動(dòng)不可壓，且不考慮源項(xiàng)時(shí)，定義k與ε的輸運(yùn)方程為

式中：Gk為湍動(dòng)能k的產(chǎn)生項(xiàng)，各常數(shù)C1ε＝1.44，C2ε＝1.92，Cμ＝0.09，σk＝1.0，σε＝1.3。

1.2 RNGk－ε模型

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

1.3 Realizable k－ε 模型

針對標(biāo)準(zhǔn)模型對時(shí)均應(yīng)變率很大時(shí)，可能導(dǎo)致負(fù)的正應(yīng)力的情況，Realizable k－ε模型［11］取消了湍動(dòng)粘度計(jì)算中Cμ的常數(shù)取值，引入相關(guān)的旋轉(zhuǎn)和曲率，使其與應(yīng)變率相關(guān)聯(lián)。Realizable k－ε模型定義k與ε的輸運(yùn)方程為

2 建模、網(wǎng)格劃分與數(shù)值求解

根據(jù)某型火炮制退機(jī)的實(shí)際結(jié)構(gòu)，建立了三維計(jì)算模型，模型參數(shù)采用真實(shí)尺寸，如圖2所示。

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

2.1 網(wǎng)格劃分

因制退機(jī)內(nèi)部結(jié)構(gòu)復(fù)雜，既有規(guī)則的圓柱形內(nèi)腔，又有傾斜的流液孔和變截面的環(huán)形縫隙。因此，合理的網(wǎng)格劃分將是提高流場計(jì)算準(zhǔn)確度的重要前提。本文中采用分塊對接的網(wǎng)格劃分方法，將制退機(jī)內(nèi)部分為若干區(qū)域，對形狀簡單規(guī)則、會產(chǎn)生大變形的區(qū)域，如工作腔、非工作腔、復(fù)進(jìn)節(jié)制腔等運(yùn)用結(jié)構(gòu)網(wǎng)格進(jìn)行劃分。對形狀復(fù)雜、不涉及變形的區(qū)域，如制退桿活塞、調(diào)速筒等運(yùn)用適應(yīng)性更強(qiáng)的非結(jié)構(gòu)網(wǎng)格進(jìn)行劃分。通過滑移網(wǎng)格界面可實(shí)現(xiàn)運(yùn)動(dòng)區(qū)域與靜止區(qū)域間的數(shù)據(jù)實(shí)時(shí)傳遞?？紤]到制退機(jī)內(nèi)部結(jié)構(gòu)為中心軸對稱，為方便計(jì)算，計(jì)算模型選為整體結(jié)構(gòu)的1／2，計(jì)算網(wǎng)格數(shù)為365 842。

2.2 網(wǎng)格運(yùn)動(dòng)方法

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

2.3 數(shù)值求解方法

利用有限體積法將控制方程離散為計(jì)算模型網(wǎng)格單元上的代數(shù)方程，擴(kuò)散項(xiàng)、壓力項(xiàng)采用二階中心差分格式離散，動(dòng)量項(xiàng)采用一階迎風(fēng)格式離散，湍動(dòng)能、湍流耗散率采用二階迎風(fēng)格式離散，近壁區(qū)域采用增強(qiáng)型壁面函數(shù)方法處理，壓力速度耦合方程采用SIMPLEC算法求解。

3 計(jì)算結(jié)果與分析

采用標(biāo)準(zhǔn)k－ε模型、RNGk－ε模型、Realizable k－ε模型分別對制退機(jī)內(nèi)部流場進(jìn)行數(shù)值計(jì)算，得到工作腔、復(fù)進(jìn)節(jié)制腔的壓力數(shù)據(jù)，并與實(shí)驗(yàn)測試的壓力結(jié)果進(jìn)行了對比，如圖3所示。

從圖3可以看出，后坐過程中制退機(jī)內(nèi)的工作腔和復(fù)進(jìn)節(jié)制腔壓力的計(jì)算結(jié)果均較好地反映出實(shí)際壓力的變化規(guī)律，采用標(biāo)準(zhǔn)k－ε模型的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果具有更好的總體符合度，其最高壓力點(diǎn)也比其余2種模型的計(jì)算結(jié)果更接近實(shí)驗(yàn)值，而RNGk－ε模型和Realizable k－ε模型的計(jì)算結(jié)果基本一致，但與實(shí)驗(yàn)結(jié)果的誤差大于標(biāo)準(zhǔn)k－ε模型的，并且相對誤差隨壓力的增減而相應(yīng)地增減。

圖3 壓力的計(jì)算值與實(shí)驗(yàn)值的比較Fig.3 Pressures calculated by different k－εmodels compared with experimental results

后坐過程中，在0～13.5ms是彈丸膛內(nèi)運(yùn)動(dòng)期，即后坐加速運(yùn)動(dòng)期，3種模型計(jì)算結(jié)果與實(shí)驗(yàn)值符合較好；在13.5～116.2ms的后效期及慣性后坐期，即后坐減速運(yùn)動(dòng)期，標(biāo)準(zhǔn)k－ε模型的計(jì)算結(jié)果最接近實(shí)驗(yàn)值，其他2種模型的計(jì)算誤差明顯增大；在116.2～172.0ms的后坐結(jié)束期，速度變化較為平穩(wěn)，3種模型的計(jì)算值與實(shí)驗(yàn)值基本吻合，誤差都較小。表1給出了工作腔、復(fù)進(jìn)節(jié)制腔壓力計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的誤差。經(jīng)對比分析可知，采用標(biāo)準(zhǔn)k－ε模型的數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合最好。

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

4 結(jié) 論

（1）利用動(dòng)網(wǎng)格與滑移網(wǎng)格技術(shù)，對后坐沖擊過程中，分別采用標(biāo)準(zhǔn)k－ε模型、RNGk－ε模型和Realizable k－ε模型對真實(shí)結(jié)構(gòu)的制退機(jī)內(nèi)部三維湍流運(yùn)動(dòng)流場進(jìn)行了數(shù)值模擬，結(jié)果與實(shí)驗(yàn)曲線的變化趨勢基本一致，較好地反映了制退機(jī)內(nèi)部壓力的變化規(guī)律。

（2）對比分析工作腔、復(fù)進(jìn)節(jié)制腔壓力計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果，得出在彈丸膛內(nèi)運(yùn)動(dòng)期與后坐結(jié)束期，3種模型的壓力計(jì)算值均與實(shí)驗(yàn)值較好吻合，誤差較小，而在后效期與慣性后坐期，標(biāo)準(zhǔn)k－ε模型計(jì)算值最接近實(shí)驗(yàn)值，誤差不超過5%，小于其他2種模型。

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

［2］高樹滋，陳運(yùn)生，張?jiān)铝?，?火炮反后坐裝置設(shè)計(jì)［M］.北京：兵器工業(yè)出版社，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］高樂南.火炮駐退機(jī)湍流流場的預(yù)測分析［D］.南京：華東工學(xué)院，1990.

［5］陳剛.駐退機(jī)內(nèi)不可壓湍流流動(dòng)研究［D］.南京：華東工學(xué)院，1991.

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

［7］鄭建國.火炮制退機(jī)流場的數(shù)值模擬［J］.力學(xué)與實(shí)踐，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］.船舶力學(xué)，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］張志榮，冉景煜，張力，等.內(nèi)燃機(jī)缸內(nèi)CFD瞬態(tài)分析中動(dòng)態(tài)網(wǎng)格劃分技術(shù)［J］.重慶大學(xué)學(xué)報(bào)：自然科學(xué)版，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.