考慮地層三維空間分布的盾構掘進模擬

張超 ，徐智文 ，劉飛香 ，陳仁朋 ?，周蘇華 ，廖金軍 ，胡偉飛

（1.湖南大學 土木工程學院，湖南 長沙 410082；2.建筑安全與節能教育部重點實驗室（湖南大學），湖南 長沙 410082；3.中國鐵建重工集團股份有限公司，湖南 長沙 410100；4.浙江大學 機械工程學院，浙江 杭州 310027）

盾構法是城市地下空間開發的主要施工方法，被廣泛應用于城市地下交通隧道建設.盾構掘進本質上是盾構機與地層動態耦合作用的過程.地質環境復雜多變導致盾構機服役性態難預測，是盾構掘進施工風險的主要誘因，極易引起地表坍塌、開挖面失穩、盾構機姿態偏移和刀盤結泥餅等工程事故.因此，有必要開展盾構掘進模擬以獲取施工過程中裝備受力狀態和地層擾動規律［1-5］.然而，現有盾構掘進模擬通常將地質環境簡化為平行成層地層［6-9］，難以反映盾構服役地質環境的三維空間變異性.

復雜地質環境的定量表征一般通過三維地質建模實現.三維地質建模可提供地層信息的幾何模型，需要借助額外的算法進一步生成有限元網格.然而，地質模型的三維空間變異性給相關算法研究帶來挑戰.目前，一些學者利用地質構造面、地層界面、軟弱夾層等特殊界面對地質模型進行分塊處理，并對獨立封閉塊體進行網格劃分，成功構建了復雜地層數值模型［10-13］.此類地層數值模型主要應用于邊坡穩定性分析、洞室穩定性分析等巖土工程問題［14-16］.目前，尚無三維地質建模技術在盾構掘進模擬中應用的相關報道.

本文提出了一種考慮地層三維空間分布的盾構掘進模擬方法.首先利用Kriging 算法基于鉆孔數據建立三維地質模型，提出一種基于地層分界面的條件判斷算法，克服了傳統方法中復雜地層曲面帶來的有限元網格畸變等問題，實現三維復雜地層的高質量網格劃分；然后，編寫了相應的復雜地層盾構掘進模擬前處理程序，實現復雜地層盾構掘進模擬的參數化建模，大幅提高建模效率；最后，以長沙地鐵四號線某區間為例，討論了所提方法的適用性.

1 三維地質建模

三維地質模型可定量表征復雜地層的地質信息，為構建復雜地層有限元模型提供幾何模型.三維地質建模包括地表及地層分界面的高程插值計算和模型可視化兩部分.

1.1 Kriging算法

在地質統計學理論中，Kriging 算法是對區域化隨機變量進行無偏最優估計的空間插值算法.該算法在礦產空間分布以及地質模型構建方面具有廣泛應用.下面以高程插值為例，介紹該算法的主要插值過程：將待插值點X0的高程定義為Z*（X0），且已知n個鉆孔樣本的數據，并將第i個鉆孔與地表或地層分界面的交點Xi的高程定義為Z（Xi），則Z*（X0）可表示為n個Z（Xi）的加權求和［17］，見式（1）.

式中：λi（i=1，2，…，n）為權重系數.Kriging 插值過程即為通過理論變異函數γ（h）構建的Kriging 方程組［18］求解最優權系數的過程，方程組見式（2）.

式中：γ（hij）是由交點Xi、Xj的水平間距hij確定的理論變異函數值；γ（hi0）是由交點Xi和插值點X0的水平間距hi0確定的理論變異函數值；μ是拉格朗日乘子.此處，理論變異函數是兩點水平間距與高程變異性的函數，可通過擬合實驗變異函數值獲取，見圖1.實驗變異函數值則由鉆孔樣本數據計算得到.球狀變異函數和實驗變異函數值分別見式（3）和式（4）.

圖1 球狀變異函數模型Fig.1 Spherical variogram model

式中：γ（h）為理論變異函數；C0為塊金常數；C為拱高；a為變程；h為樣本點間水平間距.

式中：γ*（h）為實驗變異函數；N（h）表示水平間距為h的鉆孔樣本對數量；Xi與Xi+h分別表示上述鉆孔樣本對與地表或地層分界面的交點.

