雁列斷層應力場影響因素數值模擬分析

王 時, 張環曦, 周安聘, 董 博, 章 陽, 紀春玲

(石家莊地震監測中心站, 河北 石家莊 050021)

0 引言

雁列構造是野外常見的一種間斷、不連續的經典構造形式。眾多學者對雁列構造進行了大量的研究工作,對斷層的歷史活動、相互制約和影響、地震機制以及前兆滑動過程等進行了分析和討論。唐春安等[1]為了更好研究巖石的破裂過程,在考慮材料力學非均勻特點的基礎上,對地震孕育過程進行數值模擬,并給出了應力場、變形場和微震活動的時空分布;郭玲莉等[2]通過對三種結構模型的黏滑實驗,探討了應力降與震級之間的關系,發現震源尺度是影響震級的重要因素;Dieterich[3]認為摩擦系數與時間具有一定的關聯性,隨著摩擦的增加,接觸面積也在逐步加大,在垂直壓力和接觸強度不變的情況下,摩擦系數與接觸面為正比關系,并用巖石摩擦實驗證實了這個規律;王想等[4]認為巖石地震活動幕時段背景下需排除邊界動力變化引起的區域應力場調整,再結合小震平靜區或空間圖像,判斷平靜屬于巖石活動過程中出現的間歇平靜還是震前平靜異常;朱瑤宏等[5]研究了不同溫度、動應力、初始偏應力、圍壓作用對黏土塑性應變的影響,發現影響因素的預測模型與試驗吻合性較好;羅崇亮等[6]通過研究紅層泥巖含水率對其力學特性及鄧肯-張模型參數的影響,發現含水率對泥巖黏聚力和內摩擦角產生影響,進而導致結果出現差異;許健生等[7]對巖石進行了黏滑實驗,發現巖石的黏滑錯動分為三個部分,預滑、黏滑和止滑,黏滑過程是由多次黏滑錯動過程組成;胡志平等[8]基于動三軸試驗分析了不同圍壓條件下和加載方式對壓實黃土應力應變關系的影響規律;王想等[9]利用傳統和數字地震學方法,對震群與強震孕震關系進行了分析;任雅瓊等[10]在實驗室開展巖石變形實驗,觀測相關物理量的演化特征和規律,對壓性雁列斷層進行變形實驗,根據應力時間曲線將具有前兆意義的變形階段分為強偏離線性階段、亞失穩階段和失穩階段,利用紅外熱像儀記錄各個階段的熱場變化。但是影響雁列斷層構造應力場的因素以及應力分布狀態前人并未給出足夠多的數值模擬實驗,本文結合數值模擬做了一定的分析。

1 摩擦接觸有限元數值計算方法

1.1 接觸面力條件

由牛頓第三定律可知,接觸面力應滿足:

(1)

1.2 摩擦力的計算

接觸面上的力:主要是摩擦力的計算,從節點N的法向接觸力為fs,則它的最大摩擦力值為Fy=u|fs|,其中μ為摩擦系數。

設在上一時刻tn從節點ns的摩擦力為Fn,則當前時刻tn+1可能產生的摩擦力F*

F*=Fn-kΔe

(2)

式中:k是界面剛度;Δe是從節點的位移,

(3)

當前時刻tn+1的摩擦力Fn+1根據下式計算,

(4)

最后按作用與反作用原理,計算對應節點的摩擦力[2]。

2 模型的建立

本文建立了二條未貫通的二維理論雁列斷層模型,分析了雁列走滑斷層黏滑運動特征。圖1為雁列斷層有限元網格劃分模型,在斷層兩側做了加密處理,使計算精度進一步提高,在斷層上選取了6個節點以便對比研究。該模型長160 km,寬90 km,該模型設置為均勻、各向同性,斷層面為摩擦接觸,這樣可以更精確的計算斷層黏滑運動的過程,下邊界垂向固定,水平方向自由;左邊界施加水平方向右的擠壓,擠壓速率為 1.2 cm/a,垂直方向自由;右邊界水平方向固定,垂向自由。庫倫摩擦定律被納入塑性理論的研究范疇,基于彈塑性理論的返回映射方案已被應用到有限元摩擦接觸分析,該方案使算法的收斂行為和可靠性產生了本質的提高。

圖1 雁列斷層有限元網格劃分模型Fig.1 Finite element meshing model of the echelon fault

為了對比分析幾何形態和黏滯系數對雁列斷層黏滑運動的影響,建立了5個模型,分為了兩組。第一組為基本模型1、模型2、模型3,是研究幾何間距對斷層黏滑運動的影響;第二組為模型4、基本模型1、模型5,是為了研究斷層的黏滯系數對黏滑運動的影響,具體的模型列于表1(參考了西部地區的擠壓速率及相關地質資料)。

表1 模型的相關參數Table 1 Relevant parameters of the model

3 幾何間距對斷層黏滑運動的影響

