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

(1.南京理工大學(xué) 機(jī)械工程學(xué)院, 江蘇 南京 210094; 2.中國(guó)北方工業(yè)公司 軍貿(mào)技術(shù)研究院, 北京 100053)

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

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

(1.南京理工大學(xué) 機(jī)械工程學(xué)院, 江蘇 南京 210094; 2.中國(guó)北方工業(yè)公司 軍貿(mào)技術(shù)研究院, 北京 100053)

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

兵器科學(xué)與技術(shù)； 身管； 固有頻率； 單元支持域徑向基點(diǎn)插值法

0 引言

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

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

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

1 身管振動(dòng)控制微分方程

如圖1所示：考慮火炮身管、炮口制退器、炮尾的綜合模型，炮尾通過(guò)圓柱來(lái)等效，長(zhǎng)度和質(zhì)量與炮尾的質(zhì)量一致；炮口制退器通過(guò)中空?qǐng)A柱形質(zhì)量塊來(lái)等效，長(zhǎng)度和質(zhì)量與炮口制退器相同；與身管端面固結(jié)的截面內(nèi)外徑與身管前端面的內(nèi)外徑一致。圖1中：mb和mm分別表示炮尾和炮口制退器質(zhì)量；D為身管內(nèi)徑；Ls為身管長(zhǎng)度；Lx1和Lx2分別為搖架前后襯套至身管尾端面的距離；Lb和Lm分別為炮尾和炮口制退器長(zhǎng)度；Db和Dm分別為炮尾和炮口制退器的等效尺寸；其余結(jié)構(gòu)尺寸如圖所示。

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 單元支持域徑向基點(diǎn)插值法

通常，采用徑向基點(diǎn)插值法[15]進(jìn)行平衡方程弱形式積分時(shí)，需要根據(jù)節(jié)點(diǎn)的支持域搜索每一個(gè)積分點(diǎn)對(duì)應(yīng)的鄰域節(jié)點(diǎn)，進(jìn)而根據(jù)這些鄰域節(jié)點(diǎn)構(gòu)造插值函數(shù)，這樣不僅對(duì)每個(gè)積分節(jié)點(diǎn)都需要搜索鄰域節(jié)點(diǎn)，而且對(duì)每個(gè)積分點(diǎn)都要計(jì)算一次插值函數(shù)。如圖2所示，點(diǎn)pA和pB的支持域分別為ΩA和ΩB. 不僅如此，由于徑向基函數(shù)是高階連續(xù)函數(shù)，采用數(shù)值積分獲得單元矩陣的過(guò)程中需要采用高階積分方案(即增加積分點(diǎn))，從而大大增加了計(jì)算量，這也是無(wú)網(wǎng)格方法計(jì)算效率低的一個(gè)重要因素。為此，本文提出單元支持域插值方法，單元支持域是指節(jié)點(diǎn)所在的單元及其鄰接單元構(gòu)成的區(qū)域，根據(jù)單元支持域內(nèi)的節(jié)點(diǎn)來(lái)構(gòu)造插值函數(shù)，如圖3所示，e1～e9為單元編號(hào)，點(diǎn)pA和pB的支持域均為單元e5，單元e5的支持域?yàn)棣竤.

利用支持域Ωs上的n個(gè)節(jié)點(diǎn)值和徑向基- 多項(xiàng)式基進(jìn)行插值來(lái)近似函數(shù)u(x)，u(x)的近似函數(shù)以u(píng)h(x)表示，則

(9)

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

將x的支持域上n個(gè)節(jié)點(diǎn)的坐標(biāo)值代入(9)式中，得到：

(10)

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

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

u=Ra+Pb，

(11)

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

(12)

(13)

式中：

(14)

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

Ri(rj)=Rj(ri).

(15)

(9)式中有n+m個(gè)系數(shù)，而(11)式只有n個(gè)方程，為求出n+m個(gè)系數(shù)，本文根據(jù)Golberg等[16]的建議補(bǔ)充如下m個(gè)方程：

(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中的每個(gè)矩陣都是常數(shù)矩陣，因此只需計(jì)算PTR-1一次并保存即可。

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

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

(22)

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

(23)

第k個(gè)節(jié)點(diǎn)的形函數(shù)為

(24)

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

(25)

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

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

(26)

式中:

(27)

由此可得到任一點(diǎn)x處的應(yīng)變?yōu)?/p>

(28)

式中：

(29)

任一點(diǎn)的應(yīng)力為

(30)

式中：彈性矩陣D的表達(dá)式為

(31)

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

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

(32)

(33)

(34)

(35)

(36)

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

(37)

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

(38)

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

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

4 算例

4.1 懸臂梁

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

(39)

(40)

(41)

σy=0，

(42)

(43)

式中：I為梁的截面慣性矩，I=H3/12；x和y分別為梁上的坐標(biāo)；ux和uy分別為梁x和y方向的位移；σx和σy分別為x和y方向的應(yīng)力；τxy為剪應(yīng)力。

為研究本文方法的收斂性，考慮5種不同的網(wǎng)格分布2×2，4×4，8×8，16×16，32×32. 本文方法的計(jì)算結(jié)果以及有限元結(jié)果和解析解的對(duì)比如表1所示。表1中的計(jì)算結(jié)果顯示，本文方法具有較高的收斂率和精度。

4.2 身管固有頻率及數(shù)值驗(yàn)證

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

4.3 身管固有頻率實(shí)驗(yàn)及驗(yàn)證

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

5 結(jié)論

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

)

[1] 陳光宋, 錢林方, 徐亞棟,等. 身管橫向固有振動(dòng)的半解析解法[J]. 兵工學(xué)報(bào), 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)豐, 茍文選,等. 溫度場(chǎng)下火炮身管固有頻率和振型的三維有限元分析[J]. 火炮發(fā)射與控制學(xué)報(bào), 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] 王寶元, 劉朋科, 衡剛,等. 自行火炮固有頻率特性實(shí)驗(yàn)研究[J]. 兵工學(xué)報(bào), 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] 張海航, 狄長(zhǎng)安. 某型多管火炮的振動(dòng)模態(tài)分析[J]. 兵工學(xué)報(bào), 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]. 振動(dòng)與沖擊, 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] 蘇忠亭, 徐達(dá), 楊明華,等. 基于模態(tài)試驗(yàn)的某火炮身管有限元模型修正[J]. 振動(dòng)與沖擊, 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ī)技術(shù)的某型坦克炮炮身模態(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] 龍述堯.無(wú)網(wǎng)格方法及其在固體力學(xué)中的應(yīng)用[M]. 北京：科學(xué)出版社, 2014:1-9.

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

[10] 杜超凡,章定國(guó),洪嘉振.徑向基點(diǎn)插值法在旋轉(zhuǎn)柔性梁動(dòng)力學(xué)中的應(yīng)用[J].力學(xué)學(xué)報(bào),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

國(guó)家自然科學(xué)基金項(xiàng)目(11472137)

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

錢林方(1961—), 男, 教授， 博士生導(dǎo)師。 E-mail: lfqian@vip.163.com

TJ302

A

1000-1093(2017)12-2321-07

10.3969/j.issn.1000-1093.2017.12.004