CN106897518A

CN106897518A - 车辆‑道路耦合振动***联合仿真方法

Info

Publication number: CN106897518A
Application number: CN201710098373.0A
Authority: CN
Inventors: 杨啟梁; 李星宇
Original assignee: Wuhan University of Science and Engineering WUSE
Current assignee: Wuhan University of Science and Engineering WUSE; Wuhan University of Science and Technology WHUST
Priority date: 2017-02-23
Filing date: 2017-02-23
Publication date: 2017-06-27

Abstract

本发明涉及一种车辆‑道路耦合振动***联合仿真方法，其包括如下步骤：S1、按照动力学方程建立车辆动力学模型；S2、按照结构力学有限元方程建立路面动力学模型；S3、在所述路面动力学模型上加载所述车辆动力学模型，生成轮胎‑路面接触模型，并完成车辆‑道路耦合振动***的动力学分析。本发明采用高效的工程软件Matlab和Ansys实现车辆‑道路耦合***的联合建模，以Matlab为控制软件，通过对Ansys进行二次开发，实现Matlab和Ansys的联合仿真，完成复杂道路和车辆条件下车辆‑道路耦合振动***的动力学分析。

Description

车辆-道路耦合振动***联合仿真方法

技术领域

本发明涉及车辆工程领域，尤其涉及一种车辆-道路耦合振动***联合仿真方法。

背景技术

传统的车辆动力学和路面动力学都不研究车辆与路面之间的耦合关系。在车辆动力学中，通常以路面不平度为输入来研究车辆的振动问题，即车辆行驶的平顺性，而路面动力学则是在路面模型上加载车辆静载荷或者动载荷，研究路面的动力响应，验算路面的结构强度。实际上，当车辆以一定速度在公路上行驶时，车辆与路面是一个相互作用的过程，路面的不平使得行驶的车辆产生随机振动；同时，车辆载荷作用在路面上，使路面结构发生一定的变形，路面的振动又反过来影响行驶的车辆，从而形成了车辆-道路耦合的问题。

近年来，车辆-道路耦合***动力学问题越来越受到研究人员的关注。到目前为止，车辆-道路耦合***动力学的研究方法有两种：一种是通过建立车辆-道路耦合振动***的偏微分方程，常采用直接数值积分的方法来模拟车辆-道路耦合***的动力响应，需要研究者自己动手编写动力学方程的求解程序，不易保证计算精度；另一种是在Ansys中分别建立离散质量车辆模型和有限元路面模型，通过接触单元实现车辆-道路耦合关系，可进行线性或非线性瞬态求解，这种方法虽然利用了工程软件的计算能力，但不适合采用复杂的车辆动力学模型，而且接触单元无法准确模拟车辆-道路耦合关系。

因此，有必要提供一种全新的车辆-道路耦合振动***的仿真方法，以此来实现更为复杂道路和车辆条件下的车路耦合***仿真分析。

发明内容

针对上述现有技术中存在的问题，本发明提供了一种车辆-道路耦合振动***联合仿真方法，其采用高效的工程软件Matlab和Ansys实现车辆-道路耦合***的联合建模，以Matlab为控制软件，通过对Ansys进行二次开发，实现Matlab和Ansys的联合仿真，完成复杂道路和车辆条件下车辆-道路耦合振动***的动力学分析。

本发明就上述技术问题而提出的技术方案如下：

提供了一种车辆-道路耦合振动***联合仿真方法，其包括如下步骤：

S1、按照动力学方程(1)-(2)建立车辆动力学模型；所述动力学方程(1)-(2)分别如下：

其中，m_s为簧载质量；k_s为悬架刚度，c_s为悬架阻尼，x_s为车身位移，m_u为非簧载质量，k_t为轮胎刚度，c_t为轮胎阻尼，x_u为轮胎位移，x_r为路面不平度输入与路面沉降反馈的和值；

S2、按照结构力学有限元方程(3)建立路面动力学模型；所述结构力学有限元方程(3)如下：

其中，所述[M]，[C]，[K]分别为所述路面动力学模型的质量矩阵、阻尼矩阵和刚度矩阵；{u}为节点的位移向量，{F(t)}为载荷向量；

S3、在所述路面动力学模型上加载所述车辆动力学模型，生成轮胎-路面接触模型，并完成车辆-道路耦合振动***的动力学分析。

优选的，步骤S1中，所述动力学方程(1)-(2)的状态向量为：输出向量为：

