三維模擬高溫高靜水壓下礦物顆粒邊界應(yīng)力分布

陳玉香 ，杜建國，劉紅，張建平

(1. 中國地震局 地質(zhì)研究所，北京，100029；2. 中國地震局 地震預(yù)測研究所，北京，100036；3. 湘潭大學(xué) 機(jī)械工程學(xué)院，湖南 湘潭，411105)

探明巖石、礦物的力學(xué)性質(zhì)與變形機(jī)制，有助于了解地球內(nèi)部動(dòng)力學(xué)過程[1]。巖石、礦物的力學(xué)性質(zhì)與變形機(jī)制取決于礦物成分、結(jié)構(gòu)、礦物顆粒邊界效應(yīng)和溫壓條件等多種因素。Holtzman等[2?4]沿用非靜水壓下巖石、礦物的應(yīng)變率與差應(yīng)力具有冪律關(guān)系的流變經(jīng)驗(yàn)公式，并考慮顆粒粒徑[5]、氧逸度[6]、流體[7?8]、活化體積[9]及其綜合因素[10]的影響模擬研究巖石、礦物的變形機(jī)制。隨著計(jì)算機(jī)技術(shù)的迅速發(fā)展，越來越多的研究者用數(shù)值模擬方法宏觀模擬研究巖石的力學(xué)性質(zhì)尤其是對(duì)巖石裂紋體的斷裂機(jī)制[11]，并應(yīng)用于巖土工程中。近年來，楊振濤等[12?13]采用數(shù)值模擬單礦物與雙礦物巖石在差應(yīng)力下的流變性質(zhì)。然而，人們還只是借助于數(shù)值模擬方法探討礦物集合體中顆粒邊界移動(dòng)的結(jié)果[14]、邊界濕化對(duì)礦物顆粒集合體剪應(yīng)力分布的影響[15]。事實(shí)上，礦物顆粒邊界是具有大離子半徑的元素的重要集聚地[16]，在一定的溫度、壓力條件下，礦物顆粒邊界能形成反應(yīng)邊構(gòu)造，生成新的礦物[17]或使顆粒邊界濕化[18]，從而影響巖石的變形機(jī)制以及彈性波速、電導(dǎo)率等[19]。礦物顆粒邊界效應(yīng)作為影響巖石、礦物力學(xué)性質(zhì)與變形機(jī)制的重要因素之一，人們對(duì)其研究較少。為此，本文作者應(yīng)用ANSYS有限元軟件三維模擬計(jì)算在高溫800 ℃、高靜水壓1.6 GPa條件下，上地幔、下地殼重要組成礦物鎂橄欖石與透輝石集合體在加壓瞬間完全彈性變形階段顆粒邊界層的應(yīng)力場，進(jìn)而分析探討靜水壓下2種不同礦物顆粒邊界層的應(yīng)力分布規(guī)律。

1 研究方法

1.1 計(jì)算模型

參照桑祖南等[20]對(duì)單斜輝石顆粒尺寸(長軸長度平均為0.224 mm，短軸長度平均為0.096 mm)的描述。鎂橄欖石、透輝石顆粒采用立方體結(jié)構(gòu)，立方體邊長為 0.2 mm；顆粒邊界為 2種礦物的過渡層，厚度為5 μm。

由于礦物顆粒為立方體結(jié)構(gòu)，不能將其簡化為平面應(yīng)變或帶厚度的平面應(yīng)力問題，從而直接建立三維模型。礦物顆粒排列從下至上依次為鎂橄欖石顆粒、透輝石顆粒、鎂橄欖石顆粒。根據(jù)對(duì)稱性，建立笛卡兒坐標(biāo)系，對(duì)稱中心為坐標(biāo)原點(diǎn)O，取整體模型的1/8(稱為計(jì)算模型)進(jìn)行計(jì)算。圖1所示為計(jì)算模型幾何圖，上層為鎂橄欖石，下層為透輝石，中間薄層為顆粒邊界。采用ANSYS內(nèi)置的Solid65單元類型、規(guī)則網(wǎng)格(在X軸和Z軸方向劃分20等分，在Y軸上的透輝石層劃分20等分，顆粒邊界層劃分20等分，鎂橄欖石層劃分40等分)，計(jì)算模型中共有32 000個(gè)規(guī)則單元、35 721個(gè)節(jié)點(diǎn)和107 163個(gè)自由度。

