一种无砟轨道冻结与损伤行为计算方法与流程

本发明属于轨道冻结与损伤行为的技术领域，具体涉及一种无砟轨道冻结与损伤行为计算方法。

背景技术：

我国地域辽阔，有相当大的部分地区处于严寒地带，而该区域混凝土结构面临强烈的冻融循环作用。气象统计数据结果表明，兰州地区的年平均冻融循环次数达83.6次，黑龙江哈尔滨冻融循环次数约45次作用，即便江苏南京地区的年平均冻融循环次数也达到19.9次，冻融破坏是我国北方寒冷地区混凝土结构损伤失效的主要原因之一。

无砟轨道结构设计时，虽然考虑了轨道结构的排水及抗冻性，并做了相应的排水及抗冻性设计。但由于施工和养护不当，在雨水丰富地区和排水不畅地段导致有大量的水会沿着表面裂缝或层间离缝渗入轨道结构内部，即使排水较好的区域由于低温环境下轨道板积雪，温度升高融化将导致很大比例的水分沿孔隙及微裂隙渗流进入轨道结构内部。因此，北方寒冷地区部分无砟轨道已出现较为严重的冻结与损伤行为。其中，i型双块式无砟轨道主要由钢轨、扣件、预制轨枕、道床板、混凝土底座或支承层组成，在我国铺设广泛，其中大西、西宝、兰新、宝兰等严寒地带都大量采用了i型双块式为主的轨道结构形式。双块式无砟轨道在桥梁地段采用单元结构，路基和隧道地段采用连续结构。其轨枕为工厂预制，道床板、底座板或支承层采用现场浇筑，对桥梁、路基、隧道等线下基础的适应性较好，但由于为现浇结构，浇筑养护后极易产生干缩裂缝，水分更易渗入轨道内部，因此面临更为严峻的冻结与损伤行为威胁。

powers提出的静水压理论为混凝土冻融受力的量化分析奠定了基础，此后国内外学者展开了大量研究，陆续提出了渗透压理论(powers)、结晶压理论(scherer)、微结冰理论(setzer)、粘结剥落理论(valenza)以及基于热力学原理建立的联系宏观结构与微观受力的多孔介质力学理论(coussy)。在数值计算方面，t.ueda等人(2009)基于刚体弹簧方法(rbsm)，通过引入零强度单元及塑形拉伸弹簧，建立了混凝土中尺度计算模型，分析了冻融引起的应力演化、裂纹扩展、强度退化问题以及冻融后混凝土的受弯性能；linliu，weisun等人采用水泥浆体数值微观结构，基于离散化的格构单元法，建立了含量化微观损伤的三维格构模型，分析了饱水状态、外荷载及冻结共同作用下水泥浆体微观结构的受力及微裂纹的分布；qinglidai在细观角度采用扩展有限元(xfem)研究了过冷条件下孔内结晶压力对混凝土内裂纹扩展的影响。

上述研究均通过将经验或理论预测的孔隙压力以外荷载的方式引入到模型中，并未考虑材料应力与孔隙压力的耦合关系。基于水-热-力耦合多孔介质力学物理表达式，采用毛细孔内冰水压力平衡关系，b.zuber,j.marchand在微观尺度下研究了低温冻结过程中气孔尺寸及间隔尺寸对水泥基材料变形的影响；段安分析了混凝土试块受低温冻结过程中的温度场、应变场以及孔内压力场的分布规律；曾强研究了冻结速率及气孔内压力边界条件对水泥基材料变形的影响。而混凝土冻结与损伤行为研究在无砟轨道领域中应用尚浅，低温冻结与损伤行为对无砟轨道的影响规律尚不明确。

技术实现要素：

本发明的目的在于针对现有技术中的上述不足，提供一种无砟轨道冻结与损伤行为计算方法，以解决或改善上述的问题。

为达到上述目的，本发明采取的技术方案是：

一种无砟轨道冻结与损伤行为计算方法，其包括：

