CN108021747B - 一种高铁弓网动态行为动网格不平衡力消除的仿真方法 - Google Patents

Abstract

本发明公开了一种用于高铁弓网动态行为动网格不平衡力消除的仿真方法，包括以下步骤：步骤1：根据有限元和多体动力学方法构建弓网非线性仿真模型；步骤2：每个仿真步时间t中，采用动网格方法加密接触线上接触点周围的网格；步骤3：通过迭代算法消除动网格中的不平衡力；步骤4：通过Newmark‑β算法求解仿真步时间t时的响应；步骤5：重复步骤2～4，直至受电弓遍历整个接触线完成弓网动力学仿真；本发明可消除采用动网格技术产生的不平衡力，用于非线性有限元求解中，提高求解效率和求解精度。

Description

一种高铁弓网动态行为动网格不平衡力消除的仿真方法

技术领域

本发明涉及高铁弓网动态行为的仿真方法，具体涉及一种高铁弓网动态行为动网格不平衡力消除的仿真方法。

背景技术

高速铁路受电弓和接触网系统负责为高速运行的电力机车输送电能，弓网之间通过滑动接触实现电能的持续传输；为了研究弓网之前复杂的动力学行为，通常采用有限元方法构建接触网的数学模型，通过迭代求解的方法的到弓网交互作用下的动态响应；然而，由于接触网具有大跨距、长距离等特点，弓网的动态仿真的计算量是十分可观对的；为了提升弓网的计算效率，Jimenez-Octavio、Carnicero、Sanchez-Rebollo和Such等人提出了弓网系统的动网格仿真技术，该成果《A moving mesh method to deal with cablestructures subjected to moving loads and its application to the catenary-pantograph dynamic interaction》发表于期刊《Journal of Sound and Vibration》的第349卷216-229页；该方法的核心思想是在弓网动态交互仿真过程中，局部加密接触点附近的接触线网格，起到提高计算精度、提升计算效率的目的；然而，这一方法目前只能用于线性有限元模型中，无法考虑接触线大变形产生的几何非线性；当考虑接触线的几何非线性时，新生成的动网格会产生较大的不平衡力，从而引起较大的误差，甚至无法求解。

发明内容

本发明提供一种可消除采用动网格技术产生的不平衡力的仿真方法，用于非线性有限元求解中，提高求解效率和求解精度的一种高铁弓网动态行为动网格不平衡力消除的仿真方法。

本发明采用的技术方案是：一种用于高铁弓网动态行为动网格不平衡力消除的仿真方法，包括以下步骤：

步骤1：根据有限元和多体动力学方法构建弓网非线性仿真模型；

步骤2：每个仿真步时间t中，采用动网格方法加密接触线上接触点周围的网格；

步骤3：通过迭代算法消除动网格中的不平衡力；

步骤4：通过纽马克贝塔算法求解仿真步时间t时的响应；

步骤5：重复步骤2～4，直至受电弓遍历整个接触线完成弓网动力学仿真。

进一步的，所述步骤3通过迭代算法消除动网格中的不平衡力的具体过程如下：

S1：提取动网格区域固定网格节点的位移响应；

S2：根据网格节点的位移初始化动网格节点的位移；

S3：根据弓网非线性仿真模型得到广义单元刚度矩阵和不平衡力矩阵；

S4：根据步骤S1和步骤S2中得到的动、定网格节点位移得到广义局部刚度矩阵和不平衡力向量；

S5：对步骤S3和步骤S4中得到的矩阵和向量施加局部约束条件；

S6：计算动网格位移增量，更新动网格位移响应；

S7：检验收敛性，若收敛则转入步骤4，若不收敛则转入S2。

进一步的，所述步骤1具体过程如下：

S21：采用非线性索单元离散接触网的接触线和承力索，力学平衡表达式：

式中：F₁、F₂、F₃、F₄、F₅和F₆均为两个端点x、y、z轴方向上的端点力；l_x、l_y和l_z分别为索单元两个端点在x轴方向、y轴方向和z轴方向上的相对距离；L₀为索单元初始长度；w为单位长度的索单元的自重；E为索单元杨氏模量，A为索单元的截面积；

S22：采用非线性杆单元模拟接触网吊弦的非线性伸缩行为，其力学平衡表达式为：

