李德海

（廈門(mén)地質(zhì)工程勘察院，福建廈門(mén) 361008）

隧道圍巖流變參數(shù)的粘彈性位移反演與驗(yàn)證

李德海

（廈門(mén)地質(zhì)工程勘察院，福建廈門(mén) 361008）

將解析反演和數(shù)值反演相結(jié)合，對(duì)試驗(yàn)毛洞的斷面形狀和所處應(yīng)力場(chǎng)做出適當(dāng)?shù)暮?jiǎn)化，利用對(duì)應(yīng)性原理求解出洞室位移的Kelvin粘彈性解，并以此去擬合試驗(yàn)洞不同位置處的實(shí)測(cè)位移而得到多組流變參數(shù);分別將其作為FLAC3D數(shù)值模型蠕變模式的輸入?yún)?shù)，得出各自對(duì)應(yīng)的實(shí)際斷面形狀和實(shí)際應(yīng)力場(chǎng)條件下試驗(yàn)洞的蠕變位移數(shù)值解。然后，利用BP神經(jīng)網(wǎng)絡(luò)，對(duì)各組流變參數(shù)和其對(duì)應(yīng)的試驗(yàn)洞同一位置處的蠕變位移數(shù)值解進(jìn)行訓(xùn)練，建立起兩者之間的非線性關(guān)系;利用訓(xùn)練好的網(wǎng)絡(luò)，依據(jù)實(shí)測(cè)蠕變位移值得出了圍巖流變參數(shù);最后，利用隧道實(shí)測(cè)位移數(shù)據(jù)對(duì)反演的流變參數(shù)進(jìn)行了數(shù)值正分析驗(yàn)證。

隧道;圍巖;流變;粘彈性位移;反分析;數(shù)值模擬

0 引言

由于尺寸效應(yīng)和工程因素的影響，用室內(nèi)蠕變?cè)囼?yàn)確定的巖樣的流變參數(shù)往往與隧道圍巖本身的流變參數(shù)存在一定的差別，因此不能直接將其作為數(shù)值計(jì)算的輸入?yún)?shù)。20世紀(jì)70年代發(fā)展起來(lái)的反分析法在某種意義上不僅具有消除巖樣尺寸效應(yīng)等因素影響的作用，而且由于實(shí)測(cè)物理量包含了工程因素的影響，使得以反分析所獲取的流變參數(shù)進(jìn)行實(shí)際工程巖體的流變分析能取得與現(xiàn)場(chǎng)較為一致的預(yù)測(cè)結(jié)果。位移反分析是目前反分析方法中應(yīng)用最多的一種。對(duì)于流變問(wèn)題，即是在反演的過(guò)程中，以現(xiàn)場(chǎng)流變?cè)囼?yàn)毛洞的蠕變位移為實(shí)測(cè)物理量去獲取圍巖的流變參數(shù)，其方式可分為2種:解析反演和數(shù)值反演。文獻(xiàn)［1～3］利用解析反演的方法，通過(guò)隧道蠕變位移解析解或由其推出的蠕變?nèi)崃繑M合出流變參數(shù)。文獻(xiàn)［4～6］則采用了數(shù)值反演的方法給出了圍巖流變參數(shù)值。然而兩種方法均有其缺點(diǎn):解析反演只適用于簡(jiǎn)單幾何形狀和邊界條件下的線彈性、線粘彈性和無(wú)支護(hù)隧道問(wèn)題;而數(shù)值反演則可能會(huì)因?yàn)榱繙y(cè)數(shù)據(jù)的離散性而陷入計(jì)算不收斂的困境，即便收斂，也可能會(huì)因?yàn)槟Ｐ蛷?fù)雜，節(jié)點(diǎn)多而迭代時(shí)間過(guò)長(zhǎng)。因此如何用反分析的方法迅速有效地獲取斷面形狀和應(yīng)力條件復(fù)雜的隧道圍巖流變參數(shù)是一項(xiàng)有意義的探討工作。本文提出的在簡(jiǎn)化條件下用解析反演得出的結(jié)果作為耦合BP神經(jīng)網(wǎng)絡(luò)的FLAC3D數(shù)值反演迭代初值的方法，可迅速有效地反演出圍巖的流變參數(shù)。

