RESRARCH ON THE CONTACT DYNAMIC ANALYSIS METHOD OF VEHICLE ON NONSMOOTH ROAD SURFACE
-
摘要: 车辆接触动力学研究有助于更加全面地了解整车力学性能, 对整车稳定性和座舱舒适性设计具有指导性作用. 文章基于非光滑多体系统理论, 对车辆系统展开接触动力学研究. 考虑到车辆在行驶过程中, 轮胎是与地面接触的唯一部件, 轮胎动力学模型对整车动力学分析发挥着重要的作用. 首先, 借助于张拉整体结构的思想, 建立轮胎的等效模型. 然后, 采用拉格朗日乘子方法, 将约束方程引入到拉格朗日函数中, 并基于第二类拉格朗日方程, 推导得到整车系统的一般动力学方程. 由于接触碰撞的存在, 会导致车辆的状态变量呈现出非连续性. 利用牛顿碰撞定律, 建立了碰撞响应的非光滑数学列式, 并利用光滑化的Fischer-Burmeister函数将非光滑数学列式进行等效替换, 以提高计算效率. 其次, 将摩擦响应的求解转化为优化问题进行描述, 以避免求解摩擦方向时出现奇异. 另外, 将车辆系统的一般动力学模型分解为光滑部分与非光滑部分进行求解, 并利用广义α离散策略给出了详细的求解流程. 最后, 对整车系统展开了多工况的接触碰撞数值仿真, 分析了不同工况下的接触响应对整车动力学性能的影响, 验证了所提方法的有效性.Abstract: The study of vehicle contact dynamics contributes to a more comprehensive understanding of the mechanical performance of the entire vehicle, and has a guiding role in the design of vehicle stability and cabin comfort. Therefore, this article is based on the theory of non smooth multi body system to conduct contact dynamics research on vehicle systems. Taking into account that the tire is the only component in contact with the ground during the driving process of vehicle, the tire dynamics model plays a crucial role in the overall vehicle dynamics analysis. First of all, this article establishes a multi body dynamic model of tires based on absolute coordinate modeling method. Then, by means of the Lagrange multiplier method, the constraint equation is introduced into the Lagrange function, and based on the second type of Lagrange equation, the general dynamic equation of the entire vehicle system is derived. Due to the presence of contact collisions, the state variables of the vehicle may exhibit discontinuity. In this paper, the Newton impact law is used to establish a non smooth mathematical formulation for impact response, and the smoothed Fischer-Burmeister function is used to equivalently replace the non smooth mathematical formulation to improve computational efficiency. Secondly, this article transforms the solution of friction response into an optimization problem for description, in order to avoid singularity when solving the friction direction. In addition, the general dynamic model of the vehicle system is decomposed into smooth and non smooth parts, and a detailed solution process is given based on generalized-α method. In the end, numerical simulations of contact collisions under multiple operating conditions are conducted on the entire vehicle system, and the impact of contact response under different operating conditions on the vehicle's dynamic performance is analyzed, and verified the effectiveness of the method proposed in this article.
-
Keywords:
- vehicle multibody system dynamics /
- tire modeling /
- impact /
- friction /
- nonsmooth contact
-
引 言
车辆在非光滑路面行驶过程中, 接触碰撞现象不可避免, 而接触碰撞响应会直接影响到车辆的力学性能, 对车辆展开接触动力学研究, 能够对车辆稳定性和座舱舒适性等方面设计提供指导作用. 然而, 由于接触碰撞问题具有高度的非线性, 对接触问题的求解一直以来是一个公认的难点. 同时, 由于车辆系统本身的动力学模型具有很强的非线性, 两个强非线性问题的相互耦合, 使得车辆接触动力学分析面临极大的挑战.
近年来, 随着多体系统动力学的发展, 诸多工程问题都能够利用多体动力学理论进行有效的解决[1-5]. 其中车辆是由多个部件通过运动副连接形成的机构, 其运动行为是由各个部件之间相互耦合实现, 是典型的多体系统. 因此, 利用多体系统理论对车辆系统展开动力学分析是一个有效的途径. 为了分析履带车辆传动系统扭矩在不平整路面的随机激励, 李铮等[6]结合扭转随机激励产生的机理, 推导了相关的理论计算公式, 并通过数值结果得出了路面不平度和车速是影响扭转随机激励的主要因素. 刘辉等[7]利用ADAMS软件建立了车辆系统的多自由度非线性模型, 并提出利用组合优化方法对控制系统控制参数优化分析, 通过数值算例表明了所提出的策略能够更加有效地提供准确的动力响应. 张云清等[8]提出了一种基于悬架摆臂移除的车辆多体系统动力学半递推建模技术, 推导了车辆闭环多体系统的控制方程, 并运用4阶龙格库塔方法与ABM方法对车辆动力学方程进行求解, 通过数值算例表明了所提出建模方法的高效性.
然而, 车辆在崎岖路面上行驶时, 通常会受到地面的接触碰撞响应. 在接触碰撞响应的作用下, 车辆的动力学行为会出现跳跃现象, 呈现出非光滑特性, 给车辆系统动力学方程的求解带来困难[9-23]. 同时, 地面接触力通过轮胎传递到车身时, 会对车身的稳定性带来影响. 范新秀等[24]考虑到车轮与地面的接触, 通过将车身与车轮进行简化处理, 建立了车轮与地面接触的线性互补方程, 并基于事件驱动方法对车辆纵向多体系统进行了仿真. Liu等[25]和Sun等[26]对轨道车辆系统中轮胎接触问题, 建立了轮胎-轨道接触模型, 并对轮胎-轨道接触的复杂动力学特性进行了数值研究. Peng等[27]考虑到轮胎与地面接触的摩擦系数难以估计的问题, 提出了一种新的非线性观测器, 用于估计车辆速度和轮胎-道路之间的摩擦系数, 并在不同轮胎转角下对轮胎与地面之间摩擦系数进行估计, 结果表明了所提方法的可靠性. 但是针对整车系统考虑到轮胎与地面的接触力以及接触力传递到舱体所产生的动力学特性分析较少.
为了分析轮胎与地面之间接触碰撞力对整车系统动力学特性, 给整车的设计提供必要的指导, 本文首先基于绝对坐标建模方法, 对轮胎以及舱体部分建立动力学模型, 相比于基于有限元方法的轮胎模型, 该轮胎模型既能反映出轮胎的物理特性, 又能够在实际情况中对参数进行快速迭代优化; 其次, 基于非光滑多体系统理论, 推导了车辆系统的一般非线性动力学方程, 并利用牛顿碰撞定律与库伦摩擦模型, 将轮胎与地面的接触关系表述为单边不等式约束条件, 并通过拉格朗日函数乘子法将约束条件嵌入到动力学方程中; 最后, 通过数值算例, 对轮胎模型以及整车系统进行了接触动力学仿真, 以期为车辆接触动力学分析设计提供数据支撑.
1. 轮胎建模
1.1 轮胎等效模型
轮胎是车辆与地面接触的唯一部件, 其关系到接触力从地面向舱体模块的传递, 对整车的接触动力学分析至关重要. 因此, 本部分首先对车辆系统的轮胎模型建模进行描述.
如图1(a)所示, 当轮胎行驶过程中受到挤压时, 轮胎的变形主要发生在两个方向上, 即纵向和横向的挤压变形, 两个的方向变形影响着整车的动力学分析. 考虑到轮胎的这种变形特性, 本文基于张拉整体结构的思想, 提出了一种新型的轮胎建模方法, 如图1(b)所示.
图1(b)中所提出的轮胎模型主要通过杆件与弹簧构成, 通过杆件与弹簧变形特点模拟轮胎受到挤压时的变形特性. 轮胎中的轮毂部分可通过图1(b)模型中的中间杆件组成, 胎面模型由外层的杆件围绕形成的圆面模拟, 内部气压可通过弹簧进行模拟. 由于弹簧与杆件能够承受挤压作用, 可在轴向发生变形. 因此, 当轮胎受到挤压时, 所提出轮胎模型可在纵向和横向发生变形, 在一定程度上能够等效真实轮胎的变形特性.
1.2 单元运动描述
对于上述小节中的轮胎模型, 本文采用考虑轴向变形的杆单元对杆件和弹簧进行描述, 如图2所示. 其中Ni和Ni−1是杆单元两端节点, 节点Ni和Ni−1在全局坐标系下对应的矢量分别为Ri和Ri−1, O-x-y-z代表全局坐标系. 根据位置有限元理论框架, 杆单元的广义坐标可以直接选用两端节点的广义位置坐标, 记为qr, 可以表示为
$$ {{\boldsymbol{q}}_r} = {\left[ {{\boldsymbol{R}}_{i - 1}^{\rm{T}},{\boldsymbol{R}}_i^{\rm{T}}} \right]^{\rm{T}}} = {\left[ {{x_{i - 1}},{y_{i - 1}},{z_{i - 1}},{x_i},{y_i},{z_i}} \right]^{\rm{T}}} $$ (1) 式中, xi−1, yi−1, zi−1, xi, yi和zi是杆单元两端节点在全局坐标系下的位置坐标.
根据图2所示, 杆单元的轴向矢量可以利用两端节点在全局坐标下的矢量Ri和Ri−1进行表示, 具有如下的形式
$$ \begin{split} & {{\boldsymbol{l}}_r} = {{{\boldsymbol{N}}_{i - 1}}{{\boldsymbol{N}}_i}} = {{\boldsymbol{R}}_i} - {{\boldsymbol{R}}_{i - 1}}= \\ & \qquad {\left[ {{x_i} - {x_{i - 1}},{y_i} - {y_{i - 1}},{z_i} - {z_{i - 1}}} \right]^{\rm{T}}} \end{split} $$ (2) 式中 lr为杆单元的轴向矢量.
根据式(2), 节点Ni和Ni−1中间任意一点在全局坐标系的矢量可以表示为Rξ, 且Rξ具有如下的表述形式
$$ {{\boldsymbol{R}}_\xi } = {{\boldsymbol{R}}_i} + \xi {{\boldsymbol{l}}_r} = \xi {{\boldsymbol{R}}_{i - 1}} + \left( {1 - \xi } \right){{\boldsymbol{R}}_i} = {{\boldsymbol{C}}_\xi }{{\boldsymbol{q}}_r} $$ (3) 式中, ξ$ \in $[0,1], Cξ具有如下表述形式
$$ {{\boldsymbol{C}}_\xi } = \left[ {\begin{array}{*{20}{l}} \xi &0&0&{1 - \xi }&0&0 \\ 0&\xi &0&0&{1 - \xi }&0 \\ 0&0&\xi &0&0&{1 - \xi } \end{array}} \right] $$ (4) 1.3 单元质量阵与外力列向量推导
根据式(3)中表述, 单元内部任一点的广义速度和加速度矢量可以通过式(3)分别对时间的一阶、二阶导数计算得到, 如下所示
$$ {{\boldsymbol{\dot R}}_\xi } = {{\boldsymbol{C}}_\xi }{{\boldsymbol{\dot q}}_r},{\text{ }}{{\boldsymbol{\ddot R}}_\xi } = {{\boldsymbol{C}}_\xi }{{\boldsymbol{\ddot q}}_r} $$ (5) 根据式(2), 单元的长度为
$$ {l_r} = ||{{\boldsymbol{l}}_r}|{|_2} {\text{ }} = \sqrt {{{\left( {{x_i} - {x_{i - 1}}} \right)}^2} + {{\left( {{y_i} - {y_{i - 1}}} \right)}^2} + {{\left( {{z_i} - {z_{i - 1}}} \right)}^2}} $$ (6) 假设单元的质量为mr, 单元密度沿着轴向方向线性分布, ρr = mr /lr, 单元的动能可以表示为
$$\begin{split} & {E_r} = \int_0^{{l_r}} {\frac{1}{2}{\rho _r}{\boldsymbol{\dot R}}_\xi ^{\mathrm{T}}} {{{\boldsymbol{\dot R}}}_\xi }{\mathrm{d}}l = \int_0^1 {\frac{{{m_r}}}{2}} {\boldsymbol{\dot R}}_\xi ^{\mathrm{T}}{{{\boldsymbol{\dot R}}}_\xi }{\mathrm{d}}\xi = \\ &\qquad \frac{1}{2}{\boldsymbol{\dot q}}_r^{\rm{T}}\left( {{m_r}\int_0^1 {{\boldsymbol{C}}_\xi ^{\rm{T}}{{\boldsymbol{C}}_\xi }{\mathrm{d}}\xi } } \right){{{\boldsymbol{\dot q}}}_r} \end{split} $$ (7) 根据式(7), 单元质量阵为
$$ \begin{split} & {{\boldsymbol{M}}_r} = {m_r}\int_0^1 {{\boldsymbol{C}}_\xi ^{\rm{T}}{{\boldsymbol{C}}_\xi }} {\mathrm{d}}\xi= \\ &\qquad {m_r}\left\{ {\int_0^1 {\left[ {\begin{array}{*{20}{c}} {{\xi ^2}}&{\xi \left( {1 - \xi } \right)} \\ {\xi \left( {1 - \xi } \right)}&{{\xi ^2}} \end{array}} \right]{\mathrm{d}}\xi } } \right\} \otimes {{\boldsymbol{I}}_3}= \\ &\qquad \frac{{{m_r}}}{6}\left[ {\begin{array}{*{20}{c}} {2{{\boldsymbol{I}}_3}}&{{{\boldsymbol{I}}_3}} \\ {{{\boldsymbol{I}}_3}}&{2{{\boldsymbol{I}}_3}} \end{array}} \right]\end{split}$$ (8) 基于式(6), 然后利用线弹性理论, 单元轴向内力可以通过轴向的变形计算得到, 计算公式如下
$$ {f_r} = {E_r}{\varepsilon _r}{A_r} $$ (9) 根据虚功原理, 内力所作虚功可以通过以下公式表示
$$ \begin{split} & \delta {W_r} = - {f_r}{\boldsymbol{\hat l}}_r^{\mathrm{T}}\delta {{\boldsymbol{R}}_{i - 1}} + {f_r}{\boldsymbol{\hat l}}_r^{\mathrm{T}}\delta {{\boldsymbol{R}}_i}= \\ &\qquad \left[ { - {f_r}{\boldsymbol{\hat l}}_r^{\mathrm{T}},{f_r}{\boldsymbol{\hat l}}_r^{\mathrm{T}}} \right]\delta {{\boldsymbol{q}}_r} = {{\boldsymbol{f}}_r}\delta {{\boldsymbol{q}}_r}\end{split} $$ (10) 式中, $ {\hat{\boldsymbol{l}}}_{r}^{\rm{T}} $ = $ {\boldsymbol{l}}_{r}^{\rm{T}}/l_r $. 根据式(10), 可得到单元广义力向量为 fr = [−fr$ {\hat{\boldsymbol{l}}}_{r}^{\rm{T}} $,fr$ {\hat{\boldsymbol{l}}}_{r}^{\rm{T}} $][28].
2. 整车非光滑多体动力学模型
基于多体动力学理论可知, 当车辆系统无接触作用时, 动力学方程可以描述为以下微分代数方程
$$ \left. \begin{gathered} {\boldsymbol{M}}\left( {\boldsymbol{q}} \right){\boldsymbol{\dot v}} + {\boldsymbol{\varPhi }}\left( {{\boldsymbol{q}},t} \right)_{\boldsymbol{q}}^{\rm{T}}{\boldsymbol{\lambda }} - {\boldsymbol{F}}\left( {{\boldsymbol{q}},{\boldsymbol{v}},t} \right) = {\boldsymbol{0}} \\ {\boldsymbol{\varPhi }}\left( {{\boldsymbol{q}},t} \right) = {\boldsymbol{0}} \\ \end{gathered} \right\} $$ (11) 式中, q和v为整车的广义位置和广义速度, M和F分别是系统的总质量阵和合力列向量(外力、内力和惯性力), 总质量阵和合力可以通过第1节中的单元力组装形成; ${\boldsymbol{\dot v}}$是系统的广义加速度; Φ代表双边约束; λ是约束方程的拉格朗日乘子.
当车辆在行驶过程中, 需要考虑轮胎与地面之间的接触力, 如图3(a)所示. 本文假设轮胎与地面接触过程中没有穿透的现象发生, 并用距离约束函数g (非负)表示轮胎与地面之间是否发生接触, 如图3(b)所示. 当g > 0时, 轮胎与地面之间无接触; 当g = 0时, 地面与轮胎之间发生接触, 如图3(c)所示. 通过拉格朗日乘子法, 将单边距离函数约束g引入到车辆的非线性动力学方程中, 并根据非光滑多体系统理论[29-30], 考虑接触的车辆系统动力学方程可以表示为
$$ \left.\begin{split} &{\boldsymbol{M}}\left({\boldsymbol{q}}\right){\mathrm{d}}{\boldsymbol{v}} + {\boldsymbol{\varPhi}} {\left({\boldsymbol{q}},t\right)}_{q}^{{\mathrm{T}}}{\mathrm{d}}{{\boldsymbol{i}}}_{c}-{\boldsymbol{F}}\left({\boldsymbol{q}},{\boldsymbol{v}},t\right){\mathrm{d}}t-\\ &\qquad {\boldsymbol{g}}{\left({\boldsymbol{q}}\right)}_{q}^{{\mathrm{T}}}{{\boldsymbol{\varLambda}} }_{n}-{{\boldsymbol{D}}}_{t}{{\boldsymbol{\varLambda}} }_{t} = {\boldsymbol{0}}\\ &{\boldsymbol{\varPhi}} \left({\boldsymbol{q}},t\right) = {\boldsymbol{0}}\\ &{\boldsymbol{0}}\leqslant {\boldsymbol{g}}{\left({\boldsymbol{q}}\right)}^{j}\perp {{\boldsymbol{\varLambda}} }_{n}^{j}\geqslant {\boldsymbol{0}},\text{ }j\in {N}_{A}\\ &\left|\right|{{\boldsymbol{\varLambda}} }_{t}^{j}\left|\right|\leqslant \mu {{\boldsymbol{\varLambda}} }_{n}^{j},\text{ }j\in {N}_{A}\\ &\left|\right|{{\boldsymbol{v}}}_{t}^{j}\left|\right|\cdot \left|\right|\mu {{\boldsymbol{\varLambda}} }_{n}^{j}-{{\boldsymbol{\varLambda}} }_{t}^{j}\left|\right|,\text{ }j\in {N}_{A}\end{split}\right\} $$ (12) 式中, g(q)代表单边约束; dic代表在接触作用下双边约束拉格朗日乘子的冲量; $ {\boldsymbol{v}}_{t}^{j} $代表第j个接触的切向速度; μ为动摩擦系数; Λn和Λt代表拉格朗日乘子和接触作用冲量; Dt是由切向接触方向组成的矩阵; ${\boldsymbol{g}}{\left({\boldsymbol{q}}\right)}_{q}^{{\mathrm{T}}} $表示单边约束的雅克比矩阵的转置; dt表示时间步长; 第3项表示距离函数与法向接触响应之间的互补数学列式; 第4项为由库伦摩擦表示的切向与法向之间的关系; NA表示轮胎与地面之间的接触点数量; dv表示在接触响应作用下系统速度的突变. 且方程(12)中, dic和dv具有以下表述形式
$$ \left. \begin{gathered} {\mathrm{d}}{\boldsymbol{v}} = {\boldsymbol{\dot v}}{\mathrm{d}}t + \sum\limits_i {\left[ {{\boldsymbol{v}}\left( {{t_i}} \right) - {{\boldsymbol{v}}^ - }\left( {{t_i}} \right)} \right]} \delta {t_i} \\ {\mathrm{d}}{{\boldsymbol{i}}_c} = {\boldsymbol{\lambda }}{\mathrm{d}}t + \sum\limits_i {{{\boldsymbol{\varLambda }}_b}\delta {t_i}} \\ \end{gathered} \right\} $$ (13) 式中, v(ti)代表车辆系统接触碰撞后的速度, v−(ti)代表车辆系统碰撞前的速度, Λb代表接触碰撞造成的双边约束拉格朗日乘子变化.
考虑式(13), 并将式(12)减去式(11), 可得到接触力作用下车辆系统的动力学方程, 如下式所示
$$ \left.\begin{split} &{\boldsymbol{M}}\left({\boldsymbol{q}}\right){\boldsymbol{w}} + {\boldsymbol{\varPhi}} {\left({\boldsymbol{q}},t\right)}_{q}^{{\mathrm{T}}}{{\boldsymbol{\varLambda}} }_{b}-{\boldsymbol{g}}{\left({\boldsymbol{q}}\right)}_{q}^{{\mathrm{T}}}{{\boldsymbol{\varLambda}} }_{n}-{{\boldsymbol{D}}}_{t}{{\boldsymbol{\varLambda}} }_{t} = {\boldsymbol{0}}\\ &{\boldsymbol{\varPhi}} \left({\boldsymbol{q}},t\right) = {\boldsymbol{0}}\\ &{\boldsymbol{0}}\leqslant {\boldsymbol{g}}{\left({\boldsymbol{q}}\right)}^{j}\perp {{\boldsymbol{\varLambda}} }_{n}^{j}\geqslant {\boldsymbol{0}},\text{ }j\in {N}_{A}\\ &\left|\right|{{\boldsymbol{\varLambda}} }_{t}^{j}\left|\right|\leqslant \mu {{\boldsymbol{\varLambda}} }_{n}^{j},\text{ }j\in {N}_{A}\\ &\left|\right|{{\boldsymbol{v}}}_{t}^{j}\left|\right|\cdot \left|\right|\mu {{\boldsymbol{\varLambda}} }_{n}^{j}-{{\boldsymbol{\varLambda}} }_{t}^{j}\left|\right|,\text{ }j\in {N}_{A}\end{split} \right\}$$ (14) 式中, w = v(ti) – v−(ti)代表接触碰撞前后速度跳跃. 根据式(15), 速度跳跃w是由法向接触响应和切向接触响应共同产生. 因此, 根据矢量分解过程, 可将速度跳跃在接触法向和切向进行矢量分解, 如图4所示.
根据图4中速度跳跃量的分解形式, 将速度跳跃分别在接触法向与切向进行分解, 则式(14)可以分解为下式所示形式
$$ \left.\begin{split} &\left\{\begin{aligned} &{\boldsymbol{M}}\left({\boldsymbol{q}}\right){{\boldsymbol{w}}}_{1} + {\boldsymbol{\varPhi}} {\left({\boldsymbol{q}},t\right)}_{q}^{{\mathrm{T}}}{{\boldsymbol{\varLambda}} }_{{b}_{1}}-{\boldsymbol{g}}{\left({\boldsymbol{q}}\right)}_{q}^{{\mathrm{T}}}{{\boldsymbol{\varLambda}} }_{n} = {\boldsymbol{0}}\\ &{\boldsymbol{\varPhi}} \left({\boldsymbol{q}},t\right) = {\boldsymbol{0}}\\ &{\boldsymbol{0}}\leqslant {\boldsymbol{g}}{\left({\boldsymbol{q}}\right)}^{j}\perp {{\boldsymbol{\varLambda}} }_{n}^{j}\geqslant {\boldsymbol{0}},\text{ }j\in {N}_{A}\end{aligned}\right.\\ &\left\{\begin{aligned} &\left|\right|{{\boldsymbol{\varLambda}} }_{t}^{j}\left|\right|\leqslant \mu {{\boldsymbol{\varLambda}} }_{n}^{j},\text{ }j\in {N}_{A}\\ &\left|\right|{{\boldsymbol{v}}}_{t}^{j}\left|\right|\cdot \left|\right|\mu {{\boldsymbol{\varLambda}} }_{n}^{j}-{{\boldsymbol{\varLambda}} }_{t}^{j}\left|\right|,\text{ }j\in {N}_{A}\\ &{\boldsymbol{M}}\left({\boldsymbol{q}}\right){{\boldsymbol{w}}}_{2} + {\boldsymbol{\varPhi}} {\left({\boldsymbol{q}},t\right)}_{q}^{{\mathrm{T}}}{{\boldsymbol{\varLambda}} }_{{b}_{2}}-{{\boldsymbol{D}}}_{t}{{\boldsymbol{\varLambda}} }_{t} = {\boldsymbol{0}}\\ &{\boldsymbol{\varPhi}} \left({\boldsymbol{q}},t\right) = {\boldsymbol{0}}\end{aligned}\right.\end{split}\right\} $$ (15) 式中, ${{\boldsymbol{\varLambda}} }_{b_1} $和${{\boldsymbol{\varLambda}} }_{b_2} $分别是在法向接触响应以及切向接触响应作用下的拉格朗日乘子, 且有Λb = ${\boldsymbol{\varLambda}}_{b_1}$ + ${{\boldsymbol{\varLambda}}}_{b_2}$.根据式(15)可知, 法向接触不等式约束在位置层面表述, 本文利用牛顿碰撞定律, 进一步在速度层面对法向接触进行描述. 牛顿碰撞定律定义为接触部位碰撞前后的相对速度关系, 因此下式第3项可以用下述数学列式进行表示
$$ {\boldsymbol{0}} \leqslant {\boldsymbol{g}}\left( {\boldsymbol{q}} \right)_{\boldsymbol{q}}^j{\boldsymbol{v}} + {e^j}{\boldsymbol{g}}\left( {\boldsymbol{q}} \right)_{\boldsymbol{q}}^j{{\boldsymbol{v}}^ - } \bot {\boldsymbol{\varLambda }}_n^j \geqslant {\boldsymbol{0}},{\text{ }}j \in {N_A} $$ (16) 式中, e$ \in $[0,1]代表碰撞系数. 当e = 0时, 表示完全塑性碰撞; 当e = 1时, 代表完全弹性碰撞.
考虑到式(11)、式(15)和式(16), 整车接触碰撞动力学方程的一般形式如下
$$ \left.\begin{aligned} &\boldsymbol{M}(\tilde{\boldsymbol{q}}) \dot{\tilde{\boldsymbol{v}}}+\boldsymbol{\varPhi}(\tilde{\boldsymbol{q}}, t)_q^{\mathrm{T}} \tilde{{\boldsymbol{\lambda}}}-\boldsymbol{F}(\tilde{\boldsymbol{q}}, \tilde{\boldsymbol{v}}, t)=\boldsymbol{0} \\ &\boldsymbol{\varPhi}(\tilde{\boldsymbol{q}}, t)=\boldsymbol{0}\end{aligned}\right\}\tag{17a} $$ $$ \left.\begin{aligned} &\boldsymbol{M}({\overset{\frown} {\boldsymbol{q}}}) \boldsymbol{w}_1+\boldsymbol{\varPhi}({\overset{\frown} {\boldsymbol{q}}}, t)_q^{\rm{T}} \boldsymbol{\varLambda}_{b_1}-\boldsymbol{g}({\overset{\frown} {\boldsymbol{q}}})_q^{\rm{T}} \boldsymbol{\varLambda}_n=\boldsymbol{0} \\ &\boldsymbol{\varPhi}({\overset{\frown} {\boldsymbol{q}}}, t)=\boldsymbol{0} \\ &{\boldsymbol{0}} \leqslant \boldsymbol{g}({\overset{\frown} {\boldsymbol{q}}})_{\hat{q}}^j \hat{\boldsymbol{v}}+e^j \boldsymbol{g}({\overset{\frown} {\boldsymbol{q}}})_q^j \hat{\boldsymbol{v}}^{-} \perp \boldsymbol{\varLambda}_n^j \geqslant {\boldsymbol{0}}, j \in N_A\end{aligned}\right\}\tag{17b}$$ $$ \left.\begin{aligned} &\left(\boldsymbol{\gamma}_1^j, \boldsymbol{\gamma}_2^j\right)=\underset{\sqrt{\left(\gamma_1^j\right)^2+\left(\gamma_1^j\right)^2} \leqslant \mu \boldsymbol{\varLambda}_n^j}{\arg \min } \hat{\boldsymbol{v}}^{\rm{T}}\left(\boldsymbol{D}_{u_1}^j {\boldsymbol{\gamma}}_1^j+\boldsymbol{D}_{u_2}^j {\boldsymbol{\gamma}}_2^j\right),\\ &\qquad j \in N_A \\ &\boldsymbol{M}(\hat{\boldsymbol{q}}) \boldsymbol{w}_2+\boldsymbol{\varPhi}(\hat{\boldsymbol{q}}, t)_q^{\rm{T}} \boldsymbol{\varLambda}_{b_2}-\boldsymbol{D}_{u_1} {\boldsymbol{\gamma}}_1-\boldsymbol{D}_{u_2} {\boldsymbol{\gamma}}_1=\boldsymbol{0} \\ &\boldsymbol{\varPhi}(\hat{\boldsymbol{q}}, t)=\boldsymbol{0}\end{aligned}\right\} \tag{17c}$$ 式(17c)中, 摩擦问题的求解转化为优化问题, “$ \tilde{\cdot} $”代表在光滑力作用下车辆的状态变量; “$ {\overset{\frown} \cdot} $”代表碰撞接触响应作用下车辆系统的状态变量; “$ \hat{\cdot} $”代表摩擦接触响应作用下车辆系统的状态变量. 且各部分状态变量之间具有如下的关系[29-30]
$$ \left.\begin{split} &{\overset{\frown} {\boldsymbol{q}}} = \tilde{{\boldsymbol{q}}} + \frac{1}{2}{\mathrm{d}}t\cdot{{\boldsymbol{w}}}_{1},\quad \hat{{\boldsymbol{q}}} = \tilde{{\boldsymbol{q}}} + \frac{1}{2}{\mathrm{d}}t\cdot\left({{\boldsymbol{w}}}_{1} + {{\boldsymbol{w}}}_{2}\right)\\ &{\overset{\frown} {\boldsymbol{v}}} = \tilde{{\boldsymbol{v}}} + {{\boldsymbol{w}}}_{1},\quad \text{ }\hat{{\boldsymbol{v}}} = \tilde{{\boldsymbol{v}}} + {{\boldsymbol{w}}}_{1} + {{\boldsymbol{w}}}_{2}\end{split}\right\}$$ (18) 3. 动力学分析流程
式(17)给出了整车系统的接触动力学模型, 对车辆进行接触动力学分析时, 需要对此方程进行离散求解. 本部分基于广义α离散策略, 给出具体的动力学分析流程. 广义α离散格式如下所示
$$ \left. \begin{gathered} {{{\boldsymbol{\tilde q}}}_{n + 1}} = {{\boldsymbol{q}}_n} + h{{\boldsymbol{v}}_n} + {h^2}\left( {0.5 - \beta } \right){{\boldsymbol{a}}_n} + {h^2}\beta {{\boldsymbol{a}}_{n + 1}} \\ {{{\boldsymbol{\tilde v}}}_{n + 1}} = {{\boldsymbol{v}}_n} + h\left( {1 - \gamma } \right){{\boldsymbol{a}}_n} + h\gamma {{\boldsymbol{a}}_{n + 1}} \\ \end{gathered} \right\} $$ (19) 式中, 中间变量${\boldsymbol{a}}$满足以下形式
$$ \left. \begin{gathered} \left( {1 - {\alpha _m}} \right){{\boldsymbol{a}}_{n + 1}} + {\alpha _m}{{\boldsymbol{a}}_n} = \left( {1 - {\alpha _f}} \right){{{\boldsymbol{\dot {\tilde v}}}}_{n + 1}} + {\alpha _f}{{{\boldsymbol{\dot {\tilde v}}}}_n} \\ {{\boldsymbol{a}}_0} = {{{\boldsymbol{\dot {\tilde v}}}}_0} \\ \end{gathered} \right\} $$ (20) 式中, h为时间步长, γ, β, αm和αf为算法参数, 其取值可以通过下式求得
$$\left. \begin{split} & {\alpha _m} = \frac{{2{\rho _\infty } - 1}}{{{\rho _\infty } + 1}},\quad {\alpha _f} = \frac{{{\rho _\infty }}}{{{\rho _\infty } + 1}} \\ & \beta = \frac{1}{4}{\left( {1 + {\alpha _f} - {\alpha _m}} \right)^2},\quad \gamma = \frac{1}{2} + {\alpha _f} - {\alpha _m} \end{split} \right\}$$ (21) 式中, 谱半径$ {\rho _\infty } $取值范围在0 ~ 1之间.
利用上述广义α离散格式, 可将车辆系统非线性动力学方程写为下述离散形式
$$\left. \begin{split} & {{\boldsymbol{r}}_{1}}=\boldsymbol{M}\left( {{{\boldsymbol{\tilde{q}}}}_{n+1}} \right){{{{\dot{\tilde{{\boldsymbol{v}}}}}}}_{n+1}}+\boldsymbol{\varPhi }\left( {{{\boldsymbol{\tilde{q}}}}_{n+1}},{{t}_{n+1}} \right)_{\boldsymbol{q}}^{\rm{T}}{{{\boldsymbol{\tilde{\lambda }}}}_{n+1}}- \\ &\qquad \boldsymbol{F}\left( {{{{\tilde{{\boldsymbol{q}}}}}}_{n+1}},{{{{\tilde{{\boldsymbol{v}}}}}}_{n+1}},{{t}_{n+1}} \right)=\boldsymbol{0} \\ & {{\boldsymbol{r}}_{2}}=\boldsymbol{\varPhi }\left( {{{{\tilde{{\boldsymbol{q}}}}}}_{n+1}},{{t}_{n+1}} \right)=\boldsymbol{0} \\ & {{\boldsymbol{r}}_{3}}=\boldsymbol{M}\left( {{{\overset{\frown} {\boldsymbol{q}}}}_{n+1}} \right){{\boldsymbol{w}}_{1,n+1}}+\boldsymbol{\varPhi }\left( {{{\overset{\frown} {\boldsymbol{q}}}}_{n+1}},{{t}_{n+1}} \right)_{\boldsymbol{q}}^{\rm{T}}{{\boldsymbol{\varLambda }}_{{{b}_{1}},n+1}}- \\ &\qquad \boldsymbol{g}\left( {{{\overset{\frown} {\boldsymbol{q}}}}_{n+1}} \right)_{\boldsymbol{q}}^{\rm{T}}{{\boldsymbol{\varLambda }}_{n,n+1}}=\boldsymbol{0} \\ & {{\boldsymbol{r}}_{4}}=\boldsymbol{\varPhi }\left( {{{\overset{\frown} {\boldsymbol{q}}}}_{n+1}},{{t}_{n+1}} \right)=\boldsymbol{0} \\ & {{\boldsymbol{r}}_{5}}=0\leqslant \boldsymbol{g}\left( {{{\overset{\frown} {\boldsymbol{q}}}}_{n+1}} \right)_{{\overset{\frown} {\boldsymbol{q}}}}^{j}{{{\overset{\frown} {\boldsymbol{v}}}}_{n+1}}+ \\ &\qquad {{e}^{j}}\boldsymbol{g}\left( {{{\overset{\frown} {\boldsymbol{q}}}}_{n+1}} \right)_{\boldsymbol{q}}^{j}{{\overset{\frown} {\boldsymbol{v}}}_{n}}\bot \;\boldsymbol{\varLambda }_{n.n+1}^{j}\geqslant 0,\text{ }j\in {{N}_{A}} \\ & {{\boldsymbol{r}}_{6}}=\left( \boldsymbol{\gamma }_{1,n+1}^{j},\boldsymbol{\gamma }_{2,n+1}^{j} \right)= \underset{\sqrt{{{\left( \boldsymbol{\gamma }_{1,n+1}^{j} \right)}^{2}}+{{\left( \boldsymbol{\gamma }_{1,n+1}^{j} \right)}^{2}}}\leqslant \mu \boldsymbol{\varLambda }_{n,n+1}^{j}}{ \arg \min }\,\boldsymbol{\hat{v}}_{n+1}^{\rm{T}}\cdot\\ &\qquad \left( \boldsymbol{D}_{{{\boldsymbol{u}}_{1}}}^{j}\boldsymbol{\gamma }_{1,n+1}^{j}+\boldsymbol{D}_{{{\boldsymbol{u}}_{2}}}^{j}\boldsymbol{\gamma }_{2,n+1}^{j} \right), \quad j\in {{N}_{A}} \\ & {{\boldsymbol{r}}_{7}}=\boldsymbol{M}\left( {{{\boldsymbol{\hat{q}}}}_{n+1}} \right){{\boldsymbol{w}}_{2,n+1}}+\boldsymbol{\varPhi }\left( {{{\boldsymbol{\hat{q}}}}_{n+1}},{{t}_{n+1}} \right)_{\boldsymbol{q}}^{\rm{T}}{{\boldsymbol{\varLambda }}_{{{b}_{2}},n+1}}- \\ &\qquad {{\boldsymbol{D}}_{{{\boldsymbol{u}}_{1}}}}{{\boldsymbol{\gamma }}_{1,n+1}}-{{\boldsymbol{D}}_{{{\boldsymbol{u}}_{2}}}}{{\boldsymbol{\gamma }}_{1,n+1}}=\boldsymbol{0} \\ & {{\boldsymbol{r}}_{8}}=\boldsymbol{\varPhi }\left( {{{\boldsymbol{\hat{q}}}}_{n+1}},{{t}_{n+1}} \right)=\boldsymbol{0} \end{split}\right\} $$ (22) 本文利用Newton-Raphson方法对非线性动力学方程进行迭代求解. 首先, 对式(22)中的非线性方程r1和r2, 取$ \stackrel{ ~ }{\boldsymbol{q}} $n+1和$ \stackrel{ ~ }{\boldsymbol{\lambda }} $n+1作为独立的未知变量, 将非线性方程[r1;r2]对独立未知变量求雅克比矩阵, 可得如下形式
$$ {{\boldsymbol{J}}_1} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {{\boldsymbol{r}}_1}}}{{\partial {{{\boldsymbol{\tilde q}}}_{n + 1}}}}}&{\dfrac{{\partial {{\boldsymbol{r}}_1}}}{{\partial {{{\boldsymbol{\tilde \lambda }}}_{n + 1}}}}} \\ {\dfrac{{\partial {{\boldsymbol{r}}_2}}}{{\partial {{{\boldsymbol{\tilde q}}}_{n + 1}}}}}&{\dfrac{{\partial {{\boldsymbol{r}}_2}}}{{\partial {{{\boldsymbol{\tilde \lambda }}}_{n + 1}}}}} \end{array}} \right] $$ (23) 式中, 每一项具有如下表述形式
$$ \left. \begin{aligned} & \frac{\partial {{\boldsymbol{r}}_{1}}}{\partial {{{\boldsymbol{\tilde{q}}}}_{n+1}}}=\frac{\partial \left( \boldsymbol{M}{{{\boldsymbol{\dot{\tilde{v}}}}}_{n+1}} \right)}{\partial {{{\boldsymbol{\tilde{q}}}}_{n+1}}}+\boldsymbol{M}\frac{\partial {{{\boldsymbol{\dot{\tilde{v}}}}}_{n+1}}}{\partial {{{\boldsymbol{\tilde{q}}}}_{n+1}}}+ \\ &\qquad \frac{\partial \left( \boldsymbol{\varPhi }_{\boldsymbol{q}}^{\rm{T}}{{{\boldsymbol{\tilde{\lambda }}}}_{n+1}} \right)}{\partial {{{\boldsymbol{\tilde{q}}}}_{n+1}}}-\frac{\partial \boldsymbol{F}}{\partial {{{\boldsymbol{\tilde{q}}}}_{n+1}}}-\frac{\partial \boldsymbol{F}}{\partial {{{\boldsymbol{\tilde{v}}}}_{n+1}}}\frac{\partial {{{\boldsymbol{\tilde{v}}}}_{n+1}}}{\partial {{{\boldsymbol{\tilde{q}}}}_{n+1}}} \\ & \frac{\partial {{\boldsymbol{r}}_{1}}}{\partial {{{\boldsymbol{\tilde{\lambda }}}}_{n+1}}}=\boldsymbol{\varPhi }_{\boldsymbol{q}}^{\rm{T}},\frac{\partial {{\boldsymbol{r}}_{2}}}{\partial {{{\boldsymbol{\tilde{q}}}}_{n+1}}}={{\boldsymbol{\varPhi }}_{\boldsymbol{q}}},\frac{\partial {{\boldsymbol{r}}_{2}}}{\partial {{{\boldsymbol{\tilde{\lambda }}}}_{n+1}}}=\boldsymbol{0} \end{aligned} \right\} $$ (24) 式中, $ \partial {{\boldsymbol{\dot {\tilde v}}}_{n + 1}}/\partial {{\boldsymbol{\tilde q}}_{n + 1}} = {\boldsymbol{I}}\hat \gamma $, $ \partial {{\boldsymbol{\tilde v}}_{n + 1}}/\partial {{\boldsymbol{\tilde q}}_{n + 1}} = {\boldsymbol{I}}\hat \beta $, $ \hat \beta = 1 - {\alpha _m}/\left[{h^2}\beta \left( {1 - {\alpha _f}} \right)\right] $, $ \hat \gamma = \gamma /(h\beta) $.
对于式(22)中的非线性动力学方程[r3;r4;r5], 由于非线性方程r5为互补形式, 当方程维数较高时, 将其转化为标准的互补列式求解效率较低[31]. 因而, 本文利用光滑化的Fischer-Burmeister函数, 将其转化为非线性方程组进行求解. 非线性方程r5转化后有如下的等价形式
$$ {{\boldsymbol{r}}_6} = \sqrt {{{\boldsymbol{a}}^2} + {{\boldsymbol{b}}^2} + {\zeta ^2}} - \left( {{\boldsymbol{a}} + {\boldsymbol{b}}} \right) = {\boldsymbol{0}},{\text{ }}j \in {N_A} $$ (25) 其中$ {\boldsymbol{a}} = {\boldsymbol{g}}\left( {{{{{\overset{\frown} {\boldsymbol{q}} }}}_{n + 1}}} \right)_{{\overset{\frown} {\boldsymbol{q}}}}^j{{\overset{\frown} {\boldsymbol{v}}}_{n + 1}} + {e^j}{\boldsymbol{g}}\left( {{{{\overset{\frown} {\boldsymbol{q}}}}_{n + 1}}} \right)_{\boldsymbol{q}}^j{{\overset{\frown} {\boldsymbol{v}}}_n}, {\boldsymbol{b}} = {\boldsymbol{\varLambda }}_{n.n + 1}^j $.
将非线性动力学方程r5进行上述转化后, 对于非线性方程组[r3;r4;r5], 取w1,n + 1, Λb1,n + 1和Λn,n + 1为独立未知变量, 则[r3;r4;r5]对其雅克比矩阵为
$$ {{\boldsymbol{J}}_2} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {{\boldsymbol{r}}_3}}}{{\partial {{\boldsymbol{w}}_{1,n + 1}}}}}&{\dfrac{{\partial {{\boldsymbol{r}}_3}}}{{\partial {{\boldsymbol{\varLambda }}_{b1,n + 1}}}}}&{\dfrac{{\partial {{\boldsymbol{r}}_3}}}{{\partial {{\boldsymbol{\varLambda }}_{n,n + 1}}}}} \\ {\dfrac{{\partial {{\boldsymbol{r}}_4}}}{{\partial {{\boldsymbol{w}}_{1,n + 1}}}}}&{\dfrac{{\partial {{\boldsymbol{r}}_4}}}{{\partial {{\boldsymbol{\varLambda }}_{b1,n + 1}}}}}&{\dfrac{{\partial {{\boldsymbol{r}}_4}}}{{\partial {{\boldsymbol{\varLambda }}_{n,n + 1}}}}} \\ {\dfrac{{\partial {{\boldsymbol{r}}_5}}}{{\partial {{\boldsymbol{w}}_{1,n + 1}}}}}&{\dfrac{{\partial {{\boldsymbol{r}}_5}}}{{\partial {{\boldsymbol{\varLambda }}_{b1,n + 1}}}}}&{\dfrac{{\partial {{\boldsymbol{r}}_5}}}{{\partial {{\boldsymbol{\varLambda }}_{n,n + 1}}}}} \end{array}} \right] $$ (26) 针对式(22)中的第3部分非线性动力学方程, 首先将r6进行优化求解, 然后将计算得到的切向摩擦响应带入到非线性方程r7和r8中, 对非线性方程r7和r8进行求解. 对r6中的不等式约束, 通过引入拉格朗日乘子χ, 可以将其转化为如下形式
$$ \begin{split} & L\left( {{\boldsymbol{\gamma }}_{1,n + 1}^j,{\boldsymbol{\gamma }}_{2,n + 1}^j} \right) = {\boldsymbol{\hat v}}_{n + 1}^{\rm{T}}\left( {{\boldsymbol{D}}_{{{\boldsymbol{u}}_1}}^j{\boldsymbol{\gamma }}_{1,n + 1}^j + {\boldsymbol{D}}_{{{\boldsymbol{u}}_2}}^j{\boldsymbol{\gamma }}_{2,n + 1}^j} \right) + \\ &\qquad \chi \left[ {{{\left( {{\boldsymbol{\gamma }}_{1,n + 1}^j} \right)}^2} + {{\left( {{\boldsymbol{\gamma }}_{1,n + 1}^j} \right)}^2} - {{\left( {\mu {\boldsymbol{\varLambda }}_{n,n + 1}^j} \right)}^2}} \right] \end{split} $$ (27) 利用KKT条件, 通过对式(27)求解, 可计算得到拉格朗日乘子χ
$$ \chi = \pm \frac{{\sqrt {{\wp ^2} + {\aleph ^2}} }}{{2\Re }} $$ (28) 式中, $ \wp $, $ \aleph $和$ \Re $具有如下的表达式
$$ \wp = {\left( {{\boldsymbol{D}}_{{{\boldsymbol{u}}_1}}^{j,{{{\mathrm{T}}}}}{{{\boldsymbol{\hat v}}}_{n + 1}}} \right)^{\rm{T}}},\quad \aleph = {\left( {{\boldsymbol{D}}_{{{\boldsymbol{u}}_2}}^{j,{{{\mathrm{T}}}}}{{{\boldsymbol{\hat v}}}_{n + 1}}} \right)^{\rm{T}}},\quad \Re = \mu {\boldsymbol{\varLambda }}_{n,n + 1}^j $$ (29) 进而可求得$ {\boldsymbol{\gamma }}_{1,n + 1}^j $和$ {\boldsymbol{\gamma }}_{2,n + 1}^j $
$$ {\boldsymbol{\gamma }}_{1,n + 1}^j = - \frac{\wp }{{2\chi }},\quad {\boldsymbol{\gamma }}_{2,n + 1}^j = - \frac{\aleph }{{2\chi }} $$ (30) 通过将式(30)带入到r7中, 最小的目标值即真实的$ {\boldsymbol{\gamma }}_{1,n + 1}^j $和$ {\boldsymbol{\gamma }}_{2,n + 1}^j $. 将得到的$ {\boldsymbol{\gamma }}_{1,n + 1}^j $和$ {\boldsymbol{\gamma }}_{2,n + 1}^j $带入到非线性动力学方程[r7;r8]中, 并取w2,n+1和${{\boldsymbol{\varLambda}} }_{b_{2},n+1} $作为独立未知变量, 求得[r7;r8]关于独立未知变量的雅克比矩阵如下
$$ {{\boldsymbol{J}}_3} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {{\boldsymbol{r}}_8}}}{{\partial {{\boldsymbol{w}}_{2,n + 1}}}}}&{\dfrac{{\partial {{\boldsymbol{r}}_8}}}{{\partial {{\boldsymbol{\varLambda }}_{b_{2},n + 1}}}}} \\ {\dfrac{{\partial {{\boldsymbol{r}}_9}}}{{\partial {{\boldsymbol{w}}_{2,n + 1}}}}}&{\dfrac{{\partial {{\boldsymbol{r}}_9}}}{{\partial {{\boldsymbol{\varLambda }}_{b_{2},n + 1}}}}} \end{array}} \right] $$ (31) 利用非线性动力学方程[r1;r2],[r3;r4;r5]和[r7;r8]以及式(23)、式(26)和式(31), 进行如式(32)迭代求解, 即可实现对整车的接触动力学分析
$$ \varDelta = - \frac{{\boldsymbol{r}}}{{\boldsymbol{J}}} $$ (32) 基于上述过程, 具体迭代求解流程如图5所示.
4. 数值算例
4.1 轮胎模型验证
本算例对轮胎模型的有效性进行测试, 如图6所示, 轮胎受到纵向力F的作用. 纵向力的施加过程分为两步, 分为加载过程和卸载过程. 图7给出了文献[32]中的实验结果与本文所提出的轮胎模型数值仿真结果.
从图7中的计算结果可以看出, 本文所提轮胎模型, 在加载到20 kN时, 仿真误差与实验的误差约等于4.5%, 卸载到25 kN时, 与实验的误差约等于3.3%, 卸载到10 kN时, 与实验的误差约等于3%, 且加载过程与卸载过程与实验结果比较吻合, 从而可以验证本研究提出的轮胎模型在一定程度上具有较好的效果. 图8, 给出了本文所提轮胎模型在加载20 kN和卸载到25 kN两个不同时刻的应变云图.
4.2 整车接触动力学分析
本算例对车辆在凸起路面行驶过程的接触动力学进行测试分析, 如图9所示. 凸起位置高出地面40 mm, 车辆以不同的速度通过凸起位置时, 由于轮胎与凸起位置接触响应的作用, 整车状态会发生剧烈的变化. 本文通过对车辆以不同速度驶过凸起路面的接触动力学分析, 以定性了解在崎岖路面上车速对整车的动力学影响, 具体参数信息见表1.
表 1 整车仿真参数信息Table 1. Vehicle simulation parameter informationParameter ρ∞ time steps impact
coefficientValue 0.2 10000 0.0 Parameter friction coefficient iterations error Value 0.05 50 0.0001 Parameter rod/spring elastic
modulus/Parod/spring
radius/mdensity/
(kg·m−3)Value 2.1 × 1011/1.0 × 1011 0.015/0.001 4000 图10给出了车速分别以5 km/h和10 km/h驶过凸起路面时, 舱体上面4个顶点的位置时程曲线. 其中图10(a)和图10(b)是以5 km/h驶过时的结果, 图10(c)和图10(d)是以10 km/h驶过时的结果. 从图10中的计算结果可以看出, 当车辆以较低的速度驶过凸起路面时, 舱体上部的位置曲线发生剧烈变动的幅度比较小; 当以高速驶过时, 舱体上部的位置曲线发生剧烈变动的幅度较大. 此外, 图11给出了当车速以5 km/h和10 km/h驶过凸起路面时, 舱体上面4个顶点的速度时程曲线. 其中图11(a)和图11(b)给出了以5 km/h驶过时顶点的速度曲线, 图11(c)和图11(d)是以10 km/h驶过时的速度曲线. 从速度的计算结果可以明显看出, 车辆以低速驶过凸起位置时, 舱体上面顶点处的速度抖动比较小, 随着速度的提升, 车辆驶过凸起位置时, 舱体上面顶点的速度发生了幅度较大的抖动. 计算结果可以说明, 当车辆驶过凸起路面时, 随着车速的增加, 轮胎与凸起路面之间的接触响应越大, 接触响应传递到舱体时, 对舱体的稳定性影响越大. 为了进一步验证不同速度通过凸起路面时, 接触响应对轮胎以及整车的动力学特性的影响, 在图12中给出了车辆以不同速度驶过凸起位置时的应变情况. 其中图12 (a)为车辆所在凸起位置, 图12(b)为车速为5 km/s时的整车应变图, 图12(c)为车速为10 km/s时的整车应变图. 从图12中的计算结果可以看出, 在驶过凸起位置时, 车速为10 km/s时整车的最大应变为车速为5 km/s时的3倍左右. 基于上述结果可以得出, 在崎岖路面上, 车速越大, 地面的接触响应对整车动力学行为影响越大, 此种结果与实际情况相符, 从而可说明本文所提方法的有效性.
![]()
图 10 不同速度下车辆通过凸起路面舱体顶点位置变化对比Figure 10. Comparison of changes in the vertex position of the cabin when vehicles pass through a raised road surface at different velocity![]()
图 11 不同速度下车辆通过凸起路面舱体顶点速度变化对比Figure 11. Comparison of changes in the vertex velocity of the cabin when vehicles pass through a raised road surface at different velocity![]()
图 12 不用速度下车辆通过突起路面时的整车应变情况Figure 12. The strain situation of the entire vehicle when passing through a raised road surface at different velocity图13给出了本文所提求解方法的收敛特性. 已有研究成果表明, 接触问题具有高度的非线性, 对其求解面临着收敛性差的问题. 针对车辆驶过凸起路面的复杂接触动力学问题, 本文算法能够以7个以内的迭代步完成收敛, 并且能够保持较小的收敛残差, 从而可以说明本文所提方法在分析车辆接触动力学时具有较好的优势. 图14给出了车辆在地面行驶过程中, 轮胎表面节点与地面接触距离和单边约束拉格朗日乘子的时程变化曲线. 当距离函数大于零时, 即轮胎表面节点与地面之间没有发生接触, 单边约束没有激活, 拉格朗日乘子等于0; 当距离函数与地面之间的距离函数等于0时, 节点对应的单边约束被激活, 此时拉格朗日乘子大于0. 从图14的计算结果可知, 本文所提出的求解方法能够很好地满足互补条件.
![]()
图 14 轮胎不同节点和地面的接触距离与法向单边约束拉格朗日乘子时程曲线Figure 14. Time history curves of the contact distance between different nodes of tires and the ground and the Lagrange multiplier of normal unilateral constraint图15(a)中给出了轮胎表面一节点(节点编号77)行驶过程中法向单边约束的拉格朗日乘子和切平面上的拉格朗日乘子时程变化曲线, 图15(b)中给出了该节点在x和y方向的速度时程曲线. 从图15中的计算结果可以看出, 当该节点对应的单边约束被激活时, 切平面上对应的两个方向的拉格朗日乘子大于零, 该节点在x和y方向速度均不为0. 从此计算结果可以说明, 当轮胎与地面之间的单边约束被激活时, 只要在切平面有速度移动, 所提方法能够准确计算出摩擦响应在切平面的分量. 图16中分别给出了车辆在非光滑路面行驶过程中, 前轮轮毂节点处的速度时程变化曲线.
5. 结 论
本文利用非光滑多体系统动力学理论研究了车辆在非光滑路面上行驶过程中的接触动力学问题. 轮胎是车辆与地面唯一的接触部位, 整车的接触动力学分析具有重要的影响. 因而, 本文首先基于张拉整体结构的思想, 提出了一种轮胎建模新方法. 其次, 利用第二类拉格朗日方程, 推导了整车系统的一般非线性动力学模型; 并基于接触的非贯穿假设, 将接触问题描述为不等式单边约束的形式, 利用拉格朗日乘子法, 将单边约束引入到动力学模型中, 建立了求解车辆系统接触问题的一般非光滑动力学模型. 根据接触响应对系统状态变量的影响, 将法向接触响应与切向接触响应按照矢量分解的方法进行解耦, 并利用牛顿碰撞定律对法向接触响应进行描述, 利用最大耗散功定理将摩擦转化为优化问题描述. 最后, 建立了接触解耦的车辆接触动力学一般数学列式, 并推导了数值离散求解过程. 通过数值算例, 一方面通过与实验结果的对比, 验证了本文提出的轮胎模型的有效性; 另一方面通过整车动力学仿真, 定性分析了路面与轮胎之间的接触响应对舱体状态的影响以及在非平顺路面行驶中, 车速对车辆动力学性能的影响; 并验证了本文方法对车辆接触动力学分析的正确性.
本文提出了一种轮胎建模和车辆接触动力学分析新方法, 并通过数值算例验证了轮胎模型的有效性以及求解方法的高效性, 该工作的计算结果可为轮胎模型建模方法的进一步研究以及整车接触动力学分析提供方法借鉴.
-
![]()
图 10 不同速度下车辆通过凸起路面舱体顶点位置变化对比
Figure 10. Comparison of changes in the vertex position of the cabin when vehicles pass through a raised road surface at different velocity
![]()
图 11 不同速度下车辆通过凸起路面舱体顶点速度变化对比
Figure 11. Comparison of changes in the vertex velocity of the cabin when vehicles pass through a raised road surface at different velocity
![]()
图 12 不用速度下车辆通过突起路面时的整车应变情况
Figure 12. The strain situation of the entire vehicle when passing through a raised road surface at different velocity
![]()
图 14 轮胎不同节点和地面的接触距离与法向单边约束拉格朗日乘子时程曲线
Figure 14. Time history curves of the contact distance between different nodes of tires and the ground and the Lagrange multiplier of normal unilateral constraint
表 1 整车仿真参数信息
Table 1 Vehicle simulation parameter information
Parameter ρ∞ time steps impact
coefficientValue 0.2 10000 0.0 Parameter friction coefficient iterations error Value 0.05 50 0.0001 Parameter rod/spring elastic
modulus/Parod/spring
radius/mdensity/
(kg·m−3)Value 2.1 × 1011/1.0 × 1011 0.015/0.001 4000 
下载: 导出CSV
-
[1] 吴锐, 于会龙, 董昊天等. 履带式特种车精细化动力学建模与仿真. 兵工学报, 2023, 44: 1-17 (Wu Rui, Yu Huilong Dong Tianhao, et al. Research on refined dynamics modeling and simulation of special tracked vehicles based on multibody dynamics theory. Acta Armamentarii, 2023, 44: 1-17 (in Chinese) Wu Rui, Yu Huilong Dong Tianhao, et al. Research on refined dynamics modeling and simulation of special tracked vehicles based on multibody dynamics theory. Acta Armamentarii, 2023, 44: 1-17 (in Chinese)
[2] 张影, 于波, 于慧力等. 发动机曲轴多体动力学仿真分析. 黑龙江科学, 2023, 14(12): 7-9 (Zhang Ying, Yu Bo, Yu Huili, et al. Multi-body dynamics simulation analysis of engine crankshaft. Heilongjiang Science. 2023, 14(12): 7-9 (in Chinese) Zhang Ying, Yu Bo, Yu Huili, et al. Multi-body dynamics simulation analysis of engine crankshaft. Heilongjiang Science. 2023, 14(12): 7-9 (in Chinese)
[3] 芮筱亭. 多体系统发射动力学进展与应用. 震动、测试与诊断, 2017, 37(2): 213-220 (Rui Xiaoting. New development in launch dynamics of multibody system and its application. Journal of Vibration, Measure & Diagnosis, 2017, 37(2): 213-220 (in Chinese) Rui Xiaoting. New development in launch dynamics of multibody system and its application. Journal of Vibration, Measure & Diagnosis, 2017, 37(2): 213-220 (in Chinese)
[4] 王宇, 祝小平, 周洲. 多体飞行器展开过程动力学特性研究. 西北工业大学学报, 2023, 41(3): 490-499 (Wang Yu, Zhu Xiaoping, Zhou Zhou. Dynamics modeling and flight characteristics of folded multi-body aircraft. Journal of Northwestern Polytechnical University, 2023, 41(3): 490-499 (in Chinese) doi: 10.1051/jnwpu/20234130490 Wang Yu, Zhu Xiaoping, Zhou Zhou. Dynamics modeling and flight characteristics of folded multi-body aircraft. Journal of Northwestern Polytechnical University, 2023, 41(3): 490-499 (in Chinese) doi: 10.1051/jnwpu/20234130490
[5] 姚远, 刘红, 朱浩等. 轨道车辆勾缓转置刚柔耦合多体动力学研究. 中南大学学报(自然科学版), 2023, 54(1): 390-400 (Yao Yuan, Liu Hong, Zhu Hao, et al. Research on rigid-flexible coupling multi-body dynamics of rail vehicle hook and relief device. Journal of Central South University (Science and Technology), 2023, 54(1): 390-400 (in Chinese) Yao Yuan, Liu Hong, Zhu Hao, et al. Research on rigid-flexible coupling multi-body dynamics of rail vehicle hook and relief device. Journal of Central South University (Science and Technology), 2023, 54(1): 390-400 (in Chinese)
[6] 李铮, 李敏, 王爱国等. 基于多体动力学的电动助力转向系统仿真与试验研究. 机械强度, 2022, 44(5): 1194-1200 (Li Zheng, Li Min, Wang Aiguo, et al. Simulation and experimental research of electric power steering system based on multibody dynamics. Journal of Mechanical Strength, 2022, 44(5): 1194-1200 (in Chinese) Li Zheng, Li Min, Wang Aiguo, et al. Simulation and experimental research of electric power steering system based on multibody dynamics. Journal of Mechanical Strength, 2022, 44(5): 1194-1200 (in Chinese)
[7] 刘辉, 符升平, 项昌乐. 不平整路面履带车辆动力传动系统扭转随机激励研究. 农业机械学报, 2010, 41(12): 1-6 (Liu Hui, Fu Shengping, Xiang Changle. Torsional random excitation of tracked vehicle powertrain in system caused by road roughness. Transactions of the Chinese Society for Agricultural, 2010, 41(12): 1-6 (in Chinese) Liu Hui, Fu Shengping, Xiang Changle. Torsional random excitation of tracked vehicle powertrain in system caused by road roughness. Transactions of the Chinese Society for Agricultural, 2010, 41(12): 1-6 (in Chinese)
[8] 张云清, 高斯, 李凌阳等. 基于多体力学的车辆动力学控制系统仿真及优化. 动力学与控制学报, 2007, 5(1): 68-74 (Zhang Yunqing, Gao Si, Li Lingyang, et al. Vehicle dynamic control system simulation and optimization using multibody dynamics. Journal of Dynamics and Control, 2007, 5(1): 68-74 (in Chinese) Zhang Yunqing, Gao Si, Li Lingyang, et al. Vehicle dynamic control system simulation and optimization using multibody dynamics. Journal of Dynamics and Control, 2007, 5(1): 68-74 (in Chinese)
[9] 武和全, 张家飞, 任起凡等. 不同角度正面碰撞中老年驾驶员下肢模型的响应研究. 汽车工程, 2021, 43(5): 739-745 (Wu Hequan, Zhang Jiafei, Ren Qifan, et al. Research on the response of the elderly driver’s lower limb model in frontal collisions at different angles. Automative Engineering, 2021, 43(5): 739-745 (in Chinese) Wu Hequan, Zhang Jiafei, Ren Qifan, et al. Research on the response of the elderly driver’s lower limb model in frontal collisions at different angles. Automative Engineering, 2021, 43(5): 739-745 (in Chinese)
[10] 王庚祥, 马道林, 刘洋等. 多体系统碰撞动力学中接触力模型的研究进展. 力学学报, 2022, 54(12): 3239-3266 (Wang Gengxiang, Ma Daolin, Liu Yang, et al. Research progress of contact force models in the collision mechanics of multibody system. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(12): 3239-3266 (in Chinese) Wang Gengxiang, Ma Daolin, Liu Yang, et al. Research progress of contact force models in the collision mechanics of multibody system. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(12): 3239-3266 (in Chinese)
[11] 王晓军, 王琪. 含摩擦与碰撞平面多刚体系统动力学线性互补算法. 力学学报, 2015, 47(5): 814-821 (Wang Xiaojun, Wang Qi. A LCP method for the dynamics of planar multibody systems with impact and friction. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(5): 814-821 (in Chinese) Wang Xiaojun, Wang Qi. A LCP method for the dynamics of planar multibody systems with impact and friction. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(5): 814-821 (in Chinese)
[12] 王检耀, 洪嘉振, 刘铸永. 接触碰撞动力学的多变量选取方法. 力学学报, 2014, 46(2): 318-322 (Wang Jianyao, Hong Jiazhen, Liu Zhuyong. Multi-variable selection method in contact/impact dynamics. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(2): 318-322 (in Chinese) Wang Jianyao, Hong Jiazhen, Liu Zhuyong. Multi-variable selection method in contact/impact dynamics. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(2): 318-322 (in Chinese)
[13] 邹铁方, 刘朱紫, 肖璟等. 一种降低人-地撞击损伤的车辆制动控制方法. 汽车工程, 2021, 43(1): 105-112 (Zou Tiefang, Liu Zhuzi, Xiao Jing, et al. A vehicle braking control method for reducing pedestrain-ground impact injury. Automative Engineering, 2021, 43(1): 105-112 (in Chinese) Zou Tiefang, Liu Zhuzi, Xiao Jing, et al. A vehicle braking control method for reducing pedestrain-ground impact injury. Automative Engineering, 2021, 43(1): 105-112 (in Chinese)
[14] 戴维, 潘勇军, 张志飞. 基于悬架摆臂移除的车辆多体系统动力学半递推建模. 机械工程学报, 2023, 59(16): 263-274 (Dai Wei, Pan Yongjun, Zhang Zhifei. Semi-recursive modeling of vehicle multibody system based on suspension-control-arm removal technique. Journal of Mechanical Engineering, 2023, 59(16): 263-274 (in Chinese) Dai Wei, Pan Yongjun, Zhang Zhifei. Semi-recursive modeling of vehicle multibody system based on suspension-control-arm removal technique. Journal of Mechanical Engineering, 2023, 59(16): 263-274 (in Chinese)
[15] Wang P, Wang G, Rui X, et. al. Contact dynamics analysis of the single-pin meshing pair of a tracked vehicle. Nonlinear Dynamics, 2021, 104: 1139-1155
[16] Ma J, Wang J, Peng J, et al. Data-driven modeling for complex contacting phenomena via improved neural networks considering link switches. Mechanism and Machine Theory, 2024, 191: 105521
[17] Ren D, Tao G, Yang X, et. al. Integration of a dissipative contact force model into vehicle-track dynamics for analysing wheel-rail dynamic interaction under short-wavelength irregularity. Vehicle System Dynamics, 2022, 60(12): 4317-4342
[18] Bruni S, Meijaard JP, Rill G, et al. State-of-the-art and challenges of railway and road vehicle dynamics with multibody dynamics approaches. Multibody System Dynamics, 2020, 49: 1-12
[19] Millan P, Pagaimo J, Magalaes H, et al. Clearance joints and friction models for the modeling of friction damped railway freight vehicles. Multibody System Dynamics, 2022, 58: 21-45
[20] 李韶华, 王伟达. 车辆动力学与控制研究进展. 动力学与控制学报, 2021, 19(3): 1-4 (Li Shaohua, Wang Weida. Research advance in vehicle dynamics and control. Journal of Dynamics and Control, 2021, 19(3): 1-4 (in Chinese) Li Shaohua, Wang Weida. Research advance in vehicle dynamics and control. Journal of Dynamics and Control, 2021, 19(3): 1-4 (in Chinese)
[21] 段顺昌, 许男, 石琴等. 考虑迟滞特性的轮胎纵向动力学模型. 机械工程学报, 2023, 59(12): 364-372 (Duan Shunchang, Xu Nan, Shi Qin, et al. Longitudinal dynamic model of tire considering hysteresis characteristics. Journal of Mechanical Engineering, 2023, 59(12): 364-372 (in Chinese) Duan Shunchang, Xu Nan, Shi Qin, et al. Longitudinal dynamic model of tire considering hysteresis characteristics. Journal of Mechanical Engineering, 2023, 59(12): 364-372 (in Chinese)
[22] 李韶华, 冯桂珍, 丁虎. 考虑胎路多点接触的电动汽车-路面耦合系统震动分析. 力学学报, 2021, 53(9): 2554-2568 (Li Shaohua, Feng Guizhen, Dinghu. Vibration analysis of electric vehicle-road coupling system considering tire-road multi-point contact. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(9): 2554-2568 (in Chinese) Li Shaohua, Feng Guizhen, Dinghu. Vibration analysis of electric vehicle-road coupling system considering tire-road multi-point contact. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(9): 2554-2568 (in Chinese)
[23] 任力惠, 李稳, 冷涵等. 轮胎式轨道交通车辆动力学研究现状与挑战. 交通运输工程学报, 2021, 21(6): 8-30 (Ren Lihui, Li Wen, Leng Han, et al. Research on dynamics of rail transit vehicle with tire running gears: state-of-arts and challenges. Journal of Traffic and Transportation Engineering, 2021, 21(6): 8-30 (in Chinese) Ren Lihui, Li Wen, Leng Han, et al. Research on dynamics of rail transit vehicle with tire running gears: state-of-arts and challenges. Journal of Traffic and Transportation Engineering, 2021, 21(6): 8-30 (in Chinese)
[24] 范新秀, 王琪. 车辆纵向非光滑多体动力学建模与数值算法研究. 力学学报, 2015, 47(2): 301-309 (Fan Xiuqi, Wang Qi. Research on modeling and simulation of longitudinal vehicle dynamics based on non-smooth dynamics of multibody systems. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(2): 301-309 (in Chinese) Fan Xiuqi, Wang Qi. Research on modeling and simulation of longitudinal vehicle dynamics based on non-smooth dynamics of multibody systems. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(2): 301-309 (in Chinese)
[25] Liu B, Vollebregt E, Bruni S. Review of conformal wheel/rail contact modeling approaches: towards the application in rail vehicle dynamics simulation. Vehicle System Dynamics, 2023, in press, doi: 10.1080/00423114.2023.2228438
[26] Sun Y, Shi F, Zhang S. et. al. Improving the robustness of non-Hertzian wheel-rail contact model for railway vehicle dynamics simulation. Multibody System Dynamics, 2023, in press, doi: 10. 1007/s1044-023-09903-x
[27] Peng Y, Chen J, Ma Y. Observer-based estimation of velocity and tire-road friction coefficient for vehicle control systems. Nonlinear Dynamics, 2019, 96(1): 363-387 doi: 10.1007/s11071-019-04794-0
[28] Kan Z, Peng H, Chen B, et al. Nonlinear dynamic and deployment analysis of clustered tensegrity structures using a positional formulation FEM. Composite Structures, 2018, 187: 241-258
[29] Song N, Peng H, Kan Z, et al. A novel nonsmooth approach for flexible multibody systems with contact and friction in 3D space. Nonlinear Dynamics, 2020, 102(3): 1375-1408
[30] Bruls O, Acary V, Cardona A. Simultaneous enforcement of constraints at position and velocity levels in the nonsmooth generalized-α scheme. Computer Methods in Applied Mechanics and Engineering, 2014, 281: 131-161 doi: 10.1016/j.cma.2014.07.025
[31] Kanzow C. Some nonanterior continuation methods for linear complementarity problems. Siam Journal on Matrix Analysis and Applications, 1996, 17(4): 851-868 doi: 10.1137/S0895479894273134
[32] 张志达, 李韶华, 周军魏. 重型汽车轮胎径向刚度实验研究. 石家庄铁道大学学报(自然科学版), 2018, 31(2): 35-39 (Zhang Zhida, Li Shaohua, Zhou Junwei. Experimental study on radial stiffness of heavy vehicle tire. Journal of Shijiazhuang Tiedao University (Nature Science Edition), 2018, 31(2): 35-39 (in Chinese) Zhang Zhida, Li Shaohua, Zhou Junwei. Experimental study on radial stiffness of heavy vehicle tire. Journal of Shijiazhuang Tiedao University (Nature Science Edition), 2018, 31(2): 35-39 (in Chinese)
-
期刊类型引用(0)
其他类型引用(2)
计量
- 文章访问数: 225
- HTML全文浏览量: 32
- PDF下载量: 101
- 被引次数: 2
