基于多柔体动力学理论的接触网找形方法

邱江洋，梅桂明，王江文，罗 群

(西南交通大学牵引动力国家重点实验室，成都 610031)

列车在高速运行时稳定的受流质量是高速列车设计的一个重要因素，而直接影响受流质量的就是受电弓与接触网之间的接触关系[1−2]。研究弓网系统之间的接触关系则需要进行弓网系统的动力学建模[3]。

在进行动力学计算之前，动力学模型的初始条件及其重要，它直接影响着动力学模型的计算结果[4−5]。对于弓网系统也是一样，接触网的初始条件表征了接触网系统在重力场、预张力和约束条件下的平衡构型。国内外众多学者对接触网系统的平衡构型进行了计算，利用悬链线方程[6−7]，将接触网系统离散为悬链线单元，忽略结构的抗弯刚度，对接触网的平衡构型进行计算；类似的采用抛物线单元[8]，对简单链型悬挂接触网进行离散，然后进行平衡构型的计算；意大利学者Bruni等[9]总结了全球10家致力于研究弓网关系的团队针对一个标准的弓网模型进行的动力学计算，里面也包含了每个单位对于接触网的找形方法。不难看出，从现有的接触网建模技术来看[10−19]，针对接触网找形问题，国内外学者都是针对接触网的有限元模型，对其平衡状态下的系统广义坐标(或者还包括单元初始长度)进行计算，而计算精度大多都依赖于接触网离散模型的自由度，这也是传统有限单元基于小变形假设所带来得一个不可避免的问题。

多柔体动力学理论在经过几十年的发展后，在工程领域取得广泛应用。其中由Shabana[20]提出的绝对节点坐标法(ANCF)，由于系统的广义坐标直接采用节点的绝对位置矢量和绝对梯度矢量，由于将节点广义坐标中的梯度矢量替换了传统有限单元的小转动矢量，所以没有转动幅度的限制，因此在描述铁路接触网系统、电梯升降绳以及流体建模等[19,21−23]众多柔性大变形问题具有很多优势。而基于连续介质力学[20]的ANCF在描述其有限单元时，其单元的二阶连续性能够使系统有限元模型的维度大大降低，因此能更好地进行有限元模型的动力学计算。

本文采用基于ANCF描述的多柔体动力学理论建立接触网的有限元模型，接触网由两节点的ANCF索单元进行离散，模型考虑了接触网中的弯曲变形和轴向变形，考虑到接触网中几乎不存在扭转变形，所以忽略了接触网中的扭转变形。由于接触网系统是一种张紧的工程结构，在找形分析中使单元的初始长度变化来满足接触网中的张力要求。并根据拉式方程推导了在索单元长度未知情况下的动力学模型，针对平衡构型的计算问题，对动力学模型进行退化处理。在引入未知单元长度的模型后，需要根据接触网中的张力条件进行补充方程，得到的接触网平衡构型的计算模型整体计算维度非常小，仅与每跨中的吊弦数量相关。最后通过3个数值算例，分别从单元稳态构型的计算、简单接触网稳态构型的计算以及工程应用中的接触网稳态构型的计算上验证了本文找形方法的可行性。

1 ANCF基本动力学理论

1.1 单元长度未知的ANCF索单元理论

两节点的ANCF索单元如图1所示，其中心线上任意点P的绝对位置矢量可表示为：

图1 两节点索单元变形前后示意图Fig.1 Deformed and un-deformed two-nodes cable element

式中：L为单元初始长度；l为任意点P在中心线上的弧长坐标，满足：

单元初始长度L为时变的未知量，定义单元初始长度L为：

单元的动能可表示为[24]：

单元势能由应变能和重力势能组成：

1.2 索单元动力学方程

在引入约束条件下，由带约束乘子的第一类拉格朗日方程建立索单元动力学方程：

由式(15)可推导得变长度索动力学模型的另一种表达形式：

1.3 动力学方程求解方法

