基于数值试验的风沙区散体道床内尺度比研究

张智海，肖宏，令行

（1.北京交通大学土木建筑工程学院，北京，100044；2.北京交通大学轨道工程北京市重点实验室，北京，100044）

我国先后在中西部地区建成了多条沙漠铁路。由于自然环境的制约，沙漠地区铁路一直面临着风沙危害。沙粒侵入及沉积轻则填充道砟孔隙、改变道床弹性，重则导致轨道不平顺[1]、车轮磨耗加剧，轨道结构服役性能下降[2]，进而影响列车运行的平稳性和安全性[3-4]。因此，深入研究风沙区铁路轨道结构力学特性，对于指导线路养护维修作业和保障列车平稳运营意义重大。目前，关于风沙铁路轨道结构力学特性研究，严颖等[3]开展了含沙道床有效变形模量室内试验，利用颗粒流，建立了将沙粒粒径放大20倍的桶状数值试验模型，指出低含沙量下沙石混合料的变形模量基本不变；CARRASCAL等[5]结合室内沙埋扣件试验，发现沙粒对单个扣件使用性能影响不大；ESMAEILI 等[6]针对沙粒-道砟混合材料，开展了不同沙粒比室内道砟箱载荷试验，指出含沙量与道床累积垂向位移呈正相关；FACCOLI等[7]结合车轮滚动和滑动试验，发现沙粒会增加车轮磨损。

不同尺寸散粒体材料的强度和变形特性存在着明显的差异[8]，风沙道床是沙粒与道砟的自组织散粒体系，内部沙粒粒径与道砟粒径相差悬殊，在利用离散元法研究风沙区轨道结构力学特性时，需重点关注沙粒与道砟的多尺度效应问题。在颗粒多尺度效应研究方面，尹小涛等[9]设计了12种粒径数值试验，根据颗粒结构的破坏形态，分析了内尺度比对数值试验的影响，说明了岩土材料存在内尺度比；刘海涛等[10]开展了一系列不同尺度的粗粒土三轴试验，指出颗粒粒径和试样尺寸对粗粒土的力学性质有很大影响；YANG等[11]利用平行黏结模型，对颗粒进行了单轴压缩试验，发现颗粒粒径愈大，弹性模量受粒径作用效果越明显，而单轴压缩强度无明显变化；NARDIN等[12]基于已有接触模型，提出了一种相似关系，指出颗粒粒径越小，单轴压缩强度越大；POTYONDY 等[13]建立不同颗粒尺度的三维单轴压缩数值模型，分析了不同颗粒尺度下的弹性模量和内摩擦角的变化规律，指出三维数值试验参数与颗粒尺寸有一定联系。从以上研究可知：颗粒尺寸效应会影响其力学特性，现有研究大多针对单一颗粒体系，并未涉及多元混合颗粒的尺寸效应问题，且数值试验规模相对较小。关于沙粒对轨道结构力学响应方面的研究，大多是基于室内试验，未能从细观角度给出道砟与沙粒的相互作用关系，也没有建立起细观和宏观的全尺度分析理论与方法；仅有的离散元研究虽然给出了沙粒对道砟力学特性的影响，但在模拟过程中，直接将沙粒粒径放大20倍而未考虑沙粒粒径尺度效应的影响，不能揭示沙粒与道砟之间力链传递规律及内在机理，因此，研究风沙区轨道结构力学特性时，首先要解决沙粒的合理尺度选择问题，即道床中道砟平均粒径与沙粒平均粒径之比的问题。

综上，本文作者利用虚拟试验平台，借助“胞分法”建立不同沙粒粒径尺度的轨道结构离散元仿真模型，考虑到颗粒流模型的计算效率，在开展大规模计算前，先进行数值模型试算优化分析，确定沙粒粒径合理计算区间。通过分析不同沙粒粒径下道床宏细观力学特性，利用函数模型分析预测，得出沙粒粒径尺度与道床力学特性的关系曲线，探明沙粒粒径尺度与道床受力及位移的内在联系，并给出风沙区铁路道床中沙粒最佳模拟尺度，以期为后期风沙铁路道床服役性能研究及道床脏污颗粒尺度的选择提供理论依据。

1 数值试验模型

