1. 引言

当高速列车进入隧道时，车头前方空气被压缩，形成初始压缩波在隧道内以音速向前传播，并在隧道出口处向外辐射微气压波，并且对周边环境造成噪声污染甚至结构破坏。压缩波波形会随着传播距离的增加而演化，微气压波的振幅取决于压缩波在隧道内波形的演变以及达到隧道出口处的压缩波压力梯度。因此，压缩波传播过程中的演化是预测和控制微气压波问题的关键[1]。对于压缩波在隧道内长距离传播，通常采用一维特征线法对传播过程进行预测，并且假设壁面摩擦效应为准稳态模型，但压缩波的快速瞬态特性会导致壁面摩擦力呈现非定常性，诱导壁面边界层发展，在长距离传播后，非定常摩擦带来的误差会很明显[2]，忽略非定常摩擦效应会低估压缩波的衰减速率，进而影响微气压波强度的预测精度，以往的研究常常将隧道看作光滑壁面，忽略了壁面粗糙度对非定常摩擦模型的影响，并且粗糙管道中的流动瞬态响应需结合壁面粗糙度与流体惯性力的耦合效应，用加权函数考虑粗糙管道中的非定常摩擦效应[3]。

1968年Zielke [4]导出了瞬态层流管内壁面剪切应力与瞬时平均速度和加权历史速度变化的关系式。适用于计算粘性流体中水锤现象的特征方法，用特征线法建立了瞬态流动的一维模型，并得到了粘度随频率变化的响应曲线，能较好地预测波传播过程中和频率相关的摩擦力。

1995年Vardy [5]建立了光滑管内流动瞬态摩擦的加权函数模型，假定湍流粘度在管道核心周围的厚边界层内呈线性变化，从而使其适用于高雷诺数流动，对管道中的历史流动参数信息加权函数进行了近似，预测的结果接近实验结果。

2002年Mohamed [6]推导了Vardy-Brown非定常摩擦方程的一个精确、简单、有效的近似，并将该近似应用在非定常管道流动的一维特征解中，模拟瞬态湍流管道流动中的剪切应力，并用特征线法求解得到水锤方程。提高了计算效率，降低了计算时的内存存储。

2004年Vardy [3]利用理想径向粘度分布，推导了一维管道中完全粗糙壁面流动的非定常表面摩擦的加权函数模型，假设在足够短的瞬态过程，截面上的粘度分布是恒定的，外环内的涡流粘度假设为线性变化，得到了完全粗糙管道内流动非定常摩擦瞬态加速度的计算模型。

2017年Duan [7]探讨了非定常摩擦项的相关性，分别比较了基于加权函数的模型和基于瞬态加速度的模型，通过实验研究和一维，二维数值模拟，用实验数据标定了非定常摩擦模型中的参数，并验证了各个模型所适用的领域和局限性。对于具有高频波动特性的管道瞬态流动，必须考虑非定常摩擦模型和湍流模型。

2021年Iyer [8]提出对于隧道压缩波传播的准确预测取决于隧道壁面摩擦的模拟能力，隧道壁面的摩擦损失是准稳态分量和非稳态分量的总和，并且通过实验测量证明，在接近压缩波和壁面处的表面摩擦非定常效应对传播过程有较大影响，特别是表面切应力在隧道空间内的分布影响压缩波的传播速度和波形的衰减变形。在隧道压缩波传播的情况下，不考虑壁面粗糙度和非定常摩擦带来的影响，会使得理论预测的压缩波波形和实验测量的数据有严重的偏差，并且随着传播距离的增加，误差积累得严重。对于粗糙壁面，建立引入粗糙度的加权函数和非定常摩擦因子进行求解，只考虑瞬态流动前的流动参量，而不是所有的历史流动参量，减少了大量的计算量和计算时间。

本文采用一维可压缩非定常不等熵流动模型对压缩波传播过程进行计算，对比分析不同管道摩擦模型在高速列车隧道压缩波传播中应用的准确性，并和德国Eurwang隧道实测数据[9]进行对比分析，并考虑不同的隧道壁面粗糙度和不同的非定常摩擦因子下摩擦模型的准确性，为后续高速列车预测压缩波在隧道中传播提供非定常摩擦模型，提高计算精度。

2. 数学物理模型及数值计算方法

2.1. 一维流动基本方程

