劉兵權(quán)黎立云衛(wèi)夢(mèng)希王博楠王之東

1.中國礦業(yè)大學(xué)(北京)力學(xué)與建筑工程學(xué)院,北京 100083;2.煤炭資源與安全開采國家重點(diǎn)實(shí)驗(yàn)室,北京 100083;3.北京科技大學(xué)土木與資源工程學(xué)院,北京 100083;4.北京礦冶研究總院,北京 100160;5.國電聯(lián)合動(dòng)力技術(shù)有限公司,北京 100039

隨著人類生產(chǎn)活動(dòng)的發(fā)展,地殼中各種礦產(chǎn)資源的開發(fā)不斷深入,煤礦開采深度已超過1 000 m,金屬礦開采深度甚至超過2 000 m,由深部高地應(yīng)力引起的巖爆、沖擊地壓等災(zāi)害表現(xiàn)得越來越劇烈[1],因而有必要對(duì)深部地應(yīng)力變化趨勢(shì)進(jìn)行理論探討,為深部工程地應(yīng)力預(yù)測(cè)提供參考。

為了獲得地應(yīng)力場(chǎng)的分布數(shù)據(jù),現(xiàn)有研究已提出很多高水平的測(cè)試方法[2],得到了許多寶貴的地應(yīng)力實(shí)測(cè)數(shù)據(jù),特別是深部地應(yīng)力數(shù)據(jù)。 近些年,地應(yīng)力數(shù)據(jù)大量更新,有必要建立更完善的理論模型。

關(guān)于地應(yīng)力分布的理論分析,很多學(xué)者進(jìn)行過研究[3-6]。 金尼克假說認(rèn)為,垂直地應(yīng)力為γH,水平地應(yīng)力可通過假定水平應(yīng)變?yōu)?,由胡克定理求得[7]。 Brown 等[8]基于多個(gè)國家和地區(qū)的地應(yīng)力數(shù)據(jù)(深度多在1.5 km 以內(nèi),最深約為2.9 km),給出了垂直地應(yīng)力與深度的線性回歸關(guān)系(由此得出平均重度為27 kN/m3)和側(cè)應(yīng)力系數(shù)的包絡(luò)線。 McCutchen[9]考慮了巖石自重,利用球殼模型求得了水平地應(yīng)力和徑向地應(yīng)力與深度的關(guān)系,給出了側(cè)應(yīng)力系數(shù),但沒有結(jié)合實(shí)測(cè)數(shù)據(jù)進(jìn)行擬合。也有研究利用巖石斷裂強(qiáng)度來預(yù)測(cè)水平地應(yīng)力的范圍[10]。 這些研究在工程中廣泛應(yīng)用[11-13],但大部分研究未考慮自重體積力隨深度變化的情況[14],也沒有考慮內(nèi)壓的影響。

影響地應(yīng)力分布的因素有很多,如地心引力、地殼內(nèi)部高壓力、構(gòu)造運(yùn)動(dòng)、地形地貌、巖漿侵入、地溫梯度、孔隙水壓力、地球自轉(zhuǎn)等[7]。 地殼深部完整性好于表面,淺層中由于板塊運(yùn)動(dòng)、地形地貌的不均勻性對(duì)深部地應(yīng)力的影響將逐漸減弱,應(yīng)力分布主要由地殼整體的宏觀力學(xué)條件決定。 因此,地應(yīng)力與深度的函數(shù)關(guān)系需對(duì)起主要作用的巖體自重和地殼內(nèi)部高壓力綜合分析得出,并通過大量實(shí)測(cè)數(shù)據(jù)擬合得到。 地溫梯度和孔隙水壓力引起的地應(yīng)力變化是一種規(guī)律性的靜水應(yīng)力變化,地球自轉(zhuǎn)引起的應(yīng)力變化也在全球呈現(xiàn)規(guī)律性分布[15],因此均可疊加到擬合公式中;巖漿侵入等局部非均勻因素的影響則需根據(jù)實(shí)地情況考慮。

本文基于前人的研究工作做出以下補(bǔ)充和拓展:對(duì)地殼使用球殼模型,同時(shí)考慮了隨深度變化的萬有引力體積力和地球內(nèi)部高壓力對(duì)地殼應(yīng)力變化趨勢(shì)的綜合影響,利用大量地應(yīng)力實(shí)測(cè)數(shù)據(jù)對(duì)理論公式進(jìn)行擬合,使之更好地應(yīng)用于深部工程。

1 地殼宏觀應(yīng)力場(chǎng)理論推導(dǎo)

