构造变形流变场分配多尺度数值模拟研究

向必伟， 姜大志， 谢成龙， 胡召奇

（1.安徽大学 资源与环境工程学院， 安徽 合肥 230601； 2.加拿大西安大略大学 地球科学系， 加拿大 伦敦N6A5B7； 3.合肥工业大学 资源与环境工程学院， 安徽 合肥 230009； 4. 安徽省地质调查院， 安徽 合肥230001）



构造变形流变场分配多尺度数值模拟研究

向必伟1， 姜大志2， 谢成龙3， 胡召奇4

（1.安徽大学 资源与环境工程学院， 安徽 合肥 230601； 2.加拿大西安大略大学 地球科学系， 加拿大 伦敦N6A5B7； 3.合肥工业大学 资源与环境工程学院， 安徽 合肥 230009； 4. 安徽省地质调查院， 安徽 合肥230001）

摘 要：流变场分配（flow field partitioning）现象在自然界高应变岩石中十分常见。传统的基于固体连续变形机制理论的变形分析中， 流变场分配问题往往被忽略， 致使应变带中流变场如何分配一直都缺乏深入认识。Eshelby阐述了嵌入均匀介质中的椭球体内流变场的数学方法， 为探讨流变场的分配奠定了理论基础。本文从固体连续变形机制入手， 重点介绍基于Eshelby理论的多尺度数值模拟思路和方法， 探讨流变场分配问题。模拟结果表明: 嵌入基质中的椭球体内分布的流变场主要取决于椭球体与基质间的相对流变强度， 椭球体的相对流变强度越低， 其变形越接近于简单剪切， 且有限应变积累越快。模拟还揭示， 不同流变强度的椭球体内模拟拉伸线理和面理产状的总体格局反映基质流变场特征。由此得到以下结论: （1）在流变场分配现象显著的区域， 局部小尺度上应变、涡度测量分析结果无法直接揭示区域流变场运动学边界条件， 对这些区域的构造变形分析必须是多尺度的； （2）基于Eshelby理论的以实质变形组构为约束的多尺度数值模拟分析， 能更为合理地揭示高应变岩石中流变场的分配。

关键词：流变场分配； 连续变形机制； Eshelby理论； 多尺度； 数值模拟

项目资助： 国家自然科学基金项目（41472194， 41002069， 40902064）资助。

0　引 言

构造地质学研究的一个关键目的是通过实际观察的变形构造揭示宏观区域变形的流变场及其对应的运动学边界条件。Lister和Williams （1979）引入了简单剪切和纯剪切变形的概念来描述岩石变形的流变场， Ramberg （1975）深化了这一概念并依据流变场中物质线与瞬时拉伸轴间的关系对变形形式进行了分类。同一时期研究者先后提出， 纯剪切和简单剪切是两种端员模式流变场， 自然界流变场是由这两个分量构成的， 两者的相对大小可用运动学涡度数（kinematic vorticity number）来衡量（Truesdell， 1953；Means et al.， 1980）。

Ramsay和Graham开创了岩石韧性变形的平面应变模式， 并将应变带假设为局限于刚性边界内的均匀连续体， 应变带内流变场由边界的相对运动产生（Ramsay and Graham， 1970）。在Ramsay变形分析思路的启发下， 许多学者对自然界应变带的流变场模式开展了深入研究， 流变场模式也从二维平面应变逐渐推广到三维单斜变形模式（Sanderson and Marchini， 1984； Fossen and Tikoff， 1993； Tikoff and Fossen， 1993）和一般三斜变形模式（Jiang and Williams，1998； Lin et al.， 1998）。Ramasy及其后续变形模型的理论基础是连续统内的均匀介质变形（Jiang， 2013），其中任意一物质点的运动学直接反应了连续统的宏观流变场特征。基于这一假设， 实际构造变形分析中研究者通常用从任意尺度构造（如: 露头构造、手标本构造、显微构造）获得的流变场代表宏观区域流变场。

然而， 局部流变场与区域流变场往往并不一致（Kuiper et al.， 2011）， Lister and Williams （1983）称之为流变场分配（flow field partitioning）。在野外观察的基础上人们对流变场分配取得了一些直观的认识，如自然界的应变集中带往往比围岩“软弱”， 易于变形； 简单剪切变形往往集中在“软弱”的应变局限化带内（Michibayashi and Murakami， 2007； Kuiper et al.，2011； Lee et al.， 2012； Bell et al.， 2013）。然而， 迄今为止对于流变场如何在不均匀的变形区域内分配以及局部分配的流变场与宏观区域流变场如何联系仍然缺乏深入了解。正因如此， 虽然流变场分配的概念早已被构造研究者认同， 并被用来解释一些区域内局部的不均匀变形（Wilson and Powell， 2001；Misra et al.， 2009； Shigematsu et al.， 2009； Carreras et al.， 2013； Aydin and de Joussineau， 2014）， 然而在实际研究应变带的整体运动学边界条件时却极少考虑流变场分配因素。