式中：F_g1、F_g2、F_g3为x、y、z轴方向上的端点力，l_gx、l_gy和l_gz分半为两端点在x轴方向、y轴方向和z轴方向上的相对距离；L_g0为杆单元的初始长度，E_g为杆单元杨氏模量，A_g为杆单元的截面积；

S23：根据接触网系统的拓扑结构，通过有限元方法形成接触网的刚度矩阵K_c和集中质量矩阵M_c；

S24：根据步骤S21～S23得到的式(1)～式(6)和接触网的刚度矩阵K_c和集中质量矩阵M_c构建接触网的力学模型；

S25：基于多体动力学方法构建受电弓模型；

根据多体动力学方法，生成受电弓的集中质量模型：

式中：m₁、m₂、m₃分别为弓头、上框架和下框架的质量；c₁、c₂、c₃分别为弓头、上框架和下框架的阻尼；k₁、k₂、k₃分别为弓头、上框架和下框架的质量；z₁、z₂、z₃分别为弓头、上框架和下框架的位移；f₀和f_c分别为静态抬升力和接触力；

S26：通过罚函数模拟弓网耦合接触；

通过罚函数模拟弓网之间的耦合交互作用，表达式如下：

f_c＝k_s(z₁-z_c)(8)

式中：，z_c为接触线接触点处的抬升量，k_s为接触刚度；

对接触网和受电弓系统耦合，形成弓网有限元动力学方程：

式中：M、C和K分别为整体质量、阻尼和刚度矩阵，U为整体位移向量，F为外载荷向量；

和

分别为加速度和速度向量。

进一步的，所述步骤2的具体过程如下：

根据仿真步时间t和车速v更新受电弓在接触线上的接触点：

x_c＝vt (10)

式中：x_c为受电弓在接触线上的接触点x轴坐标；

定义最小网格长度l_g，变化比率为α，动网格单侧区间长度为L_c；动网格节点横坐标为：

x_m＝[x_c-l_gαⁿ,x_c-l_gα^n-1…x_c-l_gα¹,x_c-l_g,x_c,x_c+l_g,x_c…+l_gα¹…x_c+l_gα^n-1,x_c+l_gαⁿ]^T(11)

式中：x_m为动网格节点横坐标向量，n为单侧动网格节点数，满足下述关系：

l_gαⁿ≤L_c＜l_gαⁿ⁺¹ (12)。

进一步的，所述步骤3具体过程如下：

S31：提取动网格在x_c-L_c到x_c+L_c区间内x、y、z轴方向上原有定网格位移向量x_s、y_s、z_s；

S32：根据x_s、y_s、z_s以及动网格的位置向量x_m，基于线性插值法生成动网格的初始位移向量y_m和z_m；

S33：根据弓网非线性仿真模型得到广义单元刚度矩阵和不平衡力矩阵；

分别对式(1)、(2)和(3)两端积分，得到：

式中：

和

为相应的单元柔度矩阵，dL₀为初始位移增量，dF为单元不平衡力增量。

柔度矩阵求逆得到：

其中：

和

为响应的刚度矩阵，dX^e为单元节点位移增量向量，dF^e为端点力增量向量。

S34：根据动、定网格节点位移生成广义局部刚度矩阵和不平衡力向量；

根据动网格节点位移向量x_m、y_m、z_m和定网格节点位移向量x_s、y_s、z_s，生成每个网格的单元刚度矩阵dF^e、

和

根据有限元方法生成整体刚度矩阵

和

整体静力学方程为：

dF＝K_TdX+K_GdL₀ (17)

式中：dX为动网格增量,dL₀为初始长度增量；

S34：施加局部约束条件；

在所有网格节点x方向施加虚约束，即：消除矩阵K_T和向量dX对应的列向量；

保证每个单元的总初始长度之和不变，即施加约束

在定网格点x、y、z方向上施加固定约束，即：消除矩阵K_T和向量dX对应的行、列向量；

S35：计算动网格位移增量，更新动网格位移相应；

根据式(17)计算动网格增量和初始长度增量；更新动网格位移向量：x_m、y_m、z_m；

S36：判断是否满足max[dX dL₀]＜d，如果满足则转入步骤4，如果不满足则返回步骤S32。

进一步的，所述步骤4具体过程如下：

