分离装置又称为分离器, 是一种典型的分离及连接机构, 广泛应用于航天系统中的级间分离、舱段分离、星箭分离和载荷释放等相关领域。其引导阶段为含摩擦多点接触碰撞非光滑动力学过程^[1]。以往, 工程师们多采用有限元或商业多体动力学软件对其进行建模分析, 但仿真计算时间长, 接触碰撞处理难, 参数不稳定, 难以指导实际工程设计与满足快速分析实际需求。Brogliato^[2]对非光滑动力学建模方法进行了系统总结与分析。刘才山等^[3]和赵振等^[4]在多体系统框架内, 针对含摩擦接触、碰撞等非光滑动力学行为开展了相关研究。在碰撞位形不变与常规作用力可忽略等合理假设下, Liu等^[5-6]提出在冲量与速度水平上处理多点接触、碰撞的新建模方法, 简称为LZB方法。与其他方法不同, LZB方法通过对多点碰撞接触中主要特征物理量的恰当描述以及对多尺度耦合关系的合理刻画, 能够实现短时间内对含摩擦接触与碰撞等复杂非光滑动力学问题的准确刻画。随着现代航天运载器技术的不断发展, 对运载器研制分析速度的要求也相应提高。分离装置是其中的一个关键部件, 对其含摩擦接触与碰撞等非光滑动力学快速准确的动力学建模的研究尤为重要, 也是实现快速发射、快速反应工程研制的必要条件。

本文首先基于LZB方法, 开展低冲击分离装置动力学建模, 然后进行算例验证, 证明LZB方法在工程研究分析中的可行性, 为分离装置引导阶段的动力学分析提供可行的工程模拟分析途径。

1 接触动力学描述

由于连接螺栓与分离管道都属于典型中心轴对称体, 本文动力学建模过程中将其简化为如图 1所示的平面动力学问题。连接螺栓头部半径为R, 长度为L; 螺栓螺杆部分半径为r, 长度为l。连接螺栓位于分离管道内部, 分离管道上部边界为y₁, 下部边界为y₂。连接螺栓密度为ρ, 质量为m。连接螺栓材料分布均匀, 其质心在其几何中心O点, 相对于O点绕k轴的转动惯量为J。

如图 1所示, 连接螺栓与分离管道之间存在6个潜在接触点: A, B, C, D, E, F。在分离管道左端面建立惯性坐标系{i, j, k}, 在O点建立质心固连坐标系{O, e₁, e₂, e₃}。惯性坐标下连接螺栓的广义坐标为$ \mathit{\boldsymbol{q}}{\rm{ = \{ }}{x_O}{\rm{, }}\;{y_O}{\rm{, }}\;\theta {{\rm{\} }}^{\rm{T}}} $, 其中x_O，y_O为连接螺栓质心在惯性坐标系中的坐标, θ为i与e₁之间的夹角, 初始状态时θ=0。系统的动能表示为

$ T{\rm{ = }}\frac{1}{2}m(v_{Ox}^2 + v_{Oy}^2) + \frac{1}{2}J{\omega ^2}, $

其中ω为连接螺栓的角速度, v_Ox和v_Oy为连接螺栓质心在i, j方向上的速度。取分离管道中心为零势面, 则连接螺栓的势能为

$ V{\rm{ = }}mg{y_O}。 $

惯性坐标系{i, j, k}与连接螺栓固连坐标系{O, e₁, e₂, e₃}之间的坐标变化矩阵为

$ \left[{\begin{array}{*{20}{c}} {{e_1}}\\ {{e_2}}\\ {{e_3}} \end{array}} \right] = \left[{\begin{array}{*{20}{c}} {\cos \theta }&{\sin \theta }&0\\ {-\sin \theta }&{\cos \theta }&0\\ 0&0&1 \end{array}} \right]\left[{\begin{array}{*{20}{c}} i\\ j\\ k \end{array}} \right]。 $

