周期性三维轨道结构弹性波耦合与转换

2022-04-01王树国赵才友刘伟斌

铁道学报 2022年1期

易强，王平，王树国，赵才友，刘伟斌

(1.中国铁道科学研究院集团有限公司铁道建筑研究所, 北京 100081；2.西南交通大学高速铁路线路工程教育部重点实验室，四川成都 610031)

铁路轨道结构振动噪声与结构中弹性波传播特性密切相关。轨道结构具有明显周期特征，既有研究表明弹性波在周期结构中的传播具有带隙特性，即在带隙频率范围内，周期结构中的弹性波将无法自由传播而呈现衰减特性[1]。周期结构中的弹性波带隙特性可为结构振动噪声控制提供新的研究思路和方法，因此铁路轨道结构弹性波传播特性研究得到了广泛的关注。

一般区间路段铁路轨道结构具有周期性以及近似无限长特征，目前常采用有限元方法建立轨道结构模型，为了消除边界条件影响，模型长度通常较长，轨道结构自由度较多，为弹性波分析带来困难[2]。基于周期特征，可将轨道结构弹性波分析转移至一个基本元胞内进行，从而大幅度降低轨道结构自由度，提高计算效率。在轨道结构弹性波传播分析方面，Mead等[1]率先开展了基于行波分析法的周期结构研究，将轨道结构简化为无限长周期梁结构，提出传递矩阵法、空间谐波法等计算方法，并首次给出无限周期梁结构的弹性波带隙特性。Heckl等[3]在Mead研究基础上，在模型中增加轨枕和道砟的影响，将轨枕等效为质量块，道砟层等效为弹簧，研究周期性支承Timoshenko梁中压缩波、扭转波、弯曲波的传播特性。Abe等[4]将铁路轨道考虑为离散支撑约束结构，得到轨道结构弹性波频散曲线，并通过优化轨枕间距以降低钢轨在pinned-pinned频率附近的响应幅值。Sheng等[5]采用2.5D有限元方法建立轨道结构的特征方程，得到轨道结构的振动传递系数，并分析了周期性轨道结构带隙边界频率处对应的结构模态。为了分析高频范围内轨道结构弹波传播特性，吴天行等[6]采用多梁模型模拟钢轨横向振动，其中轨头和轨底采用Timoshenko梁描述，轨道结构采用多个有限长梁进行模拟，与有限元方法结果较为一致。Thompson[7-8]采用有限元方法建立10 mm长度短钢轨模型，考虑截面形状，结合周期结构理论研究了中高频弹性波在自由钢轨中的传播特征，提取钢轨弹性波波数，分析不同类型弹性波在钢轨中的传播衰减规律，并开展现场试验验证。Gry[9]将钢轨截面进行网格划分，采用有限的截面模态模拟钢轨截面变形，进一步采用传递矩阵法建立相邻两跨钢轨之间的传递矩阵，计算周期支撑钢轨中的频散曲线与频率响应函数。

在轨道结构动力响应方面，基于结构周期特征可实现轨道结构动力响应的快速计算。马龙祥等[10]采用周期傅里叶方法建立车辆-浮置板轨道结构动力响应计算模型，在计算中只需考虑一个浮置板长度范围内的轨道结构自由度。Degrande等[11]和Chebli等[12]利用Floquet变换结合有限元、边界元方法建立轨道-土体耦合模型，计算简谐荷载作用下无限周期结构系统响应特征。Sheng等[13-14]采用广义傅里叶变换方法和2.5D有限元方法计算移动简谐荷载作用下周期性轨道结构动态响应。上述研究主要以平面半轨道结构模型为主，三维轨道结构中弹性波传播特性尚待进一步研究。

