楊文鈺，鄭俊杰，章榮軍，喬雅晴

(1.華中科技大學 土木與水利工程學院，武漢 430074；2.武漢大學 土木建筑工程學院，武漢 430072)

隨著城市的發(fā)展，盾構隧道成為充分利用地下空間最為有效的手段之一。掌子面穩(wěn)定性是盾構隧道的經典問題，合適的掌子面支護壓力不僅能夠保證盾構施工的安全，同時也保障了盾構隧道周圍建構筑物的安全。

對于盾構隧道開挖面失穩(wěn)問題，大量的學者運用理論分析法[1-3]、數值模擬法[4-6]、模型試驗法[7-9]已經得到了豐富的研究成果。然而，土體的天然變異性與試驗或場地的限制導致了人們對巖土體的認知缺乏，這決定了巖土體參數的不確定性[10]。以上這些研究都把土體視作均質且各向同性的材料，在預測極限支護壓力時可能會存在偏差。近年來，有學者將土性參數的不確定性考慮在盾構隧道掌子面穩(wěn)定性的研究中，Mollon等[11-14]基于隨機響應面等方法，研究了抗剪強度參數的變異性對掌子面穩(wěn)定性的影響，揭示了土性參數變異性對極限支護壓力的影響的基本規(guī)律。以上研究將土性參數視作隨機變量，為進一步考慮土性參數的空間變異性，Mollon等[15]首先基于極限分析理論，得到了考慮抗剪強度參數空間變異性二維掌子面失效機制，Cheng等[16-19]考慮了砂土、黏性土以及成層土抗剪強度參數的空間變異性，研究表明，抗剪強度空間變異性對掌子面穩(wěn)定性有重要影響，忽略這一特性可能高估掌子面的穩(wěn)定性；極限支護壓力與抗剪強度參數的變異系數、自相關距離密切相關。以上研究較全面地揭示了土性參數的不確定性對盾構隧道掌子面穩(wěn)定性的影響規(guī)律，然而，鮮少有研究考慮掘進參數的不確定性。

在盾構隧道施工過程中，掌子面土體經過刀盤的切削作用進入土艙內部，掌子面處土體的穩(wěn)定靠土艙壓力與掌子面處水土壓力相平衡來保持，盾構機推進油缸的推力、推進速度、螺旋輸送機的出土量[20]與刀盤開口率[21]等盾構機的參數都對土艙壓力有一定的影響。通常情況下，為保持盾構掌子面的穩(wěn)定，先根據地質條件設定土壓力值，在盾構掘進過程中再根據土壓傳感器的變化人為做出調整[22-23]。但目前土艙壓力值的計算尚無固定方法[24]，一般參考工程地質、盾構機型等因素，由施工經驗確定。工程地質條件的復雜性決定了土艙壓力的值并非定值[25]；另外，盾構機的人為控制可能存在操作不當、違規(guī)操作等問題，也會導致盾構施工中土艙壓力具有一定的變異性。因此，忽略支護壓力(本文中提及的支護壓力僅限于盾構中的土艙壓力)變異性可能無法反映工程實際的復雜性，相應地可能會高估掌子面穩(wěn)定性。基于此，筆者采用隨機場理論與數值模擬分析相結合的方法，研究了黏性土內摩擦角與黏聚力空間變異性共同作用下對盾構隧道掌子面穩(wěn)定性的影響，分析了黏聚力與內摩擦角變異性共同作用下對掌子面失穩(wěn)模式和極限支護壓力的變化規(guī)律，并結合概率分析法探討了同時考慮土性參數空間變異性與支護壓力變異性的可靠度分析方法。

1 數值模型的建立

1.1 問題描述

