一种极坐标系下的车辆动力学仿真方法

技术领域

本发明属于车辆动力学仿真技术领域，具体涉及一种极坐标系下的车辆动力学仿真方法。

背景技术

针对车辆动力学的核心问题是动力学建模和运动方程求解的问题。在车辆系统动力学理论中，常用的坐标系为笛卡尔坐标系，其定义三个相互垂直的向量为三个坐标方向，来描述空间物体的直线运动。在笛卡尔坐标系下，牛顿在1787年建立了牛顿方程解决了支点的运动学和动力学问题。在此基础上，刚体系统的动力学方程可以写成很多种形式，其中一种简单的被广泛实用的方程是1750年欧拉提出的牛顿-欧拉方程，其在笛卡尔坐标系下定义物体运动的坐标系在质心，将物体的平动和转动作作为独立的运动自由度来描述，采用旋转矢量的概念，利用对偶数形式，使牛顿-欧拉方程具有简明的表达形式，在开链和闭链空间机构的运动学和动力分析中得到广泛应用。牛顿-欧拉方程是欧拉运动定律的定量描述，是牛顿运动定律的延伸。

对于物体自身的旋转运动通常采用卡尔丹角、欧拉角或欧拉四元素来描述、卡尔丹角和欧拉角分别绕着三个相互垂直的旋转轴依次旋转，只是旋转的顺序不一样，但都存在旋转奇异点问题，即当任意一轴旋转90度的时候会导致该轴桶其他轴重合，欧拉旋转被重合的两根轴只能产生一个轴的旋转，为了解决旋转奇异点问题，欧拉提出欧拉四元素方法，其根据有限旋转理论，由物体瞬时旋转轴和旋转角度推导得到与此刻旋转运动相关的四个元素。但这四个元素不相互独立，需要满足一定的约束关系。Shabanan ect.分析了欧拉角和欧拉参数的表达式，并应用于车轮系统动力学的建模过程中，虽然欧拉四元素能很好地描述空间物体的旋转运动且不存在旋转奇异点问题，但引入的四元素间的约束方程将动力学微分方程变成微分代数方程，增加动力学方程求解的复杂性，降低了方程求解效率。国内有学者在笛卡尔坐标系下采用卡尔丹角建立了完整的车辆-轨道耦合动力学理论，由于旋转运动采用了小角度假设，有效避免了旋转奇异点问题，但对于空间物体大角度旋转问题具有一定的局限性。

Jalón和Bayo在1994年提出完全笛卡尔坐标系方法，属于绝对坐标系建模方法。这种方法的特点是避免使用一般笛卡尔方法中的欧拉角和欧拉参数，而是利用与刚体固结的若干参考点和参考矢量的笛卡尔坐标描述刚体的空间位置与姿态。参考点选择在铰的中心，参考矢量沿铰的转轴或平动方向，可通过多个刚体间共享位置来减少位置变量。但完全笛卡尔坐标所形成的动力学方程与一般笛卡尔方法本质相同，只是其雅克比矩阵统一在一个绝对坐标系下，便于求解计算。

选择合适的坐标系能够更好地描述车辆在空间中的运动，更能提高求解精度。现有的车辆动力学模型基本上都是建立在三维笛卡尔坐标系下，三维空间极坐标下的车辆动力学建模很难直接推导并直接应用于仿真计算中。Jalón给出了二维平面极坐标系下物体的运动方程，但推导过程比较复杂，且并未给出三维空间极坐标系下物体的运动方程。

发明内容

针对现有技术中的上述不足，本发明提供的极坐标系下的车辆动力学仿真方法针对现有的车辆动力学系统中车辆旋转运动过程中存在奇异点的问题。

为了达到上述发明目的，本发明采用的技术方案为：一种极坐标系下的车辆动力学仿真方法，包括以下步骤：

S1、设置车辆动力学的仿真总时长T和仿真步长dt；

