瓦屋山引水隧洞圍巖流變特性數值分析

王兆強，陳 新，張 磊，朱文凱

(四川大學水利水電學院，成都 610065)

0 引 言

在水利工程建設中，常常通過引水隧洞實現水資源跨流域、跨區域建設。引水隧洞開挖后圍巖變形受開挖空間和時間效應的影響，呈強非線性變化，建立不同類別圍巖的時效變形預測模型對工程設計和維護具有重要意義。變形預測的方法較多，但針對引水隧洞圍巖運行期圍巖變形預測的研究仍較少。金長宇等[1]建立了ANFIS模型來預測龍灘水電站地下廠房施工期的變形；張倩等[2]采用遺傳算法優化BP神經網絡，結合動態分析法建立了施工期圍巖變形預測模型。目前傳統預測模型注重對監測數據本身規律性的擬合，忽視了圍巖變形的力學機制，此類預測模型會隨著圍巖應力狀態的改變而失效，常常用于短期預測，為臨時支護結構的設計提供依據。

鑒此，本文針對砂巖的流變特性，采用Burgers模型反映其流變變形規律，基于現場位移監測資料，通過遺傳算法優化的BP神經網絡對模型中的力學參數進行反演，利用參數化程序設計語言編寫Burgers 流變方程并接入通用有限元軟件ANSYS模擬隧洞巖體充水運行條件下的流變力學行為,取得了良好的應用效果。

1 工程地質條件

瓦屋山水電站樞紐位于洪雅縣炳靈鄉，樞紐由攔河大壩、左右岸泄洪隧洞、引水系統和廠房等建筑物組成。其中引水系統包括引水隧洞、調壓井和埋藏式壓力管道。引水隧洞主洞長4 765.322 m，斷面為圓形，內徑6.2 m，襯砌厚度40～80 cm。據地勘資料，洞周圍巖類型包括：Ⅲ類圍巖，長3 260 m，樁號0+239.97～3+499.97，以巨厚層～中厚層砂巖為主，部分薄層砂巖泥質巖類；Ⅳ類圍巖，長1 425 m，樁號0+079.97～0+239.97、3+499.97～4+439.97、4+499.97～4+825.292，以中厚層泥質粉砂巖或中厚層粉砂質泥巖夾薄層粉砂質泥巖或泥質粉砂巖為主；Ⅴ類圍巖，長140 m，樁號0+000～0+079.97；4+439.97～4+499.97，風化及卸荷的中厚層泥質粉砂巖夾薄層粉砂質泥巖。Ⅲ、Ⅳ類圍巖中尚局部存在Ⅴ類圍巖，主要巖體力學參數見表1。

表1 引水隧洞圍巖分類及巖體力學參數建議值Tab.1 Classification of surrounding rock of diversion tunnel and suggested values of rock mass mechanical parameters

2 圍巖變形監測結果

圖1給出了隧洞收斂監測斷面測點、測線的布置。結合工程地質和施工條件，確定施工監測方案以洞周收斂和錨桿應力監測為主。其中變形監測斷面按照分段控制的原則，在Ⅲ、Ⅳ類圍巖洞段選擇6個變形監測斷面，在開挖近掌子面的設計斷面上埋設測點，測點按照圓形斷面5測點8測線進行布置。本文根據固定點任意三角形法計算得出引水隧洞圍巖測點的位移。固定點任意三角形法假設：周公河側測點5處于坐標原點且不發生變形；測點1和測點5處于相同水平位置，測點1只作水平方向變形；測點2、測點3和測點4的變形分解為水平和豎直兩個方向位移。

圖1 隧洞收斂監測斷面測點、測線布置圖Fig.1 Tunnel convergence monitoring section measurement point, line layout

本文選取1號支洞下游面的第二斷面(K1+90)作為研究對象，原因包括：該斷面巖石為Ⅲ類圍巖，在三類主要的圍巖中占比到達69%，具有一定的代表性；該斷面監測設備埋深時間為2005年8月16日，監測時間長達2個月，監測數據充足；期間斷面開挖活動停止后只進行臨時支護，期間圍巖應力狀態變化不大，易于模擬分析；開挖活動停止后各測線收斂值仍在增加，砂巖的流變特性得到充分表現，有利于對圍巖流變參數的反演分析。

圖2給出了1號支洞下游面的第二斷面5個測點收斂值，基于固定點任意三角形法假設測點5處于坐標原點且不發生變形故不在圖中列出。分析位移監測值可知，監測斷面整體向周公河側變形，主要由于靠山側圍巖的側向壓力和水壓力較大。

圖2 測點位移監測值(2005年)Fig.2 Measuring point displacement monitoring value

3 砂巖本構模型及力學參數

王志儉等[3]通過對砂巖巖石和巖體進行流變試驗和分析，認為Burgers 流變模型可以準確描述砂巖的流變特性。Burgers 流變模型為四元件兩單元模型，在三軸應力條件下，該模型的軸線應變表達式為：

(1)

