陳玉香

(中國地震局地質(zhì)研究所,北京 100029)

靜水壓下礦物顆粒邊界力學(xué)性質(zhì)的數(shù)值模擬研究

陳玉香

(中國地震局地質(zhì)研究所,北京 100029)

隨著地球科學(xué)研究的深入以及計算機(jī)技術(shù)的飛速發(fā)展,有限元數(shù)值模擬正成為地震機(jī)理研究中的一種強(qiáng)有力的技術(shù)方法。從震源物理力學(xué)過程、斷層類型及力學(xué)機(jī)制、區(qū)域構(gòu)造應(yīng)力場、波速比及地下水位變化、溫度效應(yīng)與應(yīng)力觸發(fā)等方面模擬研究地震孕育機(jī)理取得了很大的進(jìn)展。但是,目前的地震機(jī)理有限元模擬研究中尚存在一些關(guān)鍵的問題,如模型中塊體之間的界面效應(yīng)往往被忽視,參數(shù)(彈性模量、粘滯系數(shù))的選取不確定。

地震作為一種自然現(xiàn)象,主要是地球內(nèi)部動力作用的結(jié)果,探明巖石、礦物的力學(xué)性質(zhì)與變形機(jī)制,能為了解地球內(nèi)部動力學(xué)過程、地震物理機(jī)制與震源分布模型提供科學(xué)依據(jù)。而巖石、礦物的力學(xué)性質(zhì)與變形機(jī)制取決于礦物成分、結(jié)構(gòu)、礦物顆粒邊界效應(yīng)和溫壓條件等多種因素。邊界效應(yīng)在材料科學(xué)中,早已被證實對材料的力學(xué)性能有重要影響。邊界效應(yīng)同樣廣泛存在于不同尺度的地質(zhì)體中(地殼與地幔之間,板塊與板塊之間、地體與地體之間、巖石與巖石之間、礦物與礦物之間)。地球科學(xué)工作者也開始關(guān)注邊界效應(yīng),認(rèn)為絕大多數(shù)7級以上的地震發(fā)生在活動地塊邊界的活動斷裂帶上。汶川大地震是青藏高原東擴(kuò)、地塊邊界應(yīng)力積累和釋放的結(jié)果。在一些巖石流變實驗研究中,發(fā)現(xiàn)礦物集合體內(nèi)部由于應(yīng)力分布不均引起位錯密度不均勻。但是,邊界效應(yīng)的研究在地球動力學(xué)中依然是個薄弱環(huán)節(jié)。礦物顆粒邊界效應(yīng)作為影響巖石力學(xué)與變形機(jī)制的重要因素之一,目前還無法通過實驗?zāi)M進(jìn)行研究,而通過數(shù)值模擬研究還只側(cè)重于礦物集合體中顆粒邊界移動的結(jié)果、邊界濕化對礦物顆粒集合體剪應(yīng)力分布的影響,對于在一定溫壓條件下,礦物集合體(或巖石)中礦物顆粒邊界的應(yīng)力分布還沒有系統(tǒng)的研究。

本論文用ANSYS有限元軟件三維模擬室溫、1.6GPa靜水壓下礦物集合體的應(yīng)力場,進(jìn)而探討礦物顆粒邊界層的力學(xué)性質(zhì)。從小尺度上(介觀尺度)探討礦物顆粒邊界效應(yīng),為巖石(礦物集合體)的變形機(jī)制研究與地震機(jī)理研究以及地球動力學(xué)研究積累科學(xué)資料。

根據(jù)有限元的基本思想(是將連續(xù)的求解區(qū)域離散為一組有限個、且按一定方式相互連接在一起的單元的組合體),將邊長為1μm的立方體形態(tài)的兩種礦物組成的“三明治”式集合體(集合體結(jié)構(gòu)從下至上為礦物1-顆粒邊界層-礦物2-顆粒邊界層-礦物1),根據(jù)對稱性取1/8作為計算模型,進(jìn)行有限元三維模擬計算。在計算模型中,以對稱中心為坐標(biāo)原點(diǎn)建立坐標(biāo)系,采用規(guī)則網(wǎng)格,三個塊體X、Z方向劃分25等份,Y方向從下至上依次劃分25等份、1等份、50等份,即礦物顆粒網(wǎng)格單元尺寸為20nm×20nm×20nm,邊界層的網(wǎng)格單元尺寸為20nm×邊界層厚度×20nm;網(wǎng)格單元使用ANSYS軟件內(nèi)置的SOLID185單元網(wǎng)格后,整個計算模型,有53404個節(jié)點(diǎn),47500個單元。1μm(1000nm)正是介觀尺度(10nm～1000nm)與宏觀尺度(＞1μm)的臨界點(diǎn),在介觀尺度上,分子團(tuán)簇間的非鍵相互作用起主導(dǎo)作用,晶體顆粒在1μm時其力學(xué)性質(zhì)還是與宏觀力學(xué)性質(zhì)一樣。那么,網(wǎng)格單元中的每個節(jié)點(diǎn)可以看成是一個分子團(tuán)簇。在受壓條件下,集合體中各個單元以特定方式的相互連接作用將力傳到處于不同空間位置的分子團(tuán)簇(節(jié)點(diǎn)),從而可以模擬計算出任意分子團(tuán)簇(節(jié)點(diǎn))在整個體系受外力作用下的應(yīng)力狀態(tài)。