建立數值模型的核心問題是盾構隧道在黏性土掘進過程中的掌子面穩(wěn)定性。盾構隧道掌子面穩(wěn)定性的數值模型參考了Mollon等[12]的研究，數值模型如圖1所示。采用的數值模型隧道直徑(襯砌管片外徑)D為10 m，埋深C(隧道拱頂距地表距離)為10 m，為簡化分析，不考慮地下水的影響。整個模型的計算域為50 m×40 m×26 m(長×寬×高)，經驗證，計算域的大小能使計算結果不受邊界的影響[12]。數值模型共有52 240個單元。由于采用應力控制法分析盾構隧道掌子面穩(wěn)定性時，掌子面上會出現較高的應力梯度，所以對掌子面以及掌子面后方土體單元進行了網格的加密，掌子面上共分為了198個單元。模型的底部為固定邊界，四周為法向位移約束邊界，頂部為自由邊界。為了簡化分析，數值模型僅模擬了黏性土與襯砌管片兩種材料，本構模型采用Mohr-Coulomb模型，襯砌管片采用線彈性模型，襯砌管片的厚度為0.4 m。黏性土與襯砌管片材料參數如表1所示。

圖1 數值模型示意圖Fig.1 Schematic diagram for the numerical modelling of the problem in this

表1 黏土與襯砌材料參數Table 1 Properties for the clay and the lining segment

1.2 改進二分法

由于掌子面上的位移事先難以確定，所以采用較為常用的應力控制法進行計算。考慮到不確定性計算中的Monte-Carlo策略需要多次計算，選擇了計算效率較高的改進二分法。

改進二分法基于簡單二分法，簡單二分法的流程可簡單分為：1)確定計算上下限；2)將上下限的中間值代入計算，并根據計算結果更替上下限；3)不斷重復前兩個步驟直到計算結果達到特定的精度值。簡單二分法的本質是進行幾個穩(wěn)定與不穩(wěn)定的循環(huán)計算。然而，塑性流動往往出現在大量的計算步之后，尤其是對結果的精度要求苛刻時，將會耗費大量的時間。針對這一不足，Mollon等[12]提出了改進二分法，改進二分法的步驟如下：首先將土體的黏聚力賦予一個很大的值，這讓土體變成彈性材料；接下來，手動設置內部應力為初始值的兩倍，統(tǒng)計系統(tǒng)重新回到平衡狀態(tài)所需要的計算步數N，根據Mollon等[12]的研究，N的值接近于3 000步；N的值確定后，將初始黏聚力重置為真實值，若循環(huán)N步后系統(tǒng)仍不平衡，則可認為該工況會進入塑性流動狀態(tài)。由此可見，相比簡單二分法，改進二分法大大的縮短了塑性流動狀態(tài)工況的計算時間。所以，筆者采用改進二分法分析盾構隧道掌子面穩(wěn)定性。

1.3 數值模擬的驗證

在二分法中，選取合適的上下限對計算準確的極限支護壓力非常重要。選擇60 kPa與20 kPa作為改進二分法的上下限。圖2展示了當計算步數N為2 700時，掌子面中心點A(A點的幾何位置見圖1)位移與速度隨計算步數的變化。支護壓力為60 kPa時，在2 000步左右時，A點的水平位移已經趨于定值，A點的水平速度降為0；而當支護壓力為20 kPa時，計算2 700步后，A點的水平位移與水平速度仍在增長。說明對于表1參數的黏性土而言，分析掌子面穩(wěn)定性的改進二分法中N取2 700步，上下限取60 kPa與20 kPa是合理的。

圖2 掌子面A點水平位移與水平速度隨計算步數的變化Fig.2 Variation of the horizontal displacement and the horizontal velocity of point A on the tunnel face against the calculation

為驗證數值模型的合理性，首先運用了逐級減小支護壓力的應力控制法。支護壓力從100 kPa開始逐步減小，為使系統(tǒng)達到塑性流動的狀態(tài)，在每一支護壓力所對應的工況下計算10 000步。圖3是掌子面中心點水平位移-支護壓力變化曲線，如圖所示，曲線的拐點在34 kPa處，即極限支護壓力為34 kPa。

圖3 掌子面A點水平位移支護壓力變化曲線Fig.3 Variation of the horizontal displacement and of point A on the tunnel face against support

運用改進二分法進行計算，設置精度為1 kPa，計算結果為34.53 kPa，此時相應的速度場如圖4所示。這個結果不僅與應力控制法的計算結果相同，也與理論解[26]與數值計算結果[12]相符，這證明了模型與改進二分法的正確性。

圖4 支護壓力為34.53 kPa時的速度場Fig.4 Velocity field when the support

