基于FLAC3D的巖體初始高水平構造應力場反演方法

郭延輝，侯克鵬

(1.云南農(nóng)業(yè)大學建筑工程學院，昆明 650201；2.昆明理工大學國土資源工程學院，昆明 650093)

地應力是地質(zhì)構造和自重共同作用在地質(zhì)體內(nèi)形成的原始應力，初始地應力的確定是工程分析與設計的基礎，對于工程的穩(wěn)定與安全至關重要[1]。礦區(qū)巖體地應力測量是提供地應力場的一種最為直接的方法，然而面對廣域的工程廠址，只能開展有限測點的地應力測量，而這些數(shù)據(jù)只能反映測量點處局部地應力場的情況，且有時會由于某些原因致使測量數(shù)據(jù)存在一定誤差[2-3]。為了能夠?qū)こ淘O計和施工提供更可靠的地應力數(shù)據(jù)，根據(jù)現(xiàn)場的地質(zhì)構造特征，結合有限的地應力實測數(shù)據(jù)，對初始應力場進行反演計算和分析，獲得更為準確、適用范圍較大的地應力場十分必要[4]。

近年來，研究人員提出了多種地應力場反演方法。李飛等[5]提出批數(shù)據(jù)處理(GMDH)神經(jīng)網(wǎng)絡算法，結合地應力場分布隨埋深的非線性特征及局部地質(zhì)構造處地應力場的非連續(xù)性特點，構建出形成復雜地質(zhì)體地應力場的邊界條件模式，并基于現(xiàn)場的稀少樣本測點數(shù)據(jù)對杏山鐵礦地應力場進行反演和重構；汪吉林等[6]提出以構造指標優(yōu)化人工神經(jīng)網(wǎng)絡(ANN)模型的建模方案，并開展了基于構造控制的地應力人工神經(jīng)網(wǎng)絡反演方法研究；易達等[7]研究了巖體初始應力場的遺傳算法與有限元聯(lián)合反演法；戚藍等[8]用系統(tǒng)工程灰色理論進行了地應力分析，在處理實測資料少的小樣本問題取得了較好的效果；李欣等[9]采用改進的粒子群優(yōu)化算法，對地應力場進行了反演計算，結果表明改進后的算法計算精度和效率更高；上述這些非線性反演方法雖然能夠反演分析復雜區(qū)域地應力場，卻存在收斂速度較慢、參數(shù)選擇過于依賴經(jīng)驗等問題，并可能導致同樣的研究對象而反演結果卻并不唯一的現(xiàn)象。數(shù)值模擬方法以其三維化、可視化、成本低和能夠計算復雜問題等優(yōu)點廣泛應用于深部地應力場的模擬計算[10-12]。然而，目前關于地應力場的反演研究大多是以自重應力場為主，針對初始高水平構造應力場的反演方法方面的研究仍然不足[13-14]。為此，以獅子山礦區(qū)為研究背景，在模擬邊界討論的基礎上，采用FLAC3D反演分析各種邊界條件下計算域內(nèi)水平構造應力與垂直應力的分布特征，最終采用基于初始應變能理論的方法，對獅子山礦區(qū)初始高水平構造地應力場進行反演及驗證，期望能對其他類似工程有一定的借鑒作用。

1 關于地應力反演方法的討論

通常數(shù)值計算過程中賦地應力的方法主要有：①模型內(nèi)部部分或某些點的地應力實測值已知，采用位移邊界條件，以一定的特征和規(guī)律將已知點應力差值分布到模型中，最后對應力場反演分析并生成最終應力場；②在計算域內(nèi)根據(jù)某種回歸方程直接計算每個單元體的初始應力，并轉(zhuǎn)化為計算單元的節(jié)點力，然后根據(jù)位移和應力邊界條件，平衡計算模型系統(tǒng)內(nèi)的應力場；③在模型位移邊界條件和應力邊界條件的基礎上生成初始應力場，然后根據(jù)幾個點實測的地應力對計算域的邊界條件進行調(diào)整，并最終擬合到新的初始應力場，使之與實測點的地應力值較為吻合；④基于初始應變能理論的初始地應力場的反演方法，其基本思路為：認為將有構造應力存在的地質(zhì)體作為部分具有彈性體屬性的意義賦上去考慮是合理的，因此在數(shù)值模擬反演初始地應力場時，可以將構造應力作為內(nèi)力施加在模型內(nèi)部，即對模型內(nèi)所有單元體采用三角形應力分布初值設置，而并非在邊界上施加，以此來模擬構造應力場[15-16]。

2 計算模型構建及方案

