CN106066922A - 一种三滑块单轨柔性火箭橇的垂向响应计算方法 - Google Patents

Abstract

本发明公开了一种三滑块单轨柔性火箭橇的垂向响应的计算方法。该方法首先对火箭橇进行受力分析，将箭体离散为有限元梁模型；进而考虑轨道不平顺引起火箭橇与轨道结构振动，采用Hertz非线性接触模型，计算得到箭体各时刻受到的碰撞接触力；最终通过中心差分的数值方法，计算得到给定时间内三滑块单轨火箭橇系统全程的垂向动力学响应。本发明在计算过程中考虑了火箭橇橇体的柔性变形特性及火箭橇系统由于喷出燃料导致的质量折减效应，并且考虑了轨道不平稳引起的非线性接触问题，因此可以得到较为可靠的动力学响应结果。

Description

一种三滑块单轨柔性火箭橇的垂向响应计算方法

技术领域

本发明涉及一种三滑块单轨柔性火箭橇的垂向响应计算方法，涉及采用hertz接触模型求解碰撞接触力以及采用Newmark算法进行动力学数值计算，从而实现火箭橇在不平顺轨道上动态接触力求解和垂向动力学响应的高效和快速预测的计算方法。

背景技术

火箭橇是采用火箭发动机作为动力，沿着专门建造的轨道运行的一种试验工具。火箭橇试验受试验件体积、外形和质量的影响较小，大至数百千克的全尺寸部件、小至数十克的缩比模型都可以在火箭橇上进行试验。在轨道足够长的情况下，采用不同推力的发动机或不同数量发动机组合作为橇车的动力源时，橇车的速度可以在几十米每秒至上千米每秒的范围内进行调节，而且火箭橇试验能够很方便的观察试验结果和重复收集试验数据。因此，在先进常规武器系统研制(如导弹碰撞和侵彻试验、子母弹抛撒试验、惯性制导系统试验等)、飞机弹射座椅的开发、空气动力学试验等方面有极为广泛的应用。

第二次世界大战后，美国、英国、法国和前苏联等国家为了推动武器、航空航天的发展，相继建造了多种不同类型的火箭橇试验场，开展了包括惯性测量装置在内的一些核心部件的火箭橇综合验证试验工作。火箭橇试验的显著特点是可无损回收被试惯性测量装置，供进一步测量、检查及继续进行试验。高精度的惯性测量装置造价高，通过火箭橇试验可重复进行多类多次的测试试验，包括环境适应性试验和精度试验，增加试验样本量，确保飞行试验一次成功，减少试验次数，降低试验成本，加快研制周期。美国在应用火箭橇对惯性测量装置进行试验研究方面发展最快，应用最广。从上世纪四十年代末到七十年代初，美国相继建造了多达25条的不同类型的火箭橇试验轨道，先后开展了大量的火箭橇试验项目研究。

火箭橇约束在轨道内高速运动的过程中，受到火箭推力、碰撞接触力、空气阻力和制动力等作用，火箭橇与轨道的动态性能直接影响到车载试验件的试验环境，将直接决定动态试验结果的准确性和可靠性。由于轨道长度与轨道不平顺的数量级相差巨大，并且橇轨在轨道不平顺点接触所产生的附加接触变形量的精确解往往是通过多次迭代计算所求得，导致有限元建模和计算都十分困难。

开展高速火箭橇与轨道系统动态性能的研究，建立橇-轨动力学数学模型，有助于提高 对火箭橇、轨道动力学特性的认识，揭示系统物理本质，为研制满足要求、性能良好的试验是被提供设计依据，最终使火箭橇及试验轨道实际更为优化可靠。对于高速火箭橇-轨道耦合系统动力学仿真分析研究，具有重要的科学使用意义和工程应用价值。

目前尚未有相关技术报导。

发明内容

本发明解决的技术问题是：本发明克服现有技术建模困难、计算量巨大的不足，提供一种三滑块单轨柔性火箭橇的垂向响应计算方法，求解速度快，工程适用性强，在计算过程中考虑了火箭橇橇体的柔性变形特性及火箭橇系统由于喷出燃料导致的质量折减效应，并且考虑了轨道不平稳引起的非线性接触问题，因此可以得到较为可靠的动力学响应结果。

本发明采用的技术方案，一种三滑块单轨柔性火箭橇的垂向响应计算方法，基于Hertz接触力模型，实现步骤如下：