地球可近似視為球體,其半徑為6 371 km。 地球的圈層結(jié)構(gòu)由外至內(nèi)可分為地殼、地幔、地核,其中地殼平均厚度約為17 km,主要由巖石構(gòu)成;上地幔主要成分為巖漿,具有一定的流動(dòng)性[16]。 在地殼尺度上,可以認(rèn)為它滿足彈性力學(xué)基本假設(shè),因此參照彈性力學(xué)中球殼應(yīng)力的解法來求得地殼應(yīng)力分布。 地殼所受荷載為其中各質(zhì)點(diǎn)所受的萬有引力fr以及內(nèi)部的高壓力p[17],如圖1 所示。 以地心為原點(diǎn)建立球坐標(biāo)系,地殼區(qū)域總應(yīng)力場(chǎng)為兩部分荷載引起的應(yīng)力場(chǎng)的疊加。

圖1 地殼受載模型示意圖Fig.1 Diagram of crustal loading model

由彈性力學(xué)可知,線彈性條件下受內(nèi)壓的無體力球殼的應(yīng)力分布為[18]

式中,σr為徑向應(yīng)力;σθ為環(huán)向應(yīng)力;R為地球平均半徑;r0為地殼與地幔交界處到地心的距離;r為該地殼單元體到地心的距離。

地殼內(nèi)部具有極高的內(nèi)壓力p[17]。 由式(1)可知,若僅考慮p的作用,在整個(gè)地殼厚度上,環(huán)向地應(yīng)力都為拉應(yīng)力。 這顯然與眾多地應(yīng)力實(shí)測(cè)數(shù)據(jù)不相符合。

在彈性力學(xué)理論中,對(duì)含有體力,尤其是考慮萬有引力這種與距離有關(guān)的變體力情形鮮有討論,彈性力學(xué)書籍中沒有現(xiàn)成公式可用,因而需要重新推導(dǎo)。 對(duì)于球殼,引入自重體積力的作用,可能會(huì)使殼體中的各項(xiàng)應(yīng)力分量都成為壓應(yīng)力。

根據(jù)彈性力學(xué),球?qū)ΨQ問題的平衡方程[18]為

如圖1 所示,根據(jù)引力的高斯定理,與地心O相距為r處的地殼單元體,其所受的萬有引力等于半徑為r的球體對(duì)它的萬有引力[14]。 鑒于地球的圈層結(jié)構(gòu),密度在同一圈層內(nèi)比較接近,在不同圈層之間相差很大,因此將地殼部分的密度視為常數(shù)ρ。半徑為r的球體的質(zhì)量M可以用地球總質(zhì)量M0與該球體以外的地殼質(zhì)量的差來表示。 因此體力可以取為

式中,G為萬有引力常數(shù)。

從式(3)可看出,體積力fr與r是非線性關(guān)系。

由球?qū)ΨQ問題的幾何方程和物理方程可以得到彈性方程為[18]

再將彈性方程式(4)和體力表達(dá)式(3)代入平衡方程式(2),可得

式(5)為常微分方程中典型的歐拉方程,其通解為

式中,A1、A2、A3、A4為待定系數(shù)。

其中,A1、A3可由式(6)代入式(5),r各冪次的系數(shù)對(duì)應(yīng)相等得出

A2、A4需通過地殼的邊界條件求解。

將式(6)代入彈性方程式(4),可以得到僅由變體力自重引起的應(yīng)力場(chǎng):

邊界條件為

內(nèi)壓的作用已包含在式(1)中,這里地殼內(nèi)表面的徑向應(yīng)力應(yīng)為0。

由式(8)、(9)可得

由此,可解出A2、A4:

將僅由變體力自重引起的應(yīng)力場(chǎng)式(8)與由內(nèi)壓引起的應(yīng)力場(chǎng)式(1)疊加,即可得到地殼內(nèi)部總的應(yīng)力場(chǎng)表達(dá)式:

式(12)中各系數(shù)Ai為與地殼有關(guān)的宏觀等效物理參數(shù)。 這些參數(shù)很難通過測(cè)量得到,但可以通過對(duì)現(xiàn)有的地應(yīng)力實(shí)測(cè)數(shù)據(jù)擬合而求得。

可以看出,在地球尺度上,地殼應(yīng)力分布與深度或半徑并不是嚴(yán)格的線性關(guān)系;水平地應(yīng)力與垂直地應(yīng)力的表達(dá)式共用參數(shù)Ai,因此二者相互關(guān)聯(lián)。

2 地殼應(yīng)力場(chǎng)與深度的關(guān)系擬合

上一節(jié)基于彈性力學(xué)的推導(dǎo),按彈性力學(xué)約定,應(yīng)力σ以受拉為正。 為方便與地應(yīng)力實(shí)測(cè)數(shù)據(jù)進(jìn)行分析與比較,本文此后將地應(yīng)力用s表示。按巖土力學(xué)約定,地應(yīng)力s以受壓為正。

