冻结条件下土壤水热耦合迁移的数值模拟

步骤 4 由式( 14) 校正未冻 水含水率, 若其
值大于式( 14) , 则超出部分为本时段增加的含冰
率, 即
= k+ 1
Ii
- k+ 1
i
m+
k Ii
。
步骤 5 重复步骤 2~ 4, 直至前后两次计算
的含水率与温度之差小于 0. 01 为止。
4 结果与分析
应用上述参数及数学模型仿真分析土壤冻结 过程中的未冻水含水率及温度, 结果如图 1、2 所示。 由图可看出: 当土壤深度< 30 cm 时, 土壤温度 变化剧烈, 在周期性的环境温度和太阳辐射的作 用下, 土壤中各点 的温度亦呈周 期性变 化。50、 70、110 cm 处的土壤温度变化不 明显, 振幅趋于 消失。 在冻结阶段, 由于表层土壤未冻水含水 率较小, 所以在 0~ 10 cm 处未冻水含水率变化不 大; 10~ 50 cm 未冻水含水率变小, 主要原因为温
度为 0~ 20 cm 时:
5 108 e- 41. 379
> 0. 336
s( ) = 166 374 2- 17 403 + 4 554
( 4) 0. 336
当土层深度为 20~ 150 cm 时: 2 107 e- 35.119
s( )= 1 972 1 2- 5 749. 3 + 2 102. 9
1 g / cm3 ; L I 为冰的溶解潜热, 一般取 335 J/ cm3 ;
m 为相应土体 负温下 可能 的最大 未冻水 含量,
cm3 / cm 3 。
2. 3 边界条件 冻结土壤的水热运动状况是由土壤的上下边
界的水热条件决定的[ 4] 。热流方程上下边界条件
为第一类边界条件, 上下边界温度取实测值:
unf ro ze n w a t e r c o nt ent
土壤深度/ cm 0~ 20 20~ 60 60~ 80
a 0. 170 48 0. 214 63 0. 234 89
b 0. 162 4 0. 184 8 0. 109 2
R2 0. 739 5 0. 520 2 0. 829 3
注: R2 为相关系数。
当 + I s 时:
K1 = X / s +
1 / w+
I/ I
K2 = X s + I I + w 式中, K 1 、K 2 分别为串、并联土层的热导率; 为 串、并联土层的体积比, 一般取 0. 5; s 为饱和含 水率; X 、、I、A 分别为土颗粒、未冻水、冰、空气 所占的体积; s 、w 、I 、a 分别为土颗粒、水、冰、空 气的热传导系数, W/ ( m K) , 数值分别为 3. 244 5、 0. 562 0、2. 180 0、0. 025 0。
> 0. 304 ( 5)
0. 304
收稿日期: 2010- 01- 24, 修回日期: 2010- 03- 08 作者简介: 刘畅( 1984-) , 男, 硕士研究生, 研究方向为水土环境与生态工程, E- mail: liuchang 629@ 163. com 通讯作者: 陈晓飞( 1964- ) , 女, 教授、博导, 研究方向为水土环境与生态工程, E- mail: chenx iao20302@ vip. 163. co m
s= -
( 19)
3 计算方法与步骤
由于冻融条件下, 水热耦合运移模型的非线
性及土壤的非均质和初始边界条件的复杂性, 目
前用解析方法求解很不理想, 采用数值计算方法
最有效。有限差分法以差商近似替代微商, 将微
分方程变成差分方程, 组成可直接求解的代数方
程组。选用的中心差分格式精度较高, 但时间和空
间步长较大时, 会出现严重的数值振荡[ 3] , 因此应结
( 7)
K I( ) = K( )/ I
( 8)
式中, 1 为含水率, cm3 / cm3 。
1. 4 冻结区未冻水含水率
冻土中含有未冻水是由土壤自身性质决定。
未冻水在冻土中可自由移动。影响其含量的因素 很多, 其中影响最大、起决定性作用的是土壤的负 温。未冻水含量随土壤负温的增大而减小, 每种 土壤均有较固定的未 冻水含量与负 温的关系曲
图 1 土壤不同深度 处温度随时间的变化 Fig . 1 Va riat io n of so il t e mpe ra t ure w it h t ime
a t dif f er ent de pt hs
合实例选取时间与空间步长。本文空间步长 z 、时
间步长 t 分别取 2. 5 cm、10 min。计算步骤如下。
步骤 1 定义初始含水率、温度, 用初始值先
代替下一时刻的含水率、温度, 如
= k+ 1
i
k i
。
步骤 2 求解水分方程组, 获得未冻水含水
率预测值。
步骤 3 求解热流方程组, 获得温度预测值。
线。表 1 为根据现场试验结果拟合的未冻水含量
与负温关系, 拟合方程为[ 1] :
= aT - b
( 9)
式中, 为未冻水含水率, cm3/ cm3 ; T 为土层负温
绝对值, ; a、b 为经验常数。
表 1 温度与未冻水 含量的关系
Tab. 1 Re lat io nship o f t e mpe ra t ure and
h, q> 0 表明通量向下; q < 0 表明通量向上; z 1 、z 2 为土壤深度, cm; 1 、2 为 z 1 、z 2 处基质势; z 为
空间步长; K ( ) 为 z 处非饱和土壤导水率与基质
势的关系。
基质势 ( 单位: cm) 与土壤水吸力 s( 单位: cm) 的关系为[ 2] :
பைடு நூலகம்
- D( ) z + K ( ) = q( t) t > 0, z = 150 ( 17) 由于地下水埋深约 4 m, 下边界即 1. 5 m 处
的土壤水分处于不饱和状态, 其水分通量为[ 2] :
96
水电能源科学
2010 年
q( z ) = - K ( )
2- 1 + 1 z
( 18)
其中
z = ( z1+ z2)/2 式中, q( z ) 为边界处单位时间土壤水分通量, cm/
1 数据采集与参数测定
1. 1 试验区设计 2007 年 11 月 26 日~ 12 月 8 日在沈阳农业
大学水利学院综合试验基地采集野外数据。测区 为一个由不透水保温板围成的 3 m 3 m 26 m ( 长 宽 高) 试验区。测区埋设 T RM- Z S1 气象
及生态环境监测系统探头, 可收集气象数据及地 表至地下 1. 5 m 处土壤各 剖面的温度数据; 1. 5 m 深的冻土器可测量冻深; 在固定位置安放中子 仪测管和 T DR 测管测量总含水率和未冻水含水 率; 测区土壤为潮棕壤土, 0~ 20 cm 的土样干容 重为 1. 39 g/ cm3 , 饱和含水率 0. 48 cm3 / cm3 ; 20~
1. 5 热导率与比热容 估算土壤热导率 K h 十分困难, K h 不仅取决
于土壤中各组成物的比例, 还与各组成物的形状 等有关。本文采用的计算公式为[ 6] :
Kh =
1 /K1 + (1- )/K 2
( 10)
当 + I < s 时:
K1 = X/ s+
1 / w + I/ I + A/ a
K 2 = X s + I I+ w + A/ a
第 28 卷第 5 期
刘 畅等: 冻结条件下土壤水热耦合迁移的数值模拟
95
1. 3 冻结区水分运动参数
冻结区的扩散率和导水率很难测定, 本文采 用阻抗系数法。分别用未冻结区的扩散率和导水
率除以阻抗系数 I 表示, I 值取决于土壤中的含
冰率[ 2] :
I = 1010 I
( 6)
DI( ) = D( )/ I
150 cm 的土样干容重为 1. 47 g / cm3, 饱和含水率
0. 45 cm3 / cm3 。
1. 2 非冻结区水分运动参数 土壤非饱和扩散率 D ( ) ( 单位: cm2/ min) 通
过室内水平土柱吸渗法确定。
当土层深度为 0~ 20 cm 时:
D ( ) = 0. 098 4e7. 482 D ( ) = 0. 000 5e22. 144
土壤体积比热容具有可加性, 是土壤中各组 成物热容量之和[ 1] :
cv = cvs x s + c vw x w + cvlx vl + c va x va ( 11) 式中, x s、x w 、x vI 、x va 分别为单 位体积土壤 矿物、 水、冰、空气所占体积; c vs 、cvw 、cvI 、cva 分别为相应 组成物的体积比热容, J/ ( cm3 K ) , 数值分别为 2. 223、4. 187、1. 730、0. 001 2。
0. 38 ( 1)
> 0. 38
当土层深度为 20~ 150 cm 时:
D ( ) = 0. 000 6e18. 422
( 2)
土壤非饱和导水率 K ( ) ( 单位: cm/ m in) 与
D( ) 土壤水吸力 s( ) 的关系为:
K(
)
=-
D(
)
d ds
( 3)
s( ) ( 单位: cm) 由压力仪法测定。当土层深
( 沈阳农业大学 水 利学院, 辽宁 沈阳 110161)
摘要: 以沈阳市季节性冻土为例, 基于冻结条件下水热耦合 迁移的数 学模型, 采用 中心差 分格式并 用线性 化
迭代法求解, 对现场土壤冻结过程进行了模拟。与实测 结果对比表 明, 采 用的参数及 数学模型 可成功模拟 现
场土壤冻结过程, 并讨论了土壤冻结时水分与温 度的迁移规律, 可为后续研究奠定基础。
土壤 中 未 冻水 含 水 率 和 温 度 T 的 关 系 ( 式
( 9) ) , 水分运动方程、热流方程、联系方程分别为:
t=
z D( )
z-
K( z
)
-
I w
I
t
( 12)
cv
T t
=
z Kh
T z
+
LI
I
I
t
( 13)
m(T )
( 14)
式中, I 为冰的密度, 0. 917 g/ cm3 ; w 为水的密度,
2 数学模型
2. 1 基本假设 假设土壤骨架无膨胀和收缩变形、均质且
各向同性, 水热迁移主要发生在垂向, 可近似为一
维问题。 土壤内部满足局部热力学平衡条件。
冻土中水分迁移以液态形式运动, 认为冰固定
不动, 忽略气相的影响。
2. 2 基础方程 冻结条件下, 水热耦合运移问题必须联立求
解水分运动和热流两个基本方程。需求解 、I、 三个未知函数, 因此, 还须补充一个联系方程, 即
关键词: 冻土; 冻结条件; 冻结过程; 水热耦合迁移; 数值模拟
中图分类号: T P391. 9; S152. 9
文献标志码: A
季节性冻土是一种含冰晶的特殊土水体系。 土壤冻结过程非常复杂, 伴随物理、化学、力学现 象和过程, 包括水分、热量迁移、水分的相变和盐 分的积聚[ 1] 。国内对冻土水热耦合迁移的数值模 拟问题研究较晚, 起步于 20 世纪 80 年代末期, 但 发展较快。雷志栋等[ 2] 采用 H arlan 模型模拟了 水平、垂直土柱的冻结过程, 并定性分析了土壤初 始含水率对土壤冻胀量的影响; 尚松浩等[ 3] 模拟 了室内土柱冻结试验和田间的水热状况, 改进了 求解方法使计算速度加快; 郑秀清[ 4] 对山西省汾 河灌区土壤冻融过程水热动态进行仿真分析; 李 扬[ 5] 基于非饱和土多孔介质理论的水热迁移耦合 模型, 采用有限元法求解, 模拟了长春市季节性冻 土水分迁移, 经试验验证, 效果较好。鉴此, 本文 采用 M at lab7. 1 软件编制计算程序, 应用中心差 分格式并用线性化迭代法求解冻土一维垂直水热 耦合迁移数学模型, 对沈阳农业大学水利学院综 合试验基地土壤在冻结过程中的水、热动态进行 了模拟, 为揭示季节性冻土水热耦合迁移状况提 供依据。
Tt T( z) =
Tb
t > 0, z = 0 ( 15)
t > 0, z = 150
水流差分 方程上边界条 件为第一类 边界条
件, 地表含水率较低近似于风干含水率 c = 0. 05 cm3 / cm3; 下边界有水分通量, 为第二类边界条件,
方程如下:
( z) = c
t > 0, z = 0 ( 16)
第 28卷 第 5期 2 0 10 年5 月
文章编号: 1000- 7709( 2010) 05- 0094- 04
水电能源科学 W ater Resour ces and P ow er
V o l. 28 N o . 5 M a y. 2 0 1 0
冻结条件下土壤水热耦合迁移的数值模拟
刘 畅 陈晓飞2 苑 杰3 李晓燕3