当θ=0, EF边与j轴重合时, 连接螺栓质心O点在惯性坐标系{i, j, k}中的初始位置为$ h = \frac{{{R^2}L + {r^2}l}}{{2{r^2}}} $。

6个潜在接触点在连接螺栓固连坐标系下的位置矢量为

$ \begin{array}{l} {\mathit{\boldsymbol{r}}_A}{\rm{ = (}}L{\rm{ + }}l-h{\rm{)}}{\mathit{\boldsymbol{e}}_1}-R{\mathit{\boldsymbol{e}}_2}, \;\;{r_B}{\rm{ = (}}L{\rm{ + }}l-h{\rm{)}}{\mathit{\boldsymbol{e}}_1} + R{\mathit{\boldsymbol{e}}_2}, \\ {\mathit{\boldsymbol{r}}_C}{\rm{ = (}}L - h{\rm{)}}{\mathit{\boldsymbol{e}}_1} - R{\mathit{\boldsymbol{e}}_2}, \;\;{r_D}{\rm{ = (}}L - h{\rm{)}}{\mathit{\boldsymbol{e}}_1} + R{\mathit{\boldsymbol{e}}_2}, \\ {\mathit{\boldsymbol{r}}_E}{\rm{ = (}} - h{\rm{)}}{\mathit{\boldsymbol{e}}_1} - r{\mathit{\boldsymbol{e}}_2}, \;\;{r_F}{\rm{ = (}} - h{\rm{)}}{\mathit{\boldsymbol{e}}_1} + r{\mathit{\boldsymbol{e}}_2}。\; \end{array} $

6个潜在接触点在惯性坐标系中的位置分量为

$ \begin{array}{l} {\mathit{\boldsymbol{r}}_A} = {\mathit{\boldsymbol{r}}_O} + ((L + l-h)\cos \theta + R\sin \theta )\mathit{\boldsymbol{i}} + \\ \;\;\;\;\;\;((L + l-h)\sin \theta-R\cos \theta )\mathit{\boldsymbol{j}}{\rm{, }}\\ {\mathit{\boldsymbol{r}}_B} = {\mathit{\boldsymbol{r}}_O} + ((L + l - h)\cos \theta - R\sin \theta )\mathit{\boldsymbol{i}} + \\ \;\;\;\;\;\;((L + l - h)\sin \theta + R\cos \theta )\mathit{\boldsymbol{j}}{\rm{, }}\\ {\mathit{\boldsymbol{r}}_C} = {\mathit{\boldsymbol{r}}_O} + ((l - h)\cos \theta + R\sin \theta )\mathit{\boldsymbol{i}} + \\ \;\;\;\;\;\;((l - h)\sin \theta - R\cos \theta )\mathit{\boldsymbol{j}}{\rm{, }}\\ {\mathit{\boldsymbol{r}}_D} = {\mathit{\boldsymbol{r}}_O} + ((l - h)\cos \theta - R\sin \theta )\mathit{\boldsymbol{i}} + \\ \;\;\;\;\;\;((l - h)\sin \theta + R\cos \theta )\mathit{\boldsymbol{j}}{\rm{, }}\\ {\mathit{\boldsymbol{r}}_E} = {\mathit{\boldsymbol{r}}_O} + ( - h\cos \theta + r\sin \theta )\mathit{\boldsymbol{i}} + \\ \;\;\;\;\;\;( - h\sin \theta - r\cos \theta )\mathit{\boldsymbol{j}}{\rm{, }}\\ {\mathit{\boldsymbol{r}}_F} = {\mathit{\boldsymbol{r}}_O} + ( - h\cos \theta - r\sin \theta )\mathit{\boldsymbol{i}} + \\ \;\;\;\;\;\;( - h\sin \theta + {\rm{r}}\cos \theta )\mathit{\boldsymbol{j}}\;。 \end{array} $

6个潜在接触点的速度为