Eshelby （1957， 1959）在连续变形假设基础上，建立了均匀弹性基质与嵌入其中的均匀椭球体之间流变场的数学联系。Eshelby公式被拓展到探讨嵌入在线性粘性牛顿体材料（Newtonian）中椭球体的变形（Bilby et al.， 1975； Ellis and Watkinson， 1987）。牛顿体的流变学强度一般用粘性系数表述， 而对于非线性粘性材料（power-law）， 其流变学强度受诸多因素影响， 如流变指数、流变速率等。为了研究非线性粘性材料基质与嵌入的椭球体的变形， 非线性粘性材料的流变学强度用变形过程中的有效粘度（effect viscosity）描述， Lebensohn and Tomé （1993）、Jiang and Bentley （2012）以及Jiang （2013）对非线性粘性材料的有效粘度及其数值算法进行了深入开发。Eshelby公式奠定了研究宏观区域流变场在局部变形体内分配的理论基础。众多研究者利用给定的基质流变场代表区域流变场， 基质中嵌入的椭球体代表局部地质体， 建立了多尺度流变数值模型（Jiang and Bentley， 2012； Jiang， 2013； Xiang and Jiang，2013）。这一模型在局部变形体与区域运动学边界条件之间的鸿沟上建立了桥梁。

本文介绍基于Eshelby理论建立3D多尺度数值模型， 模拟研究流变场在不同尺度上的分配。我们将给定的基质流变场对应宏观尺度上的运动学边界条件， 基质中可变形的椭球体视为局部地质体， 椭球体内部的有限应变主轴模拟变形组构。这一简单模型可以模拟三个尺度间的构造变形关系，例如宏观尺度、露头尺度以及组构尺度， 因此这一模型是多尺度的。模拟计算揭示， 椭球体相对于基质的流变学强度越低， 其变形越倾向于简单剪切变形， 并且有限应变的积累越快； 单一椭球体内的流变场都不等同于宏观流变场， 但全部椭球体内发育的模拟面理和拉伸线理产状的总体格局能反映基质流变场的特点。模拟结果揭示天然岩石中流变场分配主要取决于岩石的相对流变强度。模拟结果还表明， 由于流变场的分配， 对自然界应变带的宏观流变场和对应的运动学边界条件研究，利用传统的以变形材料均匀为前提假设的变形分析方法是不合适的。相对于单一尺度的连续统变形理论， 基于Eshelby理论的多尺度构造变形分析思路在很大程度上避免了岩石的均一性假设， 更具普遍意义。

1　数值模拟理论背景

1.1坐标系与参考系

数值模型的建立依赖于参考坐标系， 在此我们首先定义下文所涉及的参考坐标系和相关参数。在图1a定义的右手固定坐标系x1x2x3中， 均匀高应变带的走向平行于x1轴， 任意给定边界速度v。边界速度沿坐标轴方向分解为v1， v2和v3（图1a）， 相应在应变带内产生分别平行于x1， x2， x3的拉伸应变率，，和平行于x1， x3的剪应变率和（图1b）。空间内任意一个向量e的方向可用两个角度θ和φ确定， θ为向量e在x1x2平面上的投影与x1轴正方向的夹角， φ为向量与x3轴正方向的夹角（图1c）。空间内任意一个椭球的方位由三个椭球半径轴方位确定（三对θ和φ角）， 其中只有三个角度是独立的。我们取最大椭球半径轴的θ和φ角和中间半径轴的一个方位角φ （图1d）， 当φ = 90°时， φ为椭球中间轴在x1x2平面的投影与x1轴正方向的夹角， 但当φ=90°时， φ为椭球中间半径轴与x3轴正方向的夹角。

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

图1　数值模型坐标系参数示图Fig.1 Diagrammatic drawing for the coordinate system and parameters of the numerical model

均匀介质的流变学在岩石流变学方面的应用始于简单的二维变形研究。Ramsay和Graham遵循有限应变思路阐述了刚性边界条件下变形带恒体积简单剪切型式的变形机制（平面应变）（Ramsay and Graham， 1970； Ramsay， 1980）。Ramberg （1975）详尽阐述了等体积、连续、稳态条件下岩石有限应变的二维一般剪切模式。Sanderson and Marchini （1984）提出了走滑挤压变形（transpressional zones）的三维单斜分析模式。随后大量三维分析模式陆续发表（Fossen and Tikoff， 1993； Krantz， 1995； Jones and Geoff， 1995； Teyssier et al.， 1995）。在此基础上， Jiang and Williams （1998）提出了一般的三斜变形模式， 探讨了三斜变形与以往研究的单斜变形的关系， 推演了三斜变形的速度梯度张量和位移梯度张量。

给定边界速度（图1a）均匀介质变形的流变场用速度梯度张量L来表达:

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

