基于能量泛函变分原理的无砟轨道垂向振动带隙分析

冯青松，廖宝亮，郭文杰，杨 舟，付景文，陆建飞

(1.华东交通大学 铁路环境振动与噪声教育部工程研究中心，江西 南昌 330013；2.江苏大学 土木工程与力学学院，江苏 镇江 212013)

由基本单元沿线路方向周期性排列而成是铁路轨道结构的基本构成形式。当振动以弹性波形式在轨道结构中传播时，满足Bloch理论和周期性条件，弹性波会形成相应的波带隙(带隙中的波会迅速衰减以致无法传播)。实际应用中，带隙特性对于结构减振或保障结构的可靠性有很大影响。此外，对轨道结构有关参数进行优化设计，可达到优化其带隙频率的位置与宽度的目的，提高其减振效果，从而为轨道结构减振设计提供新思路。因此，研究周期性无砟轨道结构中弹性波传播特性对其振动与噪声控制技术发展具有重要意义。

用于求解周期性结构的带隙特性的方法主要有有限元法、传递矩阵法、平面波展开法等。其中，有限元法适用性最强，一般是通过有限元软件仿真实现带隙计算，但计算过程不同于一般的特征值求解，仅有Comsol[1-3]等少数有限元软件可通过直接设置Floquet周期性边界并扫描波矢实现带隙计算。然而，这类商业软件计算模块尚未完善，以Comsol为例，对于无砟轨道结构(梁-板组合结构)需要采用实体单元建模，计算量较大。也有学者通过自编有限元程序计算周期支撑下钢轨频散特性，但尚未考虑轨道板的影响[4]。

传递矩阵法[5]对于求解简单一维结构带隙特征是十分便捷有效的，如计算考虑周期性弹簧支撑的钢轨结构带隙[6-7]。本文建立的无砟轨道模型同时考虑了钢轨和轨道板，即梁-板组合周期结构，而传递矩阵法难于对轨道板建模，因此该方法并不适用。

平面波展开法[8-11]是用于求解结构带隙的通用方法，但是受限于直接求解微分方程组的分析思路，在处理梁、板连接部位及轨道板边界条件时有一定困难。此外，也有学者提出了一些新的周期结构带隙计算方法，如谱动刚度法[12]、小波法[13]、渐进匹配展开法[14]等，这些方法虽各有特点，但处理梁-板组合模型并不适用。

此外，从振动响应频谱特征判断周期轨道结构中振动波的传播与衰减[15-21]是目前广大研究学者采用的分析其振动特性的思路与方法。但这类研究思路无法从机理上揭示带隙的形成过程，且容易受激励点和测量点位置的影响。

鉴于此，本文以我国CRTSⅢ型无砟轨道结构为研究对象，基于周期结构带隙原理，提出一种计算周期性轨道结构带隙特性的新方法，即基于能量泛函变分原理得到各组分的能量泛函，并将求解微分方程边值问题转化为泛函极值问题，易于解决组合结构的耦合问题。具体来讲，首先，利用平面波级数和布洛赫定理构造钢轨位移场函数，使其能够满足波传播条件和周期边界条件，从而得到钢轨应变能和动能；然后，结合切比雪夫级数处理轨道板位移场，得到其动能和势能，接着通过扣件弹簧建立梁与板位移场联系，并得到扣件弹性势能；最后，将所有各分部能量叠加得到总能量泛函，对所有位移场中未知系数求导后即可得到特征方程，从而求出特征频率及振型函数。本文方法解决了传统解析方法囿于直接求解微分方程组以及难于对轨道板建模的问题，为轨道结构带隙计算提供新的思路。

1 分析与计算方法

1.1 分析模型