对于三维轨道结构的既有研究以传统有限元方法为主，主要关注有限长轨道结构动力响应，无法直接应用于无限长周期轨道结构的弹性波分析及结构响应计算。近年来，国内外学者针对周期结构的分析提出传递矩阵法、平面波展开法、集中质量法等，但这些方法均有一定的限制性，如传递矩阵法容易出现矩阵病态问题[15]，且难以分析复杂三维结构模型。波有限元方法(WFE)的提出为复杂结构中弹性波分析提供了高效手段，波有限元方法可以由传统有限元方法建立元胞有限单元模型，在元胞边界施加周期边界条件以研究周期结构弹性波的传播特性[16]。波有限元方法中可建立详细的元胞三维模型，从而可应用于轨道结构弹性波传播研究与中高频振动响应分析。

对于典型的周期性轨道结构，其弹性波传播特性已有一定的研究，但主要结论均基于平面半轨道模型[17-18]。若将轨道结构简化为平面半轨道模型，将忽略相邻钢轨之间的相互作用，且轨枕模态也无法被考虑，导致无法准确分析中高频范围内结构弹性波传播。轨道结构振动衰减规律与相邻两根钢轨之间的对称模态与反对称模态密切相关，半轨道模型将过高估计部分频段内振动的衰减[19]。因此，有必要建立三维轨道结构模型，以更加准确分析弹性波在轨道结构中的传播规律。为了分析周期性三维轨道结构中弹性波耦合与转换规律，本文以有砟轨道结构为例，采用波有限元方法建立其弹性波传播分析模型，研究不同类型弹性波之间的耦合与转换特征并阐明三维轨道结构中弹性波带隙行为。

1 理论模型

采用波有限元方法建立轨道结构基本元胞模型，并提取其质量、刚度、阻尼矩阵。以有砟轨道结构为例，元胞长度等于轨枕间距，元胞结构为枕间距内轨道结构，如图1所示，周期性边界条件采用虚线表示。

图1 有砟轨道结构元胞模型

元胞结构动态矩阵方程可表示为

(K+iωC-ω2M)U=F

(1)

式中：D=K+iωC-ω2M为一个基本元胞的动刚度矩阵；K为刚度矩阵；M为质量矩阵；C为阻尼矩阵；U为广义位移向量；F为广义力向量。

根据节点位置可将U、F分为三个部分：左边界节点，内部节点与右边界节点，分别采用下标L，I和R表示，且左右节点自由度均为nc个，如图1(c)所示。力向量与位移向量之间通过动刚度矩阵联系，即

(2)

对于结构中自由波传播，即无外界作用力下，消除内部自由度，在元胞边界处力和位移满足条件

(3)

式中：“～”表示消除内部自由度后的动刚度矩阵。

(4)

根据Bloch定理可知，元胞左右两侧状态之间满足[20]

UR=eiκlUL

(5)

FR=-eiκlFL

(6)

式中：κ为Bloch波数；l为元胞长度。

结合相邻元胞之间的位移和力连续条件得

UR,ne=UL,ne+1=eiκlUL,ne

(7)

FR,ne=-FL,ne+1=-eiκlFL,ne

(8)

式中：ne为元胞编号。

因此，可构建广义特征值矩阵为

(9)

对于频率较高或结构尺寸较大的情况，以上矩阵可能出现病态问题，难以求解准确结果。因此，可根据相邻元胞向量关系直接形成标准的多项式特征问题，从而求解波数及与之对应的特征向量为

(10)

在元胞左节点位置，广义力向量与位移向量之间满足

(11)

在右节点位置，广义力向量与位移向量之间满足

(12)

根据式(10)～式(12)即可求解由有限单元建立的周期结构弹性波传播频散曲线以及所对应的特征波模态。在给定频率下，结构中可存在nc对特征波，每一对中包含了两个幅值相同但传播方向相反的弹性波。对于第j种弹性波，其波数为κj，将其正向传播的弹性波表示为[κpj,ψpj,Θpj]，负向传播的弹性波表示为[κnj,ψnj,Θnj]，其中ψ和Θ分别为力向量和位移向量，下标p和n分别表示正向传播和负向传播。

