短时冻区气候下土坡浅层水热状况的数值分析

阙 云，陈晓鹏，荣 耀

(1.福州大学 土木工程学院，福建 福州350108;2.江西农业大学 工学院，江西 南昌330045)

0 引 言

福建北部既有残积土边坡广泛在每年降雨量很少的冬季出现浅层垮塌现象，严重影响公路正常运营，但破坏原因始终存在困惑。现场调查表明，该区域内土质边坡冬季失稳原因与地处短时冻区所引发的边坡浅层短时反复冻融息息相关。因此，开展短时冻区气候下的土坡浅层水热状况研究具有重要意义。事实上，国内外众多学者针对冻土工程相关水热力变化展开了大量研究，取得了丰硕成果［1-9］，如:Harlan［1］(1973)提出了土冻结过程中水、热迁移相互作用的耦合模型并得到了广泛应用。热力学模型是Duquennoi 在1985 年首次提出［2］。何敏等［4］在饱和冻土多孔多相介质理论的基础上，建立了饱和冻土温度场、水分场、变形场耦合分析数值模型，并进行了验证。宋存牛［5］对国内外在水热力耦合作用机理、耦合作用方程方面的研究进行了综评，指出了各种研究结果之间差异、不同冻胀模式建立的依据、研究方法以及适用条件。然而以上研究多数集中于季节及多年冻土，针对短时冻土区边坡浅层水热变化研究还较少。相比与多年冻土与季节冻土，南方短时冻土无论是其深度、冻结及热融过程都不尽相同。基于此，本文依据南方短时冻区气候特点，针对短时冻区内残积土边坡浅层水热耦合进行数值分析，探究土坡浅层水热迁移状况。

1 水热耦合模型及程序编制

1.1 水热耦合模型

水热耦合模型采用Harlan 模型［6］:

式中，θ、θi分别为含水率和含冰率;D( θ) 、K(θ)分别为水分扩散率和导水系数;ρi、ρl分别为冰、水的密度;C、λ、Lf分别为比热容、导热系数及冻融潜热;T、t 分别为温度及时间。其中未冻水含量计算式通过测温法室内试验测得:θu=0.519T-1.18。

1.2 有限差分模型

基于交替隐式差分(ADI)的隐式差分方法对式(1)进行求解，采用MATLAB 语言编制程序，差分模型如图1 所示，其中坡高为7.5 m，坡底宽5 m。

1.3 边界条件

基于第二类边界条件确定坡面边界，主要考虑蒸发、太阳辐射、净辐射、对流换热、蒸发潜热。坡底坡高边界采用隔水隔温边界。其中蒸发强度采用Campbell 模型［10］，如式(2)。

式中，E 为蒸发强度;ρvs、ρva分别为地表、空气水汽密度;rv为阻力系数。

热边界主要由蒸发潜热、太阳总辐射、净辐射及对流换热累加而得。其中蒸发潜热采用式(3)［8］计算。

式中，L=2.49463×109-2.247×106Ts;Ts为土壤表层温度。

太阳辐射采用HIPERPAVⅡ模型［11］计算，如式(4)。

式中，qs为太阳辐射吸收热流，W/m2;αs为太阳辐射的吸收率;If为考虑太阳在24 h 中的辐射角度的辐射强度系数，一般假设其服从正弦分布，晚上为零;qsolar为太阳辐射峰值。

净辐射计算依据Stefan-Boltzman 定律确定［6］，其中当表层结冰时运用式(5)计算，表层不结冰时由式(6)计算。

式中，Lsp、Ls分别为雪层和土壤的表层净长波辐射量;TK为绝对温度;εac、σ 分别为热辐射系数(土壤取值为0.65，积雪取值为0.8)、常数(5.6697×10-8W/m2·K4)。

对流换热［9］采用式(7)计算，对流换热系数由式(8)确定。

式中，hc为对流换热系数;C 为常数，取决于形状和热流条件，当地表温度高于空气时取1.79，当地表温度低于空气时取0.89;w 为风速，m/s;Ta为空气温度;Ts为坡表温度。

选取2013 年12 月27 日福建省武夷山市省道303 线K364+900 处现场实测气象数据作为本文模拟基准气象边界条件，具体如表1 所示。

表1 气象边界条件参数Tab.1 Meteorological parameters of the boundary conditions

1.4 计算参数