将动网格节点位移向量x_m、y_m、z_m引入到弓网整体位移向量U中，重新生成整体质量矩阵M，整体刚度矩阵K和整体阻尼矩阵C，通过Newmark-β法求解该时刻的位移响应，并进入下一时刻t+△t进行新的迭代计算。

本发明的有益效果是：

(1)本发明通过对仿真步中划分的新的移动网格进行迭代计算可消除新网格的不平衡力；

(2)本发明的仿真技术用于弓网动力学的非线性有限元仿真中，充分保证接触线的几何非线性；

(3)本发明的仿真技术用于弓网动力学的非线性有限元仿真中，可提高其求解效率和求解精度。

附图说明

图1为本发明中接触线上的动网格划分示意图。

图2为本发明中动网格各节点上的不平衡力。

图3为本发明中经迭代算法消除后的动网格节点不平衡力。

图4为本发明动网格接触力计算结果。

具体实施方式

下面结合附图和具体实施例对本发明做进一步说明。

如图1-3所示，一种用于高铁弓网动态行为动网格不平衡力消除的仿真方法，包括以下步骤：

步骤1：根据有限元和多体动力学方法构建弓网非线性仿真模型；

其具体过程如下：

S21：采用非线性索单元离散接触网的接触线和承力索，力学平衡表达式：

式中：F₁、F₂、F₃、F₄、F₅和F₆均为两个端点x、y、z轴方向上的端点力；l_x、l_y和l_z分别为索单元两个端点在x轴方向、y轴方向和z轴方向上的相对距离；L₀为索单元初始长度；w为单位长度的索单元的自重；E为索单元杨氏模量，A为索单元的截面积；

S22：采用非线性杆单元模拟接触网吊弦的非线性伸缩行为，其力学平衡表达式为：

式中：F_g1、F_g2、F_g3为x、y、z轴方向上的端点力，l_gx、l_gy和l_gz分半为两端点在x轴方向、y轴方向和z轴方向上的相对距离；L_g0为杆单元的初始长度，E_g为杆单元杨氏模量，A_g为杆单元的截面积；

S23：根据接触网系统的拓扑结构，通过有限元方法形成接触网的刚度矩阵K_c和集中质量矩阵M_c；

S24：根据步骤S21～S23得到的式(1)～式(6)和接触网的刚度矩阵K_c和集中质量矩阵M_c构建接触网的力学模型；

S25：基于多体动力学方法构建受电弓模型；

根据多体动力学方法，生成受电弓的集中质量模型：

式中：m₁、m₂、m₃分别为弓头、上框架和下框架的质量；c₁、c₂、c₃分别为弓头、上框架和下框架的阻尼；k₁、k₂、k₃分别为弓头、上框架和下框架的质量；z₁、z₂、z₃分别为弓头、上框架和下框架的位移；f₀和f_c分别为静态抬升力和接触力；

S26：通过罚函数模拟弓网耦合接触；

通过罚函数模拟弓网之间的耦合交互作用，表达式如下：

f_c＝k_s(z₁-z_c) (8)

式中：，z_c为接触线接触点处的抬升量，k_s为接触刚度；

通过式(8)对接触网和受电弓系统耦合，形成弓网有限元动力学方程：

式中：M、C和K分别为整体质量、阻尼和刚度矩阵，U为整体位移向量，F为外载荷向量；

和

分别为加速度和速度向量。

步骤2：每个仿真步时间t中，采用动网格方法加密接触线上接触点周围的网格；

根据仿真步时间t和车速v更新受电弓在接触线上的接触点：

x_c＝vt (10)

式中：x_c为受电弓在接触线上的接触点x轴坐标；

定义最小网格长度l_g，变化比率为a，动网格单侧区间长度为L_c；动网格节点横坐标为：

x_m＝[x_c-l_gαⁿ,x_c-l_gα^n-1…x_c-l_gα¹,x_c-l_g,x_c,x_c+l_g,x_c+l_gα¹…x_c+l_gα^n-1,x_c+l_gαⁿ]^T(11)

式中x_m为动网格节点横坐标向量，n为单侧动网格节点数，α¹～αⁿ为，式中xc是否为x_c；

满足下述关系：

l_gαⁿ≤L_c＜l_gαⁿ⁺¹ (12)。

步骤3：通过迭代算法消除动网格中的不平衡力；