如此，Kriging 算法首先由鉆孔樣本數據得到離散的實驗變異函數值，并利用球狀模型等數學模型擬合離散值得到理論變異函數，再將理論變異函數代入Kriging 方程組解得插值點的最優權重系數λi，最后根據公式（1）得到該點的高程.

1.2 建模流程

三維地質建模流程具體包括：鉆孔樣本預處理、變異函數求解、地表及地層分界面高程計算、三維地質模型可視化共四個步驟.

1.2.1 鉆孔樣本預處理

鉆孔樣本預處理的目的是保證所有鉆孔的地層類別和層序一致，以便采用統一數據格式對鉆孔樣本的地層信息進行存儲和調用.根據所有鉆孔樣本提供的地層信息，對每個鉆孔揭露的地層按照沉積順序進行編號，見圖2.當鉆孔樣本缺失某地層時，在地層缺失處插入厚度為零的該地層，從而保證所有鉆孔的地層類別和層序相同.

圖2 鉆孔樣本預處理Fig.2 Borehole sample pretreatment

1.2.2 變異函數求解

基于預處理后的鉆孔樣本，采用公式（4）計算地表、各地層分界面的實驗變異函數γ*（h），并利用公式（3）對離散的實驗變異函數值進行擬合，得到理論變異函數γ（h）.

1.2.3 地表及地層分界面高程計算

對地表或地層分界面進行離散，并以離散的待插值點的水平投影作為網格節點，形成圖3 所示虛擬背景網格.以點X0高程插值為例，根據點X0、Xi和Xj的投影位置得到hi0和hij，代入理論變異函數γ（h）計算得γ（hi0）和γ（hij），再代入Kriging 方程組求解得到權重系數λi，求解出X0高程.按上述過程可得到地表及地層分界面上所待插值點的高程數據.

圖3 地層分界面的虛擬背景網格Fig.3 Virtual background grid of stratum interface

1.2.4 三維地質模型可視化

將地表及地層分界面上插值點的平面坐標和高程數據導入自主搭建的可視化平臺，得到三維地質模型，如圖4所示.

圖4 三維地質模型示意圖Fig.4 3D geological model diagram

2 復雜地層網格劃分

由于復雜地層的強空間變異性，構建的三維地質模型通常具有不規則的幾何形態，難以直接生成高質量有限元網格.因此，有必要針對復雜地層開發高效有限元網格劃分算法.

復雜地層盾構掘進模擬的有限元網格劃分難點主要體現在兩個方面.一方面，直接基于地質模型的網格劃分方法在地層尖滅處易導致網格畸變，不利于數值計算收斂，見圖5（a）；另一方面，當盾構隧道穿越圖5（b）所示的復雜地層時，隧道邊界與地層分界面相交并形成若干不規則幾何邊界，這些不規則幾何邊界易形成畸形網格.

圖5 復雜地層盾構掘進模擬網格劃分難點Fig.5 Difficulties in grid division of complex strata

針對上述兩方面問題，本文提出一種基于地層分界面的條件判斷算法.該方法從三維地質模型提取模型幾何邊界和地表及地層分界面高程數據，前者直接輸入網格劃分算法生成未賦予材料屬性的地層網格，后者輸入條件判斷算法以甄別地層網格單元的地層類別并賦予其材料屬性，最終建立復雜地層的有限元模型.

基于地層分界面的條件判斷算法流程見圖6.具體流程如下：

圖6 復雜地層網格劃分流程圖Fig.6 Flow chart of complex formation grid division

1）基于Gmsh編寫地層網格生成腳本，實現地層網格的自動化劃分.該腳本采用六面體單元對地質模型幾何邊界確定的建模區域進行均勻網格劃分，生成初始地層網格，且該網格模型與三維地質模型采用相同坐標系.在初始地層網格中，盾構掘進區域設置為網格協調區，見圖7，得到未賦予材料屬性的地層網格.

圖7 網格協調區示意圖Fig.7 Diagram of grid coordination region

2）讀取上述地層網格數據，計算單元i（i=1，2，…，n；n為單元總數）的中心點坐標（Xi，Yi，Zi）.

3）讀取三維地質模型高程數據，搜索各地層界面上與單元i中心點水平投影重合的點，并將該點的高程命名為Zlj（j為地層界面序號；j=1，2，…，N；N為地層界面總數），其中Zlj＞Zlj+1，見圖8.