在應力場反演方法討論分析的基礎上，擬通過對比分析位移邊界條件、應力邊界條件、混合邊界條件等情況下，計算模型內(nèi)垂直應力與水平應力的分布情況，從而探討模擬高水平構造應力場的最佳方法。在實際數(shù)值模擬過程中，為了盡可能消除模型尺寸過小所產(chǎn)生的邊界效應，通常模型左右邊界向外擴展一定距離，計算作為平面應變問題考慮。通過反復試算，當計算模型長3 000 m，高1 400 m，寬50 m時，在這一尺寸下，邊界效應對模型的影響可以忽略。由于地表地形對模擬結果影響較小，因此，在模型的概化過程中，將地表概化成水平地表。計算采用Mohr-Coulomb本構模型。

為了充分研究邊界條件對初始地應力場分布的影響，計算模擬采用以下邊界條件[15-16]：①位移邊界條件，如圖1所示，即在數(shù)值模擬計算中，模型左右邊界的水平方向的位移固定，模型底部的位移固定，模型頂部設置為自由表面，其中G表示自重體積力；②應力邊界條件，圖2所示，模型左右施加應力梯形；③混合邊界條件，模型底部固定豎向位移，左右邊界分別為混合的三角形分布的水平分布面力和水平向位移約束；④基于初始應變能理論的初始地應力場的反演方法，其主要為位移邊界條件的基礎上，將構造應力作為內(nèi)力施加在模型內(nèi)部，即對模型內(nèi)所有單元體采用三角形應力分布初值設置。計算所采用的力學參數(shù)為青灰色白云巖的巖體力學參數(shù)，如表1所示。

圖1 位移邊界條件示意圖

圖2 應力邊界條件示意圖

3 應力場計算反演結果對比分析

3.1 基于位移邊界條件的初始地應力場的反演

在位移邊界條件下，反演得到垂直應力和水平應力分布分別如圖3所示。計算結果表明：應力等值線水平展布，走向穩(wěn)定，等值線間間隔幾乎相等，從地表至模型底邊呈線性分布。但在計算模型內(nèi)，同一位置的自重應力明顯大于水平應力，約為3倍左右。計算的最大主應力為豎直方向，垂直應力為最大主應力，表明此種方法僅能用于計算分析以自重應力場為主導的模型，而不能模擬高水平構造應力場。

圖3 位移邊界條件下的應力分布

3.2 基于應力邊界條件的初始地應力場的反演

一般情況下，水平地應力場是在自重體積力存在的條件下由構造運動產(chǎn)生的一種應力場，尤其是在高構造應力區(qū)，巖體的原巖應力狀態(tài)特征總是與顯著的構造應力狀態(tài)特征相聯(lián)系。大量的實測結果表明，水平構造應力場一般最大主壓應力σ1與最小主壓應力σ3方向均接近水平方向，而中間主應力σ2一般為豎直方向，在數(shù)值上等于γh(γ為重度，h為埋深)。在數(shù)值模擬反演中，自重應力通常作為一種體積力進行考慮并賦值，然而對于構造應力條件的處理辦法，特別是高水平構造應力場的反演方法存在一定的分歧。初始應力場模擬是指復原其采礦工程活動前(如開挖)的應力狀態(tài)，如果將構造應力作為邊界應力加載在模型邊界上，正如2節(jié)討論(圖2)，計算時以三角形應力施加在模型的左右邊界，并且保持左右邊界三角形應力不變，通過邊界單元與模型內(nèi)部單元接觸面的節(jié)點力向內(nèi)部傳遞，并且使得模型外表面上的分布力與其內(nèi)力平衡。

上述應力邊界條件下，模型的垂直應力分布如圖4(a)所示。由圖4(a)可以看出，垂直應力等值線圖在模型中部基本上呈水平走向，且線性增加，但模型兩端自重應力場明顯大于理論值，且應力等值線分布并不理想。圖4(b)為應力邊界條件下水平應力分布圖，在模型中部，同一位置水平應力高于自重應力，但模型上部和下部均未呈水平走向分布，并且模型中間偏下部分的最大主壓應力遠低于理論值，而模型中間偏上部分的最大主壓應力又遠高于理論值?？梢?，左右側邊界上采用應力邊界條件不能夠較好地計算模擬水平高構造應力場。

圖4 應力邊界條件下的應力分布

3.3 基于混合邊界條件的初始地應力場的反演

