王丹輝

摘 要：采用有限元軟件對(duì)研究區(qū)地應(yīng)力場(chǎng)進(jìn)行了數(shù)值模擬研究。首先在ArcGIS平臺(tái)下，將模型的幾何特征點(diǎn)導(dǎo)入ANSYS，采用APDL命令流方式建立地質(zhì)模型;其次，確定了邊界約束條件和邊界外力條件;最后研究了不同深度的構(gòu)造應(yīng)力場(chǎng)變化規(guī)律。

關(guān)鍵詞：地質(zhì)模型;構(gòu)造應(yīng)力場(chǎng);數(shù)值模擬

中圖分類(lèi)號(hào)：P618 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1003-5168（2021）22-0139-04

Abstract： The in-situ stress field in the study area is numerically simulated by finite element software. Firstly， under the ArcGIS platform， the geometric feature points of the model are imported into ANSYS， and the geological model is established by APDL command flow; Then the boundary constraints and external force conditions are determined; Finally， the variation law of tectonic stress field at different depths is studied.

Keywords： geologic model;tectonic stress field;numerical simulation

通過(guò)對(duì)地應(yīng)力場(chǎng)的研究，可以得到研究區(qū)的地應(yīng)力大小及方向，再結(jié)合區(qū)域背景資料，分析地殼應(yīng)力集中區(qū)的大體位置和各個(gè)塊體及邊界斷層的差異活動(dòng)性質(zhì)，從而為區(qū)域地震研究和區(qū)域地殼穩(wěn)定性評(píng)價(jià)提供強(qiáng)有力的支撐[1]。

采用ANSYS 14.0有限元數(shù)值模擬軟件對(duì)研究區(qū)構(gòu)造應(yīng)力場(chǎng)進(jìn)行模擬分析，以獲得研究區(qū)現(xiàn)今構(gòu)造應(yīng)力場(chǎng)的分布規(guī)律。

1 地質(zhì)模型

山東省濰坊地區(qū)地質(zhì)條件復(fù)雜多樣，斷裂密集，一定程度上增加了地質(zhì)模型建立的困難程度，因此在有限元數(shù)值模擬中需要對(duì)地質(zhì)模型進(jìn)行一定程度的簡(jiǎn)化。區(qū)域計(jì)算模型以濰坊市行政界線(xiàn)為邊界，包括鄌郚-葛溝斷裂（F1）、沂水-湯頭斷裂（F2）、安丘-莒縣斷裂（F3）、昌邑-大店斷裂（F4）、景芝-賈悅斷裂（F5）、井溝-諸城斷裂（F6）、蘭底-東化山斷裂（F7）、齊河-廣饒斷裂（F9）、青州-朱劉斷裂（F10）、柴溝-膠州斷裂（F12）、百尺河-二十五里夼斷裂（F13）、紀(jì)家屯-永樂(lè)官莊斷裂（F15）、五井?dāng)嗔眩‵16）、大關(guān)-徐家河斷裂（F19）、臨朐-跋山斷裂（F22）、青州-大車(chē)溝斷裂（F27）、五圖-劉家堯斷裂（F28）、北孟-高密斷裂（F31）、安丘-辛興斷裂（F33）、吳家樓斷裂（F34）以及青島-五蓮斷裂（F38）等21條影響較大的區(qū)域斷裂[2]，模型中斷裂的走向按實(shí)際情況考慮。區(qū)域斷裂影響范圍依據(jù)斷裂活動(dòng)程度計(jì)算，本區(qū)域F1～F4斷裂活動(dòng)性較強(qiáng)，影響范圍取500 m，其他斷裂影響范圍取200 m，如圖1所示。

研究區(qū)幾何模型數(shù)據(jù)來(lái)源為1∶25萬(wàn)區(qū)調(diào)圖件和1∶5萬(wàn)地形圖，采用高斯-克呂格投影坐標(biāo)系（中央子午線(xiàn)經(jīng)度為117°）。模型前處理在ArcGIS平臺(tái)下，將模型的幾何特征點(diǎn)導(dǎo)入ANSYS，采用APDL命令流方式建立區(qū)域模型。在此過(guò)程中進(jìn)行了如下簡(jiǎn)化。

1.1 地質(zhì)體簡(jiǎn)化