2 不確定性分析

2.1 基于隨機場理論的盾構隧道掌子面穩(wěn)定性分析

運用隨機場理論對盾構掌子面的穩(wěn)定性進行不確定性分析。采用隨機場分布光滑度與連續(xù)性較好的高斯型自相關函數[27]，能夠有效的描述土性參數的空間自相關性。高斯型自相關函數主要的形式如式(1)所示。

(1)

式中：τx、τy、τz分別為空間兩點在x、y、z方向上的相對距離；δx、δy、δz分別為x、y、z方向的波動范圍。在本文中，各向異性隨機場的水平向波動范圍δx、δy為20 m，豎向波動范圍δz為2 m。中心極限定理表明[28]，受大量不確定性因素影響的因變量通常近似服從正態(tài)分布或對數正態(tài)分布，同時對數正態(tài)分布嚴格非負[29-30]，這與巖土體參數的概率分布相符合，因此，采用對數正態(tài)分布描述黏聚力與內摩擦角的不確定性。

K-L級數分解法具有運算效率較高、生成隨機場精度較好、對于各向異性隨機場的生成有較強的適應性等優(yōu)點，被廣泛應用。K-L級數分解法將土體參數隨機場H(A,θ)(其中連續(xù)坐標A∈Ω?Rn，目標空間中坐標θ∈Θ)的離散轉化為求解Fredholm積分方程的特征值問題，這一特征值問題如式(2)。

(2)

式中：A1、A2為目標離散空間Ω中的任意兩點坐標，ρ(A1,A2)為這任意兩點處隨機場特征值之間的相關函數值，λi和fi分別為與相關函數對應的特征值和特征函數。關于該特征值的求解具體可參考文獻[31]。因此，隨機場相關函數的特征值計算可以簡化為相關函數特征值的乘積，隨機場即可離散為

(x,y,z∈Ω)

(3)

式中：ξi(θ)為獨立標準正態(tài)隨機向量，θ∈Θ為外部空間坐標；(x,y,z)為隨機場區(qū)域中的任意坐標點，與A1、A2對應；μ、σ分別為隨機場的均值與標準差。為實現高效計算，通常在保證精度的前提下，截取式(3)的前n項來提高計算效率，n的選取參考文獻[31]中比率因子ε的大小進行調整，其定義為

(4)

通常ε大于0.95即可認為隨機場精度已經滿足計算要求。截取后的前n項為

(x,y,z∈Ω)

(5)

基于此，采用K-L級數分解法建立黏聚力與內摩擦角的各向異性隨機場。

采用數值模型及改進二分法，模型中僅考慮黏性土內摩擦角與黏聚力空間變異性對盾構掌子面穩(wěn)定性的影響，其他參數均為常量。保持內摩擦角與黏聚力的均值與確定性計算的參數一致，重點研究內摩擦角與黏聚力共同變異的情況下對掌子面失穩(wěn)模式、極限支護壓力的影響。分別選取黏聚力的變異系數COV(c)為0.1、0.2、0.3、0.4，內摩擦角的變異系數COV(φ)為0.05、0.01、0.15、0.2，組合得到16組工況建立隨機場進行計算分析。為方便后續(xù)分析，工況名命名如表2所示。

表2 工況名稱表Table 2 The names of numerical modelling cases

每組工況分別進行500次隨機計算。由于研究結果表明[18]，抗剪強度參數的負相關關系相對其變異系數對極限支護壓力的影響較小，因此，不考慮抗剪強度參數的負相關。

2.2 掌子面失穩(wěn)模式

以黏聚力變異系數COV(c)為0.4，內摩擦角變異系數COV(φ)為0.2這組工況為例，對掌子面失穩(wěn)的模式進行分析。