圖1 計(jì)算模型幾何圖Fig.1 Geometric plot of computed model

模型中礦物顆粒層、礦物顆粒邊界層都為各向同性彈性材料；礦物顆粒與顆粒邊界層之間接觸良好，既不考慮空隙，也不考慮兩礦物顆粒之間的摩擦。對(duì)稱中心O點(diǎn)固定，XOZ面上Y向位移為0，XOY面上Z向位移為0，YOZ面上X向位移為0；各外表面加載1.6 GPa均布?jí)毫挫o水壓力。

1.2 計(jì)算參數(shù)

在計(jì)算過程中所需要的彈性模量E、泊松比ν計(jì)算式如下：

其中：Ks為體積模量；G為剪切模量。

在溫度為800 ℃、壓力為1.6 GPa 條件下，鎂橄欖石的體積模量與剪切模量從 Anderson等[21]得出的高溫下地幔礦物的彈性常量表中直接查得；而透輝石的體積模量與剪切模量從 Bass[22]整理的立方晶體礦物的彈性模量表中查到，在室溫、常壓下，它們分別為114 GPa和64.9 GPa。計(jì)算模型中透輝石、鎂橄欖石的彈性參數(shù)見表1。

顆粒邊界層為透輝石與鎂橄欖石這2種礦物的過渡層，其組成物質(zhì)的體積模量Ks3、剪切模量G3計(jì)算式為：

表1 在800 ℃和1.6 GPa條件下透輝石和鎂橄欖石的彈性參數(shù)Table 1 Elastic parameters of diopside and forsterite at 800 ℃ and 1.6 GPa

其中：Ks1為透輝石體積模量；Ks2為鎂橄欖石體積模量；G1為透輝石剪切模量；G2為鎂橄欖石剪切模量；a為顆粒邊界層組分礦物透輝石的體積模量系數(shù)；b為顆粒邊界層組分礦物透輝石的剪切模量系數(shù)。

由于本文是計(jì)算靜水壓為1.6 GPa時(shí)顆粒邊界層的應(yīng)力，重點(diǎn)考慮體積模量的影響，故選取b為固定值0.5，a從0至1之間變化，根據(jù)式(1)～(4)計(jì)算得到相應(yīng)的顆粒邊界層的彈性模量E與泊松比ν，見表2。

表2 顆粒邊界層的彈性參數(shù)、平均等效應(yīng)力及3點(diǎn)等效應(yīng)力Table 2 Elastic parameters, equivalent stresses of average and three points at grain boundary

1.3 應(yīng)力等值線計(jì)算結(jié)果

由ANSYS軟件計(jì)算的應(yīng)力解均采用節(jié)點(diǎn)解，應(yīng)力單位為Pa。在800 ℃和1.6 GPa靜水壓下，當(dāng)b為0.5，a為0.5時(shí)，等效應(yīng)力等值線圖且根據(jù)對(duì)稱性擴(kuò)展至整個(gè)模型的等效應(yīng)力場如圖2所示。從圖2可知：在靜水壓下，透輝石與鎂橄欖石礦物顆粒之間的界面上存在應(yīng)力集中現(xiàn)象，顆粒邊界層的應(yīng)力分布不均勻。根據(jù)表2 ，依次計(jì)算顆粒邊界層具有不同彈性模量和泊松比時(shí)的節(jié)點(diǎn)等效應(yīng)力等值線圖。圖3所示是a分別為0.1和0.5時(shí)由計(jì)算模型擴(kuò)展至整個(gè)模型的邊界層等效應(yīng)力等值線圖。

1.4 結(jié)果分析

