张弘钧， 钱林方， 陈光宋， 许彬

(1.南京理工大学 机械工程学院, 江苏 南京 210094; 2.中国北方工业公司 军贸技术研究院, 北京 100053)

基于径向基点插值法的大口径火炮身管固有频率特性研究

张弘钧1， 钱林方1， 陈光宋1， 许彬2

(1.南京理工大学 机械工程学院, 江苏 南京 210094; 2.中国北方工业公司 军贸技术研究院, 北京 100053)

考虑炮口制退器和炮尾的影响，建立身管弹性动力学方程，结合有限元法与无网格方法的各自优点，基于径向基点插值法。以积分点所在的单元支持域作为积分点的插值域，提出单元支持域径向基点插值法。建立离散形式的身管弹性动力学方程，计算获得身管固有频率。通过数值算例和实验验证了所提方法求解身管模态的有效性。计算结果表明：在相同网格数的情况下，所提方法能够得到更准确的结果，而且与实验结果相吻合。

兵器科学与技术； 身管； 固有频率； 单元支持域径向基点插值法

0 引言

身管为火炮的重要部件之一，其上固接有炮口制退器、炮尾等部件。在火炮发射过程中，身管作为弹丸在膛内运动的直接载体，其动力学特性直接影响弹丸在膛内的运动状态，最终影响弹丸出炮口瞬间的状态参数。通常，在身管优化设计过程中以提高身管的1阶频率为目标来提高身管的刚度，在动力学分析过程中，身管的模态分析结果可为身管的柔性动力学响应或刚柔耦合动力学的计算提供依据。目前，对火炮身管固有频率的分析方法有解析法、半解析法、有限元方法和实验方法。陈光宋等[1]基于Chebyshev正交多项式展开，研究了身管固有振动的半解析解法，并与有限元法得到的结果进行了验证。在此基础上，Qian等[2]采用谱单元法研究了弹炮耦合问题的不确定传播。刘军等[3]考虑了温度的影响，采用有限元法研究了温度对火炮身管固有频率和振型的影响。王宝元等[4]采用实验的方法，对自行火炮固有频率的特性进行了研究。此外，文献[5-8]关于身管的固有振动也进行了理论和实验研究。上述研究成果为分析身管固有频率提供了有益的参考价值，然而当采用有限元法进行身管模态分析或动力学计算时，为提高计算精度，需要将身管划分成较多的网格，大大降低了计算效率。

利用无网格法建立系统平衡方程时，不需要在求解区域和边界上划分单元，只需在求解区域内和边界上配置适当的节点，不仅克服了有限元法复杂的前处理过程，而且避免了由于单元质量导致计算精度的降低[9]。径向基点插值法作为一种新型无网格方法，在构造形函数的过程中，为了避免单纯采用多项式基可能导致的刚度矩阵奇异，引入径向基函数并采用更多的节点插值，通常具有高阶连续性，从而提高了计算精度[10]。为了克服无网格法在单元矩阵积分过程中计算效率低的缺点，相关学者做了大量的改进，包括自然领域径向基点插值法[11]、光滑径向基点插值法[12]、谱径向基点插值法[13]、比例边界径向基点插值法[14]等。

本文基于无网格径向基点插值法，提出单元支持域作为积分点的支持域，以改进径向基点插值法形函数的构造，建立身管三维柔性动力学模型，进而以此分析身管的固有频率，并通过数值和实验验证本文方法的正确性和高效性，为身管设计和火炮发射动力学分析提供一定的技术参考。

1 身管振动控制微分方程

如图1所示：考虑火炮身管、炮口制退器、炮尾的综合模型，炮尾通过圆柱来等效，长度和质量与炮尾的质量一致；炮口制退器通过中空圆柱形质量块来等效，长度和质量与炮口制退器相同；与身管端面固结的截面内外径与身管前端面的内外径一致。图1中：mb和mm分别表示炮尾和炮口制退器质量；D为身管内径；Ls为身管长度；Lx1和Lx2分别为摇架前后衬套至身管尾端面的距离；Lb和Lm分别为炮尾和炮口制退器长度；Db和Dm分别为炮尾和炮口制退器的等效尺寸；其余结构尺寸如图所示。

1)几何方程

(1)

2)物理方程

σij=Dijklεkl，x∈Ω；

(2)

3)平衡方程

(3)

4)边界条件

面力边界条件

(4)

位移边界条件

(5)

5)初始条件

位移初始条件

