蒋龙腾, 刘贵, 石耀霖
中国科学院计算地球动力学重点实验室, 北京 100049
2016年01月18日 收稿; 2016年03月15日 收修改稿
基金项目: 国家自然科学基金重大项目（41086965）资助
通信作者: E-mail:shiyl@ucas.ac.cn

摘要: 高温高压下岩石表现出流变性质。流变参数是数值模拟地球动力过程必需的重要参数。获得流变参数需要借助于室内的高温高压蠕变实验以及野外的地球物理资料观测的反演。对单矿物的高温高压实验比较多，多矿物组成的岩石的实验则较难较少。本文探讨在单矿物流变参数已知的条件下，能否运用有限元法，计算多矿物岩石整体的流变参数。研究表明，计算结果与高温高压实验获得的参数可以对应。即使矿物力学性质可以近似为各向同性时，如果矿物排列有优势方向（例如存在面理的花岗质糜棱岩），计算表明岩石整体会显示各向异性。当加载外力与矿物面理方向垂直或平行时，岩石的等效黏滞系数较大；当所成角度为30°左右时，等效黏滞系数最小。
关键词: 有限元方法流变参数数值模拟
Calculations of rheological properties of multi-mineral foliated rocks using finite element method based on the rheological properties of minerals
JIANG Longteng, LIU Gui, SHI Yaolin
Key Laboratory of Computational Geodynamics, Chinese Academy of Sciences, Beijing 100049, China


Abstract: Rocks show rheological properties under high temperature and high pressure.Rheological parameter is an important parameter for numerical simulation of the dynamic process of the earth.The rheological parameters can be obtained by means of the high temperature and high pressure creep experiments in the laboratory and by the inversion of the field geophysical data.There have been some high temperature and high pressure experiments on single minerals.However, it is difficult to carry out multi-mineral experiments.In this study, we discuss the possibility of using the finite element method to calculate the rheological parameters of multi-mineral rocks under the condition that the rheological parameters of minerals are known.The results show that the calculated results match the parameters obtained in the high temperature and high pressure experiments.The results show that the rock as a whole will show anisotropy if the mineral arrangement has the advantage of direction (such as the existence of the granitic mylonite foliation), even though mineral mechanical properties are approximatly isotropic.When loading force is perpendicular or parallel to mineral foliation, equivalent viscous coefficient of rock is larger; and when the angle is 30 degrees, equivalent viscous coefficient is minimized.
Key words: finite element methodrheological parameternumerical simulation
岩石圈动力学研究是一个主要科研课题.对岩石力学性质的认识，提供岩石圈动力学研究的基础，前人从不同的角度出发研究了岩石的特性.比如描述岩石受力变形和撤力后变形恢复情况的弹塑性；描述岩石破裂前变形程度和变形能积累的脆韧性；岩石在短期载荷下表现出固体的弹性，而在长期载荷下，岩石会表现出一些类似于流体的性质^[1].所以，对于岩石流变学方面，也就是说在物体受力后，应变随时间变化的研究，是十分必要的.对于岩石圈流变结构的探讨将有助于对岩石圈动力学过程的理解.
岩石流变的表现主要有蠕变、应力松弛、弹性后效等.岩石蠕变是岩石在受到外力不变的情况下，应变随时间延长而增加的现象，包含暂态蠕变和稳态蠕变.本文主要讨论稳态蠕变.在稳态蠕变阶段，应变一般随时间均匀变化，也就是说应变速率为一个定值，应力和应变速率的比值为等效黏滞系数.
在建立以黏弹性岩石圈为基础的结构模型时，需要清楚不同层的流变特征和黏滞系数.到目前为止，采用的方法有以下几种.第1是，通过矿物和岩石的高温高压实验.样品在三轴高温高压实验系统上进行.在一定温度和围压下，以等应变速率加载，可测得相应应力.分别改变应变速率和温度的实验后，可得出相应岩石的幂律方程^[2].第2是，利用GPS等空间大地测量手段，可以观测到岩石圈受到加载或卸载后的地表位移变化，如冰盖融化后的岩石圈响应^[3]，或者大地震后的震后松弛变形^[4]等.通过测量数据和建立模型，能够反演出岩石圈的流变性质.第3是，我们将在这里使用的，在矿物黏滞系数已知时，通过数值计算来估算多矿物岩石整体的黏滞系数.
关于数值计算的方法，前人将其应用在很多种岩石物性参数的计算中，如岩石的热导率和电导率，取得了很好的效果^[5]，相比较实验，其花费时间较少，同时无需计较样品的尺寸.关于流变参数的有限元计算前人做了许多工作.Tullis等^[6]计算两种矿物不同比例混合所得岩石的流变参数，同时给出在已知矿物参数和所占体积比的情况下，估计岩石强度的方法.崔晓佳等^[7]计算三维情况下锗橄榄石和锗尖晶石组成的双矿物岩石的流变参数，使用Monte Carlo法得到不同比例矿物随机分布的结果.Jing等^[8]探讨火星上冰岩石混合物的流变性质.他们一般都是假定矿物颗粒随机分布的前提下进行计算，然而许多岩石具有特定组构，矿物颗粒沿一定优势方位排列，这对岩石整体流变性质会有什么影响呢？本文所做的计算是依据实验中具有片理的变质岩石样品来建立模型，同时考虑边界条件和加载方向的不同，讨论即使矿物可假定为各向同性时，矿物颗粒定向排列对岩石整体流变参数的影响，发现岩石整体会呈现各向异性，并与高温高压实验结果吻合.
1 方法原理1.1 基本原理计算中用到方程分别为：描述力之间关系的平衡方程