選取3組典型工況進行分析。圖5為3組典型工況黏聚力、內摩擦角的隨機場分布。圖6是3組典型工況的最大剪應變增量云圖?？梢园l(fā)現，考慮土性參數的空間變異性時，雖然土性參數隨機場的均值與變異系數相同，不確定性分析計算得到的開挖面失穩(wěn)模式也不盡相同，失穩(wěn)區(qū)域的大小和形式也都有一定的區(qū)別。很顯然，掌子面失穩(wěn)模式與掌子面局部區(qū)域土性參數隨機場的分布相關。幾組工況中，工況a、工況c發(fā)生整體破壞，工況b發(fā)生局部破壞。工況b的掌子面前方內摩擦角隨機場呈現明顯的分層，隧道軸線上方內摩擦角大、下方內摩擦角小，土體沿著強度軟弱處開始破壞，逐漸擴展、傳遞，導致隧道軸線下方出現局部破壞。

圖5 典型工況的黏聚力與內摩擦角隨機場分布Fig.5 Thecohesion and the friction random field

圖6 典型工況的最大剪應變增量云圖Fig.6 Maximum shear strain increment contour

為定量刻畫掌子面局部區(qū)域土性參數隨機場分布與掌子面失穩(wěn)模式間的關系，引入歸一化極限支護壓力σc_nor,掌子面局部歸一化均值cave_nor、φave_nor、標準差cdev_nor、φdev_nor與變異系數COV(c1D)、COV(φ1D)作為量化指標進行探究，這些量化指標按式(6)～式(9)計算。

(6)

(7)

(8)

(9)

式中：cave_1D、φave_1D為c、φ掌子面前方一倍直徑內的均值，cdev_1D、φdev_1D為c、φ掌子面前方一倍直徑內的標準差。cave、φave分別為7 kPa、17°。σc_det為極限支護壓力的確定性計算結果。

表3列出了3組典型工況的歸一化統(tǒng)計參數。不難發(fā)現，歸一化均值指標決定了σc_nor的大小，cave_nor與φave_nor越大，σc_nor越小。而σc_nor與歸一化標準差cdev_nor和φdev_nor、歸一化變異系數COV(c1D)和COV(φ1D)無明顯關系。為進一步說明以上結論，分析500組工況的歸一化極限支護壓力與局部歸一化黏聚力與內摩擦角隨機場的歸一化指標間的關系，繪制圖7。對比圖7(a)～(e)可以發(fā)現，極限支護壓力歸一化指標與σc_nor的相關系數均為負數，σc_nor與歸一化均值指標呈現較明顯的負相關關系，而與歸一化標準差和歸一化變異系數相關系數很小。

表3 典型工況的掌子面局部區(qū)域歸一化參數Table 3 Local normalized parameters for tunnel face of typical cases

圖7 歸一化極限支護壓力與歸一化掌子面局部區(qū)域參數間的關系散點圖Fig.7 Relationship between normalized support pressure and normalized local parameters for tunnel

結合圖5～圖7與表3分析歸一化統(tǒng)計參數對失穩(wěn)模式的影響可以發(fā)現，工況b中φdev_nor與COV(φ1D)較工況a、工況c大，掌子面發(fā)生局部破壞，說明φdev_nor，COV(φ1D)的值能夠一定程度上解釋該工況發(fā)生局部破壞的原因。這是因為，當掌子面局部隨機場分布較為分散時，可能會出現如圖5(b)中的內摩擦角局部隨機場的上下分層情況，這種情況很有可能導致局部破壞。

2.3 極限支護壓力

為進一步分析土性參數的空間變異性對極限支護壓力σc的作用，對比了16組黏聚力變異系數COV(c)與內摩擦角變異系數COV(φ)不相同的工況。

圖8繪制了工況1～工況4以及工況13中σc的概率分布直方圖。所有工況中，σc的中位值與均值均大于σc_det，很顯然土性參數的空間變異性對σc的影響不容忽視。將工況1視作基準工況進行分析。工況4與工況13的COV(φ)與COV(c)分別在基準工況的基礎上增大了3倍。相應地，相比工況1，工況4的σc中位值增大了6.27%，均值增大了5.51%，標準差與變異系數均增大了2倍左右，波動范圍從29.33～39.19 kPa增大到22.68～50.35 kPa；工況13的中位值則增大了2.26%，均值增大了1.29%，標準差與變異系數均增大了1倍左右，波動范圍從29.33～39.19 kPa增大到21.02～45.87 kPa。很顯然，隨著c、φ的變異性增大，σc的統(tǒng)計參數均相應地增大，但φ變異性增大對σc的統(tǒng)計參數造成的影響更甚。分析所有工況σc的統(tǒng)計特征值可以發(fā)現，σc的統(tǒng)計特征值與COV(c)、COV(φ)均成正比關系，c、φ的變異性共同作用時，對σc的影響較某一土性參數單獨作用時更大，會使掌子面的穩(wěn)定性更趨于不安全，所以，對于內摩擦角與黏聚力同時具有變異性的場地地質條件，應當給予重視。