在求得所有的弹性波波数以及特征向量后，结合弹性波叠加方法即可求解外界荷载作用下有限或无限周期结构响应解[16]。元胞边界位置处由自由波传播产生的位移和力响应可表示为

(13)

式中：αj为第j类弹性波所对应的广义波坐标。

对于半无限长周期结构，若在端点施加激励，则只存在向左(左-半无限周期结构)或向右(右-半无限周期结构)传播的弹性波，式(13)可进一步进行简化。

对于左-半无限周期结构，由于只能激励产生负向传播的弹性波，因此端点处结构响应可写为

(14)

通过矩阵方式可将式(14)改写成

uR=Θ-α-FR=ψ-α-

(15)

式中：下标“-”表示负方向传播弹性波所对应的特征向量和广义波坐标。

因此，左-半无限周期结构右端点处动刚度矩阵为

D-R=FRuR-1=ψ-Θ--1

(16)

对于右-半无限周期结构，由于只能激励产生正向传播的弹性波，因此左端点处结构响应为

(17)

同理，右-半无限周期结构左端点处动刚度矩阵为

D+L=FLuL-1=ψ+Θ+-1

(18)

除此之外，实际轨道结构可能存在局部缺陷，如局部刚度突变或结构局部伤损，需要考虑弹性波在局部有限结构中的传播与响应特性。对于有限结构，其结构响应则需要同时考虑两个方向传播的弹性波。当有限结构由多个元胞组成时，可利用周期性将其计算自由度大幅度缩减。有限结构中第m个元胞响应为

(19)

式中：N为有限结构包含的元胞个数。

则有限结构左右边界处力和位移为

(20)

(21)

此时，该有限结构动刚度矩阵为

(22)

值得注意的是，此时结构动刚度矩阵只包含结构左右边界处节点向量。因此，计算自由度为2nc，计算过程中的最大自由度为一个元胞自由度(用于求解波数与特征向量)，利用该方法可大大降低结构计算自由度，提升计算效率。

最后对于激励源所处的元胞，由于元胞内部自由度存在外界激励，根据式(2)，该受载元胞左右端点处力向量为

(23)

根据以上推导，结合不同类型周期结构边界处位移、力的连续条件即可组建任意有限-无限周期结构方程进行求解。如无限周期性轨道结构在某一位置受外界荷载时，整体结构可视为由左-半无限周期结构+激励单元+右-半无限周期结构组成。受载元胞左右边界处与半无限周期结构端点处满足连续条件

(24)

式中：F-inf,R、u-inf,R分别为左-半无限结构最右端力、位移向量；F0,L、u0,L分别为有限结构左端力、位移向量；F+inf,L、u+inf,L分别为右-半无限结构最左端力、位移向量；F0,R、u0,R分别为有限结构右端力、位移向量。

结合式(16)、式(19)、式(23)、式(24)可组建系统矩阵

(25)

式中：Fext为作用在结构上的外界荷载。

根据式(25)即可求解无限轨道结构动态响应。值得一提的是，上述方法中并不要求各区域周期结构基本元胞一致，因此可应用于更复杂工况分析。此外，由于结构计算自由度大量减小，可建立详细的元胞结构有限元模型。

2 弹性波传播特性分析

为了研究弹性波在三维轨道结构中传播特性，以有砟轨道结构为研究对象建立其波有限元分析模型，如图1所示。元胞模型中包含两根钢轨与一根轨枕，钢轨、轨枕均采用Timoshenko梁单元进行建模，网格尺寸为0.05 m。扣件与枕下支撑采用弹簧单元模拟并考虑六个方向刚度，重点分析垂向弯曲波在轨道结构中的传播规律。

轨道结构、扣件及道床刚度计算参数见表1、表2。

表1 有砟轨道结构参数

表2 扣件及道床刚度参数取值

2.1 与平面模型对比