1 靜水壓力下圓形洞室的Kelvin粘彈性位移解

Kelvin模型是由彈簧和粘壺并聯(lián)而成的流變?cè)Ｐ停▓D1），它是典型的可用于描述衰減蠕變的粘彈性模型，其蠕變曲線見(jiàn)圖1。鑒于試驗(yàn)洞實(shí)測(cè)蠕變位移都是隨時(shí)間趨于穩(wěn)定的衰減型曲線，因此本文選用Kelvin模型對(duì)試驗(yàn)洞的粘彈性位移進(jìn)行解析。

圖1 Kelvin流變模型及其蠕變曲線

根據(jù)對(duì)應(yīng)性原理，粘彈性問(wèn)題的解可根據(jù)相應(yīng)的彈性問(wèn)題的解通過(guò)Laplace正逆變換獲得［7］。對(duì)于靜水壓力狀態(tài)下的圓形洞室，由開(kāi)挖引起的圓形洞室無(wú)支護(hù)徑向位移為［8］:

式中:E——彈性模量;γ——泊松比;R0——圓形洞室半徑;r——計(jì)算點(diǎn)距離圓心的距離;P0——地應(yīng)力。

對(duì)于Kelvin模型，流變本構(gòu)方程為:

式中:D——對(duì)時(shí)間的微分算子;η——Kelvin模型的粘滯系數(shù)。

查L(zhǎng)aplace逆變換表，對(duì)其進(jìn)行 Laplace逆變換:

將p0=1，q0=E，q1=η代入，并考慮到彈性模量和剪切模量的關(guān)系E=2G（1+γ），得洞室相應(yīng)于Kelvin模型的粘彈性無(wú)支護(hù)徑向位移為:

對(duì)于洞壁，r=R0，代入上式即可得洞壁處的粘彈性無(wú)支護(hù)徑向位移為:

2 試驗(yàn)洞FLAC3D數(shù)值模擬方案

2.1 概述

在某隧道中開(kāi)挖兩條橫洞，分別稱為左橫洞和右橫洞，并將其作為流變?cè)囼?yàn)洞。為了獲取圍巖的蠕變變形，采用如下的FLAC3D數(shù)值模擬方案:（1）模擬初始地應(yīng)力場(chǎng)，并將相應(yīng)的初始位移場(chǎng)清零;（2）模擬試驗(yàn)洞開(kāi)挖，得到開(kāi)挖后的應(yīng)力場(chǎng)，并將相應(yīng)的由開(kāi)挖引起的位移場(chǎng)清零;（3）打開(kāi)蠕變計(jì)算模式，模擬試驗(yàn)洞的蠕變變形，獲取相應(yīng)的蠕變位移場(chǎng)。

2.2 實(shí)測(cè)初始地應(yīng)力

試驗(yàn)洞所處地應(yīng)力環(huán)境為:最大水平主應(yīng)力SH為14～17 MPa，垂直于試驗(yàn)洞軸線方向;最小水平主應(yīng)力Sh為8～9 MPa，平行于試驗(yàn)洞軸線方向;用上覆巖層重力估算的垂直應(yīng)力SV為4～5 MPa，故模擬過(guò)程中最大、最小水平主應(yīng)力和垂直應(yīng)力的方向分別為x，y和z方向，如圖2所示。

2.3 試驗(yàn)洞FLAC3D模型建立

試驗(yàn)洞為直墻拱頂形式的斷面，跨度和高度均為2 m。依據(jù)隧道開(kāi)挖的影響范圍，計(jì)算邊界取為上下左右邊界距離拱的圓心均為6 m。另外，按平面應(yīng)變計(jì)算，沿隧道軸向取為0.1 m，網(wǎng)格劃分如圖2，模型共劃分246個(gè)單元，節(jié)點(diǎn)數(shù)506。

圖2 試驗(yàn)洞FLAC3D計(jì)算模型

