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

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

0 引 言

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

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

1 工程地質(zhì)條件

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

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

2 圍巖變形監(jiān)測結(jié)果

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

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

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

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

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

3 砂巖本構(gòu)模型及力學參數(shù)

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

(1)

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

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

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

4 圍巖力學參數(shù)反演

4.1 遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡模型

通過MATLAB對遺傳算法優(yōu)化的BP神經(jīng)網(wǎng)絡模型進行設計和編程，網(wǎng)絡模型的建立主要分為以下兩個部分。

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

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

(3)

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

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

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

4.2 訓練神經(jīng)網(wǎng)絡樣本

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

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

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

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

4.3 力學參數(shù)反演及合理性驗證

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

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

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

5 隧洞運行期穩(wěn)定性預測

5.1 建立計算模型

利用參數(shù)化程序設計語言編寫B(tài)urgers 流變方程并接入通用有限元軟件ANSYS，程序每次迭代計算時均調(diào)用一次APDL，并更新應力和應變，從而實現(xiàn)砂巖流變模型的導入。

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

圖5 模型有限元網(wǎng)格Fig.5 Model finite element mesh

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

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

5.2 計算結(jié)果分析

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

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

圖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

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

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

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

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

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

6 結(jié) 語

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

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

□