高速公路路基长期THM耦合模拟性能研究

Long-Term THM Coupling Simulation Performance Study of Highway Subgrade

摘要: 为准确掌握多场耦合作用下寒区路基长期使用过程中冻胀变形规律，基于路基实际工程地质情况，通过建立THM耦合模型，使用COMSOL软件对路基长期使用性能进行预测。结果表明：经过十年左右路基结构基本达到稳定状态；路基温度变化随着深度的增加呈现明显的滞后效应，在路基达到稳定之前，路基的高低温峰值逐年减小；路基融沉位移小于冻胀位移，且两者逐年递增，另外路中位移大于左右路肩位置；路基冻结深度在非稳定期内增加了80%左右。

Abstract: In order to accurately grasp the law of freezing and deformation during the long-term use of subgrade in cold areas under the action of multi-field coupling, based on the actual engineering geology of the subgrade, the THM coupling model was established and the long-term performance of the subgrade was predicted by using COMSOL software. The results show that: the subgrade structure basically reaches a stable state after about ten years; the temperature change of the subgrade shows obvious hysteresis effect with the increase of depth, and the high and low temperature peaks of the subgrade decrease year by year before the subgrade reaches a stable state; the thawing and sinking displacement of the subgrade is smaller than the freezing and swelling displacement, and it increases year by year, and in addition, the middle of the road is displaced more than the left and right shoulders; the freezing depth of the subgrade increases by about 80% in the non- stable period.

1. 引言

冻土是指土壤或岩石中水分冻结形成的固体状态[1]。在寒冷的气候条件下，土壤中的水分会在冬季结冰并形成冻土层，这一过程称为冻结[2] [3]。随着气温升高，冻土层会逐渐解冻，土壤中的水分重新变为液态，这一过程则称为融解[4] [5]。冻土的冻结和融解过程对土壤结构、水文循环、植被生长等方面都产生重要的影响[6] [7]。研究寒冷地区路基的温度–渗流–固体(THM)耦合作用对路基长期使用性能的影响可为公路运行安全保障工作提供决策依据。

针对冻土多场耦合的THM耦合研究，有很多国内外学者对冻土的冻胀机理进行分析。Guymon、Mu等[8] [9]通过建立冻胀模型对冻土冻胀进行数值模拟，表明数值方法对解决冻胀问题有良好的适用性。毛雪松、白青波等[10] [11] [12]通过建立耦合模型研究多年冻土路基中温度场与水分场的分布和变化规律。郑洪瑞等[13]通过能量守恒方程开发了一个完全耦合的THM模型，得知温度的升高会增加土体的塑性变形。Nishimura等人[14]通过使用耦合的THM模型，用于描述涉及冻融土壤的各种地质和岩土工程过程。Li等人[15]通过对增加毛细排水板路基进行室内实验和THM耦合模拟，对路基的长期使用性能进行预测。王铁行等人[16] [17]通过冻土路基水热耦合模型分析了冻土路基的应力及应变分布。凌贤长等[18]通过非饱和膨胀土水–热–变形耦合冻胀模型，证明了膨胀土吸水膨胀特性对其冻胀特性有显著影响。

由于不同的研究对象所处的工程环境差别很大，多场耦合效应对工程的影响也不尽相同。基于G5京昆高速公路汉中至广元段A1标段实际工程情况。通过建立THM耦合方程，使用COMSOL模拟软件PDE模块进行数值建模，对该工程路基冻胀特性以及长期使用性能进行长期预测，分析长期使用过程中路基温度、位移以及冻结深度变化趋势，为路基长期使用安全提供评价依据。

2. 工程背景及水文条件

2.1. 工程概况

