張弘鈞， 錢林方， 陳光宋， 許彬

(1.南京理工大學 機械工程學院, 江蘇 南京 210094; 2.中國北方工業(yè)公司 軍貿技術研究院, 北京 100053)

基于徑向基點插值法的大口徑火炮身管固有頻率特性研究

張弘鈞1， 錢林方1， 陳光宋1， 許彬2

(1.南京理工大學 機械工程學院, 江蘇 南京 210094; 2.中國北方工業(yè)公司 軍貿技術研究院, 北京 100053)

考慮炮口制退器和炮尾的影響，建立身管彈性動力學方程，結合有限元法與無網格方法的各自優(yōu)點，基于徑向基點插值法。以積分點所在的單元支持域作為積分點的插值域，提出單元支持域徑向基點插值法。建立離散形式的身管彈性動力學方程，計算獲得身管固有頻率。通過數(shù)值算例和實驗驗證了所提方法求解身管模態(tài)的有效性。計算結果表明：在相同網格數(shù)的情況下，所提方法能夠得到更準確的結果，而且與實驗結果相吻合。

兵器科學與技術； 身管； 固有頻率； 單元支持域徑向基點插值法

0 引言

身管為火炮的重要部件之一，其上固接有炮口制退器、炮尾等部件。在火炮發(fā)射過程中，身管作為彈丸在膛內運動的直接載體，其動力學特性直接影響彈丸在膛內的運動狀態(tài)，最終影響彈丸出炮口瞬間的狀態(tài)參數(shù)。通常，在身管優(yōu)化設計過程中以提高身管的1階頻率為目標來提高身管的剛度，在動力學分析過程中，身管的模態(tài)分析結果可為身管的柔性動力學響應或剛柔耦合動力學的計算提供依據。目前，對火炮身管固有頻率的分析方法有解析法、半解析法、有限元方法和實驗方法。陳光宋等[1]基于Chebyshev正交多項式展開，研究了身管固有振動的半解析解法，并與有限元法得到的結果進行了驗證。在此基礎上，Qian等[2]采用譜單元法研究了彈炮耦合問題的不確定傳播。劉軍等[3]考慮了溫度的影響，采用有限元法研究了溫度對火炮身管固有頻率和振型的影響。王寶元等[4]采用實驗的方法，對自行火炮固有頻率的特性進行了研究。此外，文獻[5-8]關于身管的固有振動也進行了理論和實驗研究。上述研究成果為分析身管固有頻率提供了有益的參考價值，然而當采用有限元法進行身管模態(tài)分析或動力學計算時，為提高計算精度，需要將身管劃分成較多的網格，大大降低了計算效率。

利用無網格法建立系統(tǒng)平衡方程時，不需要在求解區(qū)域和邊界上劃分單元，只需在求解區(qū)域內和邊界上配置適當?shù)墓?jié)點，不僅克服了有限元法復雜的前處理過程，而且避免了由于單元質量導致計算精度的降低[9]。徑向基點插值法作為一種新型無網格方法，在構造形函數(shù)的過程中，為了避免單純采用多項式基可能導致的剛度矩陣奇異，引入徑向基函數(shù)并采用更多的節(jié)點插值，通常具有高階連續(xù)性，從而提高了計算精度[10]。為了克服無網格法在單元矩陣積分過程中計算效率低的缺點，相關學者做了大量的改進，包括自然領域徑向基點插值法[11]、光滑徑向基點插值法[12]、譜徑向基點插值法[13]、比例邊界徑向基點插值法[14]等。

本文基于無網格徑向基點插值法，提出單元支持域作為積分點的支持域，以改進徑向基點插值法形函數(shù)的構造，建立身管三維柔性動力學模型，進而以此分析身管的固有頻率，并通過數(shù)值和實驗驗證本文方法的正確性和高效性，為身管設計和火炮發(fā)射動力學分析提供一定的技術參考。

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]進行平衡方程弱形式積分時，需要根據節(jié)點的支持域搜索每一個積分點對應的鄰域節(jié)點，進而根據這些鄰域節(jié)點構造插值函數(shù)，這樣不僅對每個積分節(jié)點都需要搜索鄰域節(jié)點，而且對每個積分點都要計算一次插值函數(shù)。如圖2所示，點pA和pB的支持域分別為ΩA和ΩB. 不僅如此，由于徑向基函數(shù)是高階連續(xù)函數(shù)，采用數(shù)值積分獲得單元矩陣的過程中需要采用高階積分方案(即增加積分點)，從而大大增加了計算量，這也是無網格方法計算效率低的一個重要因素。為此，本文提出單元支持域插值方法，單元支持域是指節(jié)點所在的單元及其鄰接單元構成的區(qū)域，根據單元支持域內的節(jié)點來構造插值函數(shù)，如圖3所示，e1～e9為單元編號，點pA和pB的支持域均為單元e5，單元e5的支持域為Ωs.

利用支持域Ωs上的n個節(jié)點值和徑向基- 多項式基進行插值來近似函數(shù)u(x)，u(x)的近似函數(shù)以uh(x)表示，則

(9)

式中：n和m分別為徑向基函數(shù)和多項式基函數(shù)的個數(shù)；ai為徑向基函數(shù)Ri(x)的系數(shù)；bj為多項式基函數(shù)pj(x)的系數(shù)；a和b為對應的系數(shù)向量；R(x)和p(x)分別為徑向基向量和多項式基向量，對于二維問題，x={x,y}T.