1.3Eshelby理论

Eshelby （1957）提出了处理均匀弹性椭球嵌入体与均匀弹性基质间应变场关系的数学方法。这一方法随后被拓展到应用于线性牛顿体材料（Bilby et al.，1975）乃至非线性的粘性（power-law）材料（Lebensohn and Tomé， 1993； Jiang and Bentley， 2012； Jiang，2013）。牛顿体材料， 流变强度可由粘度表达。非线性粘性体， 其粘度由应力指数、应变率以及其他热力学性质决定， 在变形过程中呈非线性变化。设想在变形过程中的任一无限短时间内， 非线性粘性体的实际粘度与应变率间的非线性关系可以用一个线性方程近似拟合； 因此， 在这一时间段内的变形也近似地满足Eshelby理论（Lebensohn and Tomé， 1993；Jiang and Bentley， 2012）。此时， 非线性粘性体的实际粘度被称为“有效粘度”（Lebensohn and Tomé，1993）， 在Eshelby模型中嵌入的椭球体与基质间的有效粘度比值我们称为有效粘度比（reff）。Jiang and Bentley （2012）系统地探讨了非线性粘性体有效粘度的近似方法、计算精度以及数值方法。

Eshelby （1957）在数学上证明当均匀基质中的嵌入体为椭球体时， 变形过程中椭球体内部的应力和应变都是处处一致的， 但随着椭球的几何形状与方位的变化， 椭球内部流变场是非稳态的。对于各向同性非线性粘性材料， 椭球体内流变场与远离嵌入体的基质内流变场满足如下公式（Jiang， 2013）:

两式中DE和DM分别是椭球嵌入体内部和远离嵌入体的基质内的拉伸张量； WE和WM分别是椭球嵌入体内部和远离嵌入体的基质内的涡度张量。JS是四阶对称单位张量； S和Π分别是对称和反对称的Eshelby张量。A是四阶应变率分布（Partitioning） 张量， 对各向同性材料其表达式为:

式（3）中n是非线性粘性基质材料的应力指数， reff是椭球嵌入体与基质材料的有效粘度比。由于椭球嵌入体的有效粘度取决于其变形时的真实应力（或者应变率）， 嵌入体材料的有效粘度可通过Mancktelow和Jiang开发的迭代程序进行数值求解（Mancktelow，2013； Jiang， 2013）。

当给定远离嵌入体的基质中流变场速度梯度张量LM， 将LM分解为DM和WM， 我们就可以通过公式（2a）， （2b）及公式（3）求解分布于椭球嵌入体内部流变场的DE和WE， 从而求解速度梯度张量LE。在以下模拟计算中我们假设基质内的流变场是均匀且稳态的。此时， LM是恒定的， 并可由内部运动学涡度刻画。由于椭球嵌入体的形状和方向是不断变化的， 嵌入体内部的流变场虽然是均匀的但是非稳态的。由公式（2）确定的涡度张量WE， 包含有刚体旋转分量。只有内部涡度才能真正刻画流变场LE的本质， 利用Jiang （2014）开发的应用程序我们可以计算LE的内部运动学涡度数。显然， 椭球嵌入体内部流变场的运动学涡度是随时间演化的。通过求解增量应变， 我们可以求解椭球嵌入体经历任意时间段的变形后的有限应变量大小、主应变大小以及应变主轴的方位。

2　数值模拟

天然构造变形材料（岩石）组成成分往往十分多样、结构上具有多尺度特征。从变形岩石整体上（大尺度）看， 岩石由许多形状各异、流变强度不同的次级尺度元素组成。显然， 建立与天然岩石对等的理论模型来研究岩石变形行为， 是岩石流变学研究的理想目标。然而， 受限于现有的材料力学理论以及对天然岩石力学性质的了解， 这一目标暂时还难以实现。因此， 对天然岩石构成进行适当简化建立相对简单的相似模型就成为必然。目前， 材料力学领域在研究非均匀材料变形行为时， 往往将非均匀材料对等为一个力学性质十分相近的均匀介质， 这个介质被称为“均匀对等介质（HEM）”（Mura， 1987； Lebensohn et al.， 1998）。借鉴这一思路， 我们在研制岩石整体（大尺度）变形时也将其视为一个均匀介质， 在研究岩石中某一元素（次级尺度）时将其视为嵌入均匀介质中的一个非均匀体。这样， 我们就可以依据Eshelby理论建立简单的多尺度数值模型， 来模拟研究岩石的变形行为。

2.1数值模型

