一种高速铁路隧道道砟铺设量的计算方法

1.本发明涉及高速铁路隧道微气压波的缓解方法，具体为一种高速铁路隧道道砟铺设量的计算方法。


背景技术：

2.当列车高速入隧产生的压缩波传至隧道出口时，洞口处将辐射出以次声能量为主的微气压波(micro pressure wave，简称mpw)，会对周围环境产生极大危害。由于mpw能量主要取决于传至隧道出口的压缩波最大压力梯度，因此实现对压缩波隧道内传播演化过程的精确模拟非常必要
1.。已有文献表明，在隧道内部铺设碎石道砟是较为有效的微气压波缓解方式，但碎石道砟对微气压波的缓解效果与碎石道砟的铺设量有关。因此对于实际需要铺设碎石道砟的隧道，精确预测将微气压波降低至合理水平所需要的道砟数量是很有必要的。
3.目前，最常见的处理方式是将碎石道砟抑制波前陡化的作用统一视为非定常摩擦作用，通过调整摩擦系数的大小，间接模拟不同数量的碎石道砟对隧道微气压波的缓解作用
2.。但该处理方式不能反映碎石道砟对压缩波的实际作用效果，因此计算精确度较低。此外，由于摩擦系数与道砟数量之间并无直接的映射关系，在已知隧道非定常摩擦系数后，采用该方法也很难合理预测的碎石道砟的铺设量。
4.[1]vardy ae.generation and alleviation of sonic booms from rail tunnels,proceedings of the institution of civil engineers engineering and computational mechanics 161september 2008issue em3 pages 107
–
119.
[0005]
[2]fukuda t,nakamura s,miyachi t,et al.influence of ballast quantity oncompression wavefront steepening in railway tunnels,proceedings of theinstitution of mechanical engineers part f journal of rail and rapid transit,2019,234(6),607-615.


技术实现要素：