建立双层弹性支承轨道结构动力模型，轨道结构简化为由钢轨、扣件、轨道板和CA砂浆层组成的无限周期结构，见图1。扣件与CA砂浆层简化为弹簧单元，为同时考虑钢轨和轨道板的弯曲和剪切性质，钢轨简化为Timoshenko梁单元，轨道板简化为Mindlin板单元。值得一提的是，阻尼对轨道结构振动响应有较为显著的影响，尤其对振动响应频谱曲线中的峰值有明显削弱。然而，阻尼对结构振动特征频率影响相对较小，即对带隙特性影响较弱，因而很多关于轨道带隙特性的理论研究中通常选择忽略阻尼的影响[6-7]。因此，本文研究选择忽略阻尼。

图1 无限周期CRTSⅢ型无砟轨道结构

1.2 计算原理

1.2.1 钢轨和轨道板计算模型

本文基于简化的CRTSⅢ型无砟轨道模型，取一个周期的结构进行建模，模型见图2。

图2 钢轨和轨道板模型

图2分别针对钢轨和轨道板建立两个笛卡尔坐标系roB和xoy。l为单个轨道板的长度；s0为轨道板间缝隙大小，l+s0为钢轨周期，扣件周期为l0，且l0=(l+s0)/9；c为轨道板宽度的一半；y0为钢轨距轨道板外侧边缘的宽度。

本方法考虑系统总能量的组成部分有钢轨应变能与动能、轨道板应变能与动能、板下CA砂浆层势能和梁板连接处扣件弹性势能。

1.2.2 钢轨能量泛函

取一根钢轨进行分析，为全面分析钢轨对带隙特性的影响，将钢轨考虑为Timoshenko梁。根据Bloch定理及平面波级数展开，Timoshenko梁位移场可展开为

(1)

式中：v(r)为垂向位移；β(r)为梁转角，展开项均为(2η+1)项；k为波数；r为梁的横坐标；i为虚数单位；ξ(r)为关于r的试函数行向量；A1,A2为未知系数列向量，且

A1=[v-η,v-η+1,…,v0,…,vη]T

A2=[β-η,β-η+1,…,β0,…,βη]T

则(1)式可写为

(2)

双轨因形变产生的应变能可通过积分得到

(3)

式中：EB为梁的杨氏模量；IB为梁的截面惯性矩；kB为梁的剪切系数；GB为梁的剪切模量；AB为梁的截面面积。

将(1)式代入(3)式，可得

(4)

式中：A0H为A0的共轭转置向量，A0H={A1H，A2H}；KB为梁的刚度矩阵。

双轨动能为

(5)

式中：ρB为梁的线密度；ω为系统圆频率。

将(1)式代入(5)式，可得

(6)

式中：MB为梁的质量矩阵。

1.2.3 轨道板能量泛函

考虑为Mindlin板的轨道板位移场可展开为

(7)

式中：w为板的垂向位移；θx为板沿x向转角；θy为板沿y向转角；Amn、Bmn、Cmn为未知系数；φm(x)、ψn(y)为板的位移试函数，可选择具有任意性的切比雪夫级数[22-23]，其三阶微分可导、四阶微分连续的性质可以满足任意边界条件，且

φ1(x)=1

(8)

(9)

(10)

ψ1(y)=1

(11)

(12)

(13)

式(7)也可写为

(14)

式中：Amn、Bmn、Cmn均为M×N维列向量，且Amn=[A11,A12,…,Amn]T，Bmn=[B11,B12,…,Bmn]T,Cmn=[C11,C12,…,Cmn]T;λ=[φ1(x)ψ1(y),φ1(x)ψ2(y),…,φM(x)ψ1(y),…,φM(x)ψN(y)]=[Q1(x,y),Q2(x,y),…,QMN(x,y)]，对其中任意元素，有φm(x)·ψn(y)=Qa(x,y)，且a=(m-1)N+n，m∈[1,M]，n∈[1,N]。

轨道板因形变而产生弹性体运动，进而产生的应变能可由二重积分得到，可表示为

(15)

轨道板因振动产生的刚体运动而具有动能，其表达式为

(16)

式中：TP为轨道板的动能。

1.2.4 轨道板下CA砂浆层势能