$ \begin{array}{l} {\mathit{\boldsymbol{v}}_A} = {\mathit{\boldsymbol{v}}_O} + ((L + l-h)-l\;{\rm{sin}}\theta \cdot \dot \theta + R\cos \theta \cdot \dot \theta )\mathit{\boldsymbol{i}} + \\ \;\;\;\;\;\;\;((L + l-h)\cos \theta \cdot \dot \theta + R\sin \theta \cdot \dot \theta )\mathit{\boldsymbol{j}}, \\ {\mathit{\boldsymbol{v}}_B} = {\mathit{\boldsymbol{v}}_O} + ((L + l - h) - l\;{\rm{sin}}\theta \cdot \dot \theta - R\cos \theta \cdot \dot \theta )\mathit{\boldsymbol{i}} + \\ \;\;\;\;\;\;((L + l - h)\cos \theta \cdot \dot \theta - R\sin \theta \cdot \dot \theta )\mathit{\boldsymbol{j}}, \\ {\mathit{\boldsymbol{v}}_C} = {\mathit{\boldsymbol{v}}_O} + ((l - h) - l\;{\rm{sin}}\theta \cdot \dot \theta + R\cos \theta \cdot \dot \theta )\mathit{\boldsymbol{i}} + \\ \;\;\;\;\;\;((l - h)\cos \theta \cdot \dot \theta + R\sin \theta \cdot \dot \theta )\mathit{\boldsymbol{j}}, \\ {\mathit{\boldsymbol{v}}_D} = {\mathit{\boldsymbol{v}}_O} + ((l - h) - l\;{\rm{sin}}\theta \cdot \dot \theta - R\cos \theta \cdot \dot \theta )\mathit{\boldsymbol{i}} + \\ \;\;\;\;\;\;((l - h)\cos \theta \cdot \dot \theta - R\sin \theta \cdot \dot \theta )\mathit{\boldsymbol{j}}, \\ {\mathit{\boldsymbol{v}}_E} = {\mathit{\boldsymbol{v}}_O} + ( - h - l\sin \theta \cdot \dot \theta + r\cos \theta \cdot \dot \theta )\mathit{\boldsymbol{i}} + \\ \;\;\;\;\;\;( - h\cos \theta \cdot \dot \theta + r\sin \theta \cdot \dot \theta )\mathit{\boldsymbol{j}}, \\ {\mathit{\boldsymbol{v}}_F} = {\mathit{\boldsymbol{v}}_O} + ( - h - l\sin \theta \cdot \dot \theta - r\cos \theta \cdot \dot \theta )\mathit{\boldsymbol{i}} + \\ \;\;\;\;\;\;( - h\cos \theta \cdot \dot \theta - r{\rm{sin}}\theta \cdot \dot \theta )\mathit{\boldsymbol{j}}\;。 \end{array} $

考虑连接螺栓6个潜在接触点处的法线方向, A, C, E的外法线方向为n₁=-j, B, D, F的外法线方向为n₂=j。各潜在接触点对的法向相对位移为

$ \begin{array}{l} {\delta _A}(q, \;t) = {y_O} + ((L + l-h)\sin \theta-R\cos \theta )-{y_1}, \\ {\delta _B}(q, \;t) = {y_2} - {y_O} - ((L + l - h)\sin \theta + R\cos \theta ), \\ {\delta _C}(q, \;t) = {y_O} + ((l - h)\sin \theta - R\cos \theta ) - {y_1}, \\ {\delta _D}(q, \;t) = {y_2} - {y_O} - ((l - h)\sin \theta + R\cos \theta ), \\ {\delta _E}(q, \;t) = {y_O} + (( - h)\sin \theta - r\cos \theta ) - {y_1}, \\ {\delta _F}(q, \;t) = {y_2} - {y_O} - (( - h)\sin \theta + r\cos \theta )\;。 \end{array} $

对应各潜在接触点处的相对切向速度为