风沙道床中含有大量细沙，多尺度效应显著。考虑到计算耗时及模型平衡问题，本文拟采用PFC2D离散软件进行研究。根据文献[14]以及现场试验段线路实际状态，道床厚度取0.35 m，按照既有线一级碎石道砟颗粒级配生成粒径为25~63 mm的精细化不规则道砟颗粒模型。考虑到边界效应的影响，建立7根轨枕长度的风沙区多尺度轨道结构离散元模型，如图1所示。由图1可见：60 kg/m钢轨采用平行黏结模型，是由直径为0.15 mm、密度为510 kg/m3的28 个小球规则排列黏结而成[15]；新II型混凝土枕是采用颗粒重叠法建立块体，利用自编fish 语言导入PFC 中；扣件用ball 单元模拟，选用接触黏结模型；道床是利用课题组提出的“块体叠加拼装法”生成，选用线性接触模型模拟。

图1 离散元模型Fig.1 Discrete element model

“胞分法”是研究多元混合颗粒多尺度效应的一种有效方法，可以保证颗粒体系质量守恒，减少系统误差和颗粒流模型平衡时间。本文利用“胞分法”将真实沙粒粒径放大10 倍记为原胞体（一分胞体），按照等体积法进行分裂得到不同沙粒粒径的子颗粒记为胞元，不断更新道砟空隙中的沙粒胞元，得到不同沙粒粒径的离散元模型，如图2 所示。所建模型沙粒平均粒径最小为0.660 mm，最大沙粒平均粒径为1.875 mm，限于篇幅仅展示了沙粒平均粒径最大的模型，见图1。

图2 沙粒不同粒径模拟过程Fig.2 Simulation process of different particle sizes of sand

2 现场试验及接触参数标定

2.1 现场试验

在甘万铁路昌吉高勒站和巴音杭盖站之间的风沙严重区段，开展了轮轨力与钢轨动位移现场试验，具体测点位置如图3所示。现场试验轮轨垂向力及钢轨垂向位移测点各布置2个，其中轮轨垂向力A 测点布置在1 号和2 号轨枕之间，B 测点布置在2 号和3 号轨枕之间；轨枕垂向位移2 个测点均布置在跨中位置。测得试验列车C70E 速度为45 km/h时，1号位置钢轨垂向力最大为129.46 kN，最大垂向位移为0.686 mm；2 号位置垂向力为128.86 kN，最大垂向位移为0.676 mm。图4 所示为1号位实测轮轨垂向力时程曲线。

图3 现场测点布置简图Fig.3 Schematic diagram of field test

图4 实测轮轨垂向力Fig.4 Measured wheel-rail vertical force

2.2 颗粒流模型微观参数标定

为获得数值试验接触模型的微观接触参数，利用自编fish语言将图4所示列车荷载时程曲线施加到图1 所示离散元模型轮轨接触区的加载颗粒上，模拟实际列车荷载作用。参考文献[16-20]中道砟颗粒接触参数，结合文献[15]中钢轨颗粒黏结参数，得到不同接触参数下钢轨垂向位移时程曲线，并选取与实测钢轨垂向位移最为接近的垂向位移仿真曲线进行参数标定，如图5 所示。由图5可知：钢轨垂向位移仿真结果与现场实测结果数值和规律均基本一致。

图5 钢轨垂向位移现场实测值与数值试验值对比Fig.5 Comparision of measured value and numerical test value of rail vertical displacement

为进一步验证模型参数的可靠性，统计了峰值荷载下钢轨垂向位移的实测值与仿真值，如表1所示。由表1 可知：实测钢轨垂向位移为0.675~0.678 mm，数值试验位移为0.625~0.657 mm，相对误差为5.24%。由此可见，钢轨垂向位移试验值与仿真值均有很高的一致性，即数值试验模型可以用于道床内尺度比研究。在离散元模型中，扣件颗粒间抗拉强度TF为5×109Pa，抗剪强度SF为1×108Pa；钢轨颗粒间法向黏结刚度-kn为1.4×1012N/m，切向黏结刚度为5.5×1011N/m，黏结半径为0.036 m，半径乘数为0.485。此外，其他接触参数如表2所示。

表1 钢轨垂向位移实测值与仿真值统计Table 1 Statistics of measured and simulated values of rail vertical displacement

表2 模型细观接触参数Table 2 Mesoscopic contact parameters of model

3 分析方案与计算效率