第一步：根据三滑块单轨柔性火箭橇的运动形式，将运动等效为航向的刚体平动和垂向的振动；三滑块单轨柔性火箭橇受到航向的推力和垂向接触力，根据模型受力形式，自行构造一种5自由度的梁单元，即在欧拉梁单元的基础上附加一个航向刚体自由度，将三滑块单轨柔性火箭橇模型离散为3个5自由度的梁单元模型，得出三滑块单轨柔性火箭橇模型的动力学矩阵M、C和K，建立运动学方程为：

$M\stackrel{\·\·}{u}+C\stackrel{\·}{u}+Ku=F\left(t\right)$

M、C和K分别为整体质量阵、刚度阵和阻尼阵，F(t)为模型所受外力；

第二步：根据第一步求得三滑块单轨柔性火箭橇模型，然后确定模型的边界条件，所述边界条件包括载荷与约束条件；定义迭代步长ΔT，已知初始时刻即T＝0时刻的航向位移s_x、航向速度v_x及航向加速度a_x和外载荷曲线，所述外载荷曲线为三滑块单轨柔性火箭橇推力、阻力和质量折减曲线，为初始已知条件，将时刻T代入外载荷曲线，时刻T由步长ΔT和迭代步数n的乘积决定，得到对应时刻的外载荷与质量，根据Newmark算法，得到下一时刻T+ΔT的航向位移、航向速度及航向加速度，然后根据迭代步数n求得任意时刻的航向位移、航向速度及航向加速度；

第三步：根据第二步得到T时刻航向位移s_x确定三滑块单轨柔性火箭橇在轨道的距离起点即T＝0时刻所在位置，从而确定轨道不平顺信息即实测得到的轨道各测点高度与基准面高度的差，结合对应时刻T的三滑块单轨柔性火箭橇的垂向位移、速度和加速度，计算接触变形与接触相对方向速度；

第四步：滑靴相当于车轮，为火箭橇运动时与火箭橇轨道相接触的部分。确定对应T的 滑靴与轨道的接触状态，确定碰撞接触力大小，当滑靴在轨道上运动时，轨道与滑靴之间存在着三种状态：①滑靴与轨道上表面接触、②滑靴与轨道下表面接触、③滑靴与轨道不接触；接触状态不同，碰撞接触力的计算公式有所不同，首先根据轨道与滑靴的初始间隙以及当前时间步滑靴位置，确定出当前时刻的接触状态，再根据hertz接触力模型计算碰撞接触力；

第五步：根据第四步计算的碰撞接触力以及对应时刻T三滑块单轨柔性火箭橇的垂向位移、速度和加速度作为输入，施加到三滑块单轨柔性火箭橇模型上，采用Newmark方法求解出下一时刻的动力学响应，即为T+ΔT时刻的垂向位移、速度和加速度；

第六步：将第五步得到T+ΔT动力学响应作为输入，回到第四步，确定对应时刻T+ΔT的滑靴与轨道的接触状态，确定碰撞接触力大小。应用第五步中所述，将第五步得到T+ΔT动力学响应施加到三滑块单轨柔性火箭橇模型上，采用Newmark方法，从而求得T+2ΔT的动力学响应，循环至时间步长达到所要求时间步为止，最终得到三滑块单轨柔性火箭橇在给定时间范围内的动力学响应，即垂向位移、垂向速度以及垂向加速度。

本发明与现有技术相比的优点在于：本发明以轨道不平顺数值、轨道与滑靴初始间隙、火箭橇推力曲线为输入条件，将火箭橇离散为5自由度梁单元，采用Newmark数值算法，应用Hertz接触力模型，根据当前时间步的垂向位移和速率更新对应时间步的碰撞接触力，最终求得各有限元节点的垂向动力学响应，具有求解速度快，工程适用性强，在计算过程中考虑了火箭橇橇体的柔性变形特性及火箭橇系统由于喷出燃料导致的质量折减效应，并且考虑了轨道不平稳引起的非线性接触问题，因此可以得到较为可靠的动力学响应结果。

附图说明

图1为本发明实现总流程图；

图2为本发明中的5自由度梁单元模型图；

图3为本发明中滑靴与轨道接触状态种类；

图4为本发明中三滑块单轨柔性火箭模型；

图5为本发明中三滑块单轨柔性火箭橇简化模型；

图6为本发明中推力曲线；

图7为本发明中实测轨道不平顺值；

图8为本发明中节点1垂向加速度响应(正弦不平顺值下)；

图9为本发明中节点3垂向加速度响应(正弦不平顺值下)；

图10为本发明中节点5垂向加速度响应(正弦不平顺值下)；

图11为本发明中节点1垂向加速度响应；

图12为本发明中节点3垂向加速度响应；

图13为本发明中节点5垂向加速度响应。

具体实施方式

