地形復雜三維地質模型下圍巖的蠕變特性研究

張益瑄,林文凱,晏啟祥,朱洪雪,劉 陽

(西南交通大學交通隧道工程教育部重點實驗室,成都 610031)

近些年來，隨著我國在交通隧道領域和地下工程建設方面的不斷發展，TBM法具有環保、掘進速度快等優點，在許多地下工程中被采用，而在TBM法施工過程中，圍巖蠕變對隧道開挖的影響越來越大，已經成為不可忽視的因素，近些年來，許多的學者研究了圍巖蠕變對隧道開挖的影響[1-9]，蔡燕燕等[10]考慮圍巖蠕變的全過程，推導出了深埋隧洞非線性位移解，侯公羽等[11]在支護反力變化條件下提出被動支護反力連續變化下的迭代算法，并且結合工程實例對圍巖的蠕變變形進行了分析。林文凱等[12]基于Burgers的蠕變模型對圓形隧道的內力進行了研究，對比分析了地層結構模型與荷載結構模型兩種分析方法，并得出了地層結構模型分析更為有效的結論。王明年等[13]基于有限元的分析方法，結合二郎山隧道的工程實例，分析了高地應力軟巖隧道的蠕變特性。同時，也有不少學者在蠕變位移方面進行了研究，余東明等[14]為了描述深埋圓形隧道的開挖支護，基于Burgers與Drucker-parger組合的黏彈塑性模型，推導出了深埋圓形隧道考慮剪脹性能的圍巖蠕變位移解析解，結合實際算例，分析了圍巖剪脹力與支護反力對深埋圓形圍巖蠕變位移的影響規律。

目前在國內，基于Burgers本構模型，利用有限元分析軟件ANSYS等建立荷載結構模型和地層結構模型，分析圍巖蠕變作用下隧道襯砌結構的內力的研究已經很多，本文做進一步研究分析，建立地形復雜，地表崎嶇不平的三維地質模型，在初始應力作用下，實施隧道開挖數值模擬，以隧道開挖模擬的結果為基礎，分析在圍巖蠕變作用下，支護結構拱頂、拱底及左右拱腰處的位移變化，從而得到圍巖的蠕變特性。

1 圍巖蠕變理論

綜合國內外的巖體蠕變理論研究，目前常用組合模型的方式來模擬巖體的流變關系。組合模型的基本原理是根據巖體的彈性、塑性及黏滯性質來設定一些基本元件，通過這些基本元件的串并聯組合，以求能夠反映巖體的本構關系，并以此模擬巖體的應力-應變關系。常用基本元件基礎上組合的模型有:麥克斯韋體(Maxweii)模型和開爾文體(Keivin)組合模型，以及由M體和K體串聯所組合的伯格斯(Burgers)(M-K)體組合模型[6]。結合高黎貢山隧道所處的砂巖層地段開展研究，根據地層巖性特點，選用黏彈性Burgers蠕變模型開展隧道圍巖蠕變分析。Burgers模型的本構方程為

σ+p1σ1+p2σ2=q1ε1+q2ε2

(1)

式中，E1、η1分別為Maxweii模型的彈性模量、黏滯系數；E2、η2分別為Keivin模型的彈性模量、黏滯系數；σ1、σ2分別為正應力對時間的一階、二階導數；ε1、ε2分別為正應變對時間的一階、二階導數。由于Burgers模型的4個參數不能作為Prony級數輸入值，因此需將E1、E2、η1、η2轉化為剪切模量的Prony級數參數[15]。對于Burgers模型，Prony剪切模量表達式為[12]

(2)

式中，g∞為松弛時間無限大時的剪切模量；g1，τ1分別為Maxweii模型的相對剪切模量和松弛時間；g2，τ2分別為Keivin模型的相對剪切模量和剪切時間；G0為松弛時間為0的剪切模量。

將Burgers模型中的彈性模量系數進行轉換，并將獲得的表達式進行拉普拉斯變換和拉普拉斯逆變換，便可以得到剪切模量的Prony級數表達式，通過對比式(2)，便可得到Burgers模型在ANSYS中需要的4個Prony級數g1，g2，τ1，τ2[16-17]，可用于下文的蠕變分析。

2 三維地質模型與應力場回歸方法

建立高黎貢山隧道砂巖層地段三維地質模型。為了模擬地形復雜，地表崎嶇不平的模型，結合地質平面圖中的等高線資料，導入CAD中得到dxf文件，然后用GoCAD軟件調用dxf文件，提取等高線上的相關擬合點坐標數據。由于數據的不均勻性，為了便于這些提取數據在ANSYS中建模，借用Surfer軟件使數據均勻化，使數據達到均勻值，然后將處理后的數據導入ANSYS中，就可以得到地形復雜，崎嶇不平的三維地形地質模型。結合高黎貢山實際條件，建立三維地質山體模型，模型下部是一個長方體模型，長、寬、高尺寸為250 m×150 m×120 m，在山體下部模型中包含有隧道襯砌結構模型，其中模型下部位置進行隧道開挖模擬，襯砌的內、外半徑分別為4、4.5 m。從而獲得ANSYS分析的整體模型如圖1所示。

