有限元建模與實(shí)測應(yīng)變分析
——以九江地震臺為例

李 萌 盧和雄 趙愛平 嚴(yán) 麗,4

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

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

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

1.1 臺址環(huán)境

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

1.2 數(shù)字化三維地形建模

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

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

1.3 有限元三維地質(zhì)建模

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

表1 巖石物理力學(xué)參數(shù)

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

2 有限元建模應(yīng)變解算

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

2.1 熱應(yīng)變解算

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

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

圖3 施加不同溫度載荷的熱應(yīng)變

2.2 定向應(yīng)變解算

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

圖4 30 ℃溫度時(shí)的定向應(yīng)變

3 有限元建模與實(shí)測應(yīng)變分析

有限元建模的應(yīng)變解算反映的是研究區(qū)域的應(yīng)變分布,存在假設(shè)和簡化處理,而實(shí)測應(yīng)變反映的則是臺站應(yīng)變儀器處的點(diǎn)應(yīng)變狀態(tài),具有基準(zhǔn)性。如果2種應(yīng)變結(jié)果存在較好的關(guān)聯(lián)性,說明有限元建模應(yīng)變解算結(jié)果是可靠的。

3.1 模型熱應(yīng)變與實(shí)測應(yīng)變比較

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

圖5 有限元模型熱應(yīng)變與實(shí)測應(yīng)變比較

3.2 模型定向應(yīng)變與實(shí)測應(yīng)變比較

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

圖6 有限元模型定向應(yīng)變與實(shí)測應(yīng)變比較

4 結(jié) 語

本文基于九江地震臺實(shí)測溫度和應(yīng)變數(shù)據(jù),對有限元建模進(jìn)行分析,結(jié)果表明:

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

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

3)有限元建模應(yīng)變與實(shí)測應(yīng)變之間存在合乎邏輯的內(nèi)在聯(lián)系,驗(yàn)證了利用數(shù)字高程模型數(shù)據(jù)和巖性參數(shù)建立有限元三維地質(zhì)模型以分析區(qū)域應(yīng)變的方法具有可行性。