具体过程如下：

S1：提取动网格区域固定网格节点的位移响应；

根据前仿真步(t-△t)的计算结果，提取动网格在x_c-L_c到x_c+L_c区间内x、y、z轴方向上原有定网格位移向量x_s、y_s、z_s。

S2：根据网格节点的位移初始化动网格节点的位移；

根据x_s、y_s、z_s以及动网格的位置向量x_m，基于线性插值法生成动网格的初始位移向量y_m和z_m。

S3：根据弓网非线性仿真模型得到广义单元刚度矩阵和不平衡力矩阵；

分别对式(1)、(2)和(3)两端积分，得到：

式中：

和

为相应的单元柔度矩阵，dL₀为初始位移增量，dF为单元不平衡力增量。

柔度矩阵求逆得到：

其中：

和

为响应的刚度矩阵，dX^e为单元节点位移增量向量，dF^e为端点力增量向量；

S4：根据步骤S1和步骤S2中得到的动、定网格节点位移得到广义局部刚度矩阵和不平衡力向量；

根据动网格节点位移向量x_m、y_m、z_m和定网格节点位移向量x_s、y_s、z_s，生成每个网格的单元刚度矩阵dF^e、

和

根据有限元方法生成整体刚度矩阵

和

整体静力学方程为：

dF＝K_TdX+K_GdL₀ (17)

式中：dX为动网格增量,dL₀为初始长度增量。

S5：对步骤S3和步骤S4中得到的矩阵和向量施加局部约束条件；

1)在所有网格节点x方向施加虚约束，即：消除矩阵K_T和向量dX对应的列向量；

2)保证每个单元的总初始长度之和不变，即施加约束

3)在定网格点x、y、z方向上施加固定约束，即：消除矩阵K_T和向量dX对应的行、列向量。

S6：计算动网格位移增量，更新动网格位移响应；

根据式(17)计算动网格增量和初始长度增量；更新动网格位移向量：x_m、y_m、z_m。

S7：检验收敛性

判断是否满足max[dX dL₀]＜d，如果满足则转入步骤4，如果不满足则返回步骤S32。

步骤4：通过Newmark-β算法求解仿真步时间t时的响应；

将动网格节点位移向量x_m、y_m、z_m引入到弓网整体位移向量U中，重新生成整体质量矩阵M，整体刚度矩阵K和整体阻尼矩阵C，通过纽马克贝塔Newmark-β法求解该时刻的位移响应，并进入下一时刻t+△t进行新的迭代计算。

步骤5：重复步骤2～4，直至受电弓遍历整个接触线完成弓网动力学仿真。

下面以仿真列车运行速度为350km/h的京津线弓网系统动力学行为为例，按照步骤1的方法构建弓网动力学模型；具体模型构建方法可参考Yang Song、Zhigang Liu、Hongrui Wang和Xiaobing Lu等发表在《Vehicle System Dynamics》上第53卷1455-1479页的成果《Nonlinear modelling of high-speed catenary based on analyticalexpressions of cable and truss elements》；在步骤2中当仿真时间步t＝3.3s时，根据式(10)计算得弓网接触点在320.83m处；本实例中取L_c＝6.5m、l_g＝0.2m、α＝2，从t-△t时刻到t时刻，接触线上的动网格划分示意图如图1所示；在步骤2中，根据式(11)和式(12)，形成的动网格横坐标向量x_m为：

x_m＝[-6.2 -3 -1.4 -0.6 -0.2 0 0.2 0.6 1.4 3 6.2]

随后，根据前仿真步(t-△t)的计算结果，提取在动网格区间内原有定网格的位移响应x_s、y_s、z_s；其x方向上的定网格横坐标向量x_s为：

x_s＝[-7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7]

通过线性插值法确定动网格在y和z方向的位移向量y_m和z_m；此时，动网格的节点不平衡力如图2所示，可以看出，节点1和11处的不平衡力已经超过15N，继续求解会引起较大误差，甚至会导致不收敛。

然后进入步骤3，在动网格区域内，同时考虑定网格向量x_s、y_s、z_s和动网格向量x_m、y_m、z_m，根据式(13)～式(16)生成该区间内每个网格的单元刚度矩阵dF^e、

和

根据有限元方法生成整体刚度矩阵