圖1 數值模型

將圍巖視為各向同性連續介質，隧道圍巖、襯砌力學參數見表1。

表1 周圍土體與襯砌參數

圍巖蠕變參數分別為：Maxweii模型的彈性模量、黏滯系數分別為E1=5.546 GPa，η1=32.697 GPa·h，E2=3.211 GPa，η2=2.378 GPa·h，由前面公式可得到4個Prony級數，分別為g1=0.118 2，g2=0.883 6，τ1=17.265 4，τ2=0.267 6。

地應力場是隧道從開挖到支護等仿真分析的計算基礎，對合理進行地下工程設計與施工具有很大的現實意義[16-17]，初始地應力場也是影響隧道圍巖穩定性的重要影響因素[18]，以現場部分實測地應力數據為基礎，同時考慮地形、地質方面的因素，結合理論研究與數值分析進行反演[19-20]。本文考慮用多元線性回歸分析法對高黎貢山隧道砂巖層地段進行初始地應力回歸。

(3)

為了方便統計與計算，取m個相關點，由最小二乘法殘差平方和的方法得到

(4)

由最小二乘法的基本原理，可以得到回歸系數，為使得S殘最小，令S殘對Li的偏導數為0，即

(5)

整理后可以得到如下矩陣

(6)

(7)

j=1，2，3，…，6表示6個初始應力分量。

山體原地應力測試孔CS1、CS2分別位于模型x=40 m，y=75 m和x=140 m，y=75 m兩處附近，兩應力測試孔的實測數據如表2所示，其中數據已轉化為模型坐標方向的應力分量。

進行初始地應力場的反演，分別對該復雜地形三維模型施加構造應力場，應力場1、3分別為YZ面擠壓應力場與豎向的剪切應力場，應力場2、3、5分別為XZ面的擠壓應力場，水平與豎向剪切應力場。

用相關軟件分別計算上述各構造應力場和自重應力場，可以得到相關點深度的應力分量值，以表2實際測量得到的應力值為基礎，代入式(3)～式(7)進行分析，計算可以得到該處圍巖在上述假設各構造應力場下的初始狀態，具體如式(8)回歸方程

表2 測試孔應力分量

σ=1.073σ自重+6.374σ構1+1.628σ構2+
42.336σ構3+0.433σ構4+0.245σ構5

(8)

σ構1、σ構2、σ構3、σ構4、σ構5、σ自重分別對應5個構造應力場和自重應力場。對該處山體的回歸效果進行分析，可以計算得到該處山體的初始應力場回歸方程的相關系數r=0.964，接近于1.0。殘差平方和S殘=12.687，說明初始地應力場回歸具有較好的相似度。

3 結果分析

3.1 復雜地形下構造應力場

運用多元線性回歸法進行分析，得到山體在構造應力場以及自重應力場下的回歸方程，以式(8)和相關數據為基礎，可以得到相應的邊界條件，運用有限元軟件分析，復雜地形下山體地層的豎向位移和豎向應力如圖2、圖3所示。

圖2 初始地層豎向位移云圖(單位：mm)

圖3 初始地層豎向應力云圖(單位：Pa)

由圖2和圖3可知，受到山體地形復雜、崎嶇不平的影響，山頭地表處有較大的初始豎向位移，且在底層位置山體有較大的初始豎向應力。

3.2 隧道開挖支護狀態

為了得到復雜地形下隧道開挖后的支護狀態，用有限元軟件進行隧道開挖數值模擬，在上文假定的構造應力場與自重應力場作用下，隧道開挖完支護結構的豎向位移云圖見圖4。

圖4 隧道開挖支護襯砌豎向位移云圖(單位：mm)

從圖4可以看出，受到山頂復雜地形的影響，隧道縱向的中部范圍內豎向位移較大，在隧道頂部按照梭形分布，在靠近隧道兩端的位置，豎向位移逐漸趨于平緩，這種分布與圖2初始地層豎向位移較為類似。

隧道開挖支護結構豎向應力云圖如圖5所示。

圖5 隧道開挖支護襯砌豎向應力云圖(單位：Pa)

由圖5可以看出，隧道頂部沿縱向有較大的應力分布，且應力分布比較均勻，隧道頂部的應力分布與山頭處的應力分布較為相似。

3.3 圍巖蠕變下支護結構行為分析