S2、确定待仿真车辆的参数，并根据其搭建极坐标系下的车辆动力学方程；

S3、设置初始仿真时间t₀为0，对应的初始积分步数为1；

S4、根据当前仿真时间，求解搭建的车辆动力学方程，获得对应的极坐标系下的车辆自由度数据；

S5、使仿真时间增加仿真步长dt，对应的积分步数增加1，并判断增加仿真步长dt后的当前仿真时间t是否大于仿真总时长T；

若是，则进入步骤S6；

若否，则返回步骤S4；

S6、输出极坐标系下车辆动力学方程当前对应的车辆自由度数据，作为车辆动力学仿真结果。

进一步地，所述步骤S2中，待仿真车辆的参数包括车辆质量m、笛卡尔坐标系下的转动惯量J＝[J_x,J_y,J_z]^T、车辆在笛卡尔坐标系下受力点所受的力力矩车辆的空间初始位置和车辆的初始姿态R⁰＝[α₀,β₀,γ₀]^T；

其中，i为各受力点的编号，且i＝1,2,3,...M，M为受力点的总数；

α₀,β₀,γ₀分别笛卡尔坐标系下绕x,y,z轴旋转时的初始角度，且旋转顺序为先绕z轴旋转γ₀，再绕y轴旋转β₀，最后绕x轴旋转α₀；

所述步骤S2中的极坐标系包括三维空间极坐标系和车辆质心极坐标系；所述三维空间极坐标系用于描述车辆质心在三维空间中的位置，所述车辆质心极坐标系用于描述车辆在三维空间中的运动姿态。

进一步地，所述步骤S2具体为：

S21、将车辆在笛卡尔坐标系下的初始位置转换至三维空间极坐标系下，并根据其计算三维空间极坐标系下车辆位置旋转矩阵

S22、将车辆的初始姿态转换至车辆质心极坐标系中，并根据其计算车辆质心极坐标系下车辆姿态旋转矩阵

S23、根据车辆位置旋转矩阵计算三维空间极坐标系下车辆所受合力F⁰；

S24、根据车辆姿态旋转矩阵和车辆所受合力F⁰，计算车辆质心极坐标系下车辆所受合力矩M⁰；

S25、根据合力F⁰和合力矩M⁰，搭建车辆动力学方程。

进一步地，所述步骤S21中，将车辆在三维空间坐标系下的初始位置转换至三维空间极坐标系下的转换公式为：

式中，在三维空间极坐标系下，L₀为车辆初始位置P⁰沿地球半径方向的平动距离；

θ₀为车辆初始位置P⁰绕经度方向的转动角度；

为车辆初始位置P⁰绕纬度方向的转动角度；

和分别为笛卡尔坐标系中车辆初始位置的x轴、y轴和z轴坐标；

所述三维空间极坐标系下车辆位置旋转矩阵为：

所述步骤S22中，将车辆的初始姿态转换至车辆质心极坐标系中的转换公式为：

式中，在三维空间及坐标系中，为车辆绕经度方向的转动角度；ψ₀为车辆绕纬度方向的转动角度；

φ₀为车辆绕地球半径方向的转动角度；

所述车辆质心极坐标系下车辆姿态旋转矩阵

所述步骤S23中，三维空间极坐标系下车辆所受合力F⁰为：

式中，为三维空间极坐标系下车辆位置旋转矩阵；

为车辆在笛卡尔坐标系下受力点所受的力；

所述步骤S24中，车辆质心极坐标系下车辆所受合力矩M⁰：

式中，为车辆在笛卡尔坐标系下受力点的力矩。

进一步地，所述步骤S25中，搭建的车辆动力学方程为：

式中，m为车辆质量；

为三维空间极坐标系中移动自由度L的二阶导数；

为三维空间极坐标系中车辆所受合力F⁰沿移动自由度L方向的分量；

F_Q为三维空间极坐标系中车辆绕经度方向的转动θ和绕纬度方向的转动产生的离心力的合力；

