隧道圍巖流變參數的粘彈性位移反演與驗證

李德海

（廈門地質工程勘察院，福建廈門 361008）

隧道圍巖流變參數的粘彈性位移反演與驗證

李德海

（廈門地質工程勘察院，福建廈門 361008）

將解析反演和數值反演相結合，對試驗毛洞的斷面形狀和所處應力場做出適當的簡化，利用對應性原理求解出洞室位移的Kelvin粘彈性解，并以此去擬合試驗洞不同位置處的實測位移而得到多組流變參數;分別將其作為FLAC3D數值模型蠕變模式的輸入參數，得出各自對應的實際斷面形狀和實際應力場條件下試驗洞的蠕變位移數值解。然后，利用BP神經網絡，對各組流變參數和其對應的試驗洞同一位置處的蠕變位移數值解進行訓練，建立起兩者之間的非線性關系;利用訓練好的網絡，依據實測蠕變位移值得出了圍巖流變參數;最后，利用隧道實測位移數據對反演的流變參數進行了數值正分析驗證。

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

0 引言

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

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

Kelvin模型是由彈簧和粘壺并聯而成的流變元件模型（圖1），它是典型的可用于描述衰減蠕變的粘彈性模型，其蠕變曲線見圖1。鑒于試驗洞實測蠕變位移都是隨時間趨于穩定的衰減型曲線，因此本文選用Kelvin模型對試驗洞的粘彈性位移進行解析。

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

根據對應性原理，粘彈性問題的解可根據相應的彈性問題的解通過Laplace正逆變換獲得［7］。對于靜水壓力狀態下的圓形洞室，由開挖引起的圓形洞室無支護徑向位移為［8］:

式中:E——彈性模量;γ——泊松比;R0——圓形洞室半徑;r——計算點距離圓心的距離;P0——地應力。

對于Kelvin模型，流變本構方程為:

式中:D——對時間的微分算子;η——Kelvin模型的粘滯系數。

查Laplace逆變換表，對其進行 Laplace逆變換:

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

對于洞壁，r=R0，代入上式即可得洞壁處的粘彈性無支護徑向位移為:

2 試驗洞FLAC3D數值模擬方案

2.1 概述

在某隧道中開挖兩條橫洞，分別稱為左橫洞和右橫洞，并將其作為流變試驗洞。為了獲取圍巖的蠕變變形，采用如下的FLAC3D數值模擬方案:（1）模擬初始地應力場，并將相應的初始位移場清零;（2）模擬試驗洞開挖，得到開挖后的應力場，并將相應的由開挖引起的位移場清零;（3）打開蠕變計算模式，模擬試驗洞的蠕變變形，獲取相應的蠕變位移場。

2.2 實測初始地應力

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

2.3 試驗洞FLAC3D模型建立

試驗洞為直墻拱頂形式的斷面，跨度和高度均為2 m。依據隧道開挖的影響范圍，計算邊界取為上下左右邊界距離拱的圓心均為6 m。另外，按平面應變計算，沿隧道軸向取為0.1 m，網格劃分如圖2，模型共劃分246個單元，節點數506。

圖2 試驗洞FLAC3D計算模型

位移邊界條件為:按平面應變計算，隧道軸向無位移，左右邊界受水平方向的位移約束，下邊界受豎向位移約束。同時，為了模擬初始應力場，左右邊界受構造應力的面力作用，上下邊界受重力的面力作用，邊界內部受重力體力作用。

2.4 蠕變變形前的應力場模擬

試驗洞開挖前要先模擬初始地應力場。選用摩爾-庫侖彈塑性本構關系，計算參數為:圍巖密度ρ=2000 kg/m3，力學參數按照現場大剪試驗選取:彈性模量E=5.762×109Pa，泊松比μ=0.3，粘聚力c=2.1 ×105Pa，內摩擦角 φ =27.5°，抗拉強度 2.63×106Pa。在FLAC3D中，采用體積模量K和剪切模量G，按下式計算:

模擬的試驗洞初始地應力場見圖3。由圖3可知，初始垂直地應力為4～5 MPa，水平地應力為14～15 MPa，與實測值基本一致。

圖3 初始地應力云圖

在試驗洞開挖之后、蠕變變形之前，還要模擬試驗洞開挖后經過瞬時彈塑性變形后的應力場。采用相同的本構模型和計算參數，可得到模擬的開挖后應力場，見圖4。

3 試驗洞圍巖流變參數的粘彈性位移反分析

3.1 蠕變變形監測點布置

圖4 開挖后應力云圖

為了量測試驗洞圍巖的蠕變變形，在試驗洞洞壁上布置洞周收斂和拱頂下沉監測點，布置圖見圖5。對3個測點采用電子全站儀測定三維坐標，然后換算出收斂及沉降變形，且認為監測到的變形僅為蠕變變形，而不包含開挖后的瞬時變形。在建立的FLAC3D模型中，記錄相應的拱頂、拱腳位置處節點的豎向和水平位移，如圖6。

圖5 試驗洞測點布置圖

圖6 試驗洞模型中記錄位移的節點位置

