童 凱，杜圣賢，賈 超，狄勝同，馬玉驍

(1.山東大學(xué)土建與水利學(xué)院，山東 濟(jì)南 250061；2.山東大學(xué)海洋地質(zhì)與工程研究所，山東 青島 266237；3.山東省地質(zhì)科學(xué)研究院，山東 濟(jì)南 250013)

0 引 言

鹽巖因其易溶于水，所以水溶開采鹽礦是一種十分便捷經(jīng)濟(jì)的方案。又因具有極低的滲透率、良好的損傷自愈能力和塑性變形大等優(yōu)良特性，被全球各國廣泛的運(yùn)用在石油、天然氣等資源的儲(chǔ)備上。

在鹽巖力學(xué)特性研究方面，楊春和等[1]、YANG等[2]、WANG等[3]、姜德義等[4]、HOU等[5]學(xué)者開展了諸多研究，這對(duì)掌握鹽巖特性十分重要。由于鹽巖的流變特性，使得鹽腔在運(yùn)營過程中會(huì)發(fā)生收斂，進(jìn)而引發(fā)地面沉降，威脅到地表建筑物的安全[6]。由此可見研究引起地表沉降的機(jī)制對(duì)于控制由于鹽礦開采引發(fā)的次生地質(zhì)災(zāi)害十分必要。然而我國在地下礦產(chǎn)開采引起地表沉降的研究方面起步較晚，劉寶琛等[7]最先引進(jìn)并發(fā)展了概率積分法，并在煤礦地表變形預(yù)測中進(jìn)行了運(yùn)用。任松等[8]又將用于煤礦開采沉降計(jì)算中的概率積分法引入到鹽巖的沉降預(yù)測中，并根據(jù)該理論發(fā)展了新概率積分三維預(yù)測模型。李銀平等[9]利用半無限域內(nèi)球形空洞受力收縮位移解，利用疊加原理得到考慮內(nèi)壓作用下的地表沉降彈性積分解析解。陳雨等[10-11]采用Gaussian曲線表示沉降分布，結(jié)合腔體收斂函數(shù)，建立了傳遞函數(shù)法預(yù)測地表變形的理論在此基礎(chǔ)上綜合考慮影響腔體收斂的多種時(shí)效因素(腔體水溶速率、體積收斂系數(shù)以及收斂傳遞率等)，根據(jù)SCHOBER等[12]提出的兩階段收斂模型，建立了一套較為完整的變形動(dòng)態(tài)預(yù)測方法。

基于傳遞函數(shù)建立的動(dòng)態(tài)預(yù)測模型是一套較完整的用于求解鹽礦開采誘發(fā)地表沉降的計(jì)算模型，但模型中收斂系數(shù)和傳遞率對(duì)于地表沉降的動(dòng)態(tài)預(yù)測影響較大，其值估計(jì)的準(zhǔn)確性直接決定了地表沉降結(jié)果的準(zhǔn)確性。因此，本文基于傳遞函數(shù)方法和冪指函數(shù)模型，在考慮腔體收斂過程的情況下，對(duì)腔體收斂形式進(jìn)行合理的簡化，建立起一套適用于預(yù)測鹽巖開采誘發(fā)地面變形的動(dòng)態(tài)預(yù)測模型。

1 腔體收斂變形理論模型

1.1 傳遞函數(shù)模型

傳遞函數(shù)模型考慮地下腔體收斂與地表變形之間的聯(lián)系，認(rèn)為腔體臨空邊界所受應(yīng)力以水平方向?yàn)橹?，假設(shè)其受力過程中水平方向的變形遠(yuǎn)大于垂直方向的變形，從而建立起地下腔體收斂與地表變形之間的函數(shù)[10]，見式(1)。

式中：S(x,y)為地表沉降分布函數(shù)；a為考慮收斂體積經(jīng)上覆巖層傳遞折減后到達(dá)地面形成的沉降體積與溶腔體積之比；G(x,y)為聯(lián)系地下腔體收斂和地表變形的關(guān)系函數(shù)；q(z)腔體收斂函數(shù)；zt、zb分別為腔體頂部和底部埋深。

1.2 收斂函數(shù)公式推導(dǎo)

基于傳遞函數(shù)法建立的動(dòng)態(tài)預(yù)測模型可以運(yùn)用于求解地面沉降問題[10,12]，詳見式(2)～(4)。