式中：G2為黏彈性剪切模量；η1為黏性系數；η2為黏彈性系數；K為體積模量；G1為剪切模量；E1為彈性模量；η為泊松比。

其中η1、η2、E1和G2四個力學參數需要通過參數反演確定。根據地質資料、三軸試驗成果以及工程經驗綜合確定確定出了四個參數取值范圍，見表2。

表2 反演參數取值范圍Tab.2 Inversion parameter value range

4 圍巖力學參數反演

4.1 遺傳算法優化的BP神經網絡模型

通過MATLAB對遺傳算法優化的BP神經網絡模型進行設計和編程，網絡模型的建立主要分為以下兩個部分。

(1)遺傳算法優化BP神經網絡。 基于BP神經網絡，采用遺傳算法進行全局搜索，得到BP神經網絡最佳閾值和初始權值，再將優化后的閾值、權值輸入BP神經網絡，神經網絡的局部搜索功能可以快速獲得近似最優值。遺傳算法的優化解決了BP神經網絡易陷入局部極小值、穩定性差與網絡收斂慢的問題。

(2)擬定神經網絡結構。BP神經網絡的結構包括的輸入層、輸出層和隱含層，其中輸入層與輸出層節點數根據實際問題確定，而隱含層節點數N可以通過公式(3)確定。通常先擬定一個初始值用于BP神經網絡訓練，若訓練的結果不滿足要求，則增加一個隱含節點直至滿足要求。

(3)

式中：a為1～10之間的常數；m、n分別為輸入層節點數和輸出節點數。

通常先擬定一個初始值用于BP神經網絡訓練，若訓練的結果不滿足要求，則增加一個隱含節點直至滿足要求。

根據BP神經網絡結構的布置要求，以待反演力學參數作為神經網絡的目標向量，輸出節點數為4；以位移監測收斂值作為BP神經網絡的輸入向量，輸入節點數為280。根據公式(3)隱層節點數初始值擬定為20層。

4.2 訓練神經網絡樣本

訓練神經網絡樣本是一個不斷尋找最佳閥值與權值的過程，即通過訓練使得網絡的輸出誤差越來越小。本節擬利用Levenberg-Marquardt算法對神經網絡樣本進行訓練。

基于均勻設計法[4]，選用U11(117)均勻設計方案對前文四個力學參數η1、η2、E1和G2進行試驗設計組合，得到11組不同組合的參數。通過數值分析模擬計算11組力學參數分別對應的測點位移值，具體模擬過程可以參考第六節。

圖3 神經網絡訓練樣本(2005年)Fig.3 Neural network training samples

11組力學參數和對應的測點位移值共同組成人工神經網絡的訓練樣本，為力學參數反演提供初始樣本，由于篇幅有限，僅列出η1=1.03×1018、η2=1.21×1016、E1=39和G2=9.80×1010組合時的訓練樣本以供參考，見圖3。采用初次訓練好的網絡樣本進行參數反演，可能由于樣本數量偏少而不能滿足工程要求，可以繼續增加網絡訓練樣本，直至得出合理的反演參數。

4.3 力學參數反演及合理性驗證

通過MATLAB對遺傳算法優化的BP神經網絡法進行設計和編程，將K1+90監測斷面5個測點70 d監測值作為輸入向量導入初始訓練樣本中，輸出四個力學參數，結果為：η1=7.13×1017、η2=1.27×1016、E1=39、G2=1.80×1011。

將反演參數導入計算模型中，通過數值模擬得到K1+90斷面各測點70 d的流變位移，圖4給出了位移計算值與監測值得對比結果，各測點計算結果與監測數據擬合度較好，相對誤差基本控制在10%以內，符合工程需要。因此Burgers流變模型和力學參數可以較好的描述該水電站砂巖地層引水隧洞1號支洞下游面圍巖流變特性，可以用于引水隧洞圍巖運行期變形預測。

圖4 計算位移值和監測值對比分析(2005年)Fig.4 Comparison of calculated displacement values and monitored values

5 隧洞運行期穩定性預測

5.1 建立計算模型

利用參數化程序設計語言編寫Burgers 流變方程并接入通用有限元軟件ANSYS，程序每次迭代計算時均調用一次APDL，并更新應力和應變，從而實現砂巖流變模型的導入。

由于引水隧洞沿軸線方向較大范圍內圍巖巖性單一，因此，該問題可簡化為平面應變問題。為減少邊界效應，網格模型上下邊界的設定均大于洞徑的10倍，已知引水隧洞1號支洞下游面斷面為圓形，內徑6.2 m，襯砌厚度為0.6 m，因此模型計算范圍為80 m×80 m，圍巖和襯砌模擬單元均采用平面單元Plane42。模型有限元網格如圖5所示。

圖5 模型有限元網格Fig.5 Model finite element mesh

