黃德龍 宗鐘凌 劉強 岑航 湯愛平

1.江蘇海洋大學土木與港海工程學院 連云港222005

2.哈爾濱工業(yè)大學土木工程學院 150090

引言

城市地下綜合管廊的建設可以有效解決直埋管線所帶來的“馬路拉鏈”、“空中蜘蛛網(wǎng)”等“城市病”［1］。截止到2020 年，中國已有25 個試點城市開始推進管廊的建設并逐步投入使用。這些城市，有50%位于Ⅶ度設防區(qū)，有20%位于Ⅷ度設防區(qū)，因此管廊的抗震研究勢在必行。

在1995 年日本阪神大地震時，地鐵車站等地下結構遭遇嚴重毀壞，學者們才開始關注地下結構的抗震問題［2］。管廊等地下結構體系的破壞主要以外壁混凝土開裂以及內(nèi)容物破壞為主。其破壞的直接原因主要是土-結相互作用（SSI），與地面建筑物所產(chǎn)生的慣性效應不同［3，4］。因此，場地振動、地震引起土體失效以及土-結相互作用是現(xiàn)階段研究地下結構抗震的熱點問題。近些年，國內(nèi)外學者們對管廊等地下結構的抗震性能進行了大量動力試驗和數(shù)值模擬，部分學者還進行了現(xiàn)場震害調(diào)查，提出了一系列理論和技術方法。

對于試驗研究，振動臺試驗是主要采取的測試手段［5-11］，得出了均勻土體中管廊的地震響應規(guī)律以及管廊-土-地面結構（TSF）體系的響應規(guī)律，并對一致和非一致激勵下土-結相互作用機理進行了深入研究。也有一些學者利用離心機試驗進行地下結構抗震的研究［12-16］，探究了管廊和地震波之間的相互影響機制、土體致密化響應特點，并同樣對土-結相互作用進行了深入探索。對于數(shù)值模擬研究，學者們同樣主要關注均勻土或豎向成層土內(nèi)的土-結相互作用、管廊及節(jié)點的地震響應、管廊抗震設計簡化計算方法以及管廊-土-地面建筑之間的動力相互影響［4，17-25］。對于震害現(xiàn)場調(diào)查研究［1］，主要是針對地震作用下管廊的破壞形式進行總結和分析。

綜合管廊可跨河床或不同地質(zhì)單元，并且只要不是膨脹土，就可選用挖出土回填，造成管廊周圍土體的水平非一致性。上述試驗、數(shù)值模擬及現(xiàn)場震害調(diào)查中，鮮有學者考慮到穿越水平不均勻土體的管廊地震響應機理問題，并且多是運用有限元理論，無法精確模擬土顆粒與結構的接觸，更無法模擬土顆粒遷移規(guī)律。本文基于振動臺試驗以及改進離散元對上述不足進行了探究。

1 試驗與模擬方法

本文首先利用振動臺試驗得到管廊的動力響應，并將試驗與改進離散元方法的結果進行對比，以驗證改進離散元方法的適用性。

1.1 振動臺試驗

本次振動臺試驗在哈爾濱工業(yè)大學結構與抗震實驗中心進行。試驗首先預制縮尺管廊，見圖1，然后采用單向振動臺，并利用疊層剪切箱，將砂土與黏土分兩側放入剪切箱，管廊埋置深度為0.4m。針對El-Centro 波、Kobe 波和人工波，分別經(jīng)過調(diào)幅（0.1g／0.2g／0.4g／0.8g／1.0g）進行加載。試驗中分別監(jiān)測上覆土體的豎向位移、管廊外壁土壓力以及管廊的應變，其中在砂土和黏土內(nèi)各敷設7 個土壓力計，為E-1～E-14，用來測量管廊側壁（E-1～E-3 處于砂土，E-4～E-6 處于黏土）受到的土壓力增量以及管廊上（E-10～E-14）、下（E-8～E-10）的土壓力增量。具體的試驗設置方案參照黃德龍等［26］關于穿越水平非均勻場地管廊振動臺試驗。

圖1 縮尺管廊試驗前準備Fig.1 Preparation work before scale utility tunnel test

1.2 離散元模擬

1.改進離散元方法基本原理

離散元（PFC）的計算思想是使每個元素均滿足牛頓第二定律，用中心差分的方法求解各元素的運動方程，得到研究對象的整體運動形態(tài)，而地震波對結構的作用是一個能量不斷積聚的過程，和PFC的思想具有一定的差異性。因此本研究中改進的離散元方法是基于能量的思想，利用微分正交法（DQM）和Newmark-β 方法來改進離散元，提高計算精度和計算效率。

