基于響應面分析的隧洞開挖過程優(yōu)化設計研究

南 軒， 劉艷慧， 孫語晨， 張玉波， 張明飛

(1.云南農(nóng)業(yè)大學 建筑工程學院， 云南 昆明 650000； 2.云南省高校城鄉(xiāng)水安全與 節(jié)水減排重點實驗室， 云南 昆明 650000)

1 研究背景

鑒于龍開口水電站一期工程中隧洞施工存在的側壁坍塌問題，目前國內(nèi)外隧洞的開挖主要研究在隧洞的爆破和開挖方式及支護方面，貊祖國等[1]研究不同開挖方式下的隧洞引水過程中溫度應力耦合；章誠等[2]研究了連續(xù)爆破下的圍巖應力的變化；于小朋等[3]研究隧洞群的支護方式，但對于隧洞開挖過程中圍巖受力的變化情況研究較少，因此本文基于ANSYS APDL隧道分步開挖數(shù)值模擬方法對水電站二期設計中與一期的7#和9#塌方隧洞巖體實際參數(shù)相似的1#隧洞，采用D-P屈服準則描述巖體的關聯(lián)流動程度通過測定單軸抗壓、抗剪及相關參數(shù)[4-10]，應用生死單元法通過APDL編寫單元命令實現(xiàn)分部開挖(第1步實現(xiàn)隧洞頂部開挖、第2步實現(xiàn)隧洞洞身上半部分的開挖、第3步實現(xiàn)洞身下半部分的開挖)了解其開挖過程應力和位移的變化情況，并通過參數(shù)設計試驗進行優(yōu)化分析[11]，最終確定1#隧洞最優(yōu)設計方案保證側壁的應力最小,減小塌方的可能性。

2 隧洞開挖數(shù)值模擬設計

2.1 D-P準則的應用

1#隧洞長78.60 m，設計流量1.43 m3/s，加大流量1.79 m3/s，穿越地層為古生界二疊系上統(tǒng)玄武巖組上段(P2β3)致密狀玄武巖夾凝灰?guī)r，隧洞最大埋深約30 m，進、出口段屬于V類圍巖，洞身段圍巖以IV類為主，圍巖局部穩(wěn)定性差。巖體描述在之前的研究中通常采用M-C準則但其收斂性較差，通過大量的修正和計算表明D-P準則能夠更加完整和準確地描述巖體在開挖過程中的破壞和變形情況。Drunker-Prager準則是在M-C準則和塑性mises準則的基礎上擴展得到，可表示為：

(1)

式中：I1為應力張量第一不變量；J2為應力偏量第二不變量，偏應力張量第二不變量J2的大小可用來判斷物體所處的彈塑性狀態(tài)；α、k均為與材料的黏聚力和內(nèi)摩擦角φ有關的參數(shù) 。

I1，J2可通過主應力表示為：

I1=σ1+σ2+σ3

(2)

(3)

式中：σ1、σ2、σ3分別為最小主應力、中間主應力、最大主應力，MPa。

與M-C準則相比，D-P準則不僅考慮了中間主應力對材料的屈服影響，通過實際巖體的剪切強度和抗壓強度以及內(nèi)摩擦角考慮了靜水壓力和中間主應力，能更加全面地反映隧洞開挖過程中巖體的變形和應力的變化。

根據(jù)D-P的屈服流動準則在巖體塑性變形階段其變形量滿足：

(4)

在整個彈塑性屈服的階段，其總的應變量為彈性形變和塑性形變的總和:

(5)

2.2 生死單元法的應用

在隧洞施工過程模擬中，需要對模型進行合理的設置，按照龍開口水電站二期1#隧洞的實際模型進行1∶1的建模，在施工過程中分為多步開挖，施工過程中每一次的開挖就相當于去掉這個開挖單元，在建模過程能去掉有限單元的方法有很多，但由于開挖過程中每次開挖都對后面的開挖過程有較大的影響，因此需要用生死單元法，生死單元法可以記錄前一施工步結構的變形狀態(tài)、應力、位移，并可將前一施工步中的變形和應力傳遞給下一步的施工。

2.3 有限元模型的建立

輸水隧洞的模型尺寸為1.5 m×1.8 m，按照龍開口隧洞處山體形狀實測和設計施工尺寸要求進行設置，在相應的邊界處進行加密保證網(wǎng)格質(zhì)量良好，網(wǎng)格質(zhì)量為0.84時網(wǎng)格數(shù)量為1 232 987，經(jīng)計算滿足網(wǎng)格無關性的要求(網(wǎng)格繼續(xù)加密對計算結果無影響)，保證數(shù)值模擬結果的準確性。有限元模型如圖1所示。

