SPH-FEM耦合方法在彈丸膛內運動中的應用

孫玉杰，崔青春，丁宏民，徐 堅，郭俊行

(西北機電工程研究所，陜西 咸陽 712099)

由于火炮發射時間短，載荷脈沖峰值大，且彈丸在密閉的身管空間內運動，參數不易測量，試驗又難以復現火炮發射環境，發射動力學成為一種可能的分析方法，用物理模型去復現火炮發射過程中各物理參量的變化情況。文獻[1]采用自適應網格技術(Arbitrary Lagrangian-Eulerian，ALE)，解決擠進過程中彈丸網格因大變形而發生的網格畸變問題，通過仿真獲得了彈丸擠進過程中變形及殘余應力變化。文獻[2]對某大口徑榴彈炮發射條件下彈帶擠進過程中的彈帶動態擠進阻力、擠進壓力、彈帶大變形和彈丸運動規律進行了數值模擬研究。文獻[3]建立航空自動炮的熱力耦合有限元模型，采用Fortran子程序結合顯式有限元方法對彈帶擠進及內彈道過程進行了研究，獲得了次要功系數隨時間變化并存在極值，彈帶表層受熱軟化對內彈道過程有顯著影響等結論。文獻[4-5]采用顯式動力學的方法對彈帶擠進過程進行了分析，分析了坡膛結構變化對內彈道性能和坡膛受力狀態的影響。文獻[6-7]采用顯式動力學建立彈帶擠進的熱力耦合模型，將擠進時期的內彈道方程作為力學邊界條件，并采用改進庫倫摩擦模型，得到了擠進過程中擠進阻力、彈丸速度、膛壓和彈丸擺角的變化規律。文獻[8]采用光滑粒子與有限元耦合的方法對某大口徑火炮彈丸擠進過程進行數值模擬，獲得了擠進過程中彈帶材料的塑性流動及應力應變變化。文獻[9]建立了將內彈道計算與彈頭發射過程有限元計算相耦合的模型，并考慮了槍管的彎曲，獲得了彈頭受力及槍口振動情況。文獻[10]建立了預制刻槽彈帶與身管的彈炮耦合火炮發射動力學模型，模型中實現了火藥燃氣壓力隨彈丸運動而動態變化的分布。上述彈炮耦合研究中：有僅對彈丸擠進過程進行了研究；有彈帶預制刻槽，忽略擠進過程；有對小口徑彈丸全膛運動進行了研究。一方面擠進和直膛運動是個連續過程；另一方面大口徑火炮與小口徑火炮在身管內膛和彈丸結構等方面差異較大，上述研究不能直接應用于大口徑火炮全膛運動分析。

大口徑火炮彈丸全膛運動過程中，彈帶經歷侵徹、擠壓，發生大的塑性變形和材料失效。對該問題的模擬，傳統的基于網格的有限元存在單元畸變，導致計算精度和效率下降，甚至無法完成有效的計算。針對傳統FEM(Finite Element Method)的局限性，SPH(Smoothed Particle Hydrodynamics)等無網格技術應運而生，SPH法以一種拉格朗日形式的無網格粒子代替網格單元，粒子攜帶著物質的材料特性和力學量信息，具有良好的自適應性，克服了FEM由于產生大變形畸變而終止計算這一難題，在處理侵徹、穿甲等大變形問題上極具優勢。但SPH算法相比傳統有限元方法，存在計算效率低，且不易施加邊界條件。大口徑火炮彈丸全膛運動過程中，只有彈帶發生大的塑性變形和材料失效。而整個發射系統其它部分發生剛體運動和小的彈性變形。故對大變形的彈帶采用SPH粒子模擬，發射系統的其它結構采用FEM建模，通過FEM-SPH耦合算法實現結構相互間的力學行為傳遞，該方法結合了SPH方法處理大變形能力和FEM的計算效率高的互補優勢。