土水特征曲线采用滤纸法确定，并运用(Van Genuchten 模型，以下简称VG 模型)拟合，VG 模型方程［12］如式(9)所示。

式中，θ 为土壤含水率;θr、θs分别为残余含水率以及饱和含水率(由试验测得分别为0.10、0.38);H 为土壤负压;α、n、m(m=1-1/n)为拟合参数，由试验结果拟合可得α=0.0163、n=1.7445。

未冻土导水系数K(θ)与冻土导水系数Ki(θ)分别由式(10)和式(11)［3］计算。

式中，Se为有效饱和度;Ks为饱和导水率，由室内试验可得土壤饱和导水率为0.002 m/h;I 为经验常数，可近似表示为I=1010θi。

未冻土水分扩散系数由式(12)计算，冻土水分扩散系数Di(θ)由式(13)［3］计算。

土壤热参数取值:未冻土骨料颗粒比热为0.84 kJ/(kg·K);冻土骨料颗粒比热为0.77 kJ/(kg·K);土壤中水和冰的比热分别取4.182 kJ/(kg·K)、2.09 kJ/(kg·K)［11］，冻土导热系数λf与未冻土导热系数λu分布按式(14)、(15)［13］计算得到。

式中，rd为土质干密度

1.5 初始条件

模拟时设定初始含水率为0.25，由室内冻结温度试验得，该含水率下土壤冻结温度约为-1.9 ℃。由现场监测结果［14］确定初始温度分布为:

式中，x 为坡面水平深度(m)，y 为坡面垂直高度(m)。

2 计算结果分析

2.1 基准条件下边坡浅层水热状况

依据表1 边界条件，通过数值模拟可得边坡坡中(坡高4 m)浅层温度及含冰率变化如图2 所示。由图2 可见:当气温由零上逐渐降低至零下时，边坡浅层由表至内逐层冻结，当日最低气温达至-4.3℃时，土坡浅层最大冻深约为5 cm。整个冻期下，土坡浅层结冰时间段为4:00 ～8:20，原因是由于7:00以后由于辐射影响浅层冰开始消融。边坡浅层含冰率峰值约为0.14，出现在7:00 时刻，位置并非出现在边坡表层，而是在浅层2 cm 左右位置。

图2 基准模拟计算结果Fig.2 Benchmark simulation results

2.2 日最低气温对浅层水热状况的影响

依据闽北地区历史气温资料，选取日最低气温分别为-7 ℃、-5 ℃、-4.3 ℃，其他参数与气象边界参见基准条件(见表1)，计算结果如图3 至图4 所示(其中日最低气温为-4.3 ℃的计算结果见图2)。

图3 日最低气温为－7 ℃时计算结果Fig.3 Calculation results when the daily minimum temperature is －7 ℃

图4 日最低气温为－5 ℃时计算结果Fig.4 Calculation results when the daily minimum temperature is －5 ℃

由图2、图3 和图4 可见:三种日最低气温条件下，边坡浅层最大冻深分别约为11 cm、7 cm 及5 cm。当日最低气温为-7 ℃时，其冻结时段为2:00 ～11:00，其余两种工况冻结时间较为接近，均在4:00 ～8:20。三种工况下，边坡浅层含冰率峰值分别为0.17、0.15、0.14，峰值出现时刻均在7:00 ～7:20，位置依次在浅层7 cm、4 cm、2 cm，由此可见最大含冰率并非出现在表层，而是随着冻深的增大而逐渐加深。

2.3 风速对边坡浅层水热状况的影响

依据闽北地区历史气温资料，选取平均风速分别为10 m/s、5.6 m/s、2 m/s，其他参数与气象边界参照基准条件(见表1)。计算结果如图5、图6 所示(其中平均风速为5.6 m/s 的计算结果见图2)。

图5 平均风速为10 m·s－1时计算结果Fig.5 Calculation results when the average wind velocity is 10 m·s－1

图6 平均风速为2 m·s－1时计算结果Fig.6 Calculation results when the average wind velocity is 2 m·s－1

由图2、图5 和图6 可见:三种风速条件下，边坡浅层最大冻深分别约为7 cm、5 cm 及5 cm。从冻结时间段中可以看出，平均风速为10 m/s 时，冻结时段为4:00 ～9:20，平均风速为2 m/s 时，冻结时段为5:00 ～7:40，说明较强的风速有利于延长土坡浅层结冰时间。当平均风速为10 m/s 时，含冰率峰值为0.14，当平均风速为2 m/s 时，含冰率峰值为0.12，两种工况下含冰率峰值出现位置均出现在2 cm 左右，说明风速对含冰率峰值及其出现位置影响不大。