我们构建一个由单一椭球体嵌入基质材料的数值模型来模拟宏观区域构造与局部地质体。数值模型中， 椭球体与基质都是均匀的各向同性的非线性粘性（power-law）材料， 两者的有效粘度不同。以同一应变率状态为参照， 我们用椭球体材料与基质材料的流变应力指数ne和nm， 以及椭球体与基质材料的有效粘度比（reff）定义椭球体与基质材料的相对流变学性质（Jiang and Bentley， 2012）。椭球体的方向由三个球面角θ， φ和φ确定 （图1d）， 椭球几何形状由椭球半径比值（a1: a2: a3）确定（a1， a2， a3分别是最大半径、中间半径和最小半径）。

在图1a坐标系内给定恒定的基质流变场LM。除椭球嵌入体外， 基质的变形可视为均匀介质的稳态变形， LM就等同于等式1所定义的L。我们假设变形是恒体积的， 即。变形过程中跟踪计算椭球嵌入体内流变场速度梯度张量LE， 从而跟踪流变场的运动学涡度数、有限应变主轴的方位以及有限应变量大小。显然， 所有相关的参数都会影响椭球嵌入体内的流变场， 包括:

通过系统地调整上述参数， 运行大量不同椭球体初始形态和方位模拟计算， 考察基质流变场在椭球体内的分配特征以及影响流变场分配的因素。模拟过程中， 首先系统地考察椭球体在平面应变型式下的基质流变场中发生平面变形时流变场的分配特征（图2）， 然后逐步考察一般三斜的基质流变场在椭球体内的分配。本次研究中， 我们主要用椭球体中流变场的运动学涡度数和有限应变量大小来刻画流变场的分配特征。

图2　平面应变椭球的空间方位Fig.2 Orientation of the plane straining ellipsoid

2.2椭球体内的流变场特征

图3和图4分别揭示了平面应变的椭球嵌入体内流变场的运动学涡度数和有限应变量的演化历史。在相对粘度低的椭球嵌入体内（reff＜ 1）， 流变场的运动学涡度数经历了初始阶段的起伏变化后很快趋于稳定（NM≈ 2）。reff越低， wk越易于达到稳定， 且最终稳定的wk值越接近于1（图3）。本次模拟计算设定的椭球体与基质流变强度差异处于自然界岩石强度差异范围内， 低流变学强度的椭球体的变形接近于简单剪切。例如， 当基质流变场的wk=0.6， 在有效粘度是基质粘度的五分之一的椭球嵌入体内， 流变场的wk约为0.95。椭球体与基质积累有限应变的速率也有明显差异（图4）。尽管基质和椭球嵌入体内有限应变的积累都是非线性的， 但两者的比值近似线性。有效粘度越低， 椭球嵌入体积累有限应变的速度越快。当有效粘度比并不十分显著时， 低流变学强度的椭球嵌入体积累的有限应变就可达到基质材料的数倍， 如: 当reff≈0.2时， NE/ NM≈3.5（图4f）。

平面应变流变场的分配规律也适合于一般三维流变场。如图5所示， 不同型式的基质流变场在一任意初始方位的三维椭球体内的分配特点基本相同。椭球体内流变场的运动学涡度数经历短暂的波动后很快趋于稳定， 且稳定的涡度数完全取决于椭球体与基质间的相对有效粘度。当给定的基质流变场的涡度数一致（图5所示， wk=0.6）， 无论基质流变场是平面应变型式（图5a）或一般单斜型式（图5c）还是一般三斜型式（图5e）， 有效粘度比近似相等的椭球体内流变场的最终稳定涡度数近似相等。在不同类型基质流变场中， 基质与椭球体内的有限应变量仍然呈现了近似的线性关系， 并且有效粘度越低的椭球体内有限应变积累速率越快（图5b， d 和f）。

图3　平面应变椭球体内流变场运动学涡度演化图Fig.3 Evolution of the kinematical vorticity number of the flow field in plane straining ellipsoids

上述模拟计算表明， 流变场在不同流变学强度元素组成的介质中呈现差异分布。介质的相对流变学强度（粘度）越低， 其变形越接近于简单剪切变形，且有限应变积累得越快。自然界中普遍存在的应变非均匀分配（strain partition）的地质现象（Lister and Williams， 1983； Ramsay， 1980； Hudleston， 1999；Hobbs et al.， 2011）与我们的理论结果是完全吻合的。

2.3组构与区域流变场

图4　平面应变椭球有限应变演化图Fig.4 Evolution of the finite strain in plane straining ellipsoids

我们假设自然界高应变带由方向不同、形状各异的地质体构成。在图1a所示的参考坐标系下， 将 300个在三维空间内随机定向的椭球体嵌入到同一均匀基质中， 且假设任意两个椭球体之间的距离足够远， 变形过程中不会产生相互影响。椭球的长半径与中间半径的比值在1到10之间随机变化， 短半径给定为1。这些椭球能够覆盖很大范围的形态变化， 从近等轴的球状（1∶1∶1）到层状（10∶10∶1），再到杆状（10∶1∶1）。我们同样假设椭球嵌入体与基质都是非线性粘性材料， 给定材料的应力指数均为3。椭球嵌入体与基质间的相对粘度比在0.1到0.5范围内随机给定。以一般单斜形式的基质流变场为例模拟椭球嵌入体内拉伸线理与面理产状的总体演化。此时， LM可写为:

等式（4）的LM刻画了一个垂直的右行走滑应变带。其中，是平行于x2轴的缩短应变率， μ从0到1之间取值，˙可由LM的运动学涡度数wk和μ确定。

图5　三维椭球体内流变场的运动学涡度数和有限应变演化图Fig.5 Evolution of the kinematical vorticity number and finite strain in a 3D ellipsoid embedded in different imposed flow field

基质流变场特征决定了组构产状的整体格局（图6）。椭球嵌入体内有限应变最小主轴S3的方位几乎不受wk和μ的影响， 总是近平行于x2轴； 而最大主应变轴S1的方位则由wk和μ控制。S1的方位可以从以x3轴为中心的点集密经大圆环带过渡到以x1轴为中心的点集密。其代表的拉伸线理产状投影从垂直的点集密经过大圆环带过渡到水平点集密（图6b）。

图6　理论模拟变形组构演变图Fig.6 Evolution pattern of deformation fabrics from numerical modeling

3　讨 论

3.1Eshelby理论在韧性变形研究中的优势

与传统的韧性变形分析理论相比较， 利用Eshelby理论开展多尺度韧性变形分析的最主要优势在于其假设了天然岩石构成的非均匀性， 能够合理地解释自然界流变场的不均匀分配现象。自从Ramsay利用连续变形理论分析岩石构造变形以来，这一理论在岩石韧性变形分析领域一直占有统治地位。岩石材料均匀是Ramsay理论的基本假设之一。在这一前提下， 自然界构造变形区域流变场的不均匀分配现象始终难以得到合理的解释。例如在均匀性假设条件下， 人们通常利用材料弹性失稳模式来解释应变局限化带的形成机制。而按照这一观点，高应变带的尾端相容性问题（Hudleston， 1999）始终难以解决； 尽管Ramsay （1980）提出了两种解释模式，但到目前为止未见实际资料的报道。Eshelby理论的初始边界条件是均匀介质中嵌入可变形的均匀椭球体， 基质和椭球体具有不同的流变强度。这就强调了介质与嵌入其中的椭球体间的不均匀特征。随着均匀对等介质（HEM）思路的提出以及可行的数值求解方法的开发， 基于Eshelby理论的模型可以视为由流变强度不同的椭球体构成（Jiang， 2014）； 研究其中一个椭球体内的流变场时， 其他椭球体整体视为基质， 基质的流变性质用实时求解的HEM的流变性质近似（Jiang and Bentley， 2012）。显然， 通过改进的Eshelby理论模型进一步强调了变形材料的不均匀性， 与天然岩石的构成也就更加相似。上述我们利用简单的Eshelby理论模型的模拟计算， 揭示了天然岩石变形流变场的不均匀分配主要受岩石相对流变强度制约， 这一规律与已有的观察结论一致。

与其他数值模拟方法相比， 基于Eshelby理论的数值模拟的优势在于大应变和多对象同时研究。目前用于变形模拟研究的常用数值模拟方法主要有有限元法和有限差分法， 两种方法都有成熟的商业软件， 也都应用于模拟岩石的韧性变形（Mancktelow，2011）。两种方法的一个共同特点是将模型网格化，通过求解网格节点的运动学和力学参数来描述变形过程。通过对模型中不同区域内的网格设定不同的力学参数， 也可以建立由不均匀材料构成的数值模型。模拟过程中， 当任意相邻的两个网格交差产生几何破坏时， 模型材料变形的应力、应变连续性被破坏， 运算就无法继续进行。已有的模拟研究经验表明， 网格密度越大， 材料的力学性质刻画得就越精确， 计算精度也就越高， 但此时运算速率也越慢且网格越容易发生几何破坏， 模型材料经历的有限应变量就越小。因此， 无法实现模型与天然岩石应变相当的大应变是基于有限元理论和有限差分理论的模拟方法的主要困难之一， 特别是刻画具有较大流变学强度差异的非均匀材料模型尤为如此。理论上， 基于有限元理论和有限差分理论的数值模型可以无限包含非均匀元素， 但实际模拟实验过程中，这类模型能够发生的有限应变量有限而且运算量也无法由一般硬件设备承担。因此， 已有的基于这两种数学理论在岩石变形方面的研究模型材料都较为单一， 往往仅包含一个或有限的几个非均匀元素（Mancktelow， 2011）。已有研究表明（Jiang and Bentley，2012； Jiang， 2013； Xiang and Jiang， 2013）， 基于Eshelby理论的数值模型可以开展大应变模拟计算，并且可以创建由大量流变学性质不同的元素构成的数值模型。

3.2Eshelby数值模型在岩石变形中的局限性