首先进行三维模型与平面半轨道模型对比分析。在平面半轨道结构模型中，轨枕只能简化为刚体且只能考虑一根钢轨。为了分析两种模型的差异，三维轨道结构中也将轨枕考虑为刚体，即只考虑轨枕的垂向运动及转动，但三维模型中考虑了两根钢轨。为了更加直观描述结构中可传播的弹性波和衰减波，频散曲线实部给出了虚部为0的正波数，而频散曲线虚部则给出了波数虚部的最小绝对值。

两种模型中垂向弯曲波(即结构发生z方向变形)频散曲线对比结果如图2所示，图2(a)和图2(b)给出了可自由传播的弹性波频散曲线，而图2(c)和图2(d)则给出了波数虚部的最小绝对值随频率的变化，即描述弹性波在结构中的衰减特性。由频散曲线可知，垂向弯曲波在1 500 Hz频率范围内存在三个带隙[18]，如图2(c)中阴影部分所示，其带隙边界处特征模态由图2(a)中Mode A — Mode E给出，这些特征波模态均对应轨道结构中的驻波模态，且已在参考文献[17]中给出。相对于半轨道模型，三维模型中元胞边界自由度增加一倍，因此结构中可传播的弹性波类型增加一倍，对应波数个数增加一倍，在给定频率下，存在4种正向传播的垂向弯曲波。从计算结果可以看出，半轨道模型中频散曲线与三维轨道结构频散曲线部分吻合，但三维模型中存在另外一条频散曲线，对应结构中存在非对称波模态，即左右两根钢轨反相运动。由于此时将轨枕考虑为刚体，对称模态和非对称模态对应的频散曲线差异不大，如图2(a)所示。

图2 平面半轨道与三维轨道结构中垂向弯曲波频散曲线

图3给出了三维轨道结构中Mode E对应的波模态。在任意一条频散曲线起始或截止的位置附近均存在对称和反对称驻波模态，二者特征频率接近。在高频段内以钢轨变形为主，不受轨枕影响，因此对称模态与反对称模态频率一致。由图2(c)和图2(d)可以看出，虽然平面模型与三维模型得到的频散曲线结果基本一致，但三维轨道模型计算得到带隙内弹性波衰减能力明显降低。此外三维模型中包含了更多类型波模态，在外界激励下产生更多可传播弹性波。因此，有必要进一步研究三维轨道结构弹性波传播及结构响应特征。

图3 三维轨道结构中的对称与反对称模态

2.2 弹性波耦合与转换

考虑轨枕变形，将轨枕简化为Timoshenko梁建立三维有砟轨道结构模型。计算得到所有可传播弹性波频散曲线如图4所示，同样只显示了具有纯实部的频散曲线。模型中相邻元胞通过12个自由度进行耦合，因此在给定频率下，可得到12对波数。此时，结构中包含4种类型弹性波：垂向弯曲波、横向弯曲波、纵波和扭转波。从图中可以看出，此时频散曲线中出现多条平直能带，这些能带由轨枕固有模态产生。由于不同类型弹性波之间发生相互耦合与转换，从而难以根据频散曲线直观识别不同类型弹性波传播特征。

为了进一步分析不同类型弹性波在结构中的耦合与传播规律，采用波模态分离方法实现不同类型弹性波分离。在波模态分析中，采用模态置信准则中MAC值评价两种类型弹性波的相关性[21]，即

(26)

式中：ψA和ψX分别为模态A和模态X的特征向量。

MAC值表明这两个模态之间的相关性，其值处于0和1之间。当MAC>0.9时，说明这两个模态一致；而当MAC<0.6时，说明两个模态相关性较弱[22]。

图4 三维轨道结构频散曲线