${{\sigma }_{ij,j}}+{{f}_{i}}=0$

(1)

描述应变速率和变形速度之间关系的几何方程

${{{\dot{\varepsilon }}}_{ij}}=({{u}_{i,j}}+{{u}_{j,i}})/2$

(2)

以及描述差应力和应变速率之间幂函数关系的本构方程

$\bar{\sigma }=\eta {{{\bar{\dot{\varepsilon }}}}^{\frac{1}{n}}}$

(3)

其中$\bar{\sigma }=\sqrt{\frac{3}{2}{{S}_{ij}}{{S}_{ij}}},{{S}_{ij}}$为偏应力张量，S_ij=σ_ij-pδ_ij，$\bar{\dot{\varepsilon }}=\sqrt{\frac{2}{3}{{{\dot{\varepsilon }}}_{ij}}{{{\dot{\varepsilon }}}_{ij}}}$.
对于微分方程和边界条件，可以建立泛函

$\begin{align} & \prod \left( u \right)={{\int }_{V}}\frac{a}{m+1}{{{\bar{\dot{\varepsilon }}}}^{^{m+1}}}dV-{{\int }_{V}}~fudV-{{\int }_{\Gamma }}Tud\Gamma - \\ & {{\int }_{V}}\left( \Delta \cdot u \right)pdV, \\ \end{align}$

(4)

利用变分原理δΠ(u)=0，并代入

$\begin{align} & \bar{\dot{\varepsilon }}=\sqrt{\frac{2}{3}{{{\dot{\varepsilon }}}_{ij}}{{{\dot{\varepsilon }}}_{ij}}}, \\ & \delta \prod \left( u \right)={{\int }_{V}}\frac{a}{2}{{\left( \frac{2}{3}{{{\dot{\varepsilon }}}_{ij}}{{{\dot{\varepsilon }}}_{ij}} \right)}^{\frac{m-1}{2}}}\frac{4}{3}{{{\dot{\varepsilon }}}_{ij}}\delta {{{\dot{\varepsilon }}}_{ij}}dv-{{\int }_{V}}f\delta udV- \\ & {{\int }_{\Gamma }}Tud\Gamma -{{\int }_{V}}\left( \Delta \cdot u \right)pdV- \\ & {{\int }_{V}}\left( \Delta \cdot \delta u \right)pdV=0. \\ \end{align}$

(5)

定义F(u)=∫_V${a \over 2}{\left( {{2 \over 3}{{\dot \varepsilon }_{ij}}{{\dot \varepsilon }_{ij}}} \right)^{{{m - 1} \over 2}}}{4 \over 3}{\dot \varepsilon _{ij}}\delta {\dot \varepsilon _{ij}}dv - {\smallint _V}f\delta udV - {\smallint _\Gamma }Tud\Gamma - {\smallint _V}\left( {\Delta \cdot u} \right)pdV - {\smallint _V}\left( {\Delta \cdot \delta u} \right)pdV = 0$，由牛顿迭代法，F(u+Δu)=F(u)+F^′(u)Δu=0，可得F^′(u)(u+Δu)=F^′(u)u-F(u)，分别代入可建立离散化的有限元计算方程.边界条件添加的是速率条件.
1.2 模型建立模型建立依据具有面理的花岗质糜棱岩实验样品^[2]在电镜下的照片，区域面积为2mm×2mm，建立二维模型，并划分出不同矿物的区域.整个模型节点总数为6826，单元总数为13382.如图 1(a)，图 1(b).图 1(b)中绿色区域为长石，蓝色区域为云母，黄色区域为石英.
Fig. 1

	Download: JPG larger image
图 1 (a)电镜下岩石样品照片 (b)岩石模型 Fig. 1 (a) Rock sample photo under electron microscope and (b) rock model

矿物材料参数的选取依据高温高压实验测得的结果^[9-11]，具体数据如表 1.
Table 1