2.4 實驗數(shù)據(jù)確定和設計

由于ansys模擬過程中需對巖體材料屬性進行Drunker-Prager屈服準則設計，因此需要測定出龍開口1#隧洞的抗壓強度與變形量的關系和抗剪強度。基于實驗和相關文獻得到如表1的龍開口水電站1#隧洞巖體的相關參數(shù)數(shù)據(jù)[13]。

上述參數(shù)完全能夠滿足Drunker-Prager屈服準則設計下數(shù)值分析過程中巖體的變化情況。

3 數(shù)值模擬結果分析

數(shù)值模擬中采用生死單元法，這有助于詳細了解施工過程中應力、應變、變形等的變化情況，首先分析開挖過程中的巖體變形情況，開挖過程中巖體變形云圖如圖2所示。

圖1 龍開口水電站1#隧洞巖體有限元模型

表1 龍開口水電站1#隧洞巖體特性參數(shù)

圖2 龍開口水電站1#隧洞開挖過程中巖體變形情況

由圖2可以看出，在未開挖階段，整個巖體在重力加速度作用下產(chǎn)生變形，隨著上部巖體的開挖，隧洞周圍的巖體出現(xiàn)了變形不均勻的情況，隧洞頂部的巖體變形量遠大于底部的變形量，而且隧洞底部與頂部的變形在方向上也是不同的。頂部的變形主要是由重力引起的向下的位移，當巖體較松軟時頂部塌陷的機率較大，底部的巖體受到兩側的圍巖擠壓之后形成向上的位移，同時可以發(fā)現(xiàn)在側邊處的變形量較大，側邊處易于塌方，這與實際的施工中塌方主要是側壁處應力較大并且龍開口輸水隧洞巖體較松軟的情況基本一致，與相關的工程實際和相關文獻相吻合[12-13]。

由圖2可以看出，在隧洞頂部處變形出現(xiàn)峰值，監(jiān)測隧洞頂部的變形量能更加詳細地了解開挖的進程。隨著開挖步的進行，隧洞頂部的變形量逐漸增大，從未開挖時最大變形量5.894×10-5m到第1步開挖之后的最大變形量6.502×10-5m，再到第2步開挖時為1.071×10-4m，最后一步開挖后為1.098×10-4m。可以通過隧洞沿程的變形曲線了解到每一步開挖過程的變化情況，如圖3所示。

由圖3可以看出，開挖過程中第1步開挖之后與未開挖的頂部變形量相差不大，這主要是由于弧頂?shù)拈_挖量相對較小，但是在第2步開挖之后的頂部變形量遠大于第1步開挖之后，同時第3步開挖之后的頂部變形量與第2步開挖之后頂部的變形量相差不大。

隧洞的穩(wěn)定性與其隧洞圍巖處所受到的應力有密切的關系，通過應力云圖可以了解其應力的分布情況，并且得到其側邊處的應力變化情況，應力云圖如圖4所示。

圖3 開挖過程中隧洞頂部沿程變形曲線

圖4 隧洞開挖過程中應力變化云圖

由圖4可以看出，隨著開挖進程的深入，巖體逐漸出現(xiàn)應力集中的現(xiàn)象，其應力集中在隧洞的附近的圍巖之中，在隧洞的側壁面中間位置處的應力出現(xiàn)最大值，同時根據(jù)沿程應力曲線和沿程應變曲線(圖5)可以發(fā)現(xiàn)，在隧洞第2步開挖完成后，整個巖體應力峰值和應變峰值變化較大，應力峰值從8.72×105Pa增加到1.47×106Pa，應變峰值從2.955×10-5m增大到4.906 ×10-5m，結合整個過程中變形情況可以確定，在整體的開挖過程中第2步開挖對應力、應變的影響最大。

圖6為整個開挖過程中隧洞側壁的應力和變化矢量變化。由圖6可清晰地了解應力和變形參量在隧洞每一步開挖完成后的變化情況，未開挖時，隧洞壁面矢量方向沿隧洞反向方向偏下，在第1步開挖之后，其矢量方向基本沒有太大的變化，隨著第2步開挖的進行隧洞周圍的巖體的應力、變形的矢量方向發(fā)生了顯著變化，其方向已經(jīng)出現(xiàn)偏轉，隨著開挖繼續(xù)進行其矢量的方向與側壁面夾角增大，與側壁面在實際中的坍塌方向一致，在整個開挖過程中隧洞圍巖處的應力、變形的矢量方向在不斷變化最終與側壁形成夾角并指向隧洞內(nèi)部，若隧洞施工過程中沒有及時支撐，或者巖體較軟時易于形成側邊巖體塌方。