本文采用有限元和光滑粒子耦合方法，研究某大口徑火炮彈丸全膛運動過程，獲得彈帶溫度變化、影響炮口振動因素和彈丸運動姿態的變化，并與試驗數據進行對比，驗證SPH-FEM模型的有效性。該研究為大口徑火炮彈炮匹配設計及精度設計提供了參考。

1 大口徑火炮非線性有限元動力學模型

以某大口徑履帶式自行火炮為研究對象，采用SPH-FEM耦合方法建立了彈炮耦合的大口徑火炮非線性有限元動力學模型，分析其在水平地面上靜止狀態下射擊時動力響應。

1.1 各部件拓撲關系

動力學模型中包含彈帶、彈丸、身管、炮口制退器、炮尾、搖架、前后銅套、定向栓、炮塔、座圈、底盤、懸掛和履帶等。模型中各部件拓撲關系如圖1所示。hi(i=1～38)為部件間的連接關系。其中：彈丸與彈帶為綁定連接(h1)；彈丸定心部、彈帶與身管內膛間為接觸作用(h2，h3)；炮口制退器、抽氣裝置、炮尾閂體與身管間為綁定連接(h4，h5，h7)；身管與前后銅套間為接觸作用(h6，h8)；復進機和制退機的后坐部分與炮尾閂體為綁定連接(h9，h10)；炮尾閂體與定向栓間為綁定連接(h11)；前、后銅套與搖架間為綁定連接(h12，h16)；定向栓與搖架定向栓室間為接觸作用(h14)；復進機和制退機的非后坐部分與搖架為綁定連接(h13，h15)；高低機主齒輪與搖架齒弧間為接觸作用(h17)；搖架和炮塔間為耦合約束關系(h18)，即只釋放搖架繞耳軸軸線方向的旋轉自由度；平衡機、高低機、上座圈、方向機和炮塔，下座圈和底盤為綁定連接(h20～h22，h26)；上座圈和滾珠、下座圈和滾珠、方向機主齒輪與下座圈內齒圈為接觸作用(h23～h25)。扭力軸一端與車體固結，一端與平衡肘通過花鍵連接(h27，h28)；平衡肘與負重輪、誘導輪與底盤、拖帶輪與底盤、主動輪與底盤為耦合約束關系(h29～h32)，即只釋放其繞轉動中心方向的旋轉自由度；負重輪、誘導輪、拖帶輪和主動輪與履帶內表面為接觸作用(h33～h36)；履帶外表面與地面為接觸作用(h37)。

圖1 某大口徑履帶式自行火炮拓撲關系圖Fig.1 Topology relationships of large-caliber artillery

1.2 網格模型

某大口徑火炮有限元模型如圖2所示。模型中含有粒子單元、實體單元、殼單元、質量單元和連接單元等。總的單元數為1 950 346，總的節點數為2 338 318。坐標原點O位于座圈上平面圓心處，x軸沿水平方向指向炮口，y軸垂直方向向上，z軸根據右手法則確定。

身管膛線為空間曲線，且膛線斷面尺寸與身管長度尺寸存在量級差異，利用專業前處理軟件的強大網格劃分功能，除膛線起始部有楔形單元外，身管其余采用六面體單元，為準確刻畫擠進段彈帶與身管內膛的相互作用，對膛線起始部進行了網格加密，整個身管包含543 264個單元。身管局部有限元模型如圖3所示。

圖2 大口徑履帶式自行火炮有限元模型Fig.2 Finite element model of large-caliber artillery

彈帶采用光滑粒子進行離散。有兩種方法進行建模，一種為自適應算法，先建立彈帶的有限元網格，并定義單元轉變粒子準則，當某個單元滿足轉換準則，該單元轉換為粒子；第二種是分析一開始就建立彈帶的SPH粒子。由于Abaqus軟件對第一種方式的并行有限制要求[11]，為提高計算效率，本文采用第二種方式。所建立的彈帶SPH模型如圖4所示，彈帶粒子單元數為746 916。