$ \begin{array}{l} {v_A} = {v_{Ox}} + ((L + l-h)-l\sin \theta \cdot \dot \theta + R\cos \theta \cdot \dot \theta ), \\ {v_B} = {v_{Ox}} + ((L + l-h) - l\sin \theta \cdot \dot \theta - R\cos \theta \cdot \dot \theta ), \\ {v_C} = {v_{Ox}} + ((l - h) - l\sin \theta \cdot \dot \theta + R\cos \theta \cdot \dot \theta ), \\ {v_D} = {v_{Ox}} + ((l - h) - l\sin \theta \cdot \dot \theta - R\cos \theta \cdot \dot \theta ), \\ {v_E} = {v_{Ox}} + ( - h - l\sin \theta \cdot \dot \theta + r\cos \theta \cdot \dot \theta ), \\ {v_F} = {v_O} + ( - h - l\sin \theta \cdot \dot \theta - r\cos \theta \cdot \dot \theta )\;。 \end{array} $

连接螺栓动力学方程为

$ \mathit{\boldsymbol{M}}{\rm{(}}\mathit{\boldsymbol{q}}{\rm{)}}\mathit{\boldsymbol{\ddot q}} = \mathit{\boldsymbol{G}} + \mathit{\boldsymbol{H}} + \mathit{\boldsymbol{W}}{F^n}-\mathit{\boldsymbol{N}}{F^\tau }, $

式中, M(q)为质量矩阵:

$ \mathit{\boldsymbol{M}}{\rm{(}}\mathit{\boldsymbol{q}}{\rm{) = }}\left[{\begin{array}{*{20}{c}} m&0&0\\ 0&m&0\\ 0&0&J \end{array}} \right]。 $

G为重力矩阵, H为向心力矩阵, W, N分别为雅克比法向矩阵和切向矩阵。

接触点处的切向与法向相对速度为

$ \begin{array}{l} \mathit{\boldsymbol{\dot \delta }} = {\mathit{\boldsymbol{W}}^{\rm{T}}}(\mathit{\boldsymbol{q}})\mathit{\boldsymbol{\dot q}}, \\ {\mathit{\boldsymbol{v}}^\tau } = {\mathit{\boldsymbol{N}}^{\rm{T}}}(\mathit{\boldsymbol{q}})\mathit{\boldsymbol{\dot q}}。\; \end{array} $

对上式微分, 则潜在接触点处的相对加速度在法向和切向的分量为

$ \begin{array}{l} \mathit{\boldsymbol{\ddot \delta }} = {\mathit{\boldsymbol{W}}^{\rm{T}}}(\mathit{\boldsymbol{q}})\mathit{\boldsymbol{\ddot q}} + {\mathit{\boldsymbol{S}}_n}, \\ {{\mathit{\boldsymbol{\dot v}}}^\tau } = {\mathit{\boldsymbol{N}}^{\rm{T}}}(\mathit{\boldsymbol{q}})\mathit{\boldsymbol{\ddot q}} + {\mathit{\boldsymbol{S}}_\tau }。 \end{array} $

结合螺栓动力学方程, 则上式可表示为

$ \begin{array}{l} \mathit{\boldsymbol{\ddot \delta }} = {\mathit{\boldsymbol{W}}^{\rm{T}}}{\mathit{\boldsymbol{M}}^{-1}}(\mathit{\boldsymbol{q}})\mathit{\boldsymbol{W}}{F^n}-{\mathit{\boldsymbol{W}}^{\rm{T}}}{\mathit{\boldsymbol{M}}^{-1}}(\mathit{\boldsymbol{q}})\mathit{\boldsymbol{N}}{F^\tau } + {\mathit{\boldsymbol{S}}_n}, \\ {{\mathit{\boldsymbol{\dot v}}}^\tau } = {\mathit{\boldsymbol{N}}^{\rm{T}}}{\mathit{\boldsymbol{M}}^{ - 1}}(\mathit{\boldsymbol{q}})\mathit{\boldsymbol{W}}{F^n} - {\mathit{\boldsymbol{N}}^{\rm{T}}}{\mathit{\boldsymbol{M}}^{ - 1}}(\mathit{\boldsymbol{q}})\mathit{\boldsymbol{N}}{F^\tau } + {\mathit{\boldsymbol{S}}_\tau }。 \end{array} $