本文的模型中,礦物集合體中的顆粒邊界層不僅具有幾何邊界而且具有物理邊界。邊界層由兩種礦物晶體均勻混合組成,其彈性模量依據(jù)經(jīng)典的混合模型進(jìn)行計算。礦物顆粒與邊界層之間為靜摩擦接觸關(guān)系,ANSYS軟件中的面-面接觸可以轉(zhuǎn)化為兩個面上的節(jié)點(diǎn)與節(jié)點(diǎn)接觸(即面-面接觸對象集為節(jié)點(diǎn)-節(jié)點(diǎn)),也即宏觀上是面-面之間的粘結(jié)接觸,而介觀尺度上是分子團(tuán)簇-分子團(tuán)簇之間的摩擦接觸。接觸單元用ANSYS內(nèi)置的CONTA174,目標(biāo)單元為ANSYS內(nèi)置的TRAGE170。計算模型中,礦物2的頂面與邊界層的底面接觸,礦物2的頂面為目標(biāo)面,邊界層底面為接觸面;礦物1的底面與邊界層頂面接觸,礦物1的底面為目標(biāo)面,邊界層頂面為接觸面。接觸算法是用增廣拉格朗日法,接觸探測采用高斯插值點(diǎn),接觸剛度因子為0.1?？紤]到邊界層是兩種礦物均勻混合的,允許有滲透,故將滲透容忍因子設(shè)為0.1。

首先,模擬計算室溫、1.6GPa壓力下彈性模量(楊氏模量)相差38.3560GPa、泊松比相差0.0181的鎂橄欖石-透輝石集合體,在室溫、1.6GPa靜水壓下完全彈性變形時的應(yīng)力場。保持邊界層的彈性模量(鎂橄欖石、透輝石各占一半體積分?jǐn)?shù),依據(jù)Hill模型計算)不變,摩擦系數(shù)為0.65,改變邊界層的厚度,邊界層厚度從1nm～10nm之間等間距取19個值,模擬計算鎂橄欖石-透輝石集合體19個不同的邊界層厚度模型的應(yīng)力場,并選取在所有模型中處于同一幾何位置的界面層上的5個點(diǎn)(5個處在不同空間位置的分子團(tuán)簇),分析集合體中5個不同空間位置的分子團(tuán)簇的等效應(yīng)力以及系統(tǒng)中最大等效應(yīng)力與邊界層厚度的關(guān)系。結(jié)果表明,在靜水壓下,鎂橄欖石-透輝石集合體顆粒邊界層應(yīng)力集中,邊界層應(yīng)力分布不均勻,越遠(yuǎn)離對稱中心,其等效應(yīng)力越大;邊界層的等效應(yīng)力、集合體最大等效應(yīng)力與邊界層厚度呈線性負(fù)相關(guān),隨著邊界層厚度的增加,集合體最大等效應(yīng)力減小,介觀尺度上晶體界面不同空間位置的分子團(tuán)簇的等效應(yīng)力也不同程度地減少。此結(jié)果可以用來解釋前人研究中發(fā)現(xiàn)“顆粒邊界比顆粒內(nèi)部的滑動變形要大,但隨著邊界層厚度的增加,這種差距減少”的現(xiàn)象。

其次,保持邊界層物質(zhì)彈性參數(shù)不變(鎂橄欖石、透輝石各占一半體積分?jǐn)?shù),依據(jù)Hill模型計算)、邊界層厚度取1nm,改變摩擦系數(shù),摩擦系數(shù)從0.5～0.85之間等間距取15個值,并選取10個點(diǎn)(不同空間位置的10個分子團(tuán)簇),分析鎂橄欖石-透輝石集合體中10個不同空間位置的分子團(tuán)簇的等效應(yīng)力、集合體中最大等效應(yīng)力以及摩擦應(yīng)力與摩擦系數(shù)之間的關(guān)系。結(jié)果表明,分子團(tuán)簇的等效應(yīng)力以及接觸摩擦應(yīng)力在鎂橄欖石接觸界面上比在透輝石接觸界面上的大(即X、Z坐標(biāo)相同,Y坐標(biāo)不同),而且相對應(yīng)位置的分子團(tuán)簇的等效應(yīng)力以及接觸摩擦應(yīng)力隨摩擦系數(shù)變化的程度基本一致。顆粒邊界層的應(yīng)力分布不均勻,分別在與鎂橄欖石或透輝石接觸的界面上的不同空間位置的分子團(tuán)簇離對稱中心越遠(yuǎn),其等效應(yīng)力以及接觸摩擦應(yīng)力越大。顆粒邊界層的等效應(yīng)力以及接觸摩擦應(yīng)力與摩擦系數(shù)呈線性正相關(guān)變化,但變化的程度不一致。在與鎂橄欖石或透輝石相接觸的界面上,界面內(nèi)部區(qū)域的不同位置的分子團(tuán)簇距離對稱中心越近,其等效應(yīng)力以及接觸摩擦應(yīng)力隨摩擦系數(shù)變化的程度越小,而在靠近邊界層受力界面區(qū)域內(nèi)不同位置的分子團(tuán)簇的變化程度基本一致。