位移邊界條件為:按平面應(yīng)變計(jì)算，隧道軸向無(wú)位移，左右邊界受水平方向的位移約束，下邊界受豎向位移約束。同時(shí)，為了模擬初始應(yīng)力場(chǎng)，左右邊界受構(gòu)造應(yīng)力的面力作用，上下邊界受重力的面力作用，邊界內(nèi)部受重力體力作用。

2.4 蠕變變形前的應(yīng)力場(chǎng)模擬

試驗(yàn)洞開(kāi)挖前要先模擬初始地應(yīng)力場(chǎng)。選用摩爾-庫(kù)侖彈塑性本構(gòu)關(guān)系，計(jì)算參數(shù)為:圍巖密度ρ=2000 kg/m3，力學(xué)參數(shù)按照現(xiàn)場(chǎng)大剪試驗(yàn)選取:彈性模量E=5.762×109Pa，泊松比μ=0.3，粘聚力c=2.1 ×105Pa，內(nèi)摩擦角 φ =27.5°，抗拉強(qiáng)度 2.63×106Pa。在FLAC3D中，采用體積模量K和剪切模量G，按下式計(jì)算:

模擬的試驗(yàn)洞初始地應(yīng)力場(chǎng)見(jiàn)圖3。由圖3可知，初始垂直地應(yīng)力為4～5 MPa，水平地應(yīng)力為14～15 MPa，與實(shí)測(cè)值基本一致。

圖3 初始地應(yīng)力云圖

在試驗(yàn)洞開(kāi)挖之后、蠕變變形之前，還要模擬試驗(yàn)洞開(kāi)挖后經(jīng)過(guò)瞬時(shí)彈塑性變形后的應(yīng)力場(chǎng)。采用相同的本構(gòu)模型和計(jì)算參數(shù)，可得到模擬的開(kāi)挖后應(yīng)力場(chǎng)，見(jiàn)圖4。

3 試驗(yàn)洞圍巖流變參數(shù)的粘彈性位移反分析

3.1 蠕變變形監(jiān)測(cè)點(diǎn)布置

圖4 開(kāi)挖后應(yīng)力云圖

為了量測(cè)試驗(yàn)洞圍巖的蠕變變形，在試驗(yàn)洞洞壁上布置洞周收斂和拱頂下沉監(jiān)測(cè)點(diǎn)，布置圖見(jiàn)圖5。對(duì)3個(gè)測(cè)點(diǎn)采用電子全站儀測(cè)定三維坐標(biāo)，然后換算出收斂及沉降變形，且認(rèn)為監(jiān)測(cè)到的變形僅為蠕變變形，而不包含開(kāi)挖后的瞬時(shí)變形。在建立的FLAC3D模型中，記錄相應(yīng)的拱頂、拱腳位置處節(jié)點(diǎn)的豎向和水平位移，如圖6。

圖5 試驗(yàn)洞測(cè)點(diǎn)布置圖

圖6 試驗(yàn)洞模型中記錄位移的節(jié)點(diǎn)位置

3.2 BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練樣本確定

3.2.1 輸出樣本——流變參數(shù)的確定將式（10）改寫(xiě)成:其中試驗(yàn)洞的等效半徑按下式計(jì)算［9］:

式中:R0——等效隧道開(kāi)挖輪廓半徑;h——斷面高度;B——隧道開(kāi)挖輪廓跨度之半。

計(jì)算得試驗(yàn)洞等效半徑為1667 mm。

根據(jù)式（11），分別對(duì)左、右橫洞的拱頂沉降、左右拱腳沉降和拱腳水平收斂用origin進(jìn)行擬合，即可得到8組BP神經(jīng)網(wǎng)絡(luò)輸出樣本。左、右橫洞的實(shí)測(cè)時(shí)態(tài)變形曲線見(jiàn)圖7和圖8。擬合的流變參數(shù)見(jiàn)表1。

圖7 右橫洞開(kāi)挖后90天內(nèi)實(shí)測(cè)變形曲線

圖8 左橫洞開(kāi)挖后90天內(nèi)實(shí)測(cè)變形曲線

表1 擬合的流變參數(shù)

