BAT子彈藥折疊翼展開沖擊非光滑模型的隱式積分算法

張洪銘，顧曉輝，孫 麗

(1.南京理工大學 機械工程學院，南京 210094；2.江蘇科技大學 機械工程學院，江蘇 鎮江 212100)

折疊翼結構能有效縮小中近程導彈、制導炮彈的儲運空間，利于增強車輛、戰機的運載能力，提高武器系統的戰術性能，在現代飛行器的設計中獲得了廣泛的應用[1]。不同于固定翼結構近似剛性的連接，折疊翼展開機構中普遍包含鉸鏈、限位擋塊、鎖定插銷及驅動裝置?；顒蛹拇嬖诮档土私M件間的連接強度和系統的結構剛度，在迅速展開過程中，驅動加載的瞬間以及翼面展開到位時刻的碰撞、鎖定動作，均會在部件間產生沖擊載荷，存在使構件被破壞的可能性。另外，翼面的展弦比較大時，展開末段的碰撞沖擊會引起翼面振顫，不恰當的設計可能引發翼面與彈體共振，加劇結構破壞的程度或影響飛行器的運動軌跡。

因此，針對折疊翼展開的動態過程，沖擊載荷特性的分析與優化是提升展開性能的關鍵?，F有研究對折疊翼展開過程進行了廣泛研究。林三春等[2]結合隱式動力學仿真分析和試驗研究的結果，指出柵格翼展開沖擊的主要影響因素是阻尼力矩與展開到位時的角速度，氣動力和過載為次要因素；唐霄漢等[3]詳細討論了折疊翼機構中常見的鉸鏈關節在沖擊過程中的內力分析，討論了材料采用動態本構模型時，應變率對等效塑性應變的影響；甄文強等在文獻[4]中建立了考慮氣動和摩擦阻力的數學模型，并與動力學仿真結果進行了對比分析；更進一步地，其在文獻[5]中進行了發射環境下折疊翼展開過程的研究，為折疊翼展開機構的工程設計提供了參考。上述研究對折疊翼展開過程中的沖擊載荷、動力學響應進行了分析，關于沖擊載荷的優化分析，李莉[6]、蔡德詠[7]等分別建立了縱向和橫向折疊翼的剛體模型，并均以展開到位瞬時的翼片動能為優化目標，對展開系統進行了優化。趙俊鋒等[8]以剛柔耦合模型為基礎，采用序列二次規劃法(sequential quadratic programming,SQP)對驅動力-時間曲線進行優化以降低彈翼展開沖擊；韓雪峰等[9]則采用正交試驗對折疊翼展開機構進行優化。郭宏偉等[10]在提出沖擊加速度、力矩計算方法的同時，給出了兩種降低展開沖擊的設計方案。

現有研究較多關注折疊翼展開到位時，翼面對沖擊載荷的響應模擬，對于接觸/碰撞沖擊進行建模、計算的內容相對較少。實際上，準確計算折疊翼展開過程中的沖擊力，是對展開機構進行優化設計、可靠性分析的基礎和前提。折疊翼展開過程中包含的接觸/碰撞現象，屬于典型的非光滑、非線性過程，其動力學建模和求解相對困難和復雜。目前，描述接觸/碰撞過程的動力學模型主要有兩類[11-12]：基于接觸體幾何約束的非光滑模型[13-14]和基于接觸界面變形量的連續接觸力模型[15-17](或稱為罰函數模型[18])。在包含非光滑過程的動力學模型中，由不等式來描述系統中的單邊約束，對不等式的處理是建立模型的難點。針對該難點，非光滑模型通過微分包含理論，對接觸/碰撞、摩擦等單邊約束進行表達和處理；連續接觸力模型是另一種被廣泛研究的接觸/碰撞動力學模型，采用接觸力約束代替幾何約束，將接觸力模型引入到動力學方程中進行計算。相較于連續接觸力模型，非光滑模型采用局部剛性接觸避免了計算剛性微分方程的困難，計算精度較高[19]，可以處理碰撞過程中能量不協調的問題[20]。另外，基于線性互補問題(linear complementarity problem，LCP)的非光滑動力學模型，不僅使得力學問題具有更加自然、精確的數學描述，而且可以運用互補問題解的存在性和唯一性理論分析力學問題的結構模型，構建更加有效和強健的算法。同時，互補形式的非光滑動力學模型能夠自動判斷觸-離、滯-滑等狀態的切換，避免了進行事件檢測時所需的大量計算資源[21]。