研究区位于四川省东北部，北起旺苍县天星乡青峰村，经王家垭、新农村，南至旺苍县燕子乡双河村，走廊带地理坐标：东经106˚6'58.31''~106˚16'54.28''，北纬32˚25'45.59''~32˚34'50.37''，呈北东–西南向条带状展布。项目名称为G5京昆高速公路汉中至广元段A1标段，是国家运输大通道北京至昆明高速的一段，也是四川省高速网北向出川、沟通陕甘中原的重要通道，地理位置如图1所示。路线全长35.099公里，项目全线按双向六车道高速公路标准建设，设计速度采用100 km/h，路基宽度为34 m。

Figure 1. G5 Beijing-Kunming expressway Sichuan expansion line plan schematic diagram

图1. G5京昆高速公路四川境扩容线路平面示意图

2.2. 水文气候条件

根据旺苍县气象资料，年均气温16.2℃，多年平均降水量1209 mm，最大降水量为1981年的2092.4 mm，最少为1979年的728.8 mm。年降雨量在地域上分布不均，县域东南部降雨量较大，西北部较小；降雨集中在5月至9月，占全年总降水量的80.8%。其中月均降雨量最高为7月，降雨量达279.3 mm；最低为1月，降雨量仅8.9 mm。其气候条件如表1所示。

Table 1. Meteorological data for Wangcang County

表1. 旺苍县气象资料

名称	最低温度(℃)	最高温度(℃)	夏季平均气温(℃)	冬季平均气温(℃)
旺苍县	−5.3	27	22.25	−3.4

3. THM耦合理论模型

基于THM耦合的主要物理过程，提出以下几种假设：① 路基土体为不可压缩、均匀、各向同性介质。② 路基土中只有液态水的运动，忽略了水蒸气迁移对冻水和热流迁移的影响。③ 忽略了温度梯度和对流对水分迁移的影响。在此基础上，建立了季节性不饱和冻土中的THM耦合模型。

3.1. 温度场方程

在温度场理论中主要考虑温度在土质中的热传递以及当气温较低时由冰水相变引发的潜热变化，其主要方程如下所示：

$\rho c\left(\theta \right)\frac{\partial T}{\partial t}=\lambda \left(\theta \right){\nabla }^{2}T+L{\rho }_i\frac{\partial {\theta }_i}{\partial t}$ (1)