(2)

F(ζ,η,t′)=

式中：S(d,t)為距離腔體中心d處在t時(shí)刻的沉降量；Sm為腔體中心處的最大沉降量；d為計(jì)算點(diǎn)到腔體中心的距離；r～2=rt×rb，rt、rb其中分別是溶腔頂部和底部的影響半徑；T為建腔時(shí)間；ζ為收斂系數(shù)，η為傳遞率。

該模型在進(jìn)行計(jì)算時(shí)，由于ζ、η參數(shù)的取值決定了計(jì)算結(jié)果的精度，而該參數(shù)又不易獲得，為減小這一影響，本文在此基礎(chǔ)上對(duì)腔體收斂過程進(jìn)行一定合理簡化，建立一套新的計(jì)算模型。

假設(shè)腔體建成后運(yùn)行t時(shí)刻時(shí)，腔體體積滿足式(5)關(guān)系。

V(t)=ψV0(5)

式中：V0為腔體建腔時(shí)體積；V(t)為t時(shí)刻收斂后的腔體體積；ψ為體積收斂率。

假設(shè)腔體收斂時(shí)主要發(fā)生水平方向的變形，而忽略垂直方向的變形。以圓柱形、橢圓形和梨形腔體為例，建立單腔變形預(yù)計(jì)模型，見圖1。利用圖1中的模型參數(shù)，根據(jù)腔體收斂的假設(shè)，便可以建立腔體的收斂函數(shù)q(z)。

圖1 腔體收斂模型Fig.1 Convergence model of salt cavern

橢圓形腔體收斂前后的橢圓可描述為式(6)～式(7)。

(6)

式中：Δr為腔體收斂半徑；B為橢圓形短半軸；A為橢圓形長半軸；zm為橢圓形腔體中心埋深。

將橢圓形腔體體積計(jì)算公式帶入式(5)可得式(8)。對(duì)于整個(gè)橢圓形，收斂體積可表達(dá)為式(9)。聯(lián)立式(6)～(9)，可得式(10)。

(8)

同理可得圓柱形腔體和梨形腔體的收斂函數(shù),分別表達(dá)為式(11)和式(12)。

q(z)=π(1-ψ)×

式中：Zm為梨形腔體上下腔分界線埋深；h1、h2分別為梨形腔體上下腔高度；Rp為梨形腔體上下腔分界線處半徑。

1.3 沉降函數(shù)推導(dǎo)

利用高斯分布函數(shù)描述地下收斂與地表變形的傳遞函數(shù)[10]，見式(13)。

(13)

式中：R為腔體開采在的地表影響半徑，R=zcotγ；d為沉降計(jì)算點(diǎn)到腔體中心在地表的水平投影距離；γ為收斂影響角，一般小于45°[12]。

定義式(14)，并聯(lián)立式(10)～(12)和式(14),可以求得三種不同形狀腔體的地表沉降預(yù)測公式。圓柱形腔體地表沉降預(yù)測表達(dá)式見式(15)；橢圓形腔體體地表沉降預(yù)測表達(dá)式見式(16)；梨形腔體體地表沉降預(yù)測表達(dá)式見式(17)。

(14)

1.4 地表沉降與水平變形關(guān)系

地表沉降和地表水平變形之間的關(guān)系可用式(18)表示[6]。

(18)

式中，n為水平變形影響系數(shù)，可由現(xiàn)場監(jiān)測數(shù)據(jù)所得，一般簡化計(jì)算取n=1。

考慮腔體收斂過程，并由此可以得到三種不同腔形開采下地表的水平變形計(jì)算公式。定義式(19)，圓柱形腔體的地表的水平變形計(jì)算見式(20)；橢圓形腔體的地表的水平變形計(jì)算見式(21)；梨形腔體的地表的水平變形計(jì)算見式(22)。

(19)

Uo(d)=

當(dāng)式(15)～(17)和式(20)～(22)中的ψ=0時(shí)，即此時(shí)儲(chǔ)氣庫完全收斂，地表變形達(dá)到最大值。故式(15)～(17)和式(20)～(22)可統(tǒng)一表達(dá)為式(23)。

S=S0(1-ψ)(23)

式中：S0為腔體完全收斂的地表變形值；S為當(dāng)腔體收斂率為ψ時(shí)的地表變形值。