研究區(qū)地層巖性較為復(fù)雜，地層從太古代至新生代地層均有出露，三大巖類(lèi)均有出現(xiàn)，因此在建模過(guò)程中對(duì)實(shí)際地質(zhì)體進(jìn)行了一定程度的簡(jiǎn)化。由于本次模擬計(jì)算主要考慮沿深度（100～1 000 m）地應(yīng)力場(chǎng)變化情況，建模時(shí)考慮了3類(lèi)巖性：侵入巖類(lèi)、沉積巖類(lèi)和斷裂破碎帶[3]。

1.2 斷裂的簡(jiǎn)化

研究區(qū)斷裂密集、構(gòu)造活動(dòng)強(qiáng)烈，在建立模型過(guò)程中，根據(jù)地質(zhì)資料分析，結(jié)合野外調(diào)查工作，篩選了21條主干斷裂，同時(shí)去掉了一些分支斷裂，在保證計(jì)算模型合理性的同時(shí)，也可簡(jiǎn)化地質(zhì)模型，適當(dāng)減少計(jì)算量。另外，根據(jù)對(duì)活動(dòng)斷裂的研究，依據(jù)斷裂規(guī)模和活動(dòng)性的不同，對(duì)斷裂的影響范圍進(jìn)行了適當(dāng)?shù)膮^(qū)分考慮，活動(dòng)性強(qiáng)、規(guī)模較大的斷裂（沂沭斷裂帶）賦予了較寬的影響范圍。模型采用平面應(yīng)變模型，Plane183單元，離散后共剖分節(jié)點(diǎn)51 654個(gè)，單元17 265個(gè)，如圖2所示。

2 物性參數(shù)

有限元數(shù)值模擬中，巖石力學(xué)參數(shù)的選取有著非常重要的意義。參數(shù)選取是否準(zhǔn)確，從一定程度上反映了模型的正確性，決定了模擬結(jié)果的客觀現(xiàn)實(shí)性和準(zhǔn)確性。在收集研究區(qū)大量巖石力學(xué)參數(shù)資料的基礎(chǔ)上，經(jīng)過(guò)反復(fù)對(duì)比、分析、試驗(yàn)，選取數(shù)據(jù)作為模擬計(jì)算參數(shù)，如表1所示。

3 邊界條件

根據(jù)區(qū)域構(gòu)造分析結(jié)論和區(qū)域震源機(jī)制解的統(tǒng)計(jì)結(jié)果及區(qū)域的全球定位系統(tǒng)（Global Positioning System，GPS）觀測(cè)數(shù)據(jù)，區(qū)域模型在EW方向上加載最大主應(yīng)力方向，在SN方向上加載最小主應(yīng)力方向。

3.1 邊界約束條件

邊界約束的原則是不使變形受到畸變，不因約束的不合理而發(fā)生應(yīng)力集中。根據(jù)區(qū)域現(xiàn)代地殼運(yùn)動(dòng)特點(diǎn)，確定邊界條件為：模型四周采用應(yīng)力邊界條件，東西方向施加最大壓應(yīng)力，南北方向施加最小壓應(yīng)力。

3.2 邊界外力條件

區(qū)域地質(zhì)體處于構(gòu)造應(yīng)力、自重應(yīng)力和溫度熱應(yīng)力等復(fù)雜應(yīng)力體系中，在進(jìn)行有限元模擬計(jì)算時(shí)僅考慮構(gòu)造應(yīng)力。構(gòu)造應(yīng)力主要依據(jù)研究區(qū)現(xiàn)今地應(yīng)力狀態(tài)施加，最大擠壓應(yīng)力按東西向施加，施加的梯度參考山東平度水壓致裂地應(yīng)力測(cè)量結(jié)果，如圖3所示。EW向和SN向分別為：SH=0.018 4H+3.406 7，R2=0.557 6;Sh=0.014 1H+2.327 6，R2=0.692 6。

采用以上應(yīng)力邊界條件分別計(jì)算深度100 m、200 m、300 m、500 m、800 m以及1 000 m的水平主應(yīng)力SH和Sh，分析不同深度的地應(yīng)力場(chǎng)變化規(guī)律，如表2所示。

4 數(shù)值模擬計(jì)算結(jié)果與分析

模型建立完成后，在Dell的Precision系列工作站上運(yùn)行，得到100 m、200 m、300 m、500 m、800 m以及1 000 m不同深度的地應(yīng)力場(chǎng)。