式中，c(θ)为土壤的有效体积热容(J/(g∙℃)；λ(θ)为土壤的导热系数(W/(m∙k))；L为冰水相变的潜热(kJ/kg)；θ_i为孔隙中冰的体积(%)；ρ和ρ_i分别路基土和孔隙冰密度(g/cm³)。

3.2. 渗流场方程

在THM耦合过程中，水分场主要涉及路基之中自由水的迁移以及冰水相变问题。在不饱和冻土中的水分迁移定律可以用具有相变的理查兹方程来表示。

$\frac{\partial }{\partial x}\left(D\frac{\partial {\theta }_{\mu }}{\partial x}\right)+\frac{\partial }{\partial y}\left(D\frac{\partial {\theta }_{\mu }}{\partial y}\right)=\frac{\partial {\theta }_{\mu }}{\partial t}+\frac{{\rho }_i}{{\rho }_w}\frac{\partial {\theta }_i}{\partial t}$ $$ (2)

考到孔隙冰对路基土中水分迁移的阻碍作用，引入阻碍系数I。

$\begin{array}{c}D\left({\theta }_{\mu }\right)=\frac{K\left({\theta }_{\mu }\right)}{C\left({\theta }_{\mu }\right)}\cdot I\\ I={10}^{-10{\theta }_i}\end{array}$ (3)

式中，D为路基土中水分迁移扩散系数(cm²/s)。D(θμ)为未冻水扩散系数(m²/s)，K(θμ)为非饱和路基土在重力加速度方向上的渗透速度(m/s)，C(θμ)为比水容量(1/m)。

3.3. 应力场方程

在应力场中，假定路基土为各向同性的线弹性材料，使得土中的应力应变满足平衡微分方程。控制方程如下：

$\nabla ·\left[\sigma \right]+\left[F_{\nu }\right]=0$ (4)

其本构模型为

$\left\{\sigma \right\}=\left[C\right]·\left(\left\{\epsilon \right\}-\left\{{\epsilon }_0\right\}\right)$ (5)

则应力场中微分控制方程可表示为

$\nabla ·\left[\left[C\right]\left(\left[\epsilon \right]-\left[{\epsilon }_0\right]\right)\right]+\left[F_v\right]=0$ (6)

式中：$\nabla$ 为微分算子，σ为应力，F_v为物理矩阵，ε为应变，C为广义胡克定律中的弹性常数，ε₀为初始应变。

3.4. 耦合项方程

为了使应力场与温度场，渗流场进行耦合，使用相对饱和度S代替θ_μ作为水–热耦合过程的变量。路基土的相对饱和度(S)可表示为：

$S=\frac{{\theta }_u-\theta {}_r}{{\theta }_s-\theta {}_r}$ (7)

引入固液比B_I建立土体温度与含水量的动态平衡关系作为耦合方程，固液比表达式为：

$B_I=\frac{{\theta }_i}{{\theta }_{\mu }}=\{\begin{array}{l}1.1{\left(\frac{T}{T_f}\right)}^b-1,\begin{array}{cc}& T<T_f\end{array}\\ 0,\begin{array}{cc}& \begin{array}{cc}& \begin{array}{cc}& \begin{array}{cc}& T\ge T_f\end{array}\end{array}\end{array}\end{array}\end{array}$ (8)

式中，θ_r，θ_s分别为路基土剩余体积含水率和饱和体积含水率(%)。T_f为初始冻结温度(℃)。b是与含盐量和土壤类型相关的常数参数。

4. 路基长期使用性能模拟

4.1. 模型定义

根据A1标段路基实际尺寸，定义路基宽度为34 m，路基高为3 m，边坡坡度1:1.5，路面铺设18 cm厚改性沥青混凝土路基作为路面，路面以下分别为58 cm厚水泥稳定碎石基层以及224 cm厚路基填料。路基基底以下取10 m深度，地下水位线设定为路基基底以下7 m深，路基基底两侧各向外取10 m。路基填料为砂砾，地基地下水位以上7 m为粉质粘土，地下水位以下为饱和粉质粘土，路基模型简如图2所示。

(a) 路基模型示意图

(b) 数值模拟网格划分

Figure 2. Schematic diagram of subgrade model

图2. 路基模型示意图

4.2. 边界条件及参数取值

4.2.1. 边界条件

PDE模块在季节性多年冻土区的温度边界条件可表示为：

$T_n=T_0+A\sin \left(\frac{2\pi t}{365}+\frac{\pi }{2}+{\phi }_0\right)$ (9)

式中，T₀为旺苍县年平均气温(℃)，A为旺苍县气温年差(℃)，a₀为相位角。底部热边界条件选用温度梯度0.03℃/m，对应热流q为0.0375 W/m²。其温度边界条件如表2所示。

Table 2. Temperature boundary conditions

表2. 温度边界条件

温度边界参数	T0	A	φ0	Qb [w/m^2]
取值大小	12.5	17.8	7π/9	0.0375

渗流场中根据实际地质调查模拟地下水环境，饱和粉质粘土初始饱和度设为0.85，非饱和粉质粘土初始饱和度设为0.4，其余路基层初始饱和度设置为0.3。固体力学模块中应力场和应变场的边界条件设置为上边界自由约束u_x ≠ 0，u_y ≠ 0，下边界固定约束u_x = 0，u_y = 0，左、右肩边界定向约束，初始值为u_x = 0，u_y ≠ 0。

4.2.2. 参数取值

对于路面层，只考虑温度场参数对路基长期性能的影响，固体力学场以及渗流场参数直接引用填方路基层。由于路基填料层主要由泥岩，页岩和砂岩等破碎岩粒组成。所以路面层固体力学场与渗流场主要参考砂岩土。结合相关资料[19] [20]，各土层温度场、渗流场以及固体力学场的基本参数取值如表3~6所示：

Table 3. Basic parameters of temperature field

表3. 温度场基本参数

土层类型	深度(H) (m)	热容(C) (J/kg∙℃)	导热系数(λ) (W/(m∙k))
路面层	2.26~3	2450	0.669
填方层	0~10	1238	0.61
饱和粉质粘土层	−7~−3	1568	0.41

Table 4. Parameters of seepage model for unsaturated soil

表4. 非饱和土渗流模型参数

土层类型	a0 [1/m]	m	l	θs	θr	Ks [m/s]
填方路基层	0.66	0.14	0.5	0.231	0.02	1e−5
粉质粘土层	2.59	0.22	0.5	0.375	0.05	1e−8
饱和粉质粘土层	2.59	0.22	0.5	0.375	0.05	3e−8

Table 5. Mechanical basic parameters

表5. 力学基本参数

土层类型	密度(ρ) (g/cm^3)	弹性模量(e) (MPa)	泊松比(μ)
路面层	1200	16,000	0.25
填方路基层	1600	150	0.25
粉质粘土层	1540	80	0.35
饱和粉质粘土层	1540	50	0.45

Table 6. Parameters of unfrozen water content

表6. 未冻结水含量参数

土层类型	深度y (m)	起始冻结温度(Tf) (℃)	经验系数(b)	孔隙比(e)
填方路基层	0~2.26	−0.12	0.61	0.78
粉质粘土层	0~−10	−0.54	0.56	0.6

4.3. 模拟结果

4.3.1. 路基温度变化

Figure 3. Phase field of temperature change of roadbed during freezing period

图3. 冻结期路基温度变化相场图

截取路基在第一年温度降低期间温度变化情况如图3所示，在自然温度周期下第一年路基经过135天到达温度最低点。在环境温度降低期间冷能由路基表面以及路基坡面进行侵入，路基温度逐渐下降。其中，路基以外地表以及路基边坡处温度侵入速度大于路基表面处。主要为路基最上层铺设有路面层，热熔较大，在相同环境温度下温度变化较为迟缓，对路基整体温度变化起到一定阻碍作用。

4.3.2. 路基温度长期变化趋势

Figure 4. Temperature change trend graph

图4. 温度变化趋势图

通过对路基表面施加温度边界条件。模拟路基在15年中的长期使用性能。路基中轴处15年不同深度的温度变化趋势如图4所示，路基表面温度变化与边界温度变化趋势一致，路基温度波动程度较大，路基受边界温度影响较大。随着深度的增加，路基达到高低温峰值时间相比于浅层路基有明显的滞后效应。其中，低温侵入路基导致的温度变化滞后效应强于高温侵入。在距离路基表面3米深处的路基达到高低温峰值的时间相比于表层路基分别滞后了约90天与120天，随着距离路基表面越深其温度滞后效应越明显。其次，随着模拟时间的延长，路基内部达到高低温峰值的温度值在逐渐减小，经过10年左右，路基内部温度变化才趋于稳定。其中深度为3 m处的路基10年后路基最高温度相比于第一年降低了2.68℃，最低温度降低了2.21℃。而0.5 m深度处的路基经过10年最高温度下降了0.89℃，最低温度下降了1.45℃。随着深度的增加，路基能达到的高低温峰值减小的幅度越大。其主要原因为随着温度循环侵入的影响，路基孔隙中的自由水在温度较低时冻结膨胀，对路基基体产生应力影响。增强了路基骨架的连通性，结构损伤区域逐渐深入路基内部，导致路基内部对于温度变化的敏感性增强，峰值温度会逐年降低。经过10年左右，路基结构趋于稳定，路基内部温度变化也随之趋于稳定。

4.3.3. 路基位移变化

截取路基在第一年温度降低期间路基位移变化情况如图5所示，随着环境温度对路基内部的侵入，路基中的自由水在达到冻结温度时进行冻结，进而体积产生膨胀在宏观上表现出路基表面产生冻胀位移。由云图可知路基边坡首先会产生冻胀位移，其次为路基以外地表处。特别是路基边坡坡顶处位移变化随环境温度变化最为敏感。主要原因为环境温度由路基表面以及路基边坡进行侵入，路基边坡坡顶处温度变化更加敏感，当冷能侵入路基时自由水发生冻结，路基内部自由水沿横向与竖向向路基边坡坡顶处迁移，造成路基坡顶处冻胀位移值相对较大。

Figure 5. Phase field of displacement change of roadbed during freezing period

图5. 冻结期路基位移变化相场图

4.3.4. 冻胀融沉位移及冻深长期变化趋势

路基顶部的最大冻胀位移以及最大融沉位移如图6所示，在前10年时间内，路基中部以及路基的左右路肩处每年的最大冻胀位移以及最大融沉位移均缓慢增加。并且同一时间路中最大冻胀位移与最大融沉位移均大于左右路肩处。主要原因为在路基中部，当冻结锋面在路基内部移动时地下水更容易向冻结锋面进行迁移，导致位移值相对于左右路基处较大。另外，路中位移整体大于左右路肩处位移，前10年时间内，路中最大冻胀位移由3.13 mm增加至4.82 mm，平均位移为4.12 mm。路基左右肩位移变化趋势基本一致，最大冻胀位移由2.78 mm增加至4.11 mm，平均位移为3.54 mm。路中位置最大融沉位移由−0.67 mm增加至−1.68 mm，平均值为−1.18 mm。左右路肩处的最大融沉位移由−0.52 mm增加至−1.18 mm，平均位移为−0.85 mm。整体来说，融沉位移值相比于冻胀位移值较小，对路基的使用性能影响较低。经过10年时间，路基结构趋于稳定，路基受温度影响产生形变基本稳定。

Figure 6. Trend of freezing and thawing displacement as well as depth of freeze

图6. 冻胀融沉位移以及冻深变化趋势

路基在长达15年的使用过程中，路基冻结深度变化如图6所示。随着时间的增加，路基冻结深度逐渐增加，并在第10年达到稳定，这与路基的温度变化与位移变化基本一致。路基冻结深度由第一年的0.3 m增加至第十年的0.54 m，并在第10~15年基本保持不变。

5. 结论

基于实际工程概况通过COMSOL软件对路基长期使用性能进行数值模拟实验，分析路基在长期使用过程中路基温度，冻胀位移以及冻结深度的变化趋势，得到以下结论：

1) 路基以外地表以及路基边坡处温度侵入速度大于路基表面处；随着深度的增加路基温度变化有明显的滞后效应，路基3米深处达到高低温峰值的时间相比于路基表面分别滞后了约90天与120天。

2) 随着时间的增加，路基受冻融循环的影响逐渐变大，高低温峰值逐渐减小，随着深度的增加路基高低温峰值减小的幅度越大。深度为3 m处的路基最高温度10年间降低了2.68℃，最低温度降低了2.21℃。相比于0.5 m深度处路基高低温峰值降低幅值增加了1.79℃和0.76℃。

3) 随着时间增加路基冻胀位移在逐渐增大且路中位移大于左右路肩位置。路中与路肩冻胀位移经过10年时间分别增加了53.9%、47.8%，路中与路肩融沉位移10年间分别增加了1.01 mm和0.66 mm。路基冻胀效应明显强于融沉效应。路基冻结深度经过10年时间增加了80%。

4) 通过路基温度、冻胀融沉位移以及冻结深度变化趋势分析可知，THM耦合影响主要集中于前10年，10年后路基温度以及冻胀融沉位移变化基本趋于稳定。