ui(x,t0)=ui0(x)，x∈Ω；

(6)

速度初始条件

(7)

6)弹性体的总位能Πp为

(8)

2 单元支持域径向基点插值法

通常，采用径向基点插值法[15]进行平衡方程弱形式积分时，需要根据节点的支持域搜索每一个积分点对应的邻域节点，进而根据这些邻域节点构造插值函数，这样不仅对每个积分节点都需要搜索邻域节点，而且对每个积分点都要计算一次插值函数。如图2所示，点pA和pB的支持域分别为ΩA和ΩB. 不仅如此，由于径向基函数是高阶连续函数，采用数值积分获得单元矩阵的过程中需要采用高阶积分方案(即增加积分点)，从而大大增加了计算量，这也是无网格方法计算效率低的一个重要因素。为此，本文提出单元支持域插值方法，单元支持域是指节点所在的单元及其邻接单元构成的区域，根据单元支持域内的节点来构造插值函数，如图3所示，e1～e9为单元编号，点pA和pB的支持域均为单元e5，单元e5的支持域为Ωs.

利用支持域Ωs上的n个节点值和径向基- 多项式基进行插值来近似函数u(x)，u(x)的近似函数以uh(x)表示，则

(9)

式中：n和m分别为径向基函数和多项式基函数的个数；ai为径向基函数Ri(x)的系数；bj为多项式基函数pj(x)的系数；a和b为对应的系数向量；R(x)和p(x)分别为径向基向量和多项式基向量，对于二维问题，x={x,y}T.

将x的支持域上n个节点的坐标值代入(9)式中，得到：

(10)

式中：uk为第k个节点的位移；Ri(xk,yk)和pj(xk,yk)分别为径向基函数和多项式基函数在节点(xk,yk)的值。

将(10)式写成矩阵形式，有

u=Ra+Pb，

(11)

式中：u为x点的支持域上n个节点值的向量；R和P分别为径向基和多项式基的力矩矩阵，

(12)

(13)

式中：

(14)

由(13)式可知R为对称矩阵：

Ri(rj)=Rj(ri).

(15)

(9)式中有n+m个系数，而(11)式只有n个方程，为求出n+m个系数，本文根据Golberg等[16]的建议补充如下m个方程：

(16)

(16)式改写成矩阵形式表示为

PTa=0.

(17)

由(11)式和(17)式可得

a=Sau，

(18)

b=Sbu，

(19)

式中：

Sb=[PTR-1P]PTR-1，

(20)

Sa=R-1[I-PSb)]=R-1-R-1PSb，

(21)

其中I为n×n阶的单位矩阵。

上式PTR-1P中的每个矩阵都是常数矩阵，因此只需计算PTR-1一次并保存即可。

最后，将(18)式和(19)式代入(9)式中，就可得到近似函数的表达式：

uh(x)=[RT(x)Sa+pT(x)Sb]u=Φ(x)u，

(22)

式中：Φ(x)是包含n个节点的形函数矩阵，

(23)

第k个节点的形函数为

(24)

从而很容易求得形函数的导数，即

(25)

3 身管控制微分方程的离散

将身管三维空间位移场分别用插值函数近似，并写成矩阵形式有

(26)

式中:

(27)

由此可得到任一点x处的应变为

(28)

式中：

(29)

任一点的应力为

(30)

式中：弹性矩阵D的表达式为

(31)

其中λ为拉梅常数，λ=Eμ/((1+μ)(1-2μ))，E和μ分别为材料的弹性模量和泊松比。

将(26)式～(31)式代入(8)式，得到身管离散动力学方程为

(32)

(33)

(34)

(35)

(36)

如果不考虑(32)式中的阻尼和载荷的影响，则可以得到身管自由振动系统方程，即

(37)

(-ω2M+K)q=0.

(38)

(38)式为广义特征值问题，求解(38)式即可得到身管的固有圆频率ω和振型q.

本文利用有限元前处理软件划分网格，通过编程实现上述计算过程，计算流程如图4所示。

4 算例

4.1 悬臂梁

考虑如图5所示的悬臂梁，其长和高分别为L和H，梁的左端约束如图5所示，右端作用载荷F. 相关的计算参数如下：F=-1 000 kN, 弹性模量E=3.0×1010N/m2,H=12 m,L=48 m, 泊松比v=0.3. 计算过程中左端通过(39)式和(40)式进行约束，右端通过(41)式～(43)式施加载荷。则平面应力静态问题的解析解为