其中，为车身加速度，为轮胎加速度，为相对动载荷，x_u-x_s为相对位移。

且将所述动力学方程(1)-(2)转化为状态方程(4)-(5)，所述状态方程(4)-(5)分别如下：

Y＝CX+DU (5)

其中，

优选的，步骤S1中，所述路面不平度采用公式(6)进行计算，所述公式(6)如下：

Z＝z₀sin(ωt) (6)；

其中，Z为路面不平度，z₀为不平度幅值，ω为振动圆频率，ω＝2πv/L，v为车速，L为路面平整度波长，t为运行时间，t为运动时间。

优选的，步骤S2中，采用公式(7)来计算所述阻尼矩阵[C]；所述公式(7)如下：

[C]＝a₁[M]+a₂[K] (7)；

其中，a₁，a₂为阻尼系数，采用经验公式(8)-(9)确定，所述公式(8)-(9)如下：

所述公式(8)-(9)中，ω₁，ω₂为路面结构的前两阶频率，由Ansys模态分析求得；ξ₁，ξ₂为前两阶振型的阻尼比。

优选的，步骤S1中，x_r通过公式(10)-(11)计算获得，所述公式(10)-(11)分别如下：

步骤S3中，通过公式(10)-(12)完成在所述路面动力学模型上加载所述车辆动力学模型的过程，所述公式(12)如下：

其中，F_z为所述车辆动力学模型加载在所述路面动力学模型上的载荷。

优选的，步骤S3中包括如下步骤：

S31、将所述路面动力学模型中的路面进行网格划分，形成若干行车区域，每一行车区域内均设有若干车辆轮胎与所述路面的接触面，所述接触面为矩形；每一所述接触面中均设有一记录节点；

S32、确定车辆的初始路面竖向位移、初始速度以及初始路面不平度，按照所述公式(10)-(12)计算需加载到第一接触面上的车辆载荷，并计算车辆通过所述第一接触面的记录节点时的路面竖向位移以及速度，并更新路面不平度；

S33、根据车辆通过所述第一接触面的记录节点时的路面竖向位移、速度以及更新的路面不平度，按照所述公式(10)-(12)计算需加载到第二接触面上的车辆载荷，并计算车辆通过所述第二接触面的记录节点时的路面竖向位移以及速度，并更新路面不平度；

S34、循环进行所述S32-S33，模拟完成车辆在所有行车区域内的行进过程。

优选的，所述车辆车轮与接触面的接触点位于车轮轮心正下方，车辆载荷均匀分布在所述接触面内；所述记录节点均位于接触面的中心线上。

优选的，所述路面在路面宽度方向上为4-8m，在行车方向上为10-14m，在路面深度方向上为4-8m；所述接触面的在路面宽度方向上为0.2-0.24m，在行车方向上为0.3-0.32m。

优选的，所述路面动力学模型为三维线弹性模型，其从上而下依次包括各向同性的面层、基层、底基层和土基层。

本发明的技术方案具有如下技术效果：

本发明将车辆-道路耦合***分解为“车辆动力学模型”、“路面动力学模型”和“轮胎-路面接触模型”3个子***，采用高效的工程软件Matlab和Ansys实现车辆-道路耦合***的联合建模，以Matlab为控制软件，通过对Ansys进行二次开发，实现Matlab和Ansys的联合仿真，完成车辆-道路耦合振动***的动力学分析，实现复杂道路和车辆条件下的车路耦合***仿真分析，其正确性高，操作简单。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是现有技术中车辆-道路耦合***的示意图；

图2是本发明实施例一中车辆-道路耦合振动***联合仿真方法的步骤流程图；

图3是本发明实施例一中生成轮胎-路面接触模型的步骤流程图；

图4是本发明实施例一中的车辆载荷施加到路面的过程示意图；

图5是本发明实施例一中的仿真流程图；

图6是本发明实施例一中的车身加速度a_u时间历程曲线；

图7是本发明实施例一中的竖向位移u_z时程曲线；

图8是本发明实施例一中的竖向应力σ_z时程曲线。

具体实施方式

本发明针对现有技术中存在的缺陷，提供了一种车辆-道路耦合振动***联合仿真方法。

实施例一：

图1示出了现有技术中车辆-道路耦合***的示意图，图2则示出了本发明车辆-道路耦合振动***联合仿真方法，其包括如下步骤：