圖8 極限支護壓力頻率分布直方圖Fig.8 Relationship between normalized support pressure and normalized local parameters for tunnel

為進一步量化土性參數空間變異性對σc的影響，將現有的16組工況中的數據進行擬合如圖8，擬合結果表明，σc服從正態(tài)分布，未在圖8中出現的工況的擬合參數見表4。

表4 各工況極限支護壓力正態(tài)分布擬合參數Table 4 Fitting parameters for normal distribution of the critical support pressure for every cases

3 考慮支護壓力變異性的盾構掌子面可靠度分析

由于忽略支護壓力的變異性可能會高估掌子面的穩(wěn)定性，采用蒙特卡洛策略將支護壓力的變異性考慮進盾構隧道掌子面可靠度分析中，同時為將土性參數的變異性考慮進來，將實際的掌子面極限支護壓力視作服從表4中的擬合分布結果的隨機變量。對于支護壓力的概率分布，由于現有的研究中對支護壓力的統(tǒng)計規(guī)律并無統(tǒng)一定論，采用某一工程實測數據[32]的統(tǒng)計結果作為依據進行分析。該工程中土艙壓力的變異系數為0.114，并且符合正態(tài)分布，故支護壓力服從正態(tài)分布，且變異系數的上限設置為0.1。為使結果有效，蒙特卡洛模擬進行100萬次，失效概率定義為

(10)

式中：N為蒙特卡洛模擬的總計算次數，即1×106次，Nf指的是支護壓力小于極限支護壓力的次數。

圖9表示支護壓力服從正態(tài)分布、變異系數為0.1時，支護壓力均值的變化對失效概率的影響。隨著支護壓力均值的增大，失效概率逐漸減小為0%，支護壓力均值小于24 kPa或大于60 kPa，此時土性參數的變異性對失效概率的影響不大。觀察工況1～工況4或工況4、工況8、工況12、工況16的失效概率隨支護壓力均值的變化，當COV(c)與COV(φ)在增大時，失效概率隨支護壓力均值的增大而減小的速率在慢慢變緩，但減緩的路徑各不相同。這種現象在工況1～工況4中觀察得更明顯，在工況4、工況8、工況12、工況16中，隨著黏聚力的變異系數增大，失效概率-支護壓力均值的變化曲線雖然在減緩，但是減緩并不明顯。這說明在考慮支護壓力變異性時，內摩擦角的變異性對掌子面穩(wěn)定性失效概率的影響更甚。通過以上的分析可知，支護壓力的均值對失效概率的影響較大，所以在實際工程中，需要嚴格把控支護壓力的均值來控制掌子面的穩(wěn)定性。為選取合適的支護壓力均值，引入支護壓力均值特征值的概念，借鑒概率統(tǒng)計的方法，將失效概率與支護壓力均值特征值結合起來，在此規(guī)定失效概率為5%時為支護壓力均值特征值σk。在實際應用中，可在盾構的不同時段通過傳感器監(jiān)測土艙壓力，將數據進行正態(tài)分布的擬合，如果該分布中的支護壓力的均值小于支護壓力均值特征值σk，那么應當引起重視。圖9中虛線與曲線的交點代表的就是各工況支護壓力特征值σk。

圖9 支護壓力均值對失效概率的影響Fig.9 The impact of critical support pressure mean value