從采用的節(jié)點(diǎn)解的等效應(yīng)力等值線圖可知：礦物顆粒邊界層的彈性模量和泊松比不同，其應(yīng)力分布也不同。為了找出其中的變化規(guī)律，取礦物顆粒邊界層2 402個(gè)節(jié)點(diǎn)，計(jì)算其平均等效應(yīng)力，并在礦物顆粒邊界層取A，B和C3點(diǎn)，此3點(diǎn)在計(jì)算模型中的坐標(biāo)依次為(0.085 0，0.105 0，0.100 0)，(0.085 0，0.100 0，0.100 0)和(0.085 0，0.102 5，0.100 0)，單位為 mm。隨著體積模量系數(shù)a及其相應(yīng)的彈性模量和泊松比的變化，計(jì)算得到的平均等效應(yīng)力以及A，B和C3點(diǎn)的等效應(yīng)力部分解見表2。

圖2 整體模型等效應(yīng)力等值線圖Fig.2 Equivalent stress contour plot of whole model when a=0.5 and b=0.5

圖3 礦物顆粒邊界等效應(yīng)力等值線圖Fig.3 Equivalent stress contour plot at mineral grain boundaries

據(jù)表2中的數(shù)據(jù)做出散點(diǎn)圖并進(jìn)行高斯曲線擬合(見圖4(a))，平均等效應(yīng)力以及A，B和C3點(diǎn)的等效應(yīng)力σe與體積模量系數(shù)a的關(guān)系滿足高斯分布函數(shù)方程：

擬合的高斯分布函數(shù)的相關(guān)系數(shù)R和方程特征參數(shù)如下：

σAe與a的關(guān)系中，R2=0.999 77，w=1.41，A=?289.40，y0=228.12，x0=0.36；

σBe與a的關(guān)系中，R2=0.999 57，w=1.02，A=?133.80，y0=164.65，x0=0.73；

σCe與a的關(guān)系中，R2=0.998 23，w=0.90，A=?74.02，y0=121.91，x0=0.55；

σAve與a的關(guān)系中，R2=0.998 22，w=0.78，A=?63.25，y0=105.41，x0=0.51。

將平均等效應(yīng)力以及A，B和C3點(diǎn)的等效應(yīng)力與顆粒邊界層彈性模量的散點(diǎn)圖進(jìn)行擬合，擬合結(jié)果也符合高斯分布，如圖4(b)所示。擬合方程的相關(guān)系數(shù)R和特征參數(shù)如下：

σAe與E的關(guān)系中，R2=0.999 81，w=0.04，A=?4.99，y0=173.30，x0=1.56；

σBe與E的關(guān)系中，R2=0.975 10，w=0.02，A=1.32，y0=63.26，x0=1.58；

圖4 顆粒邊界平均等效應(yīng)力及A，B和C 3點(diǎn)等效應(yīng)力曲線圖Fig.4 Equivalent stress curves of point A, B, C and average at grain boundary

σCe與E的關(guān)系中，R2=0.999 91，w=0.04，A=?3.98，y0=145.27，x0=1.56；

σAve與E的關(guān)系中，R2=0.999 27，w=0.03，A=?2.23，y0=108.10，x0=1.56。

2 討論

巖石、礦物若長時(shí)間處于高溫高壓并有差應(yīng)力下會(huì)發(fā)生蠕變變形，但對(duì)于某一特定的礦物或巖石，在不同溫度和壓力下，人們還不了解應(yīng)變率與應(yīng)力的確切關(guān)系式，如橄欖石在高溫下的蠕變機(jī)制在巖石圈溫度 900～1 200 ℃范圍內(nèi)與溫度高于 1 200 ℃時(shí)不一樣[23]，后者應(yīng)變率與差應(yīng)力呈冪律關(guān)系，而前者卻呈指數(shù)關(guān)系。在一定溫度和壓力下，組成巖石的不同礦物顆粒之間可能發(fā)生化學(xué)反應(yīng)或部分熔融，在礦物顆粒接觸界面上生成新的物質(zhì)，進(jìn)而形成薄薄的礦物顆粒邊界層或直接由不同礦物以某一比例組成的過渡層。由于邊界層物質(zhì)可能不同程度地含有孔隙[24]甚至熔體[25]等，即使是在相同的溫度和壓力下，其邊界層的彈性參數(shù)也可能不一樣[26]，從而影響巖石、礦物的力學(xué)性質(zhì)與變形機(jī)制。