表 1 3种矿物的流变参数Table 1 Rheological parameters of three kinds of minerals 矿物 应力指数n 激活能E/(kJ/mol) 物质常数A/(MPa-n·s-1) 长石 3 332 2.51e3 石英 2.97 242 3.55e-3 云母 1.13 270 1.70e7

表 1 3种矿物的流变参数Table 1 Rheological parameters of three kinds of minerals

1.3 流变参数计算过程我们通过已知的矿物参数，计算出岩石整体的流变参数，即幂律方程中的应力指数和激活能.
对幂律方程等式两边取自然对数，可得

$ln\dot{\varepsilon }=nln\sigma +lnA-\frac{E}{R}\frac{1}{T}.$

(6)

为了得到应力常数，我们进行一组计算，在一定的温度下，改变边界速率条件，可以计算得到一组差应力和应变速率关系，将其取自然对数并进行线性拟合，所得斜率即为应力指数n.
等式又可以变换为

$ln\sigma =\frac{E}{nR}\frac{1}{T}+\frac{1}{n}ln\dot{\varepsilon }-\frac{1}{n}lnA.$

(7)

在同样的应变速率下，改变温度，计算得到一组差应力，将差应力取自然对数和温度倒数进行线性拟合，所得斜率再乘以nR即为激活能E.
1.4 程序可靠性验证在正式计算之前，为验证程序的可靠性，我们先对单一矿物模型进行计算.使用的是长石的流变参数.计算结果与所添加参数值相对误差在1%左右，如表 2，证明了程序的可靠性.
Table 2

表 2 流变参数值的相对误差Table 2 Relative errors in the rheological parameter values 长石流变参数 添加参数值 计算结果值 相对误差/% n 3 2.948 -1.7 E/(kJ/mol) 332 328 1.1 A/(MPa-n·s-1) 2.51e3 2.49e3 0.8

表 2 流变参数值的相对误差Table 2 Relative errors in the rheological parameter values

2 计算结果我们计算分2组情况.第1组情况，有4个模型.每个模型之间，岩石区域的速率边界条件相同，温度不同.速率边界条件是1×10^-7m/s，温度分别是700，800，850和900℃.第2组情况，5个模型，每个模型之间，温度相同，岩石区域的速率边界条件不同.温度为900℃，速率边界条件分别为1×10^-7，7.5×10^-8，5×10^-8，2.5×10^-8和1×10^-8m/s.
通过每组的结果得到岩石整体的差应力和应变速率.按照方法原理中提到的步骤，将第1组的差应力(单位为Pa)和应变速率结果取自然对数并进行线性拟合，斜率即为所求应力指数n.对第2组结果的差应力取自然对数并与温度的倒数进行线性拟合，得到斜率，可以计算出激活能E.
之后旋转模型区域，分别是原方向的30°，45°，60°和90°，每种模型按照前述2组情况计算，可以得到不同加力方向情况下的应力指数和激活能.计算结果如图 2,图 3，其中差应力单位为Pa.
Fig. 2

	Download: JPG larger image
图 2 应变速率与应力的关系 Fig. 2 Relationship between strain rate and stress

Fig. 3

	Download: JPG larger image
图 3 温度倒数与应力的关系 Fig. 3 Relationship between the reciprocal of temperature and the stress

按前面方法原理一节中谈到的步骤，进行计算，可以得到每个角度下的流变参数，结果如表 3.
Table 3

表 3 岩石流变参数Table 3 Rock rheological parameters 角度/(°) 应力指数n 激活能E/(kJ/mol) 物质常数A/(MPa-n·s-1) 0 2.83 327.7 8.87e3 30 2.47 293.9 2.45 e3 45 2.49 297.8 3.37 e3 60 2.60 307.2 4.70 e3 90 2.72 319.4 7.54 e3

表 3 岩石流变参数Table 3 Rock rheological parameters

将使用数值模拟计算得到的流变参数结果，与通过高温高压岩石实验获得的参数结果进行对比，应力指数结果很接近，激活能在数量级上也相近.其中0°和90°的实验结果n=3和激活能E=(354±54)kJ/mol ^[2]，30°、45°和60°的实验结果n=3±0.2^[12]，这些实验结果与计算模拟结果基本吻合.同时，从这几组计算结果的对比可以看出：当温度一定时，应力随着加载应变速率的增加而增加；当加载速率一定时，应力随着温度的增加而减小.
在已知流变参数的情况下，可以计算出温度为973K和应变速率为1×10^-15m/s时的等效黏滞系数^[13]，计算公式为

${{\eta }_{eff}}=\frac{{{{\dot{\varepsilon }}}^{\frac{1-n}{n}}}}{2{{A}^{\frac{1}{n}}}}{{e}^{\frac{E}{2nRT}}},$

(8)

结果如表 4.
Table 4