綜合管廊外壁的總勢能V是管廊壁的應變能U、動能K和外部土體做的功W之和，應變能可以寫成：

動能可表示為：

由于土體的壓力而產(chǎn)生的外部功可以通過下式獲得：

使用Hamilton原理推導出結構的控制方程，原理如下：

為了解決工程領域中出現(xiàn)的初始值或者邊界值問題，有很多數(shù)值方法，其中一個較好的數(shù)值方法是DQM法。與其他數(shù)值方法相比，DQM 法有以下幾個優(yōu)點：（1）是一種解決非線性微分方程導數(shù)近似的精確方法；（2）可以滿足各種邊界條件，并且需要更少的方程和編程工作；（3）準確性和收斂性都很高?；谏鲜鰞?yōu)點，近年來，DQM法在結構和動力學問題分析的數(shù)值求解中越來越受歡迎。基于此方法，結構控制方程的矩陣形式［27］可以寫成式（5）：

式中：KL、KNL、C、M、db和dd分別代表線剛度矩陣、非線性剛度矩陣、阻尼矩陣、質(zhì)量矩陣、邊界點和域點。

利用Newmark-β方法來獲得結構的時程響應序列，公式（5）可以改寫如下［28］：

基于DQM 及Newmark-β 的計算原理，對PFC2D程序進行二次開發(fā)，以求解式（6）非線性離散元運動方程。PFC2D程序中開發(fā)自定義求解原理，需先對模型設置特有的編號，避免與標準模型沖突。在C++程序中，依次修改源文件中CM_vep：：CM_vep（）、CM_vep：：Name（）、CM_vep：：PropNames（）、CM_vep：：ReturnProp（）、CM_vep：：AcceptProp（）、CM_vep：：SaveRestore（）等成員函數(shù)，完成源文件的修改。最后，將源文件編譯成動態(tài)鏈接庫文件（DLL 文件），可供PFC2D程序調(diào)用。

2.離散元方法的設置

圖2 為管廊的PFC2D離散元模型，此模型融入管廊及其土體的Newmark動力求解方法，基于黃景崎等［29］的研究，在模型外側施加等效荷載來實現(xiàn)試驗中疊層剪切箱的作用效果。對于模型的構建，Potyondy等［30］指出：由于離散元對土體粒子數(shù)目無法按照真實情況來模擬，雖然土體顆粒的尺度因子對離散元是有影響的，但是當L'／d＞32，即剪切箱的最大邊長與顆粒粒徑之比大于32時，此時宏觀響應與顆粒尺寸幾乎無關，即此時顆粒的大小對管廊的響應幾乎沒有影響。所以模擬時為了減少粒子數(shù)目和計算所需要的時間，采用比較大的球形顆粒對土粒進行建模。模型分為兩個區(qū)，接近管廊的區(qū)域顆粒小，其最大的L'／d的值為180；離管廊較遠的區(qū)域為顆粒較大區(qū)，其最大的L'／d為125，兩者都大于32，所以此模型可用來模擬管廊及周圍土體的地震響應。本研究粒子數(shù)為14733。

圖2 管廊的二維離散元模型Fig.2 2D discrete element simulation model of utility tunnel

整個布置區(qū)域與試驗振動臺的豎向截面尺寸一致，其中管廊和土箱分別采用wall來模擬。管廊周圍相近的土體孔隙比都采用0.15，砂土采用顆分的形式定義，阻尼比為0.3，黏性土的阻尼比為0.45。對于構建砂土模型來說，顆粒與顆粒之間采用線彈性模型，如圖3a。此模型需要土體的剛度系數(shù)，以土體的基床系數(shù)為依據(jù)［31］，基于本試驗的土體物理性質(zhì)，法向剛度和切向剛度取為3.3 ×104N／m，滑動摩擦系數(shù)取為0.4；構建黏土模型顆粒與顆粒之間同樣采用線彈性模型，其中法向剛度取為2.8 ×

圖3 離散元粒子與管廊接觸模型Fig.3 The contact models between discrete element particles and utility tunnel

104N／m，切向剛度取為2.0 ×104N／m，滑動摩擦系數(shù)取為0.4；對于砂土與管廊之間的接觸（即ball 與wall 接觸）采用線性平行粘結（linearpbond）模型，如圖3b。基于式（3）中的土壓力，可以計算接觸法向剛度和切向剛度都取為2 × 106N／m，摩擦系數(shù)取0.4，阻尼比取為0.2；對于黏土與管廊之間的接觸，也同樣采用線性平行粘結模型，接觸法向剛度和切向剛度都取為1.0 ×106N／m，摩擦系數(shù)取0.4，阻尼比均取為0.4。