如图1所示，本发明提出了一种基于hertz接触模型的三滑块单轨柔性火箭橇垂向响应计算方法。包括以下步骤：

(1)根据三滑块单轨柔性火箭橇的运动形式，将运动等效为航向的刚体平动和垂向的振动；三滑块单轨柔性火箭橇受到航向的推力和垂向接触力，根据模型受力形式，自行构造一种5自由度的梁单元，即在欧拉梁单元的基础上附加一个航向刚体自由度，将三滑块单轨柔性火箭橇模型离散为3个5自由度的梁单元模型，如图2所示，得出三滑块单轨柔性火箭橇模型的动力学矩阵M、C和K，建立运动学方程为：

$M\stackrel{\·\·}{u}+C\stackrel{\·}{u}+Ku=F\left(t\right)$

M、C和K分别为整体质量阵、刚度阵和阻尼阵，F(t)为模型所受外力；

对于航向刚体运动，三滑块单轨柔性火箭橇航向运动简化为杆的刚体平动，故其运动方程可由其质心运动方程表达，由牛顿第二定律得到：

$m{\stackrel{\·\·}{u}}_1=F_1+F_2$

其中，m为箭体质量，F₁和F₂分别为作用在5自由度梁单元两个节点处x方向上的合力。

对于垂向振动，梁垂向振动有限元方程可写为以下形式：

$M_1\stackrel{\·\·}{y}+C_1\stackrel{\·}{y}+K_{1}y=F\left(t\right)$

其中M₁、K₁、C₁分别为质量阵、刚度阵和阻尼阵，u＝[y₁ θ₁ y₂ θ₂]^T为节点纵向位移。F(t)外力向量。

综合考虑梁的垂向振动和航向刚体运动，三滑块单轨柔性火箭橇的运动方程可写为：

$M\stackrel{\·\·}{u}+C\stackrel{\·}{u}+Ku=F\left(t\right)$

其中M、K、C分别为整体质量阵、刚度阵和阻尼阵。

(2)根据第一步求得三滑块单轨柔性火箭橇模型，然后确定模型的边界条件，所述边界条件包括载荷与约束条件；定义迭代步长ΔT，已知初始时刻即T＝0时刻的航向位移s_x、航向速度v_x及航向加速度a_x和外载荷曲线，所述外载荷曲线为三滑块单轨柔性火箭橇推力、阻力和质量折减曲线，为初始已知条件，将时刻T代入外载荷曲线，时刻T由步长ΔT和迭代步数n的乘积决定，得到对应时刻的外载荷与质量，根据Newmark算法，得到下一时刻T+ΔT的航向位移、航向速度及航向加速度，然后根据迭代步数n求得任意时刻的航向位移、 航向速度及航向加速度；

(3)根据第二步得到T时刻航向位移确定三滑块单轨柔性火箭橇在轨道的距离起点即T＝0时刻所在位置，从而确定轨道不平顺信息即实测得到的轨道各测点高度与基准面高度的差，其中不平顺值是根据轨道实际情况每隔一定距离测定得到的一系列数据点列，要想获得任意位置的不平顺值需要对原数据进行插值(如线性插值)，轨道任意点的不平顺值与不平顺斜率计算方法如下所示：

$s=\frac{(s_{n+1}-s_n)l_1}{L}---\left(1\right)$

$s^{\′}=\frac{s}{l_1}---\left(2\right)$

其中，s为滑靴当前位置的轨道不平顺值，s′为滑靴当前位置轨道不平顺斜率，s_n+1为滑靴位置后一个轨道监测点的监测值，s_n为滑靴位置前一个轨道监测点的监测值，l₁为滑块位置距离前一个观测点的长度，L为监测点间隔。

结合对应时刻三滑块单轨柔性火箭橇的垂向位移、速度和加速度，计算接触变形与接触相对方向速度。

(4)确定对应时刻T的滑靴与轨道的接触状态，确定碰撞接触力大小。当滑靴在轨道上运动时，轨道与滑靴之间存在着三种状态(如图3所示)：①滑靴与轨道上表面接触、②滑靴与轨道下表面接触、③滑靴与轨道不接触。接触状态不同，碰撞接触力的计算公式有所不同。首先根据轨道与滑靴的初始间隙以及当前时间步滑靴位置，确定出当前时刻的接触状态，再根据hertz接触力模型计算碰撞接触力。分为如下三种情况：

A.滑靴和轨道不接触。此时得到当前碰撞接触力为0；

B.滑靴与轨道下缘接触。此时，接触变形δ由当前时间步滑靴垂向位移、不平顺值和初始间隙作差得到。接触变形速率取为当前时间步的垂向速率。则接触力有下列公式求得：

