構造變形流變場分配多尺度數值模擬研究

向必偉， 姜大志， 謝成龍， 胡召奇

（1.安徽大學 資源與環(huán)境工程學院， 安徽 合肥 230601； 2.加拿大西安大略大學 地球科學系， 加拿大 倫敦N6A5B7； 3.合肥工業(yè)大學 資源與環(huán)境工程學院， 安徽 合肥 230009； 4. 安徽省地質調查院， 安徽 合肥230001）

?

構造變形流變場分配多尺度數值模擬研究

向必偉1， 姜大志2， 謝成龍3， 胡召奇4

（1.安徽大學 資源與環(huán)境工程學院， 安徽 合肥 230601； 2.加拿大西安大略大學 地球科學系， 加拿大 倫敦N6A5B7； 3.合肥工業(yè)大學 資源與環(huán)境工程學院， 安徽 合肥 230009； 4. 安徽省地質調查院， 安徽 合肥230001）

摘 要：流變場分配（flow field partitioning）現象在自然界高應變巖石中十分常見。傳統(tǒng)的基于固體連續(xù)變形機制理論的變形分析中， 流變場分配問題往往被忽略， 致使應變帶中流變場如何分配一直都缺乏深入認識。Eshelby闡述了嵌入均勻介質中的橢球體內流變場的數學方法， 為探討流變場的分配奠定了理論基礎。本文從固體連續(xù)變形機制入手， 重點介紹基于Eshelby理論的多尺度數值模擬思路和方法， 探討流變場分配問題。模擬結果表明: 嵌入基質中的橢球體內分布的流變場主要取決于橢球體與基質間的相對流變強度， 橢球體的相對流變強度越低， 其變形越接近于簡單剪切， 且有限應變積累越快。模擬還揭示， 不同流變強度的橢球體內模擬拉伸線理和面理產狀的總體格局反映基質流變場特征。由此得到以下結論: （1）在流變場分配現象顯著的區(qū)域， 局部小尺度上應變、渦度測量分析結果無法直接揭示區(qū)域流變場運動學邊界條件， 對這些區(qū)域的構造變形分析必須是多尺度的； （2）基于Eshelby理論的以實質變形組構為約束的多尺度數值模擬分析， 能更為合理地揭示高應變巖石中流變場的分配。

關鍵詞：流變場分配； 連續(xù)變形機制； Eshelby理論； 多尺度； 數值模擬

項目資助： 國家自然科學基金項目（41472194， 41002069， 40902064）資助。

0　引 言

構造地質學研究的一個關鍵目的是通過實際觀察的變形構造揭示宏觀區(qū)域變形的流變場及其對應的運動學邊界條件。Lister和Williams （1979）引入了簡單剪切和純剪切變形的概念來描述巖石變形的流變場， Ramberg （1975）深化了這一概念并依據流變場中物質線與瞬時拉伸軸間的關系對變形形式進行了分類。同一時期研究者先后提出， 純剪切和簡單剪切是兩種端員模式流變場， 自然界流變場是由這兩個分量構成的， 兩者的相對大小可用運動學渦度數（kinematic vorticity number）來衡量（Truesdell， 1953；Means et al.， 1980）。

Ramsay和Graham開創(chuàng)了巖石韌性變形的平面應變模式， 并將應變帶假設為局限于剛性邊界內的均勻連續(xù)體， 應變帶內流變場由邊界的相對運動產生（Ramsay and Graham， 1970）。在Ramsay變形分析思路的啟發(fā)下， 許多學者對自然界應變帶的流變場模式開展了深入研究， 流變場模式也從二維平面應變逐漸推廣到三維單斜變形模式（Sanderson and Marchini， 1984； Fossen and Tikoff， 1993； Tikoff and Fossen， 1993）和一般三斜變形模式（Jiang and Williams，1998； Lin et al.， 1998）。Ramasy及其后續(xù)變形模型的理論基礎是連續(xù)統(tǒng)內的均勻介質變形（Jiang， 2013），其中任意一物質點的運動學直接反應了連續(xù)統(tǒng)的宏觀流變場特征?；谶@一假設， 實際構造變形分析中研究者通常用從任意尺度構造（如: 露頭構造、手標本構造、顯微構造）獲得的流變場代表宏觀區(qū)域流變場。