s1、根据动弹性波试验获取已浇筑混凝土试件的多孔骨架体积压缩模量，并计算得到biot系数；

s2、根据ct扫描和图像识别技术，获取混凝土试件孔结构分布参数和整体孔隙率；

s3、根据所述孔结构分布参数，计算得到混凝土试件冻结过程中温度-冻结速率的关系曲线；

s4、根据young-laplace方程冰-未冻水平衡以及冰-吸附水膜平衡关系，构建平均孔隙压力与温度、孔结构之间的函数关系；

s5、基于多孔介质力学、多孔介质内水分迁移达西定律和多孔体系热传导friour定律，构建多孔介质应力场-温度场-渗流场的耦合控制方程，得到冻融过程中无砟轨道内结构受力、热量传输和水分迁移之间的相互作用关系。

优选地，步骤s1中计算得到的biot系数b为：

其中，k0为多孔骨架压缩模量，km为固体基质压缩模量。

优选地，步骤s3中计算得到混凝土试件冻结过程中温度-冻结速率的关系曲线为：

其中，温度降低过程中水的质量结冰率，ρw为水的密度，wi为水的结冰质量；

vw-i为孔隙内冰的体积含量，其计算表达式为：

vw-i＝φ(rpeq)-vads(rpeq)

吸附层的体积含量vads为：

其中，为孔径分布函数，rn为结晶的最大孔隙半径，r为孔隙半径；

孔径rpeq为：

rpeq＝req+δ

其中，req为冰晶体尖端平衡半径；δ为孔隙水膜厚度之和，其表达式为：

其中，θ为摄氏温度；

冰晶体尖端平衡半径req为：

其中，γls为冰与水表面张力，t0为正常水冰点，rmax为冰晶体尖端最大半径。

优选地，步骤s4中构建平均孔隙压力与温度、孔结构之间的函数关系为：

p^*(θ)＝pw(θ)+x

其中，pi为孔隙冰压力，pw为孔隙水压力，γ为冰与水表面张力，req为冰晶体尖端平衡半径，πw为吸附层水压力，p^*(θ)为平均孔隙压力，pw(θ)为结冰孔，x为未结冰孔；

未结冰孔x的表达式为：

其中，n为整体孔隙率，rpeq(θ)为结晶的最小孔隙半径，req(θ)为冰晶体尖端平衡半径。

优选地，步骤s5中构建多孔介质应力场-温度场-渗流场的耦合控制方程分别为：

当不考虑体力的作用下，多孔介质微元的应力场控制方程为：

其中，为张量算子，c0多孔介质骨架的弹性刚度张量，ε为多孔骨架的应变，i为单位张量，由于温度变化是孔隙压力变化的直接原因，令：

其中，k0为多孔骨架压缩模量，α0为多孔骨架的体膨胀系数；

忽略对流对热传导的影响，并考虑相变潜热，则多孔介质的温度场控制方程为：

其中，ρm为多孔骨架密度，为变温速率，λ为热传导系数，l为水的相变潜热，多孔体系比热容为cef：

其中，cm、cw和cw分别表示多骨架、水以及冰的比热容，λef、λm、λw、λi分别为多孔体系、骨架、水以及冰的热传导系数，sw为水饱和度，ρw为水密度，si为冰饱和度，ρi为冰密度；

混凝土冻结过程中的内水分迁移，即渗流场的控制方程为：

其中，k为渗透系数，μ为水的动力粘滞系数，为体应变，s为压力源项，

其中，β为多孔弹性储水修正系数，kw为水的压缩模量，ki为冰的压缩模量；

其中，αw、αi、α0分别为水、冰、多孔骨架的体膨胀系数。

本发明提供的无砟轨道冻结与损伤行为计算方法，具有以下有益效果：