將x的支持域上n個節(jié)點的坐標值代入(9)式中，得到：

(10)

式中：uk為第k個節(jié)點的位移；Ri(xk,yk)和pj(xk,yk)分別為徑向基函數(shù)和多項式基函數(shù)在節(jié)點(xk,yk)的值。

將(10)式寫成矩陣形式，有

u=Ra+Pb，

(11)

式中：u為x點的支持域上n個節(jié)點值的向量；R和P分別為徑向基和多項式基的力矩矩陣，

(12)

(13)

式中：

(14)

由(13)式可知R為對稱矩陣：

Ri(rj)=Rj(ri).

(15)

(9)式中有n+m個系數(shù)，而(11)式只有n個方程，為求出n+m個系數(shù)，本文根據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中的每個矩陣都是常數(shù)矩陣，因此只需計算PTR-1一次并保存即可。

最后，將(18)式和(19)式代入(9)式中，就可得到近似函數(shù)的表達式：

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

(22)

式中：Φ(x)是包含n個節(jié)點的形函數(shù)矩陣，

(23)

第k個節(jié)點的形函數(shù)為

(24)

從而很容易求得形函數(shù)的導數(shù)，即

(25)

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

將身管三維空間位移場分別用插值函數(shù)近似，并寫成矩陣形式有

(26)

式中:

(27)

由此可得到任一點x處的應變?yōu)?/p>

(28)

式中：

(29)

任一點的應力為

(30)

式中：彈性矩陣D的表達式為

(31)

其中λ為拉梅常數(shù)，λ=Eμ/((1+μ)(1-2μ))，E和μ分別為材料的彈性模量和泊松比。

將(26)式～(31)式代入(8)式，得到身管離散動力學方程為

(32)

(33)

(34)

(35)

(36)

如果不考慮(32)式中的阻尼和載荷的影響，則可以得到身管自由振動系統(tǒng)方程，即

(37)

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

(38)

(38)式為廣義特征值問題，求解(38)式即可得到身管的固有圓頻率ω和振型q.

本文利用有限元前處理軟件劃分網格，通過編程實現(xiàn)上述計算過程，計算流程如圖4所示。

4 算例

4.1 懸臂梁

考慮如圖5所示的懸臂梁，其長和高分別為L和H，梁的左端約束如圖5所示，右端作用載荷F. 相關的計算參數(shù)如下：F=-1 000 kN, 彈性模量E=3.0×1010N/m2,H=12 m,L=48 m, 泊松比v=0.3. 計算過程中左端通過(39)式和(40)式進行約束，右端通過(41)式～(43)式施加載荷。則平面應力靜態(tài)問題的解析解為

(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 身管固有頻率及數(shù)值驗證

以某155 mm火炮身管為例進行數(shù)值計算對比驗證，身管的基本物理參數(shù)為：彈性模量E=2.1×1011N/m2，剪切模量G=1/2.6E，身管密度ρ=7 850 kg/m3. 考察身管邊界條件為自由狀態(tài)，將本文方法和有限元方法進行對比，其中采用本文方法時將身管離散為11 342個六面體網格，在有限元計算中分別采用11 342個網格和通過增加網格數(shù)來驗證收斂的355 744個網格。本文和有限元兩種計算方法得到的身管前10階固有頻率fi=ωi/2π結果對比見表2. 從表2可以看出：本文方法的計算結果與有限元軟件Abaqus(355 744個六面體網格)得到的結果相吻合，網格數(shù)量為其3%左右；在相同網格數(shù)的情況下，本文方法比有限元法能夠得到更準確的結果。

4.3 身管固有頻率實驗及驗證

同樣以上述某155 mm火炮身管為例，考慮炮口制退器和炮尾結構的復雜性及較多的零部件(例如炮閂等)非固接關系影響模態(tài)實驗的準確性，去除炮尾和炮口制退器，并在距身管尾端面和炮口端面1.5 m處用枕木支撐，進行身管模態(tài)實驗，得到的前10階身管固有頻率及本文的計算結果如表3所示。表3的結果顯示：本文的計算結果與實驗結果吻合較好，考慮到在進行模態(tài)實驗的過程中身管是通過枕木支撐的，并非完全自由狀態(tài)，因此實驗結果與計算結果存在誤差，但并無不合理的偏差。另外，將本文方法與有限元方法計算得到的身管固有頻率進行了比較，結果也比較吻合。

5 結論

本文利用有限元法中的網格，提出單元支持域徑向基點插值法，通過標準算例檢驗了該方法的收斂性和精度。對本文計算的身管固有頻率進行了數(shù)值驗證和實驗驗證，結果顯示在相同網格情況下，本文方法能夠得到更準確的結果，而且與實驗結果相吻合。

)

[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] 劉軍, 衛(wèi)豐, 茍文選,等. 溫度場下火炮身管固有頻率和振型的三維有限元分析[J]. 火炮發(fā)射與控制學報, 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] 張海航, 狄長安. 某型多管火炮的振動模態(tài)分析[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] 蘇忠亭, 徐達, 楊明華,等. 基于模態(tài)試驗的某火炮身管有限元模型修正[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] 吳東亞, 邢宏光, 崔軍. 基于虛擬樣機技術的某型坦克炮炮身模態(tài)分析[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