至此，針對(duì)不同形狀單腔開采考慮腔體收斂過程的地表沉降變形預(yù)測公式已經(jīng)全部給出。

1.5 地表沉降動(dòng)態(tài)預(yù)測模型

地下采空區(qū)的存在勢必會(huì)對(duì)地表的建筑物和居民造成安全隱患[13]，而鹽巖的蠕變特性又使得地表變形隨著儲(chǔ)庫的運(yùn)行持續(xù)變化。在煤田開采條件下，地表移動(dòng)過程可從6個(gè)月延續(xù)到數(shù)年，而在鹽礦開采條件下，移動(dòng)持續(xù)時(shí)間甚至達(dá)到100年以上[14]。在儲(chǔ)庫運(yùn)行初期還很難看出由于鹽礦開采對(duì)地表建筑物的影響，但隨著時(shí)間的推移，地表變形部分受壓區(qū)可能會(huì)在移動(dòng)期間遭受拉伸，反之亦然。因此在儲(chǔ)庫建設(shè)過程中除了考慮穩(wěn)定后的地表變形狀態(tài)，還須考慮地表變形隨時(shí)間的發(fā)展過程[15]。

地表變形預(yù)測公式中，體積收斂率ψ隨著儲(chǔ)庫運(yùn)行時(shí)間的不同取值也不同，若知道不同時(shí)刻腔體的體積收斂率即可求得不同時(shí)間的地表變形分布。目前確定腔體收斂體積主要有以下三種方法：現(xiàn)場聲吶測腔、數(shù)值模擬和工程經(jīng)驗(yàn)預(yù)測。考慮到現(xiàn)場聲吶測試的花費(fèi)較高，又具有一定的滯后性，而經(jīng)驗(yàn)預(yù)估又具有很大的不確定性，因此在條件有限的情況下采用數(shù)值模擬進(jìn)行求解腔體收斂體積是一個(gè)較為可行的辦法。

2 地表變形算例分析

2.1 計(jì)算模型及參數(shù)說明

為了驗(yàn)證本模型的預(yù)測效果，筆者在Matlab中編寫了相應(yīng)的計(jì)算程序，用以求解地面變形分布，同時(shí)運(yùn)用Ansys建立相應(yīng)的有限元模型，在FLAC3D中進(jìn)行計(jì)算求解，將二者的計(jì)算結(jié)果進(jìn)行對(duì)比。

本次僅針對(duì)橢球體單腔開采引起的地面變形進(jìn)行預(yù)測，根據(jù)楊春和等[16]的研究成果，腔體采用長短軸為7∶3的橢球體模型。腔體取長半軸為70 m，短半軸比為30 m，腔體中心埋深1 000 m，各地層覆存深度及有限元模型見圖2。數(shù)值模擬的瞬時(shí)穩(wěn)態(tài)性計(jì)算采用Mohr-Coulomb彈塑性模型，蠕變計(jì)算采用FLAC3D自帶的Cpower蠕變模型。參考文獻(xiàn)[17]和文獻(xiàn)[18]取各地層巖性及蠕變參數(shù)，詳見表1和表2。

圖2 巖層分布及有限元模型Fig.2 Rock strata distribution and finite element model

表1 各巖層力學(xué)參數(shù)Table 1 Mechanical parameters of each rock layer

巖性彈性模量/GPa泊松比凝聚力/MPa內(nèi)摩擦角/(°)抗拉強(qiáng)度/MPa泥巖21.780.2322.0241.5鹽巖5.150.3001.5351.0砂巖20.000.2500.5350.5

表2 鹽巖蠕變參數(shù)Table 2 Creep parameters of salt rock

模型計(jì)算中相關(guān)參數(shù)參考之前學(xué)者研究和經(jīng)驗(yàn)進(jìn)行取值，沉降預(yù)測公式中的a取0.8，影響角取37°。數(shù)值模擬直接從腔體運(yùn)營期開始計(jì)算，暫不考慮腔體運(yùn)行壓力和建腔過程對(duì)于地表沉降的影響。假設(shè)腔體運(yùn)行壓力為0，雖然這將導(dǎo)致計(jì)算結(jié)果偏大，但對(duì)于儲(chǔ)庫設(shè)計(jì)考慮是偏安全的。