$H\left(\mathrm{\δ}\right)=-{\mathrm{K\δ}}^n-{\mathrm{\μ\δ}}^n\stackrel{\·}{\mathrm{\δ}}---\left(3\right)$

C.滑靴与轨道上缘接触。此时，接触变形δ由当前时间步滑块垂向位移、不平顺值和初始间隙作差得到。接触变形速率取为当前时间步的垂向速率。则接触力有下列公式求得：

$H\left(\mathrm{\δ}\right)={\mathrm{K\δ}}^n+{\mathrm{\μ\δ}}^n\stackrel{\·}{\mathrm{\δ}}---\left(4\right)$

其中，H(δ)代表非线性接触力，Kδⁿ代表的是弹性接触部分，是冲击阻尼力部分， 为相对速度，μ是迟滞阻尼系数，其值与粘性阻尼系数、剪切和体积变形等有关。K和μ的值可由试验测定或根据数值模拟方法对其进行拟合得到。

(5)根据第四步计算的碰撞接触力以及对应时刻T三滑块单轨柔性火箭橇的垂向位移、速度和加速度作为输入，施加到三滑块单轨柔性火箭橇模型上，采用Newmark方法求解出下一时刻的动力学响应，即为T+ΔT时刻的垂向位移、速度和加速度。

(6)将第五步得到T+ΔT动力学响应作为输入，回到第四步，确定对应时刻T+ΔT的滑靴与轨道的接触状态，确定碰撞接触力大小。应用第五步中所述，将第五步得到T+ΔT动力学响应施加到三滑块单轨柔性火箭橇模型上，采用Newmark方法，从而求得T+2ΔT的动力学响应，循环至时间步长达到所要求时间步为止，最终得到三滑块单轨柔性火箭橇在给定时间范围内的动力学响应，即垂向位移、垂向速度以及垂向加速度。

实施例：

本发明实施例采用一个三滑靴的火箭橇模型，如图5所示。航向过载曲线如图6所示，前3秒为加速段。划分五个5自由度梁单元，试验件划分1个单元，发动机划分4个单元，示意图如图所示。各段截面参数以及碰撞参数如下表1所示，滑靴与轨道间的间隙取为1.5mm。

表1基本参数

利用Fortran语言按照图1的流程编写求解程序，轨道不平顺值采用两种情况：正弦不平顺值和实测获得。求解得到节点垂向加速度响应如下所述：

(1)对于正弦不平顺值，即轨道不平顺按正弦规律沿轨道坐标方向变化，幅值取1mm，波长取8m。得到各节点的加速度响应如8、9和10所示。可见三滑块单轨柔性火箭橇各节点垂向加速度响应随着时间的增加，有逐渐增大的趋势。

(2)对于实测得到的不平顺值(7)，得到各节点的加速度响应如1、图12和3所示。可见火箭橇各节点垂向加速度响应随着时间的增加，有逐渐增大的趋势，这一点与正弦不平顺值下的垂向响应相同，但是变化趋势不如正弦情况下均匀。这与实际测量值的趋势相符。

综上所述，本发明提出了一种柔性火箭橇垂向响应计算方法。首先建立火箭橇模型，将其离散为5自由度梁单元，进而组装成有限元动力学矩阵。再将不平顺值信息、推力曲线作为输入条件，对每个求解时间步进行推力与碰撞接触力的更新，并求解下一时间步的动力学响应。最终便得到全程的垂向动力学响应。

本发明未详细阐述部分属于本领域技术人员的公知技术。

Claims (3)

1.一种三滑块单轨柔性火箭橇的垂向响应计算方法，其特征在于实现步骤如下：

第一步：根据三滑块单轨柔性火箭橇的运动形式，将运动等效为航向的刚体平动和垂向的振动；三滑块单轨柔性火箭橇受到航向的推力和垂向接触力，根据模型受力形式，构造一种5自由度的梁单元，即在欧拉梁单元的基础上附加一个航向刚体自由度，将三滑块单轨柔性火箭橇模型离散为3个5自由度的梁单元模型，得出三滑块单轨柔性火箭橇模型的动力学矩阵M、C和K，建立运动学方程为：

$M\stackrel{\·\·}{u}+C\stackrel{\·}{u}+Ku=F\left(t\right)$

M、C和K分别为整体质量阵、刚度阵和阻尼阵，F(t)为模型所受外力；