其中，m_s为簧载质量；k_s为悬架刚度，c_s为悬架阻尼，x_s为车身位移，m_u为非簧载质量，k_t为轮胎刚度，c_t为轮胎阻尼，x_u为轮胎位移，x_r为路面不平度输入与路面竖向位移的和值；

优选的，该步骤S1中，根据现代控制理论，所述动力学方程(1)-(2)的状态向量为：输出向量为：

Y＝CX+DU (5)

其中，

表1为所述车辆动力学模型的参数表：

表1车辆动力学模型参数表

表2为所示路面动力学模型的参数表：

表2路面动力学模型参数表

Z＝z₀sin(ωt) (6)；

其中，Z为路面不平度，z₀为不平度幅值(可优选为0.01m)，ω为振动圆频率，ω＝2πv/L，v为车速(可优选为14m/s)，L为路面平整度波长(可优选为6m)，t为运行时间t为运动时间。

[C]＝a₁[M]+a₂[K] (7)；

其中，Z为路面不平度；

轮胎-路面接触模型采用改进的单点接触模型，具体的，如图3所示，步骤S3中包括如下步骤：

S31、如图4所示，采用solid45单元将所述路面动力学模型中的路面进行网格划分，形成若干行车区域，每一行车区域内均设有若干车辆轮胎与所述路面的接触面，所述接触面为矩形，通过每一所述行车区域的时间可为0.85s；每一所述接触面中均设有一记录节点；所述车辆车轮与接触面的接触点位于车轮轮心正下方，车辆载荷均匀分布在所述接触面内；所述记录节点均位于接触面的中心线上；优选的，针对传统车路耦合***中道路模型过于简化的问题，根据沥青路面特性和弹性层状理论体系，建立结构形式为三维线弹性模型的路面动力学模型，其从上而下依次包括各向同性的面层、基层、底基层和土基层；所述路面在路面宽度方向上为4-8m(优选6m)，在行车方向上为10-14m(优选12m)，在路面深度方向上为4-8m(优选6m)；所述接触面的在路面宽度方向上为0.2-0.24m(优选0.22m)，在行车方向上为0.3-0.32m(优选0.32m)；

由于需要Matlab调用Ansys进行联合仿真，因此须使用Ansys内部的参数化设计语言APDL编写宏文件进行道路的有限元分析。Matlab调用Ansys生成的宏文件进行道路计算的关键语法如下：

！"D:\Program Files\Ansys

Inc\v150\ANSYS\bin\intel\Ansys150.exe"-b–i

d:\inp\step_1.mac-o"d:\out\put\rst_1.txt"

其中，"！"是由Matlab提供的用以执行shell命令的操作符。"step_1"是APDL语言编写的宏文件，完成建模、加载、分析和后处理等整个仿真过程；"rst_1"是计算结果文件，保存计算过程中的信息。参数"-b"指定使用batch方式运行Ansys，"-i"指定输入文件，"-o"指定输出文件。

车辆动力学模型和路面动力学模型在simulink环境下搭建。采用Matlab作为主控程序，在batch模式下调用Ansys进行联合仿真，步骤如下(图5同时也示出了本发明中的仿真流程图)：

S32、确定车辆的初始路面竖向位移、初始速度以及初始路面不平度，并保存至外部数据交换文件，调用Ansys按照所述公式(10)-(12)计算需加载到第一接触面上的车辆载荷，并计算车辆通过所述第一接触面的记录节点时的路面竖向位移以及速度，并更新路面不平度；

具体的，Ansys运行宏文件step_1.mac，在第一接触面上进行加载和计算。Ansys计算结束后，通过*CFOPEN、*VWRITE等命令写入文件，传递计算的计算结果给外部数据交换文件。Ansys计算完成后，导出轮胎与路面接触面中心线上的记录节点的速度参数。由于*VWRITE一次最多写入19个参数，以下仅列举前11个记录节点的速度参数导出的命令流：

*DIM,z,,110,55

*DO,I,2,55

NSOL,I,17612+19*I,U,Z

VGET,z(1,I-1),I,0,

*ENDDO

！行车区域表层节点速度参数变量获取

*DIM,z_1,,55,11

*DO,I,2,12,1

NSOL,I,17612+19*I,U,Z

VGET,z_1(1,I-1),I,0,

*ENDDO

*CFOPEN,z_1,txt

*DO,I,1,55

*VWRITE,z_1(i,1),z_1(i,2),z_1(i,3),z_1(i,

4),z_1(i,5),z_1(i,6),z_1(i,7),z_1(i,8),z_1(i,9),z_1(i,

10),z_1(i,11)