根據(jù)第3節(jié)得出的FLAC3D模型（應(yīng)力場(chǎng)為試驗(yàn)洞開(kāi)挖并發(fā)生瞬時(shí)彈塑性變形后的應(yīng)力場(chǎng)，位移場(chǎng)清零），打開(kāi)蠕變計(jì)算模式，選用Burgers流變本構(gòu)模型，其中Maxwell部分的mshear和mviscosity不予賦值，Kelvin部分的kshear和kviscosity分別賦予表1中的8組G和η值，密度ρ仍取為2000 kg/m3，體積模量K仍取為4.8×109Pa，邊界條件不變，控制參數(shù)取為 lfob=1.0e-3，ufob=5.0e-3，lmul=1.01，umul=0.9，時(shí)間步取值按下述原則選取:（1）最大時(shí)間步maxdt≤η/G;（2）最小時(shí)間步mindt比maxdt小2～3個(gè)數(shù)量級(jí);（3）初始時(shí)間步dt=mindt。按此原則所取的各組樣本的時(shí)間步見(jiàn)表2。同時(shí)開(kāi)啟時(shí)間步的自動(dòng)調(diào)整功能，分別求取各組樣本對(duì)應(yīng)的第10、20和90天的拱頂沉降蠕變變形，提取圖6中的節(jié)點(diǎn)1、2、3的垂直位移和節(jié)點(diǎn)2、3的水平位移，得到各組輸出樣本對(duì)應(yīng)的輸入樣本，見(jiàn)圖9。

表2 各組樣本蠕變計(jì)算中時(shí)間步取值

圖9 輸入樣本——蠕變位移

3.2.2 輸入樣本——拱頂沉降蠕變變形的確定

限于篇幅，只給出第二組樣本對(duì)應(yīng)的拱頂位置處歷時(shí)90天后的沉降蠕變曲線，見(jiàn)圖10。

圖10 拱頂沉降蠕變90天內(nèi)蠕變曲線

3.2.3 神經(jīng)網(wǎng)絡(luò)反演圍巖的流變參數(shù)

用BP神經(jīng)網(wǎng)絡(luò)對(duì)樣本進(jìn)行訓(xùn)練以獲取輸入輸出樣本之間的高度非線性映射關(guān)系的技術(shù)已經(jīng)較為成熟。限于篇幅，本文不對(duì)BP神經(jīng)網(wǎng)絡(luò)進(jìn)行詳述，其工作原理參見(jiàn)文獻(xiàn)［10］。

比較圖7和圖8可知，右橫洞的實(shí)測(cè)變形歷時(shí)曲線更好地表現(xiàn)出了圍巖衰減蠕變的特性，其曲線形式與數(shù)值計(jì)算的結(jié)果吻合，故首先選用右橫洞的拱頂沉降實(shí)測(cè)變形反演圍巖的流變參數(shù)。

神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)計(jì)算參數(shù)對(duì)最終的訓(xùn)練結(jié)果有很大的影響。以圖9（a）數(shù)值計(jì)算的拱頂?shù)?0、20和90天的沉降值為實(shí)際輸入，以與之相應(yīng)的試驗(yàn)洞流變參數(shù)為實(shí)際輸出，分別對(duì)其進(jìn)行歸一化處理并將處理結(jié)果轉(zhuǎn)置后分別作為輸入層與輸出層，這是因?yàn)榭紤]到BP網(wǎng)絡(luò)輸入層節(jié)點(diǎn)數(shù)不至于過(guò)多以及輸出矩陣和輸入矩陣應(yīng)具有相同的列數(shù)（8列）。經(jīng)過(guò)反復(fù)調(diào)試，當(dāng)取隱含層節(jié)點(diǎn)數(shù)r=8、學(xué)習(xí)速率η=0.45、訓(xùn)練步數(shù)t=2000時(shí)，訓(xùn)練的效果較好。訓(xùn)練結(jié)束的時(shí)候，共經(jīng)歷了2000步運(yùn)算，此時(shí)訓(xùn)練誤差為3.98×10-15，符合要求。訓(xùn)練過(guò)程的誤差曲線見(jiàn)圖11。整個(gè)過(guò)程借助于MATLAB 7.0完成。網(wǎng)絡(luò)訓(xùn)練好以后，以圖7中實(shí)測(cè)的拱頂沉降為實(shí)際輸入，經(jīng)過(guò)反歸一化后即可反演試驗(yàn)洞圍巖的流變參數(shù)，反演結(jié)果見(jiàn)表3。

