岩石圈动力学研究是一个主要科研课题.对岩石力学性质的认识，提供岩石圈动力学研究的基础，前人从不同的角度出发研究了岩石的特性.比如描述岩石受力变形和撤力后变形恢复情况的弹塑性；描述岩石破裂前变形程度和变形能积累的脆韧性；岩石在短期载荷下表现出固体的弹性，而在长期载荷下，岩石会表现出一些类似于流体的性质^[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)中绿色区域为长石，蓝色区域为云母，黄色区域为石英.

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

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，证明了程序的可靠性.

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.

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

将使用数值模拟计算得到的流变参数结果，与通过高温高压岩石实验获得的参数结果进行对比，应力指数结果很接近，激活能在数量级上也相近.其中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.

从计算得到的等效黏滞系数可以看出，加载外力与矿物条带之间角度的不同，等效黏滞系数不同，其中0°和90°时等效粘滞系数值较大，而在30°时最小.这与前人所做的实验结果较一致.也印证了岩石强度受矿物层状分布的影响^[14].

3 结论

前人已经说明在矿物颗粒随机分配的情况下，可以在矿物颗粒流变性质已知的条件下，计算多矿物岩石整体的流变性质.本文通过使用有限元法，在已知矿物流变参数值已知，且近似认为矿物为各向同性的条件下，根据具有面理的花岗质糜棱岩薄片矿物沿面理层分布特征，建立岩石模型，对岩石整体的流变参数进行计算.得出的结果能与高温高压岩石实验所做的类似实验结果相互验证.说明在岩石存在片理的情况下，使用数值计算来获取岩石流变参数值的方法仍然是可行的.

通过不同方向载荷下的计算结果，可以看出岩石的流变强度存在方向性，与其矿物沿片理成层分布的方式有关.当外加载荷与矿物条带垂直或平行时，所得的等效黏滞系数值较高，也可以理解为岩石的流变强度较高.当载荷与岩石面理呈30°角时，流变强度最小.这个结果与岩石实验得到的结果一致.

当然数值计算，还需要考虑其他的实际因素，如矿物颗粒大小、含水性等，许多矿物颗粒自身也具有各向异性.各向异性矿物且空间分布具有一定方向时，如何计算它们集合形成的岩石流变性质，是下一步需要进行的工作.