流变基质中刚性椭球夹杂的运动演化机制

马连华, 李雅娆, 谭如峰

(1.河北大学质量技术监督学院，保定 071002；2.河北大学建筑工程学院，保定 071002；3.河北大学人事处，保定 071002)

弹性或黏性基体中椭球夹杂问题的研究已经有很长的历史[1-3]。针对位于地球岩石圈中的高应变岩石流变场分配问题，中外学者采用Eshelby细观力学方法[3-4]进行了充分的研究。Jeffery[5]首先推导了远场均匀流动的各向同性牛顿黏性材料中嵌入刚性椭球的角速度方程。Bilby等[6]将Eshelby细观力学理论扩展应用到椭球夹杂和宏观区域构造具有不同黏性的流变体情况。Freeman[7]基于Jeffery的理论对各种黏性流的刚性椭球夹杂的运动进行了数值研究，并且将Eshelby方程的数值解与Jeffery方程的数值解进行了比较[8]。Ježek[9]开发了FORTRAN程序，用来计算给定流变场中单个刚性椭球夹杂的旋转路径。Jiang[10-11]基于Eshelby理论并利用Mathcad计算了单个以及多个刚性和可变形椭球夹杂的演化进程。Jiang[12-13]通过引入近似算法来求解夹杂运动微分方程，提高了计算效率。向必伟等[14]从固体连续变形机制入手，重点介绍基于Eshelby理论的多尺度数值模拟思路和方法，探讨流变场分配问题，揭示了嵌入基质中的椭球体内分布的流变场主要取决于椭球体与基质间的相对流变强度，椭球体的相对流变强度越低，其变形越接近于简单剪切，且有限应变积累越快。同时还揭示，不同流变强度的椭球体内模拟拉伸线理和面理产状的总体格局反映基质流变场特征。刘稳航等[15]介绍了韧性剪切带中椭球状标志体变形的最新研究进展，探讨了基于Jeffery理论和Eshelby理论的数值模拟思路和方法，并利用Mathcad软件模拟了给定条件下的椭球夹杂的运动轨迹、变形特征。Qu等[16]基于MATLAB模拟了流变岩石中椭球夹杂的运动规律，但并未考虑椭球夹杂在有限体分比下的运动演化机制。高丽敏等[17]基于Eshelby理论,通过数值模拟,将经过两期变形叠加产生的模拟组构产状格局与只经历后期变形产生的模拟组构产状格局进行对比,揭示叠加变形组构的产状分布特征。林少泽等[18]利用适用于黏性材料的Eshelby模型，计算了椭球体内部流变场最大主应变轴和最小主应变轴方位并进行赤平投影，用以模拟线理和面理极点产状。

1 均匀介质的连续变形机制

均匀介质的流变学在岩石流变学方面的应用开始于简单的二维变形的研究。Jiang等[19]提出了一般的三维变形模式，推导出了速度梯度张量和位移梯度张量。

根据Jeffery[5]的定义可以得到，a1、a2、a3为椭球状局部地质体3个应变主轴的半轴，3个半轴应满足：a1≥a2≥a3。为了方便进行数值模拟，需要对内部和外部建立两套不同的坐标系，其中内部的参考坐标系R′是始终与椭球夹杂的3个主轴a1、a2、a3保持平行，长度用相对应的单位基本矢量η1、η2、η3表示；而外部坐标系R使用x1、x2、x3表示，使用相对应的单位基本矢量e1、e2、e3进行表示。

刚性椭球状夹杂在基体中的旋转相当于其内部坐标系R′和外部坐标系R不停发生变化，同时每个旋转位置都可以由外部坐标系R中的单位矢量ei来表示，另外也可以由内部坐标R′的单位矢量ηi来表示。在x1x2x3的坐标系统中，xi轴和x′i轴之间存在的关系为[15]

x′=Qx，x=QTx′

(1)

式(1)中：Q是从xi坐标系到xi′坐标系的过渡矩阵；QT是Q的转置矩阵。同时，也可以定义为

e′i=Qei，ei=QTe′i

(2)

均匀介质的流动场可以使用欧拉速度梯度张量可以表示为

