有限元建模與實測應變分析
——以九江地震臺為例

李 萌 盧和雄 趙愛平 嚴 麗,4

1 東華理工大學自然資源部環(huán)鄱陽湖區(qū)域礦山環(huán)境監(jiān)測與治理重點實驗室,南昌市廣蘭大道418號,330013 2 東華理工大學江西省數字國土重點實驗室,南昌市廣蘭大道418號,330013 3 江西九江揚子塊體東部地球動力學野外科學觀測研究站,江西省九江市星光大道138號,332006 4 江西省防震減災與工程地質災害探測工程研究中心,南昌市廣蘭大道418號,330013

當區(qū)域應變能積累到一定閾值時,會以地震的形式釋放[1],因此開展應變觀測,為地震趨勢判斷提供科學依據,對地震預測具有重要意義[2]。學者們通常利用全球定位系統(tǒng)(GPS)技術監(jiān)測速度場來分析區(qū)域應變[3],但江西地區(qū)GPS地殼形變監(jiān)測站分布稀疏,存在測點覆蓋不足的缺陷[4]。2007年九江地震臺開展了數字化地球物理應變觀測,彌補了這一缺陷。區(qū)域應變受氣象及區(qū)域地質構造等條件的影響[5],二維數值模型難以適應復雜的地形,因此在得不到解析解的情況下,可采用有限元分析建立區(qū)域地質模型[6],而利用數字高程模型(DEM)數據可快速建立三維地質模型,用于有限元分析,可更直觀地反映應力及應變分布[7]。本文利用DEM數據建立有限元三維地質模型,采用有限元分析軟件ANSYS模擬九江地震臺附近區(qū)域有限元應變解,同時結合九江地震臺實測溫度和應變數據,分析有限元建模與實測應變的關系。

1 地震臺區(qū)域三維建模

1.1 臺址環(huán)境

本文以九江地震臺附近區(qū)域為例進行有限元建模。九江地震臺(29.65°N,116.01°E,海拔110 m)位于江西北部廬山西北側山腳,廬山山體為NE-SW向伸展的褶皺斷塊山體,整體呈腎狀。觀測山洞內溫差較小,洞體基巖為硅質灰?guī)r[8],洞內配置有SS-Y型伸縮儀、水管儀和垂直擺等形變觀測儀器,觀測數據質量較好。

1.2 數字化三維地形建模

采用GDEMV2分辨率為30 m的DEM數據(www.gscloud.cn)進行臺站區(qū)域數字化三維地形建模,研究區(qū)域范圍為29.41°～29.66°N、115.86°～116.08°E。首先,利用DEM數據提取等高線;然后,利用等高線在AutoCAD軟件中生成三維格網;最后,將三維格網導入Rhino軟件,生成數字化三維地形模型(圖1)。數字化三維地形模型可作為有限元分析的基礎地形模型,后續(xù)還需加入區(qū)域地質巖石特性來構建有限元三維地質模型。

圖1 九江地震臺區(qū)域數字化三維地形模型

1.3 有限元三維地質建模

考慮到實際區(qū)域地質環(huán)境存在局部破裂或巖性不均一等情況,在進行有限元三維地質建模時對地形構造要素及巖性模塊進行均一化處理。根據江西省地礦局九一六大隊2003年繪制的廬山地質圖[9],將九江地震臺附近區(qū)域劃分為不同巖性模塊(圖2),并在有限元建模軟件中對不同巖性模塊設定對應的巖石物理力學參數[10-11](表1),構建有限元三維地質模型,以進行有限元模型數據解算。

表1 巖石物理力學參數

圖2 模型劃分與巖石的分布

2 有限元建模應變解算

有限元模型數據的解算需選取合理的網格單元,以避免出現因網格過密導致計算機運行崩潰或因網格過疏導致數據解算準確度不夠等情況。九江地震臺區(qū)域范圍的幾何模型規(guī)模為km級,故選取網格尺寸為1 000 m。有限元三維地質模型存在多個塊體,需明確塊體間的接觸關系,因無法確定不同塊體巖層間的摩擦系數和粘合系數,假設巖層間無滑脫現象,采用直接綁定接觸確定荷載的傳遞方式[12]。另外,有限元模型數據解算的關鍵是設置模型邊界條件,即對模型施加約束,并施加重力和熱載荷,為防止模型移動,最底端采用固定約束。施加的重力載荷設置為標準地球重力常數值9.806 6 m/s2,施加的熱荷載為溫度變化值。

2.1 熱應變解算