圖3 身管局部網格模型Fig.3 Partial mesh model of gun barrel

圖4 彈帶SPH模型Fig.4 SPH model of rotating bands

1.3 材料模型

Johnson-Cook材料模型作為一種理想剛塑性強化模型既可以反映材料在高應變速率下應變硬化、應變速率硬化，又可以反映熱軟化效應[12]。Johnson-Cook材料模型中屈服應力是塑性應變、應變率以及溫度的函數

(1)

T*=(T-Tr)/(Tm-Tr)

(2)

(3)

為了獲得更加精確的材料參數，開展彈帶材料在較寬應變速率和溫度范圍內力學性能測試，采用Gleeble3500熱模擬試驗機進行準靜態及低應變率下彈帶材料壓縮試驗；采用帶溫度控制裝置的分離式霍普金森壓桿裝置進行彈帶材料高應變率下壓縮試驗。通過對試驗數據處理，回歸得到Johnson-Cook方程的各個材料參數為：A=133 MPa，B=324 MPa，C=0.043，m=1.21和n=0.48。圖5為黃銅材料在3 000 s-1應變速率下應力-應變關系試驗測量值與Johnson-Cook模型預測值對比，從圖中可以看出，預測結果均與試驗值吻合較好。

圖5 3 000 s-1應變率下應力-應變關系試驗值與預測值對比Fig.5 Comparison of stress-strain curves between experimental measurements and predictions under 3 000 s-1 strain rate

(4)

身管為炮鋼材料，前后銅套為鑄造鋁青銅，其它部件材料為合金鋼，靶場地面為混凝土，對應的材料參數如表1所示。

表1 動力學模型中材料參數Tab.1 Material parameters of dynamical model

1.4 彈帶與身管內膛摩擦界面特性

彈帶與身管內膛的摩擦界面特性對于彈丸啟動壓力、彈丸膛內運動，身管內膛損傷有重要影響。借鑒在高速切削以及摩擦攪拌焊接研究中，切向摩擦行為多采用修正的庫倫摩擦模型來表征接觸面之間高速重載的摩擦特性。切接觸面的切向應力由式(5)決定

(5)

式中：f為切向應力，MPa；μ為靜摩擦因數，取為0.13；σn為正壓力，MPa；τs為彈帶材料的隨溫度變化的剪切極限，MPa，取為

(6)

彈丸彈帶和彈丸定心部與膛線的接觸和碰撞狀態，通過平衡主控—從屬搜索算法進行計算。

1.5 載荷、邊界條件處理

火炮發射時，作用于發射系統的載荷有膛底壓力、彈底壓力、波爾登力、制退機力、復進機力和平衡機力等。動力學模型中采用電子測壓彈獲得實彈射擊時的膛底壓力，并由文獻[14]第“3～5”節所述內彈道計算中應用的壓力換算關系，換算得到彈底壓力。由于溫度分布不均勻和重力的作用，身管在發射前向下彎曲。火炮發射時，膛內高壓火藥氣體的作用必然使彎曲的身管產生波爾登效應，即彎曲的身管在發射時受到波爾登力的作用。隨著現代火炮身管長度的增加，波爾登載荷對彈丸膛內運動和炮口振動的影響不可忽略。通過子程序VDLOAD在彈后空間施加拉格朗日分布的壓力來模擬波爾登力的影響。復進機力通過在復進機與炮尾和搖架的連接點上施加一對共線且反向的隨時間變化的集中力來模擬[15]。制退機力和平衡機力也采用相同的方法。其隨時間變化通過相應的計算書獲得。火炮靜止狀態下射擊，主動輪制動，約束其轉動自由度。

2 火炮發射動力學模型仿真分析

