BAT子弹药折叠翼展开冲击非光滑模型的隐式积分算法

张洪铭，顾晓辉，孙 丽

(1.南京理工大学 机械工程学院，南京 210094；2.江苏科技大学 机械工程学院，江苏 镇江 212100)

折叠翼结构能有效缩小中近程导弹、制导炮弹的储运空间，利于增强车辆、战机的运载能力，提高武器系统的战术性能，在现代飞行器的设计中获得了广泛的应用[1]。不同于固定翼结构近似刚性的连接，折叠翼展开机构中普遍包含铰链、限位挡块、锁定插销及驱动装置。活动件的存在降低了组件间的连接强度和系统的结构刚度，在迅速展开过程中，驱动加载的瞬间以及翼面展开到位时刻的碰撞、锁定动作，均会在部件间产生冲击载荷，存在使构件被破坏的可能性。另外，翼面的展弦比较大时，展开末段的碰撞冲击会引起翼面振颤，不恰当的设计可能引发翼面与弹体共振，加剧结构破坏的程度或影响飞行器的运动轨迹。

因此，针对折叠翼展开的动态过程，冲击载荷特性的分析与优化是提升展开性能的关键。现有研究对折叠翼展开过程进行了广泛研究。林三春等[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)

非光滑动力学角度下的互补理论是对接触/碰撞等非线性过程中物理量所遵从规律的描述。基于互补理论，对系统内满足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中点算法，文中提出的隐式积分算法具有更高的计算精度，使得基于该算法的时间步进算法在求解包含单点碰撞的非光滑过程时，具有更高的效率。文末算例检验了集成新型积分算法的时间步进算法的精度与稳定性。