初始压缩波经过的地方会引起空气压力、流速等参数的变化。由于隧道长度远大于隧道水力直径，压缩波沿整个隧道长度的传播时间也远大于在横截面上的传播时间。且压缩波的传播过程还受到隧道壁面摩擦等耗散因素的影响。因此可将初始压缩波在隧道内传播所引起的空气流动简化为一维可压缩不等熵非定常流动[10]。将空气与隧道壁面的摩擦和传热引入到动量方程和能量方程中，利用质量、动量和能量守恒定律，推导出压缩波传播过程的控制方程，其具体为：

连续性方程：

$\frac{\partial \rho }{\partial t}+\rho \frac{\partial u}{\partial x}+u\frac{\partial \rho }{\partial x}+\rho \frac{u}{F}\frac{\text{d}F}{\text{d}x}=0$ (1)

动量方程：

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\frac{1}{\rho }\frac{\partial p}{\partial x}+G\left(u\right)=0$ (2)

能量方程：

$\frac{\partial p}{\partial t}+u\frac{\partial p}{\partial x}-a^2\left(\frac{\partial \rho }{\partial t}+u\frac{\partial \rho }{\partial x}\right)=\rho \left(\kappa -1\right)\left(q+uG\left(u\right)\right)$ (3)

式中，F为隧道净空面积，ρ为空气密度，u为气流速度，a为声速， $\kappa$ 为比热容比，p为压力。q为传热项， $G\left(u\right)$ 为摩擦项。

2.2. 壁面摩擦模型

2.2.1. 定常摩擦模型

定常摩擦模型采用Darcy-Weisbach方程导出的压力损失形式，在隧道压缩波传播模型中，考虑可压缩流动的情况下，重力对流动的影响可忽略不计，因此忽略定常摩擦模型中的1/g项。则本文的定常摩擦模型为：

$G_s=\frac{\lambda }{2D}u\left|u\right|$ (4)

其中，λ为达西摩擦系数，u为隧道内空气流速，D为隧道当量直径。

2.2.2. 非定常摩擦模型

当流体在管道或隧道中流动为瞬态流动时，即流体速度快速变化，会导致壁面摩擦力的显著波动，与定常摩擦不同，非定常摩擦模型与非定常流动的雷诺数有关。定常摩擦理论根据隧道截面的平均流速和平均加速度来计算壁面带来的表面切应力，与实际隧道中压缩波传播引起的瞬变流动不同，近壁区的边界层内的切应力不能用平均切应力来表示，因此会对传播过程中的压力梯度计算带来很大误差，因此为了准确预测压缩波传播过程的流动情况，需要同时考虑定常和非定常摩擦项。本文选取两种不同的非定常摩擦模型进行对比，并分析哪种摩擦模型适用于预测初始压缩波在隧道中传播过程的衰减效应。

2002年Mohamed等人[6]在Vardy-Brown光滑管道非定常摩擦模型的基础上，推导了一个近似的，更为简单有效的加权函数，只需要存储单个时间步长下的流动变量，极大地减小了运算时的内存占用和计算时间，采用的壁面切应力为：

$\tau =\frac{4\mu }{D}{\displaystyle {\int }_0^{t}W\left(t-t^{\prime }\right)\frac{\partial u}{\partial t^{\prime }}\text{d}{t}^{\prime }}$ (5)

式中，μ为空气的动力黏度系数，D为隧道当量直径，u为当前时刻气体流速， $W\left(t-t^{\prime }\right)$ 为历史时刻流动参量的加权函数，定义为：

$w\left(t\right)=\frac{\alpha {\text{e}}^{-\beta t}}{\sqrt{\pi t}}$ (6)

式中各参数定义为：

$\alpha =\frac{D}{4\sqrt{\upsilon }};\beta =\frac{0.54\upsilon {\mathrm{Re}}^{\kappa }}{D^2};\kappa ={\log }_{10}\left(\frac{14.3}{{\mathrm{Re}}^{0.05}}\right)$ (7)

式中 $\upsilon$ 为空气的运动黏度系数，D为隧道当量直径，Re为隧道内气流的雷诺数。

动量方程中的非定常摩擦项与壁面切应力的关系为：

$G_{us}=\frac{4\tau }{\rho D}$ (8)

2021年Iyer等人采用的非定常摩擦模型用非定常摩擦因子进行了代替，表示为：

$G_{us}={\epsilon }_{us}{\displaystyle {\int }_0^{T}W\left(t\right)\frac{\partial u}{\partial t}\text{d}t}$ (9)

