三维双重孔隙-裂隙介质热-水-应力-迁移耦合模型及其有限元分析知识讲解.docx

此文档收集于网络，如有侵权请联系网站删除
三维双重孔隙-裂隙介质热-水-应力-迁移耦合模型
及其有限元分析

张玉军①*, 张维庆②
① 中国科学院武汉岩土力学研究所, 岩土力学与工程国家重点实验室, 武汉 430071;
② 西南交通大学土木工程学院, 成都 610031
* E-mail: yjzhang@
收稿日期: 2009-10-29; 接受日期: 2010-04-29
国家重点基础研究发展计划(“973”计划)(批准号: 2010CB732101)和岩土力学与工程国家重点实验室前沿探索性项目(批准号: SKLQ008)资助
摘要
建立了一种饱和-非饱和的双重孔隙-裂隙介质热-水-应力-迁移耦合三维模型,
关键词
双重孔隙-裂隙介质 热-水-应力-迁移耦合 三维模型
有限元分析
其特点是应力场和温度场是单一的, 但具有不同的孔隙渗流场、裂隙渗流场和孔隙浓
度场、裂隙浓度场, 以及可考虑裂隙的组数、间距、方向、连通率和刚度对本构关系 的影响, 并研制出相应的三维有限元程序. 通过与已有算例的对比, 验证了该模型和 程序的可靠性. 针对一个假定的高放废物地质处置库, 就岩体和缓冲层均为非饱和介 质和放射性核素泄漏的情况进行了数值分析, 考察了岩体中的温度、负孔隙水压力、 地下水流速、核素浓度和正应力的状态. 结果显示, 缓冲层中的温度、负孔隙水压力及 核素浓度呈现非线性的变化、分布; 尽管裂隙水饱和度平均仅为孔隙水饱和度的 1/9, 但因裂隙的渗透系数比孔隙的渗透系数大 4 个数量级, 故裂隙中地下水的流速约是孔 隙中相应值的 6 倍; 应力分布集中的区域位于缓冲层和处置孔壁交界两侧附近.
合课题的探讨[3~6]. 笔者也曾就国内外类似研究工作
中的缺点作了改进, 建立了一种双重孔隙-裂隙介质 热-水-应力耦合模型[7], 其特点是较严格地推导了应 力平衡、水连续性及能量守恒三大控制方程, 全耦合 求解, 用于饱和-非饱和的遍有节理的岩体, 可考虑 裂隙的组数、间距、方向、连通率、开度、刚度的影 响,并研制出相应的二维有限元程序. 但该模型也不 够完善, 即没有考虑核素迁移的耦合分析.
梁冰等人[8] 充分考虑了地下水变密度对核素浓 度分布的影响, 建立了孔隙-裂隙岩体中核素随地下 水迁移的耦合数值模型 ; Sonnenthal 等人 [9] 使用
1 引言
在自然界中, 类似于花岗岩和凝灰岩这样遍有
节理、同时具有岩石孔隙和岩体裂隙的介质可称为双 重孔隙-裂隙介质. 中国、日本、瑞典、瑞士和加拿大 等国家拟在花岗岩中进行高放射性核废物的地质处
置, 而美国已在内华达州的凝灰岩中完成了坑道规
模试验(DST), 并将在尤卡山(Yucca Mountain)地区建 成世界上第一个高放废物地质处置库[1, 2 ] . 近年来,
就双重孔隙-裂隙介质中存在的温度场、应力场、渗
流场的相互作用问题, 若干学者已开展了热-水-力耦
引用格式: Zhang Y J, Zhang W Q. 3D thermo-hydro-mechanical-migratory coupling model and femanalyses for dual-porosity medium. Sci China Tech Sci, 2010, 53:
2172−2182, doi: 10.1007/s11431-010-4031-3