3.1 幾何間距與應力降之間的關系

圖2展示了不同幾何間距(10 km、15 km、20 km)情況下雁列斷層的黏滑曲線。為了便于對比分析,選取斷層相同位置分析幾何間距對雁列斷層黏滑運動的影響。曲線都是類似于階梯式,曲線的平直段代表斷層處于閉鎖狀態(黏著狀態),斷層兩盤沒有相對運動,剪應力處于積累狀態,陡立段代表斷層的突然滑動(活動狀態),短時間內兩盤有較大的相對位錯,可以認為是一次地震的發生,形成應力降,應力減小。從曲線的變化情況可知,3個模型中斷層兩盤的相對運動都是黏著狀態與滑動狀態的交替發生。

圖2 不同幾何間距的剪應力隨時間變化曲線Fig.2 Time-dependent curves of shear stress with different geometric distances

總的來看,三個模型的都處于應力積累的狀態。但斷層幾何間距越小,斷層相互之間的控制力越強,應力場較大,容易發生黏滑運動,故而地震較多;幾何間距增大,斷層之間控制力較弱,構造應力場較弱,斷層受到的應力較小,故黏滑運動較為平穩。當斷層向前滑動時,剪應力被釋放,同時產生地震并形成相應的應力降。由此認為,沿斷裂面發生黏滑時,一次地震就是與斷裂黏滑有關的一個應力降,從而解釋了地震時應力降低的現象。地震時斷層位錯越大,產生的應力降也越大。對于不同的幾何間距,地震時的應力降也各不相同,當幾何間距為10 km(模型1)時,地震產生的應力降最大約為10 MPa,當幾何間距為15 km(模型2)時,應力降最大約為8.5 MPa,當幾何間距為20 km(模型3)時,應力降最大約為6.5 MPa。

3.2 應力與應變之間的關系

如圖3所示,當幾何形態變化從10 km到20 km逐漸增大時,斷層的應力和應變均逐漸降低,是因為斷層間距過大,斷層之間的相互作用力降低,構造應力場變小的緣故。當幾何間距為10 km時,斷層的應變和應力值最大,應變量超過3E-5有6次,幾何間距為15 km時,超過3E-5有2次,當間距為20 km時,超過3E-5有1次,由圖可以看出,幾何間距越小,應變量較大的次數較多,容易發生震級較大的地震,幾何間距越大,應變量較小,但是頻率高,容易卸載剪應力,應力不能持續積累,故而結構較為穩定,發生較大震級的地震的可能性小。

圖3 模型應力與應變隨時間的變化關系Fig.3 Variation of model stress and strain with time

3.3 幾何間距與剪應力變化云圖

圖4是斷層第3 000步、6 000步和10 000步的應力變化云圖。每一個子步即每一個階段的應力值是不同的,且斷層的應力分布是不均勻的,應力集中在貫通區域(即巖橋區),由中間向兩邊依次減小。當幾何間距為10 km、15 km、20 km時,斷層之間的最大應力分別為0.255E9、0.242E9、0.23E9。可見,隨著間距的增大,斷層之間相互的影響力在逐漸減小。

圖4 模型1～3雁列斷層剪應力隨子步變化的應力云圖Fig.4 Nephogram of shear stress variation of echelon fault with sub step in model 1～3

3.4 幾何間距與應力分量之間的關系

圖5中第一層、第二層、第三層分別為第一、第二、第三主應力的應力云圖(應力分量圖的排序與表2所示內容一致)。可以看出,第一主應力主要集中在斷層的貫通區域,與斷層面相互垂直,隨幾何間距的增大,斷層之間的應力逐漸減小,趨于平靜;第二主應力與斷層呈40°～ 45°夾角,離貫通區域越近,應力值越大,越遠離該區域,應力越小;第三主應力主要分布在斷層本身,即斷層本身主要受到第三主應力的影響,而非第一和第二主應力,斷層的受力呈現不均勻趨勢,在斷層末端受力最小,斷層的巖橋區應力最大。

表2 應力分量與幾何間距之間的數值模擬

圖5 模型1～3在幾何間距影響下應力分量的應力云圖Fig.5 Stress component nephogram of model 1～3 under the influence of geometric spacing

如圖6所示,通過對不同模型雁列斷層的橫向對比可以看出,隨著幾何間距的增大,第一和第二主應力呈現下降的趨勢。而第三主應力則呈現對稱的形式,即在15 km處斷層的第三主應力最大,斷層的結構相對不穩定,易發生滑動位移,從而產生地震。

圖6 模型1～3主應力與幾何間距之間的關系Fig.6 Relationship between principal stress and geometric distance in model 1～3

4 黏滯系數對雁列斷層黏滑運動的影響

4.1 黏滯系數與應力降之間的關系