彎曲身管將改變彈丸在膛內的運動狀態，為準確獲得火炮發射過程中的動態響應，采用隱式、顯式聯合求解的方法，充分利用兩種求解器的優勢。先采用隱式求解器獲得火炮發射動力學模型在重力載荷作用下的靜變形及等效應力，起落部分垂直方向最大變形量發生在炮口制退器處，大小約為10.6 mm，最大等效應力位置位于平衡機處，大小約為320 MPa，如圖6所示。然后將此狀態作為火炮發射動力學模型的初始構型，重新設定約束狀態、載荷等，利用顯式求解器對發射過程的動態響應進行求解。分別在引信質心和彈丸質心處設置兩個參考點，分別與少數單元設置剛體約束。定義兩參考點的連線與彈丸質心處的速度矢量之間的夾角為彈丸章動角，其單位為 °。

圖6 動力學模型在重力載荷作用下的響應Fig.6 Response of dynamics model under the gravity load

2.1 火炮發射動力學模型驗證

通過高速攝影的標記點與動力學模型中的對應位置的位移進行比較，來驗證所建立動力學模型的準確性。從圖7(a)可以得出：計算值與測量值變化趨勢相同，膛內時期處于加速后坐過程，后坐速度逐漸增大，故后坐位移曲線外凸。出炮口時刻后坐位移約為115 mm。從圖7(b)可以得出：垂直于身管方向位移的計算值與測量值變化趨勢相同，4 ms前，由于彈丸處于擠進過程，此時身管扭轉剛度很大，加上彈丸轉速低，彈炮作用不強烈，垂直于身管方向位移幾乎為零。隨著彈丸向前運動，身管扭轉剛度減小，加上彈炮相互作用強烈，垂直于身管方向位移逐漸增大，一直到彈帶出炮口前，達到最大值。當彈帶出炮口后，由于彈帶對身管扭轉作用消失，身管有一定的回復。

圖7 動力學模型計算結果與高速攝影測量值對比Fig.7 Comparisons between simulation results and high-speed video measurements

2.2 火炮發射動力學模型結果分析與討論

2.2.1 計算效率分析

顯式時間積分的穩定步長近似估計為

(7)

式中：Δtmin為穩定步長，s；lmin為單元最小邊長，mm；E為彈性模量，MPa；ρ為密度，kg/m3；ν為泊松比。

彈帶采用ALE時，由于彈帶在擠進過程中和隨后的直膛運動過程中，發生大變形，造成單元畸變，特征長度減小，造成穩定步長變小，從而大大增加求解時間。而彈帶采用SPH時，穩定步長隨彈帶發生變形而減小的程度較小，從而保證了較高的計算效率。圖8為顯式動力學模型中穩定步長隨時間變化。從圖中可以看出，采用SPH時穩定步長在4.4×10-8～4.64×10-8s。

圖3列出了連續神經網絡算法在空手道俱樂部網絡上的仿真結果。從實驗結果發現，特征向量算法、連續算法對空手道俱樂部網絡的分類結果與特征向量算法的分類結果（Q=0.3715）相同，表明本文中的算法也是合理的。

圖8 穩定步長隨時間變化曲線Fig.8 Curve of stable time increment with time

2.2.2 膛內運動彈帶溫度

獲得彈丸在膛內不同時期彈帶溫度變化，如圖9所示。擠進過程中，彈帶發生大的塑性變形，塑性功占主導，塑性功轉變為內能，彈帶溫度升高。擠進完成后，隨著彈丸速度增加，摩擦阻力做功占主導，摩擦阻力做功轉變為內能，使得彈帶與膛線緊密貼合部位溫度繼續升高。擠進初期彈帶溫度云圖如圖9(a)所示，最高溫度為102 ℃。前彈帶擠進完成時，彈帶溫度云圖如圖9(b)所示，最高溫度為457 ℃。擠進完成時，彈帶溫度云圖如圖9(c)所示，最高溫度約為585 ℃，這與文獻[16]研究相符，擠進完成，彈帶表面溫度超過再結晶溫度280 ℃，但未達到熔點。最大膛壓點時，彈帶溫度云圖如圖9(d)所示，最高溫度為636 ℃。彈丸行程一半時，彈帶溫度云圖如圖9(e)所示，最高溫度為846 ℃。彈丸進入半約束期時，彈帶溫度云圖如圖9(f)所示，最高溫度為1 033 ℃。由圖9還可以得出：彈丸膛內運動，彈帶表層溫度逐漸升高，彈帶表層甚至達到材料的熔點，彈帶與身管間建立了熔融潤滑，滑動機理變為流體動力潤滑滑動[17]，由于作用時間短，熱量來不及向彈帶內部傳導，彈帶內部仍然為初始溫度，彈帶形成軟硬結合的結構，造成摩擦因數減小。