为提取与垂向弯曲波相关的频散曲线，以波数κ=0时垂向弯曲波模态为初始参考进行波模态分离，并在波模态分离过程中随着波数增大不断更新参考波模态，得到与垂向弯曲波相关的频散曲线如图5所示。从分离后的频散曲线中可以看出，轨道结构中扭转波与弯曲波相互耦合。这是由于轨枕的作用，钢轨扭转带动轨枕产生垂向弯曲变形，从而使另一钢轨产生垂向弯曲波，即在三维轨道结构中存在弯曲波与扭转波的耦合与转换。值得注意的是，波模态转换是随着波数增大而逐渐形成的，因而在同一条频散曲线中可存在不同类型的波模态。由于不同类型弹性波之间存在耦合与转换，这为准确识别弹性波带隙范围带来困难。此外，在耦合转换位置的弹性波波模态及结构响应特性需进一步讨论。

图5 分离后得到与垂向弯曲波相关的频散曲线

从图5中可知，垂向弯曲波与扭转波频散曲线第一个交叉区域位于100～300 Hz。为了分析不同类型弹性波的耦合转换，提取300 Hz频率范围内与垂向弯曲波相关频散曲线及关键位置特征波模态如图6所示。若只从频散曲线来看，该范围内垂向弯曲波只存在一个0～127 Hz低频带隙。但由于钢轨扭转波与弯曲波的耦合作用，同一条频散曲线中可能存在不同类型的波模态。如图6中200 Hz附近黑色频散曲线，当κl=0时，结构波模态Mode S1表现为钢轨扭转，轨枕弯曲，此时垂向变形极小，表明该波模态以扭转波为主；而当κl=2.8时，结构波模态Mode S2为明显的钢轨垂向弯曲，而钢轨扭转位移较小。对于该条频散曲线，随着频率增加，其波模态由扭转波向弯曲波转变。而对于Mode S3和S4所在的频散曲线，虽然该频散曲线与垂向弯曲波相关，但在此处仍以扭转波为主。

图6 垂向弯曲波频散曲线及波模态转换

为了分析不同类型弹性波的耦合与转换过程，以Mode S4所在频散曲线为例，计算其群速度vg，得到如图7所示的结果。从发生波模态转换的位置可以看出，波模态转换均在群速度发生突变的位置。发生波模态转换时，群速度迅速降低，发生突变，如图7(b)中的波模态Mode T0，此时钢轨扭转变形与弯曲变形相当。在发生波模态转转换前，弹性波传播以扭转波为主，如Mode T1；在发生波模态转换后，则以弯曲波传播为主，扭转波成分逐渐降低，如Mode T2。发生波模态转换前，以扭转波为主，扭转波群速度较高，而波模态转换后以弯曲波为主，群速度明显显著降低，因而在群速度突变位置可视为波模态转换位置。发生转换的位置在图6(a)中以空心圆形式标出。

图7 弹性波耦合与波模态转换

结合波模态分析可以看出，轨枕发生弯曲共振，使得两根钢轨之间弹性波发生耦合，垂向弯曲波可通过轨枕弯曲运动在另一钢轨形成扭转波。因此，三维轨道结构中发生弹性波耦合与转换的根本机制是轨枕的局域共振作用。

通过以上分析可知，只有确定波模态的转换位置后，才能确定弹性波带隙特征，且在发生波模态转换位置(图6)出现带隙截止频率。在300 Hz频率范围内存在三个垂向弯曲波带隙：0～127.2 Hz、181.0～197.7 Hz、200.0～246.4 Hz。随着频率的升高，扭转波与弯曲波耦合作用逐渐降低，在1 085～1 129 Hz存在第四个垂向弯曲波带隙。

相对于平面模型，三维轨道结构中由于存在弹性波之间的耦合转换，弹性波带隙被抑制。此外，周期结构中局域共振单元的存在会形成局域共振带隙，即局域振子自身共振与结构中行波相互作用，从而抑制其传播[17]。但对于三维轨道结构，存在对称波与反对称波，轨枕局域共振只能与其中一类弹性相互作用形成局域共振，但另一类波仍然可自由传播，如图8所示。

图8 局域共振与钢轨长波形成的单一耦合