為了分析黏滯系數對雁列斷層黏滑運動的影響,本章采用了第二組模型,即模型4(1E22)、基本模型1(1E23)、模型5(1E24),選取了相同的區域進行對比,三組雁列斷層的主要模型參數相同,具體參數見表1。

由圖7可知黏滯系數越大,斷層之間的應力值越大,當發生黏滑運動時,產生地震,形成的應力降也較大,當黏滯系數為1E22(模型4)時,地震產生的應力降最大約為4.5 MPa,當黏滯系數為1E23(模型1)時,應力降最大約為8.5 MPa,當黏滯系數為1E24(模型5)時,應力降最大約為9.5 MPa。

圖7 不同黏滯系數的剪應力隨時間變化曲線Fig.7 Time varying curves of shear stress with different viscosity coefficients

4.2 應力與應變之間的關系

如圖8所示,當黏滯系數從1E22到1E24逐漸增大時,斷層的應力和應變是兩種不同的情況,應力是逐漸上升而應變則呈現下降的趨勢。當黏滯系數較大時,斷層結構較為穩定,應變較小,而應力一致處于積累狀態,故應力值較大、應變量小。當黏滯系數較小時,結構松散,容易發生較大的應變,形成應力降,應力能被卸載故而不能累積,故應力值較小、應變量較大。

由圖可以看出當黏滯系數為1E22時,應變量超過3E-5的有4次;當黏滯系數為1E23時應變量超過最大值3E-5的有4次;當黏滯系數為1E24時,應變量超過最大值3E-5的有0次。當黏滯系數越小,應變越大,越容易發生震級較大的地震。黏滯系數越大,斷層較為堅固,不易發生形變,但應力不斷持續累積,平時容易發生較小的變形,發生震級較小的地震,當應力值積累到一定的程度,超過巖石的承受極限時,發生大的應變,就會發生震級較大的地震。

4.3 黏滯系數與剪應力變化云圖

圖9為各雁列斷層剪應力在第3 000步、6 000步和10 000步的剪應力變化云圖,隨著子步的變化,斷層的剪應力也在發生改變,每一子步即是一年,可以看出黏滯系數的改變不會影響斷層的應力分布,即斷層的空間區域不同,受到的應力也不同。

圖9 各類模型雁列斷層剪應力隨時間變化的應力云圖Fig.9 Shear stress nephogram of echelon fault with time in various models

當黏滯系數為1E22～1E24時,斷層的最大應力值分別為0.126E9、0.255E9、0.281E9。當其他參數條件不變,隨著黏滯系數的增大,斷層之間的應力不斷積累,應力值逐漸增大。

4.4 黏滯系數與應力分量之間的關系

圖10為主應力分量的應力云圖,其排列順序與表3所示內容相同。從圖可以看出,隨著黏滯系數的增大,第一、第二、第三主應力逐漸增大且第一主應力與斷層垂直;第二主應力與斷層呈40°～ 45°夾角;第三主應力與斷層平行。

表3 應力與黏滯系數之間的數值模擬

圖10 在黏滯系數影響下各應力分量的應力云圖Fig.10 Nephogram of each stress component under the influence of viscosity coefficient

如圖11所示,通過對不同模型雁列斷層的橫向對比可以看出,隨著黏滯系數的增大,斷層的第一主應力、第二主應力和第三主應力逐漸增大,符合斷層的剪應力變化和應力-應變之間規律相符合。

圖11 各黏滯系數與應力之間的關系Fig.11 Relationship between viscosity coefficient and stress

5 結論

通過對比的方法可以看出,斷層的幾何形態和黏滯系數是影響雁列斷層黏滑活動的重要因素,但這兩個因素不會影響斷層應力分布不均勻的狀態。

(1) 斷層的幾何間距在10～20 km之間時,隨著間距的增大,斷層之間的相互控制力在逐漸減小(應力和應變減小),不易發生震級較大的地震,斷層之間的活動相對平緩。

(2) 當黏滯系數在1E22～1E24范圍之間波動時,黏滯系數越小,結構越不穩定,易發生錯動從而產生地震,并伴隨應力降的產生;當黏滯系數較大時,結構較為堅固,不易發生較大變形,但應力在不斷地積累,積蓄的應力值逐漸增大,一旦積累到一定程度超過巖石的承受極限,巖石就會發生破碎形成震級較大的地震,但這需要一個較長的地震周期。

(3) 斷層的空間區域不同,受到的應力也不同,越接近巖橋區域,受到的應力越大,由巖橋區域向斷層的兩側逐漸延伸,受到的應力依次減小,體現了斷層空間區域受力不均勻。應力在斷層的巖橋區域最為集中,兩端的應力較小,摩擦較為平緩。

(4) 第一、第二、第三主應力與斷層呈一定的角度,第一主應力控制斷層的巖橋區域,第二主應力控制斷層的應力向外延伸,與斷層呈40°～ 45°夾角,第三主應力主要控制斷層本身。