基于上述研究結果，以某型智能反坦克(brainpower anti-tank,BAT)子彈藥的折疊翼機構為研究背景，建立該縱向折疊翼展開過程的非光滑動力學模型，并推導模型的線性互補形式；文中提出了一種新型隱式積分算法用于求解線性互補形式的非光滑動力學模型，并使用基于所提隱式積分算法的時間步進算法，來計算翼面展開到位時碰撞沖擊力量值的大小。將模型計算結果與數值仿真結果進行了對比，驗證文中非光滑模型的計算精度與求解效率。結果表明，線性互補形式的非光滑動力學模型能夠有效描述接觸/碰撞過程中的單邊約束，并且文中提出的隱式積分算法相較于Moreau中點法具有更高的計算精度。

1 折疊翼展開過程的非光滑動力學模型

1.1 折疊翼構型及展開特性

某型BAT子彈藥外形如圖1所示，由運載平臺釋放后展開折疊翼和尾翼，配合紅外/聲傳感器復合導引，滑翔至目標上空實施攻擊。彈體中部裝有四片與彈體縱軸垂直、展開后呈十字型分布的折疊式彈翼。

圖1 某型BAT子彈藥外形圖Fig.1 Shape diagram of BAT submunitions

本文的研究內容主要關注縱向彈載折疊翼展開過程中，碰撞沖擊的建模與計算，因此，忽略聲傳感器支撐桿部分。單組折疊翼組件的裝配示意如圖2所示，展開系統由彈體基座、限位擋塊、轉軸、翼面組件、連桿、滑塊以及驅動彈簧共同組成。圖2中，Ox為彈軸方向，Oy為彈徑方向，A、B、C為三處鉸接點。BAT子彈藥折疊翼展開動作過程為：被運載平臺拋散釋放后，在展開過程的初始時刻，翼面處于收攏狀態，圖2中φ1為0°；隨后滑塊在壓縮彈簧的驅動下沿滑軌向彈頭方向運動，通過連桿對翼面組件施加轉動力矩，使得翼面繞鉸接點A逆時針轉動；翼面展開到90°位時與限位擋塊之間發生碰撞，同時鎖定機構啟動，對展開到位的翼面進行鎖定防止回轉，展開動作結束。

圖2 單組折疊翼展開機構Fig.2 Schematic diagram of folding wing deployment mechanism

依據設計指標，折疊翼的展開動作需要在0.2 s內完成，即翼面組件需要在短時間內完成大范圍轉動，轉動速度高、動能大，在展開動作末段翼面與限位擋塊、鎖定銷之間會發生相對劇烈的碰撞沖擊。沖擊過程持續時間短、作用力大，使得包含碰撞沖擊的折疊翼的展開動作是一個典型的非光滑動力學過程。準確計算碰撞力量值對展開機構優化、機構強度校核均具有重要意義。后續內容建立能夠準確描述碰撞過程的非光滑動力學模型，和用于解算非光滑模型數值計算方法，用于計算翼面與限位擋塊首次碰撞力的大小。

1.2 折疊翼非光滑動力學建模