为研究沙粒多尺度效应对道床受力的影响，量化道床内部颗粒粒径尺度，提出以量纲一参数内尺度比反映道床颗粒体系的跨尺度效应。道床内尺度比符号记为k，表征道砟平均粒径与沙粒平均粒径的比值，其定义式如下：

图6 内尺度比与胞元数的关系Fig.6 Relationship between intrinsic scale ratio and cell number

散体材料利用离散单元法模拟时，都会面临研究问题的规模和计算耗时的问题。为提高计算效率，参考文献[3]中的沙粒尺度，通过试算分析确定合理的沙粒尺度计算区间。胞元个数与模型计算效率的试算结果如图7 所示。由图7（a）可知：随着胞元个数增加，道床中颗粒数量快速增长，计算耗时呈指数型增长；当x>6时，颗粒总数与计算耗时曲线斜率明显变大，计算时间急剧上涨，计算效率低。由图7（b）可知：当x<3时，道床内部峰值应力曲线急剧上升，波动很大；但当x>6 时，曲线斜率较小，数值逐渐收敛。这与道床内部颗粒接触密度、力链传递途径等密切相关。综合考虑沙粒尺度敏感性分析试验、单机计算效率和研究问题需要，建议道床颗粒内尺度比为57.52左右较合理。

为进一步表征沙粒尺度对道床峰值应力的影响，对图7 中峰值应力关系曲线进行非线性拟合，得出胞元个数与峰值应力函数关系式：

式中：σ 为峰值应力，σ0为当x →+∞时的极限应力，实际研究中x为有限值。拟合优度R2=0.997 7。拟合结果表明，随着胞元个数增加，沙粒粒径逐渐减小，峰值应力总会趋于临界值。

图7 计算效率优化分析Fig.7 Optimization analysis of computational efficiency

在上述数值试验模型分析的基础上，设计8种沙粒粒径的计算工况，进行道床宏细观力学特性数值分析，通过可变沙粒粒径形成的道砟空隙填充和物质结构离散性差异，了解风沙区道床是否具有内尺度特性，以确定合理内尺度比，见表3。

表3 不同沙粒尺度计算方案Table 3 Calculation scheme of different sand size

4 风沙道床沙粒内尺度比特性

4.1 道床细观力学特性

4.1.1 细观接触特性

散粒体道床内部道砟颗粒之间的接触力直接会影响其劣化及服役寿命[17-18]。为研究不同内尺度比下道床内部荷载传递机制及颗粒之间的接触力，更加直观反映道床内部接触力的分布规律，选取加载点的主影响区域进行分析，利用自编fish语言提取颗粒接触信息，绘制峰值荷载作用下不同内尺度比的颗粒法向接触力与切向接触力分布云图，分别如图8和图9所示。

由图8可知：道床中强力键都主要分布在轨枕正下方，且沿45°角向下扩散[19]，而轨枕两侧以弱力键为主。经对比可知：随着沙粒胞元个数增加，内尺度比逐渐增大，道床内部接触力较大区域逐渐减小，接触力较小区域逐渐增大，道床内部法向接触力趋于均匀化，弱力键的作用逐渐增强；随着道床深度增加，接触力较大区域面积逐渐减小；当内尺度比k≥57.52时，道床内部法向接触力较大区域面积基本保持不变。这表明沙粒粒径减小可以减小道床内部法向接触力，自五分胞体后，内尺度比对道床法向接触力分布影响减弱。

图8 不同内尺度比下道床内部法向接触力分布Fig.8 Normal contact force distribution in the ballast bed with different internal scale ratios

由图9可知：轨枕正下方切向接触力较大，接触力扩散角大于45°，且随着道床深度递增，切向接触力较大区域面积递减；随着沙粒胞元个数增加，内尺度比逐渐增大，切向接触力分布较大区域（强力键）逐渐减小，较小区域（弱力键）逐渐增大，道床中约90%以上切向接触力小于0.5 kN，这与道床的受载方向有密切关系。从强力键分布面积可知：内尺度比k≥57.52时，切向接触力较大区域面积低于5%，道床内部颗粒之间的切向接触力趋于均匀化。

综合图8和图9可知：峰值荷载作用下道床中颗粒之间的法向接触力是切向接触力的5倍，接触力传递规律基本一致，这也说明了沙粒粒径变化对法向接触力的影响更为显著。