和

整体静力学方程可以写成：

其中，[K_T K_G]为69×91的矩阵，不是方阵，无法直接求解；采用步骤3中的1)和2)虚约束条件，使得该刚度矩阵转化为69×69的方阵；再施加3)中的固定约束条件，进一步减少其自由度，使该刚度矩阵减少为24×24的方阵，从而能够进行求解。

本实施例中设定阈值d＝10^-5，通过迭代计算，当步骤4中的收敛条件满足时，动网格节点的不平衡力如图3所示；可以看出，最大的节点不平衡力仅为3×10^-4N，通过该算法可以成功地消除动网格节点的不平衡力。

通过步骤5中的迭代求解计算，可以最终输出弓网间的接触压力曲线；图4给出了本发明的计算结果和传统定网格计算结果的比较；可以看出，无论是接触力的波形、波动范围，本发明的计算结果和传统定网格计算结果具有高度的一致性；接触压力统计量如表1所示。

表1.接触压力统计量表

从表1可以看出，本发明的接触力计算结果和传统定网格计算结果相差不大，尤其是接触力标准差、均值等较为重要的统计量具有高度的一致性，证明了本发明的有效性；此外，本发明的计算时间仅需要2200s，相比于传统的定网格计算方法节省了超过40％的时间；体现了本发明的高效性。

本发明仿真方法能够在保证几何非线性的同时有效消除新网格的不平衡力，并且将本发明仿真方法扩展到非线性有限元求解中，可提升求解效率和求解精度。

Claims (5)

1.一种用于高铁弓网动态行为动网格不平衡力消除的仿真方法，其特征在于，包括以下步骤：

步骤1：根据有限元和多体动力学方法构建弓网非线性仿真模型；

步骤2：每个仿真步时间t中，采用动网格方法加密接触线上接触点周围的网格；

步骤3：通过迭代算法消除动网格中的不平衡力，具体为：

S1：提取动网格区域固定网格节点的位移响应；

S2：根据网格节点的位移初始化动网格节点的位移；

S3：根据弓网非线性仿真模型得到广义单元刚度矩阵和不平衡力矩阵；

S4：根据步骤S1和步骤S2中得到的动、定网格节点位移得到广义局部刚度矩阵和不平衡力向量；

S5：对步骤S3和步骤S4中得到的矩阵和向量施加局部约束条件；

S6：计算动网格位移增量，更新动网格位移响应；

S7：检验收敛性，若收敛则转入步骤4，若不收敛则转入S2；

步骤4：通过纽马克贝塔算法求解仿真步时间t时的响应；

步骤5：重复步骤2～4，直至受电弓遍历整个接触线完成弓网动力学仿真。

2.根据权利要求1所述的一种用于高铁弓网动态行为动网格不平衡力消除的仿真方法，其特征在于，所述步骤1具体过程如下：

S21：采用非线性索单元离散接触网的接触线和承力索，力学平衡表达式：

式中：F₁、F₂、F₃、F₄、F₅和F₆均为两个端点x、y、z轴方向上的端点力；l_x、l_y和l_z分别为索单元两个端点在x轴方向、y轴方向和z轴方向上的相对距离；L₀为索单元初始长度；w为单位长度的索单元的自重；E为索单元杨氏模量，A为索单元的截面积；

S22：采用非线性杆单元模拟接触网吊弦的非线性伸缩行为，其力学平衡表达式为：

式中：F_g1、F_g2、F_g3为x、y、z轴方向上的端点力，L_gx、L_gy和L_gz分别为两端点在x轴方向、y轴方向和z轴方向上的相对距离；L_g0为杆单元的初始长度，E_g为杆单元杨氏模量，A_g为杆单元的截面积；

S23：根据接触网系统的拓扑结构，通过有限元方法形成接触网的刚度矩阵K_c和集中质量矩阵M_c；

S24：根据步骤S21～S23得到的式(1)～式(6)和接触网的刚度矩阵K_c和集中质量矩阵M_c构建接触网的力学模型；

S25：基于多体动力学方法构建受电弓模型；

根据多体动力学方法，生成受电弓的集中质量模型：