圖11 BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練誤差曲線

類似地，可以右橫洞的左拱腳沉降、右拱腳沉降和拱腳水平收斂為輸入樣本，輸出樣本仍取表1的值，對(duì)其進(jìn)行訓(xùn)練，將實(shí)測(cè)蠕變位移輸入各自訓(xùn)練好的網(wǎng)絡(luò)，即可得出反演的流變參數(shù)，見(jiàn)表3。

表3 反演的流變參數(shù)

反演結(jié)果表明，對(duì)于Kelvin模型，可認(rèn)為其剪切模量G在1.19×108Pa附近，粘滯系數(shù)η在（2.5～6）×108Pa·s之間。

4 反演結(jié)果的數(shù)值正分析驗(yàn)證

某隧道軸線與試驗(yàn)洞軸線垂直，斷面形式為曲墻拱頂仰拱式，其跨度為7.66 m，高度為8 m，依據(jù)隧道開(kāi)挖的影響范圍，計(jì)算邊界取為:左右邊界距斜井?dāng)嗝尕Q直中軸線20 m，上下邊界距斜井?dāng)嗝婀叭A心所在水平面20 m。另外，按平面應(yīng)變計(jì)算，沿隧道軸向取為1 m。網(wǎng)格劃分如圖12，模型共劃分1680個(gè)單元，節(jié)點(diǎn)數(shù)2541。

圖12 隧道FLAC3D計(jì)算模型

在建立的FLAC3D模型中，記錄圖13所示節(jié)點(diǎn)的水平位移，計(jì)算相應(yīng)的拱腳及墻腰位置處的水平收斂，并與實(shí)測(cè)的斜井拱腳及墻腰水平收斂對(duì)比，以驗(yàn)證反演流變參數(shù)的合理性。

采用相同的數(shù)值模擬方案，僅需在應(yīng)力邊界中將水平構(gòu)造應(yīng)力改為Sh=8 MPa。在獲得隧道開(kāi)挖后經(jīng)過(guò)彈塑性變形的應(yīng)力場(chǎng)后，打開(kāi)蠕變計(jì)算模式，根據(jù)上節(jié)反演的流變參數(shù)，取流變參數(shù)均值，即kshear=1.19 ×108Pa，kviscosity=3.76 ×1013Pa·s，并考慮隨著計(jì)算蠕變時(shí)間的不同，采用變時(shí)間步dt，見(jiàn)圖14。其他計(jì)算參數(shù)和控制參數(shù)不變，計(jì)算所得的拱腳和墻腰蠕變收斂值與實(shí)測(cè)值對(duì)比見(jiàn)圖15和圖16。

圖13 隧道模型中記錄位移的節(jié)點(diǎn)位置

圖14 蠕變時(shí)間步取值

圖15 拱腳水平收斂實(shí)測(cè)值與數(shù)值正分析值對(duì)比圖

5 結(jié)論

圖16 墻腰水平收斂實(shí)測(cè)值與數(shù)值正分析值對(duì)比圖

根據(jù)對(duì)應(yīng)性原理，推導(dǎo)了圓形洞室在靜水壓力狀態(tài)下的Kelvin粘彈性位移解析解，用此解析解初步反演得出簡(jiǎn)化條件下的圍巖流變參數(shù)，并將其作為耦合BP神經(jīng)網(wǎng)絡(luò)的FLAC3D數(shù)值反演迭代初值，一方面避免了迭代初值選取的隨意性，給出了一種有效的神經(jīng)網(wǎng)絡(luò)樣本確定方法，另一方面也加快了數(shù)值反演的速度和反演結(jié)果的精確性。將反演的結(jié)果用于數(shù)值正分析驗(yàn)證，結(jié)果與實(shí)測(cè)值吻合的很好，表明了該方法是合理的，能夠用于斷面形式和應(yīng)力條件復(fù)雜的隧道圍巖流變參數(shù)反分析。