图9 不同沙粒内尺度比下道床内部切向接触力分布Fig.9 Distribution of tangential contact forces in the ballast bed with different internal scale ratios of sand grains

除了接触力分布外，配位数也是表征颗粒细观接触特性的重要指标。配位数决定颗粒的受载形式和微观接触各向异性等。当配位数较大时，颗粒各个方向均有接触力键，颗粒的动力响应也变得复杂，趋于各向同性加载；当配位数小于等于3 时，颗粒极不稳定，接触各向异性显著增强。图10所示为8种数值试验的道砟颗粒平均配位数。由图10 可知：胞元个数愈多，内尺度比愈大，道砟颗粒的平均配位数逐渐增加，最大平均配位数为26.3，最小平均配位数为10.7，极差为15.6。这表明沙粒粒径变化会引起道砟颗粒的配位数变化，影响道砟颗粒受载形式及微观接触分布规律。

图10 道砟颗粒平均配位数Fig.10 Average coordination number of ballast

为进一步揭示胞元个数与道砟颗粒配位数的关系，建立非线性Logistic回归模型进行分析，如式（6）所示。

式中：P 为道砟颗粒平均配位数；f 为下收敛值，取10.04；g 为上收敛值，取30.24；P0为道砟平均配位数中值，与g和f相关；x为胞元个数，取离散值，但在研究连续粒径问题时，可取连续值。拟合优度R2=0.998 2。拟合结果表明：胞元个数与道砟平均配位数呈良好的非线性对应关系，也说明沙粒愈小，内尺度比作用越弱。由Logistic回归曲线的走势可知，道砟平均配位数存在收敛点，可预测在真实沙粒作用下，道砟颗粒平均配位数接近30.24。

颗粒内尺度比的差异会导致道床内部接触力分布和道砟颗粒平均配位数均存在差异，说明沙粒粒径的改变亦会影响道床内部颗粒的微观结构。

4.1.2 细观位移概率密度

为揭示内尺度比对道床内部颗粒细观位移的影响，利用LIU 等[21-22]提出的概率密度的统计方法，计算颗粒细观位移。细观位移的概率密度可以衡量颗粒位移出现在某个范围的概率。本文利用自编fish语言统计了颗粒细观位移，采用式（7）计算颗粒位移的概率密度。

式中：y 为颗粒位移，其取值范围为-∞＜y<＋∞。对散粒体系，其位移是离散且有限的，y取值视实际研究的问题而定；p（y）为概率密度函数。颗粒位移概率密度计算结果如图11所示。

颗粒细观位移可以反映道床宏观动位移，细观位移概率集中分布区决定了道床动位移的整体水平。由图11 可知：颗粒位移概率密度曲线大致呈“人”字型，出现较大位移的概率较小，颗粒位移大多集中在位移0点附近。道砟颗粒水平位移主要集中在-1.0~1.0 mm，垂向位移主要集中在-0.5~0.5 mm；沙粒水平位移主要集中在-2.0~2.0 mm，垂向位移主要集中在-1~1 mm，即沙粒位移约为道砟颗粒位移的2倍。道床内部颗粒水平位移关于“0”位移线对称分布，而道砟颗粒垂向位移概率密度峰值向与道床受载方向相反的正向移动，沙粒垂向位移也向正向移动。随着胞元个数增加，概率密度主峰值呈“V”字型变化，主峰值先增大后减小，说明沙粒粒径越小，道砟颗粒向下移动的概率越小，道床的力与位移非线性越强，道床宏观动位移会逐渐减小，并趋于稳定。

此外，当道床颗粒内尺度比k≥52.38 时，颗粒水平位移的概率密度曲线变化不大；当道床颗粒内尺度比k≥57.52时，道砟颗粒垂向位移概率密度曲线变化很小；但沙粒垂向位移概率密度曲线稳定性较差，变异性较强。这表明颗粒越小，越容易受扰动，道砟颗粒质量较大，其惯性越大，受载不易改变其力学性质，即道砟颗粒比沙粒稳定。改变沙粒粒径导致道砟的微观接触状态发生变化，进而引起道砟颗粒位移分布规律发生变化，道砟颗粒概率密度极值减小，趋于均匀化。这与道床内部道砟接触数量的增加及接触力的减小有一定关系。