式中：m₁、m₂、m₃分别为弓头、上框架和下框架的质量；c₁、c₂、c₃分别为弓头、上框架和下框架的阻尼；k₁、k₂、k₃分别为弓头、上框架和下框架的质量；z₁、z₂、z₃分别为弓头、上框架和下框架的位移；f₀和f_c分别为静态抬升力和接触力；

S26：通过罚函数模拟弓网耦合接触；

通过罚函数模拟弓网之间的耦合交互作用，表达式如下：

f_c＝k_s(z₁-z_c) (8)

式中：z_c为接触线接触点处的抬升量，k_s为接触刚度；

对接触网和受电弓系统耦合，形成弓网有限元动力学方程：

式中：M、C和K分别为整体质量、阻尼和刚度矩阵，U为整体位移向量，F为外载荷向量；

和

分别为加速度和速度向量。

3.根据权利要求2所述的一种用于高铁弓网动态行为动网格不平衡力消除的仿真方法，其特征在于，所述步骤2的具体过程如下：

根据仿真步时间t和车速v更新受电弓在接触线上的接触点：

x_c＝vt (10)

式中：x_c为受电弓在接触线上的接触点x轴坐标；

定义最小网格长度l_g，变化比率为α，动网格单侧区间长度为L_c；动网格节点横坐标为：

x_m＝[x_c-l_gαⁿ,x_c-l_gα^n-1，...，x_c-l_gα¹,x_c-l_g,x_c,x_c+l_g,x_c+l_gα¹，...，x_c+l_gα^n-1,x_c+l_gαⁿ]^T (11)

式中：x_m为动网格节点横坐标向量，n为单侧动网格节点数，满足下述关系：

l_gαⁿ≤L_c＜l_gαⁿ⁺¹ (12)。

4.根据权利要求3所述的一种用于高铁弓网动态行为动网格不平衡力消除的仿真方法，其特征在于，所述步骤3具体过程如下：

S31：提取动网格在x_c-L_c到x_c+L_c区间内x、y、z轴方向上原有定网格位移向量x_s、y_s、z_s；

S32：根据x_s、y_s、z_s以及动网格的位置向量x_m，基于线性插值法生成动网格的初始位移向量y_m和z_m；

S33：根据弓网非线性仿真模型得到广义单元刚度矩阵和不平衡力矩阵；

分别对式(1)、(2)和(3)两端积分，得到：

式中：

和

为与dF和dL₀对应的单元柔度矩阵，dL₀为初始位移增量，dF为单元不平衡力增量；

柔度矩阵求逆得到：

其中：

和

为响应的刚度矩阵，dX^e为单元节点位移增量向量，dF^e为端点力增量向量；

其中，x₁、y₁、z₁和x₂、y₂、z₂分别是索单元两个端点在x、y、z轴方向上的位移；

S34：根据动、定网格节点位移生成广义局部刚度矩阵和不平衡力向量；

根据动网格节点位移向量x_m、y_m、z_m和定网格节点位移向量x_s、y_s、z_s，生成每个网格的单元刚度矩阵dF^e、

和

根据有限元方法生成整体刚度矩阵

和

整体静力学方程为：

dF＝K_TdX+K_GdL₀ (17)

式中：dX为动网格增量,dL₀为初始长度增量；

S34：施加局部约束条件；

在所有网格节点x方向施加虚约束，即：消除矩阵K_T和向量dX对应的列向量；

保证每个单元的总初始长度之和不变，即施加约束

其中，

为单个动网格单元长度；

在定网格点x、y、z方向上施加固定约束，即：消除矩阵K_T和向量dX对应的行、列向量；

S35：计算动网格位移增量，更新动网格位移相应；

根据式(17)计算动网格增量和初始长度增量；更新动网格节点位移向量：x_m、y_m、z_m；

S36：判断是否满足max[dX，dL₀]＜d，其中，d为预设允许误差阈值；如果满足则转入步骤4，如果不满足则返回步骤S32。

5.根据权利要求4所述的一种用于高铁弓网动态行为动网格不平衡力消除的仿真方法，其特征在于，所述步骤4具体过程如下：

将动网格节点位移向量x_m、y_m、z_m引入到弓网整体位移向量U中，重新生成整体质量矩阵M，整体刚度矩阵K和整体阻尼矩阵C，通过纽马克贝塔算法Newmark-β法求解仿真步时间t时的位移响应，并进入下一时刻t+△t进行新的迭代计算。

Images