(3)

(4)

(5)

速度梯度张量L可以分解拉伸张量D和涡度张量W[18]：

L=D+W

(6)

式(6)中：

(7)

(8)

根据坐标转换，有：

D′=QDQT

(9)

(10)

只要给定流变场的速度梯度张量L，就可以求解任一时刻流变场的运动学涡度，并可以求解经历任意时间段的变形后的有限应变量大小、主应变大小以及应变主轴的方位。

2 等效远场法

为了应用经典的、适用于无限大基体的Eshelby等效夹杂理论，将含椭球夹杂的代表体积单元(RVE)置于无限大的基体中并且承受远场速度梯度L∞边界条件，假设在远场的边界处施加的速度梯度张量L∞在RVE体积V内产生的流变场与直接在RVE边界处施加的速度梯度张量L所产生的流变场是一致的，即在RVE边界处施加L与在远场处施加的L∞是等效的，该种方法就是等效远场法[20]。同时，速度梯度张量L(L∞)可以分解拉伸张量D(D∞)和涡度张量W(W∞)，如图1所示。

图1 旋转的等效远场模型Fig.1 Rotational equivalent far field model

当椭球体为刚性夹杂时，仅发生刚性旋转变形，其拉伸张量为0。

根据向必伟等[14]给出的椭球体内流变场与基体内的流变场的关系，可得到椭球夹杂的涡度张量WE为

WE=W∞-Π:S-1:D∞

(11)

式(11)中：Π为反对称的Eshelby张量；S-1为对称的Eshelby张量的逆；WE和W∞分别为椭球体夹杂和RVE边界的涡度张量；D∞为RVE边界处的拉伸张量。

根据等效远场法，可分别建立拉伸和涡度张量与远场拉伸与涡度张量的关系：

D=(1-f)D∞

(12)

W=(1-f)W∞+fWE

(13)

式中:f为椭球夹杂所占的体积分数。

根据上述公式，可得到：

(14)

W∞=W+f·Π:S-1:D∞

(15)

将式(13)代入式(11)，并利用式(12)，可得到：

WE=W-Π:S-1:D

(16)

由上式可以看出，椭球夹杂的涡度张量与其体分比无关。

刚性椭球夹杂的旋转演化方程可以表示为[13]

(17)

式(17)中：Q为方位角矩阵；Θ为椭球体夹杂的角速度张量。

为提高计算效率，采用Basar等[21]给出的近似方法来求解上述微分方程：

(18)

(19)

对于可变形椭球夹杂的变形，其形状随着时间变化的关系式为[13]

(20)

3 数值模拟

基于上述细观力学模型，采用MATLAB数值软件，模拟了刚性椭球夹杂的演化进程。以文献[16]中的简单剪切流变场中给出的参数为例，研究刚性椭球夹杂的演化进程。采用的计算数据如下。

(1)椭球夹杂的3个主半轴分别为：a1=5，a1=3，a1=1。

(2)初始方位为(100,60,30)。

(4)每个步骤的时间增量δt=0.01。

(5)模拟运行的总步数为1 000，每隔20步取一次结果。

经过计算，可以得到刚性椭球夹杂的情况。刚性椭球体夹杂的演化进程如图2、图3所示。

“x”是起始点，“*”是终止点图2 单个刚体椭球夹杂Fig.2 Single rigid body ellipsoid inclusion

Step表示不同的时间步刚性椭球体夹杂的演化图3 刚体椭球夹杂的演化进程Fig.3 Evolution of rigid body ellipsoid inclusion

4 结论

基于等效远场法研究了黏性基体中刚性椭球夹杂的运动演化情况，得到如下结论。

(1)利用等效远场法建立了黏性基质中刚性椭球夹杂的细观力学理论，并基于MATLAB发展了数值计算程序。

(2)在给定速度梯度张量条件下，嵌入黏性材料的刚性椭球夹杂运动演化与体分比无关。

(3)基于发展的MATLAB程序，模拟了嵌入黏性基体中刚性椭球夹杂的运动情况，揭示了椭球夹杂的演化机制。