轨道板由其板下CA砂浆层均布支承，同时结合轨道板垂向位移可以通过二重积分得到支承层势能，其表达式为

(17)

式中：U1为支承层势能；ks为板下CA砂浆层平均刚度。

1.2.5 扣件弹性势能

由扣件垂向刚度和扣件处钢轨与轨道板的垂向位移差值可求出扣件弹性势能。模型包含对称的两根钢轨，因此，扣件的总弹性势能U2可表示为

(18)

轨-板位移差值为

(19)

将式(19)代入式(18)，可得

(20)

式中：AH为A的共轭转置向量，AH={A1H,A2H,AmnT,BmnT,CmnT}为未知系数列向量；j为扣件序列编号；j0为单轨扣件个数；kz为扣件垂向刚度；01为 (2η+1)×(2η+1)维零矩阵；02为(2η+1)×MN维零矩阵；03为MN×(2η+1)维零矩阵；04为MN×MN维零矩阵。

1.2.6 总能量泛函

系统总能量泛函Π为

Π=UB+UP+U1+U2-TB-TP

(21)

对未知系数求极值，即

(22)

于是，结构振动问题转化为求解特征值问题，可表示为

(K-ω2M)A=0

(23)

式中：K为系统总刚度矩阵；M为系统总质量矩阵；A为未知的系数列向量；ω为圆频率。

2 算例分析

2.1 收敛性分析

本文将位移进行级数展开，为得到频率结果，计算时需将位移级数进行截断。截断项个数的取值不同，将影响计算精度和效率，因此有必要对不同截断项数下的频率进行收敛性分析。

本节采用控制变量法研究文中涉及的3个截断项数M、N、η的取值对带隙频率的影响。我国CRTSⅢ型轨道板长为5 600 mm，板宽为2 500 mm，板厚为210 mm，板缝宽为70 mm，钢轨为60 kg/m钢轨，其余轨道结构参数取值见表1。其中，CA砂浆刚度为换算后的平均支承刚度。从频散曲线中随机取3条曲线波数位于π/l处的频率作为研究对象，其结果见图3。从图3可以看出，当M=20，N=η=8时，带隙频率基本收敛。因此，将其作为后面的算例分析中截断项数取值。

表1 无砟轨道结构参数

图3 带隙频率收敛性分析

2.2 将轨道板考虑为质量块的垂向振动频散特性分析

国内外现有涉及无砟轨道振动特性的研究中，或忽略轨道板的影响[7]，或将轨道板考虑为梁[24]，鲜有研究将其作为板进行建模。因此，在运用本文方法之前，有必要在原方法的基础上，不考虑轨道板形变将其简化为具有无限刚度的有限长质量块，对模型垂向振动频散特性进行分析。将结果作为对照组，以探索本文将轨道板考虑为Mindlin板的分析方法的优越性，其带隙结果见图4。为便于分析计算，引入参数m，用于代替波矢k，且满足

(24)

式中：k表示被限制在第一不可约Brillioun区内的波矢。

轨道结构可看成沿轨道纵向的一维周期结构，对于一维周期结构，其第一不可约Brillioun范围即k的取值范围为[-π/l,π/l]，由式(24)可知，参数m的取值为[-1,1]。

从图4中可以看出，将轨道板考虑为质量块时，在0～1 200 Hz范围内，结构共存在三阶垂向振动带隙，频率范围分别为0～89.7、92.7～217、1 002～1 046 Hz。

图4 轨道板考虑为质量块时结构带隙结果

2.3 将轨道板考虑为板的垂向振动频散特性分析

将轨道板考虑为Mindlin板时，运用理论法计算得到0～1 200 Hz范围内无砟轨道结构垂向振动频散曲线见图5。

由图5可见，无砟轨道结构垂向振动会产生四阶带隙，前三阶在低频产生(0～89.2、104.8.5～115.5、119.5～195.4 Hz)，第四阶在高频产生(1 001.7～1 032.0 Hz)，其中，第一阶、第三阶带隙带宽较大，分别为89.2、75.9 Hz。