鎂橄欖石、透輝石的熔點(diǎn)都高于1 000 ℃，本研究中采用800 ℃和1.6 GPa靜水壓條件，在上地幔—下地殼溫度、壓力條件范圍內(nèi)，兩礦物集合體在加靜水壓瞬間必然會(huì)發(fā)生彈性變形。通過三維模擬計(jì)算，在完全彈性變形階段(如圖2 所示)，兩礦物顆粒邊界層有應(yīng)力集中現(xiàn)象。為了找出其邊界層應(yīng)力分布規(guī)律，設(shè)定顆粒邊界層為兩礦物的過渡層，邊界層(過渡層)的體積模量為兩礦物的體積模量乘以1個(gè)變化系數(shù)后之和，依此計(jì)算得到邊界層不同的泊松比、彈性模量，進(jìn)行多次模擬計(jì)算，部分結(jié)果見表2。圖4中擬合的高斯曲線方程，相關(guān)系數(shù)R2都在0.95以上。在靜水壓力作用下，當(dāng)立方體結(jié)構(gòu)的鎂橄欖石與透輝石集合體在其顆粒內(nèi)部與顆粒邊界都處于完全彈性變形時(shí)，顆粒邊界層的平均等效應(yīng)力及其空間任意點(diǎn)的等效應(yīng)力與顆粒邊界層組分礦物的體積模量系數(shù)、邊界層組成物質(zhì)的彈性模量之間具有高斯分布規(guī)律，并且不同位置點(diǎn)擬合的高斯分布函數(shù)特征參數(shù)不同，說明在靜水壓下，處于彈性變形階段的礦物顆粒邊界層內(nèi)部應(yīng)力分布不均勻。

不同礦物顆粒邊界的應(yīng)力分布規(guī)律有助于從力學(xué)上解釋巖石、礦物在高溫高壓下的一些微觀變形特征，如位錯(cuò)、擴(kuò)散、邊界遷移、滑移等以及由礦物顆粒邊界引起巖石、礦物的彈性波速的變化。巖石中的顆粒邊界對(duì)應(yīng)力敏感而導(dǎo)致在鉆孔周圍存在彈性波速與應(yīng)力之間的非線性變化[27]；靜水壓試驗(yàn)測出P 波速度隨著主應(yīng)力(靜水壓)的增大呈非線性上升[28]；Ji等[29]通過對(duì)大別—蘇魯超高壓變質(zhì)巖層的研究，得出P波速度?壓力非線性曲線可以用1個(gè)有4個(gè)參數(shù)的指數(shù)方程進(jìn)行描述，并賦予了這4個(gè)參數(shù)的物理意義。盡管本文的模擬計(jì)算是在理想狀態(tài)下進(jìn)行的，只是模擬計(jì)算靜水壓力下立方體結(jié)構(gòu)的鎂橄欖石與透輝石集合體在完全彈性變形時(shí)的邊界應(yīng)力，但可以類推為靜水壓下2種立方體結(jié)構(gòu)的不同礦物的集合體。計(jì)算結(jié)果表明：在礦物顆粒邊界應(yīng)力集中，邊界層應(yīng)力分布不均勻，顆粒邊界層的等效應(yīng)力與邊界層組成物質(zhì)的彈性模量之間具有高斯分布函數(shù)關(guān)系。這對(duì)于具有確切物質(zhì)組成的礦物顆粒邊界，根據(jù)其擬合的高斯分布函數(shù)進(jìn)行化簡，可以得到與彈性模量相對(duì)應(yīng)的等效應(yīng)力，再根據(jù)彈性模量與彈性波速的計(jì)算公式[30]，可以推導(dǎo)出靜水壓力下巖石、礦物內(nèi)部的等效應(yīng)力與彈性波速的關(guān)系式。

在高溫高靜水壓下，立方體結(jié)構(gòu)的2種不同礦物顆粒邊界應(yīng)力的這種高斯分布規(guī)律為礦物顆粒邊界效應(yīng)的研究提供了依據(jù)。在今后的研究中，有必要考慮非靜水壓條件、礦物顆粒內(nèi)部及邊界有蠕變變形、礦物顆粒之間的接觸有摩擦和間隙等情況，并結(jié)合高溫高壓實(shí)驗(yàn)、分子動(dòng)力學(xué)模擬精確確定計(jì)算模型中的彈性模量、泊松比等參數(shù)，進(jìn)而綜合模擬研究礦物顆粒邊界效應(yīng)對(duì)巖石的力學(xué)性質(zhì)與變形機(jī)制的影響。