然而， 局部流變場與區(qū)域流變場往往并不一致（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）。然而， 迄今為止對于流變場如何在不均勻的變形區(qū)域內分配以及局部分配的流變場與宏觀區(qū)域流變場如何聯系仍然缺乏深入了解。正因如此， 雖然流變場分配的概念早已被構造研究者認同， 并被用來解釋一些區(qū)域內局部的不均勻變形（Wilson and Powell， 2001；Misra et al.， 2009； Shigematsu et al.， 2009； Carreras et al.， 2013； Aydin and de Joussineau， 2014）， 然而在實際研究應變帶的整體運動學邊界條件時卻極少考慮流變場分配因素。

Eshelby （1957， 1959）在連續(xù)變形假設基礎上，建立了均勻彈性基質與嵌入其中的均勻橢球體之間流變場的數學聯系。Eshelby公式被拓展到探討嵌入在線性粘性牛頓體材料（Newtonian）中橢球體的變形（Bilby et al.， 1975； Ellis and Watkinson， 1987）。牛頓體的流變學強度一般用粘性系數表述， 而對于非線性粘性材料（power-law）， 其流變學強度受諸多因素影響， 如流變指數、流變速率等。為了研究非線性粘性材料基質與嵌入的橢球體的變形， 非線性粘性材料的流變學強度用變形過程中的有效粘度（effect viscosity）描述， Lebensohn and Tomé （1993）、Jiang and Bentley （2012）以及Jiang （2013）對非線性粘性材料的有效粘度及其數值算法進行了深入開發(fā)。Eshelby公式奠定了研究宏觀區(qū)域流變場在局部變形體內分配的理論基礎。眾多研究者利用給定的基質流變場代表區(qū)域流變場， 基質中嵌入的橢球體代表局部地質體， 建立了多尺度流變數值模型（Jiang and Bentley， 2012； Jiang， 2013； Xiang and Jiang，2013）。這一模型在局部變形體與區(qū)域運動學邊界條件之間的鴻溝上建立了橋梁。

本文介紹基于Eshelby理論建立3D多尺度數值模型， 模擬研究流變場在不同尺度上的分配。我們將給定的基質流變場對應宏觀尺度上的運動學邊界條件， 基質中可變形的橢球體視為局部地質體， 橢球體內部的有限應變主軸模擬變形組構。這一簡單模型可以模擬三個尺度間的構造變形關系，例如宏觀尺度、露頭尺度以及組構尺度， 因此這一模型是多尺度的。模擬計算揭示， 橢球體相對于基質的流變學強度越低， 其變形越傾向于簡單剪切變形， 并且有限應變的積累越快； 單一橢球體內的流變場都不等同于宏觀流變場， 但全部橢球體內發(fā)育的模擬面理和拉伸線理產狀的總體格局能反映基質流變場的特點。模擬結果揭示天然巖石中流變場分配主要取決于巖石的相對流變強度。模擬結果還表明， 由于流變場的分配， 對自然界應變帶的宏觀流變場和對應的運動學邊界條件研究，利用傳統(tǒng)的以變形材料均勻為前提假設的變形分析方法是不合適的。相對于單一尺度的連續(xù)統(tǒng)變形理論， 基于Eshelby理論的多尺度構造變形分析思路在很大程度上避免了巖石的均一性假設， 更具普遍意義。

1　數值模擬理論背景

1.1坐標系與參考系

數值模型的建立依賴于參考坐標系， 在此我們首先定義下文所涉及的參考坐標系和相關參數。在圖1a定義的右手固定坐標系x1x2x3中， 均勻高應變帶的走向平行于x1軸， 任意給定邊界速度v。邊界速度沿坐標軸方向分解為v1， v2和v3（圖1a）， 相應在應變帶內產生分別平行于x1， x2， x3的拉伸應變率，，和平行于x1， x3的剪應變率和（圖1b）?？臻g內任意一個向量e的方向可用兩個角度θ和φ確定， θ為向量e在x1x2平面上的投影與x1軸正方向的夾角， φ為向量與x3軸正方向的夾角（圖1c）?？臻g內任意一個橢球的方位由三個橢球半徑軸方位確定（三對θ和φ角）， 其中只有三個角度是獨立的。我們取最大橢球半徑軸的θ和φ角和中間半徑軸的一個方位角φ （圖1d）， 當φ = 90°時， φ為橢球中間軸在x1x2平面的投影與x1軸正方向的夾角， 但當φ=90°時， φ為橢球中間半徑軸與x3軸正方向的夾角。