图5 垂向振动频散曲线

从结果可以看出，将轨道板考虑为Mindlin板和质量块所得到的带隙结果基本吻合，区别在于考虑为质量块时得到的二阶带隙频率为92.7～217 Hz，而考虑为Mindlin板时在这一频率范围内生成了两条带隙，分别为104.8～115.5、119.5～195.4 Hz。这是因为将轨道板考虑为Mindlin板后，轨道板不仅存在刚体运动，还存在因形变而产生的弹性体运动，多种运动耦合后在该频段范围内产生通带，从而将一条带隙分割为两条，因此该频段带隙的产生受轨道板影响较大。此外，这一结果也表明，更为细致地考虑轨道板在结构振动中的表现，可更全面地展现轨道结构的带隙特性。

为进一步验证前述理论分析的正确性，建立单个周期无砟轨道结构有限元模型，见图6。利用对称性建模，将钢轨与轨道板均考虑为实体单元，扣件与道砟简化为连接弹簧，其中，钢轨与轨道板通过扣件弹簧连接，轨道板与基础通过道砟弹簧连接。将有限元模型划分为有限个通过节点连接的单元，单元类型为Langrage-quadratic，为满足计算精度需求，将其划分为521 655个域单元、138 250个边界元以及24 065个边单元，求解总自由度为2447 210个。选用多物理场仿真软件Comsol Multiphysics中的固体力学模块进行求解，在图6中板的内侧边界添加对称边界，其余为自由边界，同时钢轨右端截面定义为源边界，左端截面定义为目标边界，在钢轨源边界添加floquet周期性边界条件，计算机配置为2.7 GHz CPU，192 GB内存，计算时长约14 h。理论法采用Matlab进行编程计算，耗时约32 s，表明本文方法计算效率较高。图7～图9给出了解析法与有限元法得到的垂向振动频散曲线局部对比分析图，不同方法计算得到的带隙频率对比见表2。

图6 双层无砟轨道结构有限元模型

图7 轨道板考虑为板时第一、二阶带隙对比分析

图8 轨道板考虑为板时第三阶带隙对比分析

图9 轨道板考虑为板时第四阶带隙对比分析

表2 不同方法计算得到的带隙频率对比

由表2可见，比较四阶带隙频率处频散曲线可以发现，理论解和有限元法计算得到的带隙频率位置基本一致，整体频散曲线也基本吻合，说明本文方法不仅计算效率高，能保证计算结果的可靠性，还表明0～1 200 Hz范围内，采用Timoshenko梁与Mindlin板耦合模型能够较为准确地表征双层周期性无砟轨道结构的带隙特性。

3 带隙形成机理分析

本文2.3节从无砟轨道结构的带隙计算结果中选出了位于0～1 200 Hz范围内的四条带隙，其中前三阶带隙为局域共振带隙，第四阶带隙为Bragg带隙。

钢轨和周期性排列的局域共振单元是局域共振带隙形成的主要条件。局域共振带隙的起始与截止频率对应振动模态可以用“质量-弹簧”模型描述[25-27]。本文建立了与双层无砟轨道结构相对应的“质量-弹簧”模型，以进一步分析带隙中局域共振带隙的产生机理，见图10。其中，k1为CA砂浆支承总刚度；k2为扣件垂向总刚度；m1和m2分别为轨道板和钢轨质量，且满足关系式：k1=2kslc，m1=2ρPlch，k2=18kz，m2=2ρBAB·(l+s0)。

图10 双层无砟轨道结构局域共振带隙分析模型

由图10(a)可得局域共振带隙的起始频率fs为

fs1=0

(25)

(26)

由图10(b)可得其特征频率ω满足关系式为

(27)

从而可以得到局域共振带隙截止频率fc为

(28)

(29)

(30)

(31)