對圖2所示折疊翼多體系統，選取最大坐標數目，在絕對笛卡爾坐標下描述各部件的位形。在未發生接觸/碰撞現象的自由運動階段，系統內僅有光滑的等式雙邊約束，基于拉格朗日乘子法將約束方程與運動微分方程結合，對含有n個自由度、m個雙邊約束(m

(1)

(2)

(3)

通常取α=β∈[5,50]系統響應最快，在每個時間步內，系統的位置約束與速度約束同時被滿足，使得積分誤差得到部分抑制。式(3)即為光滑多體系統的完備動力學方程。

多體系統內部構件間發生接觸/碰撞時，系統拓撲構型發生改變，除原有雙邊約束外，系統內增加了額外的單邊約束使構件運動狀態發生突變。對于多體系統內所有潛在發生正、斜碰撞的接觸面，借助Lagrange乘子法向多體系統的運動微分方程中添加法向、切向約束反力得到

(4)

式中：λN、λT分別為未知的法向碰撞力和切向摩擦力；WN與WT分別為法向、切向距離約束函數對廣義坐標q的Jacobian矩陣的轉置

(5)

(6)

式中,gN(q,t)、gT(q,t)分別為接觸面法向、切向距離函數列陣。

由于碰撞過程發生時間短、作用力大，碰撞前后速度發生突變，考慮將式(6)中力-加速度關系在碰撞過程時間Δt上積分，寫成速度-沖量形式。在碰撞起始時刻ta至結束時刻tb上，式(4)兩側同時關于時間積分得到

(7)

(8)

式中，ΛN和ΛT分別為法向、切向上的單邊約束力沖量。將系統的幾何約束方程對時間微分，代入tb時刻的廣義速度，得到目標構件在碰撞結束時刻的速度約束方程

(9)

將碰撞結束時刻接觸面間法向、切向上的相對速度用廣義速度表示為

(10)

(11)

非光滑動力學角度下的互補理論是對接觸/碰撞等非線性過程中物理量所遵從規律的描述?；诨パa理論，對系統內滿足gN(q,t)=0條件的接觸對，法向碰撞力沖量ΛN和切向摩擦余量的沖量ΛTS，分別與各自對應方向上，碰撞前后相對速度的合成量構成互補關系

0≤ΛN⊥ξN≥0

(12)

0≤ΛTS⊥ξT≥0

(13)

式中：符號⊥表示兩個向量垂直；切向摩擦余量的沖量ΛTS=μΛN-ΛT；ξN、ξT分別為法向、切向合成速度

(14)

(15)

式中,eN、eT分別為法向、切向上的恢復系數對角矩陣。

由式(7)、式(9)、式(12)～式(15)共6個方程共同構成了基于互補接觸定律的非光滑動力學模型

(16)

2 非光滑模型的數值計算方法

形如式(16)的微分代數方程組(differential algebraic equations,DAEs)從速度-沖量角度描述了多剛體系統內部構件之間發生接觸/碰撞現象時，多體系統的運行規律。式(16)闡明了單邊約束力沖量與合成速度之間存在互補關系，但未顯式給出兩者間的數量關系，無法直接求解。本章推導式(16)的LCP形式以求解非光滑過程，求解主要涉及兩部分內容，即成對互補量的計算和DAEs的迭代求解。通過確定相鄰時間點上廣義速度和廣義坐標的隱式積分算法，從而對形如式(16)的DAEs采用時間步進算法進行求解。

2.1 非光滑模型的線性互補形式

采用庫倫模型描述接觸面件的動/靜摩擦力，包含庫倫摩擦模型的非光滑動力學模型是一個微分包含[22]，在切向相對速度為零時，摩擦力是一個多值函數。為了便于判斷接觸狀態的轉換，對于法向和切向上相對距離、相對速度為零的接觸對，將摩擦余量的沖量分解為左、右分量

ΛL=μΛN-ΛT

(17)

ΛR=μΛN+ΛT

(18)

則

ΛR=2μΛN-ΛL

(19)

相應地將切向合成速度分解為左、右分量之差

ξT=ξR-ξL

(20)

將式(17)代入式(16)第一式、式(20)代入式(16)第四式，聯立式(14)、式(19)后寫成矩陣形式有

(21)

將式(21)等號兩側內容對換并按分塊矩陣簡寫為

(22)

(23)

將式(23)代入式(22)第二個方程，得到Λ的表達式

(24)

將式(24)代入式(22)第三個方程中，得到互補量y和x之間的標準線性互補形式

(25)

式中：

2.2 非光滑模型的隱式積分算法

(26)

基于上述討論，構造非光滑動力學模型在速度-沖量層面的時間步進算法，求解描述折疊翼展開系統中含有接觸/碰撞過程的微分代數方程組。時間步進算法的流程如圖3所示，已知展開系統在t時刻的狀態，判斷該時刻是否有接觸對在切向、法向上的距離為零：若非零，則采用Newmark隱式迭代算法求解系統平穩運行至下一時刻的坐標和速度；若單邊約束為零，表明系統內部構件間發生接觸/碰撞，動力學方程組由式(5)變為式(18)，解出一對互補的單邊約束力沖量與合成速度后，按式(23)中的迭代算法計算非光滑過程下一時刻的廣義坐標和速度，而后重新判斷單邊距離約束是否為零，從而進入迭代，直至總計算時間的終點時刻tF。

3 BAT子彈藥折疊翼展開沖擊分析

前兩章內容建立了縱向折疊翼展開過程的非光滑動力學模型，并構造了非光滑模型的線性互補形式，該形式使用文中提出的隱式積分算法進行求解時，相較于Moreau中點法能夠獲得更加精確的碰撞力值。本章通過一個工程算例，驗證基于隱式積分算法的時間步進算法對非光滑動力學模型進行數值計算的精度和計算效率，并對計算結果進行分析。

某型BAT子彈藥折疊翼展開動作原理如圖4所示。翼面解鎖后展開的主要動力來源于滑塊底部的彈簧組件，處于壓縮狀態的彈簧推動滑塊沿滑桿向前運動，滑塊帶動連桿推動翼面繞鉸點A轉動，翼面與限位擋塊碰撞后展開動作完成。在折疊翼從收攏至展開的過程中，驅動彈簧始終處于壓縮狀態。圖4系統中各量的取值如表1所示。表1中：K為驅動彈簧的剛度系數；L0和L分別為驅動彈簧的初始長度和壓縮剩余長度；H為0.075 m。圖4展開機構中滑移鉸和柱鉸均處于良好潤滑條件，摩擦力量值相較于碰撞力很小，故而不考慮系統中摩擦力的影響。翼面以及擋塊使用的材料為7A04-T6鋁合金，硬度HB150。

圖4 展開機構原理Fig.4 Schematic diagram of folding wing expansion principle

表1 BAT子彈藥折疊翼機構參數Tab.1 Parameters of folding wing mechanism of BAT submunitions

圖4中彈簧的驅動力-時間關系如圖5所示。力的量值為負表示彈簧處于壓縮狀態，為正表示拉伸狀態。

圖5 彈簧驅動力-時間曲線Fig.5 Spring force-time curve

對圖4中的折疊翼展開系統，建立與彈體固連的慣性笛卡爾坐標系Oxy，在圖示二維平面內選取最大坐標數目，令q=[x1,y1,θ,x2,y2,φ,x3,y3,γ]T為廣義坐標，建立折疊翼展開過程的非光滑動力學模型。折疊翼展開系統的雙邊約束函數為

由于該系統為定常系統，雙邊約束函數中不顯含時間，則可將系統的雙邊約束方程寫成Φ(q)=0；在翼面展開到位并與限位擋塊之間發生碰撞階段的單邊約束函數為

將對應項代入式(21)～式(22)中，得到折疊翼展開過程的非光滑動力學模型。使用文中提出的隱式積分算法，取碰撞過程的時間步Δt=0.005 s，解算從收攏至展開到位的完整過程中，翼面自由端的角位移和角速度在不同時刻的值，結果如圖6所示。

圖6 翼面角位移-時間曲線和角速度-時間曲線Fig.6 Angular displacement-time curve and angular velocity-time curve of folded wing

為了驗證2.2節中基于隱式積分算法的數值算法的計算精度與穩定性，分別考察文中所提算法的計算結果與連續力模型計算結果的吻合程度，以及計算結果對幾何約束的相對誤差。采用與數值計算相同的參數，使用連續接觸力模型對折疊翼展開過程進行計算(連續力模型的計算在RecurDyn?軟件中實現)。將不同展開時間下，兩種模型計算得到的翼面角位移、角速度分別進行對比，結果如圖7所示。

圖7中曲線對比結果表明：數值計算結果與仿真結果吻合；首次碰撞發生后，基于隱式積分算法的數值計算方法所得出的最大反彈角位移與反彈角速度，略微大于連續接觸力模型的結果。原因在于：連續力模型假設物體界面在相互接觸/碰撞的過程中可以發生微小的穿透，且碰撞力由接觸力和阻尼力兩部分組成，兩者均為界面相互侵入量的函數。連續力模型中的阻尼力在接觸/碰撞過程中起耗能作用，因此使用罰函數模型算得的仿真結果略小于數值計算的結果。

(a) 角位移對比

圖4所示系統在展開過程中，翼面的質心與鉸點A之間的距離理論上應該為固定值，但由于數值計算存在積分誤差，使得該固定值的計算結果在理論值附近有微小偏移。采用該固定值數值計算結果與理論值之間的相對誤差，來度量文中所提算法的精度具有較高的合理性?；陔[式積分算法的數值計算結果與理論值之間的相對誤差如圖8所示。相對誤差值的數量級在10-8量級，且整個展開過程中誤差均值穩定并且分布均勻，表明文中所提的基于隱式積分算法求解的非光滑動力學模型能夠準確、穩定地描述含有碰撞過程的折疊翼展開動作。

圖8 文中所提方法求解結果的相對誤差Fig.8 The relative error between the solution result of the proposed method and the exact solution

基于隱式積分算法的時間步進算法能夠更準確地計算接觸力。使用文中提出方法、Moreau中點算法、連續力模型算法，分別對圖4中的折疊翼展開過程進行計算，依據力與沖量之間的關系，將文中提出的隱式積分法和Moreau中點算法的碰撞沖量結果換算成碰撞力，三種算法在所有時間步上的碰撞力隨展開時間的變化如圖9中曲線所示。觀察圖9中的放大圖，比較三種方法算得的首次碰撞力的大小，文中提出的隱式積分算法結果與RecurDyn計算結果吻合，Moreau中點算法的計算結果偏大；將文中提出的隱式積分法計算結果與連續接觸力模型的計算結果(由RecurDyn實現)相比較，結果表明：隱式積分法基于非光滑動力學模型，碰撞過程中接觸對表面之間不存在相互侵入過程，高速、低載荷情形下計算結果更接近真實物理過程，得到的碰撞力量值更加準確。

圖9 三種方法對折疊翼展開過程碰撞力的計算結果Fig.9 The calculation results of impact forces of folding wings by three methods

4 結 論

對含有接觸非線性的折疊翼展開系統內部構件間存在的接觸/碰撞過程，在現有的速度-沖量形式非光滑動力學模型基礎上，推導了模型的線性互補形式?；诰€性互補形式的非光滑動力學模型，借助互補原理能夠以較少的計算量檢測和判斷觸-離、滯-滑狀態的切換，節省事件檢測所需的計算資源。

文中提出了新的隱式積分算法，用于迭代求解線性互補形式的非光滑動力學模型，并基于該隱式積分算法構造了基于速度和沖量的時間步進算法。文中所提積分算法基于隱式歐拉法，用前一時刻的真實速度和后一時刻的近似速度，加權構造后一時刻的位移，更新后一時刻各狀態量，算出后一時刻的真實速度。相較于Moreau中點算法，文中提出的隱式積分算法具有更高的計算精度，使得基于該算法的時間步進算法在求解包含單點碰撞的非光滑過程時，具有更高的效率。文末算例檢驗了集成新型積分算法的時間步進算法的精度與穩定性。