圖9 彈丸運動不同位置彈帶溫度分布云圖Fig.9 Contour of rotating band temperature distribution at different locations

2.2.3 影響炮口運動因素分析

火炮發射過程中，受力復雜，各種因素對炮口運動相互影響，且同時發生，采用試驗測試很難區分其單獨的影響，發射動力學成為一種可能的方法去區分各個因素對炮口點運動的影響。采用所建立的模型分別計算后坐運動、后坐運動中計及波爾登力、后坐運動中計及彈炮耦合和后坐運動中計及波爾登力和彈炮耦合四種工況。為描述炮口振動，定義炮口中心點與炮口制退器前端面耦合，輸出其在整個發射過程中的位移變化，耦合示意圖如圖10所示。各工況下炮口點的運動如圖11所示，圖11(a)為炮口點垂直方向位移，圖11(b)為炮口點水平方向位移。

圖10 炮口點耦合示意圖Fig.10 The illustration of the coupling between reference point and muzzle

垂直方向：該火炮后坐部分質心在炮膛軸線下方，膛內期，由于動力偶臂的作用，起落部分發生俯仰運動，炮口點位移向上；由于重力作用，身管初始構型向下彎曲，則波爾登力有使身管變直的趨勢，后坐運動也使垂直方向位移向上，兩者綜合作用使得炮口點位移向上。曲線的振蕩是由于膛壓移動導致；9 ms之前由于彈丸與身管的作用不強烈，和后坐運動計算結果基本相同。此后彈丸定心部與身管內膛的接觸碰撞加劇，影響炮口點的振動。11.6 ms前，計及波爾登力和彈炮耦合與考慮波爾登力時兩者變化一致，此后由于彈丸前定心部與身管內膛的碰撞加劇，影響炮口位移。波爾登力和彈炮作用對垂直方向位移影響顯著。

圖11 不同因素對炮口點位移影響Fig.11 The effect of different factors on muzzle displacement

垂直方向：由于身管初始構型在水平方向變形量很小，波爾登力作用不如垂直方向顯著，引起的位移很小；4 ms之前，位移為零，此后由于彈丸和身管的相互作用引起位移在零點附近跳動，隨著彈丸前定心部與身管內膛的碰撞加劇，水平方向位移變化增大。水平方向位移主要受彈炮作用影響，且水平方向位移小于垂直方向位移。

2.2.4 彈丸不平衡因素的影響

由第“2.2.3”節分析可知，彈炮作用對于炮口點垂直方向和水平方向位移都有顯著影響。膛內期，彈丸約逆時針旋轉2.2圈(從炮口看向炮尾)，出炮口時彈丸質量偏心方位與卡膛時質量偏心方位恒相差約72°。當卡膛時質量偏心方位位于象限點1上，出炮口時質量偏心方位位于左上方1*處，箭頭所示方向為彈丸旋轉方向，如圖12所示。炮口點在Y-Z面的位移軌跡分布在第一、二象限內，膛內運動前期，彈丸轉速較低，彈炮作用不強烈，曲線重合，隨著轉速增大，質量偏心大的，慣性離心力也大，曲線分離。圖13為膛內彈丸章動角隨時間變化，膛內運動前期，章動角幅值較小，但變化較快，隨后的直膛內運動，章動角幅值變大，出炮口時達到最大值，質量偏心影響明顯，零偏心時0.068 37 °增大至偏心0.2 mm時0.138 33 °，章動角幾乎增大兩倍。