(39)

(40)

(41)

σy=0，

(42)

(43)

式中：I为梁的截面惯性矩，I=H3/12；x和y分别为梁上的坐标；ux和uy分别为梁x和y方向的位移；σx和σy分别为x和y方向的应力；τxy为剪应力。

为研究本文方法的收敛性，考虑5种不同的网格分布2×2，4×4，8×8，16×16，32×32. 本文方法的计算结果以及有限元结果和解析解的对比如表1所示。表1中的计算结果显示，本文方法具有较高的收敛率和精度。

4.2 身管固有频率及数值验证

以某155 mm火炮身管为例进行数值计算对比验证，身管的基本物理参数为：弹性模量E=2.1×1011N/m2，剪切模量G=1/2.6E，身管密度ρ=7 850 kg/m3. 考察身管边界条件为自由状态，将本文方法和有限元方法进行对比，其中采用本文方法时将身管离散为11 342个六面体网格，在有限元计算中分别采用11 342个网格和通过增加网格数来验证收敛的355 744个网格。本文和有限元两种计算方法得到的身管前10阶固有频率fi=ωi/2π结果对比见表2. 从表2可以看出：本文方法的计算结果与有限元软件Abaqus(355 744个六面体网格)得到的结果相吻合，网格数量为其3%左右；在相同网格数的情况下，本文方法比有限元法能够得到更准确的结果。

4.3 身管固有频率实验及验证

同样以上述某155 mm火炮身管为例，考虑炮口制退器和炮尾结构的复杂性及较多的零部件(例如炮闩等)非固接关系影响模态实验的准确性，去除炮尾和炮口制退器，并在距身管尾端面和炮口端面1.5 m处用枕木支撑，进行身管模态实验，得到的前10阶身管固有频率及本文的计算结果如表3所示。表3的结果显示：本文的计算结果与实验结果吻合较好，考虑到在进行模态实验的过程中身管是通过枕木支撑的，并非完全自由状态，因此实验结果与计算结果存在误差，但并无不合理的偏差。另外，将本文方法与有限元方法计算得到的身管固有频率进行了比较，结果也比较吻合。

5 结论

本文利用有限元法中的网格，提出单元支持域径向基点插值法，通过标准算例检验了该方法的收敛性和精度。对本文计算的身管固有频率进行了数值验证和实验验证，结果显示在相同网格情况下，本文方法能够得到更准确的结果，而且与实验结果相吻合。

)

[1] 陈光宋, 钱林方, 徐亚栋,等. 身管横向固有振动的半解析解法[J]. 兵工学报, 2012, 33(10):1168-1172.

CHEN Guang-song, QIAN Lin-fang, XU Ya-dong, et al. Semi-analytic solution of nature frequency of transverse vibration of a barrel[J]. Acta Armamentarii, 2012, 33(10): 1168-1172.(in Chinese)

[2] Qian L F, Chen G S. The uncertainty propagation analysis of the projectile-barrel coupling problem[J]. Defence Technology，2017,13(4):1-5.

[3] 刘军, 卫丰, 苟文选,等. 温度场下火炮身管固有频率和振型的三维有限元分析[J]. 火炮发射与控制学报, 2004, 25(4):5-7.

LIU Jun, WEI Feng, GOU Wen-xuan, et al. 3D FEM analysis of natural frequency and mode of gun tube in temperature field[J]. Journal of Gun Launch & Control, 2004, 25(4):5-7.(in Chinese)

[4] 王宝元, 刘朋科, 衡刚,等. 自行火炮固有频率特性实验研究[J]. 兵工学报, 2011, 32(11):1315-1319.

WANG Bao-yuan, LIU Peng-ke, HENG Gang, et al. Experimental research on natural frequency characteristics of self-propelled guns[J]. Acta Armamentarii, 2011, 32(11):1315-1319.(in Chinese)

[5] 张海航, 狄长安. 某型多管火炮的振动模态分析[J]. 兵工学报, 2008, 29(12):1514-1517.

ZHANG Hai-hang, DI Chang-an. Analysis of vibration modal of a multibarrel cannon [J]. Acta Armamentarii, 2008, 29(12):1514-1517.(in Chinese)

[6] 陈光宋, 钱林方, 吉磊. 身管固有频率高效全局灵敏度分析[J]. 振动与冲击, 2015, 34(21):31-36.