为三维空间极坐标系中第一旋转自由度θ的二阶导数；

为三维空间极坐标系中车辆所受合力F⁰沿第一旋转自由度θ方向的分量；

F_Cθ为三维空间极坐标系中车辆绕经度方向的转动θ产生的科氏力；

为三维空间极坐标系中第二旋转自由度的二阶导数；

为三维空间极坐标系中车辆所受合力F⁰沿第二旋转自由度方向的分量；

为三维空间极坐标系中车辆绕纬度方向的转动产生的科氏力；

为车辆质心极坐标系中车辆绕第三旋转自由度的转动惯量；

为车辆质心极坐标系中第三旋转自由度的二阶导数；

为车辆质心极坐标系中车辆所受合力矩M⁰沿第三旋转自由度方向的分量；

J_ψ为车辆质心极坐标系中绕第四旋转自由度ψ的转动惯量；

为车辆质心极坐标系中第四旋转自由度ψ的二阶导数；

车辆质心极坐标系车辆所受合力矩M⁰沿第四旋转自由度ψ方向的分量；

J_φ为车辆质心极坐标系中绕第五旋转自由度φ的转动惯量；

为车辆质心极坐标系中第五旋转自由度φ的二阶导数；

为车辆质心极坐标系中车辆所受合力矩M⁰沿第五旋转自由度ψ方向的分量；

为三维空间极坐标系中第一旋转自由度θ的一阶导数；

为三维空间极坐标系中第二旋转自由度的一阶导数；

为三维空间极坐标系中移动自由度L的一阶导数；

为三维空间极坐标系中第一旋转自由度θ上一积分步数时的值；

为三维空间极坐标系中第二旋转自由度上一积分步数时的值。

进一步地，所述步骤S4中，根据当前仿真时间，通过四阶龙格库塔法对车辆动力学方程进行求解；

所述车辆自由度数据包括当前仿真时间t对应的积分步数下，三维空间极坐标系下车辆移动及转动自由度数据及车辆质心极坐标系下车辆转动自由度数据。

进一步地，所述三维空间极坐标系下车辆移动及转动自由度数据包括L_j,θ_j,由四阶龙格库塔法所建立的车辆动力学方程进行迭代求解得到；

其中，j为车辆动力学仿真过程中，积分步数的编号，且j＝1,2,3,...J，J为积分步数总数；

L_j为第j积分步数下，三维空间极坐标系下移动自由度L的值；

θ_j为第j积分步数下，三维空间极坐标系下第一旋转自由度θ的值；

为第j积分步数下，三维空间极坐标系下第二旋转自由度的值；

为第j积分步数下，移动自由度L的一阶导数值；

为第j积分步数下，第一旋转自由度θ的值；

为第j积分步数下，第二旋转自由度的一阶导数值；

为第j积分步数下，移动自由度L的二阶导数值；

为第j积分步数下，第一旋转自由度θ的二阶导数值；

为第j积分步数下，第二旋转自由度的二阶导数值；

所述车辆质心极坐标系下车辆转动自由度数据包括ψ_j,φ_j、由四阶龙格库塔法对所建立的车辆动力学方程进行迭代求解得到；

其中，为第j积分步数下，车辆质心极坐标系下第三旋转自由度的值；

ψ_j为第j积分步数下，车辆质心极坐标系下第四旋转自由度ψ的值；

φ_j为第j积分步数下，车辆质心极坐标系下第五旋转自由度φ的值；

为第j积分步数下，车辆质心极坐标系下第三旋转自由度的一阶导数值；

为第j积分步数下，车辆质心极坐标系下第四旋转自由度ψ的一阶导数值；

为第j积分步数下，车辆质心极坐标系下第五旋转自由度φ的一阶导数值；

为第j积分步数下，车辆质心极坐标系下第三旋转自由度的二阶导数值；