模型邊界條件及初始條件：左右兩側、底部均采用法向約束；初始應力場以自重應力場為主；外水壓力對引水隧洞的圍巖和襯砌結構的長期穩定性有較大的影響，因此應當考慮外水壓力對隧洞的影響，根據工程地質報告，對隧洞模型頂部施加80 m的壓力水頭，底面施加160 m的壓力水頭，側面施加沿重力方向梯度變化的壓力水頭。

給定初始應力場和外水壓力的模擬完成后，依次進行隧洞全斷面開挖模擬及隧洞支護模擬(初期采用噴混凝土、掛網、錨桿支護)。由于開挖與襯砌施作之間的間隔時間較長，圍巖變形量較大，開挖與襯砌施作期間及竣工之后還需要進行蠕變模擬。經過上述施工模擬建立運行期的初始應力場，在此基礎上施加運行期的內水壓力，最后分別對隧洞充水運行初期(流變時間為T=0 d)、T=10 d、T=30 d、T=60 d、T=120 d、T=240 d，T=360 d不同時段的流變情況進行模擬，每次流變模擬均調用Burgers 流變模型，并更新應力和應變，進而對引水隧洞長期的穩定性進行評價。建立圍巖流變模型所需的力學參數η1、η2、E1和G2通過上述參數反演已經確定，其余計算所需的力學參數見表3。

5.2 計算結果分析

圖6、圖7給出了隧洞充水運行期流變時間T=60 d和T=360 d的位移場云圖。隧洞運行T=60 d后，從位移的Y分量和位移的X分量云圖可以看出：隧洞頂拱變形方向向下，Y方向最大位移分量為3.316 1 mm；底拱變形方向向上，Y方向最大位移分量為2.131 7 mm；山體側隧洞水平位移分量要大于周公河側，均向左側發展變形，X方向最大位移分量為5.661 8 mm。

隧洞運行一年T=360 d，與隧洞運行60 d的位移情況相比發現：隧洞頂拱繼續向下變形，Y方向最大位移分量達到3.836 1 mm；底拱發生少量回彈變形，變形方向向下，Y方向最大位移分量減小為2.001 2 mm，位移方向向上；山體側隧洞水平位移分量仍大于周公河側，繼續向左側發展變形，最大位移分量為達到5.661 8 mm；位移增量較小，流變變形處于收斂狀態。

圖6 T=60 d圍巖及襯砌位移云圖Fig.6 T=60 d displacement cloud chart of surrounding rock and lining

圖7 T=360 d 圍巖及襯砌位移云圖Fig.7 T=360 d displacement cloud chart of surrounding rock and lining

為了系統分析引水隧洞1號支洞下游面圍巖的流變規律，對不同時刻的位移云圖進行整理，得到對應時刻圍巖特征點的流變位移，見表4。從量值來分析，圍巖總體變形量較小，在流變計算360 d后，圍巖特征點位移約為4.5～6.9 mm。從流變趨勢分析，隧洞圍巖各特征點流變趨勢都較為一致，初期變形較大，前10 d的變形量達到2.6～4.2 mm，之后變形速率逐漸減小。流變在10～30 d進入流變衰減期，在第30 d之后變形趨于穩定，進入穩定蠕變階段，平均變形速率為0.005～0.006 mm/d。

分析結果表明：該斷面圍巖變形能在30 d內趨于穩定，穩定變形速率很低，且隨著流變時間的推移洞周未出現新的塑形區，因此，運行期圍巖流變情況滿足引水隧洞圍巖長期穩定性的要求。

表4 不同時刻隧洞圍巖特征點位移Tab.4 The characteristic point displacement of surrounding rock at different moments

表5給出了不同流變時刻洞周圍巖應力最值的統計結果。從量值上分析，運行期圍巖和襯砌結構主要出于受壓狀態，拉應力與流變時間成反比，壓應力與流變時間成正比。各時刻大應力的最大壓應力均出現在洞底右側及洞頂左側區域，洞周圍巖未出現拉應力區。各時刻小主應力最大壓應力均出現在圍巖兩側3～4 m區域內，拉應力主要出現在襯砌結構底板右側上，應力范圍隨與流變時間的推移而減小。因此，運行期圍巖應力情況滿足引水隧洞圍巖長期穩定性的要求。

表5 流變計算過程中應力最值Tab.5 The maximum stress in the rheological calculation

6 結 語

(1)本文針對砂巖的流變特性，采用Burgers模型反映其流變變形規律。利用現有的監測資料，通過遺傳算法優化BP神經網法反演得到模型中的力學參數。利用APDL編寫Burgers 流變方程并接入ANSYS軟件以模擬隧洞運行期穩定性分析。該預測模型能夠適應隧洞充水運行條件下圍巖應力狀態的改變，研究成果可為引水隧洞運行期的穩定性分析提供參考。

(2)計算結果表明，該洞段引水隧洞整體穩定性良好。位移值變化較小，變形收斂速度快，穩定變形速率低;應力值趨于均勻化，應力集中程度漸少，塑性區面積較小，且進入流變穩定期后未出現新增塑性區。支護的受力情況也是安全可靠的。

□