彈丸裝填時，質量偏心方位在整個圓周范圍內隨機分布。分析彈丸質量偏心為0.2 mm，質量偏心分別位于四個象限點上時炮口點振動和章動角的變化。從圖15可以看出：彈丸質量偏心方位相位差為90°，則彈炮作用的方位也相差90°，出炮口時位移也基本上相差90°，質量偏心方位對炮口點的振動影響顯著，如圖14所示。當質量偏心方位位于2位置時，彈丸膛內時期章動角幅值最大，質量偏心方位位于4位置時，彈丸膛內時期章動角幅值最小，如圖15所示。

2.2.5 初始裝填姿態的影響

上述模型中，彈丸幾何軸線與身管軸線平行。由于裝填誤差等原因，裝填到位后，彈丸軸線與身管軸線不完全重合，也就是兩者之間有一定的夾角。分析質量偏心為0.1 mm，卡膛時質量偏心位于象限點1上，卡膛角取為仰4′、仰2′、0′、俯2′和俯4′(仰為鉛垂面內，彈軸在身管軸線上方，并與之相交；反之定義為俯)。由于鉛垂面內，身管向下彎曲，身管內膛約束彈丸沿其軸線運動，當卡膛角為仰時，炮口點的振動變化相同，當卡膛角為俯時，也有相似變化，如圖16所示。隨著裝填角度的增大，彈丸的章動角波動范圍變大，出炮口時，仰4′相對于0′，章動角增大2.08倍。當卡膛角為俯時，出炮口時，仰4′相對于0′，章動角增大1.7倍。相同卡膛角時，卡膛角為俯時，無論章動角的最大值還是出炮口時相應值都大于卡膛角為仰時，如圖17所示。

圖12 相同偏心方位不同質量偏心對炮口點振動的影響Fig.12 The effect of mass eccentricity on the vibration of muzzle under the same orientation of mass eccentricity

圖13 相同偏心方位不同質量偏心對章動角的影響Fig.13 The effect of mass eccentricity on the nutation angle of projectile under the same orientation of mass eccentricity

圖14 相同質量偏心不同質量偏心方位對炮口點振動的影響Fig.14 The effect of the orientation of mass eccentricity on the vibration of muzzle under the same mass eccentricity

圖15 相同質量偏心不同質量偏心方位對章動角的影響Fig.15 The effect of the orientation of mass eccentricity on the nutation angle of projectile under the same mass eccentricity

圖16 不同初始裝填姿態對炮口點振動的影響Fig.16 The effect of the initial loading posture on the vibration of muzzle

圖17 不同初始裝填姿態對章動角的影響Fig.18 The effect of the initial loading posture on the nutation angle of projectile

3 結 論

(1)針對大口徑火炮彈丸全膛運動過程中因為彈帶高速侵徹和擠壓導致網格畸變而使計算效率低的問題，彈帶采用SPH模擬，火炮其它結構采用FEM模擬，通過FEM-SPH耦合的方法實現結構相互間的力學行為傳遞。即利用了SPH處理大變形優勢，又充分利用了FEM計算精度高的特性。

(2)彈丸膛內運動過程中，由于彈帶發生塑性變形和摩擦力做功，擠進完成后彈帶表面超過再結晶溫度；隨后的膛內運動中，彈帶表面溫度接近熔點。

(3)波爾登力和彈炮作用對炮口點垂直方向位移影響顯著；水平方向位移主要受彈炮作用影響，且水平方向位移小于垂直方向位移。

(4)質量偏心和質量偏心方位對于炮口點運動影響顯著，裝填姿態對炮口點振動有一定的影響；裝填姿態對章動角影響最大，質量偏心影響次之，質量偏心方位影響最小。