为第j积分步数下，车辆质心极坐标系下第四旋转自由度ψ的二阶导数值；

为第j积分步数下，车辆质心极坐标系下第五旋转自由度φ的二阶导数值。

进一步地，获得所述车辆在三维空间极坐标系下，车辆移动及转动自由度数据的方法具体为：

A1、根据当前仿真时间，确定当前积分步数j下的车辆位置旋转矩阵为；

式中，为第j积分步数时，车辆的第一旋转自由度θ_j和第二旋转自由度变化产生的旋转矩阵变化量；

为第j-1积分步数时在三维空间极坐标系下车辆位置的旋转矩阵；

所述表示为：

A2、根据三维空间极坐标系中车辆位置旋转矩阵确定当前积分步数下车辆所受合力F^j为：

式中，当j＝1时，

A3、将每一积分步数下的合力F^j代入车辆动力学方程，并通过四阶库塔求解器对车辆动力学方程进行求解，输出三维空间极坐标系下车辆移动及转动自由度数据。

进一步地，获得所述车辆在车辆质心极坐标系下的车辆转动自由度数据的方法具体为：

B1、根据当前仿真时间，确定当前步长下的车辆姿态旋转矩阵

式中，为第j积分步数时，车辆转动产生的角度自由度变化产生的旋转矩阵变换量；

为第j-1积分步数时的车辆姿态旋转矩阵；

所述表示为：

式中，

B2、根据车辆位置姿态矩阵确定当前积分步数下车辆质心极坐标系下车辆所受合力矩M^j为；

式中，为车辆质心极坐标系下车辆姿态旋转矩阵；

为车辆在笛卡尔坐标系下的受力点；

F^j为第j积分步数时车辆所受合力；

B3、将每一积分步数下的合力矩M^j代入车辆动力学方程，并通过四阶龙格库塔求解器求解车辆动力学方程，输出车辆质心极坐标系下的车辆转动自由度数据。

进一步地，所述步骤S6中，当需要输出当前积分步数j下，车辆各受力点在笛卡尔坐标系下的绝对坐标位置时，通过以下公式计算：

式中，为在笛卡尔坐标系中，受力点的位置矩阵；

为车辆在笛卡尔坐标系下的受力点转换到车辆质心极坐标系中时，对应的旋转矩阵；

为第j积分步下，车辆质心极坐标系中车辆姿态转置矩阵；

为第i个受力点的初始方向矢量，其中，为第i个受力点在三维空间极坐标系下的移动自由度值，上标T为转置运算符；

为第j积分步下，三维空间极坐标系中车辆位置旋转矩阵；

O_j为第j积分步数下，车辆中心点O的初始矢量方向，且O_j＝[L_j,0,0]，L_j为第j积分步数下，三维空间极坐标系下的移动自由度L的值；

其中，车辆在笛卡尔坐标系下的受力点转换到车辆质心极坐标系中时的转换公式为：

式中，在车辆质心极坐标系中，为受力点沿地球半径方向的平动距离；

为受力点绕经度方向转动的角度；

为受力点绕纬度方向转动的角度；

和分别为车辆质心所在的笛卡尔坐标系中，受力点在x轴、y轴和z轴上的坐标；

所述车辆质心极坐标系下受力点的旋转矩阵为：

本发明的有益效果为：

本发明提供的极坐标系下的车辆动力学仿真方法，针对车辆动力学系统中车辆旋转运动存在奇异点的问题，采用空间极坐标系结合牛顿定律、离心力和科氏力计算方法，给出了车辆在极坐标系下的六自由度动力学仿真方法，用于描述车辆在空间中的运动位置和姿态。利用车辆在转动过程中各瞬时转轴会产生变化、空间旋转具有继承性的原理，将传统动力学模型中的三个平动和三个转动自由度转换为一个平动五个转动自由度，能够很好地描述车辆在空间中的运动，同时能够有效解决车辆动力学系统中车体等部件空间旋转奇异点的问题。