在地應(yīng)力的工程測(cè)量中,通常用sv表示垂直主應(yīng)力,sH表示水平方向較大的主應(yīng)力,sh表示水平方向較小的主應(yīng)力,三者在多數(shù)情況下為正值,表現(xiàn)為壓應(yīng)力(在彈性理論中應(yīng)為負(fù)值)。 由于構(gòu)造等非均勻性因素的作用,這3 個(gè)主應(yīng)力的方向并不嚴(yán)格與鉛垂方向和水平方向相同,通常會(huì)有一個(gè)較小的傾角。 本文引用世界地應(yīng)力圖(WSM2016)[19]共404 條數(shù)據(jù)的3 向主應(yīng)力都有說明傾角,其中2 個(gè)水平主應(yīng)力傾角都不超過10°的數(shù)據(jù)共288 條,占71.3% 。 本文主要研究均勻性因素的作用,且多數(shù)傾角的數(shù)值較小、影響不大,在此忽略這個(gè)傾角的影響,將sH和sh取算術(shù)平均值sa,作為水平方向地應(yīng)力值,即

將式(12)中的多個(gè)常數(shù)合并,并結(jié)合式(13),得

Bi和Ai的關(guān)系可以通過對(duì)比式(12) ～式(14)得到。

不同地區(qū)、不同深度下的地應(yīng)力測(cè)量結(jié)果[19-31]見圖2、圖3 及附表1(篇幅所限,附表1 見于本刊網(wǎng)站)。 由于陸地地殼與海洋地殼的天然差異,并且現(xiàn)階段能搜集到的海洋地應(yīng)力數(shù)據(jù)較少,本文只引用陸地地應(yīng)力數(shù)據(jù)并對(duì)其擬合。 由圖1 可知

式中,H為深度,m。

通過式(15),可將測(cè)點(diǎn)到地心的距離r與深度H相互轉(zhuǎn)化。

以式(14)作為擬合公式,利用附表1 中大量的地應(yīng)力實(shí)測(cè)數(shù)據(jù),應(yīng)用最小二乘原理,通過MATLAB 計(jì)算得到參數(shù)Bi的值為

B1=-0.031 226,B2=1.083 0×107,B3=-8.372 5×1010,

B4=-4.628 0×1017,B5=0.119 45,B6=-2.670 4×1010。

最終得到垂直地應(yīng)力sv和水平向地應(yīng)力sa的數(shù)學(xué)表達(dá)式為

若將式(16)中r代換為H,會(huì)增加公式的復(fù)雜程度,為工程應(yīng)用方便,只需得到地應(yīng)力隨測(cè)點(diǎn)半徑r變化的表達(dá)式即可。

擬合結(jié)果表明,垂直地應(yīng)力sv和水平地應(yīng)力sa都隨r的變小而增大,即兩者都隨深度H的增大而增大。 垂直地應(yīng)力sv擬合曲線的標(biāo)準(zhǔn)偏差為6.44 MPa,擬合相關(guān)系數(shù)為0.898 4;水平地應(yīng)力sa的擬合曲線標(biāo)準(zhǔn)偏差為8.19 MPa,擬合相關(guān)系數(shù)為0.867 3。

3 擬合結(jié)果對(duì)比分析

地應(yīng)力的擬合結(jié)果對(duì)比如圖2 和圖3 所示。由水平地應(yīng)力與垂直地應(yīng)力的比值可以得出側(cè)應(yīng)力系數(shù)(圖4)。

圖3 水平地應(yīng)力對(duì)比Fig.3 Comparison of horizontal in-situ stress

3.1 地應(yīng)力擬合結(jié)果分析

圖2 展示了各文獻(xiàn)中垂直地應(yīng)力實(shí)測(cè)值svm、本文擬合的sv曲線和依據(jù)金尼克假說得到的svγ曲線(重度取γ=22 kN/m3)的對(duì)比。 可以發(fā)現(xiàn),在5 km 的深度之內(nèi)(r>6 366 km),垂直地應(yīng)力基本呈線性變化,與金尼克假說非常吻合;而Brown 和Hoek 使用線性擬合得到的平均重度為γ=27 kN/m3[8],這可能是擬合時(shí)引用的數(shù)據(jù)差異引起的結(jié)果。

圖2 垂直地應(yīng)力對(duì)比Fig.2 Comparison of vertical in-situ stress