下文研究地形復雜三維地質模型下，隧道開挖后，圍巖的蠕變對支護結構的影響。為了更好地得到結果進行分析，在建模過程中，取隧道模型的縱向長度為300 m，在模型的縱向方向，分別選取位置不同，埋深不同的5個數值采集面，并按照順序進行編號，數值檢測斷面的位置與埋深見表3。

表3 各數值采集斷面分布

其中，斷面1和斷面5位于模型縱向的兩側，位置上位于隧道埋設較淺處，斷面2、斷面3、斷面4位于模型縱向中部，相對于斷面1、斷面5，埋深較大。各數值采集斷面分布如圖6所示。

支護結構的拱頂部位也設置相應隧道縱向的數值采集斷面，各支護結構數值采集斷面示意見圖6。支護結構橫斷面監測點布置見圖7，如監測斷面2的拱頂為A2，拱底為B2，左側拱腰為C2，右側拱腰為D2。

圖6 數值監測點斷面示意

圖7 監測點橫向示意

以5個隧道橫向數值采集斷面和隧道縱向數值采集斷面的數據為基礎，分析研究隧道在圍巖蠕變荷載作用下，支護結構的響應特點。

3.3.1隧道各監測斷面Ai處的位移

以隧道5個橫向數值采集斷面支護結構Ai處位移隨圍巖蠕變變化的數據為基礎，得到Ai處位移隨時間變化的趨勢如圖8所示。

圖8 隧道各數值采集斷面Ai位移曲線

由圖8可知，在圍巖蠕變荷載作用下，位于隧道中部的斷面2、斷面3、斷面4橫向數值采集斷面的Ai處位移不斷增加，而位于隧道兩端的1、5數值采集斷面隨著時間的增加位移不斷減小。

圖8中，圍巖蠕變產生的初始階段，A3處沉降為38.71 mm，是產生最大沉降位置，A1處沉降為28.87 mm，是所有拱頂位置處的最小沉降值，這與在初始應力場作用下支護結構的位移分布是相似的，因為A3相比于A1埋深較大。由于A1相對于A5有更大的埋深，A1處的位移減小率更大，在斷面2、斷面3、斷面4中，在圍巖蠕變作用下，可以發現A3處位移增加的速率是最快的，原因是斷面3相對于其他斷面埋深較大。

3.3.2隧道各數值采集斷面Bi處的位移

通過對5個橫向數值采集斷面得到的數據進行匯總分析，得到B1處位移隨時間變化的趨勢如圖9所示。

圖9 隧道各數值采集斷面Bi處位移曲線

由圖9可知，受到蠕變荷載的影響，1與5數值采集斷面支護結構Bi的位移隨著天數增加不斷減小， 隧道斷面2、斷面3、斷面4橫向數值采集斷面的B1處位移是不斷增加的。Bi處的位移變化趨勢與Ai處的位移變化趨勢是類似的。

在圖9中，圍巖蠕變產生的初始階段，B3處位移的大小為22.01 mm，是產生位移的最大位置，B1沉降為15.98 mm，是所有拱頂位置處沉降的最小值。同時可以發現，在斷面2、斷面3、斷面4中，斷面3擁有較大的埋深，因此A3處位移最大，在2、3、4斷面中，在圍巖蠕變隨時間作用下，埋深越大的位置處的位移增長速率也越大。

3.3.3隧道各監測斷面Ci與Di的收斂位移

通過對5個橫向數值采集斷面得到的數據進行匯總分析，得到Ci與Di的收斂位移隨時間變化的趨勢如圖10所示。

圖10 隧道各監測斷面Ci與Di的收斂位移曲線

由圖10可以看出，隧道橫向數值采集斷面2、斷面3、斷面4隨著時間增加，左右拱腰收斂位移呈減小趨勢，其中埋深較大的斷面3處收斂位移減小最大，斷面1和斷面5處左右拱腰收斂位移呈現增大趨勢，由于斷面1和斷面5埋深接近，增長趨勢也較為接近。在圍巖蠕變作用下，通過對比分析Ai與Bi處的位移，可以發現Ci與Di之間的收斂位移的變化趨勢與拱頂和拱底處的位移變化趨勢恰好是相反的。

在圖10中，圍巖蠕變產生的初始階段，Ci與Di位置的收斂位移4.62 mm是收斂位移最大值處，相對地，C3與D3位置的收斂位移大小為3.89 mm，是收斂位移的最小值處。由圖中各曲線的變化趨勢，在圍巖蠕變的作用下，隨著斷面埋深的不斷增大，相應左右拱腰的收斂位移不斷減小，與上述各斷面Ai、Bi處位移隨蠕變變化相比，Ci與Di位置的收斂位移數值變化幅度較小，并且變化的數值也比較小。

3.3.4 隧道縱向各數值采集斷面拱頂位移分析