[0006]
现有的研究在分析碎石道砟对隧道压缩波陡化的抑制作用时，往往将碎石道砟对压缩波的作用效果归结为摩擦作用，通过修正隧道壁面的非定常摩擦系数来间接考虑道砟的影响。但实际碎石道砟大多为孔隙率等于0.6的多孔结构，其对压缩波波前的气体有着很强的吸收作用，上述方法并不能考虑该作用对压缩波波前的影响，因此并不能准确的模拟碎石道砟对压缩波对波前陡化的抑制作用。此外，由于摩擦系数与道砟数量之间并无直接的映射关系，在已知隧道非定常摩擦系数后，采用该方法也很难合理预测碎石道砟的铺设量。
[0007]
针对现有方法很难合理预测碎石道砟铺设量的问题，本发明提出了一种高速铁路隧道道砟数量的预测方法，该方法可将质量项主导的碎石道砟和动量项主导的壁面摩擦进行有效辨识和区分，进而使道砟铺设量与储气容积进行直接关联。
[0008]
本发明是采用如下的技术方案实现的：一种高速铁路隧道道砟铺设量的计算方法，包括以下步骤：
[0009]
(1)给定需加设道砟隧道的长度、横截面积以及气室阵列体积比、微气压波的辐射角、隧道外部环境的测点到隧道出口的距离；
[0010]
(2)通过实车试验测得隧道入口的初始压缩波以及出口位置压缩波的波前形状，并通过该试验数据修正一维特征线法的定常摩擦系数f以及非定常摩擦系数ε
us
；
[0011]
(3)通过试验测量需要铺设的碎石道砟的孔隙率以及阻力系数k1和k2；
[0012]
(4)将上述得到的数据作为已知条件，采用考虑道砟影响的特征线法数值计算程序，可算得隧道出口位置的微气压波大小；结合工程中许用的微气压波限值，判断计算所得隧道出口位置的微气压波是否小于许用微气压波大小，若满足，直接输出气室阵列体积比；若不满足，调整气室阵列体积比，再次进行计算，直到计算得到满足隧道出口的微气压波值；
[0013]
特征线法数值计算程序为：气室阵列中每个气室作为一个计算节点，每个计算节点设置双网格点，即设定当前气室与隧道的连接器的左右两端为a和b，两个特征线方程为沿la方向(p
ab-p
l
)+(ρc)
l
(u
a-u
l
)＝-(ρc)
l
δt(ua|ua|f+ε
us
yq)和沿rb方向(p
ab-pr)-(ρc)r(u
b-ur)＝(ρc)rδt(ua|ua|f+ε
us
yq)，其中，(ρc)
l
为当前气室左侧节点在当前时刻和上一时刻的密度和声速乘积的插值，(ρc)r为当前节点上一时刻密度与声速的乘积和当前气室右侧节点上一时刻密度与声速的乘积的插值，p
l
表示当前节点左侧的隧道压力，pr表示当前节点右侧的隧道压力，u
l
表示当前节点左侧的气体速度，ur表示当前节点右侧的气体速度，δt为时间步长；连续性方程为气室方程为式中v为气室体积，特征线方程可通过连续性方程与气室方程进行耦合，同时，连续性方程可将当前节点左右两点的速度差u
a-ub与通过连接器流向气室的流速相关联，即：式中，a为隧道的横截面积，表示通过计算网格点的总质量流；
[0014]
在每个时间步长中，依次在每个计算节点上求解上述方程，可得隧道压力p
ab
、气室压力p
ch
、当前节点的左端点气体速度ua、右端点气体速度ub、隧道与气室间的质量流变化速率五个未知量的值；该求解过程是逐步迭代的，最后得到隧道出口气室位置处的隧道压力p
ab
，此时，将隧道压力p
ab
对时间求导，即可得到隧道出口压缩波的最大压力梯度，通过公式即可求得隧道出口位置的微气压波大小，式中，a为隧道的横截面积；ω为隧道出口向外辐微气压波的辐射角；z为测点到隧道出口的距离；c为声速；为隧道出口压缩波的最大压力梯度；
[0015]
(5)结合隧道横截面积、碎石道砟孔隙率即可求得每米隧道所需铺设的碎石道砟体积，
[0016]
本发明的数值计算方法可以较好的预测碎石道砟对波前陡化的抑制作用。此外，对于工程上需要加装碎石道砟去改善出口微气压波的隧道，使用本发明的数值计算方法可以很好的预测出隧道内合理的道砟铺设量，因此在工程设计和分析中具有更好的实用性。
附图说明
[0017]
图1为碎石道砟与气室阵列作用原理示意图。
[0018]
图2为特征线法的数值网格结构图。
[0019]
图3为整个数值计算程序求解的网格结构图。
[0020]
图4为计算所需道砟数量的流程图。
[0021]
图5为与道砟实测结果
2.对比图。
具体实施方式
[0022]
为缓解高速铁路隧道出口微气压波，本发明提出了一种高速铁路隧道道砟铺设量的计算方法。与现有技术方法相比，本发明可将质量项主导的碎石道砟和动量项主导的壁面摩擦进行有效辨识和区分，进而使道砟铺设量与储气容积进行直接关联，因而在工程设计和分析中具有更好的实用性。由于采用了1d模型，因此在保证计算速度的同时，还可以较好的预测碎石道砟对波前陡化的抑制作用，并且可以较好的反应道砟铺设量的影响。其具体方案如下：
[0023]
压缩波在隧道内的传播的数值计算方法
[0024]
压缩波沿隧道向前传播时，所包含的控制方程如下：
[0025]
连续性方程：
[0026][0027]
动量方程：
[0028][0029]
式中，ρ气体的密度；u表示沿隧道轴向的气体流速；t为时间；x表示沿隧道轴向的传播距离；p为气体的压力；g为摩擦力项，可表征为隧道定常摩擦gs和非定常摩擦g
us
两项之和。
[0030]
对上述两方程采用特征线法简化，可将其转化为更易处理的常微分方程：
[0031][0032]
其特征线方程为：
[0033][0034]
其中，f为定常摩擦项的常数，q为非定常摩擦项的常数，y表示非定常摩擦项的积分部分，c为声速。
[0035]
基于以上公式，通过设置初始边界条件，使用插值、迭代等方法即可对压缩波在隧道内的传播过程进行求解。
[0036]
道砟储气的模拟
[0037]
由于压缩波在隧道内传播时，道砟的主要作用是存储多余空气，且气流主要沿道砟深度方向流动，本发明将碎石道砟简化为沿隧道纵向分布的小型气室阵列，气室阵列中的气室彼此间互不连通，气室通过连接孔和隧道连通。
[0038]
压缩波沿隧道传播时，气室的主要作用是通过连接孔吸收波前位置的空气。隧道与气室间的压差δp取决于隧道与气室间的质量流及其变化率因此连接处动量方程为：
[0039][0040]
其中，k1和k2为连接处的阻力系数，a
acc
为等效惯性系数。当隧道中的气体流向气室时，气室内部的空气密度ρ
ch
可由下式获得：
[0041][0042]
其中，v为每单位长度隧道对应的气室体积，且假设气室内压力的变化在所考虑的较小时间尺度上几乎是等熵的，则气室内压力的变化率为：
[0043][0044]
式中，c
ch
为声音在气室中的等熵速度，且气室内的空气被视为理想气体，声速可以表示为：
[0045][0046]
式中，γ为比热容比。
[0047]
在比本发明数值方法大的多的时间尺度上，气室中的气体与壁面将会产生一定的热传递，但这一过程是逐渐发生的，它对早期气室内部压力快速升高的影响可以忽略不计。通过设置双网格点，可将上述气室阵列的储气作用融入特征线法的数值计算程序。其中，气室阵列间的距离相同，且均设置在双网格点上。此外，由于考虑了气室阵列的储气作用，因此双网格点前后位置的速度会发生突变，但该处压力和密度被视为是连续不变的。
[0048]
数值积分与迭代求解
[0049]
为研究气室阵列对压缩波波前演化的影响，数值计算中所采用的网格结构如图2所示。由于隧道与气室阵列间存在气流流动，为考虑气室阵列对压缩波演化的影响，在每个计算节点上都有设置一个双网格点。在数值求解时，计算节点处的速度会发生突变，但该处压力和密度被视为是连续不变的。
[0050]
求解过程以时间推进的方式进行，在求解t1+δt时刻前，需先得到t1时刻所有计算节点处的结果。在图2中，假设初始压缩波沿隧道的传播方向是从左到右的，则求解过程按照图2中(1)、(2)、(3)、(4)计算节点的顺序依次展开。此外，即使使用固定大小的网格结构，特征线的梯度也可能不是常数，因此在使用特征线法求解t1+δt时刻的物理参数时，需要对上一时刻的结果进行插值。其中，l位置处的数值结果可沿时间线通过l进行插值，同理，r处的数值结果可在时间t1通过空间插值获得。
[0051]
道砟铺设量的预测方法
[0052]
通过上述数学模型得到压缩波沿隧道传播的数值计算方法后，即可求得隧道内部所需铺设的道砟数量，其流程图如图4所示，具体执行方法如下：
[0053]
(1)给定需加设道砟隧道的长度、横截面积以及铺设方案(道砟在隧道内的铺设位置)、气室阵列体积比、微气压波辐射角、测点到隧道出口距离；
[0054]
(2)通过实车试验测得隧道入口的初始压缩波以及出口位置压缩波的波前形状，并通过该实验数据修正一维特征线法的定常摩擦系数f以及非定常摩擦系数ε
us
。
[0055]
(3)通过实验测量需要铺设的碎石道砟的孔隙率以及阻力系数k1和k2。
[0056]
(4)将上述得到的数据作为已知条件，采用考虑道砟影响的特征线法数值计算程序，可算得隧道出口位置的微气压波大小；结合工程中许用的微气压波限值，判断计算所得的微气压波是否略小于许用微气压波大小，若满足，直接输出气室阵列的体积比；若不满足，调整气室阵列的体积比，再次使用上述方法进行计算，直到满足隧道出口的微气压波限值。
[0057]
特征线法数值计算程序为：
[0058]
沿la方向
[0059]
(p
ab-p
l
)+(ρc)
l
(u
a-u
l
)＝-(ρc)
l
δt(ua|ua|f+ε
us
yq)
ꢀꢀꢀ
(9)
[0060]
沿rb方向
[0061]
(p
ab-pr)-(ρc)r(u
b-ur)＝(ρc)rδt(ua|ua|f+ε
us
yq)
ꢀꢀꢀꢀꢀꢀ
(10)
[0062]
其中，(ρc)
l
为当前气室左侧节点当前时刻和上一时刻的密度和声速乘积的插值，(ρc)r为当前节点上一时刻密度与声速的乘积和当前气室右侧节点上一时刻密度与声速的乘积的插值，p
l
表示当前节点左侧的隧道压力，pr表示当前节点右侧的隧道压力，u
l
表示当前节点左侧的气体速度，ur表示当前节点右侧的气体速度；
[0063]
方程(5)和(7)描述了气室的对波前演化的影响，其在数值上近似写为：
[0064][0065]
以及
[0066][0067]
其中公式(11)和(12)中为当前节点上一个时刻的质量流，k1和k2为连接处的阻力系数。
[0068]
特征线方程(9)和(10)可通过连续性方程与气室方程(11)、(12)进行耦合。同时，连续性方程(11)可将当前节点左右两点的速度差u
a-ub与通过连接器流向气室的流速相关联，即：
[0069][0070]
上式中，a为隧道的横截面积，表示通过计算网格点的总质量流。
[0071]
在每个时间步长中，依次在每个网格点上求解方程(9)—(13)，可得隧道压力p
ab
、气室压力p
ch
，当前节点的左端点速度ua、右端点速度ub、隧道与气室间的质量流变化速率
五个主要未知量的值。
[0072]
气室阵列体积比与每单位长度隧道对应的气室体积v之间存在如下关系：
[0073]
在整个数值计算程序过程如下所示：
[0074]
(1)已知隧道长度l
tun
，且设定计算的网格长度为δx。
[0075]
则沿隧道长度方向划分的网格节点个数沿时间方向划分的时间节点个数nn＝n
j-1。
[0076]
(2)根据划分好的计算节点，可以进一步划分数值计算的网格，如图3所示。图中，横轴表示沿隧道长度方向的网格划分(从1到nj)，纵轴表示沿时间方向的网格划分(从0到n
n-1)。
[0077]
假设在中，x代表计算过程中的各关键参数，如p
ab
、ua、ub，n表示时间节点，j表示空间节点。此时，则表示在n
n-1时刻，出口位置气室(nj处)的各关键参数。
[0078]
(3)设置计算过程中的初始条件和边界条件。首先，设置初始条件，即已知其中，j从1到nj。接着设置隧道入口边界条件，即已知其中，n从1到n
n-1。接着根据(9)-(13)5个方程，依次按如下顺序对未知节点依次求解。
[0079]
即先按从左到右的顺序求解n＝1时刻的所有未知节点的关键参数x，即求解
[0080]
接着，按从左到右的顺序依次求解n＝2时刻的所有未知节点的关键参数x，即
[0081]
然后，依次求解n＝3时刻的所有未知节点的关键参数x，即然后，依次求解n＝3时刻的所有未知节点的关键参数x，即
[0082]
以时间节点n推进的顺序，按上述的求解方式一直求解，直到求解到n
n-1时刻的各位置节点的关键参数
[0083]
得到以上所有参数后，在隧道出口位置气室处，各个关键参数随时间节点的变化关系(其中，n＝0,1,2,
…
,n
n-1)即为所求。此时，将关键参数中的隧道压力p
ab
对时间求导，即可得到隧道出口压缩波的最大压力梯度。
[0084]
(4)最后，通过公式(15)即可求得隧道洞口的微气压波大小。
[0085][0086]
式中，a为隧道的横截面积；
[0087]
ω为隧道出口向外辐射微气压波的辐射角；
[0088]
z为测点到隧道出口的距离；
[0089]
为隧道出口压缩波的最大压力梯度。
[0090]
(5)结合隧道横截面积、碎石道砟孔隙率即可求得每米隧道所需铺设的碎石道砟体积，具体公式为：
[0091]