此文档仅供学习和交流
中国科学: 技术科学
2010 年 第 40 卷 第 12 期
TOUGHREACT 程序模拟了尤卡山坑道规模试验中
双重介质内部的热-水-化学过程. 不过这些研究还局 限于水-迁移耦合、或热-水-迁移(化学)耦合. 武文华、 刘泽佳等人[10,11]曾提出和发展了多孔介质中的化学- 热-水力-力学耦合本构模型; 笔者也建立过一个饱和
-非饱和多孔介质中热-水-应力-迁移耦合模型[12]. 但 该类模型只适用于单重孔隙介质和平面应变问题.
为此, 针对上述研究中的不足之处, 笔者在自己
以及前人已有研究成果的基础上, 首先在理论上发 展了双重孔隙-裂隙介质中热-水-应力-迁移耦合的数 学模型, 其次研制出相应的三维弹性有限元程序, 然 后从方法论研究的角度, 以一个假定的高放废物地
质处置库为算例, 在一定的初始温度、孔隙水压力、
岩体应力和核素释放强度条件下, 就岩体是非饱和 双重孔隙-裂隙介质的情况, 考察了处置库近场的温 度、孔隙水压力、水流速、核素浓度和正应力的分布 与变化, 得出了若干的认识.
图 2 局部与整体坐标系
dσ = D ⎡ dε − m ⎛ C
⎞
− (s + D p
1
dp
) w1
⎜
⎝
⎟
⎢
1
w1 s1 w1
dt ⎣ dt
3Ks ⎠
dt

dT ⎤
⎛
⎞
β
1
dp
w2 − m S
−m ⎜ C
⎝
− ⎟ (s
+ D p )
⎥ , (1)
2
w2 s 2 w2
3K
dt
3 dt ⎦
s ⎠
这里, σ, ε 分别是总应力和总应变; D = (C1+C2)−1 为弹
性矩阵; mT =[1 1 1 0 0 0]为法向应力单位列阵; Ks, βS,
T 依次是岩质固体的体积模量、综合热膨胀系数和温 度; sw1, pw1, Ds1, C1 和 sw2, pw2, Ds2, C2 分别是孔隙岩体
及裂隙介质中的饱和度、水压力、湿气容量和柔度矩 阵; t 为时间. 不考虑负的孔隙水压力和裂隙水压力
对应力平衡的影响[14].
而且
2
双重孔隙-裂隙介质热-水-应力-迁移耦
合方程

对于图１所示的双重孔隙-裂隙介质, 认为其中
存在着孔隙水压力和裂隙水压力、孔隙浓度和裂隙浓
度, 但应力场和温度场是单一的, 从而建立了一种热-
∂s
∂s
D =
w1 ,
D =
w2 .
(2)
s1
s 2
水-应力-迁移耦合三维模型.
给出控制方程如下.
略去繁琐的推导过程,
∂p
∂p
w1
w2
对于三维问题, 有
−μ
1
−μ
0
0
0
−μ
−μ
1
0
0
0
2.1 应力平衡方程
借鉴了 Leiws 提出的在非饱和条件下有效应力 的计算方法[13], 并假定孔隙岩体中发育有 n 组裂隙 (见图 2), 在整体坐标系中有
⎡ 1
0
0
0
2(1 + μ)
0
0
0
0
0
0
2(1 + μ)
0
0
0
0
0
0
⎤
⎥
⎢−μ
⎢
⎥
= 1 ⎢−μ
⎥ , (3)
C
⎢
⎥
⎥
⎥
1
E ⎢ 0
⎢ 0
⎢
⎢⎣ 0
⎥
2(1 + μ)⎥⎦
= L C ′ L T ,
C
(4)
2i i 2 i i
⎡1/ Kni
0
0
0
0
0

0
0
0
0
0
0

0
0
0
0
1 / K
0

0
0
0
0
0
1 / Ksy ′i
0
0 ⎤
⎢
0 ⎥
0
0
0
0

0
⎢
⎥
A ⎢
0 ⎥
C ′ =
i ⎢
⎥ ,
(5)
2i
0
S
i ⎢
⎢
⎢
⎢⎣
⎥
sx′i
0 ⎥
⎥
0 ⎥⎦
图 1 孔隙-裂隙介质
1427