附图说明

图1为本发明中极坐标坐标系下的车辆动力学仿真方法流程图。

图2为本发明中构建车辆动力学方程方法流程图。

图3为本发明提供的实施例中笛卡尔坐标系下受力点的位移示意图。

图4为本发明提供的实施例中笛卡尔坐标系中旋转角度示意图。

图5为本发明提供的实施例中极坐标系下受力点的位移示意图。

图6为本发明提供的实施例中极坐标系下旋转角度示意图。

图7为本发明提供的实施例中笛卡尔坐标系下受力点的轨迹示意图。

图8为本发明提供的实施例中极坐标系下受力点的轨迹示意图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

如图1所示，一种极坐标系下的车辆动力学仿真方法，包括以下步骤：

S1、设置车辆动力学的仿真总时长T和仿真步长dt；

S2、确定待仿真车辆的参数，并根据其搭建极坐标系下的车辆动力学方程；

S3、设置初始仿真时间t₀为0，对应的初始积分步数为1；

S4、根据当前仿真时间，求解搭建的车辆动力学方程，获得对应的极坐标系下的车辆自由度数据；

S5、使仿真时间增加仿真步长dt，对应的积分步数增加1，并判断增加仿真步长dt后的当前仿真时间t是否大于仿真总时长T；

若是，则进入步骤S6；

若否，则返回步骤S4；

S6、输出极坐标系下车辆动力学方程当前对应的车辆自由度数据，作为车辆动力学仿真结果。

上述步骤S2中，待仿真车辆的参数包括车辆质量m、笛卡尔坐标系下的转动惯量J＝[J_x,J_y,J_z]^T、车辆在笛卡尔坐标系下受力点所受的力力矩车辆的空间初始位置和车辆的初始姿态R⁰＝[α₀,β₀,γ₀]^T；

其中，i为各受力点的编号，且i＝1,2,3,...M，M为受力点的总数；

α₀,β₀,γ₀分别笛卡尔坐标系下绕x,y,z轴旋转时的初始角度，且旋转顺序为先绕z轴旋转γ₀，再绕y轴旋转β₀，最后绕x轴旋转α₀；

步骤S2中的极坐标系包括三维空间极坐标系和车辆质心极坐标系；三维空间极坐标系用于描述车辆质心在三维空间中的位置，车辆质心极坐标系用于描述车辆在三维空间中的运动姿态。以地球为例，单个质点在空间极坐标系中的运动包含沿地球半径方向的平动、沿经度和纬度方向的转动。

如图2所示，步骤S2具体为：

S21、将车辆在笛卡尔坐标系下的初始位置转换至三维空间极坐标系下，并根据其计算三维空间极坐标系下车辆位置旋转矩阵

其中，将车辆在三维空间坐标系下的初始位置转换至三维空间极坐标系下的转换公式为：

式中，在三维空间极坐标系下，L₀为车辆初始位置P⁰沿地球半径方向的平动距离；

θ₀为车辆初始位置P⁰绕经度方向的转动角度；

为车辆初始位置P⁰绕纬度方向的转动角度；

和分别为笛卡尔坐标系中车辆初始位置的x轴、y轴和z轴坐标；

三维空间极坐标系下车辆位置旋转矩阵为：

S22、将车辆的初始姿态转换至车辆质心极坐标系中，并根据其计算车辆质心极坐标系下车辆姿态旋转矩阵

在笛卡尔坐标系下定义的绕x,y,z轴的旋转的初始角度与极坐标系下定义的初始角度一致；因此，将车辆的初始姿态转换至车辆质心极坐标系中的转换公式为：

式中，在三维空间及坐标系中，为车辆绕经度方向的转动角度；

ψ₀为车辆绕纬度方向的转动角度；

φ₀为车辆绕地球半径方向的转动角度；

车辆质心极坐标系下车辆姿态旋转矩阵