3 結(jié)論

(1) 在高溫高靜水壓下完全彈性變形時(shí)，2種立方體結(jié)構(gòu)的不同礦物顆粒邊界具有應(yīng)力集中現(xiàn)象，顆粒邊界層內(nèi)部應(yīng)力分布不均勻。

(2) 2種立方體結(jié)構(gòu)的不同礦物集合體在高溫高靜水壓下完全彈性變形時(shí)，過渡邊界層的等效應(yīng)力與其組分礦物的體積模量系數(shù)及其組成物質(zhì)的彈性模量之間具有高斯分布規(guī)律。

致謝：本文得到了湘潭大學(xué)機(jī)械學(xué)院龔曙光教授、中國科學(xué)院地質(zhì)與地球物理所伍向陽研究員、中國地震局地質(zhì)研究所何昌榮研究員、湘潭大學(xué)材料與光電物理學(xué)院毛衛(wèi)國副教授的指導(dǎo)，在此一并表示衷心感謝！

[1] Holtzman B K, Kohlstedt D L, Zimmerman M E, et al. Melt segregation and strain partitioning: implications for seismic anisotropy and mantle flow[J]. Science, 2003, 301(5637):1227?1230.

[2] Holtzman B K, Kohlstedt D L, Morgan J P. Viscous energy dissipation and strain partitioning in partially molten rocks[J]. J Petrol, 2005, 46(12): 2569?2592.

[3] JI Shao-cheng, ZHAO Ping-lao, XIA Bin. Flow laws of multiphase materials and rocks from end-member flow laws[J].Tectonophysics, 2003, 370(1/4): 129?145.

[4] HE Chang-rong. ZHOU Yong-sheng, SANG Zu-nan. An experimental study on semi-brittle and plastic rheology of Panzhihua gabbro[J]. Science in China: Series D, 2003, 46(7):730?742.

[5] Dimanov A, Dresen G, Wirth R. High-temperature creep of partially molten plagioclase aggregates[J]. J Geophys Res, 1998,103(B5): 9651?9664.

[6] Jaoul O, Raterron P. High-temperature deformation of diopside crystal 3. influences ofpO2and SiO2precipitation[J]. J Geophys Res, 1994, 99(B5): 9423?9439.

[7] MEI Sheng-hua, Kohlstedt D L. Influence of water on plastic deformation of olivine aggregates 1. diffusion creep regime[J]. J Geophys Res, 2000, 105(B9): 21457?21469.

[8] MEI Sheng-hua, Kohlstedt D L. Influence of water on plastic deformation of olivine aggregates 2. dislocation creep regime[J].J Geophys Res, 2000, 105(B9): 21471?21481.

[9] Hier-Majumder S, Anderson I M, Kohlstedt D L. Influence of protons on Fe-Mg interdiffusion in olivine[J]. J Geophys Res,2005, 110, B02202, doi: 10.1029/2004JB003292.

[10] Hier-Majumder S, Mei Sheng-hua, Kohlstedt D L. Water weakening of clinopyroxenite in diffusion creep[J]. J Geophys Res, 2005, 110, B07406, doi: 10.1029/2004JB003414.

[11] 趙延林, 曹平, 文有道, 等. 滲透壓作用下壓剪巖石裂紋損傷斷裂機(jī)制[J]. 中南大學(xué)學(xué)報(bào): 自然科學(xué)版, 2008, 39(4):838?845.ZHAO Yan-lin, CAO Ping, WEN You-dao, et al. Damage fracture failure mechanism of compressive-shear rock cracks under seepage pressure[J]. Journal of Central South University:Science and Technology, 2008, 39(4): 838?845.

[12] 楊振濤, 寧杰遠(yuǎn), 趙永紅. 多礦物巖石流變性的有限元計(jì)算[J]. 巖石學(xué)報(bào), 2004, 20(6): 1469?1476.YANG Zhen-tao, NING Jie-yuan, ZHAO Yong-hong. Finite element computation of rheology for multi-mineral rocks[J].Acta Petrologica Sinica, 2004, 20(6): 1469?1476.

