馬連華, 李雅嬈, 譚如峰

(1.河北大學(xué)質(zhì)量技術(shù)監(jiān)督學(xué)院，保定 071002；2.河北大學(xué)建筑工程學(xué)院，保定 071002；3.河北大學(xué)人事處，保定 071002)

彈性或黏性基體中橢球夾雜問題的研究已經(jīng)有很長的歷史[1-3]。針對位于地球巖石圈中的高應(yīng)變巖石流變場分配問題，中外學(xué)者采用Eshelby細觀力學(xué)方法[3-4]進行了充分的研究。Jeffery[5]首先推導(dǎo)了遠場均勻流動的各向同性牛頓黏性材料中嵌入剛性橢球的角速度方程。Bilby等[6]將Eshelby細觀力學(xué)理論擴展應(yīng)用到橢球夾雜和宏觀區(qū)域構(gòu)造具有不同黏性的流變體情況。Freeman[7]基于Jeffery的理論對各種黏性流的剛性橢球夾雜的運動進行了數(shù)值研究，并且將Eshelby方程的數(shù)值解與Jeffery方程的數(shù)值解進行了比較[8]。Je?ek[9]開發(fā)了FORTRAN程序，用來計算給定流變場中單個剛性橢球夾雜的旋轉(zhuǎn)路徑。Jiang[10-11]基于Eshelby理論并利用Mathcad計算了單個以及多個剛性和可變形橢球夾雜的演化進程。Jiang[12-13]通過引入近似算法來求解夾雜運動微分方程，提高了計算效率。向必偉等[14]從固體連續(xù)變形機制入手，重點介紹基于Eshelby理論的多尺度數(shù)值模擬思路和方法，探討流變場分配問題，揭示了嵌入基質(zhì)中的橢球體內(nèi)分布的流變場主要取決于橢球體與基質(zhì)間的相對流變強度，橢球體的相對流變強度越低，其變形越接近于簡單剪切，且有限應(yīng)變積累越快。同時還揭示，不同流變強度的橢球體內(nèi)模擬拉伸線理和面理產(chǎn)狀的總體格局反映基質(zhì)流變場特征。劉穩(wěn)航等[15]介紹了韌性剪切帶中橢球狀標志體變形的最新研究進展，探討了基于Jeffery理論和Eshelby理論的數(shù)值模擬思路和方法，并利用Mathcad軟件模擬了給定條件下的橢球夾雜的運動軌跡、變形特征。Qu等[16]基于MATLAB模擬了流變巖石中橢球夾雜的運動規(guī)律，但并未考慮橢球夾雜在有限體分比下的運動演化機制。高麗敏等[17]基于Eshelby理論,通過數(shù)值模擬,將經(jīng)過兩期變形疊加產(chǎn)生的模擬組構(gòu)產(chǎn)狀格局與只經(jīng)歷后期變形產(chǎn)生的模擬組構(gòu)產(chǎn)狀格局進行對比,揭示疊加變形組構(gòu)的產(chǎn)狀分布特征。林少澤等[18]利用適用于黏性材料的Eshelby模型，計算了橢球體內(nèi)部流變場最大主應(yīng)變軸和最小主應(yīng)變軸方位并進行赤平投影，用以模擬線理和面理極點產(chǎn)狀。

1 均勻介質(zhì)的連續(xù)變形機制

均勻介質(zhì)的流變學(xué)在巖石流變學(xué)方面的應(yīng)用開始于簡單的二維變形的研究。Jiang等[19]提出了一般的三維變形模式，推導(dǎo)出了速度梯度張量和位移梯度張量。

根據(jù)Jeffery[5]的定義可以得到，a1、a2、a3為橢球狀局部地質(zhì)體3個應(yīng)變主軸的半軸，3個半軸應(yīng)滿足：a1≥a2≥a3。為了方便進行數(shù)值模擬，需要對內(nèi)部和外部建立兩套不同的坐標系，其中內(nèi)部的參考坐標系R′是始終與橢球夾雜的3個主軸a1、a2、a3保持平行，長度用相對應(yīng)的單位基本矢量η1、η2、η3表示；而外部坐標系R使用x1、x2、x3表示，使用相對應(yīng)的單位基本矢量e1、e2、e3進行表示。

剛性橢球狀夾雜在基體中的旋轉(zhuǎn)相當于其內(nèi)部坐標系R′和外部坐標系R不停發(fā)生變化，同時每個旋轉(zhuǎn)位置都可以由外部坐標系R中的單位矢量ei來表示，另外也可以由內(nèi)部坐標R′的單位矢量ηi來表示。在x1x2x3的坐標系統(tǒng)中，xi軸和x′i軸之間存在的關(guān)系為[15]

x′=Qx，x=QTx′