连接螺栓与分离管道之间的6个潜在接触点可分为闭合接触点(δ_j=0)和非闭合接触点(δ_j > 0)。对于非闭合接触点处, 连接螺栓与分离管道之间的相互作用力为零。对于闭合接触点(δ_j=0), 需根据接触点的相对法向运动速度判断接触点状态。若$ {\dot \delta _j} $>0, 说明在该接触点的两接触体在该时刻分离, 相互作用力为零。若存在任意一个闭合接触点处的法向相对速度$ {\dot \delta _j} $ < 0, 说明该接触点将出现一个碰撞过程, 从而使连接螺栓与分离管道之间的接触状态转入碰撞过程。如果连接螺栓与分离管道之间不存在任何一个碰撞接触点(此时$ {\dot \delta _j} $=0), 说明对应的接触点将保持接触状态。接触状态可分为滑动与黏滞两种状态, 滑动状态下接触点的法向与切向作用力需满足$ F_j^\tau = \mu F_j^n $; 黏滞状态下需满足$ F_j^\tau \le \mu F_j^n $。

2 碰撞动力学描述

由于连接螺栓与分离管道多为硬性金属材料, 连接螺栓分离后的速度也小于10 m/s, 连接螺栓与分离管道之间一旦进入碰撞状态, 其动力学过程具有以下特征: 1) 在很短时间内产生很大的冲击力, 但冲击力的冲量为有限值; 2) 冲击前后速度发生突变, 但位形几乎不变。

因此, 分离装置连接螺栓碰撞过程符合LZB方法的假设前提: 1) 碰撞前后连接螺栓与分离管道位形不发生变化; 2) 碰撞过程中由于冲击力相对其他力很大, 碰撞动力学分析中忽视重力与向心力等常规作用力。

设法向作用力与切向作用力对应的冲量为

$ {\rm{d}}{P^n} = {F^n}{\rm{d}}t, \;\;{\rm{d}}{P^\tau } = {F^\tau }{\rm{d}}t, $

则连接螺栓与分离管道之间的接触动力学可转化为碰撞动力学:

$ \mathit{\boldsymbol{M}}{\rm{d}}\dot q = \mathit{\boldsymbol{W}}(\mathit{\boldsymbol{q}}){\rm{d}}{P^n} + \mathit{\boldsymbol{N}}(\mathit{\boldsymbol{q}}){\rm{d}}{P^\tau }。 $

实际工程设计中, 连接螺栓各潜在接触点处都做倒圆角处理, 则各潜在接触点处法向接触力F_jⁿ与局部弹性变形刺穿深度δ_j的关系可近似地按Hertz接触定义为$ F_j^n = {k_j}\delta _j^{\frac{3}{2}} $, 对应碰撞过程中接触点处的弹性势能变化量为

$ {\rm{d}}{E_j} =-F_j^n{\rm{d}}{\delta _j} =-\frac{{{\rm{d}}P_j^n}}{{{\rm{d}}t}}{\rm{d}}{\delta _j} =-{\dot \delta _j}{\rm{d}}P_j^n。 $

Liu等^[5]提出可利用各点存储势能获得碰撞过程中不同点之间法向冲量的关系式:

$ \frac{{{\rm{d}}P_i^n}}{{{\rm{d}}P_j^n}} = {\left( {\frac{{{k_i}}}{{{k_j}}}} \right)^{\frac{2}{5}}}{\left( {\frac{{{E_i}}}{{{E_j}}}} \right)^{\frac{1}{5}}} = R_{i, j}^n。 $

同时, 各碰撞点处的切向冲量与法向冲量之间仍然满足库仑摩擦定律, 即滑动状态下j点法向冲量与切向冲量满足$ {\rm{d}}P_j^\tau =-\mu {\rm{d}}P_j^n{\rm{sign}}(v_j^\tau ) $, 黏滞状态下满足$ \left| {\frac{{{\rm{d}}P_j^\tau }}{{{\rm{d}}P_j^n}}} \right| \le {\mu _{\rm{s}}} $。

