冻土水热运动的数值模拟
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
附表 实验初始温度分布
பைடு நூலகம்
深度 Z ( cm)
0
2
6
12
20
25
温度 T ( °C)
2. 3
3. 9
6. 1
5. 9
6. 0
5. 5
2 数学模型
假定土壤介质不可压缩、均质且各向同性; 冻土中只有液态水的运动, 认为冰固定不 动, 忽略汽相影响; 冻土中未冻水含量与土壤负温处于动态平衡中; 水、热的运动主要发生 在垂向, 可近似为一维问题。 2. 1 水流方程
3 表面半控制容积的水量平衡图。1——( 1+
1/ 2) 半控制容积的水量平衡方程为:
(R-
q1+
1 2
)
△
t=
(
i△ il /
W+
△
ul
)
△Z 2
1
(
7)
其中表面通量 R= 0
q1+
1 2
=
-
D(
u)
u
Z
+
K
(
u)
将式( 7) 简化后得:
图 1 半控制容积水量平衡图
E 1 2uo+ F 1 u 12+ G1 2u2= H 1
第 19 卷 第 1 期 V ol. 19 No . 1
内 蒙 古 农 牧 学 院 学 报 Journ al of Inn er M ongolia In st it ut e of A gricult ure & A nimal Husb andry
冻土水热运动的数值模拟
王海丽
( 广东省水利科学研究所 广州)
i
t
( 2)
其中 T 为土壤温度, L 为融化潜热, CV 为土壤热容量, 为土壤热传导率。
2. 3 联系方程
由( 1) 、( 2) 组成的方程组, 有 u、T 及 i 三个未知数, 却只有两个方程, 求解的关键需
补充一联系方程, 才可使方程组解。冻土中, 未冻水含量与负温保持动态平衡, 这一关系表
Key words soil f reezing , coupled moisture and heat t ransfer, num erical sim ulatio n, implicit f init e dif ference
引 言
全球陆地面积的 70% 存在冻土, 其中 14% 为永久冻土, 56% 为季节性冻土〔6〕; 在中
件的, 但考虑到其精确性, 起初作者采用了中心差分格式计算, 结果出现强烈振荡, 如温度
可以从 60°C 振荡到- 17°C, 而含水量有时会出现负值。改用全隐式差分法后, 振荡现象大
为减轻, 且收敛较好, 因此最后采用全隐式差分格式。
对水流方程离散得全隐式差分方程:
- k+ 1
ui
△t
k ui
对冻土中水热迁移规律的研究已成为国际前沿学科关注热点, 但对工程冻土的研究 较多, 对农业冻土的研究较少。
Nakano 等( 1982～1984) 在试验基础上研究了等温条件下的水分迁移; M illy 采用了 以土壤基质势和温度为变量的土壤水、汽、热耦合方程, 适用于非均质土壤, 并用有限单元 法模拟了等温、非等温条件下的土壤水分运动〔10〕; Jame 和 N orum ( 1980) 利用有限差分法 模拟 了水 平土柱 冻结情 况下温 度、含 水量、含冰量 的动 态变 化; F ukuda 和 Nakag aw a ( 1985) 对自然条件下土壤的冻结过程进行了模拟。
国内此方面的研究起步较晚, 但发展很快。杨诗秀、雷志栋等模拟了水平、垂直土柱冻 结过程〔4〕; 叶佰生, 陈肖柏( 1990) , 胡和平等引入了 Clapeyron 方程来研究冻土中水热迁 移问题; 尚松浩( 1997) 根据冻土水热基本方程导出了水热耦合方程, 对不同条件下整个越 冬期长时间的土壤冻融过程中的水热迁移进行了模拟。
1 2
-
K
k-
1 2
)
-
i(
- k+ 1
ji
W
kj i)
P i=
Zi+
△t 1- Zi-
)
1
同理; 将热传导方程列全隐式差分方程并整理得:
E
iT
k+ i-
11+
F
iT
k+ i
1+
G
iT
k+ i+
1 1
=
H
i
( 6)
其中,
E i=
-
2P
i
i-
1 2
Z i -Z i -
)
1
Gi =
1 998
整理后得:
E i
+ k+ 1
ui- 1
Fi
+ k+ 1
ui
Gi
= k+ 1
ui+ 1
Hi
( 5)
其中,
E i=
-
2P
i
D i-
1 2
Zi- Zi-
)
1
Gi =
-
2Pi
D i+ Zi+ 1 -
1 2
Z
i
)
F i= 1- Ei - Gi
H i =
k ui
-
2P i( K k+
Abstract Based o n t he m at hemat ical model for solv ing t he problem of coupled moist ure and heat t ransf er during soil f reezing, t he numerical calculation w as carr ied o ut in t his paper. An impl icit f init e dif ference scheme is used t o solve t he mo del. An ex periment w as conduct ed t o v er if y t he model . T he calculat ed result s agreed w ell w ith t he o bserved ones, which show ed that t he model w as reliable. U sing t his mo del, Sim ulat ions on o ther condit ions ( freezing dur ing 30d, et c. ) w er e conducted.
示了冻土中水、热运动间的相互关系:
u ≤ m( T )
( 3)
其中 m( T ) 为相应壤负温下可能的最大未冻水含量( 土壤冻结特性曲线) 。
方程( 1) ～( 3) 构成了冻结条件下土壤水热迁移的数学模型。
3 数值计算方法
隐式差分格式和中心差分格式的差分方程及求解方程组形式上完全一样, 差别仅在
在以上假定情况下, 冻土中水分运动规律与非饱和土壤水分运动规律类似, 土壤水分 运动可用含有相变项的 Rechar ds 方程表示:
第 1 期 王海丽: 冻土水热运动的数值模拟
101
u
t
=
Z[ D(
u)
u
Z
-
K
(
u) ] -
i W
i t
(
1
于方程系数与常数项的计算取值稍有不同。求解抛物型方程 u = D
2
u
2
(
其中
D
为常数)
t
x
时, 在数学上虽可以证明隐式差分格式和中心差分格式是无条件收敛和稳定的, 而对非线
性微分方程数值计算的收敛性和稳定性, 在数学上很难进行严格分析, 目前并无一套成熟
的理论来选取时间、空间步长。一般认为中心差分精度较高, 但其稳定是有条件的, 对步长
1998 年 3 月 M ar . 1998
摘要 本文在建立冻土水热耦合迁移的数学模型基础上, 采用全隐式有限差 分格式并用线性
化迭代方法求解。边界条件处理采用半控制 容积法, 时间步长根据 迭代情况自动调节, 求解过 程中为防止振荡过大采用松驰迭代 等措施。通过计算与实验结果对比, 取得较好效果, 表明采 用的数学模型及数值计算方法可靠。用此模型模拟了非实验状态下( 历时 30d 及表面温度降至 - 20°C) 的土壤冻结过程。模拟显示温度分布随冻结历时的加长而变缓, 最后趋于线性分布, 随 上、下表面温差的加大而变陡; 水分分布基本 形状并不随时间及上、下表面温差而变, 只是冻深 及上迁水量随历时加长、温差加大而增大。
收稿日期: 1997- 07- 31 本项研究系日本文部省资助的中日联合科研项目 作者: 女, 24 岁, 工学硕士, 从事农业水土工程的研究
10 0
内 蒙 古 农 牧 学 院 学 报
1 998
国, 多年冻土区的分布面积约为 2. 086×106km2, 季节性冻土区的 分布面积为 5. 137× 106km2, 两者合计占全国面积的 75% ( 徐学祖等, 1983) 。季节性冻土的研究是内蒙古自治 区盐渍化防治及工程防冻研究的目标。
4 土壤特性参数
4. 1 土壤冻结特性曲线
由于土壤持水特性的滞后效应, 不同初始含水量及冻融过程中, 土壤未冻水含量与负
温的关系曲线 m( T ) 并非单值。但对于同一种土壤在负温一定时, 不同条件下未冻水含量
相差不 大, 为了 简化 可将 m ( T ) 视 为单 值曲 线〔5〕。图 2 给出沙壕渠土的冻结特性曲线。
选取要求很严格, 隐式差分精度虽然较中心差分差些, 但是无条件稳定并收敛, 故在满足
精度要求前提下, 隐式差分在数值计算中常被采用。史、陈也曾考虑用中心差分进行计算
饱和 非饱和流水分传输, 在比较小时稳定性及收敛性较好, △t 较大时出现数值跳动 或不收敛 。 〔1〕、〔2〕
结合到本文实例, 方程( 1) 、( 2) 均属非线性微分方程, 虽中心差分格式的稳定是有条
( 8)
其中,
E 1= 0. 0
G1 =
-
D 1+
1 2
△t △Z1
2
F 1=
1 2
-
G1
H 1 =
1 2
1 1-
1 2
i(
W
il 2-
il 1) -
K 1+
1 △t 2 △Z1
第 1 期 王海丽: 冻土水热运动的数值模拟
103
对下边界采取同样的方法处理。 时间步长采取自动控制, 即计算中收敛情况较好时, 选取较大△t, 较差时选取较小 △t , 既提高了计算精度, 又提高了计算效率。空间步长开始时曾考虑用变化的△Z, 在冻结 锋面处取△Z 最小, 然后按等比级数向两端逐渐加大, 但由于冻结锋面不断下移, 节点坐 标变亦随之改变, 涉及到新节点处水量及温度的取值等一系列复杂的问题, 处理不好反而 适得其反, 最终还是采用了恒定△Z。迭代中为防振荡过大, 还采取了松驰迭代等措施。
)
其中 t、Z 分别为时间、空间坐标( 以垂直向下为正) ; u、i 分别为土壤未冻水、冰的体积含
量; i、W 分别为冰、水的密度; D ( u ) 、K ( u) 分别为非饱和冻土水分扩散率、导水率。
2. 2 热传导方程
以相变潜热作为内热源的热传导方程可表示为:
CV
T t
=
Z[
TZ] + L i
关键词 土壤冻结 水热耦合迁移 数值模拟 隐 式差分 中图资料分类号 S152. 7
NUMERICAL SIMULATION OF MOISTURE AND HEAT TRANSFER DURING SOIL FREEZING
Wang Haili
( Depar tment o f Wat er Conser vancy Engineer ing , Inner M ongo lia Instit ute o f A gr icult ur e and A nimal Husbandr y )
=
( Di+
1 2
- k+ 1
ui+ 1
Zi+ 1 -
k+ 1
ui
Zi
-
D i-
1 2
k+ ui
Z
1i- Z
k+ 1
) ui- 1
i- 1
1 2
(
Zi+
1-
Zi- 1 )
-
K i+
1 2
-
K i-
1 2
1 2
(
Zi+
1-
Zi- 1 )
-
i(
- k+ 1
i
W △t
k i
)
( 4)
10 2
内 蒙 古 农 牧 学 院 学 报
综上所述, 国内外对冻土水热的研究已有了很大进展, 但有的模拟缺乏充分的实验背 景( 如参数等) , 有的则未进行充分展开。本文在室内实验基础上进行数学模型检验, 对非 实验状态进行了模拟预测。
1 实验
实验在日本冈山大学农学部农地整备实验室进行, 试料采用内蒙古河套灌区沙壕渠 试验场的土壤, 初始含水量均一为 12. 87cm 3/ cm3, 初始温度分布如附表, 上表面温度按一 定规律降至- 6. 0°C, 下表面温度降至 4. 5°C, 土柱上、下两端密封。实验开始时, 在上表面 瞬时加一负温, 冻结历时 2880min。
-
2Pi
Zi+
i+
1 2
1-
Zi
)
F i= CVi - Ei - Gi
H i = CViT ki + L i
i(
- k+ 1
ji
k ji
)
对于边界条件处理, 由实验所述可知温度上、下边界条件均为已知, 属第一类边界, 水
分上、下边界为通量为零的第二类边界, 计算
时采用半控制容积法。如对上表面, 图 1 列出
பைடு நூலகம்
深度 Z ( cm)
0
2
6
12
20
25
温度 T ( °C)
2. 3
3. 9
6. 1
5. 9
6. 0
5. 5
2 数学模型
假定土壤介质不可压缩、均质且各向同性; 冻土中只有液态水的运动, 认为冰固定不 动, 忽略汽相影响; 冻土中未冻水含量与土壤负温处于动态平衡中; 水、热的运动主要发生 在垂向, 可近似为一维问题。 2. 1 水流方程
3 表面半控制容积的水量平衡图。1——( 1+
1/ 2) 半控制容积的水量平衡方程为:
(R-
q1+
1 2
)
△
t=
(
i△ il /
W+
△
ul
)
△Z 2
1
(
7)
其中表面通量 R= 0
q1+
1 2
=
-
D(
u)
u
Z
+
K
(
u)
将式( 7) 简化后得:
图 1 半控制容积水量平衡图
E 1 2uo+ F 1 u 12+ G1 2u2= H 1
第 19 卷 第 1 期 V ol. 19 No . 1
内 蒙 古 农 牧 学 院 学 报 Journ al of Inn er M ongolia In st it ut e of A gricult ure & A nimal Husb andry
冻土水热运动的数值模拟
王海丽
( 广东省水利科学研究所 广州)
i
t
( 2)
其中 T 为土壤温度, L 为融化潜热, CV 为土壤热容量, 为土壤热传导率。
2. 3 联系方程
由( 1) 、( 2) 组成的方程组, 有 u、T 及 i 三个未知数, 却只有两个方程, 求解的关键需
补充一联系方程, 才可使方程组解。冻土中, 未冻水含量与负温保持动态平衡, 这一关系表
Key words soil f reezing , coupled moisture and heat t ransfer, num erical sim ulatio n, implicit f init e dif ference
引 言
全球陆地面积的 70% 存在冻土, 其中 14% 为永久冻土, 56% 为季节性冻土〔6〕; 在中
件的, 但考虑到其精确性, 起初作者采用了中心差分格式计算, 结果出现强烈振荡, 如温度
可以从 60°C 振荡到- 17°C, 而含水量有时会出现负值。改用全隐式差分法后, 振荡现象大
为减轻, 且收敛较好, 因此最后采用全隐式差分格式。
对水流方程离散得全隐式差分方程:
- k+ 1
ui
△t
k ui
对冻土中水热迁移规律的研究已成为国际前沿学科关注热点, 但对工程冻土的研究 较多, 对农业冻土的研究较少。
Nakano 等( 1982～1984) 在试验基础上研究了等温条件下的水分迁移; M illy 采用了 以土壤基质势和温度为变量的土壤水、汽、热耦合方程, 适用于非均质土壤, 并用有限单元 法模拟了等温、非等温条件下的土壤水分运动〔10〕; Jame 和 N orum ( 1980) 利用有限差分法 模拟 了水 平土柱 冻结情 况下温 度、含 水量、含冰量 的动 态变 化; F ukuda 和 Nakag aw a ( 1985) 对自然条件下土壤的冻结过程进行了模拟。
国内此方面的研究起步较晚, 但发展很快。杨诗秀、雷志栋等模拟了水平、垂直土柱冻 结过程〔4〕; 叶佰生, 陈肖柏( 1990) , 胡和平等引入了 Clapeyron 方程来研究冻土中水热迁 移问题; 尚松浩( 1997) 根据冻土水热基本方程导出了水热耦合方程, 对不同条件下整个越 冬期长时间的土壤冻融过程中的水热迁移进行了模拟。
1 2
-
K
k-
1 2
)
-
i(
- k+ 1
ji
W
kj i)
P i=
Zi+
△t 1- Zi-
)
1
同理; 将热传导方程列全隐式差分方程并整理得:
E
iT
k+ i-
11+
F
iT
k+ i
1+
G
iT
k+ i+
1 1
=
H
i
( 6)
其中,
E i=
-
2P
i
i-
1 2
Z i -Z i -
)
1
Gi =
1 998
整理后得:
E i
+ k+ 1
ui- 1
Fi
+ k+ 1
ui
Gi
= k+ 1
ui+ 1
Hi
( 5)
其中,
E i=
-
2P
i
D i-
1 2
Zi- Zi-
)
1
Gi =
-
2Pi
D i+ Zi+ 1 -
1 2
Z
i
)
F i= 1- Ei - Gi
H i =
k ui
-
2P i( K k+
Abstract Based o n t he m at hemat ical model for solv ing t he problem of coupled moist ure and heat t ransf er during soil f reezing, t he numerical calculation w as carr ied o ut in t his paper. An impl icit f init e dif ference scheme is used t o solve t he mo del. An ex periment w as conduct ed t o v er if y t he model . T he calculat ed result s agreed w ell w ith t he o bserved ones, which show ed that t he model w as reliable. U sing t his mo del, Sim ulat ions on o ther condit ions ( freezing dur ing 30d, et c. ) w er e conducted.
示了冻土中水、热运动间的相互关系:
u ≤ m( T )
( 3)
其中 m( T ) 为相应壤负温下可能的最大未冻水含量( 土壤冻结特性曲线) 。
方程( 1) ～( 3) 构成了冻结条件下土壤水热迁移的数学模型。
3 数值计算方法
隐式差分格式和中心差分格式的差分方程及求解方程组形式上完全一样, 差别仅在
在以上假定情况下, 冻土中水分运动规律与非饱和土壤水分运动规律类似, 土壤水分 运动可用含有相变项的 Rechar ds 方程表示:
第 1 期 王海丽: 冻土水热运动的数值模拟
101
u
t
=
Z[ D(
u)
u
Z
-
K
(
u) ] -
i W
i t
(
1
于方程系数与常数项的计算取值稍有不同。求解抛物型方程 u = D
2
u
2
(
其中
D
为常数)
t
x
时, 在数学上虽可以证明隐式差分格式和中心差分格式是无条件收敛和稳定的, 而对非线
性微分方程数值计算的收敛性和稳定性, 在数学上很难进行严格分析, 目前并无一套成熟
的理论来选取时间、空间步长。一般认为中心差分精度较高, 但其稳定是有条件的, 对步长
1998 年 3 月 M ar . 1998
摘要 本文在建立冻土水热耦合迁移的数学模型基础上, 采用全隐式有限差 分格式并用线性
化迭代方法求解。边界条件处理采用半控制 容积法, 时间步长根据 迭代情况自动调节, 求解过 程中为防止振荡过大采用松驰迭代 等措施。通过计算与实验结果对比, 取得较好效果, 表明采 用的数学模型及数值计算方法可靠。用此模型模拟了非实验状态下( 历时 30d 及表面温度降至 - 20°C) 的土壤冻结过程。模拟显示温度分布随冻结历时的加长而变缓, 最后趋于线性分布, 随 上、下表面温差的加大而变陡; 水分分布基本 形状并不随时间及上、下表面温差而变, 只是冻深 及上迁水量随历时加长、温差加大而增大。
收稿日期: 1997- 07- 31 本项研究系日本文部省资助的中日联合科研项目 作者: 女, 24 岁, 工学硕士, 从事农业水土工程的研究
10 0
内 蒙 古 农 牧 学 院 学 报
1 998
国, 多年冻土区的分布面积约为 2. 086×106km2, 季节性冻土区的 分布面积为 5. 137× 106km2, 两者合计占全国面积的 75% ( 徐学祖等, 1983) 。季节性冻土的研究是内蒙古自治 区盐渍化防治及工程防冻研究的目标。
4 土壤特性参数
4. 1 土壤冻结特性曲线
由于土壤持水特性的滞后效应, 不同初始含水量及冻融过程中, 土壤未冻水含量与负
温的关系曲线 m( T ) 并非单值。但对于同一种土壤在负温一定时, 不同条件下未冻水含量
相差不 大, 为了 简化 可将 m ( T ) 视 为单 值曲 线〔5〕。图 2 给出沙壕渠土的冻结特性曲线。
选取要求很严格, 隐式差分精度虽然较中心差分差些, 但是无条件稳定并收敛, 故在满足
精度要求前提下, 隐式差分在数值计算中常被采用。史、陈也曾考虑用中心差分进行计算
饱和 非饱和流水分传输, 在比较小时稳定性及收敛性较好, △t 较大时出现数值跳动 或不收敛 。 〔1〕、〔2〕
结合到本文实例, 方程( 1) 、( 2) 均属非线性微分方程, 虽中心差分格式的稳定是有条
( 8)
其中,
E 1= 0. 0
G1 =
-
D 1+
1 2
△t △Z1
2
F 1=
1 2
-
G1
H 1 =
1 2
1 1-
1 2
i(
W
il 2-
il 1) -
K 1+
1 △t 2 △Z1
第 1 期 王海丽: 冻土水热运动的数值模拟
103
对下边界采取同样的方法处理。 时间步长采取自动控制, 即计算中收敛情况较好时, 选取较大△t, 较差时选取较小 △t , 既提高了计算精度, 又提高了计算效率。空间步长开始时曾考虑用变化的△Z, 在冻结 锋面处取△Z 最小, 然后按等比级数向两端逐渐加大, 但由于冻结锋面不断下移, 节点坐 标变亦随之改变, 涉及到新节点处水量及温度的取值等一系列复杂的问题, 处理不好反而 适得其反, 最终还是采用了恒定△Z。迭代中为防振荡过大, 还采取了松驰迭代等措施。
)
其中 t、Z 分别为时间、空间坐标( 以垂直向下为正) ; u、i 分别为土壤未冻水、冰的体积含
量; i、W 分别为冰、水的密度; D ( u ) 、K ( u) 分别为非饱和冻土水分扩散率、导水率。
2. 2 热传导方程
以相变潜热作为内热源的热传导方程可表示为:
CV
T t
=
Z[
TZ] + L i
关键词 土壤冻结 水热耦合迁移 数值模拟 隐 式差分 中图资料分类号 S152. 7
NUMERICAL SIMULATION OF MOISTURE AND HEAT TRANSFER DURING SOIL FREEZING
Wang Haili
( Depar tment o f Wat er Conser vancy Engineer ing , Inner M ongo lia Instit ute o f A gr icult ur e and A nimal Husbandr y )
=
( Di+
1 2
- k+ 1
ui+ 1
Zi+ 1 -
k+ 1
ui
Zi
-
D i-
1 2
k+ ui
Z
1i- Z
k+ 1
) ui- 1
i- 1
1 2
(
Zi+
1-
Zi- 1 )
-
K i+
1 2
-
K i-
1 2
1 2
(
Zi+
1-
Zi- 1 )
-
i(
- k+ 1
i
W △t
k i
)
( 4)
10 2
内 蒙 古 农 牧 学 院 学 报
综上所述, 国内外对冻土水热的研究已有了很大进展, 但有的模拟缺乏充分的实验背 景( 如参数等) , 有的则未进行充分展开。本文在室内实验基础上进行数学模型检验, 对非 实验状态进行了模拟预测。
1 实验
实验在日本冈山大学农学部农地整备实验室进行, 试料采用内蒙古河套灌区沙壕渠 试验场的土壤, 初始含水量均一为 12. 87cm 3/ cm3, 初始温度分布如附表, 上表面温度按一 定规律降至- 6. 0°C, 下表面温度降至 4. 5°C, 土柱上、下两端密封。实验开始时, 在上表面 瞬时加一负温, 冻结历时 2880min。
-
2Pi
Zi+
i+
1 2
1-
Zi
)
F i= CVi - Ei - Gi
H i = CViT ki + L i
i(
- k+ 1
ji
k ji
)
对于边界条件处理, 由实验所述可知温度上、下边界条件均为已知, 属第一类边界, 水
分上、下边界为通量为零的第二类边界, 计算
时采用半控制容积法。如对上表面, 图 1 列出