研究區現今構造應力場數值模擬分析

王丹輝

摘 要：采用有限元軟件對研究區地應力場進行了數值模擬研究。首先在ArcGIS平臺下，將模型的幾何特征點導入ANSYS，采用APDL命令流方式建立地質模型;其次，確定了邊界約束條件和邊界外力條件;最后研究了不同深度的構造應力場變化規律。

關鍵詞：地質模型;構造應力場;數值模擬

中圖分類號：P618 文獻標識碼：A 文章編號：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

通過對地應力場的研究，可以得到研究區的地應力大小及方向，再結合區域背景資料，分析地殼應力集中區的大體位置和各個塊體及邊界斷層的差異活動性質，從而為區域地震研究和區域地殼穩定性評價提供強有力的支撐[1]。

采用ANSYS 14.0有限元數值模擬軟件對研究區構造應力場進行模擬分析，以獲得研究區現今構造應力場的分布規律。

1 地質模型

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

研究區幾何模型數據來源為1∶25萬區調圖件和1∶5萬地形圖，采用高斯-克呂格投影坐標系（中央子午線經度為117°）。模型前處理在ArcGIS平臺下，將模型的幾何特征點導入ANSYS，采用APDL命令流方式建立區域模型。在此過程中進行了如下簡化。

1.1 地質體簡化

研究區地層巖性較為復雜，地層從太古代至新生代地層均有出露，三大巖類均有出現，因此在建模過程中對實際地質體進行了一定程度的簡化。由于本次模擬計算主要考慮沿深度（100～1 000 m）地應力場變化情況，建模時考慮了3類巖性：侵入巖類、沉積巖類和斷裂破碎帶[3]。

1.2 斷裂的簡化

研究區斷裂密集、構造活動強烈，在建立模型過程中，根據地質資料分析，結合野外調查工作，篩選了21條主干斷裂，同時去掉了一些分支斷裂，在保證計算模型合理性的同時，也可簡化地質模型，適當減少計算量。另外，根據對活動斷裂的研究，依據斷裂規模和活動性的不同，對斷裂的影響范圍進行了適當的區分考慮，活動性強、規模較大的斷裂（沂沭斷裂帶）賦予了較寬的影響范圍。模型采用平面應變模型，Plane183單元，離散后共剖分節點51 654個，單元17 265個，如圖2所示。

2 物性參數

有限元數值模擬中，巖石力學參數的選取有著非常重要的意義。參數選取是否準確，從一定程度上反映了模型的正確性，決定了模擬結果的客觀現實性和準確性。在收集研究區大量巖石力學參數資料的基礎上，經過反復對比、分析、試驗，選取數據作為模擬計算參數，如表1所示。

3 邊界條件

根據區域構造分析結論和區域震源機制解的統計結果及區域的全球定位系統（Global Positioning System，GPS）觀測數據，區域模型在EW方向上加載最大主應力方向，在SN方向上加載最小主應力方向。

3.1 邊界約束條件

邊界約束的原則是不使變形受到畸變，不因約束的不合理而發生應力集中。根據區域現代地殼運動特點，確定邊界條件為：模型四周采用應力邊界條件，東西方向施加最大壓應力，南北方向施加最小壓應力。

3.2 邊界外力條件

區域地質體處于構造應力、自重應力和溫度熱應力等復雜應力體系中，在進行有限元模擬計算時僅考慮構造應力。構造應力主要依據研究區現今地應力狀態施加，最大擠壓應力按東西向施加，施加的梯度參考山東平度水壓致裂地應力測量結果，如圖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。

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

4 數值模擬計算結果與分析

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

4.1 最大水平主壓應力分布特征

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

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

4.2 最小水平主壓應力分布特征

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

4.3 最大剪應力分布特征

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

5 結語

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

參考文獻：

[1]白世偉，李光煜.二灘水電站壩區巖體應力場研究[J].巖石力學與工程學報，1982（1）：45-56.

[2]陳慶宣，王維襄，孫葉，等.巖石力學與構造應力場分析[M].北京：地質出版社，1998：24.

[3]鄧起東，張培震，冉康勇，等.中國活動構造基本特征[J].中國科學D輯，2002（12）：1020-1031.