将各项轨道参数代入式(26)和式(28)中，得到结构局域共振带隙起始频率和截止频率，从而可估算结构局域共振带隙频率范围为：0～89.3 Hz，111.4～193.7 Hz。估算结果与表2基本一致。需要强调的是，考虑为Mindlin板后，轨道板不仅存在刚体运动，还存在因形变而产生的弹性体运动，多种运动耦合后在该频段范围内产生通带，从而将其中一条低频带隙分割为两条，证明了2.3中两种方法得到的前三阶带隙为局域共振带隙，再次证明本文方法求解轨道结构带隙特性的准确性。

根据Bragg散射机理[28]，当弯曲波以波长λ在周期性轨道结构中传播时，此时晶格常数为l0(即扣件间距)，则Bragg条件可表示为：2l0=nλ(n=1,2,…)。当n取最小值1时，对应的弯曲波频率为第一阶“pinned-pinned”频率的起始频率，此时位于扣件处的钢轨垂向位移为零，而跨中处钢轨垂向位移达到最大值[29]，振型见图11，且可设钢轨的位移场满足关系：

图11 第一阶“pinned-pinned”频率振型

v=H1sin(nπx/l0)

(32)

β=H2cos(nπx/l0)

(33)

式中：H1和H2为未知系数，且不为零。

钢轨的振动微分方程为[30]

将式(32)和式(33)代入方程组(34)，可得

(35)

要使式(35)成立，则有

(36)

将式(36)行列式展开，代入n=1及各项轨道结构参数可解得：f1=1 001.7 Hz，f2=6 598.3 Hz。其中，f1为以弯曲位移为主对应的“pinned-pinned”频率，f2为以剪切位移为主对应的“pinned-pinned”频率。因此，第一阶“pinned-pinned”频率的起始频率为1 001.7 Hz，与表2中各方法所求最后一阶带隙频率起始频率相吻合，证明第四阶带隙为Bragg带隙。从式(36)可以看出，Bragg带隙的起始频率仅由扣件间距和钢轨物理性质所决定，这也可以解释将轨道板由考虑为质量块转变为考虑成Mindlin板后，第四阶带隙起始频率基本不变的现象。此外，Bragg带隙的带宽主要由扣件刚度决定[31]。

4 结论

本文以我国CRTSⅢ型无砟轨道结构为研究对象，提出了一种基于能量泛函变分原理和平面波级数展开的混合方法，用于计算周期性无砟轨道结构振动带隙。理论计算结果与有限元仿真结果和简化模型结果进行了对比分析，得到如下结论：

(1)本文提出的基于能量泛函变分原理的新方法在保证计算准确性的前提下，极大地提高了计算效率(同模型下，计算速度约为有限元法的上千倍)。本文方法克服了传统解析方法需要直接求解微分方程组以及难以对完整轨道板进行建模的问题，可更为全面地分析轨道板模态对结构振动特性的影响。

(2) 0～1200 Hz频率范围内，用Timoshenko梁和Mindlin板能较准确地表征周期无砟轨道结构带隙特性，且垂向振动产生0～89.2、104.8.5～115.5、119.5～195.4、1 001.7～1 032.0 Hz四阶带隙。

(3) 与将轨道板考虑为质量块相比，考虑为Mindlin板后，因轨道板存在多种模态，使得其不仅存在刚体运动，还存在因形变而产生的弹性体运动，多种运动耦合后产生通带，从而将在该频段范围内一条带隙分割为多条，可更为细致地考虑轨道板在结构振动中的表现，更全面地展现轨道结构的带隙特性。

(4) 文中求得带隙中前三阶为局域共振带隙，第四阶为Bragg带隙。这两类带隙频率均可通过估算得到其大致范围，其中，局域共振带隙频率范围可由“质量-弹簧”模型估算得到，Bragg带隙频率范围可由不同的“pinned-pinned”频率振型分析计算得到。此外，Bragg带隙起始频率仅由扣件间距和钢轨物理性质决定。