本发明通过多孔介质内的渗流、传热过程、应力场本构及相互耦合作用，考虑孔隙结构对水分冻结温度的影响，并基于弹性多孔介质力学，研究了冻融过程中无砟轨道内结构受力、热量传输和水分迁移之间的相互作用规律，为无砟轨道混凝土冻结与损伤行为问题的应对及防治工作提供理论指导。

附图说明

图1为孔结构分布曲线。

图2为双块式无砟轨道横截面示意图。

图3为pastei试验与数值结果。

图4为pasteii试验与数值结果。

图5为水压力空间分布。

图6为支承层孔内水压力演化。

图7为道床板孔内水压力演化。

图8为道床板和支承层结晶压力。

图9为水压力源项组成及演化。

图10为支承层结冰速率。

图11为道床板结冰速率。

图12为道床板和支承层冰饱和度。

图13为温度空间分布。

图14为支承层降温速率。

图15为道床板降温速率。

图16为tresca应力空间分布规律。

图17为支承层tresca应力演化规律。

图18为道床板tresca应力演化规律。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

根据本申请的一个实施例，参考图1，本方案的无砟轨道冻结与损伤行为计算方法，包括：

s1、根据动弹性波试验获取已浇筑混凝土试件的多孔骨架体积压缩模量，并计算得到biot系数；

s2、根据ct扫描和图像识别技术，获取混凝土试件孔结构分布参数和整体孔隙率；

s3、根据所述孔结构分布参数，计算得到混凝土试件冻结过程中温度-冻结速率的关系曲线；

s4、根据young-laplace方程冰-未冻水平衡以及冰-吸附水膜平衡关系，构建平均孔隙压力与温度、孔结构之间的函数关系；

s5、基于多孔介质力学、多孔介质内水分迁移达西定律和多孔体系热传导friour定律，构建多孔介质应力场-温度场-渗流场的耦合控制方程，得到冻融过程中无砟轨道内结构受力、热量传输和水分迁移之间的相互作用关系。

以下对上述步骤进行详细说明

步骤s1、根据动弹性波试验获取已浇筑混凝土试件的多孔骨架体积压缩模量，并计算得到biot系数b为：

其中，k0为多孔骨架压缩模量，km为固体基质压缩模量。

步骤s2、通过ct(x射线或核磁共振)扫描及图像识别技术，获取混凝土试件孔结构分布参数及整体孔隙率n。

步骤s3、计算得到混凝土试件冻结过程中温度-冻结速率的关系曲线，其具体包括：

根据gibbs-thomson方程，孔径与结晶温度的关系表示为：

其中，tf为孔隙水冰点(k)，t0为正常水冰点(273.15k)，γls为冰与水表面张力(39×10^-3n/m)，ρw为水的密度(1000kg/m³)，l为水的相变潜热(333.5kj/kg)，req为温度降至tf时的冰晶体尖端平衡半径(nm)。

在一般环境下，

将其进行tylor展开，可简化为：

其中，θ-孔隙水冰点(℃)，带入相关数值后可得：

其中，64为计算常数(℃·nm)。

在已结冰的孔隙内，孔隙壁与冰晶之间式中存在一层水膜，孔隙水膜厚度可表示为：

其中，1.93为计算常数℃^1/3·nm，因此，温度为θ时，孔径大于rpeq(nm)的孔隙才会结冰：

rpeq＝req+δ

假设孔隙中只有水和冰，混凝土完全饱水，则

si+sw＝1

其中，si、sw分别表示冰和水的饱和度，当温度为θ时，假设孔隙为半径r的理想圆柱形，吸附层所占孔隙体积含量为2δ/r，因此结晶孔隙中，吸附层的体积含量为：

孔隙内冰的体积含量vw-i(m³/m³)表示为：

式中，为孔径分布函数，rn为结晶的最大孔隙半径。由此可得，温度降低过程中水的质量结冰率可表示为：

步骤s4，根据young-laplace方程冰-未冻水平衡以及冰-吸附水膜平衡关系，构建平均孔隙压力与温度、孔结构之间的函数关系，其具体步骤包括：