圖3 為各文獻(xiàn)中水平地應(yīng)力實(shí)測(cè)平均值sam、本文擬合的sa曲線和依據(jù)金尼克假說得到的saγ曲線(重度取γ=22 kN/m3,側(cè)應(yīng)力系數(shù)取λ=0.9)的對(duì)比。 可以發(fā)現(xiàn),水平地應(yīng)力也大致呈線性變化。 本文所得平均水平地應(yīng)力與金尼克假說比較接近,淺部相對(duì)差異較大,深部相對(duì)差異較小。

3.2 側(cè)應(yīng)力系數(shù)對(duì)比分析

圖4 為各文獻(xiàn)中側(cè)應(yīng)力系數(shù)實(shí)測(cè)值λm、本文得出的平均側(cè)應(yīng)力系數(shù)λ與Brown 和Hoek 給出的側(cè)應(yīng)力系數(shù)的包絡(luò)線λmax、λmin的對(duì)比。 其中,側(cè)應(yīng)力系數(shù)是根據(jù)公式λ=sv/sa計(jì)算得出的。

圖4 側(cè)應(yīng)力系數(shù)對(duì)比Fig.4 Comparison of lateral pressure coefficient

由圖4 可以看出,深度超過約3 km(r＜6 368 km)時(shí),本文得到的平均側(cè)應(yīng)力系數(shù)已高于Brown和Hoek 給出的側(cè)應(yīng)力系數(shù)上限,更接近于1;實(shí)測(cè)數(shù)據(jù)也大多超出此上限,分布在本文得出的λ曲線兩側(cè)。 這可能是Brown 和Hoek 原文缺乏超過3 km 深度的實(shí)測(cè)數(shù)據(jù)而導(dǎo)致的。

深度小于70 m(r>6 370.3 km)時(shí),本文得到的平均側(cè)應(yīng)力系數(shù)會(huì)低于Brown 和Hoek 給出的側(cè)應(yīng)力系數(shù)下限,此時(shí)自重因素減弱而構(gòu)造等非均勻性因素顯著,水平應(yīng)力分布離散,平均側(cè)應(yīng)力系數(shù)不具有代表性。 深度處于70 m ～3 km(6 368 km＜r＜6 370.3 km)時(shí),平均側(cè)應(yīng)力系數(shù)落在上下限之間,和實(shí)際結(jié)果比較吻合。

圖2 ～圖4 中仍有一些較為離散的地應(yīng)力實(shí)測(cè)值和側(cè)應(yīng)力系數(shù)值,說明在實(shí)際工程中該數(shù)據(jù)來源區(qū)域的構(gòu)造、溫度、孔隙水壓力或其他因素比較顯著,需要疊加這些因素的影響。

由式(16)可以得到地殼應(yīng)力分布變化趨勢(shì)圖(圖5)。 可以看出,在地殼更深處,地應(yīng)力具有一定的非線性趨勢(shì),還需要基于大量深部地應(yīng)力數(shù)據(jù)進(jìn)一步研究。 內(nèi)壓會(huì)使得球殼內(nèi)部產(chǎn)生環(huán)向拉應(yīng)力,但通過擬合結(jié)果來看,自重等因素的作用使得平均地應(yīng)力始終為壓應(yīng)力。

圖5 地殼應(yīng)力變化趨勢(shì)示意圖Fig.5 Diagram of crustal stress trends

4 結(jié) 論

(1) 在地球尺度上,地殼應(yīng)力隨深度的分布不再是線性關(guān)系,其理論關(guān)系符合式(12),水平地應(yīng)力與垂直地應(yīng)力相互關(guān)聯(lián)。

(2) 垂直地應(yīng)力sv和水平地應(yīng)力sa都隨深度H的增加而增大,隨測(cè)點(diǎn)位置變化的擬合關(guān)系見式(16)。

(3) 在5 km 的深度之內(nèi),地應(yīng)力基本隨深度線性變化,地應(yīng)力的大小與金尼克假說在γ=22 kN/m3、λ=0.9 時(shí)的結(jié)果最為接近,僅淺部的平均水平地應(yīng)力差異較為明顯。

(4) 當(dāng)深度超過70 m 時(shí),側(cè)應(yīng)力系數(shù)逐漸趨近于1,比較符合實(shí)際數(shù)據(jù);當(dāng)深度超過3 km 時(shí),平均側(cè)應(yīng)力系數(shù)比Brown 和Hoek 給出包絡(luò)線更接近實(shí)測(cè)數(shù)據(jù);在近地表區(qū)域,側(cè)應(yīng)力系數(shù)分布離散,難以預(yù)估。

(5) 地殼內(nèi)部的高壓會(huì)引起水平方向產(chǎn)生拉應(yīng)力,而自重引起的壓應(yīng)力更高,使得整個(gè)地殼厚度內(nèi)的垂直地應(yīng)力和水平地應(yīng)力都為壓應(yīng)力。