第28卷 增1
岩石力学与工程学报 V ol.28 Supp.1
2009年5月 Chinese Journal of Rock Mechanics and Engineering May ，2009
层状盐岩体非线性蠕变本构模型
王安明
1，2
，杨春和，陈剑文3，黄志全1，李银平2
3. 中国地质大学，湖北 武汉 430074)
2
(1. 华北水利水电学院，河南 郑州 450011；2. 中国科学院武汉岩土力学研究所 岩土力学与工程国家重点实验室，湖北 武汉 430071；
摘要：针对由不同岩层交替而组成的层状盐岩体，建立泥岩夹层和盐岩复合体代表单元，根据应变协调原理，从细观力学分析角度通过考虑泥岩夹层弹性性质、盐岩弹性及蠕变力学特性以及两相体积含量建立层状盐岩体宏观各向异性非线性蠕变增量型本构模型，分析层状盐岩在蠕变过程中因保持细观应变协调而产生的应力重分布问题，给出新本构模型ABAQUS 有限元二次开发增量迭代算法实现方法，对一层状盐岩体简单试样算例进行初步计算分析验证，分析计算与试验结果吻合较好，表明该本构模型能反映细观应变协调时层状盐岩体宏观蠕变力学特性，该模型将为层状盐岩体内硐室长期稳定性分析提供理论计算基础。
关键词：岩石力学；层状盐岩；应变协调；细观力学；非线性蠕变；应力重分布
中图分类号：TU 45 文献标识码：A 文章编号：1000–6915(2009)增1–2708–07
NONLINEAR CREEP CONSTITUTIVE MODEL OF BEDDED SALT ROCK
WANG Anming12，YANG Chunhe2，CHEN Jianwen3，HUANG Zhiquan1，LI Yinping2
，
(1. North China Institute of Water Conservancy and Hydroelectric Power，Zhengzhou ，Henan 450011，China ；2. State Key Laboratory of Geomechanics and Geotechnical Engineering，Institute of Rock and Soil Mechanics，Chinese Academy of Sciences，
Wuhan ，Hubei 430071，China ；3. China University of Geosciences，Wuhan ，Hubei 430074，China )
Abstract ：As for alternating layers of mudstone and salt rock，representative volume element(RVE) in the composite is established. Based on strain compatibility principle and micromechanical method，considering elastic behavior of mudstone interbeds，elastic and creep behavior of salt，an anisotropic incremental nonlinear composite creep constitutive model of bedded salt rock is derived；stress redistribution is analyzed during creep deformation of salt. The new model is introduced into finite element software ABAQUS. A simple bedded salt rock sample is modeled as case study to verify the composite model. The result agrees well with the experiment，so it is shown that the constitutive model can reflect mechanical behavior of bedded salt rock considering micro-strain compatibility. The model can be employed to analyze long-term stability of bedded salt rock cavern in the future. Key words：rock mechanics；bedded salt rock；strain compatibility；micromechanical method；nonlinear creep；stress redistribution
1 引 言
盐岩具有良好的密封性、低渗透性、良好的蠕
收稿日期：2008–03–13；修回日期：2008–07–09
变行为和损伤自愈合性，使得盐岩构造已被国际上公认为石油、天然气储存和高放射性废物处置的理想场所。我国层状盐岩的基本赋存特点是：盐岩层数多，单层厚度薄，盐岩体中一般含有众多夹层(如
基金项目：国家重点基础研究发展规划(973)项目(2009CB724600)；国家杰出青年科学基金项目(50725414)；河南省高校杰出科研人才创新工程(2005KYCX015)；华北水利水电学院高层次人才基金项目(200901)；
作者简介：王安明(1975–) ，男，博士，1999年毕业于中南工业大学交通土建专业，现为讲师，主要从事层状盐岩力学特性，地下工程稳定性方面的教学与研究工作。E-mail ：[email protected]
第28卷 增1 王安明，等. 层状盐岩体非线性蠕变本构模型 • 2709 •
泥岩层等) ，在这种特殊地层中建腔进行长期能源储存，国际上没有先例可循，国外储气库的建设经验也不能直接借鉴到我国的储气库的建设中。M. Langer 和S. Hensermann[1]认为储藏库建设的核心是围岩的模拟，这样的模拟应该是尽可能的接近现场的条件包括地质条件，应力条件和本构模型。因此针对我国层状盐岩矿层的基本特点，杨春和和李银平[2]指出建立反映层状盐岩体本构模型及长期强度方程是层状盐岩废弃溶腔油(气) 地下可储性研究中的关键理论问题。
对层状盐岩溶腔进行稳定性分析时，由于盐岩体中含夹层众多，要考虑并模拟其中的每一小层是相当困难的，为避免给建模及计算带来困难，常常采用非均匀介质均匀化等效概化方法，即将厚度很小的夹层与盐岩层概化在一起作为一种等效均匀介质参与分析计算，而且这种概化后等效均匀介质的宏观力学特性与其细观相的力学特性和相对含量是密切相关的。杨春和等
[2，3]
图1 层状盐岩复合体代表单元
Fig.1 Representative volume element for bedded salt rock
(1) 互层层状盐岩介质是由力学性质不同的介质黏结而成，在界面处完美连接，没有滑动、张开或者嵌入；
(2) 细观上夹层和盐岩为各向同性连续均匀介
质，这样考虑下的复合体宏观上为均匀各向异性体(横观各向同性体) 。 2.2 应变协调原理
以σij ，εij 分别表示代表单元宏观平均应力和平
s s ，εij 分别表示盐岩的平均应力和平均均应变；以σij
m m
应变，以σij ，εij 分别表示代表泥岩夹层的平均应
建立的层状盐岩体宏观
平均意义下Cosserat 介质扩展本构模型即体现出这个特点，此模型在层状岩盐溶腔的静力分析时为层状盐岩体变形破坏提供了一条新思路；同理在对层状岩盐容腔进行长期稳定性分析时，也需要对层状盐岩建立等效概化的蠕变本构模型，最新的试验研究表明，层状盐岩的蠕变呈现出明显的非线性性质[4]，因此建立的本构模型尽可能的反应出这种非线性性质。刘 江等[5
，6]
力和平均应变，以V s ，V m 分别表示盐岩和泥岩夹层的体积含量，其中V s +V m = 1，由应变协调原理和多相混和物本构关系的一般原理[7形式之间的关系：
s m
d σij =V s d σij +V m d σij (ij = 11，33，13) (1) s m d εij =d εij =d εij (ij = 11，33，13) (2) s m d εij =V s d εij +V m d εij (ij = 22，32，12) (3) s m d σij =d σij =d σij (ij = 22，32，12) (4)
～9]
可以列出单元体
宏观平均应力–应变与细观上各相应力–应变增量
用黏弹塑性蠕变理论来建立层状(互
层) 盐岩蠕变本构方程，这些本构模型多半为经验型本构模型，没有反映出层状(互层) 盐岩本构模型的宏细观特性，也没有体现出层状盐岩蠕变的非线性性质。
本文用细观力学分析方法通过考虑泥岩夹层弹性性质、盐岩弹性及蠕变力学特性和各自相对体积含量建立层状盐岩体等效非线性蠕变增量型本构模型，最后给出了本构模型的二次开发实现方法并进行了初步的试验验证。
2.3 层状盐岩复合体中各相应力–应变关系
为了求出层状盐岩复合体应力–应变关系，必须先知道各相应力–应变关系，杨春和等[2
，3]
在建立
层状盐岩复合弹性本构模型时假定泥岩夹层和盐岩均为弹性，本文在此基础上建立复合体蠕变本构关系时引入盐岩层的蠕变项。试验研究[5]表明，泥岩夹层和盐岩均具有较为明显的蠕变特性，但泥岩的稳态蠕变率比盐岩层低一个数量级左右，层状盐岩的蠕变性质主要来自于盐岩部分。因此这里假定泥岩夹层仅为弹性变形，忽略泥岩夹层的蠕变，而盐岩层的变形由弹性加蠕变复合而成，这样一种假定大大简化了层状盐岩蠕变性质的分析。
对泥岩夹层有：
2 层状盐岩体增量型复合蠕变本构
模型
2.1 基本假定
选取如图1所示盐岩和泥岩层组成的层状盐岩复合体代表单元(RVE)，由连续介质力学理论推导复合体等效本构关系时采用以下假设：
• 2710 • 岩石力学与工程学报 2009年
m m
d σij =D ijkl d εkl (5)
对盐岩层有：
s s scr d εij =D ijkl (dεkl −d εkl ) (6)
⎛2E µs ⎞2E s 3n −1
s 11q n −1+s q ∆tA 1⎜s 2⎟ (10) 233
41µ1µ−−s s ⎝⎠
s m
d σ22=d σ22=d σ22 (11)
由式(5)得
m d σij µm m m m m
d εij =(dσ11+dσ33) δij (7) −+d σ22
m m 2G E
由式(6)得弹性变形：
s d σij µs s s s s
d εij =(dσ11+dσ33) δij −+d σ22
s s 2G E 式中：D ijkl 为刚度矩阵；G m ，E m 分别为泥岩夹层的剪切模量和弹性模量；G s ，E s 分别为盐岩层的剪切模量和弹性模量；µm ，µs 分别为泥岩夹层和盐岩层的泊松比。对于δij ，当i =j 时，δij = 1；当i ≠j 时，δij = 0。
研究结果[10]表明，盐岩蠕变变形主要来自于其稳态蠕变，盐岩的稳态蠕变率可以表示成
ij =ε
s ij 3
A 1q n
q 2
m
d σ11=
µm E m E m µm
d ε+d ε+d σ22 (12) 113322
1−µm 1−µm 1−µm
µm E m µm E m
d ε+d ε+d σ22 (13) 113322
1−µm 1−µm 1−µm
s s s
d σ13=2G s d ε13−3G s 13∆tAq n =2G s d ε13−
q
s
3G s 13∆tA 1q n (14)
q
m d σ33=
m m d σ13=2G m d ε13=2G m d ε13 (15)
s m
d σij =d σij =d σij (ij = 32，12) (16)
以上各式表明了层状盐岩中细观上盐岩和泥岩分担应力增量与层状岩盐宏观应力应变增量之间的关系，显然随时间步不同，盐岩和泥岩分担应力增量也在不断调整当中，因此发生了应力重分布。
式中：A 1，n 均为材料参数；q 为Mises 应力；s ij 为偏应力。
稳态蠕变为
cr
=d εij
s ij 3
A 1q n ∆t
q 2
s
d σij
4 层状盐岩体增量型蠕变本构模型矩
阵形式
消去岩盐层和泥岩夹层的应力–应变增量之后，可以得到层状盐岩体宏观等效应力–应变关系，写成矩阵形式为 ⎧d ε11⎫⎡M 11⎪d ε⎪⎢M ⎪22⎪⎢21⎪d ε33⎪⎪⎢M 31⎪⎨⎬=⎢d ε⎪12⎪⎢0⎪d ε23⎪⎢0⎪⎪⎢d ε⎢0⎪13⎭⎪⎣⎩
M 12
M 22M 32000
M 13M 23M 33000
000M 4400
0000M 550
0⎤0⎥⎥0⎥⎥⋅ 0⎥0⎥⎥M 66⎦⎥
盐岩总应变为
d ε=d ε+d ε=
s
ij
e ij
cr ij
−(dσ+d σ+ 2G s E s 3s ij s
∆tA 1q n (8) d σ33) δij +2q
s 11
s 22
µs
式中：∆t 为时间增量。
3 层状盐岩复合体中盐岩层和泥岩夹
层分担应力
联立上面的式(1)～(8)并经过推导得到
E s E s µs µs s
d σ11=d +d +d σ22− εε1133
1−µs 21−µs 21−µs
⎛2E s ⎞2E s µs 3n −1n −1∆tA 1⎜+s q s q ⎟ (9) 211233
−−41µ1µs s ⎝⎠
⎧d σ11⎫⎧F 11⎫
⎪d σ⎪⎪F ⎪⎪22⎪⎪22⎪⎪⎪d σ33⎪⎪⎪⎪F 33⎪⎪
=⎨⎬⎨⎬ (17)
⎪d σ12⎪⎪F 12⎪⎪d σ23⎪⎪F 23⎪⎪⎪⎪⎪d σ⎪⎩13⎪⎭⎪⎩F 13⎪⎭
其中，
s
d σ33=
E s µs E s
d ε+d ε3311
1−µs 21−µs 2
s
d σ33=
µs E s µs E s
εεd +d +d σ22−1133
1−µs 21−µs 21−µs
⎛⎞
⎜⎟1⎜11⎟ +M 11=
2⎜s s s s m m m m ⎟++⎜1−µ1−µ1+µs 1+µm ⎟s m ⎝⎠
第28卷 增1 王安明，等. 层状盐岩体非线性蠕变本构模型 • 2711 •
V s µs V µ
+m m
1−µs 1−µm
M 12=
s s +m m
1−µs 1−µm
⎛⎞⎜⎟1⎜11⎟ M 31=−2⎜s s ⎟s s m m
+m m ⎟⎜1−µ+1−µ1+µs 1+µm ⎠s m ⎝
V s µs V µ
+m m
1−µs 1−µm =−
s s +m m
1−µs 1−µm
⎛⎞
⎟1⎜11⎜⎟ M 13=−
s s m m m m ⎟2⎜s s
++⎜1−µ1−µ1+µs 1+µm ⎟s m ⎝⎠⎛1
(q n −1s 11−q n −1s 33) ⎜1+µs 3
+ F 11=E s V s A 1∆t ⎜
E V 4⎜s s
+m m
⎜1+µs 1+µm ⎝
1⎞(q n −1s 11+q n −1s 33) ⎟1−µs
⎟
s s m m ⎟
+⎟1−µs 1−µm ⎠V s µs V µ+m m
1−µs 1−µm
M 21=−
s s +m m
1−µs 1−µm M 22=
M 32
⎛⎞
⎟1⎜11⎟ M 33=⎜+
s s V m E m V m E m ⎟2⎜s s
+⎜1−µ+1−µ11+µm ⎟µ+s m s ⎝⎠⎡1
(q n −1s 11+q n −1s 33) ⎢1−µs 3
F 33=E s V s A 1∆t ⎢− s s m m 4⎢+⎢+µ11+µm s ⎣
1⎤
(q n −1s 11−q n −1s 33) ⎥1+µs
⎥
s s m m ⎥+⎥1−µs 1−µm ⎦
M 44=
V s V m ⎛µs V m µm V m
+−⎜−E s E m ⎝E s E m ⎞
⎟⋅ ⎠
V s V m
+ 2G s 2G m
3s
F 12=V s 12∆tA 1q n
2q M 55=M 66
V s V m
+ 2G s 2G m
1
= 2V s G s +2V m G m
⎡E m ⎛2V s µs 2V m µm ⎞⎤
+⎢⎥⎜⎟
⎢1−µm ⎝1−µs 1−µm ⎠−2µm ⎥ ⎢s s m m 1−µm ⎥
+⎢⎥11µµ−−s m ⎢⎥⎣⎦V s µs V µ+m m
1−µs 1−µm =−
s s +m m
1−µs 1−µm
M 23
3s
F 23=V s 32∆tA 1q n
2q s 13
∆tA 1q n q
F 13=
2V s G s +2V m G m
3V s G s
M 11=M 33
F 22=A 1∆t ⋅ ⎡⎛µs V m µm V m
−⎢⎜E E m ⎢⎝s ⎢⎢⎢⎣
⎤
⎥3n −1
V s q s 22⎥ 2⎥
⎥⎦
⎞3E s E m V s
(q n −1s 11+q n −1s 33) ⎟
⎠2(1−µs )(1−µm ) + V s E s V m E m
+
1−µs 1−µm
M 12=M 21M 13=M 31M 12=M 32M 21=M 23M 44=M 55
式(17)即为层状盐岩复合体宏观增量型应力–应变关系，简写成下面的统一形式：
d ε=M d σ+F (18)
• 2712 • 岩石力学与工程学报 2009年
式中：M 为常数项，与泥岩夹层和盐岩的弹性常数及其各自体积含量有关；F 为蠕变增量项，与泥岩夹层和盐岩的弹性常数、各自体积含量、盐岩层当前应力及时间步长有关。
显然宏观本构模型中含有其细观成分盐岩层当前应力及时间步长，从式(9)～(16)可以看出，随蠕变时间变化，盐岩层分担应力在不断调整之中，可以看出此本构模型能够反映层状盐岩在蠕变过程中因保持细观应变协调而产生的应力重分布问题，体现了蠕变的非线性性质
[11]
(4) 第四步：ABAQUS 的平衡迭代开始于新时刻t n +1，ABAQUS 对于最大迭代次数的默认设置为
9，最小容差TOL 为5×106，当通过小于9次的迭
－
代后，如果误差小于TOL ，则为收敛，n = n + 1时进行下一步增量计算；如果误差大于TOL ，ABAQUS 将会缩减该步时间增量∆t 直到计算结果收敛为止。
6 算例验证与讨论
为验证新本构模型的有效性，考虑由泥岩夹层和盐岩组成的层状岩盐复合体，在ABAQUS 软件中建立圆柱体复合体模型，采用C3D8I 单元，底边施加固定边界条件，按照上述二次开发算法实现方法将本文建立的非线性蠕变本构方程通过UMAT 子程序植入ABAQUS 有限元软件。对某层状盐岩试样的三轴压缩蠕变试验[5](围压为5 MPa，差应力为10 MPa)进行初步的计算模拟，盐岩和泥岩夹层力学参数和体积含量如下面给出：取盐岩弹性模量
；另外对于互层盐岩
体，其宏观上呈现横观各向同性性质，显然上述建立的本构模型也体现了这样一种各向异性性质；如果不考虑∆t ，即不考虑蠕变项的话，上述本构模型退化为增量形式的弹性本构模型，这时与李银平和杨春和建立的层状盐岩体三维Cosserat 介质扩展本构模型(去掉偶应力项) 本质上是一致的，这也反映了本文本构模型的合理性。
[3]
5 新本构模型增量迭代算法及其在
ABAQUS 中的实现
在给定远场外载荷的情况下，层状盐岩体宏观应力应变关系可以通过上述建立的增量型应力应变关系逐次进行增量跌代方法得到，参考叶 峥等[12
，13]
E = 5 150 MPa，泊松比µs = 0.31，夹层弹性模量E =
21 780 MPa，泊松比µm = 0.23，夹层含量为44.1%，其盐岩层稳态蠕变本构关系符合Norton 蠕变模型，参数[5]取为A 1=1.12×108/h，n = 3.8，图2给出了
－
理论预测结果与试验结果对比曲线。
从图2中可以看出，数值预测与试验结果比较吻合，这说明本文非线性蠕变本构模型可较好地模拟层状岩盐的蠕变变形。当然蠕变计算值与试验值相比也有一定误差，其原因是多方面的[2
，3]
的方法，将上述本构模型和迭代算法编写成子程序
UMAT 嵌入大型有限元软件ABAQUS 进行二次开发，下面简述其主要步骤。将节4建立的本构模型式(18)简写成下面的形式：
1d σ=(dε−F ) (19)
M
，首先
是岩体材料力学特性的离散性，参数选取难度较大，然后是理论和试验的模型以及边界条件是有区别的，对于层状岩盐试样，现场钻孔时为圆柱形试样，而本文理论模型建立时以长方体形状的代表单元为起点进行推导的，复合体为横观各向同性性质，表明试验试样与理论模型存在不一致性，另对于三轴试验，虽然侧向压力相同，由于各层力学特性不同，其侧向变形是不相同的，这表明试验条件并不完全满足理论分析中的假设(1)，而理论分析时假设层面方向是无限展开的，故而各层侧向位移是相同的，但其侧向压力是不相同的。当然，刘 江[5]只是做了一个初步的试验，试验岩样是现场钻取的含小夹层的盐岩岩芯，百分比含量不易较大范围变化，系统地开展层状盐岩蠕变试验工作量大等，当前只验证的是层状盐岩其垂直于层面方向的蠕变力学性质，
(1) 第一步：由层状盐岩初始应力场，代入式(9)～(16)求得没有蠕变时盐岩和泥岩的初始分担应力。
(2) 第二步：在平衡时刻t n ，ABAQUS 提供总应变增量∆εtotal (t n ) ，总应变εtotal (t n ) 和时间增量∆t ，由盐岩和泥岩的当前应力(在前一步已经计算出来) ，进而可以求出时间步长∆t 内蠕变应变增量F 项的值，代入上述本构关系式(19)之后可以求出∆t 时段内总应力增量d σ的值。
(3) 第三步：由总应变增量、总应力增量d σ及当前分析步的当前应力代入式(9)～(16)，得到当前时步盐岩和泥岩的应力增量d σs ，d σm ，与盐岩和泥岩当前应力进行叠加后可以得到盐岩和泥岩下一
时步初始应力σs (t n +1) ，从而进行应力更新。
第28卷 增1 王安明，等. 层状盐岩体非线性蠕变本构模型 • 2713 •
累，即把非线性蠕变计算进行线性化处理[11
，15]
，运
用非线性问题增量迭代算法将本文本构模型写成子程序与ABAQUS 软件连接，为层状盐岩体容腔长
蠕变/10
－4
期稳定性分析奠定了基础。
15
1 3 5 7时间/(104 s)
91113
参考文献(References)：
[1] LANGER M ，HEUSERMANN S. Geomechanical stability and
integrity of waste disposal mines in salt structures[J]. Engineering Geology ，2001，21(2/3)：155–161. [2]
杨春和，李银平. 互层盐岩体的Cosserat 介质扩展本构模型[J]. 岩石力学与工程学报，2005，24(23)：4 226–4 232.(YANG Chunhe，LI Yinping. Expanded Cosserat medium constitutive model for laminated salt rock[J]. Chinese Journal of Rock Mechanics and Engineering ，2005，24(23)：4 226–4 232.(in Chinese)) [3]
李银平，杨春和. 层状盐岩体的三维Cosserat 介质扩展本构模型[J]. 岩土力学，2006，27(4)：509–513.(LI Yinping，YANG Chunhe. Three-dimensional expanded Cosserat medium constitutive model for laminated salt rock[J]. Rock and Soil Mechanics，2006，27(4)：509–513.(in Chinese)) [4]
郤保平，赵阳升，赵金昌，等. 层状盐岩温度应力耦合作用蠕变特性研究[J]. 岩石力学与工程学报，2008，27(1)：90–96.(XI Baoping ，ZHAO Yangsheng，ZHAO Jinchang. Study on coupled thermo-mechanical creep properties of bedded rock salt[J]. Chinese Journal of Rock Mechanics and Engineering，2008，27(1)：90–96.(in Chinese)) [5]
刘 江. 层状盐岩力学特性试验研究及理论分析[硕士学位论文][D]. 武汉：中国科学院武汉岩土力学研究所，2006.(LIU Jiang. Experimental investigation and theoretic analysis on the mechanical properties of layered rock salt[M. S. Thesis][D]. Wuhan：Institute of Rock and Soil Mechanics，Chinese Academy of Sciences，2006.(in Chinese)) [6]
郤保平. 含高盐份泥岩夹层的盐岩蠕变特性及油气储库稳定性研究[硕士学位论文][D]. 太原：太原理工大学，2006.(XI Baoping. Study on creep property of rock salt with salt-mudstone interlayer and stability of oil and gas storage cavern[M. S. Thesis][D]. Taiyuan：Taiyuan University of Technology，2006.(in Chinese)) [7]
范镜泓，高芝晖，刘天模. 具层状细观结构的PD3珠光体钢的细/宏观力学分析[J]. 重庆大学学报，1999，22(5)：37–46.(FAN Jinghong ，GAO Zhihui，LIU Tianmo. The micro/macro mechanical
图2 层状岩盐蠕变理论预测值与试验值对比曲线 Fig.2 Comparison between simulated values and experimental
results based on creep theory of bedded salt rock
别的方向的试验验证尚需钻取更多与理论模型比较一致的试样来进行进一步的试验。
此外，对于本文本构模型，当t = 0时，层状盐岩复合体呈弹性响应，盐岩层在任意时间步t +∆t 发生蠕变，但对于深部地下层状盐岩容腔加内压后，其围岩力学响应盐岩层往往是弹塑性的，但由于盐岩弹塑性变形相对其蠕变变形来说很小，况且在容腔设计中，与盐岩的高应变率弹塑性响应相比，盐岩的黏弹性响应是最为重要的[14]，从这一点来说本文模型因没有考虑塑性变形导致的误差影响不大，但是要建立能考虑塑性变形的复合本构模型仍需进一步深入研究。
7 结 论
本文针对由不同岩层交替而组成的层状盐岩体，建立泥岩夹层和盐岩复合体代表单元，根据应变协调原理，用细观力学分析方法通过考虑泥岩夹层弹性性质、盐岩弹性及蠕变力学特性以及两者体积含量建立层状盐岩体宏观非线性蠕变增量型本构模型，将细观上不均匀的两相介质在宏观上进行了均匀化处理，分析层状盐岩在蠕变过程中因保持细观位移协调而产生的应力重分布问题。该模型既能考虑由于泥岩夹层的存在而导致的层状盐岩构造体宏观平均意义下力学特性的各向异性，又能很好地反映层状盐岩蠕变的非线性性质，与试验研究结论非常一致[4]。
在数值计算中把蠕变经历的时间分成有限间隔，假定在每一个时间步内应力不变，可按某一时刻初的应力水平进行蠕变，这样就可以把变应力情况求蠕变量问题看作是求逐段常应力蠕变量的积
• 2714 • 岩石力学与工程学报 2009年
analysis of PD3 pearlitic steel with layered microstructures[J]. Journal of Chongqing University(Natural Science Edition)，1999，22(5)：37–46.(in Chinese))
[8] HAJ-ALI R M ，MULIANA A H. A multi-scale constitutive
formulation for the nonlinear viscoelastic analysis of laminated composite materials and structures[J]，International Journal of Solids and Structures，2004，41(13)：3 461–3 490. [9]
HAJ-ALI R M，MULIANA A H. A micromechanical constitutive framework for the nonlinear viscoelastic behavior of pultruded composite materials[J]. International Journal of Solids and Structures ，2003，40(5)：1 037–1 057.
[10] 陈 锋. 盐岩力学特性及其在储气库建设中的应用研究[博士学
位论文][D]. 武汉：中国科学院武汉岩土力学研究所，2006.(CHEN Feng. The mechanical behavior of salt rock and its application to the natural gas storage construction[Ph. D. Thesis][D]. Wuhan：Institute of Rock and Soil Mechanics，Chinese Academy of Sciences，2006.(in Chinese))
[11] 彭 凡，彭 献，傅依铭. 纤维增强复合材料非线性蠕变流动特
性[J]. 湖南大学学报(自然科学版) ，1999，26(5)：6–9.(PENG Fan，
PENG Xian，FU Yiming. Behavioral analyses of creep flow for fiber composites[J]. Journal of Hunan University(Natural Sciences)，1999，22(5)：6–9.(in Chinese))
[12] 叶 峥. 基于微观滑移的蠕变本构模型及应用研究[硕士学位论
文][D]. 北京：北京航空航天大学，2007.(YE Zheng. Slip theory based modeling of creep deformation and application in Ni-based single crystal super alloys[M. S. Thesis][D]. Beijing：Beijing University of Aeronautics and Astronautics，2007.(in Chinese)) [13] 罗海安，翟 军. 兼具基体蠕变和轻度界面蠕变的复合材料的本
构特性[J]. 力学学报，1996，28(5)：550–557.(LUO Haian，ZHAI Jun. Constitutive behavior of the composites with both matrix creep and slight interface creep[J]. Acta Mechanica Sinica，1996，28(5)：550–557.(in Chinese))
[14] LIANG W ，YANG C，ZHAO Y，et al. Experimental investigation of
mechanical properties of bedded salt rock[J]. International Journal of Rock Mechanics and Mining Sciences，2007，44(3)：400–411. [15] 穆霞英. 蠕变力学[M]. 西安：西安交通大学出版社，1990.(MU
Xiaying. Creep mechanics[M]. Xi′an ：Xi ′an Jiaotong University Press ，1990.(in Chinese))
第28卷 增1
岩石力学与工程学报 V ol.28 Supp.1
2009年5月 Chinese Journal of Rock Mechanics and Engineering May ，2009
层状盐岩体非线性蠕变本构模型
王安明
1，2
，杨春和，陈剑文3，黄志全1，李银平2
3. 中国地质大学，湖北 武汉 430074)
2
(1. 华北水利水电学院，河南 郑州 450011；2. 中国科学院武汉岩土力学研究所 岩土力学与工程国家重点实验室，湖北 武汉 430071；
摘要：针对由不同岩层交替而组成的层状盐岩体，建立泥岩夹层和盐岩复合体代表单元，根据应变协调原理，从细观力学分析角度通过考虑泥岩夹层弹性性质、盐岩弹性及蠕变力学特性以及两相体积含量建立层状盐岩体宏观各向异性非线性蠕变增量型本构模型，分析层状盐岩在蠕变过程中因保持细观应变协调而产生的应力重分布问题，给出新本构模型ABAQUS 有限元二次开发增量迭代算法实现方法，对一层状盐岩体简单试样算例进行初步计算分析验证，分析计算与试验结果吻合较好，表明该本构模型能反映细观应变协调时层状盐岩体宏观蠕变力学特性，该模型将为层状盐岩体内硐室长期稳定性分析提供理论计算基础。
关键词：岩石力学；层状盐岩；应变协调；细观力学；非线性蠕变；应力重分布
中图分类号：TU 45 文献标识码：A 文章编号：1000–6915(2009)增1–2708–07
NONLINEAR CREEP CONSTITUTIVE MODEL OF BEDDED SALT ROCK
WANG Anming12，YANG Chunhe2，CHEN Jianwen3，HUANG Zhiquan1，LI Yinping2
，
(1. North China Institute of Water Conservancy and Hydroelectric Power，Zhengzhou ，Henan 450011，China ；2. State Key Laboratory of Geomechanics and Geotechnical Engineering，Institute of Rock and Soil Mechanics，Chinese Academy of Sciences，
Wuhan ，Hubei 430071，China ；3. China University of Geosciences，Wuhan ，Hubei 430074，China )
Abstract ：As for alternating layers of mudstone and salt rock，representative volume element(RVE) in the composite is established. Based on strain compatibility principle and micromechanical method，considering elastic behavior of mudstone interbeds，elastic and creep behavior of salt，an anisotropic incremental nonlinear composite creep constitutive model of bedded salt rock is derived；stress redistribution is analyzed during creep deformation of salt. The new model is introduced into finite element software ABAQUS. A simple bedded salt rock sample is modeled as case study to verify the composite model. The result agrees well with the experiment，so it is shown that the constitutive model can reflect mechanical behavior of bedded salt rock considering micro-strain compatibility. The model can be employed to analyze long-term stability of bedded salt rock cavern in the future. Key words：rock mechanics；bedded salt rock；strain compatibility；micromechanical method；nonlinear creep；stress redistribution
1 引 言
盐岩具有良好的密封性、低渗透性、良好的蠕
收稿日期：2008–03–13；修回日期：2008–07–09
变行为和损伤自愈合性，使得盐岩构造已被国际上公认为石油、天然气储存和高放射性废物处置的理想场所。我国层状盐岩的基本赋存特点是：盐岩层数多，单层厚度薄，盐岩体中一般含有众多夹层(如
基金项目：国家重点基础研究发展规划(973)项目(2009CB724600)；国家杰出青年科学基金项目(50725414)；河南省高校杰出科研人才创新工程(2005KYCX015)；华北水利水电学院高层次人才基金项目(200901)；
作者简介：王安明(1975–) ，男，博士，1999年毕业于中南工业大学交通土建专业，现为讲师，主要从事层状盐岩力学特性，地下工程稳定性方面的教学与研究工作。E-mail ：[email protected]
第28卷 增1 王安明，等. 层状盐岩体非线性蠕变本构模型 • 2709 •
泥岩层等) ，在这种特殊地层中建腔进行长期能源储存，国际上没有先例可循，国外储气库的建设经验也不能直接借鉴到我国的储气库的建设中。M. Langer 和S. Hensermann[1]认为储藏库建设的核心是围岩的模拟，这样的模拟应该是尽可能的接近现场的条件包括地质条件，应力条件和本构模型。因此针对我国层状盐岩矿层的基本特点，杨春和和李银平[2]指出建立反映层状盐岩体本构模型及长期强度方程是层状盐岩废弃溶腔油(气) 地下可储性研究中的关键理论问题。
对层状盐岩溶腔进行稳定性分析时，由于盐岩体中含夹层众多，要考虑并模拟其中的每一小层是相当困难的，为避免给建模及计算带来困难，常常采用非均匀介质均匀化等效概化方法，即将厚度很小的夹层与盐岩层概化在一起作为一种等效均匀介质参与分析计算，而且这种概化后等效均匀介质的宏观力学特性与其细观相的力学特性和相对含量是密切相关的。杨春和等
[2，3]
图1 层状盐岩复合体代表单元
Fig.1 Representative volume element for bedded salt rock
(1) 互层层状盐岩介质是由力学性质不同的介质黏结而成，在界面处完美连接，没有滑动、张开或者嵌入；
(2) 细观上夹层和盐岩为各向同性连续均匀介
质，这样考虑下的复合体宏观上为均匀各向异性体(横观各向同性体) 。 2.2 应变协调原理
以σij ，εij 分别表示代表单元宏观平均应力和平
s s ，εij 分别表示盐岩的平均应力和平均均应变；以σij
m m
应变，以σij ，εij 分别表示代表泥岩夹层的平均应
建立的层状盐岩体宏观
平均意义下Cosserat 介质扩展本构模型即体现出这个特点，此模型在层状岩盐溶腔的静力分析时为层状盐岩体变形破坏提供了一条新思路；同理在对层状岩盐容腔进行长期稳定性分析时，也需要对层状盐岩建立等效概化的蠕变本构模型，最新的试验研究表明，层状盐岩的蠕变呈现出明显的非线性性质[4]，因此建立的本构模型尽可能的反应出这种非线性性质。刘 江等[5
，6]
力和平均应变，以V s ，V m 分别表示盐岩和泥岩夹层的体积含量，其中V s +V m = 1，由应变协调原理和多相混和物本构关系的一般原理[7形式之间的关系：
s m
d σij =V s d σij +V m d σij (ij = 11，33，13) (1) s m d εij =d εij =d εij (ij = 11，33，13) (2) s m d εij =V s d εij +V m d εij (ij = 22，32，12) (3) s m d σij =d σij =d σij (ij = 22，32，12) (4)
～9]
可以列出单元体
宏观平均应力–应变与细观上各相应力–应变增量
用黏弹塑性蠕变理论来建立层状(互
层) 盐岩蠕变本构方程，这些本构模型多半为经验型本构模型，没有反映出层状(互层) 盐岩本构模型的宏细观特性，也没有体现出层状盐岩蠕变的非线性性质。
本文用细观力学分析方法通过考虑泥岩夹层弹性性质、盐岩弹性及蠕变力学特性和各自相对体积含量建立层状盐岩体等效非线性蠕变增量型本构模型，最后给出了本构模型的二次开发实现方法并进行了初步的试验验证。
2.3 层状盐岩复合体中各相应力–应变关系
为了求出层状盐岩复合体应力–应变关系，必须先知道各相应力–应变关系，杨春和等[2
，3]
在建立
层状盐岩复合弹性本构模型时假定泥岩夹层和盐岩均为弹性，本文在此基础上建立复合体蠕变本构关系时引入盐岩层的蠕变项。试验研究[5]表明，泥岩夹层和盐岩均具有较为明显的蠕变特性，但泥岩的稳态蠕变率比盐岩层低一个数量级左右，层状盐岩的蠕变性质主要来自于盐岩部分。因此这里假定泥岩夹层仅为弹性变形，忽略泥岩夹层的蠕变，而盐岩层的变形由弹性加蠕变复合而成，这样一种假定大大简化了层状盐岩蠕变性质的分析。
对泥岩夹层有：
2 层状盐岩体增量型复合蠕变本构
模型
2.1 基本假定
选取如图1所示盐岩和泥岩层组成的层状盐岩复合体代表单元(RVE)，由连续介质力学理论推导复合体等效本构关系时采用以下假设：
• 2710 • 岩石力学与工程学报 2009年
m m
d σij =D ijkl d εkl (5)
对盐岩层有：
s s scr d εij =D ijkl (dεkl −d εkl ) (6)
⎛2E µs ⎞2E s 3n −1
s 11q n −1+s q ∆tA 1⎜s 2⎟ (10) 233
41µ1µ−−s s ⎝⎠
s m
d σ22=d σ22=d σ22 (11)
由式(5)得
m d σij µm m m m m
d εij =(dσ11+dσ33) δij (7) −+d σ22
m m 2G E
由式(6)得弹性变形：
s d σij µs s s s s
d εij =(dσ11+dσ33) δij −+d σ22
s s 2G E 式中：D ijkl 为刚度矩阵；G m ，E m 分别为泥岩夹层的剪切模量和弹性模量；G s ，E s 分别为盐岩层的剪切模量和弹性模量；µm ，µs 分别为泥岩夹层和盐岩层的泊松比。对于δij ，当i =j 时，δij = 1；当i ≠j 时，δij = 0。
研究结果[10]表明，盐岩蠕变变形主要来自于其稳态蠕变，盐岩的稳态蠕变率可以表示成
ij =ε
s ij 3
A 1q n
q 2
m
d σ11=
µm E m E m µm
d ε+d ε+d σ22 (12) 113322
1−µm 1−µm 1−µm
µm E m µm E m
d ε+d ε+d σ22 (13) 113322
1−µm 1−µm 1−µm
s s s
d σ13=2G s d ε13−3G s 13∆tAq n =2G s d ε13−
q
s
3G s 13∆tA 1q n (14)
q
m d σ33=
m m d σ13=2G m d ε13=2G m d ε13 (15)
s m
d σij =d σij =d σij (ij = 32，12) (16)
以上各式表明了层状盐岩中细观上盐岩和泥岩分担应力增量与层状岩盐宏观应力应变增量之间的关系，显然随时间步不同，盐岩和泥岩分担应力增量也在不断调整当中，因此发生了应力重分布。
式中：A 1，n 均为材料参数；q 为Mises 应力；s ij 为偏应力。
稳态蠕变为
cr
=d εij
s ij 3
A 1q n ∆t
q 2
s
d σij
4 层状盐岩体增量型蠕变本构模型矩
阵形式
消去岩盐层和泥岩夹层的应力–应变增量之后，可以得到层状盐岩体宏观等效应力–应变关系，写成矩阵形式为 ⎧d ε11⎫⎡M 11⎪d ε⎪⎢M ⎪22⎪⎢21⎪d ε33⎪⎪⎢M 31⎪⎨⎬=⎢d ε⎪12⎪⎢0⎪d ε23⎪⎢0⎪⎪⎢d ε⎢0⎪13⎭⎪⎣⎩
M 12
M 22M 32000
M 13M 23M 33000
000M 4400
0000M 550
0⎤0⎥⎥0⎥⎥⋅ 0⎥0⎥⎥M 66⎦⎥
盐岩总应变为
d ε=d ε+d ε=
s
ij
e ij
cr ij
−(dσ+d σ+ 2G s E s 3s ij s
∆tA 1q n (8) d σ33) δij +2q
s 11
s 22
µs
式中：∆t 为时间增量。
3 层状盐岩复合体中盐岩层和泥岩夹
层分担应力
联立上面的式(1)～(8)并经过推导得到
E s E s µs µs s
d σ11=d +d +d σ22− εε1133
1−µs 21−µs 21−µs
⎛2E s ⎞2E s µs 3n −1n −1∆tA 1⎜+s q s q ⎟ (9) 211233
−−41µ1µs s ⎝⎠
⎧d σ11⎫⎧F 11⎫
⎪d σ⎪⎪F ⎪⎪22⎪⎪22⎪⎪⎪d σ33⎪⎪⎪⎪F 33⎪⎪
=⎨⎬⎨⎬ (17)
⎪d σ12⎪⎪F 12⎪⎪d σ23⎪⎪F 23⎪⎪⎪⎪⎪d σ⎪⎩13⎪⎭⎪⎩F 13⎪⎭
其中，
s
d σ33=
E s µs E s
d ε+d ε3311
1−µs 21−µs 2
s
d σ33=
µs E s µs E s
εεd +d +d σ22−1133
1−µs 21−µs 21−µs
⎛⎞
⎜⎟1⎜11⎟ +M 11=
2⎜s s s s m m m m ⎟++⎜1−µ1−µ1+µs 1+µm ⎟s m ⎝⎠
第28卷 增1 王安明，等. 层状盐岩体非线性蠕变本构模型 • 2711 •
V s µs V µ
+m m
1−µs 1−µm
M 12=
s s +m m
1−µs 1−µm
⎛⎞⎜⎟1⎜11⎟ M 31=−2⎜s s ⎟s s m m
+m m ⎟⎜1−µ+1−µ1+µs 1+µm ⎠s m ⎝
V s µs V µ
+m m
1−µs 1−µm =−
s s +m m
1−µs 1−µm
⎛⎞
⎟1⎜11⎜⎟ M 13=−
s s m m m m ⎟2⎜s s
++⎜1−µ1−µ1+µs 1+µm ⎟s m ⎝⎠⎛1
(q n −1s 11−q n −1s 33) ⎜1+µs 3
+ F 11=E s V s A 1∆t ⎜
E V 4⎜s s
+m m
⎜1+µs 1+µm ⎝
1⎞(q n −1s 11+q n −1s 33) ⎟1−µs
⎟
s s m m ⎟
+⎟1−µs 1−µm ⎠V s µs V µ+m m
1−µs 1−µm
M 21=−
s s +m m
1−µs 1−µm M 22=
M 32
⎛⎞
⎟1⎜11⎟ M 33=⎜+
s s V m E m V m E m ⎟2⎜s s
+⎜1−µ+1−µ11+µm ⎟µ+s m s ⎝⎠⎡1
(q n −1s 11+q n −1s 33) ⎢1−µs 3
F 33=E s V s A 1∆t ⎢− s s m m 4⎢+⎢+µ11+µm s ⎣
1⎤
(q n −1s 11−q n −1s 33) ⎥1+µs
⎥
s s m m ⎥+⎥1−µs 1−µm ⎦
M 44=
V s V m ⎛µs V m µm V m
+−⎜−E s E m ⎝E s E m ⎞
⎟⋅ ⎠
V s V m
+ 2G s 2G m
3s
F 12=V s 12∆tA 1q n
2q M 55=M 66
V s V m
+ 2G s 2G m
1
= 2V s G s +2V m G m
⎡E m ⎛2V s µs 2V m µm ⎞⎤
+⎢⎥⎜⎟
⎢1−µm ⎝1−µs 1−µm ⎠−2µm ⎥ ⎢s s m m 1−µm ⎥
+⎢⎥11µµ−−s m ⎢⎥⎣⎦V s µs V µ+m m
1−µs 1−µm =−
s s +m m
1−µs 1−µm
M 23
3s
F 23=V s 32∆tA 1q n
2q s 13
∆tA 1q n q
F 13=
2V s G s +2V m G m
3V s G s
M 11=M 33
F 22=A 1∆t ⋅ ⎡⎛µs V m µm V m
−⎢⎜E E m ⎢⎝s ⎢⎢⎢⎣
⎤
⎥3n −1
V s q s 22⎥ 2⎥
⎥⎦
⎞3E s E m V s
(q n −1s 11+q n −1s 33) ⎟
⎠2(1−µs )(1−µm ) + V s E s V m E m
+
1−µs 1−µm
M 12=M 21M 13=M 31M 12=M 32M 21=M 23M 44=M 55
式(17)即为层状盐岩复合体宏观增量型应力–应变关系，简写成下面的统一形式：
d ε=M d σ+F (18)
• 2712 • 岩石力学与工程学报 2009年
式中：M 为常数项，与泥岩夹层和盐岩的弹性常数及其各自体积含量有关；F 为蠕变增量项，与泥岩夹层和盐岩的弹性常数、各自体积含量、盐岩层当前应力及时间步长有关。
显然宏观本构模型中含有其细观成分盐岩层当前应力及时间步长，从式(9)～(16)可以看出，随蠕变时间变化，盐岩层分担应力在不断调整之中，可以看出此本构模型能够反映层状盐岩在蠕变过程中因保持细观应变协调而产生的应力重分布问题，体现了蠕变的非线性性质
[11]
(4) 第四步：ABAQUS 的平衡迭代开始于新时刻t n +1，ABAQUS 对于最大迭代次数的默认设置为
9，最小容差TOL 为5×106，当通过小于9次的迭
－
代后，如果误差小于TOL ，则为收敛，n = n + 1时进行下一步增量计算；如果误差大于TOL ，ABAQUS 将会缩减该步时间增量∆t 直到计算结果收敛为止。
6 算例验证与讨论
为验证新本构模型的有效性，考虑由泥岩夹层和盐岩组成的层状岩盐复合体，在ABAQUS 软件中建立圆柱体复合体模型，采用C3D8I 单元，底边施加固定边界条件，按照上述二次开发算法实现方法将本文建立的非线性蠕变本构方程通过UMAT 子程序植入ABAQUS 有限元软件。对某层状盐岩试样的三轴压缩蠕变试验[5](围压为5 MPa，差应力为10 MPa)进行初步的计算模拟，盐岩和泥岩夹层力学参数和体积含量如下面给出：取盐岩弹性模量
；另外对于互层盐岩
体，其宏观上呈现横观各向同性性质，显然上述建立的本构模型也体现了这样一种各向异性性质；如果不考虑∆t ，即不考虑蠕变项的话，上述本构模型退化为增量形式的弹性本构模型，这时与李银平和杨春和建立的层状盐岩体三维Cosserat 介质扩展本构模型(去掉偶应力项) 本质上是一致的，这也反映了本文本构模型的合理性。
[3]
5 新本构模型增量迭代算法及其在
ABAQUS 中的实现
在给定远场外载荷的情况下，层状盐岩体宏观应力应变关系可以通过上述建立的增量型应力应变关系逐次进行增量跌代方法得到，参考叶 峥等[12
，13]
E = 5 150 MPa，泊松比µs = 0.31，夹层弹性模量E =
21 780 MPa，泊松比µm = 0.23，夹层含量为44.1%，其盐岩层稳态蠕变本构关系符合Norton 蠕变模型，参数[5]取为A 1=1.12×108/h，n = 3.8，图2给出了
－
理论预测结果与试验结果对比曲线。
从图2中可以看出，数值预测与试验结果比较吻合，这说明本文非线性蠕变本构模型可较好地模拟层状岩盐的蠕变变形。当然蠕变计算值与试验值相比也有一定误差，其原因是多方面的[2
，3]
的方法，将上述本构模型和迭代算法编写成子程序
UMAT 嵌入大型有限元软件ABAQUS 进行二次开发，下面简述其主要步骤。将节4建立的本构模型式(18)简写成下面的形式：
1d σ=(dε−F ) (19)
M
，首先
是岩体材料力学特性的离散性，参数选取难度较大，然后是理论和试验的模型以及边界条件是有区别的，对于层状岩盐试样，现场钻孔时为圆柱形试样，而本文理论模型建立时以长方体形状的代表单元为起点进行推导的，复合体为横观各向同性性质，表明试验试样与理论模型存在不一致性，另对于三轴试验，虽然侧向压力相同，由于各层力学特性不同，其侧向变形是不相同的，这表明试验条件并不完全满足理论分析中的假设(1)，而理论分析时假设层面方向是无限展开的，故而各层侧向位移是相同的，但其侧向压力是不相同的。当然，刘 江[5]只是做了一个初步的试验，试验岩样是现场钻取的含小夹层的盐岩岩芯，百分比含量不易较大范围变化，系统地开展层状盐岩蠕变试验工作量大等，当前只验证的是层状盐岩其垂直于层面方向的蠕变力学性质，
(1) 第一步：由层状盐岩初始应力场，代入式(9)～(16)求得没有蠕变时盐岩和泥岩的初始分担应力。
(2) 第二步：在平衡时刻t n ，ABAQUS 提供总应变增量∆εtotal (t n ) ，总应变εtotal (t n ) 和时间增量∆t ，由盐岩和泥岩的当前应力(在前一步已经计算出来) ，进而可以求出时间步长∆t 内蠕变应变增量F 项的值，代入上述本构关系式(19)之后可以求出∆t 时段内总应力增量d σ的值。
(3) 第三步：由总应变增量、总应力增量d σ及当前分析步的当前应力代入式(9)～(16)，得到当前时步盐岩和泥岩的应力增量d σs ，d σm ，与盐岩和泥岩当前应力进行叠加后可以得到盐岩和泥岩下一
时步初始应力σs (t n +1) ，从而进行应力更新。
第28卷 增1 王安明，等. 层状盐岩体非线性蠕变本构模型 • 2713 •
累，即把非线性蠕变计算进行线性化处理[11
，15]
，运
用非线性问题增量迭代算法将本文本构模型写成子程序与ABAQUS 软件连接，为层状盐岩体容腔长
蠕变/10
－4
期稳定性分析奠定了基础。
15
1 3 5 7时间/(104 s)
91113
参考文献(References)：
[1] LANGER M ，HEUSERMANN S. Geomechanical stability and
integrity of waste disposal mines in salt structures[J]. Engineering Geology ，2001，21(2/3)：155–161. [2]
杨春和，李银平. 互层盐岩体的Cosserat 介质扩展本构模型[J]. 岩石力学与工程学报，2005，24(23)：4 226–4 232.(YANG Chunhe，LI Yinping. Expanded Cosserat medium constitutive model for laminated salt rock[J]. Chinese Journal of Rock Mechanics and Engineering ，2005，24(23)：4 226–4 232.(in Chinese)) [3]
李银平，杨春和. 层状盐岩体的三维Cosserat 介质扩展本构模型[J]. 岩土力学，2006，27(4)：509–513.(LI Yinping，YANG Chunhe. Three-dimensional expanded Cosserat medium constitutive model for laminated salt rock[J]. Rock and Soil Mechanics，2006，27(4)：509–513.(in Chinese)) [4]
郤保平，赵阳升，赵金昌，等. 层状盐岩温度应力耦合作用蠕变特性研究[J]. 岩石力学与工程学报，2008，27(1)：90–96.(XI Baoping ，ZHAO Yangsheng，ZHAO Jinchang. Study on coupled thermo-mechanical creep properties of bedded rock salt[J]. Chinese Journal of Rock Mechanics and Engineering，2008，27(1)：90–96.(in Chinese)) [5]
刘 江. 层状盐岩力学特性试验研究及理论分析[硕士学位论文][D]. 武汉：中国科学院武汉岩土力学研究所，2006.(LIU Jiang. Experimental investigation and theoretic analysis on the mechanical properties of layered rock salt[M. S. Thesis][D]. Wuhan：Institute of Rock and Soil Mechanics，Chinese Academy of Sciences，2006.(in Chinese)) [6]
郤保平. 含高盐份泥岩夹层的盐岩蠕变特性及油气储库稳定性研究[硕士学位论文][D]. 太原：太原理工大学，2006.(XI Baoping. Study on creep property of rock salt with salt-mudstone interlayer and stability of oil and gas storage cavern[M. S. Thesis][D]. Taiyuan：Taiyuan University of Technology，2006.(in Chinese)) [7]
范镜泓，高芝晖，刘天模. 具层状细观结构的PD3珠光体钢的细/宏观力学分析[J]. 重庆大学学报，1999，22(5)：37–46.(FAN Jinghong ，GAO Zhihui，LIU Tianmo. The micro/macro mechanical
图2 层状岩盐蠕变理论预测值与试验值对比曲线 Fig.2 Comparison between simulated values and experimental
results based on creep theory of bedded salt rock
别的方向的试验验证尚需钻取更多与理论模型比较一致的试样来进行进一步的试验。
此外，对于本文本构模型，当t = 0时，层状盐岩复合体呈弹性响应，盐岩层在任意时间步t +∆t 发生蠕变，但对于深部地下层状盐岩容腔加内压后，其围岩力学响应盐岩层往往是弹塑性的，但由于盐岩弹塑性变形相对其蠕变变形来说很小，况且在容腔设计中，与盐岩的高应变率弹塑性响应相比，盐岩的黏弹性响应是最为重要的[14]，从这一点来说本文模型因没有考虑塑性变形导致的误差影响不大，但是要建立能考虑塑性变形的复合本构模型仍需进一步深入研究。
7 结 论
本文针对由不同岩层交替而组成的层状盐岩体，建立泥岩夹层和盐岩复合体代表单元，根据应变协调原理，用细观力学分析方法通过考虑泥岩夹层弹性性质、盐岩弹性及蠕变力学特性以及两者体积含量建立层状盐岩体宏观非线性蠕变增量型本构模型，将细观上不均匀的两相介质在宏观上进行了均匀化处理，分析层状盐岩在蠕变过程中因保持细观位移协调而产生的应力重分布问题。该模型既能考虑由于泥岩夹层的存在而导致的层状盐岩构造体宏观平均意义下力学特性的各向异性，又能很好地反映层状盐岩蠕变的非线性性质，与试验研究结论非常一致[4]。
在数值计算中把蠕变经历的时间分成有限间隔，假定在每一个时间步内应力不变，可按某一时刻初的应力水平进行蠕变，这样就可以把变应力情况求蠕变量问题看作是求逐段常应力蠕变量的积
• 2714 • 岩石力学与工程学报 2009年
analysis of PD3 pearlitic steel with layered microstructures[J]. Journal of Chongqing University(Natural Science Edition)，1999，22(5)：37–46.(in Chinese))
[8] HAJ-ALI R M ，MULIANA A H. A multi-scale constitutive
formulation for the nonlinear viscoelastic analysis of laminated composite materials and structures[J]，International Journal of Solids and Structures，2004，41(13)：3 461–3 490. [9]
HAJ-ALI R M，MULIANA A H. A micromechanical constitutive framework for the nonlinear viscoelastic behavior of pultruded composite materials[J]. International Journal of Solids and Structures ，2003，40(5)：1 037–1 057.
[10] 陈 锋. 盐岩力学特性及其在储气库建设中的应用研究[博士学
位论文][D]. 武汉：中国科学院武汉岩土力学研究所，2006.(CHEN Feng. The mechanical behavior of salt rock and its application to the natural gas storage construction[Ph. D. Thesis][D]. Wuhan：Institute of Rock and Soil Mechanics，Chinese Academy of Sciences，2006.(in Chinese))
[11] 彭 凡，彭 献，傅依铭. 纤维增强复合材料非线性蠕变流动特
性[J]. 湖南大学学报(自然科学版) ，1999，26(5)：6–9.(PENG Fan，
PENG Xian，FU Yiming. Behavioral analyses of creep flow for fiber composites[J]. Journal of Hunan University(Natural Sciences)，1999，22(5)：6–9.(in Chinese))
[12] 叶 峥. 基于微观滑移的蠕变本构模型及应用研究[硕士学位论
文][D]. 北京：北京航空航天大学，2007.(YE Zheng. Slip theory based modeling of creep deformation and application in Ni-based single crystal super alloys[M. S. Thesis][D]. Beijing：Beijing University of Aeronautics and Astronautics，2007.(in Chinese)) [13] 罗海安，翟 军. 兼具基体蠕变和轻度界面蠕变的复合材料的本
构特性[J]. 力学学报，1996，28(5)：550–557.(LUO Haian，ZHAI Jun. Constitutive behavior of the composites with both matrix creep and slight interface creep[J]. Acta Mechanica Sinica，1996，28(5)：550–557.(in Chinese))
[14] LIANG W ，YANG C，ZHAO Y，et al. Experimental investigation of
mechanical properties of bedded salt rock[J]. International Journal of Rock Mechanics and Mining Sciences，2007，44(3)：400–411. [15] 穆霞英. 蠕变力学[M]. 西安：西安交通大学出版社，1990.(MU
Xiaying. Creep mechanics[M]. Xi′an ：Xi ′an Jiaotong University Press ，1990.(in Chinese))