圖8 條件判斷示意圖Fig.8 Diagram of conditional judgment

4）對比各地層界面高程值Zlj與單元i中心點高程值Zi，進而確定單元i的地層類別.由地表（Zlj，j=1）開始，判斷Zlj≤Zi是否成立，若成立則判斷單元i位于第j-1 層地層，繼續第5 步；否則令j=j+1，再次判斷，直到滿足Zlj≤Zi.

5）重復步驟2）～4），對地層網格模型中所有單元進行逐一判斷，形成各地層的網格單元集.

6）對各地層單元集賦予相應的土體物理力學參數，完成復雜地層的有限元模型構建.

通過上述步驟可得到圖9 所示的復雜地層有限元模型，該模型包含了三維地質模型的地質信息.

圖9 融合三維地質信息的復雜地層有限元模型Fig.9 Finite element model of complex strata incorporating 3D geological information

3 復雜地層盾構掘進模擬前處理程序

為實現盾構掘進快速參數化建模，基于上述方法采用C#語言編寫了復雜地層盾構掘進模擬前處理程序.該程序包括3 個模塊：復雜地層有限元建模模塊、盾構隧道掘進建模模塊和計算求解模塊，見圖10.復雜地層有限元建模模塊具備三維地質建模功能和地層有限元建模功能，可利用鉆孔數據文件快速建立融合三維地質信息的地層有限元模型.在盾構隧道掘進建模模塊設置隧道幾何與掘進參數后，可自動生成盾構掘進模型的inp文件.計算求解模塊通過批處理命令調用通用有限元求解器ABAQUS，并讀取盾構掘進模型的inp文件進行計算，最終實現復雜地層條件下盾構掘進數值模擬.

圖10 復雜地層盾構掘進模擬前處理程序Fig.10 Pre-processing program for simulation of shield tunneling in complex strata

在復雜地層盾構掘進模擬前處理程序輸入隧道幾何和掘進參數，即可實現盾構掘進的快速參數化建模.隧道幾何設置涵蓋隧道埋深、隧道直徑、管片尺寸等參數.隧道掘進參數設置包括掘進步長、掌子面支護系數、管片和土體材料等參數.

基于上述輸入參數，盾構隧道掘進建模模塊可自動完成以下四方面設置，并生成inp文件：

1）分段開挖設置.在每個掘進步中，采用“生死”單元法移除掌子面前方土體單元，同時激活盾構管片單元，實現分段開挖過程的模擬.

2）掌子面支護力設置.將掌子面的靜止土壓力作為支護力默認值，沿高度方向呈梯度分布.每個掘進步可根據實際施工資料設置掌子面支護系數，將該系數乘以默認支護力得到實際支護力.

3）地層損失設置.盾構掘進引起的地層損失是導致地表沉降的主要因素，常以地層損失率［19］表示，并可通過整環擾動層近似模擬［20］.通常盾體與圍巖間的間隙具備不均勻性，即拱頂間隙大于兩旁和底部間隙［21］.為近似考慮該不均勻性，本文將擾動層分為上下兩半環分別設置材料參數，見圖11.上半環擾動層參數取上部圍巖參數折減值；下半環擾動層采用鄰近圍巖的材料參數.

圖11 地層損失設置示意圖Fig.11 Diagram of stratum loss setting

4）材料設置.復雜地層盾構掘進模擬前處理程序包括線彈性本構模型和摩爾庫倫本構模型，前者用于管片和上半環擾動層，后者用于土體和下半環擾動層.

4 工程案例驗證

本節以長沙地鐵四號線區間隧道施工工程為例，利用前處理程序構建盾構地層掘進模型，并開展掘進模擬.將沿隧道軸線的地表沉降模擬結果與實測數據對比，探討本文所提方法的有效性.

4.1 工程概況

長沙地鐵四號線桐梓坡—望月湖區間隧道全長1 512 m，該區間段采用土壓平衡盾構機掘進施工，盾構機外徑為6.28 m.該區段為雙線隧道.本文選取先期施工的左線隧道掘進引起的地表沉降數據進行驗證分析.隧道位置與鉆孔平面分布見圖12.依據鉆孔樣本數據，利用所編寫的前處理程序構建該區域的三維地質模型，見圖13.

圖12 隧道與鉆孔分布平面圖Fig.12 Plan view of tunnels and borehole distribution

圖13 三維地質模型Fig.13 3D geological model