此时，在车辆质心极坐标系下车辆的转动惯量其计算公式为：

式中，J车辆在笛卡尔坐标系下的转动惯量，且J＝[J_x,J_y,J_z]^T；

S23、根据车辆位置旋转矩阵计算三维空间极坐标系下车辆所受合力F⁰；

其中，三维空间极坐标系下车辆所受合力F⁰为：

式中，为三维空间极坐标系下车辆位置旋转矩阵；

为车辆在笛卡尔坐标系下受力点所受的力；

S24、根据车辆姿态旋转矩阵和车辆所受合力F⁰，计算车辆质心极坐标系下车辆所受合力矩M⁰；

其中，车辆质心极坐标系下车辆所受合力矩M⁰：

式中，为车辆在笛卡尔坐标系下受力点的力矩。

S25、根据合力F⁰和合力矩M⁰，搭建车辆动力学方程。

根据牛顿运动定律，结合合力F⁰和合力矩M⁰得到搭建的车辆动力学方程为：

式中，m为车辆质量；

为三维空间极坐标系中移动自由度L的二阶导数；

为三维空间极坐标系中车辆所受合力F⁰沿移动自由度L方向的分量；

F_Q为三维空间极坐标系中车辆绕经度方向的转动θ和绕纬度方向的转动产生的离心力的合力；

为三维空间极坐标系中第一旋转自由度θ的二阶导数；

为三维空间极坐标系中车辆所受合力F⁰沿第一旋转自由度θ方向的分量；

F_Cθ为三维空间极坐标系中车辆绕经度方向的转动θ产生的科氏力；

为三维空间极坐标系中第二旋转自由度的二阶导数；

为三维空间极坐标系中车辆所受合力F⁰沿第二旋转自由度方向的分量；

为三维空间极坐标系中车辆绕纬度方向的转动产生的科氏力；

为车辆质心极坐标系中车辆绕第三旋转自由度的转动惯量；

为车辆质心极坐标系中第三旋转自由度的二阶导数；

为车辆质心极坐标系中车辆所受合力矩M⁰沿第三旋转自由度方向的分量；

J_ψ为车辆质心极坐标系中绕第四旋转自由度ψ的转动惯量；

为车辆质心极坐标系中第四旋转自由度ψ的二阶导数；

车辆质心极坐标系车辆所受合力矩M⁰沿第四旋转自由度ψ方向的分量；

J_φ为车辆质心极坐标系中绕第五旋转自由度φ的转动惯量；

为车辆质心极坐标系中第五旋转自由度φ的二阶导数；

为车辆质心极坐标系中车辆所受合力矩M⁰沿第五旋转自由度ψ方向的分量；

为三维空间极坐标系中第一旋转自由度θ的一阶导数；

为三维空间极坐标系中第二旋转自由度的一阶导数；

为三维空间极坐标系中移动自由度L的一阶导数；

为三维空间极坐标系中第一旋转自由度θ上一积分步数时的值；

为三维空间极坐标系中第二旋转自由度上一积分步数时的值。

上述步骤S4具体为：根据当前仿真时间，通过四阶龙格库塔法对车辆动力学方程进行求解；

车辆自由度数据包括当前仿真时间t对应的积分步数下，三维空间极坐标系下车辆移动及转动自由度数据及车辆质心极坐标系下车辆转动自由度数据；

其中，三维空间极坐标系下车辆移动及转动自由度数据包括L_j,θ_j,由四阶龙格库塔法所建立的车辆动力学方程进行迭代求解得到；

其中，j为车辆动力学仿真过程中，积分步数的编号，且j＝1,2,3,...J，J为积分步数总数；

L_j为第j积分步数下，三维空间极坐标系下移动自由度L的值；

θ_j为第j积分步数下，三维空间极坐标系下第一旋转自由度θ的值；

为第j积分步数下，三维空间极坐标系下第二旋转自由度的值；