(F12.10,2X,F12.10,2X,F12.10,2X,F12.10,2

X,F12.10,2X,F12.10,2X,F12.10,2X,F12.10,2X,F1

2.10,2X,F12.10,2X,F12.10)

*ENDDO

*CFCLOSE

具体的，Ansys运行结束后把第一接触面上的记录节点的路面竖向位移和速度结果传给外部文件，Simulink车辆动力学模型读取更新的路面不平度，用来计算需加载到第二接触面上的车辆载荷并导出。

用Matlab处理Ansys的输出文件，将第一接触面上的由Ansys计算的路面竖向位移和速度输入车辆动力学模型，计算需加载到第二接触面上的车辆载荷大小，并将计算的载荷大小写入外部文件，以供下一次计算时Ansys进行调用。Matlab中导出载荷文件的命令如下：

fid＝fopen('ftop.txt','wt')；

fprintf(fid,'％f\n',ftop_1)；

fclose(fid)；

具体的，调用Ansys进行计算，此时调用宏文件step_2.mac，在第二行车区域的第一接触点上进行加载和计算，加载数据由上一步生成的ftop.txt文件提供，如此循环，以此模拟模拟完成车辆在所有行车区域内的行进过程。加载命令如下：

*DO,I,1,55

TIME,0.0156*I

NSEL,R,LOC,Z,6

NSEL,R,LOC,X,-0.2,0.2

NSEL,R,LOC,Y,0.2,0.45

F,ALL,FZ,-FTOP(I,1)/10

ALLSEL

SOLVE

*ENDDO

仿真结果：

(1)车身加速度时间历程如图6所示，本发明中计算的车身加速度最大值为3.2m/s2,最小值为-3.1m/s²。引文“陈静云,孙依人.基于接触非线性的车-路耦合***动力响应分析[J].振动与冲击,2013；32(19):119-124”(以下简称引文)中最大加速度分别为2.7m/s 2和-2.8m/s²。与引文的计算结果相比，车身加速度变化规律基本一致，但加速度数值大小存在较大误差，这是由于车速为14m/s时，路面对车辆的激励频率接近车辆自振频率，因而，加速度有逐渐变大的趋势。本文计算较引文提前1/4周期，因此，加速度变化的趋势较快。

(2)竖向位移是路面设计最重要的参数，是验算路面结构的强度和刚度的关键指标。沥青路面各层竖向位移时间历程计算结果如图7所示。面层最大位移为0.39mm，与引文的结果(0.37mm)误差控制在5.4％。面层以外其余各层的位移沉降由上至下逐渐变小，其中土基层顶部的最大位移为0.35mm，由土基层引起的位移占路面位移的89.7％，由此可见对土基层的压实可有效降低路面沉降。与引文的计算结果相比，位移曲线的趋势和各层之间位移的比例关系基本相符，其误差在合理范围内。

(3)沥青路面各层中部竖向应力计算结果如图8所示。当行车载荷离开和接近计算点时，出现了短暂微小的拉应力，但路面结构主要承受的是压应力。最大压应力为-0.21MPa，与引文的结果误差控制在5％，进一步说明仿真方法的准确性。面层以下各层所受最大压应力逐层递减，至最底部的土层中部所受最大压应力为-0.018MPa。可见路面的作用就是降低路基所受的垂直压力，路面表层结构最容易发生压密形变进而形成车辙。

综上所述，本发明通过Ansys和Matlab协同仿真，进行车辆-道路耦合振动***动力学计算的方法。采用Matlab作为主控程序，以batch模式调用Ansys进行道路计算，通过读取和写入外部文件，实时传递轮胎动载和路面沉降，实现车辆-道路耦合***动力学联合仿真。通过与引文的计算结果进行对比，证明了本文提出的Ansys和Matlab联合仿真方法的正确性，为车辆-道路耦合***的动力学研究提供了一种新的方法，使用Matlab，实现复杂道路和车辆条件下的车路耦合***仿真分析，其正确性高，操作简单，为以后基于车路耦合环境下的整车平顺性研究奠定基础。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，且所述实施例中的技术特征均可进行任意组合，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims

1.一种车辆-道路耦合振动***联合仿真方法，其特征在于，包括如下步骤：

{\overset{\cdot\cdot}{m}}_{s} {\overset{\cdot\cdot}{x}}_{s} + k_{s} (x_{s} - x_{u}) + c_{s} ({\overset{\cdot}{x}}_{s} - {\overset{\cdot}{x}}_{u}) = 0 - - - (1)