2.4 辐射对边坡浅层水热状况的影响

依据闽北地区现场辐射监测结果，分别选取辐射峰值为200 W/m2、539 W/m2、800 W/m2，其他参数与气象边界参照基准条件(表1)。计算结果如图7、图8 所示(其中辐射峰值为539 W/m2的计算结果见图2)。

图7 辐射峰值为200 W·m－2计算结果Fig.7 Calculation results when the radiation peak is 200 W·m－2

图8 辐射峰值为800 W·m－2计算结果Fig.9 Calculation results when the radiation peak is 800 W·m－2

由图2、图7 和图8 可见:三种辐射条件下，土坡浅层冻深依次为6 cm、5 cm、5 cm，说明辐射对土坡浅层冻深影响不大。冻结时间上，三种工况下土坡浅层冻结时间段依次为:4:00 ～11:00、4:00 ～8:20、4:00 ～7:40，说明虽然辐射对起冻时间几乎无影响，但对土坡浅层冰保存时间有较大影响。三种辐射条件下，含冰率峰值均为0.14，且均出现在浅层2 cm 位置，说明辐射对于土坡浅层冻结主要影响其冻结时间，对含冰率峰值并没有太大影响。

3 结 论

①福建典型短时冻区气候下，土坡浅层最大的冻深约为10 cm 以内，浅层一般约在4:00 开始冻结，含冰率峰值一般出现在表层以下。短时冻区土质边坡冻融失稳分析时重点关注浅层10 cm 左右稳定性。

②日最低气温对土坡浅层冻深有较明显影响，风速及太阳辐射对土坡浅层冻深影响不大。

③日最低气温及风速对于土坡浅层冻结时间影响较大，二者的增大均会导致土坡浅层冻结时间的延长，此外，低辐射也能够促进土坡浅层冻结时间的延长。

④日最低气温对于含冰率峰值有较显著影响;风速及太阳辐射对于土坡浅层含冰率峰值影响不大。

［1］ HARLAN R L.Analysis of coupled heat-fluid transport in partially frozen soil［J］.Water Resoure Reserch，1973，9(5):1314-1323.

［2］ TAYLOR G S，LUTHIN J N.A model for coupled heat and moisture transfer during soil freezing［J］.Canadian Geotechnical Journal，1978，15(4):548-555.

［3］ 毛雪松，李宁，王秉纲，等.多年冻土路基水-热-力耦合理论模型及数值模拟［J］.长安大学学报:自然科学版，2006，26(4):16-20.

［4］ 何敏，李宁，刘乃飞.饱和冻土水热力耦合模型解析及验证［J］.岩土工程学报，2012，34(10):1058-1065.

［5］ 宋存牛.冻融过程中土体水热力耦合作用理论和模型研究进展［J］.冰川冻土，2010，32(5):982-987.

［6］ 郑秀清，樊贵盛，邢述彦.水分在季节性非饱和冻融土壤中的运动［M］.北京:地质出版社，2002.

［7］ 李双洋，张明义，黄志军，等.地铁双线隧道施工人工冻结水热力数值分析［J］.交通运输工程学报，2009，9(1):67-73.

［8］ 沈科伟，荣传新，余文文.查干淖尔矿副井冻结温度场有限元分析［J］.广西大学学报:自然科学版，2011，36(4):647-652.

［9］ 程学磊，李顺群，孙世娟，等.围压和负温对冻土强度和微结构的影响研究［J］.广西大学学报:自然科学版，2014，39(1):95-104.

［10］BRISTOW K L，HORTON R.Modeling the impact of partial surface mulch on soil heat and water flow［J］.Theoretical and Applied Climatology，1996，54(1-2):85-98.

［11］陈祖鑫.福建省花岗岩残积土路堤温湿耦合效应研究［D］.福州:福州大学土木工程学院，2012.

［12］VAN GENUCHTEN M T.A closed-form equation for predicting the hydraulic conductivity of unsaturated soils［J］.Soil Science Society of America Journal，1980，44(5):892-898.

［13］徐斅祖，王家澄，张立新.冻土物理学［M］.北京:科学出版社，2010.

［14］陈晓鹏.短时冻土区残积土边坡水热耦合机理的研究［D］.福州:福州大学土木工程学院，2014.