因九江地震臺區(qū)域環(huán)境溫度四季變化顯著,月均溫度變化在0～30 ℃之間,而廬山地區(qū)海拔為1 474 m,溫度受地形高度影響較大,底部至頂端有約9 ℃的溫差。因此,僅需設置有限元模型底部溫度,再將模型從低到頂均分為4個層次,每上升一個層次溫度遞降3 ℃。

在有限元模型數據解算過程中,設置4種模型底部溫度條件(0 ℃、10 ℃、20 ℃、30 ℃),并對不同高度施加不同溫度的熱荷載。4種溫度條件下有限元模型熱應變的解算結果(圖3)反映了九江地震臺區(qū)域環(huán)境因溫度造成的應變差異:溫度升高時,熱應變增大,山體應變不斷增大,處于拉張狀態(tài);溫度降低時,熱應變減小,山體應變不斷減小,處于壓縮狀態(tài);山體頂端溫度最低,其熱應變總是低于山體底部。

圖3 施加不同溫度載荷的熱應變

2.2 定向應變解算

有限元模型除了能求解熱應變外,還能求解區(qū)域任意方向的應變。設定模型底部的熱載荷溫度為30 ℃,NS向和EW向的應變分別如圖4(a)和4(b)所示。因地質構造和溫度的差異,不同位置及高度的應變不同。NS向與EW向應變的最大值總是出現在山體底部,最小值總是出現在山體頂端,其余部分的應變情況保持在一定范圍內,有限元模型定向應變均值的量級為10-5。

圖4 30 ℃溫度時的定向應變

3 有限元建模與實測應變分析

有限元建模的應變解算反映的是研究區(qū)域的應變分布,存在假設和簡化處理,而實測應變反映的則是臺站應變儀器處的點應變狀態(tài),具有基準性。如果2種應變結果存在較好的關聯性,說明有限元建模應變解算結果是可靠的。

3.1 模型熱應變與實測應變比較

有限元模型熱應變隨溫度及方位的變化而變化。根據2011～2021年九江地震臺實測的外界溫度(圖5(a)),均一化處理得到每月對應溫度,求解每月對應溫度下的熱應變(圖5(b)),并與實測NS向及EW向應變(圖5(c)和5(d ))進行比較。結果發(fā)現,有限元模型熱應變均值與實測應變值的量級大小皆為10-5;模型熱應變曲線與實測應變曲線存在基本相同的年周期變化,實測應變曲線的相位有所滯后。實測應變曲線年變周期不明顯且幅度較小,原因是實測應變的儀器處于山洞之中,存在消除溫度干擾的措施,洞內年變溫度小于1 ℃,儀器受溫度的影響很小。如果儀器與有限元模型具有相同的溫度環(huán)境,實測應變的年變周期形態(tài)會顯著接近模型熱應變,說明有限元模型熱應變與實測應變具有較好的一致性。

圖5 有限元模型熱應變與實測應變比較

3.2 模型定向應變與實測應變比較

利用有限元模型解算每月對應溫度下的應變,并與實測應變進行比較,結果如圖6所示。在有限元模型構建過程中,因模型簡化及巖性均一化處理,有限元建模應變計算結果僅為概略值,模型中每個點的應變各有不同,故有限元分析的應變解算取值為區(qū)間值。圖6中,有限元模型定向應變的均值量級為10-5,與實測應變觀測結果量級一致,且實測NS向及EW向應變完全在有限元分析的定向應變區(qū)間內;有限元模型定向應變解算區(qū)間值存在明顯的年周期變化,而實測應變客觀上也存在年周期變化,只是相位滯后較大。實測應變曲線的年變周期不明顯,且幅度較小,其原因與有限元模型熱應變和實測應變比較中的闡述相同,由此說明有限元模型的定向應變解算在一定程度上能反映區(qū)域應變狀態(tài)。

圖6 有限元模型定向應變與實測應變比較

4 結 語

本文基于九江地震臺實測溫度和應變數據,對有限元建模進行分析,結果表明:

1)有限元模型解算的熱應變均值和定向應變均值與九江地震臺實測應變結果具有相同的量級,皆為10-5。

2)有限元模型解算的熱應變和定向應變與實測應變的變化趨勢基本一致,均呈現出較為明顯的年周期變化特征。由于山洞內存在保溫措施,使實測應變的年變周期不明顯,且變化幅度較小;實測應變的相位相對有限元建模應變有所滯后,原因是外界溫度傳導至山洞內需要一定時間。

3)有限元建模應變與實測應變之間存在合乎邏輯的內在聯系,驗證了利用數字高程模型數據和巖性參數建立有限元三維地質模型以分析區(qū)域應變的方法具有可行性。