张玉军等: 三维双重孔隙-裂隙介质热-水-应力-迁移耦合模型及其有限元分析
l 2
l 2
l 2
⎡
2l l
2l l
2l l
⎤
⎥
⎥
1
2
3
1 2
2 3
3 1
⎢
2
2
2
⎢ m1
m2
m3
2m1 m2
2m2 m3
2m3 m1
⎢ n2
⎥
n2
n2
2n n
2n n
2n n
L = ⎢
1
2
l2 m2
3
l3 m3
1 2
2 3
3 1
⎥ .
(6)
i
⎢ l1 m1
l1 m2 + l2 m1
l2 m3 + l3 m2
l3 m1 + l1 m3 ⎥
⎢m n
m n + m n ⎥
m n
m n
m n + m n
m n + m n
⎢
3 1 1 3 ⎥
n3l1 + n1l3 ⎥⎦
1 1
2 2
3 3
n3l3
1 2 2 1
2 3 3 2
n2 l3 + n3l2
⎢⎣ n1l1
n2 l2
n1l2 + n2 l1
上述式子中, E, μ 分别为孔隙岩体的弹性模量和
泊桑系数; C2′i 为局部坐标系中的裂隙柔度矩阵; Li 为 坐标变换矩阵; Kni, Ksx′i, Ksy′i, Ai, Si 依次为裂隙的法向 刚度、沿 x′方向的剪切刚度、沿 y′方向的剪切刚度、 连通率和间距; {li, mi, ni, i = 1, 2, 3}为局部坐标轴与 整体坐标轴之间夹角的余弦.
⎛
⎞
m ⎟ ,
⎠
1
+ D p ) ⎜ C m −
s 2 w 2 2
3K
⎝
s
⎡
β
⎤
1
F1 = −ϕ1 sw1 βw + sw1 ⎢(ϕ1 −1)βs + 3K m D1m 3 ⎥
T
s
⎣
⎦
s
⎛ m mT −
T
D1 ⎞ (β m − C Dmβ ),
+ 1 s
w1 ⎜
⎝
⎟
s 1 S
3
3K
s ⎠
2.2 水连续性方程
根据质量守恒原理, 在 dt 时段内流入某一物体 的水量应等于其内部储水量的增加. 设定水的渗流 可以用达西定律来描述, 则对于孔隙介质有
式中, ϕ1, βs 分别为孔隙岩体的孔隙率和热膨胀系数;
Kw 和βw 分别是水的体积模量与热膨胀系数; D1 是孔 隙岩体的弹性矩阵.
对于裂隙介质,
水的连续性方程为
经过与上述类似的推导, 可得出
−∇T ⎧k
krw1 ∇（p
⎫ α k k
+ γ z) + 1 rw1 ( p − p )
α k k
⎨ 1 μ
⎬
⎭
⎧
k
⎫
∂ε
w1 w
w1 w2
μ
−∇T ⎨k rw2 ∇（p + γ z ⎬ −
1 rw1
( p − p
+ A
)
w1 w2 )
⎩
w
w
2
w2
w
2 ∂t
μ
μ
⎩ w
⎭
w
+ A ∂ε + B ∂pw1 + E
∂pw2
∂p
∂p
∂T
1 ∂t
1
1
∂t
∂t
w2 + E w1 + F
− ∇T D ∇T = 0 ,
+ B
(8)
2
2
∂t 2 ∂t
t 2
∂t
+ F ∂T − ∇T D ∇T = 0,
(7)
式中, k2, krw2 分别是裂隙介质的固有渗透率张量和比
渗透率; A2, B2, E2, F2 可由 A1, B1, E1, F1 表达式中的下 标“1”, “2” 分别改为“2” , “1”得到.
在本模型中, 不考虑应力场对裂隙作用所引起
的裂隙渗透张量的变化.
1 ∂t
t1
这里, k1, krw1 分别是孔隙岩体的固有渗透率张量和比
渗透率; ρw, μw, γw 依次为水的密度、粘滞系数和容重;
z 为位置水头; α 为取决于裂隙的开度与几何形状的
参数; Dt1 为孔隙岩体温度梯度水分扩散系数.
并且有
2.3 能量守恒方程
根据能量守恒原理, 在 dt 时段内流入某一物体 的热量应等于其内能的增加. 这里认为双重介质中 温度场是单一的, 可得
⎛
m T D ⎞
A = s m T −
1
C D ,
1 w1 ⎜
⎝
⎟
1
3K
s ⎠
⎛ 1 − ϕ1 −
⎞
m T D m (s
sw1 + s
1
B = D ϕ + ϕ
−∇T λ∇T + (s ϕ V a + s ϕ V a )ρ C (∇TT )
1 s1 1 1 K w1 ⎜
⎟
⎠
1
w1
w1 1 1 w2 2 2 w w
(3K )2
K
w ⎝ s
s
⎡
ρ
ρ
⎤
+ ⎢(1 − ϕ1 )CsT + ϕ1CwT
K K
s
w
⎥ (sw1 + Ds1 pw1 )
⎛
m T D ⎞
+ D p ) − s m T −
1
(s + D p )
⎣
w ⎦
s1 w1 w1 ⎜
⎝
⎟
s
w1 s1 w1
3K
s ⎠
∂p
× w1 − {(1 − ϕ )CT ρ β + (ϕ + ϕ )C T ρ β
1 s s s 1 2 w w w
⎡ ⎛
× C D C m −
⎞⎤
,
∂t
⎞ ⎛
m +
1
1
m − C1 m ⎟⎥
⎢ ⎜
⎟
⎜ 3K
∂T
∂t
1 1
3K
⎢⎣
⎝
⎠ ⎝ s
⎠⎥⎦
−[(1 − ϕ1 )ρs Cs + (ϕ1 + ϕ2 )ρwCw ]}
s
mT D ⎞
⎛
+ 1 (1 − ϕ )β T ∂ (u + u )δ = 0,
T
E1 = −sw1 ⎜ m −
⎝
1
⎟ C1 D(sw 2
(9)
3K
1 s
i , j j ,i ij
s ⎠
2
∂t
1428

中国科学: 技术科学
2010 年 第 40 卷 第 12 期
这里, Cw, Cs 分别是水及双重介质的比热; ρs, λ分别是
3 有限元格式

3.1 空间域离散
使用 Galerkin 方法, 对前述的应力平衡方程、水 连续性方程、能量守恒方程和渗透迁移方程进行空间
a a
双重介质的密度和导热系数矩阵; V1, V2 分别是孔隙
水和裂隙水的平均表观流速; ui, uj 是位移分量; δij 为 克罗内克符号.
2.4 渗透迁移方程
对参考文献[15]中渗透迁移方程进行了改进,
域离散,
整理可得
K du + C dpw1 + I dpw 2 + J dT = df ,
其
新意在于增加了孔隙岩体和裂隙介质之间因浓度差
产生的溶质交换, 其形式如下
dt
dt
dp
dt
dt dt
du dp

w1 + H
w 2 + F p + Mp − Mp
E
+ G
1
1
1
1 w1 w1 w 2
R θ ρ ∂ci = ∇Tθ ρ D ∇c − θ ρ V ∇c − R θ ρ
dt
dt
dt
χ c
i i w
i w i i i w i i i i w i
dT
∂t
+U1T + V1 = 0,
dt
dp
+(−1) ωθi ρ w D1 (c1 − c2 ) − Qc i ,
i +1
(10)
du
dp

E
+ G
w1 + H
w 2 + F p − Mp + Mp
2
2
2
2 w 2 w1 w 2
dt
dt
dt
式中, i = 1, 2 分别对应于孔隙岩体和裂隙介质;
+U T + V dT = X,
V ⎛ ρ
R = i = 1 + di K
⎞
,
(12)
其为阻滞系数;
V 为地下水
2 2
⎜
⎝
di ⎟
⎠
dt
i
i
∗
θi
Vi
du
dp
dT
A + P + RT + Q = Y , dt dt
w1

流表观速度; Vi*为溶质输运速度; ρdi 为孔隙岩体或裂
隙介质的干密度; Kdi 为对饱和介质的分配系数; θi 为 体积含水量; ρw 为流体密度; Di 为扩散张量; ci 为溶质 浓度; Vi 为地下水表观流速矢量; χ为衰减常数; ω 为
dt

dc
W c + L 1 + Sc − Sc = Z,
1 1 1
1 2
dt
W c + L dc2 + Sc − Sc = 0,
2 2 2
2 1
dt
取决于裂隙的开度与几何形状的参数; Qci 为源汇项.
而扩散张量可表示为
上述式子中, 未知量为 u, pw1 , pw2 , T , c1 , c2 ,
其前面的黑体大写字母为常数矩阵; f 为取决于作用
在单元上的体力、面力与初应力的力荷载; X, Y, Z 分 别为由源汇项和边界流量项求得的水荷载、热荷载和
溶质浓度“荷载”.
V V
) iα i β
iαβ = αiT
i δαβ + αiL − αiT
+ αimτ iδαβ ,
D
V
(
(11)
V
i
式中, αiT 为横向弥散度; αiL 为纵向弥散度; |Vi|为表观
流速的绝对值; αim 为分子扩散系数; τi 为曲折率; δαβ
为克罗内克符号.
在上述式子中, (1), (7)~(10)式即是热-水-应力-迁 移耦合控制方程.
3.2 时间域离散
认为在每个时步Δt 内 u, pw1 , pw2 , T , c1 , c2 是
线性变化的, 故可将(12)式变为如下的迭代求解方程
⎡ K
G
G1 + η (F1 + M )Δtk
G2 −η M Δtk
P
I
H1 −η M Δtk
H 2 + η (F2 + M )Δtk
0
J
⎤
⎧ U ⎫
⎪ ⎪
⎢
⎥
V1 + ηU1 Δtk ⎥
⎢ E1
⎪ Pw1 ⎪
⎨ ⎬
⎢
⎥
E2
⎢
V2 + ηU2 Δtk
⎥
⎪Pw2 ⎪
⎪⎩ T ⎪⎭
Q + η RΔt
⎣ A

⎡ K
k ⎦
t ,η
tk +Δtk
⎧ ∂f ⎫
G
G1 − (1 −η )(F1 + M )Δtk
G2 + (1 + η ) M Δtk
P
I
H1 + (1 −η ) M Δtk
H 2 − (1 −η )(F2 + M )Δtk
0
J
⎤ ⎧ U ⎫
⎪ ∂t ⎪
⎢
⎥ ⎪ ⎪
= ⎢ E1
V1 − (1 −η )U1 Δtk ⎥
⎪ Pw1 ⎪ + ⎪ 0 ⎪
Δt ,
(13)
⎨ ⎬ ⎨ ⎬
⎢
V2 − (1 −η )U2 Δtk
⎥
⎥
k
E2
⎢
⎪Pw2 ⎪ ⎪ X ⎪
⎩⎪ T ⎭⎪ ⎪ Y ⎪
Q − (1 −η ) RΔt
⎣ A
k ⎦
t ,η
⎪⎩ ⎪⎭
tk
t ,η
1429

张玉军等: 三维双重孔隙-裂隙介质热-水-应力-迁移耦合模型及其有限元分析
⎡ L1 + η (W1 + S)Δtk
−η SΔtk
⎤ ⎧ c1 ⎫
⎨ ⎬
⎢
⎣
⎥
−η SΔtk
L2 + η (W2 + S)Δtk ⎦t ,η ⎩c2 ⎭t +Δt
k k
⎡ L − (1 −η )(W + S)Δt
= 1 1 k
(1 + η )SΔtk
⎤ ⎧ c1 ⎫ + ⎧Z ⎫
(14)
Δt .
⎨ ⎬ ⎨ ⎬
⎢
⎣
⎥
k
(1 + η )SΔtk
L2 − (1 −η )(W2 + S)Δtk ⎦t ,η ⎩c2 ⎭t
⎩ 0 ⎭t ,η
k
在每一时步的迭代计算中, 若前后两次迭代所
得的孔隙水压力、裂隙水压力和相应的渗透系数的相 对误差均小于 0.01, 则转入下一时步的迭代, 直到完
成整个确定时段的分析.
4 双重孔隙-裂隙介质耦合模型的验证

对于上述的双重孔隙-裂隙介质热-水-应力-迁移 耦合模型, 笔者开发了相应的三维有限元程序. 耦合 模型及有限元程序的合理性验证如下.
由于目前尚无热-水-应力-迁移耦合模型的解析 解可供对比, 加之一时也难以找到合适的现场试验
数据来模拟, 故笔者试图用已有的数值解来验证本 模型. Elsworth 曾提出过一个双重孔隙介质的水-应力 耦合模型. 作为考题, 其假定: 一饱和柱状体中分布
有间距为 0.1 m 的正交裂隙, 在轴向承受均布荷载 P0,
材料的物理力学参数见表 1. 从而使用相应的有限元 程序对该柱状体进行了流-固耦合数值计算模拟, 其
结果符合客观规律[16]. 笔者为了将所建立的双重孔
隙-裂隙介质热-水-应力-迁移耦合模型与 Elsworth 的 成果相比较, 也针对上述柱状体进行了流-固耦合(温 度为常量, 暂不计迁移)的有限元分析. Elsworth 和笔 者得出的柱状体中孔隙流体压力和端面沉降的分布 及随时间变化的曲线见图 3 和 4. 看到 2 种程序的模 拟结果相当接近, 这表明了笔者所建立的模型是可 靠的.
图 3 4000 s 时柱状体中孔隙流体压力分布
图 4 柱状体表面沉降-时间曲线
5 算例

根据目前国外所设计的核废物地质处置概念库
的构成[14],
立足于方法论的定性研究, 笔者采用了一
个简单的计算模型, 见图 5. 该模型的尺寸为长×宽×
高=10 m×10 m×20 m, 其内部缓冲材料的尺寸为长×
宽×高=4 m×4 m×4 m, 核废物贮存罐位于模型的中心.
由于对称性, 取沿垂向剖开的 1/4 模型作为分析域, 对应的有限元网格见图 6, 其中有 4891 个节点, 4000 个八节点等参数单元. 假定核废物贮存罐的埋深为
500 m, 岩体中的初始温度为 20.0℃. 取边界条件为:
顶、底面水头及温度固定, 其余面为绝热及不透水面,
四个垂直面及底面的法线方向位移约束, 而顶面是
Table 1 Column coefficients
Young’s modulus, E (MPa)

Poisson’s ratio, μ
Fissure stiffness, kn( MPa/m) Fluid bulk modulus, kf (MPa) Matrix porosity, φ1
Fissure porosity, φ2
Matrix permeability, k1/μ (m2/Pa·s) Fissure permeability, k2/μ (m2/Pa·s) Fissure spacing, s (m)
1.0

0.15
0.1

0.1

0.1
0.05
0.01×10−9
0.1×10−6
0.1
位移自由的
, 其上并作用有上部岩体自重转化的节
点荷载. 岩体是非饱和介质, 其中发育有一组水平裂
隙及两组正交的垂直裂隙. 笔者参考了美国 DST 试
1430

中国科学: 技术科学
2010 年 第 40 卷 第 12 期
验[17]和日本 BMT1 试验[14]的情况, 取岩石、缓冲层、
固化体和裂隙组的主要计算参数见表 2 和 3.
孔隙岩体、裂隙介质和缓冲层的初始饱和度依次 为 0.9, 0.1 和 0.2, 其水分特性曲线符合 Van Genuchten
模型, 即
n −m
sw = (sws − swr )(1+ αψ
) + swr ,
(15)
式中, α, n, m = 1−1/n 为材料参数; ψ为水势; sws, swr 分
别为最大、最小饱和度. 对于孔隙岩体: α = 2.25×10−6
−1
m , n = 1.33, sws, swr 分别为 1.0 和 0.18, 从而孔隙岩
体中的初始孔隙水压力为−0.34 MPa; 对于裂隙介质:
α = 9.74 × 10−5 m−1, n = 1.97, sws, swr 分别为 1.0 和 0.01, 故裂隙介质中的初始孔隙水压力为−0.14 MPa; 对于 缓冲层: α = 8.0×10−3 m−1, n = 1.6, sws, swr 分别为 1.0 和
0.0, 所以缓冲层中的初始孔隙水压力为−18.12 MPa.
比渗透率与饱和度的关系为
图 5 核废物地质处置模型
(a) 立体图; (b) 剖面图(其中点号由内向外依次 为 2422, 2423, 2424, 2425)
= s2.0 .
k
(16)
rw w
取孔隙岩体、裂隙介质和缓冲层的温度梯度水分
扩散系数均为
D = 2.0 ×10−10 m2/s ℃.
(17)
t
假定核废料贮存罐埋置完毕后, 由于某种原因
贮存罐有所损坏, 并以恒定的浓度及速率均匀地向
近场泄漏放射性核素,
以此时为计算起点, 贮存罐
图 6 三维有限元网格
Table 2 Main computation parameters for rock, buffer and vitrified waste
Property
Rock mass
Buffer
Vitrified waste
Density, ρ ( kN/m3)
Porosity, φ1
Permeability, k1/μw (m2/Pa.s) Young’s modulus, E (MPa) Poisson’s ratio, μ
Specific heat, C (kJ/kg°C) Thermal expan. coeff., β (1/°C) Thermal conductivity, λ (W/m °C)
25.3
0.11
1.24×10−10
1.48×104
0.21
0.95
8.8×10−6
2.0
16.0
0.33
4.0×10−13
5.0×102
0.3
0.34
1.0×10−6
0.44
25.0
0.0
1.0×10−27
5.3×104
0.21
0.7
1.0×10−5
5.3
Table 3 Parameters for fracture sets used in calculation
Parameter
Horizontal fracture
Vertical fracture
Spacing, S (m)
Continuity ratio, A
Dip angle, θ (°)
Normal stiffness, kn (MPa/m) Shearing stiffness, ks (MPa/m) Porosity, φ2
Constant, α (m−2)
0.3
1
0
1000.0
500.0
0.01
100.0
1.0×10−6
0.3
1
90
2000.0
1000.0
0.01
100.0
1.0×10−6
Permeability, k2/μw (m2/Pa·s)
Note: the properties of two vertical fracture sets are the same, and the shearing stiffnesses in any direction on a fracture surface are the same.
1431

张玉军等: 三维双重孔隙-裂隙介质热-水-应力-迁移耦合模型及其有限元分析
为孔隙岩体中的源项, 某种放射性核素的释放强度
为 Qc = 1.44 × 10−6 mol·kg/m3·s−1. 设与核素的渗透 迁移有关且在计算中不变的参数为: 孔隙曲折率τ1, τ2 分别为 0.4, 0.8, 纵向弥散度α1L, α2L 分别为 1.0 和
2.0 m, 横向弥散度αiT = αiL/10, 分子扩散系数α1m,
α2m 分别为 1.0×10−9 和 2.0×10−9 m2/s, 分配系数 Kd1,
Kd2 分别为 8.0 和 5.3 mL/g, 干密度ρd1, ρd2 分别为 23.0
和 21.0 kg/m3 , 参数 ω 为 100.0 m-2, 放射性核素的衰
减常数χ = ln2/Thalf, 其中 Thalf 是半衰期, 取为 1000 a.
核废物以 300 W 的不变功率释放热量, 取时间步长为
100000 s, 经历了 4 a. 计算的主要结果及分析如下.
图 5(b)显示的剖面即为有限元网格的正面, 上面 标注了 4 个计算结果输出点的位置, 其中 2422 点位
于贮存罐的外缘, 2423 点、2424 点位于缓冲层的内部,
2425 点位于缓冲层和处置孔壁的交界处. 图 7 是上述
4 个点的温度随时间变化曲线, 从中看到, 各点的温 度在开始阶段快速上升, 离贮存罐越近的点其温度 越高, 在前 0.1 a 温度增加最明显, 之后温度增长缓
慢, 到 4 a 时 4 个点的温度依次为 72.5℃(2422 点)、
45.0℃(2423 点)、32.3℃(2424 点)和 24.8º℃(2425 点).
图 8 是 4 a 时在模型高度之半的水平剖面上贮存罐周 围 2.5 m×2.5 m 范围内的温度等值线分布.
图 9 是 2422, 2423, 2424 和 2425 各点处的孔隙水 压力随时间的变化曲线. 从中得知, 2425 点位于缓冲
层和处置孔壁的交界处, 该点得到地下水的快速浸
润, 其负孔隙水压力在约 0.46 年时降到了−0.25 MPa.
而位于缓冲层内部的 3 个点的负孔隙水压力开始均 有所上升, 其原因在于: 一是在开始阶段热-水-应力
耦合现象强烈, 缓冲材料受到围岩的挤压作用后其 含水量减少(特别是在离围岩较近之处, 如 2424 点),
二是在开始阶段缓冲层中距离贮存罐近的部位的温
度升高较快, 温度梯度引起了水分扩散(向外), 而水
图 8 4 a 时缓冲层水平剖面上温度等值线(℃)
图 9 缓冲层中孔隙水压力-时间曲线
力梯度产生的水分迁移(向内)尚未到达该部位, 使得
该部位的干燥程度不断增加. 一定时间后随着地下 水分的到达和浸润作用的增强, 缓冲层中各点的负 孔隙水压力才逐渐变小. 到 4 a 时 4 个点的负孔隙水
压力依次为−0.53 MPa(2422 点)、−0.44 MPa(2423 点)、
−0.22 MPa(2424 点)和−0.13 MPa(2425 点).
图 10 是 0.8 a, 2.4 a 和 4 a 时在模型高度之半的 水平剖面上缓冲层区域中的孔隙水压力等值线分布.
图 11 是 0.8 a 和 2.4 a 时在模型高度之半的水平 剖面上围岩区域中的孔隙水压力和裂隙水压力等值
线分布. 从中看到, 随着时间的推移, 孔隙水压力和
裂隙水压力的变化趋于平缓,其差值逐渐减小. 如在 沿水平方向距贮存罐外缘 2.0 m 处,在 0.8 a 和 2.4 a
时的孔隙水压力和裂隙水压力分别为: −0.130 MPa,
−0.128 MPa; −0.120 MPa, −0.121 MPa. 图 12 是 0.8 a
时该水平剖面上的孔隙水流速和裂隙水流速矢量分 布(后者的比例尺为前者的 1/200000). 如在沿水平方
图 7 温度-时间曲线
1432

中国科学: 技术科学
2010 年
第 40 卷
第 12 期
图 10 缓冲层水平剖面上孔隙水压力等值线 (MPa)
(a) 0.8 a 时; (b) 2.4 a 时; (c) 4 a 时
向距贮存罐外缘 2.0 和 3.0 m 处, 孔隙岩体和裂隙介
质中的流速分别为 3.77×10−9, 1.34×10−3, 5.54×10−10
和 1.55×10−4 m/s, 可见裂隙中的流速要比裂隙中的流 速大 6 个数量级左右.
图 13 是上述缓冲层中 4 个点处的核素浓度随时 间的变化曲线. 看到位于贮存罐外缘的 2422 点处的
核素浓度最大, 在 4 a 的时间内基本呈直线上升, 而
其余 3 个点处的核素浓度依点距贮存罐的近、远而逐 次降低, 并且随时间显示了非线性的的变化, 这表明
近场的热-水-应力耦合过程是复杂的, 核素迁移要受
其影响. 到 4 a 时 2422, 2423, 2424 和 2425 各点的核 素浓度依次为: 1.26, 0.28, 0.07 和 0.06 mol/m3. 图 14
是 4 a 时在模型高度之半的水平剖面上贮存罐周围
2.5 m × 2.5 m 范围内的核素浓度等值线分布.
在上述水平剖面上 4 a 时的正应力σx, σy 和σz 的等 值线见图 15. 从中看到, 应力分布集中的区域位于缓 冲层和处置孔壁交界两侧附近, 该面上最大的σx, σy
和σz 值依次为−2.47, −2.47 和−16.6 MPa(均为压性).
6 结束语

对于饱和-非饱和的遍有节理岩体, 建立了一种 双重孔隙-裂隙介质的热-水-应力-迁移耦合三维模型, 其中岩体的应力场和温度场是单一的, 但孔隙岩体 和裂隙介质具有不同的渗流场和浓度场, 并开发了 相应的三维有限元计算程序.
对一个假定的高放废物地质处置库进行了三维 数值模拟. 其计算结果表明: 缓冲层中温度在开初的
0.1 a 内快速上升, 之后变化很小, 计算终了(4 a)时,
缓冲层内的温度可达到