1.2均勻介質的連續(xù)變形機制

圖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）詳盡闡述了等體積、連續(xù)、穩(wěn)態(tài)條件下巖石有限應變的二維一般剪切模式。Sanderson and Marchini （1984）提出了走滑擠壓變形（transpressional zones）的三維單斜分析模式。隨后大量三維分析模式陸續(xù)發(fā)表（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）系統(tǒng)地探討了非線性粘性體有效粘度的近似方法、計算精度以及數值方法。

Eshelby （1957）在數學上證明當均勻基質中的嵌入體為橢球體時， 變形過程中橢球體內部的應力和應變都是處處一致的， 但隨著橢球的幾何形狀與方位的變化， 橢球內部流變場是非穩(wěn)態(tài)的。對于各向同性非線性粘性材料， 橢球體內流變場與遠離嵌入體的基質內流變場滿足如下公式（Jiang， 2013）:

兩式中DE和DM分別是橢球嵌入體內部和遠離嵌入體的基質內的拉伸張量； WE和WM分別是橢球嵌入體內部和遠離嵌入體的基質內的渦度張量。JS是四階對稱單位張量； S和Π分別是對稱和反對稱的Eshelby張量。A是四階應變率分布（Partitioning） 張量， 對各向同性材料其表達式為:

式（3）中n是非線性粘性基質材料的應力指數， reff是橢球嵌入體與基質材料的有效粘度比。由于橢球嵌入體的有效粘度取決于其變形時的真實應力（或者應變率）， 嵌入體材料的有效粘度可通過Mancktelow和Jiang開發(fā)的迭代程序進行數值求解（Mancktelow，2013； Jiang， 2013）。

當給定遠離嵌入體的基質中流變場速度梯度張量LM， 將LM分解為DM和WM， 我們就可以通過公式（2a）， （2b）及公式（3）求解分布于橢球嵌入體內部流變場的DE和WE， 從而求解速度梯度張量LE。在以下模擬計算中我們假設基質內的流變場是均勻且穩(wěn)態(tài)的。此時， LM是恒定的， 并可由內部運動學渦度刻畫。由于橢球嵌入體的形狀和方向是不斷變化的， 嵌入體內部的流變場雖然是均勻的但是非穩(wěn)態(tài)的。由公式（2）確定的渦度張量WE， 包含有剛體旋轉分量。只有內部渦度才能真正刻畫流變場LE的本質， 利用Jiang （2014）開發(fā)的應用程序我們可以計算LE的內部運動學渦度數。顯然， 橢球嵌入體內部流變場的運動學渦度是隨時間演化的。通過求解增量應變， 我們可以求解橢球嵌入體經歷任意時間段的變形后的有限應變量大小、主應變大小以及應變主軸的方位。

2　數值模擬

天然構造變形材料（巖石）組成成分往往十分多樣、結構上具有多尺度特征。從變形巖石整體上（大尺度）看， 巖石由許多形狀各異、流變強度不同的次級尺度元素組成。顯然， 建立與天然巖石對等的理論模型來研究巖石變形行為， 是巖石流變學研究的理想目標。然而， 受限于現有的材料力學理論以及對天然巖石力學性質的了解， 這一目標暫時還難以實現。因此， 對天然巖石構成進行適當簡化建立相對簡單的相似模型就成為必然。目前， 材料力學領域在研究非均勻材料變形行為時， 往往將非均勻材料對等為一個力學性質十分相近的均勻介質， 這個介質被稱為“均勻對等介質（HEM）”（Mura， 1987； Lebensohn et al.， 1998）。借鑒這一思路， 我們在研制巖石整體（大尺度）變形時也將其視為一個均勻介質， 在研究巖石中某一元素（次級尺度）時將其視為嵌入均勻介質中的一個非均勻體。這樣， 我們就可以依據Eshelby理論建立簡單的多尺度數值模型， 來模擬研究巖石的變形行為。

2.1數值模型

我們構建一個由單一橢球體嵌入基質材料的數值模型來模擬宏觀區(qū)域構造與局部地質體。數值模型中， 橢球體與基質都是均勻的各向同性的非線性粘性（power-law）材料， 兩者的有效粘度不同。以同一應變率狀態(tài)為參照， 我們用橢球體材料與基質材料的流變應力指數ne和nm， 以及橢球體與基質材料的有效粘度比（reff）定義橢球體與基質材料的相對流變學性質（Jiang and Bentley， 2012）。橢球體的方向由三個球面角θ， φ和φ確定 （圖1d）， 橢球幾何形狀由橢球半徑比值（a1: a2: a3）確定（a1， a2， a3分別是最大半徑、中間半徑和最小半徑）。