式中， ${\epsilon }_{us}$ 为非定常摩擦因子， $W\left(t\right)$ 为历史时刻流动参量的加权函数，定义为：

$w\left(t\right)=\frac{{\text{e}}^{-B\psi }}{2\sqrt{\pi \psi }}$ (10)

式中各参数定义为：

$B=\frac{{\mathrm{Re}}^{\kappa }}{12.86};\kappa ={\log }_{10}\left(\frac{15.29}{{\mathrm{Re}}^{0.0567}}\right);\psi =\frac{4\upsilon \theta }{D^2};\theta =t-t^{\prime }$ (11)

式中 $\upsilon$ 为空气的运动黏度系数，D为隧道当量直径，Re为隧道内气流的雷诺数。

2.3. 数值计算方法

本文采用一维可压缩非定常不等熵流动模型，并用特征线法对控制方程进行求解，使得其转化为方便计算的常微分方程组，特征线法分为正步进特征线法和逆步进特征线法，本文采用改进的逆步进特征线法对网格进行计算，采用修正的迭代算法确定特征线在网格中的插值点位置，插值点和网格的相对位置如图1所示。

Figure 1. Discrete grids and interpolation points

图1. 离散网格和插值点

针对控制方程中的定常和非定常摩擦项进行数值计算，可将非定常摩擦模型中的积分项转变为用网格点表示的离散求解公式，Mohamed文献采用的非定常摩擦模型离散公式第n + 1步的切应力 ${\tau }_u^{n+1}$ 可用第n步的切应力 ${\tau }_u^n$ 计算得到，离散公式为：

${\tau }_u^{n+1}=\frac{\rho \sqrt{\nu }}{\Delta t\sqrt{\beta }}\text{erf}\sqrt{\beta \Delta t}\left[Q^{n+1}-Q^n\right]+\xi \exp \left(-\beta \Delta t\right)\left[\sqrt{\frac{n+1}{n}}+\frac{1}{2n}\ln \frac{\sqrt{\frac{n+1}{n}}-1}{\sqrt{\frac{n+1}{n}}+1}\right]{\tau }_u^n+\left(1-\zeta \right){\tau }_u^n$ (12)

Iyer文献采用的非定常摩擦模型离散公式为各时间步保存的流动参量的代数和， $t+\Delta t$ 时刻受到的壁面非定常摩擦力表示为：

$f_{us}\left(t+\Delta t\right)=\frac{16\upsilon }{D^2}{\displaystyle \underset{i=1}{\overset{N}{\sum }}W\left(t\right){\text{e}}^{-\frac{4n_i\upsilon }{D^2}}}+\frac{m_i{D}^2}{4n_i\upsilon }\left(\frac{1-{\text{e}}^{-\frac{4n_i\upsilon }{D^2}}\partial u}{\partial t}\right)$ (13)

3. 压缩波传播摩擦模型探究

3.1. 定常摩擦模型参数探究

由定常摩擦模型公式得知，压缩波传播过程中受到的衰减效应与定常摩擦系数λ和空气流速u有关，而空气流速与初始压缩波的最大压力幅值有关。本节在只考虑定常摩擦项的情况下对比德国文献中的实测数据[9]和西班牙文献初始压缩波[11]在隧道中的传播现象，采用的初始压缩波特征参数如表1所示。

Table 1. Maximum pressure and gradients for different references

表1. 不同文献的最大压力幅值和压力梯度

分析定常摩擦模型的两个参数对传播衰减激化过程的影响，选取定常摩擦系数0.005和0.015进行对比在隧道中传播3000 m后的压缩波压力幅值。由图2可知，初始压缩波压力幅值越大，压缩波传播过程中受到的衰减现象越明显，且随着传播的距离越远，衰减越大，定常摩擦系数越大，传播过程的衰减效应越明显，因此，对于不同的隧道，可以采用不同的定常摩擦系数进行计算，根据实验数据确定合适的定常摩擦系数。

Figure 2. Comparison of parameters of steady friction model

图2.定常摩擦模型参数对比

3.2. 非定常摩擦模型参数探究

由非定常摩擦模型公式可知，非定常摩擦力由上一时刻流体的流动参量以及非定常摩擦因子决定，因此对比不同初始压缩波在不同非定常摩擦因子的情况下计算得到的压力幅值衰减效应，分别选取非定常摩擦因子 ${\epsilon }_{us}$ 为20和15进行对比，并分析同一初始压缩波在不同非定常摩擦因子下传播5000 m后的压力幅值变化规律，探究其对压缩波传播过程中非定常摩擦项的影响，并且将压缩波传播到隧道不同位置处的压力幅值与初始压缩波压力幅值的比值定义为压缩波衰减率 $\xi$ ，表示为：