3.2 BP神經網絡訓練樣本確定

3.2.1 輸出樣本——流變參數的確定將式（10）改寫成:其中試驗洞的等效半徑按下式計算［9］:

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

計算得試驗洞等效半徑為1667 mm。

根據式（11），分別對左、右橫洞的拱頂沉降、左右拱腳沉降和拱腳水平收斂用origin進行擬合，即可得到8組BP神經網絡輸出樣本。左、右橫洞的實測時態變形曲線見圖7和圖8。擬合的流變參數見表1。

圖7 右橫洞開挖后90天內實測變形曲線

圖8 左橫洞開挖后90天內實測變形曲線

表1 擬合的流變參數

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

表2 各組樣本蠕變計算中時間步取值

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

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

限于篇幅，只給出第二組樣本對應的拱頂位置處歷時90天后的沉降蠕變曲線，見圖10。

圖10 拱頂沉降蠕變90天內蠕變曲線

3.2.3 神經網絡反演圍巖的流變參數

用BP神經網絡對樣本進行訓練以獲取輸入輸出樣本之間的高度非線性映射關系的技術已經較為成熟。限于篇幅，本文不對BP神經網絡進行詳述，其工作原理參見文獻［10］。

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

神經網絡的結構計算參數對最終的訓練結果有很大的影響。以圖9（a）數值計算的拱頂第10、20和90天的沉降值為實際輸入，以與之相應的試驗洞流變參數為實際輸出，分別對其進行歸一化處理并將處理結果轉置后分別作為輸入層與輸出層，這是因為考慮到BP網絡輸入層節點數不至于過多以及輸出矩陣和輸入矩陣應具有相同的列數（8列）。經過反復調試，當取隱含層節點數r=8、學習速率η=0.45、訓練步數t=2000時，訓練的效果較好。訓練結束的時候，共經歷了2000步運算，此時訓練誤差為3.98×10-15，符合要求。訓練過程的誤差曲線見圖11。整個過程借助于MATLAB 7.0完成。網絡訓練好以后，以圖7中實測的拱頂沉降為實際輸入，經過反歸一化后即可反演試驗洞圍巖的流變參數，反演結果見表3。

圖11 BP神經網絡訓練誤差曲線

類似地，可以右橫洞的左拱腳沉降、右拱腳沉降和拱腳水平收斂為輸入樣本，輸出樣本仍取表1的值，對其進行訓練，將實測蠕變位移輸入各自訓練好的網絡，即可得出反演的流變參數，見表3。

表3 反演的流變參數

反演結果表明，對于Kelvin模型，可認為其剪切模量G在1.19×108Pa附近，粘滯系數η在（2.5～6）×108Pa·s之間。

4 反演結果的數值正分析驗證

某隧道軸線與試驗洞軸線垂直，斷面形式為曲墻拱頂仰拱式，其跨度為7.66 m，高度為8 m，依據隧道開挖的影響范圍，計算邊界取為:左右邊界距斜井斷面豎直中軸線20 m，上下邊界距斜井斷面拱圈圓心所在水平面20 m。另外，按平面應變計算，沿隧道軸向取為1 m。網格劃分如圖12，模型共劃分1680個單元，節點數2541。

圖12 隧道FLAC3D計算模型

在建立的FLAC3D模型中，記錄圖13所示節點的水平位移，計算相應的拱腳及墻腰位置處的水平收斂，并與實測的斜井拱腳及墻腰水平收斂對比，以驗證反演流變參數的合理性。

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

圖13 隧道模型中記錄位移的節點位置

圖14 蠕變時間步取值

圖15 拱腳水平收斂實測值與數值正分析值對比圖

5 結論

圖16 墻腰水平收斂實測值與數值正分析值對比圖

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

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

［1］ 薛琳.圓形隧道圍巖蠕變柔量的確定及粘彈性力學模型的識別［J］.巖石力學與工程學報，1993，12（4）:338-344.

［2］ 劉世君，徐衛亞，邵建富.巖石粘彈性模型辨識及參數反演［J］.水力學報，2002，（6）:101-105.

［3］ 趙旭峰，孫鈞.擠壓性軟巖流變參數反演與本構模型辨識［J］.鐵道工程學報，2008，（5）:5-8.

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

［5］ 陳炳瑞，馮夏庭，丁秀麗，等.基于模式–遺傳–神經網絡的流變參數反演［J］.巖石力學與工程學報，2005，24（4）:553-558.

［6］ 楊有海，王長虹.考慮時空效應的隧道工程黏彈性位移反分析［J］.地下空間與工程學報，2005，5（3）:468-472.

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

［8］ 羅佑榮，唐輝明.巖體力學［M］.湖北武漢:中國地質大學出版社，1999.

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

［10］ 崔志盛，金磊，趙凱.雙向八車道連拱隧道圍巖力學參數反演分析［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-），男（漢族），福建漳平人，廈門地質工程勘察院工程師、國家注冊巖土工程師，巖土工程專業，從事工程勘察、巖土設計的工程技術和管理工作，福建省廈門市蓮前西路192號，282588020@qq.com。