第二步：根据第一步求得三滑块单轨柔性火箭橇模型，确定该模型的边界条件，所述边界条件包括载荷和约束条件；定义迭代步长ΔT，已知初始时刻即T＝0时刻的航向位移、航向速度及航向加速度和外载荷曲线，所述外载荷曲线为三滑块单轨柔性火箭橇推力、阻力和质量折减曲线，为初始已知条件，将时刻T代入外载荷曲线，时刻T由步长ΔT和迭代步数n的乘积决定，得到对应时刻的外载荷与质量，根据Newmark算法，得到下一时刻T+ΔT的航向位移、航向速度及航向加速度，然后根据迭代步数n求得任意时刻的航向位移、航向速度及航向加速度；

第三步：根据第二步得到T时刻航向位移，确定三滑块单轨柔性火箭橇在轨道的距离起点即初始时刻T＝0时所在位置，从而确定轨道不平顺信息即实测得到的轨道各测点高度与基准面高度的差，结合对应时刻T的三滑块单轨柔性火箭橇的垂向位移，计算接触变形与接触相对方向速度；

第四步：确定对应时刻T的滑靴与轨道的接触状态，确定碰撞接触力大小；所述滑靴相当于车轮，为三滑块单轨柔性火箭橇运动时与三滑块单轨柔性火箭橇轨道相接触的部分；当滑靴在轨道上运动时，轨道与滑靴之间存在着三种状态：①滑靴与轨道上表面接触；②滑靴与轨道下表面接触；③滑靴与轨道不接触；接触状态不同，碰撞接触力的计算公式有所不同，首先根据轨道与滑靴的初始间隙以及当前时间步滑靴位置，确定出当前时刻的接触状态，再根据hertz接触力模型计算碰撞接触力；

第五步：根据第四步计算的碰撞接触力以及对应时刻T三滑块单轨柔性火箭橇的垂向位移、速度和加速度作为输入，施加到三滑块单轨柔性火箭橇模型上，采用Newmark方法求解出下一时刻的动力学响应，即为T+ΔT时刻的垂向位移、速度和加速度；

第六步：将第五步得到T+ΔT动力学响应作为输入，回到第四步，确定对应时刻T+ΔT的滑靴与轨道的接触状态，确定碰撞接触力大小，应用第五步中所述，将第五步得到T+ΔT动力学响应施加到三滑块单轨柔性火箭橇模型上，采用Newmark方法，从而求得T+2ΔT的动力学响应，循环至时间步长达到所要求时间步为止，最终得到三滑块单轨柔性火箭橇在给定时间范围内的动力学响应，即垂向位移、垂向速度以及垂向加速度。

2.根据权利要求书1中一种三滑块单轨柔性火箭橇的垂向响应计算方法，其特征在于：所述第一步中得出三滑块单轨柔性火箭橇模型的动力学矩阵M、C和K过程如下：

将三滑块单轨柔性火箭橇模型离散为5自由度单元，5自由度梁单元的自由度为两端节点的挠度和转角y₁，θ₁，y₂，θ₂以及航向的刚体平动自由度u₁，对于航向刚体平动，三滑块单轨柔性火箭橇航向运动方程可由其质心运动方程表达，由牛顿第二定律得到：

$m{\stackrel{\·\·}{u}}_1=F_1+F_2$

其中，m为箭体质量，F₁和F₂分别为作用在5自由度梁单元两个节点处x方向上的合力；

对于垂向振动，5自由度梁单元垂向振动有限元方程写为以下形式：

$M_1\stackrel{\·\·}{y}+C_1\stackrel{\·}{y}+K_{1}y=F\left(t\right)$

其中M₁、K₁、C₁分别为质量阵、刚度阵和阻尼阵，u＝[y₁ θ₁ y₂ θ₂]^T为节点纵向位移，F(t)外力向量，

综合考虑梁的垂向振动和航向刚体运动，三滑块单轨柔性火箭橇的运动方程为：

$M\stackrel{\·\·}{u}+C\stackrel{\·}{u}+Ku=F\left(t\right)$

其中M、K、C分别为整体质量阵、刚度阵和阻尼阵。

3.根据权利要求书1中一种三滑块单轨柔性火箭橇的垂向响应计算方法，其特征在于：第四步中，所述hertz接触力模型为：

$H\left(\mathrm{\δ}\right)={\mathrm{K\δ}}^n+{\mathrm{\μ\δ}}^n\stackrel{\·}{\mathrm{\δ}}$

其中，H(δ)代表非线性接触力，Kδⁿ代表的是弹性接触部分，是冲击阻尼力部分，为相对速度，μ是迟滞阻尼系数，μ值与粘性阻尼系数、剪切和体积变形有关，K和μ的值由试验测定或根据数值模拟方法进行拟合得到。