同传统的岩石韧性变形分析理论一样， Eshelby理论也以连续变形假设为前提， 这一假设对于天然岩石变形可能过于严格。天然岩石在构造变形过程中， 不同块体（或元素）之间不仅发生连续变形， 还可能会发生相对滑移的非连续变形。显然， 如果岩石变形过程中非连续变形十分显著， 用Eshelby理论建立的连续变形模式来研究其变形行为就不合适了。然而， 在考虑引入新的理论建立包含非连续变形基质的模型之前， 我们必须对岩石变形行为开展系统的研究， 揭示岩石发生非连续变形的标志， 非连续变形对岩石整体变形的贡献大小等等。如果包含非连续变形的理论模型是研究天然岩石变形的最终理想模型， 那么基于Eshelby理论的数值模型只是人们迈向最终目标的一个中间过程。

Eshelby理论对嵌入基质中的次级元素的假设条件可能过于严格。一方面， 地质体的初始几何形态可能与标准椭球体有一定出入， 这样在实际变形过程中其内部的流变场并非严格的处处均一。另一方面， 利用均匀椭球体模拟地质体也是简化地近似；天然地质体的流变性质往往并不完全均匀， 变形过程中， 流变学性质的不均匀性也会导致实际流变场比理想的模拟流变场复杂得多。因此， 实际天然变形岩石中记录流变场特征的组构可能比数值模拟预测结果更为复杂， 数值模拟结果更多地揭示了岩石的总体变形行为（Xiang and Jiang， 2013）。此外， 天然地质体还可能具有各向异性的流变学性质。尽管对一些相对简单的各向异性体的变形已开展了一些初步的模拟研究工作， 如云母残斑晶的变形和定向（Chen et al.， 2014）， 但是多数情况下地质体的各向异性特征并不被人们所掌握， 对构建具有各向异性特征的数值模型缺乏基本的约束。因此， 地质体的各向异性特征也会影响地质体内部的流变场以及变形组构发育。

3.3韧性变形分析思路

本次模拟揭示， 构造变形区域内分配在不同流变学强度的局部岩石内的流变场往往不同。因此， 传统的以均匀性假设为前提的变形分析方法就难以胜任对流变场不均匀分配区域的宏观流变场特征的研究了。首先， 以往的韧性构造变形分析中， 拉伸线理方向或优选方位直接代表宏观流变场的剪切方向。本次模拟研究表明， 变形区域内拉伸线理产状的总体分布格局是由宏观流变场形式决定的， 随着流变场特征的变化， 拉伸线理优选方位可以从与剪切方向一致渐变为与剪切方向垂直。模拟实验表明， 可以利用变形区域内拉伸线理产状的总体分布格局来揭示宏观流变场形式。其次， 传统研究方法中通过运动学涡度测量来直接指示区域上的变形程度和运动学边界条件， 在流变场不均匀分配的区域就不合适了。本次模拟研究表明， 局部岩石内流变场的运动学涡度数不仅由宏观流变场的运动学涡度决定， 而且还受局部岩石的相对流变学强度影响。对于如何利用从不同局部变形岩石中测量的运动学涡度数来揭示宏观区域流变场的运动学涡度数， 从而揭示宏观运动学边界条件， 还需要进一步研究， 我们将在后续研究中进一步探讨。

4　结 论

本次工作系统地综述了基于Eshelby理论的多尺度数值模拟研究思路和方法。通过多尺度数值模拟， 我们对天然岩石的构造变形研究取得了两点初步认识: （1）构造变形区域， 流变场的分配主要受岩石的相对流变学强度影响， 流变强度越低的岩石内流变场的剪应变越集中且有限应变累积速度越快；（2）在流变场不均匀分配的区域开展构造变形分析，传统的以变形岩石均匀性假设为前提的韧性变形分析方法难以胜任； 变形区域上宏观流变场可以通过面理和拉伸线理组构产状的总体分布格局并结合可能的变形形式来推断； 构造变形区域上流变场的运动学涡度数无法从测量局部变形岩石中的运动学涡度数直接获取。

致谢： 两位审稿专家对本文的最终成稿提出了科学严谨的修改意见， 让笔者获益匪浅， 在此表示由衷的感谢。

参考文献（References）:

Aydin A and de Joussineau G. 2014. The relationship between normal and strike-slip faults in Valley of Fire State Park， Nevada， and its implications for stress rotation and partitioning of deformation in the east-central Basin and Range. Journal of Structural Geology， 63: 12-26.

Bell T H， Rieuwers M T， Cihan M， Evans T P， Ham A P and Welch P W. 2013. Inter-relationships between deformation partitioning， metamorphism and tectonism. Tectonophysics， 587: 119-132.

Bilby B A， Eshelby J D and Kundu A K. 1975. The change of shape of a viscous ellipsoidal region embedded in a slowly deforming matrix having a different viscosity. Tectonophysics， 28（4）: 265-274.