3.離散元與試驗結果對比

基于DQM 計算方法，加速度、位移和應變的響應結果與試驗對比如圖4 所示?？梢园l(fā)現(xiàn)，加速度和位移的計算結果吻合較好。

圖4 離散元模擬和試驗結果對比（El-Centro，PGA ＝1.0g）Fig.4 Comparison of numerical simulation and experimental results（El-Centro，PGA ＝1.0g）

在模擬中，由于是二維模型并且土體顆粒尺寸與真實顆粒相差較大，導致顆粒之間的剛度和阻尼無法按照真實情況模擬。這種情況下，粒子的宏觀運動特性可以準確實現(xiàn)，微應變等微觀變形響應的時程曲線將很難與試驗準確對應，但是其幅值均在可接受范圍內(nèi)。

2 試驗與模擬結果討論

2.1 豎直位移

圖5 為輸入PGA 從0.1g 調(diào)幅至1.0g（增幅0.2g）的三種地震波后，利用激光位移計所測得的砂土和黏土豎向位移。可以看出砂土數(shù)值增大明顯，砂土區(qū)發(fā)生了剪沉現(xiàn)象（沉降），而黏土區(qū)數(shù)值呈減小趨勢，發(fā)生了剪脹現(xiàn)象。隨著加速度的增大，剪沉和剪脹增加量逐漸減小，即剪沉和剪脹雖然逐漸累積，但累積速度在放緩。到最后加載PGA ＝1.0g的人工波后，發(fā)現(xiàn)砂土依然在剪沉，黏土同樣也發(fā)生了輕微剪沉現(xiàn)象，但是其量值很小。說明隨著加載的進行，砂土的密度逐漸增大，而黏土的密度逐漸減小到一定值后還有增加的趨勢。對于三種波在PGA ＝0.8g時出現(xiàn)了與其他地震動強度不同的情況，即砂土出現(xiàn)了明顯的變形上浮，主要是因為三種波在PGA ＝0.4g均發(fā)生了明顯的剪沉。

圖5 不同土體表面的豎向位移Fig.5 Vertical displacement of different soil surfaces

土體變形將會導致土壓力發(fā)生變化，而管廊所受的土壓力又是管廊發(fā)生變形的直接因素。圖6 為試驗中砂土發(fā)生了剪沉現(xiàn)象，可以發(fā)現(xiàn)砂土流向與剪切箱接觸邊的中心，并逐漸向下流動。

圖6 試驗觀測砂土沉降Fig.6 Sand settlement observation

圖7 為土體發(fā)生豎向位移的改進離散元二維截面模擬情況。此程序無法實現(xiàn)內(nèi)部管廊的轉動作用，但是從圖7a 中可以發(fā)現(xiàn)左側砂土區(qū)發(fā)生了剪沉。土體與管廊之間出現(xiàn)了明顯的空隙，之后砂土區(qū)左側土體補充此空隙，但是有明顯的滯后效應。這與圖6 砂土區(qū)剪切箱邊砂土發(fā)生的流動變化相一致。并且圖7b 土顆粒的運動也可以說明砂土主要是向下沉降，而黏土是向側壁擴展，進而導致剪脹現(xiàn)象的發(fā)生。

圖7 土體豎向移動模擬Fig.7 Simulation of vertical movement of the soil

2.2 管廊上部土壓力增量

圖8 為El-Centro 波，PGA ＝0.4g 和1.0g，分別對應砂土和黏土區(qū)20cm和40cm深度處（管廊上部）的土壓力增量ΔσE時程。可以發(fā)現(xiàn)深度越深，土壓力增量越大；小震時砂土區(qū)的土壓力較大，而大震時，黏土區(qū)的土壓力較大；其中部分土壓力增量迅速達到一個平衡階段，然后在該位置震蕩直到結束，見圖8a PGA ＝0.4g 的曲線，這與Cilingir 和Madabhushi［32］發(fā)現(xiàn)的規(guī)律一致。Hushmand等［14］給出的解釋是土體致密化會導致土壓力增量增大，土壓力最后可能有殘余值。結合圖5也可以發(fā)現(xiàn)，隨著加載的進行，在PGA 增大的同時，砂土和黏土各自發(fā)生剪沉和剪脹現(xiàn)象，其殘余內(nèi)力在增大，表現(xiàn)為土壓力增量在增大。圖9 為利用改進的離散元程序對加載1.0g El-Centro波的土體內(nèi)部不平衡力變化過程進行模擬，此不平衡力可以反映土壓力的大小。由于El-Centro波在3s左右加速度達到最大峰值，因此對模擬的前3s中6 個時刻的土壓力進行研究（限于篇幅，圖9 中僅給出3 個代表時刻）?？梢园l(fā)現(xiàn)隨著加載時刻的進行，加速度逐漸增大，砂土區(qū)管廊與土體之間的空隙逐漸增大。由于上部土體的補充，此空隙逐漸減小?？梢钥闯鐾翂毫φw在逐漸增大。開始時，由于砂土區(qū)下部土體沉降，此區(qū)域土體內(nèi)部不平衡力逐漸增大，即土壓力增大。隨著加速度的增大，管廊上覆土體的土壓力也在增大，達到加速度峰值1.0g時，管廊上部土壓力超過下部，在邊界以及兩種土體分界區(qū)土壓力最大。