在圖1a坐標系內給定恒定的基質流變場LM。除橢球嵌入體外， 基質的變形可視為均勻介質的穩(wěn)態(tài)變形， LM就等同于等式1所定義的L。我們假設變形是恒體積的， 即。變形過程中跟蹤計算橢球嵌入體內流變場速度梯度張量LE， 從而跟蹤流變場的運動學渦度數、有限應變主軸的方位以及有限應變量大小。顯然， 所有相關的參數都會影響橢球嵌入體內的流變場， 包括:

通過系統(tǒng)地調整上述參數， 運行大量不同橢球體初始形態(tài)和方位模擬計算， 考察基質流變場在橢球體內的分配特征以及影響流變場分配的因素。模擬過程中， 首先系統(tǒng)地考察橢球體在平面應變型式下的基質流變場中發(fā)生平面變形時流變場的分配特征（圖2）， 然后逐步考察一般三斜的基質流變場在橢球體內的分配。本次研究中， 我們主要用橢球體中流變場的運動學渦度數和有限應變量大小來刻畫流變場的分配特征。

圖2　平面應變橢球的空間方位Fig.2 Orientation of the plane straining ellipsoid

2.2橢球體內的流變場特征

圖3和圖4分別揭示了平面應變的橢球嵌入體內流變場的運動學渦度數和有限應變量的演化歷史。在相對粘度低的橢球嵌入體內（reff＜ 1）， 流變場的運動學渦度數經歷了初始階段的起伏變化后很快趨于穩(wěn)定（NM≈ 2）。reff越低， wk越易于達到穩(wěn)定， 且最終穩(wěn)定的wk值越接近于1（圖3）。本次模擬計算設定的橢球體與基質流變強度差異處于自然界巖石強度差異范圍內， 低流變學強度的橢球體的變形接近于簡單剪切。例如， 當基質流變場的wk=0.6， 在有效粘度是基質粘度的五分之一的橢球嵌入體內， 流變場的wk約為0.95。橢球體與基質積累有限應變的速率也有明顯差異（圖4）。盡管基質和橢球嵌入體內有限應變的積累都是非線性的， 但兩者的比值近似線性。有效粘度越低， 橢球嵌入體積累有限應變的速度越快。當有效粘度比并不十分顯著時， 低流變學強度的橢球嵌入體積累的有限應變就可達到基質材料的數倍， 如: 當reff≈0.2時， NE/ NM≈3.5（圖4f）。

平面應變流變場的分配規(guī)律也適合于一般三維流變場。如圖5所示， 不同型式的基質流變場在一任意初始方位的三維橢球體內的分配特點基本相同。橢球體內流變場的運動學渦度數經歷短暫的波動后很快趨于穩(wěn)定， 且穩(wěn)定的渦度數完全取決于橢球體與基質間的相對有效粘度。當給定的基質流變場的渦度數一致（圖5所示， wk=0.6）， 無論基質流變場是平面應變型式（圖5a）或一般單斜型式（圖5c）還是一般三斜型式（圖5e）， 有效粘度比近似相等的橢球體內流變場的最終穩(wěn)定渦度數近似相等。在不同類型基質流變場中， 基質與橢球體內的有限應變量仍然呈現了近似的線性關系， 并且有效粘度越低的橢球體內有限應變積累速率越快（圖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組構與區(qū)域流變場

圖4　平面應變橢球有限應變演化圖Fig.4 Evolution of the finite strain in plane straining ellipsoids

我們假設自然界高應變帶由方向不同、形狀各異的地質體構成。在圖1a所示的參考坐標系下， 將 300個在三維空間內隨機定向的橢球體嵌入到同一均勻基質中， 且假設任意兩個橢球體之間的距離足夠遠， 變形過程中不會產生相互影響。橢球的長半徑與中間半徑的比值在1到10之間隨機變化， 短半徑給定為1。這些橢球能夠覆蓋很大范圍的形態(tài)變化， 從近等軸的球狀（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軸為中心的點集密經大圓環(huán)帶過渡到以x1軸為中心的點集密。其代表的拉伸線理產狀投影從垂直的點集密經過大圓環(huán)帶過渡到水平點集密（圖6b）。