4.2 宏观动位移

图11 荷载峰值作用下道床内部颗粒细观位移概率密度分布Fig.11 Probability density distribution of particle micro displacement in ballast bed under peak load

为建立细观颗粒位移与道床表层宏观动位移的联系，揭示沙粒粒径变化对道床宏细观位移的影响，选取与加载颗粒临近的2根轨枕，提取峰值荷载作用下轨枕的最大位移及平均位移，绘制荷载峰值下不同胞元个数的轨枕动位移变化三维柱状图，如图12所示。

图12 荷载峰值下轨枕动位移Fig.12 Sleeper dynamic displacement under load peak

由图12 可知：随着胞元个数增加，道床颗粒内尺度比增大，轨枕动位移都逐渐减小。左轨枕最大动位移为0.591 mm，最小动位移为0.321 mm；右轨枕最大动位移为0.512 mm，最小动位移为0.295 mm，由此可知：沙粒粒径变化会影响道床表层位移变化，进而引起轨枕动位移改变，且由于散体材料的离散性强，左右轨枕的下部接触个数和接触面积各不相同，导致不同沙粒粒径下轨枕动位移最大值存在差异。为更好地表征道床宏观位移随沙粒粒径的变化规律，拟合轨枕平均动位移，建立胞元个数与平均动位移的非线性Logistic 相关性分析模型，根据平均动位移变化速率将平均动位移曲线划分为非敏感区（A~B），敏感区（B~C~D）和收敛区（D~E），如图13所示。

由图13 可知：轨枕平均动位移与胞元个数呈负相关，胞元个数愈多，平均动位移愈小。此外，对平均动位移曲线进行拟合得到的Logistic相关性分析函数如下：

图13 轨枕平均动位移与胞元个数的内在联系Fig.13 Average dynamic displacement of sleeper under load peak

式中：S 为平均动位移；S0为敏感区C 点平均动位移；x为胞元个数，与道床内尺度比及沙粒粒径相关；a 为平均动位移下收敛值，取0.302 mm；b 为平均动位移上收敛值，取0.533 mm，b - a 为沙粒粒径主影响区。拟合优度R2=0.998 2，拟合结果表明沙粒胞元个数与平均动位移之间存在很强的非线性相关性，a反映了真实沙粒粒径下轨枕的平均动位移，b反映了沙粒粒径的放大倍数的限值，b-a代表了沙粒粒径对道床力学特性影响的合理研究区间。

此外，由图13 的Logistic 模型可知：当1≤x≤3时，平均动位移随沙粒粒径变化速率较慢，处于非敏感区；当3

5 结论

1）风沙铁路道床数值模拟中内尺度比对其计算结果有显著影响，内尺度比k≥57.52 是合理可靠的。

2）随着胞元个数增加，道床中颗粒数量快速增长，计算耗时呈指数型增长；当沙粒胞元个数x>6 时，颗粒总数与计算耗时曲线斜率明显变大，计算时间急剧增长，计算效率低。

3）沙粒粒径会影响道床内部微观力键分布。无论是道床内部颗粒法向接触力还是切向接触力均随内尺度比增大，强力键分布面积逐渐减小，且沙粒粒径对法向接触力的影响更为直观；当内尺度比k≥57.52时，道床内部接触力分布面积趋于稳定。道砟平均配位数存在收敛点，利用Logistic模型可预测真实沙粒作用下，道砟颗粒平均配位数接近30.24。

3）沙粒粒径变化对颗粒细观位移有一定影响。道砟水平位移主要集中在-1.0~1.0 mm，垂向位移主要集中在-0.5~0.5 mm；沙粒水平位移主要集中在-2.0~2.0 mm，垂向位移主要集中在-1.0~1.0 mm，即沙粒位移约为道砟颗粒位移的2倍。随着沙粒粒径减小，颗粒水平位移概率密度峰值逐渐减小，主概率密度分布区逐渐增大；垂向位移概率密度峰值向正位移方向移动。

4）道床宏观动位移幅值与沙粒粒径密切相关。轨枕平均动位移与内尺度比呈负相关，当x≥6 时，轨枕平均动位移过渡到收敛区，沙粒粒径对轨枕动位移的影响减弱。仅从宏观位移角度分析，真实沙粒粒径作用下轨枕动位移平均值不超过0.302 mm。