為了具體分析研究縱向拱頂位移隨蠕變的變化趨勢，這里沿著隧道的縱向方向，在拱頂位置取一系列采集點，采集點的范圍為10～240 m，每隔10 m共取24個點，研究蠕變作用天數下位移的變化趨勢，結果如圖11所示。

圖11 隧道橫向數值采集斷面拱底位移曲線

結合圖11進行分析，在地形崎嶇復雜的山體中間部位，沿著隧道縱向進行分析，相對埋深較大處，支護結構的拱頂位移大于山體兩邊處的拱頂位移。隨著圍巖蠕變作用不斷增加，將圖11根據變化趨勢分為3個區間，其中區間I，在隧道縱向x=10 m到x=60 m范圍內，支護結構拱頂位移隨著蠕變時間的增加不斷減小；區間Ⅱ，沿隧道縱向區間x=60 m到x=200 m，支護結構拱頂位移在蠕變荷載作用下，隨著時間是不斷增加的；區間Ⅲ，沿隧道縱向區間x=200 m到x=240 m范圍內，隨著荷載作用時間的增加，拱頂位移不斷減小，區間Ⅱ范圍內形成了一個拱形曲線，觀察圖2發現，支護結構拱頂位移隨著時間逐漸形成一個山頭，與地形崎嶇復雜的三維地質模型的地表形狀是類似的。

其中，區間Ⅰ處的拱頂位移隨著蠕變時間增加，變化幅度大于區間Ⅲ的位移變化，在圖11中表現為區間Ⅰ處的曲線比區間Ⅲ的曲線更為稀疏，在三維地質模型中表現為區間Ⅰ處的地表地形更加復雜陡峭。

3.3.5同一數值采集斷面下支護結構位移隨時間變化分析

為了研究圍巖蠕變對支護結構的響應，選取上文的5個監測斷面，分別提取它們在Ai，Bi處的位移以及Ci，Di兩處的收斂位移進行分析。

數值采集斷面1支護結構的位移隨時間變化如圖12所示。

圖12 數值采集斷面1支護結構位移曲線

由圖12可以看出，數值采集斷面1在整個蠕變分析過程中，拱頂處的位移要大于拱底處的位移，拱頂處和拱底處的位移都隨蠕變時間的增加而減小，但減小的速率不斷增大，左右拱腰處的位移隨時間緩慢增加，增加速率較為穩定。

數值采集斷面3支護結構的位移隨時間變化曲線如圖13所示。

圖13 數值采集斷面3支護結構位移曲線

由圖13可以看出，數值采集斷面3在蠕變時間增加的過程中，拱頂的位移要大于拱底的位移，拱頂和拱底位移是不斷增加的，增加的趨勢較為穩定，而左右拱腰的位移隨蠕變時間增加不斷減小，整體上監測斷面3左右拱腰的位移較小。

數值采集斷面5支護結構的位移隨時間變化如圖14所示。

圖14 數值采集斷面5支護結構位移曲線

由圖14可以看出，在蠕變時間增加的過程中，拱頂位移數值要大于拱底位移，但拱頂位移和拱底位移都隨時間不斷減小，左右拱腰位移隨蠕變時間呈平緩增加的趨勢。

4 結論

通過建立地形復雜，崎嶇不平的三維地質模型，在初始應力作用下，實施隧道開挖數值模擬，分析隧道開挖模擬后在圍巖蠕變作用下，支護結構不同部位的位移變化，對支護結構5個橫向數值采集斷面的Ai，Bi，Ci，Di處，以及1個縱向數值采集斷面進行相關數據提取，從而得出以下結論。

(1)在山體中間部位，支護結構的拱頂Ai處的位移在蠕變荷載作用下，位移隨時間變化的趨勢是不相同的，其中，埋深相對較淺的斷面A1和A5處拱頂位移在圍巖蠕變荷載作用下，位移呈逐漸減小的趨勢，而位于隧道中部的A2，A3，A4隨著蠕變時間的增加，位移不斷增大，原因是上述3個斷面埋深相對較大。

(2)支護結構Bi處的拱底位移受到圍巖蠕變荷載的影響，隨著時間的增加，Bi處的位移變化趨勢類似于Ai處的拱頂位移變化，同時Bi處的位移數值要小于Ai處。

(3)Ci與Di之間的收斂位移的變化趨勢也隨著斷面位置的變化而不同，其中C2與D2，C3與D3，C4與D4之間的收斂位移在蠕變荷載作用下，隨時間呈現減小的趨勢，C1與D1，C5與D5之間的收斂位移則隨著時間的增大不斷增大。

(4)通過分析同一數值采集斷面下支護結構位移隨時間變化曲線，發現拱頂處與拱底處的位移變化趨勢類似，左右拱腰處的收斂位移變化趨勢相反，位于隧道兩端的斷面1和5各位置的位移變化趨勢是相似的，而位于中部的數值采集斷面各位置處則與兩端位置有著相反的變化趨勢。