式(16)是一组典型的微分-代数方程(DAEs)，本文采用经典的增广法[26]进行求解，为了使解趋于稳定，引入闭环系统：

对于多单元系统，采用传统有限元共用节点组装的方法，根据式(20)，即可求出系统的广义坐标q、广义速度q˙和广义加速度q¨。

2 动力学模型退化的找形方法

2.1 平衡状态下的动力学模型

第1节推导了单元长度未知情况下的索单元动力学模型，而平衡状态就是动力学模型中一个特殊情况。对应的条件是与时间相关的一阶、二阶导数均为0，即：

而式(23)在找形问题中需要求解的变量只有单元的广义坐标，所以可以根据式(23)得到：

式中：DOFs为单元自由度，对于本模型的ANCF索单元，单元自由度为12；另外定义：

值得注意的是，在式(25)中，未知的只有由动能附加项引起的广义力附加项，动能附加项表达式为：

而在应用拉式方程建立索单元动力学模型时，动能附加项参与计算时，会引起两项广义力附加项：

但当单元处于平衡状态时，由于式(28)和式(29)中均含有与时间相关的导数项，根据式(22)可知：

由此可知，由动能附加项引起的广义力附加项：

所以，式(25)可改写为：

2.2 引入张力条件的找形方法

根据式(24)可知，由动力学模型退化得到一组独立的非线性方程组的维度只是单元自由度数，但还需要求解的还有单元初始长度，所以需要再添加一个独立的方程，对于张紧的单元，可根据张力条件补充方程：

式中，T为给定的张力条件。联立式(24)和式(32)可得到一组非线性方程组：

3 接触网找形问题

3.1 接触网有限元模型

简单链型悬挂接触网如图2所示，以3跨简单链型悬挂接触网为例，吊弦作为连接承力索和接触线的结构，使得接触网成为一个弹性均匀分布的系统，同时，定位装置将接触线和承力索并不是沿轨道呈直线分布，而是呈“之”形分布，如图2中俯视图所示，从而使受电弓滑板能够均匀磨耗。另外，接触线和承力索的下锚处安装了棘轮补偿装置，使得接触线和承力索一直处于张紧状态，在图中用Tcw、Tmw表示。

两个吊弦之间的接触线/承力索、接触线/承力索与定位装置之间，还有定位器均采用一个变长度索单元。接触网系统离散为N个变长度ANCF索单元：

式中：ns为跨数；nd为一跨内的吊弦数量。

3.2 接触网中的约束处理