使用Stronge能量恢复系数e_j对碰撞过程中各个碰撞点处能量的耗散进行刻画。碰撞过程中各点能量变化量随冲量步长变化关系为

$ {\rm{d}}{E_i}(P_j^n) =-\frac{1}{{{\eta _i}(\mathop {{\delta _i}}\limits^. )}}\mathop {{\delta _i}}\limits^. R_{i, j}^n{\rm{d}}P_j^n, $

$ {\eta _i}(\mathop {{\delta _i}}\limits^. ) = \left\{ {\begin{array}{*{20}{l}} {1, }&{当{{\dot \delta }_i} \le 0, }\\ {e_i^2, }&{当{{\dot \delta }_i} > 0。\;} \end{array}} \right. $

由此可以得出, 在碰撞点之间相对法向速度不大于零的压缩阶段, 接触点处的弹性势能不断增加。反之, 在碰撞点之间的相对法向速度大于零的恢复阶段, 接触点处的弹性势能不断减少, 若减少到零, 此对应点的碰撞结束。

3 算例

取螺帽半径R=14 mm, 厚度L=12 mm, 螺栓半径r=10 mm, 长度l=34 mm。分离管道内径20mm, 即y₁=10 mm, y₂=-10 mm。螺栓密度r=7.8×10^-3 g/mm³, 重力加速度g=9.8 m/s。各潜在接触点处的Stronge恢复系数为e=0.97, 对应接触点处的摩擦系数为μ=0.3, μ_s=0.5。螺栓初始位形如图 1所示, 初始速度为$ {{\mathit{\boldsymbol{\dot q}}}^0} = {[1\;{\rm{m/s}}, \;-0.2\;{\rm{m/s}}, -50\;{\rm{rad/s}}]^{\rm{T}}} $。

图 2为螺栓在管道运动过程中6个潜在接触点法向位置随时间变化曲线。螺栓自由运动后, 先在A点与分离管道发生第一次碰撞, 然后依次在F点、B点、E点、C点发生4次碰撞。由每次碰撞时各潜在接触点与分离管道内壁法向距离可知, 每次碰撞过程皆为单点碰撞过程。

5次碰撞过程中, 每次碰撞螺栓能量都发生一次衰减, 如图 3所示。碰撞前后螺栓都为自由状态, 与分离管道内壁无接触, 对应自由状态下螺栓能量守恒不变。如图 4所示, 每次碰撞前后, 碰撞点的法向速度也发生变化。以A点碰撞为例, 碰撞前A点法向速度为负值, 碰撞后变为正值。潜在接触点A由碰撞点变为自由状态点。碰撞过程中A, B, E, F处的弹性势能如图 5所示。由于每次碰撞都为单点碰撞, 即每次碰撞过程中除碰撞点外, 其余点处弹性势能为零。每次碰撞时, 碰撞点处的势能先在压缩阶段增加, 再在释放阶段减少。

4 结论

本文基于LZB方法, 对分离装置引导阶段的动力学过程进行系统研究。在分离装置引导阶段, 连接螺栓与分离管道之间存在含摩擦、接触、碰撞、自由等多个状态, 且运动过程中状态之间还会发生切换, 属于复杂非光滑动力学过程。使用常规电脑, 在MATLAB软件平台上, 针对连接螺栓与分离管道之间的4次碰撞动力过程, 仿真计算仅需不到一分钟, 远小于有限元软件服务器仿真小时级计算时间。算例的数值计算结果证明, LZB方法能够较好地描述分离装置引导阶段非光滑动力学过程, 除接触、自由等常规状态外, 还能对碰撞过程进行详细的分析。分离装置引导阶段涉及许多复杂的力学问题, 还需进一步发展相关理论并完善数值计算方法。