采用高斯函数表示孔结构分布：

其中，为大于半径r的累积孔隙比，r为孔隙半径，f0、ai、bi及ri为待定系数。

孔隙压力，假设在混凝土内孔隙为圆柱形，冰-未冻水接触面为球面，通过young-laplace方程冰-未冻水平衡以及冰-吸附水膜平衡关系分别表示为：

施加在混凝土微元孔隙内的孔隙压力可表示为微元内结晶孔隙内吸附水膜压力以及未冻结孔隙内的水压力：

步骤s5、基于多孔介质力学、多孔介质内水分迁移达西定律和多孔体系热传导friour定律，构建多孔介质应力场-温度场-渗流场的耦合控制方程，得到冻融过程中无砟轨道内结构受力、热量传输和水分迁移之间的相互作用关系，其具体步骤包括：

s5.1、应力场控制方程；

对于各项同性的多孔体系，其应力可分为两部分：多孔体系骨架的有效应力和孔内压力：

σ＝σ^*-bp^*i

σ^*＝c0：ε^θ

其中，σ^*为骨架有效应力，p^*为孔内压力，i为单位张量，ε^e多孔介质弹性应变，c0为多孔介质骨架的弹性刚度张量，b为boit系数表征孔隙压力对变形的贡献，对于各项同性材料，一般表示为：

其中，k0和km分别表示多孔骨架以及固体基质压缩模量。

假设骨架没有塑形应变，考虑到热膨胀效应，多孔骨架的应变可分为温度应变ε^th以及弹性应变ε^e：

ε＝ε^e+ε^th

而骨架的温度应变可表示为：

联立式可得：

由于温度变化是孔隙压力变化的直接原因，令

则：

当不考虑体力的作用下，多孔介质微元的力学平衡方程可表示为：

s5.2、渗流场控制方程；

假设多孔体系内水封迁移服从darcy渗流定律，由质量守恒及冰、水、骨架本构方程，混凝土冻结过程中的内水分迁移表示为：

其中，k为渗透系数(m²)，μ为水的动力粘滞系数(pa/s)；εv为体应变；kw、ki分别为水和冰的压缩模量(pa)；αw、ai、α0分别为水、冰以及多孔骨架的体膨胀系数。

s5.3、温度场控制方程；

忽略对流对热传导的影响，并考虑相变潜热，则多孔介质的热传导方程表示为：

其中，cef、cm、cw和ci分别表示多孔体系、骨架、水以及冰的比热容(j/(kg·k))；λef、λm、λm和λi分别表示多孔体系、骨架、水以及冰的热传导系数(w/(m·k))；l为水的相变潜热(j/kg)。

根据本申请的一个实施例，本发明基于构建的多孔介质应力场-温度场-渗流场的耦合控制方程，构建i型双块式无砟轨道水-热-力耦合平面模型，模拟大气温度在10h内由10℃等速率降低至-20℃然后在-20℃持续2h过程中，正常状态(h＝12w/(m²·k))和做简单保温隔热处理(h＝6w/(m²·k))工况下的冻结与损伤行为特征。

参考图1，图中位于上方的折线为pastei，下方的折线为pasteii，模型中支承层和道床板分别采用pastei和pasteii的孔结构，其中pastei的三个孔隙分布峰为327000nm、15.49nm、1.03nm，pasteii的三个孔隙分布峰为11.35nm、6.99nm、3.24nm，因此pastei孔结构更分散，大孔更多，孔隙率更大。

参考图2，模型中各设置5个数值监测点，下文均表示为点i(i＝1～10)。

数值计算时，模型中采用的参数详见下表2-1：

表2-1数值模拟参数

本发明通过有限元方法对建立的混凝土冻结与损伤行为水-热-力耦合理论模型进行求解。计算通过模型验证，得到了整个冻结过程中不同温度状态下的混凝土的应变，并于试验结果相对比。