为第j积分步数下，移动自由度L的一阶导数值；

为第j积分步数下，第一旋转自由度θ的值；

为第j积分步数下，第二旋转自由度的一阶导数值；

为第j积分步数下，移动自由度L的二阶导数值；

为第j积分步数下，第一旋转自由度θ的二阶导数值；

为第j积分步数下，第二旋转自由度的二阶导数值；

车辆质心极坐标系下车辆转动自由度数据包括ψ_j,φ_j、由四阶龙格库塔法对所建立的车辆动力学方程进行迭代求解得到；

其中，为第j积分步数下，车辆质心极坐标系下第三旋转自由度的值；

ψ_j为第j积分步数下，车辆质心极坐标系下第四旋转自由度ψ的值；

φ_j为第j积分步数下，车辆质心极坐标系下第五旋转自由度φ的值；

为第j积分步数下，车辆质心极坐标系下第三旋转自由度的一阶导数值；

为第j积分步数下，车辆质心极坐标系下第四旋转自由度ψ的一阶导数值；

为第j积分步数下，车辆质心极坐标系下第五旋转自由度φ的一阶导数值；

为第j积分步数下，车辆质心极坐标系下第三旋转自由度的二阶导数值；

为第j积分步数下，车辆质心极坐标系下第四旋转自由度ψ的二阶导数值；

为第j积分步数下，车辆质心极坐标系下第五旋转自由度φ的二阶导数值。

具体地，获得车辆在三维空间极坐标系下，车辆移动及转动自由度数据的方法具体为：

A1、根据当前仿真时间，确定当前积分步数j下的车辆位置旋转矩阵为；

式中，为第j积分步数时，车辆的第一旋转自由度θ_j和第二旋转自由度变化产生的旋转矩阵变化量；

为第j-1积分步数时在三维空间极坐标系下车辆位置的旋转矩阵；

上述公式体现了车辆自由度数据求解过程中旋转矩阵的继承性；

表示为：

式中，

A2、根据三维空间极坐标系中车辆位置旋转矩阵确定当前积分步数下车辆所受合力F^j为：

式中，当j＝1时，

A3、将每一积分步数下的合力F^j代入车辆动力学方程，并通过四阶库塔求解器对车辆动力学方程进行求解，输出三维空间极坐标系下车辆移动及转动自由度数据。

具体地，获得车辆在车辆质心极坐标系下的车辆转动自由度数据的方法具体为：

B1、根据当前仿真时间，确定当前步长下的车辆姿态旋转矩阵

式中，为第j积分步数时，车辆转动产生的角度自由度变化产生的旋转矩阵变换量；

为第j-1积分步数时的车辆姿态旋转矩阵；

上述计算公式体现了旋转矩阵的继承性；

表示为：

式中，

B2、根据车辆位置姿态矩阵确定当前积分步数下车辆质心极坐标系下车辆所受合力矩M^j为；

式中，为车辆质心极坐标系下车辆姿态旋转矩阵；

为车辆在笛卡尔坐标系下的受力点；

F^j为第j积分步数时车辆所受合力；

B3、将每一积分步数下的合力矩M^j代入车辆动力学方程，并通过四阶龙格库塔求解器求解车辆动力学方程，输出车辆质心极坐标系下的车辆转动自由度数据。

在上述计算三维空间极坐标系中的车辆转动及移动自由度数据和计算车辆质心极坐标系中的车辆转动自由度数据时，将第j积分步数得到的车辆所受合力F^j和受合力矩M^j，代入车辆动力学方程中，通过四阶龙格库塔求解器积分得到下一积分步数j+1下的车辆各移动和转动自由度及其一阶、二阶导数值。又重新更新j+1积分步数得到的车辆所受合力F^j+1和受合力矩M^j+1，重新进行迭代运算。最终达到第J积分步数，则停止迭代运算，输出最终的自由度数据，作为车辆动力学的仿真结果。