圖14 為沿隧道縱軸線方向的地質剖面圖.該區段內地層復雜多變，自上而下分別為素填土、硬塑黏土、強風化板巖、中風化板巖和微風化板巖.各地層的物理力學參數見表1.

表1 地層物理力學參數表Tab.1 Physical and mechanical parameters of soils

圖14 地質剖面圖Fig.14 Longitudinal geological profile

4.2 盾構掘進模擬與驗證

為研究地層空間變異性對地表沉降的影響，采用表2 所示盾構隧道幾何和掘進參數，將其輸入復雜地層盾構掘進模擬前處理程序建立圖15 所示盾構掘進有限元模型，并開展盾構掘進模擬.

表2 盾構隧道幾何和掘進參數表Tab.2 Geometry and excavation parameter of shield tunnel

圖15 盾構掘進模擬有限元模型Fig.15 Shield tunneling simulation finite element model

為驗證有限元模型的準確性，選取掘進斷面內極軟巖占比為50%與100%的兩種工況，分別對應711 環和730 環地表沉降監測點，并對比兩種工況下地表沉降變化過程的模擬結果與實測結果，見圖16.根據盾構施工工法，盾構掘進對監控斷面地表沉降的影響分為兩個階段，盾構掘進抵達監控斷面前稱為第一階段，盾構通過監控斷面后稱為第二階段.在地表沉降規律方面，兩種工況下，兩階段的模擬結果與實測結果趨勢均一致.在第一階段期間，地表沉降均隨掌子面與監控斷面的距離減小而逐漸增大，而在第二階段期間，地表沉降增長變緩且最終趨于穩定.在沉降幅值方面，兩種工況下地表產生了超過10 mm 的差異沉降量，具體為：當盾構機由掘進斷面部分極軟巖地層進入全斷面極軟巖地層時，地表最終沉降量模擬與實測分別增加了11.67 mm 和18.02 mm，模擬與實測的誤差為35%.誤差主要來源于730環，根據施工資料，盾構機開挖至730環停機約10 h，而有限元模型無法考慮停機時未及時注漿引起的地層損失，因此730 環最終沉降模擬值小于實測值.綜上，有限元模擬結果與實測結果趨勢一致，吻合程度好，表明有限元模型計算準確.

圖16 監測點地表沉降變化過程模擬與實測結果對比Fig.16 Comparison of simulation and measured results of surface subsidence change process at monitoring points

為研究地層空間變異性對地表沉降的影響，將各監測點最終沉降的模擬結果和實測數據繪制于圖17 中.模擬結果顯示，地表最終沉降隨掘進斷面內極軟巖占比增大而顯著增加，與實測數據規律一致.綜合圖16～17 結果，掌子面接近監測斷面前，開挖卸載作用引起掌子面前方極軟巖地層發生變形，導致地表發生第一階段沉降.掌子面通過監測斷面后，極軟巖地層向隧道斷面內發生位移，導致第二階段地表沉降.當極軟巖占比增加時，上述兩階段的地表沉降均增大，進而導致地表最終沉降顯著增加.

圖17 沿隧道軸線地表最終沉降模擬與實測結果對比圖Fig.17 Comparison of simulation and measured results of final settlement along tunnel axis

綜上，復雜地層盾構掘進模擬前處理程序可實現盾構掘進快速參數化建模，所建有限元模型能較準確地反映復雜地層條件下盾構掘進過程，并能有效揭示地層空間變異性對地表沉降的影響規律.

5 結 論

本文提出了考慮地層三維空間分布的盾構掘進模擬方法，編寫了相應的復雜地層盾構掘進模擬前處理程序，并以長沙地鐵四號線區間隧道工程為例，驗證了所提方法和相應程序的有效性.主要結論如下：

1）基于地層分界面的條件判斷算法可以規避復雜地層曲面的有限元網格畸變問題，并且有效解決了復雜地層曲面與隧道斷面的網格協調問題.

2）復雜地層盾構掘進模擬前處理程序集成了三維地質建模功能和復雜地層網格劃分算法，實現了復雜地層掘進模擬的參數化建模，大幅提高了建模效率.

3）盾構掘進模擬結果表明，地表沉降量隨盾構掘進斷面內軟巖占比增大而顯著增加，與工程實測數據規律一致，驗證了本文所提方法的有效性.