计算结果并对其分析

1、水压力分布特征

混凝土冻结过程中，未结冰孔隙水迁移滞阻将产生孔隙内部水压力，较大的孔隙内压力将引起混凝土内部微元解体。因此，孔隙水压力是研究混凝土冻结与损伤行为的关键指标。不同工况下无砟轨道内部孔内水压力分布如图5所示。

由图5可知，最大水压力在道床板板角向内部一定深度处，且道床板水压力远大于支承层内部水压力，支承层和道床板各观测点在冻结过程中，水压力的演变规律如图。

由图6和图7可知，由于冻结阻滞效应，轨道板冻结区渗透系数迅速降低，随深度增加水压力有巨大的提升。道床板的水压力增长峰值分别在-5℃及-10℃附近，其分别对应的孔分布峰为11.35nm、6.99nm。支承层的水压力增长峰值在0℃和-5℃附近，其对应的孔分布峰为327000nm、15.49nm。相较于道床板，支承层内水压力随深度增大明显提早，由点1处的0℃提早到点6处的3℃，这并不是由于冻结温度提高，而是支承层渗透系数较小冻结区产生的水压力快速传递到轨道内部引起的。

支承层冻结区域内各点的水压力在h＝6w/(m2·k)工况下相较于h＝12w/(m2·k)约有0.5mpa的减小；而道床板冻结区域内各点的水压力在h＝6w/(m2·k)工况下相较于h＝12w/(m2·k)也约有1mpa的减小。因此，保温隔热处理可以降低轨道内部水压力，最主要的是有效延缓冻结过程，减少无砟轨道的冻结深度，从而有效降低无砟轨道内的水压力空间分布大小。

由公式可知，结晶压为温度和孔结构的函数，支承层和道床板内结晶压与温度的关系如图9所示

由图可知，道床板在-10～15℃区间内结晶压力有增长较快。支承层在0～-5℃区间内结晶压增长较为缓慢且明显小于道床板，-5℃以后开始迅速增大。

式中的s表示水压力的来源包含四个相关部分：第一项与结冰速率相关，第二项与热膨胀效应相关，第三项平均孔压相关，第四项和结晶压力有关。其中，支承层内水压力源s的四项随着温度下降的变化如图9所示。

由图9可知，在冻结过程中，结冰速率占据水压力源的主导项，尤其是在冻结初期，随着温度逐渐降低，结冰速率所占比重逐渐降低，但仍然占据主导地位。因此，有必要研究支承层和道床板在冻结过程中结冰速率的演化规律。支承层和道床板各观测点在冻结过程中，结冰速率的演变规律如图10及图11所示。

由图10及图11可知，随深度增大，点内结冰速率逐渐减小，结合图6、图7中水压力随冻结深度的增加而变大的关系，揭示出冰阻导致渗透系数减小对水压力的巨大影响。在冻结过程中，在0℃很小温度区间内结冰速率最大。此后支承层在-3℃附近有一个冻结峰，道床板在-6℃和-10℃同样存在冻结峰。在冻结峰内道床板和支承层最大结冰速率分别为0.0035、0.013，结合图6、图7可知道床板水压力更大，由此也反映出渗透系数对水压力的影响大于结冰速率。支承层和道床板在h＝6w/(m2·k)工况下冻结区域内的结冰速率约为h＝12w/(m2·k)的工况的50％，因此保温隔热处理可有效减小轨道内结冰速度率。

由图12可见支承层内孔隙水在-10℃后结晶速率迅速降低，最终冰饱和度约为0.35，而道床板由于孔结构较为致密，在-10～-15℃区间内冰饱和度仍有显著增大，其最终冰饱和度约为0.55的试件，这是结冰速度更快。

2、温度分布特征

温度是引起混凝土冻结与损伤行为的最直接因素，温度分布规律可以间接判断混凝土的冻结与损伤行为深度。对于不同的表面热传导系数条件下，无砟轨道的最终温度场如图13所示。

