Fast Simulation Method for Pantograph and Overhead Conductor Rail System
-
摘要:
针对当前刚性接触网-受电弓系统有限元模型仿真速度慢、计算时间成本高的问题,对采用三维接触算法的弓网仿真方法与流程进行改进. 首先,采用中心差分思想,将求解弓网接触副相对运动速度时需迭代计算的方程转换为可直接计算的显式方程;然后,将刚性接触网在静平衡处进行线性化处理,以避免刚度矩阵组装耗时,并加快刚性网内力计算;其次,对弓网接触状态进行惰性判断以减少计算量;最后,对本文所提快速仿真方法在不同情况下的计算效率与精度进行分析. 研究结果表明:在30跨8 m跨距的刚性接触网-受电弓仿真算例中,快速仿真方法相比标准仿真方法节省97.67%的仿真时间,且接触力结果最大偏差仅为0.48%;随着模型规模的增大,其节省的时间迅速增加,计算效率优势愈发显著,同时接触力结果偏差均小于1.0%;随着运行速度的提高,所节省的时间占比基本不变,接触力结果偏差略有增大趋势,在230 km/h以下的速度工况中,接触力标准差偏差均小于1.0%.
Abstract:Finite element model simulation of the pantograph and overhead conductor rail (OCR) system is slow, and the time cost of calculation is high. Therefore, the simulation method and process of the pantograph and OCR system using three-dimensional contact formulation were improved. Firstly, the equation that needs to be calculated iteratively when solving the relative velocity of the contact pair between the pantograph and the OCR was replaced with an explicit equation that can be calculated directly based on the central difference method. Then, the OCR model was linearized at the static equilibrium state to avoid the time-consuming stiffness matrix assembly procedure and increase the efficiency when calculating the internal force of the OCR. Next, a lazy judgment strategy was used to estimate the contact state of the pantograph and the OCR to reduce the computational load. Finally, the computational efficiency and accuracy of the fast simulation method in different cases were analyzed. The results show that compared with the standard simulation method, the proposed fast simulation method can save 97.67% of computational time in the example of pantograph and OCR with 30 spans of 8 m, and the maximum error of contact force results is only 0.48%. With the increase in the model scale, the time saved by the fast simulation method increases sharply, and its computational efficiency advantage becomes more and more significant. Meanwhile, the errors of contact force results are all less than 1.0%. With the increase in the operation speed, the proportion of time saved by the fast simulation method remains stable, and the errors of contact force results increase slightly. At speeds under 230 km/h, the standard deviation errors of contact force are all less than 1.0%.
-
受电弓与接触网之间的稳定接触是电气化铁路安全稳定运行的重要前提. 刚性接触网具有载流量大、结构简单、易维护、可靠性高等优点,被广泛应用于城际轨道交通与山区隧道内电气化铁路. 目前,国内刚性接触网线路最高运行速度为160 km/h,如北京大兴机场线[1]、广州18号线[2]等. 近年来,各地刚性接触网线路纷纷提出提速需求,给目前刚性接触网-受电弓系统的稳定安全运行带来挑战. 因此,研究刚性接触网-受电弓系统动态性能是刚性接触网线路提速的必要前提条件.
受限于隧道线路条件以及高昂的试验成本,刚性接触网-受电弓系统动力学研究主要是基于动力学仿真开展的[3-4]. 目前,有限元方法是接触网动力学建模领域应用最广泛的方法[5-6],其在高速工况下能够得到更为精确的接触力结果[7-10]. Bautista等[11]基于梁单元建立了刚性接触网的有限元模型,并进行了刚性接触网-受电弓系统的动态接触仿真. Vera等[12]基于ANSYS软件开发了刚性接触网有限元模型. Chover等[13]在SIMPACK中建立了受电弓与刚柔过渡结构的耦合模型. Chen等[14]基于绝对节点坐标法(ANCF)建立了刚性接触网模型.
弓网动态接触力是评估弓网受流质量的重要依据[15],弓网动态接触行为的准确模拟则是弓网系统建模与仿真的关键. 部分研究采用罚函数法[16-19]模拟弓头与接触线之间的接触. 罚函数法在弓网接触位置引入了接触刚度与虚拟穿透的概念. 合适的接触刚度取值是保障罚函数法计算结果精度的关键,但目前刚性接触网-受电弓仿真中接触刚度的选取并没有统一的理论支撑,一部分沿用柔性接触网系统
50000 N/m的接触刚度[16-17],一部分则采用了更大的接触刚度[18-19]. 除罚函数法外,拉格朗日乘子约束法也被广泛用于模拟弓网接触[20-22]. 拉格朗日乘子约束法将弓网接触视为一组约束条件,将其引入弓网运动方程并进行求解,从而避免了接触刚度的选取以及接触副虚拟穿透的问题[23]. 采用罚函数法模拟接触时,一般采用Newmark法对模型进行求解[24-25]. 而拉格朗日乘子约束法则一般采用中心差分法进行求解[22,26]. 采用上述2种方法时的系统方程及求解流程均不相同.标准的刚性接触网-受电弓系统仿真方法求解效率较低. 在进行设计以及性能优化分析时,通常需要进行大批量的仿真工作,计算时间成本极大. 另外,硬件在环仿真设备也对接触网仿真程序的计算效率提出了极高挑战. 因此,弓网系统模型仿真计算效率的研究受到了相关学者的关注. 文献[27-29]建立的硬件在环仿真平台中均采用模态叠加法来模拟虚拟的接触网部分. 通过对接触网模型进行简化并缩减所计算的模态阶数,能够实现硬件在环的实时仿真,但不可避免地损失了计算精度. 文献[30-31]采用移动网格的思想,降低了弓网系统运动方程数量,从而加快仿真计算速度,但其在每个时间步均需要重新划分网格并进行插值处理,操作较为复杂,且存在一定偏差. Gregori等[32]提出了一种两阶段仿真方法,预先计算系统运动方程在单位载荷下的解,并在动态仿真阶段采用线性叠加的思想,改进仿真计算效率. 以上研究均是针对柔性接触网系统开展的,关于刚性接触网的仿真效率研究较少. Chen等[14]对刚性接触网模型进行了线性化处理,基于线性化的模型提出一种快速迭代格式,能够在几乎不损失精度的情况下,大幅减少计算时间.
以上关于弓网模型仿真效率的研究主要是针对柔性接触网进行的,且都是基于罚函数法的隐式求解格式开展的. 但目前尚无针对拉格朗日乘子约束法开发的快速仿真方法. 因此,本文在采用拉格朗日乘子约束法的弓网三维接触模型基础上,针对刚性接触网-受电弓系统的动态仿真,从优化仿真流程、加速矩阵运算2个方面对标准计算方法进行改进,提出一种显式的快速仿真方法,并分析了该方法的计算效率与计算精度.
1. 弓网系统三维接触及求解算法
文献[22]中已经详细阐述了采用ANCF单元建立刚性接触网模型的具体流程以及系统仿真流程的详细推导过程. 本文目的是在此模型与仿真流程基础上,对刚性接触网-受电弓系统仿真流程与方法进行改进,提高仿真计算效率. 因此,此处仅对标准仿真方法进行简要介绍.
1.1 三维接触算法
图1为刚性接触网接触线所受接触力在三维空间中的示意. 图中:${{\boldsymbol{v}}_{\text{L}}}$为车辆运行速度, ${\boldsymbol{F}}$为弓网间总接触力(式(1)所示),${{\boldsymbol{F}}_{\text{n}}}$为法向接触力,${{\boldsymbol{F}}_{\text{t}}}$为切向接触力,${f_x}$、${f_y}$和${f_{\textit{z}}}$分别为总接触力在3个全局坐标轴上的分量.
F=Fn+Ft=[fxfyfz]T. (1) 设受电弓上接触点的全局坐标为${\text{(}}{x_{\text{P}}},{y_{\text{P}}},{{\textit{z}}_{\text{P}}}{\text{)}}$,刚性接触网上接触点的全局坐标为${\text{(}}{x_{\text{C}}},{y_{\text{C}}},{{\textit{z}}_{\text{C}}}{\text{)}}$. 考虑到在实际运行过程中受电弓始终位于刚性接触网的下方,弓网间的接触状态可归于如式(2)所示的2种情况.
{yC−yP=0,接触,yC−yP>0,离线. (2) 在不考虑弓网接触副时,刚性接触网-受电弓系统整体运动方程为
M¨e+C˙e+Q=P, (3) 式中:${\boldsymbol{M}}$为系统质量矩阵,${\boldsymbol{C}}$为系统阻尼矩阵,${\boldsymbol{Q}}$为系统弹性内力向量,${\boldsymbol{P}}$为系统所受外部载荷向量,${\boldsymbol{e}}$为系统位置向量,${\boldsymbol{\dot e}}$为系统速度向量,${\boldsymbol{\ddot e}}$为系统加速度向量.
当弓网处于接触状态时,需对式(1)施加额外约束以计算弓网间接触力. 此时,考虑接触约束的弓网系统运动方程为
{M¨e+C˙e+Q=P+GTF,TGe=0, (4) 式中:${\boldsymbol{T}}$为用于计算弓网接触状态的辅助矩阵;${\boldsymbol{G}}$为接触约束矩阵,具体形式可参考文献[22].
1.2 标准计算方法
采用拉格朗日乘子法模拟接触行为时,中心差分法的求解效果较好[26]. 将中心差分假设代入式(4)可得
{(MΔt2+C2Δt)ei+1−fiGTini=Pi−Qi+MΔt2(2ei−ei−1)+C2Δtei−1,TGi+1ei+1=0, (5) 式中:下标$i$为迭代步索引,${f_i}$为总接触力幅值,$\Delta t$为时间步长,${{\boldsymbol{n}}_i}$为总接触力方向向量.
在计算$i + 1$步信息时,${{\boldsymbol{n}}_i}$是未知的,因为弓头与刚性接触网的相对速度是未知的. 相对速度的计算式为
vi=−12ΔtGi(ei+1−ei−1)+vL. (6) 可以看出,在计算弓头与接触网相对速度时,需要$i + 1$步的坐标信息. 因此,在每个时间步需要对总接触力方向进行牛顿迭代求解. 式(5)的增量格式可写为式(7).
JΔU=ΔˆP, (7) J=[MΔt2+C2Δt−∂(fiGTini)∂ei+1−GTiniTGi+10], (8) ΔˆP=[Pi−Qi+MΔt2(2ei−ei−1)+C2Δtei−1+(fiGTini−(MΔt2+C2Δt)ei+1)−TGi+1ei+1]T, (9) ΔU=[Δei+1Δfi]T, (10) 式中:$\Delta {{\boldsymbol{e}}_{i + 1}} $、$ \Delta f_i$分别为${{\boldsymbol{e}}_{i + 1}} $、$f_i $的增量.
当弓网处于离线状态时,受电弓与刚性接触网的运动互不干涉,则式(5)中的接触约束不存在,系统状态可由式(11)直接进行显式计算.
(MΔt2+C2Δt)ei+1=Pi−Qi+MΔt2(2ei−ei−1)+C2Δtei−1. (11) 在计算过程中,先假设弓网独立运动,并对弓网接触状态进行判断,判断式为
{TGi+1ei+1,P⩽ (12) 式中:${\boldsymbol{e}}_{i + 1,{\text{P}}}$为通过式(11)计算所得的弓网系统在$i + 1$步的预测坐标向量.
刚性接触网-受电弓系统整体计算流程如图2所示,图中,$\varepsilon $为迭代计算的收敛容限.
刚性接触网在重力载荷下的初始形态计算方法可参考文献[14],此处不再赘述. 后文中称图2所示为标准方法. 可以看出,整体仿真流程中包含2个嵌套的迭代循环,其中一个是为精确求解接触力方向进行的牛顿迭代,另一个是表示受电弓随车辆运动的中心差分逐步时间积分.
2. 计算方法改进
本文建立了30跨8 m跨距的刚性接触网模型,以及DSA380受电弓的三质量块模型,并采用上述拉格朗日乘子约束法模拟弓网间动态接触,模型中设定线路沿全局直角坐标系x轴布置,垂向方向为y轴,横向为z轴. 刚性接触网最大拉出值200 m;汇流排-接触线抗拉刚度为1.76 × 108 N,抗弯刚度为 2.77 × 105 N·m2,线密度为 7.29 kg/m;悬挂结构等效刚度为6.7 × 107 N/m,等效质量为 2.8 kg. 受电弓参数如表1所示,弓网摩擦系数设定为0.3[33],运行速度为200 km/h. 仿真时间步长为5 × 10−5 s. 本文计算效率分析均基于此模型开展,仿真平台为Intel-i7-7700 (CPU),16 G,
2400 MHz (DRAM). 鉴于文献[22]中已经通过将仿真结果与实测数据对比,验证了第1节中所述模型及标准方法的准确性,在后文中采用快速仿真方法与标准方法结果的相对偏差来评价快速仿真方法的准确性.表 1 受电弓参数Table 1. Parameters of pantograph自由度 质量/kg 刚度/(N·m−1) 阻尼/(N·s·m−1) 自由度 1 7.12 9430.00 20 自由度 2 6.00 14100.00 20 自由度 3 5.80 0.01 70 2.1 接触力方向计算
由图2可以看出,在计算弓网接触力的过程中,采用牛顿迭代法对接触力进行精确求解. 原因是式(6)在计算第$i$步速度向量时需使用第$i + 1$步的坐标信息. 而在每次迭代过程中,都需要进行求偏导、组装${\boldsymbol{J}}$和$\Delta {\boldsymbol{\hat P}}$以及矩阵求逆等运算,需要消耗大量计算时间. 因此,本文对接触力方向向量求解过程进行改进.
中心差分方法基本假设为
\left\{ {\begin{array}{*{20}{l}} {{{{\boldsymbol{\ddot e}}}_i} = \dfrac{{{{\boldsymbol{e}}_{i + 1}} - 2{{\boldsymbol{e}}_i} + {{\boldsymbol{e}}_{i - 1}}}}{{\Delta {t^2}}}}, \\ {{{{\boldsymbol{\dot e}}}_i} = \dfrac{{{{\boldsymbol{e}}_{i + 1}} - {{\boldsymbol{e}}_{i - 1}}}}{{2\Delta t}}} . \end{array} } \right. (13) 根据中心差分的思想,假设在时刻$(i - 1) / 2$为$i - 1$和$i$的中间时刻,则时刻$(i - 1) / 2$的速度为
{{\boldsymbol{\dot e}}_{(i - 1) / 2}} = \frac{{{{\boldsymbol{e}}_i} - {{\boldsymbol{e}}_{i - 1}}}}{{\Delta t}}. (14) 设系统在$i - 1$至$i$时段内加速度保持不变,则时刻$i$的速度近似为
{{\boldsymbol{\dot e}}_i} = 2{{\boldsymbol{\dot e}}_{(i - 1) / 2}} - {{\boldsymbol{\dot e}}_{i - 1}}. (15) 则式(6)可替换为
{{\boldsymbol{v}}_i} = {{\boldsymbol{G}}_i}\left[\frac{{2{\text{(}}{{\boldsymbol{e}}_i} - {{\boldsymbol{e}}_{i - 1}}{\text{)}}}}{{\Delta t}} - {{\boldsymbol{\dot e}}_{i - 1}}\right] + {{\boldsymbol{v}}_{\text{L}}}. (16) 显然,式(16)中计算相对速度所需的信息均为已知的,可以直接求解,不需迭代. 采用此种方式计算相对速度时,${{\boldsymbol{n}}_i}$与${{\boldsymbol{e}}_{i + 1}}$无关,式(5)可以写作
{{\boldsymbol{J}}_{\text{s}}}{\boldsymbol{U}} = {{\boldsymbol{P}}_{\text{s}}}, (17) 式中:
{{\boldsymbol{J}}_{\text{s}}} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{\boldsymbol{M}}}{{\Delta {t^2}}} + \dfrac{{\boldsymbol{C}}}{{2\Delta t}}}&{ - {\boldsymbol{G}}_i^{\text{T}}{{\boldsymbol{n}}_i}} \\ {{\boldsymbol{T}}{{\boldsymbol{G}}_{i + 1}}}&{\boldsymbol{0}} \end{array}} \right], (18) {\boldsymbol{U}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{e}}_{i + 1}}} & {{f_i}} \end{array}} \right]^{\mathrm{T}}, (19) {{\boldsymbol{P}}_{\text{s}}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{P}}_i} - {{\boldsymbol{Q}}_i} + \dfrac{{\boldsymbol{M}}}{{\Delta {t^2}}}{\text{(}}2{{\boldsymbol{e}}_i} - {{\boldsymbol{e}}_{i - 1}}{\text{)}} + \dfrac{{\boldsymbol{C}}}{{2\Delta t}}{{\boldsymbol{e}}_{i - 1}}} \quad 0 \end{array}} \right]^{\mathrm{T}}. (20) 可以看出,式(17)为显式格式,可直接求解出${{\boldsymbol{e}}_{i + 1}}$和${f_i}$. 此外,在仿真过程中,${{\boldsymbol{J}}_{\text{s}}}$矩阵中只有最后一行与最后一列是变化的. 因此,在逐步时间积分过程中,仅更新${{\boldsymbol{J}}_{\text{s}}}$的最后一行与最后一列即可,可以大大减少矩阵组装所需时间.
采用式(16)~(20)的计算方法替换图2中的牛顿迭代部分(后称替换方法)进行仿真,并将结果与标准方法所得结果进行对比,如图3所示. 图中,局部放大图为绝对偏差最大处. 由图3可以看出:采用式(16) ~ (20)的显式计算替换掉牛顿迭代之后所得结果与标准方法结果几乎重合,整体曲线中观察不到两者区别,两者相对偏差绝对值最大仅为
0.00035 %;标准方法耗时2092.89 s,替换方法耗时960.13 s,仅为标准方法的45.88%. 替换后,能够在不损失计算精度的基础上,大幅提高计算效率.2.2 刚性接触网线性化
通过对仿真程序逐行计时,统计各环节所耗时间,各环节耗时占比如图4所示. 在替换方法中,由于${{\boldsymbol{J}}_{\text{s}}}$仅需更新最后一行与最后一列,其更新所耗时间可忽略不计. 但组装${{\boldsymbol{Q}}_i}$与计算$ {\boldsymbol{U}} $仍然占据了绝大部分的计算时间,其中,组装${{\boldsymbol{Q}}_i}$消耗了55.51%的计算时间. 在标准方法中,弓网系统内力为
{\boldsymbol{Q}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{K}}_{\text{C}}}}& {\boldsymbol{0}} \\ {\boldsymbol{0}} &{{{\boldsymbol{K}}_{\text{P}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{e}}_{\text{C}}}} \\ {{{\boldsymbol{e}}_{\text{P}}} - {\boldsymbol{e}}_{\text{P},{\text{s}}}} \end{array}} \right], (21) 式中:${{\boldsymbol{e}}_{\text{C}}}$和${{\boldsymbol{e}}_{\text{P}}}$分别为刚性接触网和受电弓上的节点坐标向量,${\boldsymbol{e}}_{\text{P},{\text{s}}}$为受电弓的初始节点坐标,${{\boldsymbol{K}}_{\text{C}}}$和${{\boldsymbol{K}}_{\text{P}}}$分别为刚性接触网和受电弓的刚度矩阵.
${{\boldsymbol{K}}_{\text{P}}}$为固定矩阵,而${{\boldsymbol{K}}_{\text{C}}}$与坐标向量${{\boldsymbol{e}}_{\text{C}}}$有关[34]. 在逐步时间积分的每个时间步内,由于刚性接触网单元坐标发生变动,均需要重新组装${{\boldsymbol{K}}_{\text{C}}}$并计算单元内力. 改进系统内力的计算与组装方法将能够进一步缩短仿真时间. 适当地对刚性接触网进行线性化,将能够大幅提高系统仿真效率. 因此,对刚性接触网进行线性化处理可以提高计算速度.
刚性接触网在重力作用下存在一定的垂度,此种状态称为初始静平衡状态,此时的系统状态可由静态找形计算获得. 一般来说,刚性接触网跨中位置垂向刚度最小,其振动幅度最大. 取弓网系统动态运行过程中刚性接触网156 m处(即第20跨跨中点)振动轨迹,如图5所示. 图中,将跨中点的位移按照受电弓运行位置展开,红色虚线表示跨中点在初始平衡态时的位移量基准,蓝色三维曲线表示跨中点位移轨迹,紫色曲线表示跨中点垂向位移,红色曲线表示跨中点水平方向位移,黄色曲线表示跨中点水平位移和垂向位移叠加的摆动轨迹,紫色、红色和黄色曲线与蓝色轨迹呈投影关系.
由图5可以看出:刚性接触网垂向位移整体较小,在受电弓接近所观测跨中点的过程中,跨中点位移逐渐增大,当受电弓运行至观测点正下方时达到最大值,但最大值也仅为3.66 mm;接触网横向位移幅值极小,均小于0.01 mm,几乎可以忽略. 整体看来,在弓网系统运行过程中,刚性接触网基本在静平衡位置附近振动,偏离程度较小. 因此,为尽量减小线性化处理对仿真精度的影响,本文采用刚性接触网处于初始静平衡态时的切线刚度矩阵对刚性接触网弹性内力进行计算.
将刚性接触网在初始平衡位置进行线性化之后,式(21)可改写为
{\boldsymbol{Q}} = {{\boldsymbol{Q}}_0} + {{\boldsymbol{K}}_{\text{t}}}\Delta {\boldsymbol{e}} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{K}}_{\text{C},\text{s}}{\boldsymbol{e}}_{\text{C},\text{s}}} \\ {\boldsymbol{0}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{K}}_{\text{t},\text{C}}}&{\boldsymbol{0}} \\ {\boldsymbol{0}}&{{{\boldsymbol{K}}_{\text{P}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{e}}_{\text{C}}} - {\boldsymbol{e}}_{\text{C},\text{s}}} \\ {{{\boldsymbol{e}}_{\text{P}}} - {\boldsymbol{e}}_{\text{P},\text{s}}} \end{array}} \right], (22) 式中:${{\boldsymbol{Q}}_{\text{0}}}$为初始平衡态时的系统内力,${{\boldsymbol{K}}_{\text{t}}}$为系统处于初始平衡态时的整体切线刚度矩阵,$\Delta {\boldsymbol{e}}$为系统运动过程中相对于初始平衡态的坐标增量,${\boldsymbol{K}}_{\text{C},{\text{s}}}$为刚性接触网在初始平衡态时的刚度矩阵,${\boldsymbol{e}}_{\text{C},{\text{s}}}$为刚性接触网在初始平衡态时的坐标向量,${\boldsymbol{K}}_{\text{t},{\text{C}}}$为刚性接触网在初始平衡态时的切线刚度矩阵.
通过将刚性接触网线性化,在逐步时间积分迭代中,不再需要反复组装系统刚度矩阵,只需在开始迭代之前组装一次系统切线刚度矩阵即可. 结合2.1节中所述的速度方向计算方法与此处的线性化处理,对模型进行仿真,并将结果与标准方法结果进行对比. 结果对比表明,线性化的模型结合替换方法所得的结果与标准方法结果基本重合. 但由于线性化处理不可避免地引入了额外偏差,接触力时域相对偏差绝对值最大值上升至0.48%,仍处于可接受范围内. 计算效率方面,线性化模型结合替换方法仅需496.64 s,为标准方法的23.73%,为替换后方法的51.73%. 可见,在采用中心差分法的显式格式计算刚性接触网-受电弓的三维动态接触行为时,对刚性接触网进行线性化处理,能够显著提高系统仿真计算效率.
2.3 其他改进措施
由图4可以看出,除了计算组装系统弹性内力向量${{\boldsymbol{Q}}_i}$耗时较多外,计算式(17)也耗费了大量时间,求解涉及大矩阵的求逆,导致计算效率较低.
刚性接触网属于细长结构,采用梁单元进行模拟时,1个节点一般只与2个单元相连,导致刚性接触网的质量矩阵和阻尼矩阵较为稀疏,非零元素基本集中在主轴附近. 进一步地,系统矩阵${\boldsymbol{J}}$与${{\boldsymbol{J}}_{\text{s}}}$也较为稀疏. 经统计可知,矩阵${{\boldsymbol{J}}_{\text{s}}}$的稀疏度为99.19%. 因此,对${{\boldsymbol{J}}_{\text{s}}}$进行稀疏处理将能够有效减小计算式(17)时所占的时间与内存. 同理,对系统切线刚度矩阵${{\boldsymbol{K}}_{\text{t}}}$进行稀疏处理,能够提高式(22)计算弹性内力的效率.
此外,矩阵LU分解能够减小式(17)的计算量,对${{\boldsymbol{J}}_{\text{s}}}$分解后,式(17)改写为
{\boldsymbol{U}} = {\boldsymbol{J}}_{\text{U}}^{ - 1}{\boldsymbol{J}}_{\text{L}}^{ - 1}{{\boldsymbol{P}}_{\text{s}}}, (23) 式中:${{\boldsymbol{J}}_{\text{U}}}$和${{\boldsymbol{J}}_{\text{L}}}$分别为${{\boldsymbol{J}}_{\text{s}}}$分解所得上三角和下三角矩阵,满足$ {{\boldsymbol{J}}_{\text{L}}} {{\boldsymbol{J}}_{\text{U}}} = {{\boldsymbol{J}}_{\text{s}}} $.
一般情况下,刚性接触网-受电弓系统离线状态的时间远小于正常接触的时间. 因此,绝大多数时候,式(12)的判断结果均为接触. 当弓网处于离线状态时,若强行施加接触约束,求解得到的接触力也将为负值. 因此,可采用惰性计算的思想,将判断接触状态这一环节调整为在计算出接触力之后,再根据接触力正负判断接触状态,从而削减判断接触状态所占用的计算时间成本支出.
2.4 快速仿真方法整体流程
本文提出的刚性接触网-受电弓系统快速仿真方法流程如图6所示. 可以看出,除了逐步时间积分之外,没有任何其他迭代循环.
快速仿真方法与标准方法结果对比如图7所示. 图中,局部放大图为绝对偏差最大处. 可以看出,快速仿真方法所得结果基本与标准方法所得结果重合,总接触力的相对偏差绝对值最大值仍为0.48%,表明稀疏矩阵处理、LU分解以及接触状态的惰性判断对计算精度没有影响. 采用标准方法计算所需时间为
2092.89 s,而快速仿真方法所需时间为48.75 s,仅为标准方法的2.33%,节省了97.67% 的计算时间成本. 通过对程序逐行计时,可得快速仿真方法各环节耗时占比如图8所示. 可以看出,在标准方法中,计算${{\boldsymbol{Q}}_i}$、$\Delta {\boldsymbol{U}}$和${\boldsymbol{U}}$等环节耗时较多,而在快速仿真方法中,其耗时占比几乎可以忽略不计.本文所采用的改进措施共有5种:① 替换速度向量计算公式、② 模型线性化、③ 稀疏处理、④ LU分解、⑤ 惰性判断接触. 为直观展示各种措施对计算结果及仿真效率的影响,表2列出了上文中各个阶段的接触力结果统计值和仿真耗时对比. 本文最终的快速仿真方法采用了上述全部5种措施. 可以看出,快速仿真方法所得结果偏差极小,其中,大部分偏差是由模型线性化引起的,替换速度向量计算公式仅引起了极其微小的偏差,而稀疏处理、LU分解、惰性判断接触则没有引起额外偏差. 各措施均能有效改善仿真效率,但随着整体仿真时间越来越短,效率改善存在边际递减效应. 这是由于仿真程序基础框架存在不可避免的基础时间成本.
表 2 各改进措施效果Table 2. Effect of different improvement measures改进方法 总接触力 仿真时间 最大相对偏差/% 平均相对偏差/% 标准差相对偏差/% 耗时/s 节省时间比/% 标准方法 2092.89 ① 3.50 × 10−4 4.80 × 10−5 5.90 × 10−5 960.13 54.12 ① + ② 0.48 0.10 0.06 496.64 76.27 ① + ② + ③ 0.48 0.10 0.06 116.83 94.42 ① + ② + ③ + ④ 0.48 0.10 0.06 56.53 97.30 ① + ② + ③ + ④ + ⑤ 0.48 0.10 0.06 48.75 97.67 3. 计算效率与精度分析
本节对所提出的刚性接触网-受电弓快速仿真方法在不同情况下的计算效率和精度进行分析. 以标准方法结果为基准,采用快速仿真方法结果与标准方法结果间的相对偏差表示快速仿真方法的计算精度.
3.1 模型规模的影响
在本文所采用的刚性接触网模型中,每个ANCF节点具有6个自由度. 因此,对于具有N个节点的刚性接触网有限元模型,其刚度矩阵和切线刚度矩阵均为6N × 6N的方阵. 理论上,随着模型自由度的增加,整体计算耗时将呈现显著的几何增长. 本文测试了200 km/h速度下单跨内单元数目分别为1~8,即模型总单元数为30~240时,快速仿真方法与标准方法的计算效率与结果偏差,如图9、10所示.
由图9可以看出,随着模型总单元数的增大,使用快速仿真方法所节省的计算时间呈几何增长. 其主要原因是增加模型单元数增大了矩阵维度,而快速仿真方法简化了模型求解过程中的矩阵计算. 模型规模越大,单元数越多,则快速仿真方法效率优势越显著.
计算精度方面,随着模型规模的增大,快速仿真方法与标准方法所得的标准差之间偏差基本保持不变,偏差最大值略微增大.
3.2 运行速度的影响
为验证本文所提出的快速仿真算法在不同速度工况下的适用性与高效性,本文分析了180~230 km/h时不同速度对计算效率与结果精度的影响,结果如图11和图12所示.
由图11可以看出,随着运行速度的提高,快速仿真方法所节省的总时间减少. 运行速度提高,导致模型所模拟的运行时间缩短,而仿真步长固定不变,故逐步时间积分的迭代步数减少,整体仿真计算时间也随之缩短. 但在每个迭代步中,快速仿真算法所节省的时间比例基本不变. 整体来看,不同速度下,快速仿真算法节省时间比基本不变.
由图12可以看出,随着运行速度的提高,快速仿真方法与标准方法的结果偏差逐渐增大. 偏差呈增大趋势的主要原因是,随着运行速度提高,刚性接触网振动幅度增大,线性化的刚性接触网模型与原始模型之间的差异逐渐表现出来. 速度在230 km/h以下时,接触力平均时域相对偏差和接触力标准差的相对偏差均小于1.0%,接触力最大时域相对偏差小于2.5%,结果偏差基本可以忽略. 并且,国内刚性接触网运营速度目前均在160 km/h及以下,因此,本文所提出的快速仿真算法精度完全满足目前的仿真需求.
4. 总 结
针对刚性接触网-受电弓系统仿真效率低下、计算时间成本高昂的问题,基于采用三维拉格朗日接触算法的弓网有限元模型,提出了一种快速仿真方法,并对快速仿真方法的计算效率与精度进行了分析. 本文改进了弓网接触副相对运动速度计算公式,并采用了惰性判断接触状态的策略,以优化仿真计算流程. 首次将刚性接触网在静平衡处进行线性化处理,并引入至采用拉格朗日乘子约束法模拟弓网三维接触的仿真中,从而避免了刚度矩阵组装耗时并加快刚性网内力计算. 本文研究结论如下:
1) 提出的快速仿真计算方法能够显著提高刚性接触网-受电弓系统仿真计算效率,同时几乎不影响计算结果精度. 在30跨8 m跨距的刚性接触网-受电弓的仿真算例中,能够节省97.67%的计算时间成本,同时接触力相对偏差最大值仅为0.48%.
2) 随着模型规模的提高,快速仿真方法相对于标准方法节省的时间比例将进一步增大. 接触力时域相对偏差最大值呈略微增大趋势,接触力标准差偏差及时域偏差的平均值则基本不变. 模型规模越大,快速仿真方法的效率优势越显著,同时能保持计算精度不变.
3) 随着车辆运行速度的提高,快速仿真方法相对于标准方法节省的时间比例基本保持不变,快速仿真方法的结果偏差呈现增大趋势. 结果对比表明,运行速度在230 km/h以下时,接触力时域偏差的平均值和接触力标准差偏差均小于1.0%,在目前刚性接触网的主要运行速度下,完全满足仿真精度需求,同时能够节省大量计算时间成本.
本文所提出的快速仿真方法仍需要进一步完善. 随着运行工况的提高,刚性接触网的振动幅度增大. 线性化模型与原始非线性模型之间的差异逐渐显现,进而导致快速仿真方法的结果偏差逐渐增大. 考虑未来可能存在的更高速度工况的仿真需求,需寻找线性化处理的改进或替代方法,保障更高速度下快速仿真方法的结果精度.
-
表 1 受电弓参数
Table 1. Parameters of pantograph
自由度 质量/kg 刚度/(N·m−1) 阻尼/(N·s·m−1) 自由度 1 7.12 9430.00 20 自由度 2 6.00 14100.00 20 自由度 3 5.80 0.01 70 表 2 各改进措施效果
Table 2. Effect of different improvement measures
改进方法 总接触力 仿真时间 最大相对偏差/% 平均相对偏差/% 标准差相对偏差/% 耗时/s 节省时间比/% 标准方法 2092.89 ① 3.50 × 10−4 4.80 × 10−5 5.90 × 10−5 960.13 54.12 ① + ② 0.48 0.10 0.06 496.64 76.27 ① + ② + ③ 0.48 0.10 0.06 116.83 94.42 ① + ② + ③ + ④ 0.48 0.10 0.06 56.53 97.30 ① + ② + ③ + ④ + ⑤ 0.48 0.10 0.06 48.75 97.67 -
[1] 杜智恒. 北京大兴国际机场线160 km/h刚性接触网系统弓网耦合受流质量分析[J]. 城市轨道交通研究,2019,22(12): 157-160.DU Zhiheng. Current collection quality of pantograph-catenary coupling for 160 km/h rigid catenary system on Beijing Daxing international airport express[J]. Urban Mass Transit, 2019, 22(12): 157-160. [2] 何国军. 160 km/h架空刚性接触网技术研究进展及展望[J]. 工程建设与设计,2022,19:131-134.HE Guojun. Research progress and prospect of 160 km/h overhead rigid catenary[J]. Construction & Design for Engineering, 2022, 19: 131-134. [3] 梅桂明,张卫华. 刚性悬挂接触网动力学研究[J]. 铁道学报,2003,25(2): 24-29. doi: 10.3321/j.issn:1001-8360.2003.02.006MEI Guiming, ZHANG Weihua. Study on dynamics of rigid suspension catenary[J]. Journal of the China Railway Society, 2003, 25(2): 24-29. doi: 10.3321/j.issn:1001-8360.2003.02.006 [4] 关金发,田志军,吴积钦. 基于弓网动力学仿真的160 km/h刚柔过渡系统方案研究[J]. 铁道学报,2018,40(9): 48-56. doi: 10.3969/j.issn.1001-8360.2018.09.007GUAN Jinfa, TIAN Zhijun, WU Jiqin. Research of 160 km/h transition structure proposal between overhead conductor rail and contact line based on dynamic simulation[J]. Journal of the China Railway Society, 2018, 40(9): 48-56. doi: 10.3969/j.issn.1001-8360.2018.09.007 [5] SONG Y, ANTUNES P, POMBO J, et al. A methodology to study high-speed pantograph-catenary interaction with realistic contact wire irregularities[J]. Mechanism and Machine Theory, 2020, 152: 103940.1-103940.18. [6] GREGORI S, GIL J, TUR M, et al. Analysis of the overlap section in a high-speed railway catenary by means of numerical simulations[J]. Engineering Structures, 2020, 221: 110963.1-110963.14. [7] 周宁,邹欢,邹栋,等. 城市轨道交通弓网系统仿真模型适应性研究[J]. 西南交通大学学报,2017,52(2): 408-415,423. doi: 10.3969/j.issn.0258-2724.2017.02.026ZHOU Ning, ZOU Huan, ZOU Dong, et al. Investigation on the applicability of pantograph and catenary model urban railway system[J]. Journal of Southwest Jiaotong University, 2017, 52(2): 408-415,423. doi: 10.3969/j.issn.0258-2724.2017.02.026 [8] SONG Y, ZHANG M J, ØISETH O, et al. Wind deflection analysis of railway catenary under crosswind based on nonlinear finite element model and wind tunnel test[J]. Mechanism and Machine Theory, 2022, 168: 104608.1-104608.22. [9] 吴凡平,徐钊,刘志刚,等. 山区峡谷地形风场下柔性接触网风振特性研究[J]. 铁道学报,2021,43(5): 47-61. doi: 10.3969/j.issn.1001-8360.2021.05.006WU Fanping, XU Zhao, LIU Zhigang, et al. Research on wind-induced vibration of flexible catenary under canyon terrain wind[J]. Journal of the China Railway Society, 2021, 43(5): 47-61. doi: 10.3969/j.issn.1001-8360.2021.05.006 [10] 熊嘉铭,徐钊,鲁小兵,等. 曲线区段受电弓-接触网系统建模及动态性能分析[J]. 铁道学报,2022,44(8): 25-35. doi: 10.3969/j.issn.1001-8360.2022.08.003XIONG Jiaming, XU Zhao, LU Xiaobing, et al. Modeling and dynamic performance analysis of pantograph-catenary system in curved railway track[J]. Journal of the China Railway Society, 2022, 44(8): 25-35. doi: 10.3969/j.issn.1001-8360.2022.08.003 [11] BAUTISTA A, MONTESINOS J, PINTADO P. Dynamic interaction between pantograph and rigid overhead lines using a coupled FEM—multibody procedure[J]. Mechanism and Machine Theory, 2016, 97: 100-111. doi: 10.1016/j.mechmachtheory.2015.10.009 [12] VERA C, SUAREZ B, PAULIN J, et al. Simulation model for the study of overhead rail current collector systems dynamics, focused on the design of a new conductor rail[J]. Vehicle System Dynamics, 2006, 44(8): 595-614. doi: 10.1080/00423110500165499 [13] CHOVER J A, SUAREZ B, RODRIGUEZ P. Simulation techniques for design of overhead conductor rail lines for speeds over 140 km/h[C]//Proceedings of the 22nd International Symposium on Dynamics of Vehicles on Roads and Tracks (IAVSD 2011), Manchester: [s.n.], 2011: 14-19. [14] CHEN L, DUAN F C, SONG Y, et al. Assessment of dynamic interaction performance of high-speed pantograph and overhead conductor rail system[J]. IEEE Transactions on Instrumentation and Measurement, 2021, 71: 9001914.1-9001914.14. [15] 周宁,王俊东,刘跃平,等. 基于图像处理技术检测弓网接触力的新方法[J]. 西南交通大学学报,2023,58(1): 1-8,57. doi: 10.3969/j.issn.0258-2724.20210509ZHOU Ning, WANG Jundong, LIU Yueping, et al. Image processing based method for measuring contact force in pantograph-catenary system[J]. Journal of Southwest Jiaotong University, 2023, 58(1): 1-8,57. doi: 10.3969/j.issn.0258-2724.20210509 [16] 关金发. 受电弓与刚性接触网动力相互作用研究[D]. 成都: 西南交通大学,2016. [17] 唐浩. AC25kV刚性接触网若干关键技术研究[D]. 成都: 西南交通大学,2018. [18] 田珂,高仕斌,于金鑫. 基于ADAMS的地铁弓网耦合仿真分析[J]. 电气化铁道,2018,29(5): 51-54. [19] 代洪宇. 200 km/h交流刚性接触网方案研究[D]. 成都: 西南交通大学,2019. [20] MASSAT J P, LAINE J P, BOBILLOT A. Pantograph-catenary dynamics simulation[J]. Vehicle System Dynamics, 2006, 44(S1): 551-559. [21] KULKARNI S, PAPPALARDO C M, SHABANA A A. Pantograph/catenary contact formulations[J]. Journal of Vibration and Acoustics, 2017, 139(1): 011010.1-011010.12. [22] CHEN L, DUAN F C, SONG Y, et al. Three-dimensional contact formulation for assessment of dynamic interaction of pantograph and overhead conductor rail system[J]. Vehicle System Dynamics, 2023, 61(9): 2432-2455. doi: 10.1080/00423114.2022.2112607 [23] PAPPALARDO C M, PATEL M D, TINSLEY B, et al. Contact force control in multibody pantograph/catenary systems[J]. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-Body Dynamics, 2016, 230(4): 307-328. doi: 10.1177/1464419315604756 [24] SONG Y, LIU Z G, WANG H R, et al. Nonlinear modelling of high-speed catenary based on analytical expressions of cable and truss elements[J]. Vehicle System Dynamics, 2015, 53(10): 1455-1479. doi: 10.1080/00423114.2015.1051548 [25] COLLINA A, BRUNI S. Numerical simulation of pantograph-overhead equipment interaction[J]. Vehicle System Dynamics, 2002, 38(4): 261-291. doi: 10.1076/vesd.38.4.261.8286 [26] CARPENTER N J, TAYLOR R L, KATONA M G. Lagrange constraints for transient finite element surface contact[J]. International Journal for Numerical Methods in Engineering, 1991, 32(1): 103-128. doi: 10.1002/nme.1620320107 [27] ZHANG W H, MEI G M, WU X J, et al. Hybrid simulation of dynamics for the pantograph-catenary system[J]. Vehicle System Dynamics, 2002, 38(6): 393-414. doi: 10.1076/vesd.38.6.393.8347 [28] FACCHINETTI A, GASPARETTO L, BRUNI S. Real-time catenary models for the hardware-in-the-loop simulation of the pantograph-catenary interaction[J]. Vehicle System Dynamics, 2013, 51(4): 499-516. doi: 10.1080/00423114.2012.748920 [29] GIL J, TUR M, CORRECHER A, et al. Hardware-in-the-loop pantograph tests using analytical catenary models[J]. Vehicle System Dynamics, 2022, 60(10): 3504-3518. doi: 10.1080/00423114.2021.1962538 [30] JIMENEZ-OCTAVIO J R, CARNICERO A, SANCHEZ-REBOLLO C, et al. A moving mesh method to deal with cable structures subjected to moving loads and its application to the catenary-pantograph dynamic interaction[J]. Journal of Sound and Vibration, 2015, 349: 216-229. doi: 10.1016/j.jsv.2015.03.051 [31] SONG Y, LIU Z G, XU Z, et al. Developed moving mesh method for high-speed railway pantograph-catenary interaction based on nonlinear finite element procedure[J]. International Journal of Rail Transportation, 2019, 7(3): 173-190. doi: 10.1080/23248378.2018.1532330 [32] GREGORI S, TUR M, NADAL E, et al. Fast simulation of the pantograph-catenary dynamic interaction[J]. Finite Elements in Analysis and Design, 2017, 129: 1-13. doi: 10.1016/j.finel.2017.01.007 [33] MEI G M. Tribological performance of rigid overhead lines against pantograph sliders under DC passage[J]. Tribology International, 2020, 151: 106538.1-106538.9. [34] BERZERI M, SHABANA A A. Development of simple models for the elastic forces in the absolute nodal co-ordinate formulation[J]. Journal of Sound and Vibration, 2000, 235(4): 539-565. doi: 10.1006/jsvi.1999.2935 -