2.2 數(shù)值模擬與預(yù)測模型結(jié)果對(duì)比

通過數(shù)值模擬可以求得不同運(yùn)行時(shí)間的腔體體積收縮量，再帶入式(23)即可求得地表沉降變形的動(dòng)態(tài)過程。

2.2.1 地表沉降動(dòng)態(tài)曲線對(duì)比

根據(jù)數(shù)值模擬求出不同蠕變時(shí)間的腔體體積收縮量，再利用預(yù)測模型求取地表中心沉降值。以數(shù)值模擬的中心沉降值為標(biāo)準(zhǔn)，將預(yù)測模型和數(shù)值模擬的誤差統(tǒng)計(jì)在表3中。

表3 鹽腔地表中心沉降值計(jì)算表Table 3 Central surface subsidence of salt cavern

將數(shù)值模擬和預(yù)測模型的計(jì)算的不同蠕變時(shí)間的地表中心沉降值繪制在圖3中，可以更直觀地比較二者之間的關(guān)系。

從圖3和表3中可以看出，預(yù)測模型和數(shù)值模擬的發(fā)展趨勢基本相同。蠕變前10年二者相對(duì)誤差較大，由于該階段沉降值較小，即便絕對(duì)差值只有0.84 mm的第一年，相對(duì)誤差也達(dá)到了18.8%，隨著蠕變時(shí)間的增加，誤差逐漸減小，蠕變50年時(shí)，相對(duì)誤差僅有1.3%，整個(gè)計(jì)算周期內(nèi)，最大沉降差值為1.4 mm，說明本預(yù)測模型具有較高的精度。

2.2.2 地表變形空間分布對(duì)比

為了進(jìn)一步檢驗(yàn)預(yù)測模型的地表變形在空間中的分布情況，計(jì)算儲(chǔ)氣庫運(yùn)行50年時(shí)，以數(shù)值模擬的地表變形值為標(biāo)準(zhǔn)，將預(yù)測模型和數(shù)值模擬的誤差統(tǒng)計(jì)在表4中，將其與數(shù)值模擬的計(jì)算結(jié)果繪制在圖4中。

圖3 數(shù)值模擬和預(yù)測模型地表中心沉降對(duì)比Fig.3 Comparison of surface subsidence between numerical simulation and prediction model

表4 地表變形分布曲線對(duì)比Table 4 Comparison of surface deformation value

計(jì)算結(jié)果水平距離/m預(yù)測模型/mm數(shù)值模擬/mm相對(duì)誤差/%0-19.89-19.641.3100-19.54-19.142.1200-18.50-17.972.9300-16.93-16.403.2400-14.93-14.453.3地表沉降/mm500-12.70-12.342.9600-10.43-10.172.6700-8.26-8.735.4800-6.32-7.019.8900-4.66-5.6617.71 000-3.32-4.2622.100.000.20-1001.961.970.52003.713.904.93005.095.466.84005.996.457.1水平變形/mm5006.366.877.46006.276.676.07005.796.338.58005.065.7912.69004.205.0316.51 0003.324.1019.0

從圖4和表4中可以看出，數(shù)值模擬和預(yù)測模型計(jì)算的地表中心沉降最大值計(jì)算結(jié)果分別為19.64 mm和19.89 mm，相對(duì)誤差僅為1.3%。在整個(gè)計(jì)算區(qū)域內(nèi)，當(dāng)徑向距離小于600 m時(shí)，數(shù)值模擬的計(jì)算結(jié)果稍小于預(yù)測模型，但隨著距離的增大，數(shù)值模擬的結(jié)果會(huì)偏大一些。當(dāng)距離達(dá)到1 000 m時(shí)，由于地表沉降值已經(jīng)很小，即使只有0.94 mm的差值，相對(duì)誤差也達(dá)到了22.1%。

分析水平變形分布曲線和誤差統(tǒng)計(jì)可以得出，地表水平變形最大值發(fā)生在徑向距離為500 m的位置，二者計(jì)算的最大值分別為6.36 mm和6.87 mm，相對(duì)誤差7.4%，計(jì)算結(jié)果很接近。徑向距離小于700 m時(shí)，二者的相對(duì)誤差均小于10%，但隨著徑向距離的增大，很明顯數(shù)值模擬的計(jì)算結(jié)果較預(yù)測模型稍大，最大誤差達(dá)到了19%。