通過模擬與分析隧洞開挖過程中圍巖的形變、應力、應變等的變化，可以了解到在開挖過程中隧洞的塌方主要是由于側壁處的應力、應變、變形等增大并且其矢量方向從側壁面向外逐漸轉變?yōu)閺谋诿嫦騼?nèi)，如何在施工中減小或改善這種狀況，尤其是在開挖的巖體所對的Drunker-Prager屈服準則以及非關聯(lián)流動與各向同性硬化確定的情況下，如何減小其坍塌尤為重要。目前可靠性分析開始被大量運用在結構和工程的失效、疲勞、破壞等方面[14-18]，通過前人在工程可靠性方面的研究分析，本文將進一步運用響應面分析和優(yōu)化對龍開口水電站1#輸水隧洞側壁應力和變形提出改善意見。

圖5 隧洞開挖過程中應力及應變沿程變化曲線

圖6 整個開挖過程中隧洞側壁應力和變形的矢量變化

4 可靠性及優(yōu)化設計分析

輸水隧洞的可靠性分析和優(yōu)化設計可以保證隧洞耐用程度和降低坍塌機率。隨著計算機仿真技術的發(fā)展，試驗設計的范圍和結果更加準確，可通過評價試驗因素、確定最佳試驗參數(shù)組合、分析輸入?yún)?shù)和輸出參數(shù)之間的關系、構建回歸方程和近似模型等方法提高工程設計的穩(wěn)定性，本試驗選用中心組合設計(central composite design)由2k全因子設計、軸點設計與零水平的中心點重復試驗3部分構成，應用各影響因子的二次多項式來預測其對評價指標的作用。該方法擴展了試驗設計空間并得到高階信息，能夠給響應面近似模型提供樣本數(shù)據(jù)，具有設計簡單、試驗次數(shù)少、預測性好等優(yōu)點，在優(yōu)化設計中由于山體中的巖體和土體為確定的參數(shù)，因此只能對隧洞的埋深、斷面的長寬尺寸、隧洞的長度進行相應響應面優(yōu)化和修改。

響應面分析法，即響應曲面設計方法(Response Surface Methodology，RSM)，是利用合理的試驗設計方法并通過試驗得到一定數(shù)據(jù)，采用多元二次回歸方程來擬合因素與響應值之間的函數(shù)關系，通過對回歸方程的分析來尋求最優(yōu)工藝參數(shù)，解決多變量問題的一種統(tǒng)計方法，響應變量y與各輸入變量xi的關系可用簡單的代數(shù)形式表達出來：

(6)

當通過中心組合試驗獲得M個樣本點對應的響應量y=(y(1),y(2),…,y(M))后，利用最小二乘法得到系數(shù)矩陣：

β=(θTθ)-1(θTy)

(7)

式中：θ為響應面樣本矢量點，可以表示為：

(8)

結合公式(6)～(8)可得二階響應面近似模型的表達式。

以龍開口水電站1#隧洞尺寸設計為優(yōu)化對象，將ANSYS數(shù)值模擬計算和響應面方法相結合，對側壁應力和變形進行優(yōu)化分析。

4.1 輸水隧洞設計參數(shù)敏感性評價

輸水隧洞設計主要參數(shù)包括隧洞的埋深、斷面的長寬尺寸、隧洞的長度等，要了解這些參數(shù)能對應力、變形產(chǎn)生多大的影響，靈敏性分析方法(定量分析某一組研究數(shù)據(jù)對某些因數(shù)響應值變化的分析技術，實質(zhì)是逐步改變相關變量值以了解關鍵指標受這些因數(shù)影響的規(guī)律)能準確分析各因素的權重，是確定因數(shù)對結果影響的重要評價方法。在多目標優(yōu)化(MOGA)方法中很難求出解析解，這里借鑒故障樹分析法中計算基本事件臨界重要度的方法來進行參數(shù)的敏感性評價分析，用S表示敏感度[19-21]，定義為：

(9)

式中：Δyi為應力峰值的變化量，Pa;y0為設計模型下的隧洞側壁應力峰值,Pa;Δxi為試驗中某因數(shù)改變量;x0為設計模型對應因數(shù)參數(shù)。S越大表明應力對該因素的敏感度越高，因此借助敏感度可以了解各因素對應力的相對影響程度，以選出最大的不確定因數(shù)。

表2為通過計算得出的隧洞幾何參數(shù)敏感度：

從表2各因素對應力峰值的敏感度可以看出，隧洞的長度對側壁的應力峰值無影響，隧洞埋深對應力峰值有較大的影響，隧洞寬度對應力峰值的影響為負相關，因此隧洞長度為無關變量，在設計時不考慮隧洞長度的影響。

4.2 中心復合設計方法的應用