圖5(a)、圖5(b)分別為分布應力施加在模型左側時混合邊界條件下垂直和水平應力分布，圖6(a)、圖6(b)分別為分布應力施加在模型右側時混合邊界條件下垂直和水平應力分布。從圖5(a)、圖6(a)可以看出，遠離模型施加三角形分布面力的端部，自重應力等值線呈線性增加，一旦靠近該端部，自重應力等值線便發(fā)生變異，尤其是在施加分布面力一側，豎直應力等值線起伏較大，自重應力遠大于理論值。圖5(b)、圖6(b)中模型施加分布面力一側的較小范圍內(nèi)，模擬結果與預期結果較好，一旦離開該端部，水平應力等值線迅速變異。具體表現(xiàn)為在水平位移約束側淺部等值線過于密集，但在深部等值線又過于稀疏。礦體附近的水平應力低于預期結果，水平應力等值線越往深部越偏斜，等值線間隔在不同位置變化大，表明模擬得到的水平應力梯度變化不均勻。綜上分析可知，混合邊界條件既不適于模擬自重應力場，也不適于模擬水平構造應力場。

圖5 混合邊界條件下應力場分布(左側應力邊界、右側位移邊界)

圖6 混合邊界條件下應力場分布左位移邊界、右應力邊界

3.4 基于初始應變能理論的初始地應力場的反演分析

在地質(zhì)力學中，構造應力場通常是指地殼內(nèi)一定范圍內(nèi)某一瞬時的應力狀態(tài)。構造應力場中的應力分布和變化是連續(xù)并且有規(guī)律的。根據(jù)傳統(tǒng)的彈性理論和彈塑性理論，在自重體積力和構造運動的共同作用下，產(chǎn)生了位移并且做了相應功，在這一過程中，除塑性變形功不能轉(zhuǎn)化為應變能外，地質(zhì)體主要以彈性變形的方式儲存應變能。近年來，諸多學者從能量的角度和觀點出發(fā)，對巖石變形破壞過程中發(fā)生的能量積聚，耗散與釋放進行了分析[15]，使得能量觀點在巖體變形失穩(wěn)和破壞分析中的應用愈加廣泛。

假若將巖體視為彈性體，那么巖體中的某一點的微小單元的單位體積具有的應變能可表示為[16-17]

(1)

式(1)可改寫為

Gεijεij=WV+WF

(2)

式中:E為彈性模量；ν為泊松比；G為剪切模量；W為應變能；WV為體積變形應變能；WF為形狀改變應變能；θ為應變，θ=ε11+ε22+ε33。

如式(1)、式(2)所示，工程開挖后，初始應力場發(fā)生變化，二次應力場重新分布形成過程中，受到開挖影響的地質(zhì)體中某單元體將會同時出現(xiàn)體積改變和形狀改變。即微小單元體體積改變彈性比能WV和形狀改變彈性比能WF一般都會隨之變化。也就是說，巖體發(fā)生的變形破壞，其本質(zhì)是能量轉(zhuǎn)化的過程[13-15]，在工程開挖過程中，受影響的地質(zhì)體中所儲存的彈性應變能必然會隨之變化。

通過前面的討論，將有構造應力存在的地質(zhì)體作為部分具有彈性體屬性的意義上去考慮，在理論上是比較合理的。因此，在數(shù)值法反演初始地應力場時，可以將構造應力作為內(nèi)力施加在模型內(nèi)部，即對模型內(nèi)所有單元體采用梯形應力分布初值設置，以此來模擬高水平構造應力場，并采用位移邊界條件的情況下進行初始高水平構造應力場的反演計算。

圖7(a)為初始地應力場中垂直應力等值線分布。從圖7(a)可以看出，垂直應力從地表向下呈線性增加，應力跡線水平展布，走向穩(wěn)定，自重應力計算值與理論值較為吻合。圖7(b)為初始地應力場中水平壓應力等值線分布，水平應力從地表向下呈線性增加，應力跡線較平直，水平展布，等值線間隔也很均勻，走向穩(wěn)定，水平應力模擬結果與預期吻合較好。

圖7 初始應地應力場中的應力分布

圖8為初始地應力場主壓應力跡線分布。從局部放大圖上看到最大主應力均為水平方向，鉛直方向為中間主應力，滿足預期要實現(xiàn)的水平構造應力為最大主壓應力的初始應力場特征量的條件。

圖8 模型中部的主壓應力邊跡線分布

4 基于FLAC3D的獅子山礦區(qū)高水平構造應力場的反演與驗證

4.1 獅子山礦區(qū)地應力場測試結果

分別在礦區(qū)4個不同中段對應的4個測點進行了實測，通過對實測數(shù)據(jù)綜合分析，得出如下結果。