(1)

式(1)中：Q是從xi坐標系到xi′坐標系的過渡矩陣；QT是Q的轉(zhuǎn)置矩陣。同時，也可以定義為

e′i=Qei，ei=QTe′i

(2)

均勻介質(zhì)的流動場可以使用歐拉速度梯度張量可以表示為

(3)

(4)

(5)

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

L=D+W

(6)

式(6)中：

(7)

(8)

根據(jù)坐標轉(zhuǎn)換，有：

D′=QDQT

(9)

(10)

只要給定流變場的速度梯度張量L，就可以求解任一時刻流變場的運動學(xué)渦度，并可以求解經(jīng)歷任意時間段的變形后的有限應(yīng)變量大小、主應(yīng)變大小以及應(yīng)變主軸的方位。

2 等效遠場法

為了應(yīng)用經(jīng)典的、適用于無限大基體的Eshelby等效夾雜理論，將含橢球夾雜的代表體積單元(RVE)置于無限大的基體中并且承受遠場速度梯度L∞邊界條件，假設(shè)在遠場的邊界處施加的速度梯度張量L∞在RVE體積V內(nèi)產(chǎn)生的流變場與直接在RVE邊界處施加的速度梯度張量L所產(chǎn)生的流變場是一致的，即在RVE邊界處施加L與在遠場處施加的L∞是等效的，該種方法就是等效遠場法[20]。同時，速度梯度張量L(L∞)可以分解拉伸張量D(D∞)和渦度張量W(W∞)，如圖1所示。

圖1 旋轉(zhuǎn)的等效遠場模型Fig.1 Rotational equivalent far field model

當橢球體為剛性夾雜時，僅發(fā)生剛性旋轉(zhuǎn)變形，其拉伸張量為0。

根據(jù)向必偉等[14]給出的橢球體內(nèi)流變場與基體內(nèi)的流變場的關(guān)系，可得到橢球夾雜的渦度張量WE為

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

(11)

式(11)中：Π為反對稱的Eshelby張量；S-1為對稱的Eshelby張量的逆；WE和W∞分別為橢球體夾雜和RVE邊界的渦度張量；D∞為RVE邊界處的拉伸張量。

根據(jù)等效遠場法，可分別建立拉伸和渦度張量與遠場拉伸與渦度張量的關(guān)系：

D=(1-f)D∞

(12)

W=(1-f)W∞+fWE

(13)

式中:f為橢球夾雜所占的體積分數(shù)。

根據(jù)上述公式，可得到：

(14)

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

(15)

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

WE=W-Π:S-1:D

(16)

由上式可以看出，橢球夾雜的渦度張量與其體分比無關(guān)。

剛性橢球夾雜的旋轉(zhuǎn)演化方程可以表示為[13]

(17)

式(17)中：Q為方位角矩陣；Θ為橢球體夾雜的角速度張量。

為提高計算效率，采用Basar等[21]給出的近似方法來求解上述微分方程：

(18)

(19)

對于可變形橢球夾雜的變形，其形狀隨著時間變化的關(guān)系式為[13]

(20)

3 數(shù)值模擬

基于上述細觀力學(xué)模型，采用MATLAB數(shù)值軟件，模擬了剛性橢球夾雜的演化進程。以文獻[16]中的簡單剪切流變場中給出的參數(shù)為例，研究剛性橢球夾雜的演化進程。采用的計算數(shù)據(jù)如下。

(1)橢球夾雜的3個主半軸分別為：a1=5，a1=3，a1=1。

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

(4)每個步驟的時間增量δt=0.01。

(5)模擬運行的總步數(shù)為1 000，每隔20步取一次結(jié)果。

經(jīng)過計算，可以得到剛性橢球夾雜的情況。剛性橢球體夾雜的演化進程如圖2、圖3所示。

“x”是起始點，“*”是終止點圖2 單個剛體橢球夾雜Fig.2 Single rigid body ellipsoid inclusion

Step表示不同的時間步剛性橢球體夾雜的演化圖3 剛體橢球夾雜的演化進程Fig.3 Evolution of rigid body ellipsoid inclusion

4 結(jié)論

基于等效遠場法研究了黏性基體中剛性橢球夾雜的運動演化情況，得到如下結(jié)論。

(1)利用等效遠場法建立了黏性基質(zhì)中剛性橢球夾雜的細觀力學(xué)理論，并基于MATLAB發(fā)展了數(shù)值計算程序。

(2)在給定速度梯度張量條件下，嵌入黏性材料的剛性橢球夾雜運動演化與體分比無關(guān)。

(3)基于發(fā)展的MATLAB程序，模擬了嵌入黏性基體中剛性橢球夾雜的運動情況，揭示了橢球夾雜的演化機制。