Carreras J， Cosgrove J W and Druguet E. 2013. Strain partitioning in banded and/or anisotropic rocks: Implications for inferring tectonic regimes. Journal of Structural Geology， 50: 7-21.

Chen Y， Jiang D， Zhu G and Xiang B. 2014. The formation of micafish: A modeling investigation based on micromechanics. Journal of Structural Geology，68（Part B）: 300-315.

Elliott D. 1972. Deformation paths in structural geology. Geological Society of America Bulletin， 83（9）: 2621-2638.

Ellis M and Watkinson A J. 1987. Orogen-parallel extension and oblique tectonics: The relation between stretching lineations and relative plate motions. Geology， 15（11）: 1022-1026.

Eshelby J D. 1957. The determination of the elastic field of an ellipsoidal inclusion， and related problems. Proceedings of the Royal Society of London. Series A，Mathematical and Physical Sciences， 241（1226）: 376-396.

Eshelby J D. 1959. The elastic field outside an ellipsoidal inclusion. Proceedings of the Royal Society of London. Series A， Mathematical and Physical Sciences，252（1271）: 561-569.

Fossen H and Tikoff B. 1993. The deformation matrix for simultaneous simple shearing， pure shearing and volume change， and its application to transpressiontranstension tectonics. Journal of Structural Geology，15（3-5）: 413-422.

Goscombe B D and Gray D R. 2008. Structure and strain variation at mid-crustal levels in a transpressional orogen: A review of Kaoko Belt structure and the character of West Gondwana amalgamation and dispersal. Gondwana Research， 13（1）: 45-85.

Hobbs B E， Ord A and Regenauer-Lieb K. 2011. The thermodynamics of deformed metamorphic rocks: A review. Journal of Structural Geology， 33（5）: 758-818.

Hudleston P. 1999. Strain compatibility and shear zones: Is there a problem？ Journal of Structural Geology，21（8-9）: 923-932.

Jiang D. 2007. Sustainable transpression: An examination of strain and kinematics in deforming zones with migrating boundaries. Journal of Structural Geology，29（12）: 1984-2005.

Jiang D. 2013. The motion of deformable ellipsoids in power-law viscous materials: Formulation and numerical implementation of a micromechanical approach applicable to flow partitioning and heterogeneous deformation in Earth's lithosphere. Journal of Structural Geology， 50: 22-34.

Jiang D. 2014. Structural geology meets micromechanics: A self-consistent model for the multiscale deformationand fabric development in Earth's ductile lithosphere. Journal of Structural Geology， 68（Part B）: 247-272.

Jiang D and Bentley C. 2012. A micromechanical approach for simulating multi-scale fabrics in large-scale high-strain zones: Theory and application. Journal of Geophysical Research， 117（12）: 12201-12216.

Jiang D and Williams P F. 1998. High-strain zones: A unified model. Journal of Structural Geology， 20（8）: 1105-1120.

Jones R and Geoff T P W. 1995. Strain partitioning in transpression zones. Journal of Structural Geology，17（6）: 793-802.

Krantz R W. 1995. The transpressional strain model applied to strike-slip， oblique-convergent and obliquedivergent deformation. Journal of Structural Geology，17（8）: 1125-1137.

Kuiper Y D， Lin S and Jiang D. 2011. Deformation partitioning in transpressional shear zones with an along-strike stretch component: An example from the Superior Boundary Zone， Manitoba， Canada. Journal of Structural Geology， 33（3）: 192-202.

Lebensohn R A and Tomé C N. 1993. A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: Application to zirconium alloys. Acta Metallurgica et Materialia， 41（9）: 2611-2624.

Lebensohn R A， Turner P A， Signorelli J W， Canova G R and Tomé C N. 1998. Calculation of intergranular stresses based on a large-strain viscoplastic self-consistent polycrystal model. Modelling and Simulation in Materials Science and Engineering， 6（4）: 447-465.

Lee P E， Jessup M J， Shaw C A， Hicks Iii G L and Allen J L. 2012. Strain partitioning in the mid-crust of a transpressional shear zone system: Insights from the Homestake and Slide Lake shear zones， central Colorado. Journal of Structural Geology， 39: 237-252.

Lin S， Jiang D and Williams P F. 1998. Transpression （or transtension） zones of triclinic symmetry: Natural example and theoretical modelling. Geological Society，London， Special Publications， 135（1）: 41-57.

Lister G S and Williams P F. 1979. Fabric development in shear zones: Theoretical controls and observed phenomena. Journal of Structural Geology， 1（4）: 283-297.

Lister G S and Williams P F. 1983. The partitioning of deformation in flowing rock masses. Tectonophysics， 92（1-3）: 1-33.

Mancktelow N S. 2011. Deformation of an elliptical inclusion in two-dimensional incompressible power-law viscous flow. Journal of Structural Geology， 33（9）: 1378-1393.