技术特征：
1.一种高速铁路隧道道砟铺设量的计算方法，其特征在于：包括以下步骤：(1)给定需加设道砟隧道的长度、横截面积以及气室阵列体积比、微气压波的辐射角、测点到隧道出口的距离；(2)通过实车试验测得隧道入口的初始压缩波以及出口位置压缩波的波前形状，并通过该试验数据修正一维特征线法的定常摩擦系数f以及非定常摩擦系数ε
us
；(3)通过试验测量需要铺设的碎石道砟的孔隙率以及阻力系数k1和k2；(4)将上述得到的数据作为已知条件，采用考虑道砟影响的特征线法数值计算程序，可算得隧道出口位置的微气压波大小；结合工程中许用的微气压波限值，判断计算所得的微气压波是否略小于许用微气压波大小，若满足，直接输出气室阵列体积比；若不满足，调整气室阵列体积比，再次使用上述方法进行计算，直到满足隧道出口的微气压波限值；特征线法数值计算程序为：气室阵列中每个气室作为一个计算节点，每个计算节点设置双网格点，即设定当前气室与隧道的连接器的左右两端为a和b，两个特征线方程为沿la方向(p
ab-p
l
)+(ρc)
l
(u
a-u
l
)＝-(ρc)
l
δt(u
a
|u
a
|f+ε
us
yq)和沿rb方向(p
ab-p
r
)-(ρc)
r
(u
b-u
r
)＝(ρc)
r
δt(u
a
|u
a
|f+ε
us
yq)，其中，(ρc)
l
为当前气室左侧节点当前时刻和上一时刻的密度和声速乘积的插值，(ρc)
r
为当前节点上一时刻密度与声速的乘积和当前气室右侧节点上一时刻密度与声速的乘积的插值，p
l
表示当前节点左侧的隧道压力，p
r
表示当前节点右侧的隧道压力，u
l
表示当前节点左侧的气体速度，u
r
表示当前节点右侧的气体速度；连续性方程为气室方程为特征线方程可通过连续性方程与气室方程进行耦合，同时，连续性方程可将当前节点左右两点的速度差u
a-u
b
与通过连接器流向气室的流速相关联，即：式中，a为隧道的横截面积，表示通过计算网格点的总质量流；在每个时间步长中，依次在每个计算节点上求解上述方程，可得隧道压力p
ab
、气室压力p
ch
、当前节点的左端点气体速度u
a
、右端点气体速度u
b
、隧道与气室间的质量流变化速率五个未知量的值；该求解过程是逐步迭代的，最后得到隧道出口气室位置处的隧道压力p
ab
，此时，将隧道压力p
ab
对时间求导，即可得到隧道出口压缩波的最大压力梯度，通过公式即可求得隧道洞口的微气压波大小，式中，a为隧道的横截面积；ω为隧道出口向外辐微气压波的辐射角；z为测点到隧道出口的距离；c为声速；为隧道出口压缩波的最大压力梯度；(5)结合隧道横截面积、碎石道砟孔隙率即可求得每米隧道所需铺设的碎石道砟体积，

技术总结
本发明涉及高速铁路隧道微气压波的缓解方法，具体为一种高速铁路隧道道砟铺设量的计算方法。本发明可将质量项主导的碎石道砟和动量项主导的壁面摩擦进行有效辨识和区分，进而使道砟铺设量与储气容积进行直接关联，因而在工程设计和分析中具有更好的实用性。由于采用了1D模型，因此在保证计算速度的同时，还可以较好的预测碎石道砟对波前陡化的抑制作用，并且可以较好的反应道砟铺设量的影响。此外，对于工程上需要加装碎石道砟去改善出口微气压波的隧道，使用本发明的数值计算方法可以很好的预测出隧道内合理的道砟铺设量，因此在工程设计和分析中具有更好的实用性。设计和分析中具有更好的实用性。设计和分析中具有更好的实用性。


技术研发人员：刘峰 雷海洋 邸娟 卫梦杰 马健斌 陈大伟
受保护的技术使用者：太原理工大学
技术研发日：2023.06.21
技术公布日：2023/9/13