4.1 最大水平主壓應(yīng)力分布特征

圖4和圖5是數(shù)值模擬計(jì)算得到研究區(qū)在100 m和1 000 m深度位置最大水平主應(yīng)力值分布云圖。從圖4和圖5可直觀地看出，在100 m和1 000 m深度位置，研究區(qū)均處于壓應(yīng)力狀態(tài)，最大水平主應(yīng)力值分別為2.21～9.09 MPa和9.47～37.59 MPa。

最大值均分布在諸城市、臨朐市以南的山區(qū)內(nèi)，最大值分別為9.09 MPa和37.59 MPa;各斷裂帶內(nèi)及附近最大水平主應(yīng)力值相對(duì)較低，分別為2.21～4.17 MPa和9.47～15.10 MPa。

4.2 最小水平主壓應(yīng)力分布特征

圖6和圖7是數(shù)值模擬計(jì)算得到研究區(qū)在100 m和1 000 m深度位置最小水平主應(yīng)力值分布云圖。由圖6和圖7可直觀地看出，在100 m和1 000 m深度位置，研究區(qū)最小水平主應(yīng)力均處于壓應(yīng)力狀態(tài)，最小水平主應(yīng)力值分別為0.90～2.63 MPa和3.82～11.07 MPa;最大值均分布在諸城市、臨朐市以南的山區(qū)內(nèi)，最大值分別為2.63 MPa和11.07 MPa;各斷裂帶內(nèi)及附近最小水平主應(yīng)力值相對(duì)較低，分別為0.90～1.50 MPa和3.82～6.25 MPa。

4.3 最大剪應(yīng)力分布特征

依據(jù)摩爾-庫(kù)侖（Mohr-Coulomb）強(qiáng)度理論，地殼的破裂現(xiàn)象與最大剪應(yīng)力有關(guān)。根據(jù)數(shù)值模擬計(jì)算結(jié)果，可以得到研究區(qū)不同深度（100 m、200 m、300 m、500 m、800 m、1 000 m）范圍的最大剪應(yīng)力分布，如圖8和圖9所示。由圖8和圖9可以看出，在100 m和1 000 m深度位置，研究區(qū)最大剪應(yīng)力值分別為0.61～3.25 MPa和2.60～13.35 MPa;最大值均分布在諸城市以南的山區(qū)內(nèi)，最大值分別為4.45 MPa和13.35 MPa;各斷裂帶內(nèi)及附近最大剪應(yīng)力值相對(duì)較低，分別為0.61～0.98 MPa和2.60～4.80 MPa。

5 結(jié)語(yǔ)

從研究區(qū)各深度范圍最大水平主應(yīng)力、最小水平主應(yīng)力、最大剪應(yīng)力大小隨深度變化特征來(lái)看，最大水平主應(yīng)力、最小水平主應(yīng)力、最大剪應(yīng)力的分布形態(tài)具有相似的特征，均受區(qū)域構(gòu)造格架的控制。在南部臨朐、諸城山區(qū)應(yīng)力值偏高，北部壽光、寒亭、昌邑一帶的平原區(qū)和東部高密地區(qū)應(yīng)力值偏低，中部昌樂(lè)、寒亭南以及昌邑中部地區(qū)應(yīng)力值較高，各斷裂帶影響范圍內(nèi)應(yīng)力值明顯偏低。值得注意的是，在各斷裂帶交匯區(qū)尤其是沿沂沭斷裂帶四條主干斷裂易形成應(yīng)力集中區(qū)域，在沂沭斷裂帶之內(nèi)局部地區(qū)（尤其是在鄌郚-葛溝斷裂和沂水-湯頭斷裂中段以及寒亭區(qū)南部地區(qū)）有應(yīng)力集中現(xiàn)象。

參考文獻(xiàn)：

[1]白世偉，李光煜.二灘水電站壩區(qū)巖體應(yīng)力場(chǎng)研究[J].巖石力學(xué)與工程學(xué)報(bào)，1982（1）：45-56.

[2]陳慶宣，王維襄，孫葉，等.巖石力學(xué)與構(gòu)造應(yīng)力場(chǎng)分析[M].北京：地質(zhì)出版社，1998：24.

[3]鄧起東，張培震，冉康勇，等.中國(guó)活動(dòng)構(gòu)造基本特征[J].中國(guó)科學(xué)D輯，2002（12）：1020-1031.