为了简化接触网中的约束，在接触网有限元模型中的连接方式均采用球铰约束[19]，在绝对坐标系下的3单元球铰连接如图3所示。其数学模型表示为：{

图3 3个单元球铰连接示意图Fig.3 Schematic diagram of spherical hinge connection of 3 elements

在接触网中，吊弦与接触线单元、吊弦与承力索单元均采用这种约束。

另外，接触网系统的边界条件也可简化为接触网系统的边界单元铰接在固定点：

式中：rb表示边界处的位置矢量，即包含承力索和接触线的两端点的位置矢量；rfix表示边界点约束的位置矢量。

3.3 接触网的平衡构型

接触网的找形问题就是确定接触网有限元模型在平衡状态下系统的广义坐标，由于接触网采用变长度索单元进行有限元划分，平衡构型还需要确定的是单元的初始长度，用x表示为平衡构型：

针对接触网的找形问题，不仅只是确定接触网系统在边界条件下的平衡构型，吊弦还有预分配的约束条件：

式中：rd为吊弦下端的位置矢量；rdesign为吊弦下端点的预设条件。

3.4 接触网找形应用

根据吊弦预配要求，那么接触线就可从接触网系统中分离出来单独求解，再根据接触网的离散方式，那么每个接触线单元就可处理为两端铰接在已知点的张紧悬链线模型

3.4.1接触线找形分析

将接触线的每个单元处理为两端悬挂在已知点的悬链线模型，而对于接触网中的张力主要分量是水平分量，所以这里简化处理为：

式中：Tcw为接触线的张力条件：上标“x”为张力的水平分量。

利用水平张力作为张力补充条件，那么应用式(33)对每个接触线单元进行找形分析时，则需要改写为：

式中，Ec、Ac分别为接触线的弹性模量和截面积。

完成接触线的平衡构型求解之后，根据相应接触线单元的稳态构型，提取在定位点处的合力作为定位器单元的张力补充条件，可以得到定位器稳态构型的求解方程：

式中：i=[0,1,0]T；EcAc分别为定位器的弹性模量和截面积；Ts为反求得到的定位器所受张力。

3.4.2吊弦承力索系统的找形应用

在确定了接触线的平衡构型后，与定位器类似的，每根吊弦单元的下端张力也可由相应的接触线单元的平衡构型条件反求出：

式中，j=[0,1,0]T，补充了吊弦的张力方程，对于承力索，张力补充方程与接触线的张力处理方式一致：式中：Tmw为承力索的张力条件；上标“x”表示承力索张力的水平分量。

则对于承力索与吊弦组成的系统，其动力学方程由经典的有限元填装方法可得到，相应的，对于由多单元组成的承力索吊弦系统，需要补充相应单元数的张力条件。利用式(42)可补充吊弦单元的张力条件，再由承力索的水平张力条件，那么可以得到吊弦承力索系统的平衡构型求解方程：

式中：n=nd+1；Em、Am、Ed、Ad分别为承力索和吊弦的弹性模量和截面积。张力补充方程的前n=nd+1个方程由各个承力索单元提供，而最后nd个方程由吊弦补充方程提供；对应了每个单元的补充方程。

3.4.3小结

总结前两部分的找形应用过程，可以看出本文的方法是基于分模法[27]的思路，先将接触网进行分模处理，将接触线的平衡构型求出后，利用其平衡构型求出相应吊弦点处的吊弦力，以此补充相应缺失的张力控制方程，最后完成整个接触网的找形分析；整个流程如图4所示。

图4 接触网找形流程示意图Fig.4 Schematic diagram of form-finding steps of the railway catenary

4 算例验证

本节用3个算例对本文提出的找形方法进行验证，包括悬链线模型、两吊弦简单接触网模型和EN50318中的简单链型悬挂接触网模型。本文在针对接触网进行找形分析时，针对接触线的找形时是将每个接触线单元视作一个悬链线模型(包含等高与不等高悬链线模型)，所以引入了悬链线模型进行初步验证本文提出方法的可行性。

4.1 悬链线模型

悬链线示意图如图5所示，两端铰接在固定点处，同时端部还满足张力控制要求。悬链线的密度为ρ=7800 kg/m3，截面直径D=0.01 m，弹性模量为E=2.01 GPa，端部张力控制为T=15 kN。

图5 悬链线示意图Fig.5 Schematic diagram of catenary

对该模型进行找形分析，跨距d=20 m，悬挂高度差分别取h=0 m、0.1 m、0.2 m、0.3 m、0.4 m、0.5 m。找形结果如下图6所示。

图6 不同悬挂高度的悬链线找形结果Fig.6 Form-finding resultsof catenary at different suspension height

随着悬挂高度的增大，悬链线稳态构型的最大静挠度逐渐变小，稳态悬链线构型也逐渐逼近其对应悬挂点之间的连线。通过提取施加张力段的竖直张力分量分别为54.88 N、129.88 N、204.87 N、279.85 N、354.81 N、429.75 N，而张力的主要分力是水平分量，所以在端部或者悬链线中任意点偏置水平位置角度是逐渐增大的。

为了验证找形结果的正确性，建立相应的动力学模型，再代入找形结果进行动力学仿真，不进行重复性的工作，只做图6中高度为h=0 m,0.3 m的验证，提取端部张力控制时程曲线如图7所示。由张力的时程曲线可以看出，在使用本文中的方法进行找形的结果在代入到动力学模型中时，悬链线的端部张力的几乎是不变的，尽管图7(b)不等高悬链线的端部张力最后的稳态值与目标张力在数值上有一点差异(20 N)，但是其相对误差仅为0.13%，这完全可以认为模型是稳定的。

图7 悬链线端部张力时程曲线Fig.7 Time history of tension at the tip of catenary

4.2 两吊弦简单接触网模型

文献[28]提出的两吊弦单跨简单接触网如图8所示，用于验证弓网接触的动力学响应。材料参数由表1给出，几何参数由表2给出。作为一个经典的标准模型(benchmark catenary)，文献[6,29]也对该模型进行了找形分析，也不乏基于绝对节点坐标法的找形方法[29]。

图8 两吊弦接触网示意图Fig.8 Schematic diagram of two dropperscatenary

表1接触网材料参数Table 1 Material parameters for the benchmark catenary

应用动力学退化方法进行找形分析，得到接触网的稳态构型如图9所示。

图9 两吊弦接触网找形结果Fig.9 Form-finding result of the benchmark catenary

提取吊弦点处的吊弦力均为75.69 N，吊弦长度均为0.9552 m；从仿真结果也体现出了结构的对称性；另外对于吊弦长度的计算结果来看，与文献中给出的计算结果分别为：0.9579 m[6]、0.95 m[28]及0.9540 m[29]非常接近，相对误差分别为0.55%、0.13%、0.28%。而在文献[6]的模型中，接触网的稳态构型是通过悬链线方程推导得到的解析形式，并没有考虑接触网的弯曲和轴向变形，但其他两个模型以及本文中的模型均是将接触网离散为有限单元模型，同时考虑了单元的轴向刚度和抗弯刚度，这说明针对该模型的轴向和弯曲变形在找形分析中是可以忽略的。另外值得注意的是，本文中接触网的有限元模型仅采用8个ANCF索单元，但应用本文的找形方法得到的结果精度却没有降低，这是由于ANCF单元的插值函数是基于三次多项式插值推导得到的，具有二阶连续性。

表2接触网几何参数Table 2 Geometrical parameters for the benchmark catenary

同样的，将求解到的平衡构型代入到相应的动力学模型中进行验证，提取控制端张力时程曲线如下：

由端部张力时程曲线图10(a)接触线端部和图10(b)承力索端部，可以看出，由本文提出的找形方法得到的平衡构型代入到动力学模型中，动力学模型几乎是不震荡的，端部的张力的稳态值与设计张力的相对误差分别为：接触线0.0023%，承力索0.067%。可以看出，对于该模型，本文提出的找形方法能够精确的求解到其平衡构型。

图10 接触网端部张力时程曲线Fig.10 Time history of tension of the benchmark catenary

4.3 铁路接触网系统

在工程应用上，本文采用欧洲标准EN50318[30]中的简单链型悬挂接触网为例，分别进行单跨和多跨接触网找形分析。单跨接触网如图11所示。

图11 简单链型悬挂接触网示意图Fig.11 Schematic diagram of a samplecatenary in EN50318

这个模型也在文献[9]中作为一个标准模型进行研究，有全球10个致力于弓网研究的团队针对本模型进行找形分析和弓网动力学性能评估。这种简单链型悬挂结构应用在法国的LN2和意大利的C270接触网系统中。接触网的材料参数如表3所示，几何参数如表4所示；另外由于存在吊弦的预分配条件，吊弦的预分配位置如表5所示。

表5 吊弦预配参数Table 5 Pre-allocated parametersfor droppers

应用本文的动力学方程退化方法进行找形分析，得到单跨接触网的平衡构型如图12所示。

图12 简单链型悬挂接触网平衡构型Fig.12 Equilibrium configuration of thesample catenary

10家单位都做过该模型的找形研究，由于方法的不同和模型约束考虑的不同，每家单位得到的结果都有差异，但大体上近似[9]。这里采用西班牙的一家单位PACDIN[19]的结果进行对比，提取接触线的平衡构型对比如图13所示。

表3接触网材料参数Table 3 Material parameters for the EN50318 catenary

表4接触网几何参数Table 4 Geometrical parametersfor the EN50318 catenary

图13 接触线平衡构型的对比结果Fig.13 Comparison of equilibrium configuration of contact wire

吊弦是接触网的重要部件，对于吊弦长度和吊弦力的计算，也与参考文献进行对比，如表6所示。

从表6可以看出，在量化之后的找形结果对比中，对于吊弦长度的计算，两种方法得到的结果相对误差不大于1%，虽然吊弦力的计算上相对

表6 吊弦长度与吊弦力计算结果对比Table 6 Comparison of the dropper's length and force by different methods

误差较大，但也没有超过2%。这也说明了本文找形方法的正确性，但本文中接触网有限元模型的自由度却非常小，相比于文献[9]中10家单位的接触网模型的自由度具有极大的优势。

除此之外，从动力学的角度进行验证，将得到的平衡构型代入到动力学模型中，验证模型是否震荡，提取得到的接触线和承力索的端部张力如图14所示。

图14 简单链型悬挂接触网端部张力时程曲线Fig.14 Time history of tension at thesample catenary’s tips

通过提取接触线张力时程图14(a)和承力索端部张力时程图14(b)曲线可以看出，动力学模型在该初始条件下是基本处于稳态的，虽然一开始会有轻微的震荡，但是可以看出张力的振幅并不是很大(20 N左右)，并且最终趋于稳定，稳定的端部张力数值与控制张力的相对误差均小于1%。通过验证动力学模型是否震荡，更加直接地证明了平衡构型的正确性，进一步证明了本文提出的找形方法的可行性。

最后，对于该模型，进行多跨建模，需要考虑定位器的影响，得到找形结果如图15所示。

图15 10跨简单链型悬挂接触网平衡构型Fig.15 Equilibrium configuration of the sample catenary with ten spans

同样的，建立多跨接触网动力学模型，将找形结果代入动力学模型中作为仿真初始条件，提取接触线和承力索的端部张力时程曲线如图16所示。

与单跨接触网的验证结果一致，根据接触线张力时程图16(a)和承力索的端部张力时程图16(b)曲线可以说明系统是基本处于稳定的，开始震荡的赋值相比也设计张力来说非常小,最后收敛的张力数值与设计张力的相对误差不超过1%。

5 结论

本文提出了一种基于动力学退化模型的接触网找形方法，根据平衡条件将未知单元长度的索单元动力学模型退化，再补充求解单元的初始长度的张力控制方程，推导了张紧索单元的找形方程。基于分模法的思想，将接触线的平衡构型预先求出，进而得到吊弦单元的张力补充方程，再根据有限元填装技术，推导了承力索吊弦系统平衡构型的计算方程，最后根据三个算例验证，可以得到以下结论：

(1)基于ANCF理论和未知单元初始长度条件，根据平衡条件对动力学模型进行退化处理，推导了张紧单元的静态构型的计算模型；该模型适用于长大张紧结构的静态构型的计算，如铁路接触网，输电线等。

(2)应用该模型进行EN50318中的简单链型悬挂接触网的找形分析，对比Benchmark中的某单位的计算结果，得到的接触线的静态构型几乎一致，并且在吊弦长度的计算上，相对误差均小于1%，吊弦力的计算上，相对误差小于2%。满足工程应用的要求。并且通过动力学模型的验证，得到的张力控制的相对误差小于1‰，验证了本文方法的准确性。

(3)由于模型中能够考虑稳态构型下的吊弦预分配条件，所以该模型适用于所有简单链型悬挂接触网的找形计算，能够指导接触网的设计和施工。

(4)由于该找形方法是基于动力学模型推导得到的，所以找形结果可应用于接触网系统的动力学计算，动力学模型能够快速地趋于稳定。