實際工程中，支護壓力受工程地質條件與人為因素的限制而非定值，具有一定的變異性，接下來分析支護壓力變異系數對掌子面穩(wěn)定性的影響。圖10表示支護壓力變異系數對支護壓力均值特征值σk的影響?？傮w而言，σk隨著支護壓力變異性增強而增大，近似線性關系。對于工況1～工況4，σk都隨支護壓力的變異系數的增大而增大，支護壓力變異系數從0.01增大到0.1，4組工況的σk分別增大了12.88%、9.57%、8.06%、6.49%，增幅隨內摩擦角變異系數增大而逐步減小。這說明，隨著內摩擦角的變異性增強，支護壓力的變異性對σk的影響會被減弱。對于工況4、工況8、工況12、工況16，σk都隨支護壓力的變異系數的增大而增大，但增幅并不隨工況的黏聚力變異系數變化而變化，支護壓力的變異系數從0.01增大到0.1，4組工況的σk分別增大了6.49%、6.15%、5.98%、6.00%。這說明，支護壓力的變異性對σk的影響主要受到內摩擦角變異性的限制，而受黏聚力的變異性影響較小。

圖10 支護壓力變異系數對支護壓力特征值的影響Fig.10 The impact of critical support pressure coefficient of variation on thecritical support pressure characteristic

綜上所述，考慮支護壓力變異性的盾構隧道掌子面穩(wěn)定性分析中，確定合理的支護壓力均值特征值是保證掌子面穩(wěn)定性的關鍵。σk受土性參數與支護壓力變異系數的影響，將σk的值與確定性計算的結果、土性參數、支護壓力的變異系數結合起來，能夠得到式(11)～式(13)。

σk=[k·COV(σ)+b]·σc_det

(11)

k=a1+a2·COV(c)+a3·COV(φ)·ln(COV(φ))

(12)

b=a4+a5·COV(c)+a6·COV(φ)·ln(COV(φ))

(13)

圖10表明σk與支護壓力的變異系數呈線性關系，其截距k、斜率b與土性參數變異系數COV(c)、COV(φ)相關，可以由式(12)與式(13)表示，其中常量a1=1.828、a2=-0.638、a3=2.234、a4=0.844、a5=0.297、a6=-1.251。圖11展示了通過式(11)預測的支護壓力均值特征值σk與通過數值模擬與蒙特卡洛模擬得到的σk的對比，能夠發(fā)現大部分數據落在了1∶1線附近，擬合精度為0.89，表示擬合公式的預測效果不錯。

圖11 σk的數值計算值與式(11)的預測值對比Fig.11 Comparison between the values of σk obtained from numerical modelling and Eq.

4 結論

運用隨機場理論與數值模擬結合的盾構隧道開挖面穩(wěn)定性分析方法，同時考慮了黏性土黏聚力與內摩擦角的空間變異性，進一步考慮了服從正態(tài)分布的支護壓力的變異性，研究了土性參數與支護壓力變異性對盾構隧道掌子面穩(wěn)定性的影響，得到以下結論：

1)同時考慮黏性土黏聚力與內摩擦角變異性時，盾構隧道掌子面的失穩(wěn)模式與掌子面局部區(qū)域的土性參數隨機場分布相關，土性參數在掌子面前方一倍直徑處的局部隨機場均值越小，掌子面越不穩(wěn)定，而土性參數在掌子面前方一倍直徑處的局部隨機場的標準差或變異系數越大，發(fā)生局部破壞的可能性越大。

2)黏性土黏聚力與內摩擦角的空間變異性對盾構隧道掌子面穩(wěn)定性有重要的影響，其中內摩擦角對掌子面穩(wěn)定性的影響更甚。黏聚力與內摩擦角對極限支護壓力有正向的作用，隨著黏聚力與內摩擦角的變異性增大，對掌子面穩(wěn)定性的影響逐漸增強，掌子面越不穩(wěn)定。

3)支護壓力的不確定性對盾構隧道掌子面穩(wěn)定性也有一定的影響，支護壓力均值越大，變異系數越小，掌子面越穩(wěn)定。然而，支護壓力的變異系數對掌子面穩(wěn)定性的影響受到內摩擦角變異性的限制，內摩擦角的強變異性會減弱支護壓力變異性對掌子面穩(wěn)定性的影響。

4)提出了掌子面支護壓力均值特征值的概念，結合掌子面的失效概率、極限支護壓力的確定性結果、土性參數的變異系數以及支護壓力的變異系數，對掌子面支護壓力均值特征值給出了初步的確定方法。