圖8 不同土體的土壓力增量（El-Centro 波）Fig.8 Earth pressure increments for different soils（El-Centro）

圖10 表示利用改進的離散元程序對PGA ＝1.0g時三種波作用下的土壓力峰值，可以發(fā)現(xiàn)Kobe波的土壓力峰值最小。

圖10 三種地震波PGA ＝1.0g 時土壓力峰值模擬結果（單位：kPa）Fig.10 Simulation of peak earth pressure at 1.0g for three waves（unit：kPa）

2.3 管廊側壁土壓力增量

圖11 分別為砂土區(qū)和黏土區(qū)管廊側壁沿縱向受到的土壓力增量ΔσE時程，可以發(fā)現(xiàn)砂土區(qū)土體對側壁土壓力增量略大。并且側壁中間（砂土E-2 和黏土E-5）土壓力增量較大，其原因是處于側壁中間處的土體受到側壁的約束，殘余內(nèi)力無法自由釋放，而導致對側壁土壓力的增大。同樣可以發(fā)現(xiàn)砂土區(qū)E-2 發(fā)生了土壓力增量的負向殘余值，其原因是縱向運動的過程中，砂土與管廊有垂直側壁方向的脫離趨勢。

圖11 沿豎直方向分布土壓力增量響應（El-Centro，PGA ＝1.0g）Fig.11 Response of earth pressure increments along vertically distributed（El-Centro，PGA ＝1.0g）

柔度比F表示結構周圍土體與結構之間的相對剛度［2］。Wang［3］和Hashash等［4］都對矩形地下結構的柔度比進行了深入研究，其計算公式如下：

式中：Gm為在自由場中沿管廊高度土體的平均應變所對應的剪切模量；K 為管廊的剛度系數(shù)，可通過結構力學的方法確定，即在限制管廊底部水平運動的同時，在管廊頂部施加單位力所引起的管廊頂部橫向位移的倒數(shù)；W為管廊截面的寬度；H為管廊截面的高度。

本研究的工況只有砂土和黏土兩種，只涉及兩個柔度比，即砂土區(qū)柔度比F1和黏土區(qū)柔度比F2。由于本研究砂土的剪切模量Gm-sand近似為7.7MPa，黏土的剪切模量Gm-clay近似為3.5MPa，而兩個工況柔度比所涉及的其他參數(shù)均相同，因此F1＞F2。

圖12 為砂土和黏土區(qū)沿管廊高度分布土壓力增量峰值變化曲線?？梢园l(fā)現(xiàn)隨著PGA 的增大，土壓力增量峰值有增大的趨勢，但是由于土體和管廊有分離的趨勢，部分土壓力增量減小。小震PGA ＝0.1g時，在側壁中間高度處，土壓力增量都要比PGA ＝0.2g時大，這是由于地震波以PGA ＝0.1g開始加載，當突然加載或波形（頻率）發(fā)生改變時，土壓力均將發(fā)生顯著變化。大震時砂土區(qū)土壓力增量要比黏土區(qū)大，即柔度比在一定范圍內(nèi)越大時，管廊相對于土體越柔，受到的土壓力越大。還可以發(fā)現(xiàn)隨著加速度的增大，土壓力增量分布近似由底部（頂部）較大向中間位置較大變化。Hushmand 等［14］指出隨著地震作用的增大，土壓力增量的分布近似由三角形變?yōu)閽佄锞€形，柔度比可以決定側壁土壓力的形狀［2］。本文得出的結果與上述結論相符，但是兩者不同點在于砂土區(qū)管廊下部的土壓力增量峰值要大于上部，而黏土區(qū)則相反。因此砂土區(qū)要著重關注管廊下部，而黏土區(qū)上部較危險，這與2.1 節(jié)土體豎向變形有密切關系。