Mancktelow N S. 2013. Behaviour of an isolated rimmed elliptical inclusion in 2d slow incompressible viscous flow. Journal of Structural Geology， 46: 235-254.

Means W D， Hobbs B E， Lister G S and Williams P F. 1980. Vorticity and non-coaxiality in progressive deformations. Journal of Structural Geology， 2（3）: 371-378.

Michibayashi K and Murakami M. 2007. Development of a shear band cleavage as a result of strain partitioning. Journal of Structural Geology， 29（6）: 1070-1082.

Misra S， Burlini L and Burg J P. 2009. Strain localization and melt segregation in deforming metapelites. Physics of The Earth and Planetary Interiors， 177（3-4）: 173-179.

Mura T. 1987. Micromechanics of Defects in Solids. Martinus Nijhoff Publishers: 388.

Ramberg H. 1975. Particle paths， displacement and progressive strain applicable to rocks. Tectonophysics，28（1-2）: 1-37.

Ramsay J G. 1980. Shear zone geometry: A review. Journal of Structural Geology， 2（1-2）: 83-99.

Ramsay J G and Graham R H. 1970. Strain variation in shear belts. Canadian Journal of Earth Sciences， 7（3）: 786-813.

Sanderson D J and Marchini W R D. 1984. Transpression. Journal of Structural Geology， 6（5）: 449-458.

Shigematsu N， Fujimoto K， Ohtani T， Shibazaki B， Tomita T，Tanaka H and Miyashita Y. 2009. Localisation of plastic flow in the mid-crust along a crustal-scale fault: Insight from the Hatagawa Fault Zone， NE Japan. Journal of Structural Geology， 31（6）: 601-614.

Teyssier C， Tikoff B and Markley M. 1995. Oblique plate motion and continental tectonics. Geology， 23（5）: 447-450.

Tikoff B and Fossen H. 1993. Simultaneous pure and simple shear: The unifying deformation matrix. Tectonophysics，217（3-4）: 267-283.

Truesdell C A. 1953. Two measures of vorticity. Journal of Rational Mechanics and Analysis， 2: 173-217.

Truesdell C A and Toupin R A. 1960. The classic field theories // Flügge S. Encyclopedia of Physics， Vol. II: Principles of Classical Mechanics and Field Theory. Springer-Verlag， Berlin: 226-793.

Wilson C J L and Powell R. 2001. Strain localisation and high-grade metamorphism at Broken Hill， Australia: A view from the Southern Cross area. Tectonophysics，335（1-2）: 193-210.

Xiang B and Jiang D. 2013. Small-scale ductile shear zones as transposed rheologically weak domains: A numerical modeling investigation and practical application. Journal of Structural Geology， 54: 184-198.

Multi-scale Numerically Modeling Flow Field Partitioning During Structural Deformation

XIANG Biwei1， JIANG Dazhi2， XIE Chenglong3and HU Zhaoqi4
（1. School of Resources and Environmental Engineering， Anhui University， Hefei 230601， Anhui， China； 2. Department of Earth Sciences， University of Western Ontario， London N6A5B7， Canada； 3. School of Resources and Environmental Engineering， Hefei University of Technology， Hefei 230009， Anhui， China； 4. Geological Survey of Anhui Province， Hefei 230001， Anhui， China）

Abstract：Flow field partitioning is common in inhomogeneous rocks in high strain area. However， the classic structural analysis based on the continuum mechanics usually ignores the ways of flow field partition. Eshelby's formulation deals with the flow field in an ellipsoid embedded in a homogeneous matrix， which established the foundation for the study of flow field partitioning. In this paper， we begin with the continuum mechanism and focus on the multi-scales structural analysis concept and method by numerical modeling based on the Eshelby's formulation. Numerical modeling results show that the flow field in embedding ellipsoids depends on the relative rheology between the ellipsoid and matrix. Under an imposed matrix flow field， the weaker the ellipsoid is， the more noncoaxial the partitioned flow field is and the more dramatic finite strain increase. Modeling results also show that the imposed flow field can be characterized by the whole pattern composed of fabrics in all embeddings with different rheology. Thereafter， we got two conclusions: （1） in a flow field partitioning area， the traditional structural analysis by the finite and vorticity number measurements can not be directly used to explain the regional kinematical boundary conditions； （2） the traditional single-scale method is not enough to explain flow field partitioning and a multi-scale modeling limited by natural fabrics should be a feasible way to understand the kinematic boundary conditions of a regional deformation.

Keywords:flow field partitioning； continuum mechanism； Eshelby's theory； multi-scale； numerical modeling

中图分类号：P542

文献标志码：A

文章编号：1001-1552（2016）01-0001-013

收稿日期：2014-09-27； 改回日期： 2015-05-08

第一作者简介：向必伟（1976-）， 男， 副教授， 从事造山带构造变形分析等方面研究。Email: xbw1977@163.com