第三,保持邊界層厚度為1nm、摩擦系數(shù)為0.65,依據(jù)Reuss模型改變邊界層物質(zhì)的體積模量與剪切模量,進(jìn)而計算得到模擬所需要的楊氏模量與泊松比,對鎂橄欖石-透輝石集合體分別模擬計算應(yīng)力場21次,分析10個不同空間位置的分子團(tuán)簇的應(yīng)力分布,并且以同樣邊界條件和有限元模型模擬計算室溫、1.6GPa壓力下彈性模量(楊氏模量)相差64.9220 GPa、泊松比相差0.0393的透輝石-鈣長石集合體和室溫、1.6GPa壓力下楊氏模量相差3.2577GPa、泊松比相差0.1618的鈣長石-石英集合體的應(yīng)力場各21次,對比分析邊界層物質(zhì)的彈性參數(shù)對集合體應(yīng)力分布的影響。研究結(jié)果表明,礦物顆粒邊界層的等效應(yīng)力與邊界層物質(zhì)的彈性參數(shù)呈高斯函數(shù)關(guān)系,利用此函數(shù)關(guān)系式結(jié)合彈性參數(shù)與波速的計算公式,可以幫助研究應(yīng)力與波速的關(guān)系;當(dāng)邊界層物質(zhì)的彈性參數(shù)為兩礦物的彈性參數(shù)的平均值時,集合體邊界層應(yīng)力集中程度最小;不同彈性參數(shù)的礦物組成的集合體,礦物之間的彈性參數(shù)尤其是泊松比相差越大,界面層應(yīng)力集中程度越大,這說明泊松比相差越大的集合體越容易發(fā)生塑性變形或破壞。此結(jié)果可以用來解釋一些巖石力學(xué)現(xiàn)象,也可用于復(fù)合材料的設(shè)計。

通過三維有限元模擬計算、分析,本文的主要結(jié)論如下:

(1)立方體形態(tài)的兩種不同礦物集合體,在室溫、1.6GPa靜水壓下完全彈性變形時,礦物集合體內(nèi)部應(yīng)力分布不均,應(yīng)力集中在顆粒邊界,邊界層的應(yīng)力分布也不均勻。說明礦物顆粒邊界層是弱化帶,在外力作用下邊界層容易發(fā)生塑性變形或是破裂。

(2)邊長為1μm的立方體形態(tài)的兩種不同礦物集合體,在室溫、1.6GPa靜水壓下完全彈性變形時,介觀尺度上,礦物晶體內(nèi)部的分子團(tuán)簇與晶體界面的分子團(tuán)簇的應(yīng)力狀態(tài)不同,晶體界域的分子團(tuán)簇的等效應(yīng)力大,且在界面域內(nèi),離對稱中心越遠(yuǎn),其分子團(tuán)簇的等效應(yīng)力越大。

(3)立方體形態(tài)的兩種不同礦物集合體,在室溫、1.6GPa靜水壓下完全彈性變形時,礦物顆粒邊界層的等效應(yīng)力、集合體最大等效應(yīng)力與邊界層厚度呈線性負(fù)相關(guān),隨著邊界層厚度的增加,集合體最大等效應(yīng)力減小,介觀尺度上晶體界面層的分子團(tuán)簇的等效應(yīng)力也隨著離對稱中心的距離的不同而不同程度地減少。

(4)在給定的模擬條件下,顆粒邊界層的等效應(yīng)力以及接觸摩擦應(yīng)力與摩擦系數(shù)呈線性正相關(guān)變化,但變化的程度不一致,界面內(nèi)部區(qū)域的不同位置的分子團(tuán)簇距離對稱中心越近,其等效應(yīng)力以及接觸摩擦應(yīng)力隨摩擦系數(shù)變化的程度越小,而在靠近邊界層受力界面區(qū)域內(nèi)不同位置的分子團(tuán)簇的變化程度基本一致。

(5)在室溫、1.6GPa靜水壓下完全彈性變形時,不同彈性參數(shù)的礦物組成的集合體或巖石,礦物之間的彈性參數(shù)尤其是泊松比相差越大,界面層應(yīng)力集中程度越大。

(6)在室溫、1.6GPa靜水壓下完全彈性變形時,礦物顆粒邊界層的等效應(yīng)力與邊界層物質(zhì)的彈性參數(shù)呈高斯函數(shù)關(guān)系,當(dāng)邊界層物質(zhì)的彈性參數(shù)為兩礦物的彈性參數(shù)的平均值時,集合體邊界層應(yīng)力集中程度最小。

邊界效應(yīng);等效應(yīng)力;應(yīng)力集中;礦物顆粒邊界;數(shù)值模擬

(作者電子信箱,陳玉香:chyuxiang@tom.com)

P574.1;

A;

10.3969/j.issn.0235-4975.2011.07.013