(1)4個測點的最大主應力均位于近水平方向，各測點的最大主應力方向均接近為NNW-SSE向(近垂直于礦體走向方向)，說明礦區(qū)以近水平構造應力場為主。

(2)各測點最大主應力、垂直主應力和最小主應力的大小均隨深度基本成線性增長，并且垂直主應力和最小主應力的隨埋深的增加速率要遠小于最大主應力的增加速率，進一步證實了礦區(qū)最大主應力的影響具有明顯的主導性。通過實測結果分析，得到最大主應力、中間主應力和最小主應力三個方向地應力(單位：MPa)隨深度H(單位：m)的線性回歸方程：

σh,max=-0.016 3+0.051 1H

(3)

σh,min=-0.388 3+0.021H

(4)

σv=-0.601 8+0.030 7H

(5)

式中：σh,max為最大主應力；σh,min為最小主應力；σv為垂直主應力。

(3)結合室內(nèi)巖石力學試驗結果和三維地應力實測結果，按照高地應力定量標準[18]，獅子山銅礦屬于高地應力區(qū)。

4.2 獅子山銅礦地應力場反演數(shù)值模型的構建與計算參數(shù)

計算模型以礦體為中心，計算域上邊界取至地表。計算模型X方向(垂直礦體走向方向)長度1 400 m；模型Y方向(沿礦體走向方向)，長度1 500 m；模型Z方向為豎直方向，模型底部標高687 m，模型最大高度1 430 m。模型共包括394 713個單元體和411 825個節(jié)點，計算模型和礦體形態(tài)分別如圖9、圖10所示。采用Mohr-Coulomb彈塑性本構模型[15]。計算模型X方向兩端約束X方向位移，模型Y方向兩端約束Y方向位移，模型底部固定位移，模型頂部為自由邊界。地應力按3.4節(jié)討論的高水平構造應力場的施加方法進行反演，即在位移邊界條件下，對模型內(nèi)所有單元體采用三角形應力分布初值設置，最大主應力為X方向，最小主應力為Y方向，豎直方向為中間主應力。礦區(qū)各巖性宏觀巖體力學參數(shù)如表1所示。

圖9 礦區(qū)三維有限差分數(shù)值模型

圖10 礦體形態(tài)

圖11分別為初始應力平衡后模型最大主應力和最小主應力分布。初始應力場反演計算表明，X方向(垂直于礦體走向方向)為最大主應力方向，Y方向(沿礦體走向方向)為最小主應力，Z方向(豎直方向)為中間主應力方向，與現(xiàn)場三維地應力實測主應力方向一致。根據(jù)三維地應力實測的4個測點位置，在模型內(nèi)部相應位置布設4個監(jiān)測點，以分析初始應力平衡后各監(jiān)測點的主應力數(shù)值與地應力實測值之間的關系。

圖11 初始應力平衡后模型主應力分布

為檢驗獅子山礦區(qū)地應力反演結果的可靠性，繪制實測地應力與反演地應力對比曲線，如圖12所示。各測點最大主應力反演值與實測地應力值之間的最大相對誤差為-8.10%，為測點D3；4個測點中，D1測點最小主應力反演值與實測地應力值之間的相對誤差最大，為-18.29%。各測點中，中間主應力反演值與實測地應力之間最大相對誤差為9.25%，為測點D3。綜上可見，該地應力場的反演方法具有一定的可靠性。

圖12 獅子山礦區(qū)4個測點實測地應力與反演地應力對比

5 結論

(1)在地應力邊界條件討論的基礎上，采用FLAC3D對4種地應力場反演方法和邊界條件進行數(shù)值模擬對比分析，結果表明位移邊界條件在模擬自重應力場中效果較好，但不能較好地模擬高水平構造應力場?；旌线吔鐥l件既不適于模擬自重應力場，也不適于模擬水平構造應力場。

(2)采用基于初始應變能理論的初始地應力場的反演方法，即在位移邊界條件下，對計算域內(nèi)所有單元體采用梯形應力分布初值設置，即設定初始彈性應變能狀態(tài)，使?jié)M足或擬合應力隨深度變化的實測值，模擬效果與實際吻合最好，能夠較好地模擬初始高水平構造應力場。

(3)建立了獅子山礦區(qū)復雜三維力學模型，以初始應變能理論的初始地應力場的反演方法進行設置，數(shù)值反演表明，各測點主應力的方向與實測地應力方向一致，各測點主應力的反演值與實測地應力回歸值之間的誤差在允許范圍內(nèi)，表明礦區(qū)水平構造應力場反演具有一定的可靠性。