[13] Madi K, Forest S, Cordier P, et al. Numerical study of creep in two-phase aggregates with a large rheology contrast:implications for the lower mantle[J]. Earth Planet Sci Lett, 2005,237(1/2): 223?238.

[14] Becker J K, Bons P D, Jessell M W. A new front-tracking method to model anisotropic grain and phase boundary motion in rocks[J]. Comput Geosci, 2008, 34(3): 201?212.

[15] Hier-Majumder S, Leo P H, Kohlstedt D L. On grain boundary wetting during deformation[J]. Acta Mater, 2004, 52(12):3425?3433.

[16] Hiraga T, Anderson I M, Kohlstedt D L. Grain boundaries as reservoirs of incompatible elements in the Earth’s mantle[J].Nature, 2004, 427: 699?703.

[17] de Ronde A A, Stünitz H. Deformation-enhanced reaction in experimentally deformed plagioclase-olivine aggregates[J].Contrib Mineral Petrol, 2007, 153(6): 699?717.

[18] Yoshino T, Takei Y, Wark D A, et al. Grain boundary wetness of texturally equilibrated rocks, with implications for seismic properties of the upper mantle[J]. J Geophys Res, 2005,110(B08205), doi: 10.1029/2004JB003544.

[19] Takei Y. Deformation-induced grain boundary wetting and its effects on the acoustic and rheological properties of partially molten rock analogue[J]. J Geophys Res, 2005, 110(B12203),doi: 10.1029/2005JB003801.

[20] 桑祖南, 夏斌, 周永勝, 等. 含礦輝長巖液態(tài)不混溶作用實(shí)驗(yàn)研究[J]. 中國科學(xué): D輯, 2003, 33(4): 353?361.SANG Zu-nan, XIA Bin, ZHOU Yong-sheng, et al. Experimental study of ore gabbro liquid immiscibility[J]. Science in China:Series D, 2003, 33(4): 353?361.

[21] Anderson O L, Isaak D G. Elastic constants of mantle minerals at high temperature[C]//Ahrens T J. Mineral physics &crystallography: A handbook of physical constants. Washington:American Geophysical Union, 1995: 64?97.

[22] Bass J D. Elastic of minerals, glasses, and melts[C]//Ahrens T J.Mineral physics & crystallography: A handbook of physical constants. Washington: American Geophysical Union, 1995:45?63.

[23] Demouchy S, Schneider S E, Mackwell S J, et al. Experimental deformation of olivine single crystals at lithospheric temperatures[J]. Geophys Res Lett, 2009, 36(L04304), doi:10.1029/2008GL036611.

[24] Ciz R, Siggins A F, Gurevich B, et al. Influence of microheterogeneity on effective stress law for elastic properties of rocks[J]. Geophysics, 2008, 73(1): E7?E14.

[25] JIN Zhen-ming, Green H W, ZHOU Yi. Melt topology in partially molten mantle peridotite during ductile deformation[J].Nature, 1994, 372: 164?167.

[26] WANG Qian, JI Shao-cheng, SUN Sheng-si, et al. Correlations between compressional and shear wave velocities and corresponding Poisson’s ratios for some common rocks and sulfide ores[J]. Tectonophysics, 2009, 469(1/4): 61?72.

[27] Sayers C M. Effect of borehole stress concentration on elastic wave velocities in sandstones[J]. Int J Rock Mech Min Sci, 2007,44(7): 1045?1052.

[28] Scott T E. The effects of stress paths on acoustic velocities and 4D seismic imaging[J]. The Leading Edge, 2007, 26(5):602?608.

[29] JI Shao-cheng, WANG Qian, Marcotte D, et al. P wave velocities,anisotropy and hysteresis in ultrahigh-pressure metamorphic rocks as a function of confining pressure[J]. J Geophys Res,2007, 112(B9), doi: 10.1029/2006JB004867.

[30] 謝鴻森. 地球深部物質(zhì)科學(xué)導(dǎo)論[M]. 北京: 科學(xué)出版社,1997:55?57.XIE Hong-sen. Introduction to materials science of the earth’s interior[M]. Beijing: Science Press, 1997: 55?57.