圖6　理論模擬變形組構演變圖Fig.6 Evolution pattern of deformation fabrics from numerical modeling

3　討 論

3.1Eshelby理論在韌性變形研究中的優(yōu)勢

與傳統(tǒng)的韌性變形分析理論相比較， 利用Eshelby理論開展多尺度韌性變形分析的最主要優(yōu)勢在于其假設了天然巖石構成的非均勻性， 能夠合理地解釋自然界流變場的不均勻分配現象。自從Ramsay利用連續(xù)變形理論分析巖石構造變形以來，這一理論在巖石韌性變形分析領域一直占有統(tǒng)治地位。巖石材料均勻是Ramsay理論的基本假設之一。在這一前提下， 自然界構造變形區(qū)域流變場的不均勻分配現象始終難以得到合理的解釋。例如在均勻性假設條件下， 人們通常利用材料彈性失穩(wěn)模式來解釋應變局限化帶的形成機制。而按照這一觀點，高應變帶的尾端相容性問題（Hudleston， 1999）始終難以解決； 盡管Ramsay （1980）提出了兩種解釋模式，但到目前為止未見實際資料的報道。Eshelby理論的初始邊界條件是均勻介質中嵌入可變形的均勻橢球體， 基質和橢球體具有不同的流變強度。這就強調了介質與嵌入其中的橢球體間的不均勻特征。隨著均勻對等介質（HEM）思路的提出以及可行的數值求解方法的開發(fā)， 基于Eshelby理論的模型可以視為由流變強度不同的橢球體構成（Jiang， 2014）； 研究其中一個橢球體內的流變場時， 其他橢球體整體視為基質， 基質的流變性質用實時求解的HEM的流變性質近似（Jiang and Bentley， 2012）。顯然， 通過改進的Eshelby理論模型進一步強調了變形材料的不均勻性， 與天然巖石的構成也就更加相似。上述我們利用簡單的Eshelby理論模型的模擬計算， 揭示了天然巖石變形流變場的不均勻分配主要受巖石相對流變強度制約， 這一規(guī)律與已有的觀察結論一致。

與其他數值模擬方法相比， 基于Eshelby理論的數值模擬的優(yōu)勢在于大應變和多對象同時研究。目前用于變形模擬研究的常用數值模擬方法主要有有限元法和有限差分法， 兩種方法都有成熟的商業(yè)軟件， 也都應用于模擬巖石的韌性變形（Mancktelow，2011）。兩種方法的一個共同特點是將模型網格化，通過求解網格節(jié)點的運動學和力學參數來描述變形過程。通過對模型中不同區(qū)域內的網格設定不同的力學參數， 也可以建立由不均勻材料構成的數值模型。模擬過程中， 當任意相鄰的兩個網格交差產生幾何破壞時， 模型材料變形的應力、應變連續(xù)性被破壞， 運算就無法繼續(xù)進行。已有的模擬研究經驗表明， 網格密度越大， 材料的力學性質刻畫得就越精確， 計算精度也就越高， 但此時運算速率也越慢且網格越容易發(fā)生幾何破壞， 模型材料經歷的有限應變量就越小。因此， 無法實現模型與天然巖石應變相當的大應變是基于有限元理論和有限差分理論的模擬方法的主要困難之一， 特別是刻畫具有較大流變學強度差異的非均勻材料模型尤為如此。理論上， 基于有限元理論和有限差分理論的數值模型可以無限包含非均勻元素， 但實際模擬實驗過程中，這類模型能夠發(fā)生的有限應變量有限而且運算量也無法由一般硬件設備承擔。因此， 已有的基于這兩種數學理論在巖石變形方面的研究模型材料都較為單一， 往往僅包含一個或有限的幾個非均勻元素（Mancktelow， 2011）。已有研究表明（Jiang and Bentley，2012； Jiang， 2013； Xiang and Jiang， 2013）， 基于Eshelby理論的數值模型可以開展大應變模擬計算，并且可以創(chuàng)建由大量流變學性質不同的元素構成的數值模型。

3.2Eshelby數值模型在巖石變形中的局限性

