接触刚度的计算
step(time,0,0d,0.68,-12000d)+step(time,0.68,0d,1.77,0d)+step(time,1.77,0d,2.45,12000d)
3.2.3定义齿轮啮合的接触碰撞力
为了保证仿真分析的真实性，齿轮之间的啮合运动关系没有被定义成理想化的几何约束关系，而是被定义为基于接触碰撞的力约束关系，即齿轮之间只能通过接触碰撞力（法向）和摩擦力（切向）相互约束，而不存在其他的约束关系。
在ADAMS 中有两种接触碰撞的计算模型，一种是基于Hertz 理论的Impact 函数模型，一中是基于恢复系数（Coefficient of restitution）的泊松（POISSON ）模型。两种力模型都来自于法向接触约束的惩罚函数。ADAMS/C++Solver使用惩罚因子来转换所有的接触约束。
采用Impact 函数来计算各啮合齿轮轮齿之间的接触碰撞力。Impact 函数模型将实际中物体的碰撞过程等效为基于穿透深度的非线性弹簧—阻尼模型，其计算表达式为：
. ⎧n ⎧⎫⎪Max ⎨K (x 1-x )-step (x , x 1-d , C max , x 1, 0)x , 0⎬F =⎨⎩⎭⎪0⎩x x 1 （3-1）
其中 K ——接触刚度系数；
x 1——位移开关量，用于确定单侧碰撞是否起作用；
x ——接触物体之间的实测位移变量；
d ——阻尼达到最大时两接触物体的穿透深度；
C max ——最大接触阻尼； . x ——穿透速度；
n ——非线性弹簧力指数。
当x >x 1时，两物体不发生接触，接触力为0，当x
（1）弹性分量K (x 1-x ) n ，相当于一个非线性弹簧；
（2）阻尼分量step (x , x 1-d , C max , x 1, 0)x ，其方向与运动方向相反，为了避免阻尼分量突变而使得函数变得不连续，采用了阶跃函数step ()来定义阻尼，step ()函数是利用三次多项式逼近海赛（Heacisde ）阶跃函数，具有连续的一阶导数，但在起始点处二阶导数不连续。在ADAMS 中的表达形式为： .
step (x , x 0, h 0, x 1, h 1)
x ≤x 0⎧h 0⎪2=⎨h 0+(h 1-h 0) [(x -x 0) /(x 1-x 0) ]{3-2[(x -x 0) /(x 1-x 0) ]}x 0
其中，x 为自变量，当x 小于x 0时，因变量的值为初始值h 0，当x 大于x 1时，因变量的值为终止值h 1；当x 在初始值和终止值之间变化时，因变量根据一定规律光滑过渡，避免出现数值过渡突变、微分值不连续。由此定义的碰撞阻尼系数，如图3-5所示。
图3-5碰撞过程中的阻尼变化曲线
（1） 弹性分量计算
物体接触刚度系数与物体的材料属性和接触表面的几何形状有关，在此根据文献[1]提供的接触刚度计算式来计算各齿轮啮合的齿廓面接触刚度，计算式为：
⎡R 1R 2⎤4 K =⎢⎥3π(h 1+h 2) ⎣R 1+R 2⎦1/2 （3-3）
其中，R 1，R 2——相啮合的两齿廓面在啮合点处的曲率半径。对于渐开线
齿轮，其工作过程中啮合点在齿廓面上的位置是不断变化
的，其曲率半径也是不断变化的，主动轮齿廓面上啮合点
处的曲率半径由小变大，从动轮尺阔面上啮合点处的曲率
半径由大变小，因此以两轮齿在节点处啮合作为计算点，mz 2，m 为模数，z 为啮合齿轮的齿则R i =i sin α；i =1，2
数，α为节圆压力角，对于标准啮合，节圆压力角等于分
度圆压力角20o 。
1-u i 22，u 为泊松比，E 为h 1，h 2——材料参数，定义为：h i =；i =1，πE i
弹性模量
非线性指数n ，根据Hertz 理论，一般取1.5较合适[2]。
（2） 阻尼分量计算
采用文献[2]给出的非线性阻尼模型来计算轮齿啮合接触阻尼系数，该阻尼模型认为物体表面接触—碰撞过程中的能量损失是由接触阻尼引起的，在基于等效能量损失的基础上给出接触阻尼的计算式：
3K (1-e 2) n δ （3-4） C = 4U
其中，K ——接触刚度；
e ——弹性恢复系数，一般定义为碰撞前法向速度差值与碰撞后法向速
度差值的比值。跟物体的材料、碰撞表面曲率半径、碰撞速度、
以及润滑介质的粘度有关，一般通过实验测定。
δ——穿透深度，对应ADAMS 取最大阻尼系数时的穿透深度，在此取
δ=0. 1mm ；
n ——非线性指数，一般取1.5[2]。
U ——碰撞速度，以相啮合的两个齿轮在节点处的线速度的差值代替。
在齿廓面间的动态碰撞力的作用下，相啮合的两个齿轮在节点处的
线速度是不等的，且随时间变化，即碰撞速度也不是定值
在阻尼分量中，弹性恢复系数和碰撞速度与实际的工况有关，其具体的取值在各仿真工况下确定。
（3） 摩擦力计算
接触体之间的摩擦力采用库仑摩擦模型，考虑静摩擦和动摩擦，有润滑时，取静摩擦系数为0.08，动摩擦系数为0.05，无润滑时，取静摩擦系数为0.15，动摩擦系数为0.1。在ADAMS 中摩擦力采用下面的函数表达式计算：
F = -N * step(v, -Vs, -1, Vs, 1) * step(ABS(v), Vs, Cst, Vtr, Cdy) （3-5） 其中：N ——法向力；
v ——表面相对滑移速度；
Vs ——最大静摩擦对应的相对滑移速度；
Cst ——静摩擦系数；
Vtr ——动摩擦对应的相对滑移速度；
Cdy ——动摩擦系数；
摩擦系数与相对滑移速度的关系，如图3-5所示。
图3-5 摩擦系数与相对滑移速度的关系曲线
4发射装置动力学仿真
4.1 ADAMS动力学分析算法及求解器的选择
ADAMS 对动力学微分方程，根据机械系统特性，选择不同的积分算法：对刚性系统，采用变系数的BDF （Backwards Differentiation Formulation）刚性积分程序，它是自动变阶、变步长的预估校正法（PECE ，Predict-Evaluate-Correct-Evaluate ），并分别有Index3、SI2、SI1三种积分格式，在积分的每一步采用了修正的Newton-Raphson 迭代算法，对应的求解器有GSTIFF 、WSTIFF 、DSTIFF 、CONSTANT_BDF；对高频系统，采用坐标分块法将微分—代数（DAE ）方程简化为常微分（ODE ）方程，分别利用ABAM （Adams-Bashforth and Adams-Moulton)方法和龙格—库塔（RKF45）方法求解，对应的求解器有ABAM 求解器和RKF45求解器。
4.1.1 BDF刚性积分法步骤
（1）预估阶段
用Gear 预估—校正算法可以有效地求解微分—代数方程。首先，根据当前时刻的系统状态变量，用泰勒级数预估下一时刻系统的状态矢量值：
y n +1∂y n 1∂2y n 2=y y +h +h + (4-1) ∂t 2！∂t 2
其中，时间步长h =t n +1-t n 。这种预估算法得到的新时刻的系统状态矢量值通常不准确，可以用Gear k+1阶积分求解程序（或其他向后差分积分程序）来
校正。
y n +1=-h β0y n +1+∑αi y n -i +1 （4-2） i =1. k
其中，y n +1为y (t ) 在t =t n +1时的近似值；β0和αi 为Gear 积分程序系统值。 上式经过整理，可表示为：
y n +1. k -1=[y n +1-∑αi y n -i +1] （4-3） i =1h β0
（2）校正阶段
① 求解系统方程G ，如G (y , y , t ) =0，则方程成立，此时的y 为方程的解，否则继续。
② 求解Newton-Raphson 线性方程，得到∆y ，以更新y ，使系统方程G 更接近于0。J ∆y =G (y , y , t n +1) ，其中J 为系统的雅可比矩阵。
③ 采用Neeton-Raphson 迭代，更新y :y k +1=y k +∆y k
④ 重复以上步骤直到∆y 足够小。
（3）误差控制阶段
① 预估计积分误差并与误差精度比较，如积分误差过大，则舍弃此步。 ② 计算优化的步长h 和阶数n 。
③ 如达到仿真结束时间，则停止，否则t =t +∆t ，重新进入第一步。 . .
4.1.2 坐标缩减的微分方程求解过程步骤
（1）坐标分离
将坐标的约束方程进行矩阵的满秩分解，可将系统的广义坐标列阵{q }分解成独立坐标列阵q i 和非独立坐标列阵q d ，即{q }=q i
（3）预估
用Adams-Bashforth 显示公式，根据独立坐标前几个小时步长的值，预估t n +1时刻的独立坐标值q i
（4）校正
用Adams-Moulton 隐示公式对上面的预估值，根据给定的收敛误差限进行校正，以得到独立坐标的校正值q i
（5）确定相关坐标
确定独立坐标的校正值后，可由相应公式计算出非独立坐标和其他系统状态变量值。
（6）积分误差控制
与上面预估—校正算法积分误差控制过程相同，如果预估值与校正的差值小{}{}{q d }。 T {}，P 表示预估值。 p C {}，C 表示校正值。
于给定的积分误差限，接受该解，进行下一时刻的求解。否则减小积分步长，重复第二步开始的预估步骤。
采用刚性积分算法Gstiff 对该发射装置动力学模型进行求解，Gstiff 求解器为刚性稳定算法，采用多步、变阶、变步长、固定系数算法，可直接求解DAE （微分—代数方程）。其特点是计算速度快，位移精度高。积分格式采用SI2格式，该积分格式考虑了速度约束方程，对速度和加速度的解比较精确，而且在步长很小时仍能保持计算的稳定。