张晓东，张培林，傅建平，王 成，杨玉栋

（军械工程学院一系，河北 石家庄 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.