同傳統(tǒng)的巖石韌性變形分析理論一樣， Eshelby理論也以連續(xù)變形假設為前提， 這一假設對于天然巖石變形可能過于嚴格。天然巖石在構造變形過程中， 不同塊體（或元素）之間不僅發(fā)生連續(xù)變形， 還可能會發(fā)生相對滑移的非連續(xù)變形。顯然， 如果巖石變形過程中非連續(xù)變形十分顯著， 用Eshelby理論建立的連續(xù)變形模式來研究其變形行為就不合適了。然而， 在考慮引入新的理論建立包含非連續(xù)變形基質的模型之前， 我們必須對巖石變形行為開展系統(tǒng)的研究， 揭示巖石發(fā)生非連續(xù)變形的標志， 非連續(xù)變形對巖石整體變形的貢獻大小等等。如果包含非連續(xù)變形的理論模型是研究天然巖石變形的最終理想模型， 那么基于Eshelby理論的數值模型只是人們邁向最終目標的一個中間過程。

Eshelby理論對嵌入基質中的次級元素的假設條件可能過于嚴格。一方面， 地質體的初始幾何形態(tài)可能與標準橢球體有一定出入， 這樣在實際變形過程中其內部的流變場并非嚴格的處處均一。另一方面， 利用均勻橢球體模擬地質體也是簡化地近似；天然地質體的流變性質往往并不完全均勻， 變形過程中， 流變學性質的不均勻性也會導致實際流變場比理想的模擬流變場復雜得多。因此， 實際天然變形巖石中記錄流變場特征的組構可能比數值模擬預測結果更為復雜， 數值模擬結果更多地揭示了巖石的總體變形行為（Xiang and Jiang， 2013）。此外， 天然地質體還可能具有各向異性的流變學性質。盡管對一些相對簡單的各向異性體的變形已開展了一些初步的模擬研究工作， 如云母殘斑晶的變形和定向（Chen et al.， 2014）， 但是多數情況下地質體的各向異性特征并不被人們所掌握， 對構建具有各向異性特征的數值模型缺乏基本的約束。因此， 地質體的各向異性特征也會影響地質體內部的流變場以及變形組構發(fā)育。

3.3韌性變形分析思路

本次模擬揭示， 構造變形區(qū)域內分配在不同流變學強度的局部巖石內的流變場往往不同。因此， 傳統(tǒng)的以均勻性假設為前提的變形分析方法就難以勝任對流變場不均勻分配區(qū)域的宏觀流變場特征的研究了。首先， 以往的韌性構造變形分析中， 拉伸線理方向或優(yōu)選方位直接代表宏觀流變場的剪切方向。本次模擬研究表明， 變形區(qū)域內拉伸線理產狀的總體分布格局是由宏觀流變場形式決定的， 隨著流變場特征的變化， 拉伸線理優(yōu)選方位可以從與剪切方向一致漸變?yōu)榕c剪切方向垂直。模擬實驗表明， 可以利用變形區(qū)域內拉伸線理產狀的總體分布格局來揭示宏觀流變場形式。其次， 傳統(tǒng)研究方法中通過運動學渦度測量來直接指示區(qū)域上的變形程度和運動學邊界條件， 在流變場不均勻分配的區(qū)域就不合適了。本次模擬研究表明， 局部巖石內流變場的運動學渦度數不僅由宏觀流變場的運動學渦度決定， 而且還受局部巖石的相對流變學強度影響。對于如何利用從不同局部變形巖石中測量的運動學渦度數來揭示宏觀區(qū)域流變場的運動學渦度數， 從而揭示宏觀運動學邊界條件， 還需要進一步研究， 我們將在后續(xù)研究中進一步探討。

4　結 論

本次工作系統(tǒng)地綜述了基于Eshelby理論的多尺度數值模擬研究思路和方法。通過多尺度數值模擬， 我們對天然巖石的構造變形研究取得了兩點初步認識: （1）構造變形區(qū)域， 流變場的分配主要受巖石的相對流變學強度影響， 流變強度越低的巖石內流變場的剪應變越集中且有限應變累積速度越快；（2）在流變場不均勻分配的區(qū)域開展構造變形分析，傳統(tǒng)的以變形巖石均勻性假設為前提的韌性變形分析方法難以勝任； 變形區(qū)域上宏觀流變場可以通過面理和拉伸線理組構產狀的總體分布格局并結合可能的變形形式來推斷； 構造變形區(qū)域上流變場的運動學渦度數無法從測量局部變形巖石中的運動學渦度數直接獲取。

致謝： 兩位審稿專家對本文的最終成稿提出了科學嚴謹的修改意見， 讓筆者獲益匪淺， 在此表示由衷的感謝。

參考文獻（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