綜合地表沉降和水平變形分布曲線的對(duì)比結(jié)果，可以得出以下結(jié)論，本文所提出的預(yù)測模型計(jì)算的地表變形空間分布規(guī)律和數(shù)值模擬一致，但在距離腔體中心較遠(yuǎn)時(shí)，預(yù)測結(jié)果存在一定誤差。但主要關(guān)心的變形較大的集中區(qū)域擬合程度較好，因此在一定條件下能滿足工程精度需求。

圖4 儲(chǔ)氣庫運(yùn)行50年地表變形對(duì)比Fig.4 Comparison of surface deformation in operation for 50 years of gas reservoir

2.3 群腔地面變形動(dòng)態(tài)預(yù)測

為了保證地表變形在可接受范圍內(nèi)，盡可能多地開采鹽巖，充分利用鹽礦資源，使有限空間內(nèi)的經(jīng)濟(jì)效益達(dá)到最大化，實(shí)行群腔開采。運(yùn)用本文提出的地表變形動(dòng)態(tài)預(yù)測模型，再利用疊加原理進(jìn)行群腔的地面沉降動(dòng)態(tài)預(yù)測，假設(shè)各腔體沉降計(jì)算的參數(shù)相互獨(dú)立[12]，則地表任一點(diǎn)p(xi,yi)的變形值可由式(24)表達(dá)。

(24)

式中：Sp(xi,yi)(t)為點(diǎn)p在t時(shí)刻的地表沉降值；di為點(diǎn)p到i號(hào)腔體的徑向距離。

各地層分布情況及腔體形狀和圖2相同，開采形式如圖5所示，鄰腔夾角60°，腔距為2.0D[19]?；诏B加原理的假設(shè)可知，各腔體之間的變形互不影響，故以表3中的腔體收斂率作為該算例中各腔體的收斂動(dòng)態(tài)過程，將地表沉降結(jié)果繪制在圖6中。

圖5 群腔布置平面圖Fig.5 Plan position of salt caverns field

圖6 群腔開采地面沉降曲線Fig.6 The subsidence curve of salt caverns

采用圖5所示鹽群布置時(shí)，地面沉降最大值發(fā)生在腔群中心地表位置，在50年時(shí)達(dá)到了132.5 mm，是相同單腔開采條件下的6.66倍。由此可見，群腔開采對(duì)地面沉降的影響很大，因此需要對(duì)鹽巖開采規(guī)模和腔群布置方式進(jìn)行有效的控制。

3 結(jié) 論

1) 本文提出的預(yù)測模型所含未知參數(shù)較少，計(jì)算方便，在確定某工程的埋深及腔體形狀的前提下，可以通過預(yù)估或者參考類似工程儲(chǔ)庫運(yùn)行年份的腔體收斂情況，求得地表變形的動(dòng)態(tài)過程。根據(jù)不同的腔群布置形式可以很方便的求解地表沉降分布，從而進(jìn)行優(yōu)化設(shè)計(jì)。

2) 通過和數(shù)值模擬的計(jì)算結(jié)果對(duì)比，本模型在變形動(dòng)態(tài)預(yù)測和空間分布上具有較高的準(zhǔn)確性，地表中心最大沉降值的動(dòng)態(tài)預(yù)測差值最大為1.4 mm。當(dāng)徑向距離小于700 m時(shí)，地表變形的空間分布誤差均小于10%。腔群開采的沉降計(jì)算表明，多腔儲(chǔ)庫建設(shè)對(duì)地表沉降的影響顯著，所以在工程設(shè)計(jì)之初有必要針對(duì)腔群布置形式和建設(shè)規(guī)模進(jìn)行進(jìn)一步的研究和優(yōu)化。

由于缺少實(shí)際監(jiān)測資料，計(jì)算結(jié)果的準(zhǔn)確性有待驗(yàn)證，同時(shí)本模型未考慮運(yùn)行壓力以及造腔過程的影響，這也必然導(dǎo)致計(jì)算結(jié)果與實(shí)際情況存在一定偏差。模型計(jì)算中的參數(shù)還需根據(jù)實(shí)際儲(chǔ)庫運(yùn)行情況和腔體收斂結(jié)果進(jìn)行進(jìn)一步校核，以提高預(yù)測精度。