圖12 沿管廊側壁高度土壓力增量峰值變化Fig.12 Variations of earth pressure increment peaks along the height of utility tunnel

在地震作用下，地基土體會發(fā)生沉降變形。管廊地基土體的沉降變形是破壞管廊穩(wěn)定性的直接原因。土體一旦出現(xiàn)沉降變形就會使管廊受力發(fā)生變化。受力超過允許值，管廊將遭到破壞。由地震作用所產(chǎn)生的管廊側壁水平向動土壓力為［33］：

這表明擋板上的土壓力水平分量py隨深度增加趨于一個極值，而傳統(tǒng)的Coulomb 土壓力理論和Rankine土壓力理論都認為土壓力沿深度線性增加，這說明地震動力作用下的擋板土壓力分布與Coulomb土壓力理論和Rankine 土壓力理論有較大差別。這與圖12 管廊側壁土壓力呈現(xiàn)中間大兩端小的極值曲線相接近。

分析式（8）可以發(fā)現(xiàn)，土壓力的極值與側壓力系數(shù)K、土體重度γ、土體內(nèi)摩擦角φ、擋板與土體摩擦角δ、地震加速度系數(shù)kv、kh和剪切箱-管廊側壁凈距L有關。

2.4 地表豎向位移與管廊側壁土壓力的關系

圖13 為試驗中地表豎向位移增量與管廊側壁土壓力增量的關系?？梢园l(fā)現(xiàn)，砂土區(qū)隨著豎向沉降增量峰值的增大，側壁土壓力增量在開始階段較大，砂土區(qū)側壁中間E-2 的土壓力增量在沉降增量為0.4mm～0.8mm 附近出現(xiàn)極小值。說明在發(fā)生沉降初期，側壁土壓力變化明顯，側壁中間土壓力增量最大，并且存在極小值。而黏土區(qū)的剪脹與土壓力增量關系存在極大值，當剪脹量在0.2mm～0.3mm 附近時，側壁土壓力增量出現(xiàn)極大值。

圖13 管廊側壁土壓力增量峰值與地表豎向位移之間的關系Fig.13 Relationship between the vertical displacement of the ground surface and peak value of the earth pressure on the side wall of the tunnel

圖14 為利用改進的離散元程序對PGA ＝1.0g時El-Centro波作用下砂土和黏土區(qū)管廊側壁土壓力峰值的模擬。可以發(fā)現(xiàn)砂土區(qū)土壓力整體要比黏土區(qū)土壓力大，并且側壁中部土壓力大于側壁上下兩側的土壓力，與圖12 曲線相一致。還可以發(fā)現(xiàn)管廊側壁的土壓力要比管廊上下壁土壓力大，管廊的破壞主要受側壁土壓力控制。

圖14 El-Centro 波管廊側壁土壓力峰值模擬（PGA ＝1.0g）（單位：kPa）Fig.14 Simulation of peak earth pressure corresponding to El-Centro wave（PGA ＝1.0g）（unit：kPa）

3 結論

本文基于振動臺試驗以及改進的離散元數(shù)值方法對水平非均勻場地綜合管廊地震響應進行探究，通過利用改進的離散元程序對管廊及其周圍土體進行模擬，并與試驗進行對比，得出如下結論：

1.基于能量思想，并利用微分正交法以及Newmark-β方法得到模型體系的運動控制方程，以此改進離散元程序，建立PFC2D管廊模型。模擬結果與試驗結果的誤差均在20%以內(nèi)，即二次開發(fā)的離散元數(shù)值模擬方法是準確的；

2.砂土區(qū)土體發(fā)生了剪沉，黏土區(qū)發(fā)生了剪脹現(xiàn)象，并且伴隨著土體密度的變化，管廊發(fā)生了輕微轉動，砂土區(qū)剪沉現(xiàn)象的發(fā)生是為了補充管廊下部砂土與管廊的空隙；

3.土體致密化會導致土壓力增大，土壓力增量可能會出現(xiàn)殘余值，土壓力的大小與所輸入PGA有關，不均勻土體對土壓力的影響較大。管廊側壁的土壓力增量呈中間大兩端小的極值形，并且在一定范圍內(nèi)，柔度比越大，管廊所受到的土壓力越大，管廊側壁變形也越大；

4.砂土區(qū)在發(fā)生沉降初期，側壁土壓力變化明顯，且之后存在極小值；黏土區(qū)發(fā)生剪脹，隨著剪脹量的減小，土壓力增量存在極大值。