在600 Hz附近，存在轨枕弯曲共振，此时轨枕弯曲共振模态为对称模态，可以与钢轨中对称波相互作用形成局域共振带隙，但不影响非对称波的传播。单一的局域共振模态将无法在三维轨道结构中形成完全局域共振带隙，因而该不完全带隙可称之为对称型局域共振带隙，只有在对称荷载激励下才能表现为带隙。只有当轨枕局域共振模态与钢轨中的对称波和反对称波同时产生耦合作用，才能形成完全的局域共振带隙。

3 振动响应及传递规律

为了进一步分析和验证弹性波传播特性，基于波有限元方法(WFE)计算无限长三维轨道结构振动响应。与此同时，建立由70个元胞组成的轨道结构有限元模型(FE)，在扣件上方一侧钢轨施加垂向单位简谐荷载，加载位置处钢轨位移响应见图9。

图9 有砟轨道结构响应(单侧非对称荷载)

由图9可见，两者吻合良好，因此可以验证上述理论分析的正确性。除了利用波有限元方法可大大降低计算时间以外，由于轨道结构具有无限长特征，传统有限元模型难以消除边界弹性波的反射作用，从而使得结构响应在通带频率范围内形成波动，无法得到准确结果，如图9(a)中通带位置。

由于只在单侧钢轨施加荷载，因此可同时激励产生轨道结构中的对称与反对称波模态，从而在带隙边界出现两个峰值，对应结构的对称模态与反对称模态，如图9(b)和图9(c)所示。但当在两根钢轨施加对称荷载时，则只能激励产生对称波模态，从而在带隙边界表现为单一峰值，如图10所示。与此同时，在外界荷载激励下，结构频响函数在波模态转换位置(f=197.7 Hz)出现了峰值。虽然该频率并不对应任何驻波模态，但由于波模式转换位置弹性波群速度降低至谷值，形成近似驻波。

图10 有砟轨道结构响应(对称荷载,WFE)

为了验证弹性波的传播规律，在单侧钢轨处施加垂向单位简谐荷载以模拟非对称激励，在两根钢轨同时施加同相位的垂向单位简谐荷载以模拟对称激励。以激励位置为参考点，提取距激励点6 m处结构响应并计算不同激励下振动传递系数，如图11所示。

图11 有砟轨道结构振动传递规律

从结果可以看出，第一阶垂向弯曲振动带隙为0～127 Hz，而对于181～254 Hz范围内局域共振带隙则较为复杂。由于弯扭耦合的存在，在局域共振带隙内振动传递系数曲线形成多个峰值(如f=188、204 Hz)，这些峰值是由于垂向激励作用使得钢轨产生扭转共振(对应图6(a)中κ=0的情况)，从而增强垂向弹性波沿结构传递能力，但并不会形成通带。轨枕的一阶弯曲共振使得垂向弯曲波与扭转波产生转换，在198～200 Hz区间形成窄通带，通带起始频率对应波模式转换位置，且在该通带范围内无明显的振动传递峰值。根据前文分析可知，该通带由扭转波与弯曲波之间的转换所形成，在结构中无法形成驻波模态，因此振动传递曲线在该通带内无法形成明显的峰值。此外，对称激励与非对称激励可得到不同的带隙范围，对称激励产生的带隙范围更宽。如图11(c)所示，在非对称荷载激励下，可同时产生对称模态峰值与非对称模态峰值，带隙截止频率由254 Hz降低至246 Hz，因而带隙宽度降低，且带隙内对弹性波的衰减能力也大幅度减小。

此外，轨枕对称/非对称型局域共振带隙表现行为也与激励类型密切相关。如600 Hz附近轨枕对称共振模态，该局域共振模态只与结构中对称弹性波产生耦合作用形成对称型局域共振带隙。因此，在非对称荷载激励下，产生的非对称波可自由传播，而在对称激励下，则在614～619 Hz范围内形成局域共振带隙，如图12所示。