另外，在使用FLAC3D蠕變模式模擬隧道蠕變變形時(shí)，計(jì)算所取的時(shí)間步與蠕變時(shí)間為線性關(guān)系時(shí)模擬的效果較好，可為以后的研究提供借鑒。

［1］ 薛琳.圓形隧道圍巖蠕變?nèi)崃康拇_定及粘彈性力學(xué)模型的識(shí)別［J］.巖石力學(xué)與工程學(xué)報(bào)，1993，12（4）:338-344.

［2］ 劉世君，徐衛(wèi)亞，邵建富.巖石粘彈性模型辨識(shí)及參數(shù)反演［J］.水力學(xué)報(bào)，2002，（6）:101-105.

［3］ 趙旭峰，孫鈞.擠壓性軟巖流變參數(shù)反演與本構(gòu)模型辨識(shí)［J］.鐵道工程學(xué)報(bào)，2008，（5）:5-8.

［4］ 李云鵬，王芝銀.粘彈性位移反分析的邊界元法［J］.西安礦業(yè)學(xué)院學(xué)報(bào)，1989，9（1）:17 –24.

［5］ 陳炳瑞，馮夏庭，丁秀麗，等.基于模式–遺傳–神經(jīng)網(wǎng)絡(luò)的流變參數(shù)反演［J］.巖石力學(xué)與工程學(xué)報(bào)，2005，24（4）:553-558.

［6］ 楊有海，王長(zhǎng)虹.考慮時(shí)空效應(yīng)的隧道工程黏彈性位移反分析［J］.地下空間與工程學(xué)報(bào)，2005，5（3）:468-472.

［7］ 周光泉，劉孝敏.粘彈性理論［M］.安徽合肥:中國(guó)科學(xué)技術(shù)大學(xué)出版社，1996.

［8］ 羅佑榮，唐輝明.巖體力學(xué)［M］.湖北武漢:中國(guó)地質(zhì)大學(xué)出版社，1999.

［9］ 王中文，方建勤，夏才初，等.考慮圍巖蠕變特性的隧道二襯合理支護(hù)時(shí)機(jī)確定方法［J］.巖石力學(xué)與工程學(xué)報(bào)，2010，29（S1）:3241-3246.

［10］ 崔志盛，金磊，趙凱.雙向八車(chē)道連拱隧道圍巖力學(xué)參數(shù)反演分析［J］.探礦工程（巖土鉆掘工程），2011，38（5）:65-69.

Back Analysis and Verification on Rheological Parameters of Tunnel Surrounding Rock through Visco-elastic Displacement

LI De-hai（Xiamen Geological Engineering Investigation Institute，Xiamen Fujian 361008，China）

Based on the correspondence principles，the analytical visco-elastic solution of round tunnel displacement under hydrostatic pressure is given using Kelvin constitutive law.After simplifying the testing tunnel as a round one subjected to hydrostatic pressure，the preliminary rheological parameters are obtained by using the analytical solution to fit the monitoring displacement data.A FLAC3Dmodel is then established in which the testing tunnel is of the real shape and is under the real ground stress.The preliminary rhelogical parameters are then adopted as parameters required under the creep option of FLAC3D，and the numerical solution of the test tunnel displacement is then given.With the help of BP neural network，the mapping relationship between the rheological parameters and the corresponding numerical solution of displacement are established.And the rheological parameters are then back analyzed by applying monitoring data to this trained network.Finally，the verification of the usefulness and reliability of the back-analyzed rheological parameters is given by implementing a normal calculation process to a tunnel.

tunnel;surrounding rock;rheology;visco-elastic displacement;back analysis;numerical simulation

U45

A

1672-7428（2012）02-0074-06

2011-08-03

李德海（1977-），男（漢族），福建漳平人，廈門(mén)地質(zhì)工程勘察院工程師、國(guó)家注冊(cè)巖土工程師，巖土工程專業(yè)，從事工程勘察、巖土設(shè)計(jì)的工程技術(shù)和管理工作，福建省廈門(mén)市蓮前西路192號(hào)，282588020@qq.com。