$\xi =\frac{\Delta p_x}{\Delta p_0}$ (16)

由图3可知，考虑非定常摩擦因子使得压缩波传播过程的衰减现象延后了，并且衰减效应大于只考虑定常摩擦模型的情况，非定常摩擦因子越大，压缩波的衰减现象越滞后，并且衰减后的压力幅值更低，而初始压缩波压力幅值越大的压缩波传播过程中受到非定常摩擦效应带来的衰减现象越明显。由图4可知，不同非定常摩擦因子带来的衰减效应不同，非定常摩擦因子越大，衰减现象越明显，并且，传播距离越远，衰减幅度越大， ${\epsilon }_{us}=15$ 时传播到5000 m处压缩波的衰减比为0.84， ${\epsilon }_{us}=20$ 时传播到5000 m处压缩波的衰减比为0.77。

3.3. 压缩波传播中的应用

根据Mohamed文献和Iyer文献的非定常摩擦模型，对比不同非定常摩擦模型应用于高速列车隧道压缩波传播计算程序的正确性和精确性，并且根据德国Eurwang隧道实测数据[9]，确定对于不同的非定常摩擦模型下的定常摩擦系数和非定常摩擦因子，采用德国文献初始压缩波，在传播距离为7300 m后的压缩波波形以及梯度进行对比，由图5可知，定常摩擦系数为0.005时与实测数据压力幅值的误差为4.70%，定常摩擦系数为0.015时与实测数据压力幅值的误差为6.89%，因此定常摩擦系数取0.005更适应德国Eurwang隧道壁面的粗糙程度，能很好地反映出初始压缩波在隧道内的摩擦效应。当考虑Iyer非定常摩擦模型时，由图6可知，非定常摩擦因子 ${\epsilon }_{us}=0.5$ ，定常摩擦系数 $f=0.015$ 时，与实测数据压力幅值的误差为5.40%，取 ${\epsilon }_{us}=0.1$ ， $f=0.02$ 时，误差为3.70%，而取 ${\epsilon }_{us}=0.1$ ， $f=0.015$ 时，误差为0.71%，能够很好模拟Eurwang隧道壁面的粗糙程度。

Figure 3. Comparison of unsteady friction factors

图3. 非定常摩擦因子对比

Figure 4. Comparison of unsteady friction factor attenuation ratio

图4. 非定常摩擦因子衰减比对比

Figure 5. Mohamed parameter comparison of unsteady friction model

图5. Mohamed非定常摩擦模型参数对比

4. 结论

本文通过一维可压缩非定常不等熵流动模型探究了定常摩擦模型各参数对于压缩波传播过程衰减效应的影响，以及不同文献中的非定常摩擦模型各参数对于压缩波传播过程衰减效应的影响，并探究了

Figure 6. Iyer parameter comparison of unsteady friction model

图6. Iyer非定常摩擦模型参数对比

不同非定常摩擦模型在隧道压缩波传播计算程序中的应用，确定了不同摩擦模型下参数。当压缩波传播模型只考虑定常摩擦效应时，定常摩擦系数越大，传播过程中的衰减效应越明显，初始压缩波最大压力幅值越大时，传播过程受到的衰减效应也越明显；当压缩波传播模型考虑定常摩擦模型和非定常摩擦因子时，衰减效应更加明显，也能更精确地模拟传播过程中的衰减现象，非定常摩擦因子越大，传播过程的衰减现象越明显，并且对于初始压缩波最大压力幅值越大的波形衰减越明显， ${\epsilon }_{us}=15$ 时传播到5000 m处压缩波的衰减比为0.84， ${\epsilon }_{us}=20$ 时传播到5000 m处压缩波的衰减比为0.77；考虑Mohamed摩擦模型时，定常摩擦系数取0.005可以很好地模拟Eurwang隧道壁面的粗糙程度，而考虑Iyer非定常摩擦模型时，非定常摩擦因子取0.1，定常摩擦系数取0.015可以很好地模拟Eurwang隧道壁面的粗糙程度，对于不同隧道，应采取不同的系数建立隧道内的摩擦模型，能够精确地预测长隧道传播下的压缩波衰减效应。

参考文献