由图13可知，在道床板和支承层板角处冻结最深，保温隔热处理后轨道表面的温度略有提高。h＝6w/(m2·k)和h＝12w/(m2·k)工况下轨道表面温度分别为-10℃和-15℃。其中各监测点的温度变化速率分别如图14及图15所示。

由图14及图15可知，随着深度的增加，轨道内降温速率也越来越小。由于孔隙内水化热的释放，在0℃附近，变温速率有非常大的波动。同时支承层和道床分别在-3℃和-10℃附近区间内降温速率变小，分别对应于图10和图11的结晶峰。保温隔热处理后轨道内部降温速率有显著的减小。

3、应力分布特征

冻结过程中，混凝土孔内压力作用于固体骨架，为分析轨道内的应力分布特征，不同工况下无砟轨道的tresca应力由图16所示。

由图16可知，由于支承层弹模较低，密实度较小，boit系数更大，孔内水压力更多的反映到混凝土骨架上，道床板内最大剪应力出现在道床板表层和深度约为5cm处。该结果很好的解释了现场双块式无砟轨道表面冻结与损伤行为剥落情况。而支承层的最大tresca应力则出现在板角表面和距板角约5cm处。各观测点在不同工况下的tresca应力演化规律分别如图17及图18所示。

由图17及图18可知，在0℃以下时，随着温度的降低支承层内tresca应力显著增大，距离表面越深tresca应力增长越快。h＝12w/(m2·k)和h＝6w/(m2·k)工况下，支承层最大tresca应力分别为5.5mpa和2mpa，而道床板最大tresca应力分别为5mpa和4mpa，保温隔热处理可减弱冻结与损伤行为对无砟轨道的破坏力度。由于支承层内水压力传递速度较快，tresca应力也体现出明显的提前性。

根据以上描述，综合可得：

严寒地区无砟轨道受冻结与损伤行为影响较大，长期服役过程中不断出现表面剥蚀的状况。本发明基于多孔介质力学，建立了i型双块式无砟轨道水-热-力耦合模型，分析了冻结作用下无砟轨道应力-应变分布特征、温度以及孔隙压力的分布状态。得出主要结论如下：

(1)结冰速率占据混凝土内冻结水压力源的主导地位，由于水分结冰的阻止效应致使渗透系数降低，距离轨道表面越深的观测点在相同冻结温度下水压力峰值越大。

(2)混凝土的孔结构显著影响轨道内部水分的结冰速率及水压力演化规律。道床板孔径主要集中于3～15nm，冻结峰温度为0℃、-6℃及-11℃，结冰速率较为平缓。而支承层内部大孔结构较多，孔结构分布峰值对应的孔径为300000nm、15nm及1nm，结晶峰温度为0℃及-4℃，冻结速率更快。

(3)支承层点1在-5℃结冰峰结冰速率为0.003kg/(m3·s)，水压力为1mpa。而道床板孔隙率小，渗透系数为支承层的1/30，点6在-12℃结冰峰的结冰速率为0.0007kg/(m3·s)，水压力为8.5mpa。孔隙率及渗透系数对水压力的影响更显著。

(4)虽然道床板水压力更大，但支承层孔隙率更高，boit系数更大，孔内水压力更多的反映到混凝土骨架上，应力更大，受冻结与损伤行为最大。

(5)保温隔热层处理可有效减小轨道的冻结速率及深度，在h＝6w/(m2·k)和h＝12w/(m2·k)工况下，轨道表面温度分别为-14℃、-18℃，最大水压力分别为6mpa、2mpa。因此可通过无砟轨道两侧堆道砟的方法减轻无砟轨道的冻结与损伤行为。

虽然结合附图对发明的具体实施方式进行了详细地描述，但不应理解为对本专利的保护范围的限定。在权利要求书所描述的范围内，本领域技术人员不经创造性劳动即可做出的各种修改和变形仍属本专利的保护范围。