在響應面的試驗中應用Design Expert軟件采用中心復合試驗設計方法，對隧洞幾何尺寸因素進行試驗設計，見表3。試驗總次數(shù)為20次,其中14個為外部分析點，6個區(qū)域中心點，能夠在重復計算中降低誤差[22-26]。

表3 中心復合設計試驗各因素編碼水平

4.3 響應及優(yōu)化設計結果

中心復合設計試驗結果見表4。

將試驗結果的數(shù)據(jù)代入二階響應面試驗近似模型擬合后得到應力峰值基于實際因子水平的回歸方程,可表示為：

(10)

(11)

相關系數(shù)為完全擬合的度量值，反映了響應面與試驗數(shù)據(jù)的相符合程度，要在0.9以上才能保證擬合效果。

(12)

表4 中心復合設計試驗結果

4.4 影響應力峰值各參數(shù)間交互作用結果分析

試驗將隧道側壁應力作為研究隧洞坍塌問題的主要考察指標，影響側壁應力峰值的各參數(shù)之間交互作用下的側壁應力峰值變化趨勢見圖7。參數(shù)包括隧洞截面長度H1、截面寬度H2和埋深V1。

從圖7中的各因素面的斜率情況可以清楚的了解到，對于隧洞側邊的應力峰值大小，從影響效果上來看是隧洞埋深最大、隧洞截面寬度次之、隧洞截面長度影響較小；從響應面分析結果可以看出，所有的參數(shù)與應力峰值之間均處于線性，僅僅是響應程度不同，這也表明在隧洞的設計過程中，在設計許可的范圍內(nèi)盡可能減少截面長度、增大截面寬度，從而保證最終形成的側壁應力最小并且降低龍開口水電站1#隧洞松軟巖體的坍塌機率，使輸水隧洞保持良好的穩(wěn)定性。

為了綜合考慮實際施工過程中的參數(shù)設計，在優(yōu)化過程中沿對角線選取外部分析點，由于隧洞的埋深與山體的位置以及水庫的位置相關且不易改變，因而在不改變埋深的情況下選取：隧洞截面長度H1=1.35 m和隧洞截面寬度H2=2m、隧洞埋深V1=30 m、隧洞長度L1=78.6 m。

在ANSYS中按照優(yōu)化前的開挖步驟和方法求解并與優(yōu)化后的應力與變形曲線進行對比，可以發(fā)現(xiàn)整個隧洞側壁的沿程應力峰值在初始設計尺寸下為1.4159×106Pa，優(yōu)化尺寸設計后的側壁沿程應力峰值為1.2741×106Pa，側壁沿程應力峰值下降了10.01%；同時對沿程應變峰值進行對比,初始設計的應變峰值為4.7053×105，優(yōu)化后的沿程應變峰值為4.2313×105，側壁沿程應變峰值下降了10.07%，整體的變化情況與之基本相同。具體如圖8所示。

表5 應力峰值回歸方程的方差分析

圖7 隧洞各幾何參數(shù)之間交互作用下的應力峰值變化趨勢

圖8 中心復合試驗優(yōu)化前后應力及應變沿程變化曲線

5 結 論

在分析隧洞開挖過程中隧洞內(nèi)部的巖體受力和變形時，采用生死單元和數(shù)值分析相結合，得到內(nèi)部的巖體受力情況；同時將敏感性評價分析和響應面法應用于隧洞內(nèi)部壁面處沿程受力變化情況的優(yōu)化，選取對受力情況有影響的幾何參數(shù)，分析各影響隧洞沿程受力參數(shù)間交互作用的結果，確定符合隧洞設計要求的最優(yōu)化參數(shù)組合。得到如下結論：

(1)在生死單元法開挖的過程中，從第1步開挖完成到第2步開挖完成，這個階段隧洞側壁巖體受到的應力在大小和方向上均發(fā)生了極大的變化，應力值從9.7×105Pa增加到1.467×106Pa，其方向從沿側壁向外轉變?yōu)檠貍缺谙騼?nèi)，此變化極大地增加了隧洞坍塌的可能性。

(2)應用中心復合設計的方法對選取的3個敏感度較高的參數(shù)進行試驗設計，建立了隧洞沿程應力應變的多元回歸模型，同時采用方差分析驗證了模型的擬合度。結果表明：回歸模型能夠準確地反映影響參數(shù)與響應值之間的變化關系，并且對實際情況擬合良好。

(3)通過建立的多元回歸模型求解，得到了截面尺寸在隧洞設計要求內(nèi)的最優(yōu)組合設計參數(shù)，在ANSYS中驗證之后，發(fā)現(xiàn)比原始模型設計得到的沿程應力峰值降低10.01%，應變降低10.07%。