在上述计算过程中，需要说明的是，与其他任何方法不同的是，本发明中当前积分步数的旋转矩阵，不仅与当前步长下各旋转自由度变化量有关，而且与上一步旋转矩阵有关。即在上一步旋转矩阵的基础上，叠加当前各旋转自由度变化量的旋转矩阵，只要通过控制积分步数，使当前各旋转自由度变化量小于90度，即可避免旋转奇异点的问题。

在本发明的一个实施例中，当需要输出当前积分步数j下，车辆各受力点在笛卡尔坐标系下的绝对坐标位置时，通过以下公式计算：

式中，为在笛卡尔坐标系中，受力点的位置矩阵；

为车辆在笛卡尔坐标系下的受力点转换到车辆质心极坐标系中时，对应的旋转矩阵；

为第j积分步下，车辆质心极坐标系中车辆姿态转置矩阵；

为第i个受力点的初始方向矢量，其中，为第i个受力点在三维空间极坐标系下的移动自由度值，上标T为转置运算符；

为第j积分步下，三维空间极坐标系中车辆位置旋转矩阵；

O_j为第j积分步数下，车辆中心点的矢量方向矩阵，且O_j＝[O_j,0,0]，L_j为第j积分步数下，三维空间极坐标系下的移动自由度L的值；

其中，车辆在笛卡尔坐标系下的受力点转换到车辆质心极坐标系中时的转换公式为：

式中，在车辆质心极坐标系中，为受力点沿地球半径方向的平动距离；

为受力点绕经度方向转动的角度；

为受力点绕纬度方向转动的角度；

和分别为车辆质心所在的笛卡尔坐标系中，受力点在x轴、y轴和z轴上的坐标；

车辆质心极坐标系下受力点的旋转矩阵为：

在本发明的一个实施例中，提供了验证本发明中极坐标系无奇异点六自由度车辆动力学仿真方法正确性的实例：

令在笛卡尔坐标系下，车辆的质量为车辆的质量m＝100kg，转动惯量J_x＝J_y＝J_z＝120kg·m²，受到两个方向的力，受力点坐标为(3m，4m，3m)，力坐标为(0m，1m，1m)，力车辆初始位置和姿态均为0，同时约束车辆的移动自由度仿真时间为40s，采用龙格库塔积分算法分别计算车辆的运动，并与SIMPACK软件计算结果进行对比，结果如图3-图8所示。

从图3-图8中可以看出，本发明中的无奇异点六自由度车辆动力学仿真方法能够很好地描述车辆在空间中的运动，同时能够准确地求解车辆在空间中的运动轨迹。求解的结果与在传统笛卡尔坐标系建立的车辆动力学模型完全一致。

但极坐标系下三个旋转自由度与笛卡尔坐标系的三个旋转自由度结果完全不同。笛卡尔坐标系下的旋转自由度(α,β,γ)的角度值均较大，存在多个自由度的旋转轴同时旋转90°的风险，引起旋转奇异问题。而在极坐标系下，由于当前积分步的旋转矩阵只与当前自由度的旋转角度变化量和上一步的旋转矩阵有关。只要满足当前自由度的旋转角度变化量小于90°，就不会存在旋转奇异问题。

本发明的有益效果为：

本发明提供的极坐标系下的车辆动力学仿真方法，针对车辆动力学系统中车辆旋转运动存在奇异点的问题，采用空间极坐标系结合牛顿定律、离心力和科氏力计算方法，给出了车辆在极坐标系下的六自由度动力学仿真方法，用于描述车辆在空间中的运动位置和姿态。利用车辆在转动过程中各瞬时转轴会产生变化、空间旋转具有继承性的原理，将传统动力学模型中的三个平动和三个转动自由度转换为一个平动五个转动自由度，能够很好地描述车辆在空间中的运动，同时能够有效解决车辆动力学系统中车体等部件空间旋转奇异点的问题。