表 4 等效黏滞系数Table 4 Equivalent viscous coefficient 角度/(°) 等效黏滞系数/(Pa·s) 0 1.66e20 30 4.39e19 45 4.79e19 60 7.22e19 90 1.16e20

表 4 等效黏滞系数Table 4 Equivalent viscous coefficient

从计算得到的等效黏滞系数可以看出，加载外力与矿物条带之间角度的不同，等效黏滞系数不同，其中0°和90°时等效粘滞系数值较大，而在30°时最小.这与前人所做的实验结果较一致.也印证了岩石强度受矿物层状分布的影响^[14].
3 结论前人已经说明在矿物颗粒随机分配的情况下，可以在矿物颗粒流变性质已知的条件下，计算多矿物岩石整体的流变性质.本文通过使用有限元法，在已知矿物流变参数值已知，且近似认为矿物为各向同性的条件下，根据具有面理的花岗质糜棱岩薄片矿物沿面理层分布特征，建立岩石模型，对岩石整体的流变参数进行计算.得出的结果能与高温高压岩石实验所做的类似实验结果相互验证.说明在岩石存在片理的情况下，使用数值计算来获取岩石流变参数值的方法仍然是可行的.
通过不同方向载荷下的计算结果，可以看出岩石的流变强度存在方向性，与其矿物沿片理成层分布的方式有关.当外加载荷与矿物条带垂直或平行时，所得的等效黏滞系数值较高，也可以理解为岩石的流变强度较高.当载荷与岩石面理呈30°角时，流变强度最小.这个结果与岩石实验得到的结果一致.
当然数值计算，还需要考虑其他的实际因素，如矿物颗粒大小、含水性等，许多矿物颗粒自身也具有各向异性.各向异性矿物且空间分布具有一定方向时，如何计算它们集合形成的岩石流变性质，是下一步需要进行的工作.
参考文献

[1]	石耀霖, 张贝, 张斯奇. 地震数值预报[J].物理, 2013, 42(4):237–255.
[2]	Liu G, Zhou Y, He C h, et al. An experimental study of effect of pre-existing fabric on deformation of foliated mylonite at high temperature and pressure[J].Geological Journal, 2014, 51(1):92–112.
[3]	Burgmann R, Dresen G. Rheology of the lower crust and upper mantle: evidence from rock mechanics, geodesy, and field observations[J].Annual Review of Earth and Planetary Sciences, 2008, 36:531–567.DOI:10.1146/annurev.earth.36.031207.124326
[4]	张晁军, 曹建玲, 石耀霖. 从震后形变探讨青藏高原下地壳黏滞系数[J].中国科学 D辑:地球科学, 2008, 38(10):1250–1257.
[5]	刘善琪, 李永兵, 田会全, 等. 含湿孔隙岩石有效热导率的数值分析[J].地球物理学报, 2012, 55(12):4239–4248.
[6]	Tullis T, Horowitz F G, Tullis J. Flow laws of polyphase aggregates from end-member flow laws[J].Journal of Geophysical Research, 1991, 96(B5):8081–8096.DOI:10.1029/90JB02491
[7]	崔晓佳, 张怀, 石耀霖. 岩石宏观流变性质的数值模拟[J].岩石学报, 2008, 24(6):1417–1424.
[8]	Jing H, Zhang H, Li H, et al. Parallel numerical analysis on the rheology of the martian ice-rock mixture[J].Journal of Earth Science, 2011, 22(2):176–181.DOI:10.1007/s12583-011-0170-0
[9]	Offerhaus L J, Wirth R, Dresen G. High-temperature creep of polycrystalline albite[J].Deformation Mechanisms, Rheology and Tectonics, 2001, 124:107–131.
[10]	Rutter E H, Brodie K H. Experimental intracrystalline plastic flow in hot-pressed synthetic quartzite prepared from Brazilian quartz crystals[J].Journal of Structural Geology, 2004, 26(2):259–270.DOI:10.1016/S0191-8141(03)00096-8
[11]	Mariani E, Brodie K H, Rutter E H. Experimental deformation of muscovite shear zones at high temperatures under hydrothermal conditions and the strength of phyllosilicate-bearing faults in nature[J].Journal of Structural Geology, 2006, 28:1569–1587.DOI:10.1016/j.jsg.2006.06.009
[12]	刘贵, 石耀霖, 周永胜. 各向异性花岗岩的强度和变形机制实验研究[J].吉林大学学报:地球科学版, 2015, 45(7):1506–1514.
[13]	石耀霖, 曹建玲. 中国大陆岩石圈等效粘滞系数的计算和讨论[J].地学前缘, 2008, 15(3):82–95.
[14]	刘贵, 周永胜, 石耀霖. 先存组构对各向异性岩石流变强度的影响[J].地震地质, 2014, 36(3):918–928.