\begin{matrix} {\overset{\cdot\cdot}{m}}_{u} {\overset{\cdot\cdot}{x}}_{u} + k_{s} (x_{s} - x_{u}) + c_{s} ({\overset{\cdot}{x}}_{s} - {\overset{\cdot}{x}}_{u}) + k_{t} (x_{u} - x_{r}) + \\ c_{t} ({\overset{\cdot}{x}}_{u} - {\overset{\cdot}{x}}_{r}) = 0 \end{matrix} - - - (2);

[M] {\overset{\cdot\cdot}{u}} + [C] {\overset{\cdot}{u}} + [K] {u} = {F (t)} - - - (3);

2.如权利要求1所述的方法，其特征在于，步骤S1中，所述动力学方程(1)-(2)的状态向量为：输出向量为：

其中，为车身加速度，为轮胎加速度，为相对动载荷，x_u-x_s为相对位移；。

\overset{\cdot}{X} = A X + B U - - - (4);

Y＝CX+DU (5)

A = [\begin{matrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ - \frac{k_{s} + k_{t}}{m_{u}} & \frac{k_{s}}{m_{u}} & - \frac{c_{s} + c_{t}}{m_{u}} & \frac{c_{s}}{m_{u}} \\ \frac{k_{s}}{m_{s}} & - \frac{k_{s}}{m_{s}} & \frac{c_{s}}{m_{s}} & - \frac{c_{s}}{m_{s}} \end{matrix}];

C = [\begin{matrix} \frac{k_{s}}{m_{s}} & - \frac{k_{s}}{m_{s}} & \frac{c_{s}}{m_{s}} & - \frac{c_{s}}{m_{s}} \\ - \frac{k_{s} + k_{t}}{m_{u}} & \frac{k_{s}}{m_{u}} & - \frac{c_{s} + c_{t}}{m_{u}} & \frac{c_{s}}{m_{u}} \\ \frac{k_{t}}{m_{s} g} & 0 & 0 & 0 \\ - 1 & 1 & 0 & 0 \end{matrix}];

其中，

3.如权利要求1所述的方法，其特征在于，步骤S1中，所述路面不平度采用公式(6)进行计算，所述公式(6)如下：

Z＝z₀sin(ωt) (6)；

其中，Z为路面不平度，z₀为不平度幅值，ω为振动圆频率，ω＝2πv/L，v为车速，L为路面平整度波长，t为运行时间。

4.如权利要求1所述的方法，其特征在于，步骤S2中，采用公式(7)来计算所述阻尼矩阵[C]；所述公式(7)如下：

[C]＝a₁[M]+a₂[K] (7)；

a_{1} = \frac{2 ω_{1} ω_{2} (ξ_{1} ω_{2} - ξ_{2} ω_{1})}{ω_{2}^{2} - ω_{1}^{2}} - - - (8);

a_{2} = \frac{2 (ξ_{2} ω_{2} - ξ_{1} ω_{1})}{ω_{2}^{2} - ω_{1}^{2}} - - - (9);

5.如权利要求3所述的方法，其特征在于，步骤S1中，x_r通过公式(10)-(11)计算获得，所述公式(10)-(11)分别如下：

x_{r} = Z + Σ_{i = 1}^{55} z_{i} (t) - - - (10)

{\overset{\cdot}{x}}_{r} = \overset{\cdot}{Z} + Σ_{i = 1}^{55} v_{i} (t) - - - (11);

F_{z} (t) = k_{t} (x_{r} - x_{u}) + c_{t} ({\overset{\cdot}{x}}_{r} - {\overset{\cdot}{x}}_{u}) - - - (12);

6.如权利要求5所述的方法，其特征在于，步骤S3中包括如下步骤：

7.如权利要求6所述的方法，其特征在于，所述车辆车轮与接触面的接触点位于车轮轮心正下方，车辆载荷均匀分布在所述接触面内；所述记录节点均位于接触面的中心线上。

8.如权利要求6所述的方法，其特征在于，所述路面在路面宽度方向上为4-8m，在行车方向上为10-14m，在路面深度方向上为4-8m；所述接触面的在路面宽度方向上为0.2-0.24m，在行车方向上为0.3-0.32m。

9.如权利要求1所述的方法，其特征在于，所述路面动力学模型为三维线弹性模型，其从上而下依次包括各向同性的面层、基层、底基层和土基层。