CHEN Guang-song, QIAN Lin-fang, JI Lei. An effective global sensitivity analysis method for natural frequencies of a barrel[J]. Journal of Vibration and Shock, 2015, 34(21):31-36.(in Chinese)

[7] 苏忠亭, 徐达, 杨明华,等. 基于模态试验的某火炮身管有限元模型修正[J]. 振动与冲击, 2012, 31(24):54-59.

SU Zhong-ting, XU Da, YANG Ming-hua, et al. Finite-element model updating for a gun barrel based on modal test[J]. Journal of Vibration and Shock, 2012, 31(24):54-59.(in Chinese)

[8] 吴东亚, 邢宏光, 崔军. 基于虚拟样机技术的某型坦克炮炮身模态分析[J]. 火力与指挥控制, 2011, 36(7):65-67.

WU Dong-ya, XING Hong-guang, CUI Jun. Modal analysis on a certain tank gun barrel based on virtual prototyping technology[J]. Fire Control and Command Control, 2011, 36(7): 65-67.(in Chinese)

[9] 龙述尧.无网格方法及其在固体力学中的应用[M]. 北京：科学出版社, 2014:1-9.

LONG Shu-yao. Meshless methods and its applications on solid mechanics[M]. Beijing: Science Press, 2014:1-9. (in Chinese)

[10] 杜超凡,章定国,洪嘉振.径向基点插值法在旋转柔性梁动力学中的应用[J].力学学报,2015,47(2):279-288.

DU Chao-fan, ZHANG Ding-guo, HONG Jia-zhen. A meshfree method based on radial point interpolation method for the dynamic analysis of rotating flexible beams[J]. Chinese Journal of Theoretical and Applied Mechanics, 2015,47(2):279-288. (in Chinese)

[11] Dinis L M J S, Jorge R M N, Belinha J. Analysis of 3D solids using the natural neighbour radial point interpolation method[J]. Computer Methods in Applied Mechanics & Engineering, 2007, 196(13/14/15/16):2009-2028.

[12] Feng S Z, Cui X Y, Chen F, et al. An edge/face-based smoothed radial point interpolation method for static analysis of structures[J]. Engineering Analysis with Boundary Elements, 2016, 68:1-10.

[13] Shivanian E. A new spectral meshless radial point interpolation (SMRPI) method: a well-behaved alternative to the meshless weak forms[J]. Engineering Analysis with Boundary Elements, 2015, 54:1-12.

[14] Hajiazizi M, Graili A. A scaled boundary radial point interpolation method for 2D elasticity problems[J]. International Journal for Numerical Methods in Engineering, 2017, 112(7): 832-851.

[15] Wang J G, Liu G R. A point interpolation meshless method based on radial basis functions[J]. International Journal for Numerical Methods in Engineering, 2002, 54(11): 1623-1648.

[16] Golberg M A, Chen C S, Bowman H. Some recent results and proposals for the use of radial basis functions in the BEM[J]. Engineering Analysis with Boundary Elements, 1999, 23(4): 285-296.

ResearchontheNaturalFrequenciesofaLarge-caliberHowitzerBarrelBasedonRadialPointInterpolationMethod

ZHANG Hong-jun1， QIAN Lin-fang1， CHEN Guang-song1， XU Bin2

(1.School of Mechanical Engineering, Nanjing University of Science and Technology, Nanjing 210094，Jiangsu，China； 2.Technical Research Institute of Military Trade, China North Industries Corporation, Beijing 100053，China)

Considering the influence of muzzle brake and breech ring, an elastic dynamics equation of the barrel is established. Combining the advantages of the finite element method and the meshless method, an element supported domain radial point interpolation method is proposed based on the radial basis point interpolation method, in which the element supported domain is treated as the interpolation domain of the integration point to establish the discrete elastic dynamics equation. The natural frequencies of barrel are estimated by solving the dynamic equation. Finally, the effectiveness of the proposed method is demonstrated by numerical examples and experiments, and the results show that the proposed method can provide more accurate results with the same element number，which can coincide with the experimental results.

ordnance science and technology; barrel; natural frequency; element supported domain radial point interpolation method

2017-05-05

国家自然科学基金项目(11472137)

张弘钧(1977—), 男， 副研究员, 博士研究生。 E-mail: zhanghj@njust.edu.cn

钱林方(1961—), 男, 教授， 博士生导师。 E-mail: lfqian@vip.163.com

TJ302

A

1000-1093(2017)12-2321-07

10.3969/j.issn.1000-1093.2017.12.004