丁祖德，陳譽升，資 昊

（昆明理工大學建筑工程學院，云南昆明 650500）

引言

動力時程分析法能夠反映巖土體－結(jié)構(gòu)動力相互作用過程以及材料的非線性行為，已廣泛應用于隧道及地下結(jié)構(gòu)的地震響應分析中［1］。在開展地下結(jié)構(gòu)動力時程分析時，需要考慮無限域地層與結(jié)構(gòu)的相互作用，但受限于計算資源，通常只截取有限域進行研究。為使有限域地層具備模擬無限域地層的輻射阻尼［2］，需要在人工截斷的地層周圍施加人工邊界。此時，地震動輸入方法及動力人工邊界的選擇，將直接影響地下結(jié)構(gòu)動力響應結(jié)果精度和可靠性。

目前，國內(nèi)外學者針對動力人工邊界以及地震動輸入方法開展了較豐富的研究工作，發(fā)展了如粘性邊界［3］、粘彈性邊界［4］、等位移邊界［5］等多種人工邊界模型。劉晶波等［6］推導了與波動輸入法相匹配的二維粘彈性邊界。楊喻聲等［7］驗證了等位移邊界在彈性基巖假設條件下的有效性，簡化了常規(guī)波動分析方法。李嘉瑞等［8］對比分析了采用不同人工邊界對地下結(jié)構(gòu)動力計算結(jié)果的影響，建議側(cè)邊界采用綁定邊界或豎向約束邊界，根據(jù)下臥基巖場地條件選擇底部邊界類型。趙武勝等［9］和馬笙杰等［10］分別研究了無限元邊界、粘彈性邊界與地震動輸入方法之間的影響機制，結(jié)果表明加速度時程輸入將導致底部人工邊界失效，且導致位移場與理論值產(chǎn)生嚴重偏差。

已有研究基本上都是從單一人工邊界的地震動輸入原理展開闡述，且大都基于自由場模型，從位移角度來驗證人工邊界及地震動輸入的有效性［9－11］，而從結(jié)構(gòu)力學響應角度進行驗證［8，12］，對地震動輸入方法與動力人工邊界的匹配性問題的探討較少。此外，人工邊界的施加較為繁瑣，在隧道動力響應有限元模型前處理過程中的工作量巨大?；诖?，文中依托ABAQUS 有限元軟件平臺，結(jié)合粘彈性邊界、無限元邊界和等位移邊界等動力人工邊界理論，基于Python 編程開發(fā)相應人工邊界下地震動輸入程序，以實現(xiàn)模型邊界及地震動的自動施加。進一步結(jié)合自由場模型和土－隧道相互作用模型，對比振動輸入和波動輸入方法及常用人工邊界的性能，驗證文中開發(fā)程序的有效性，最后探討2種地震動輸入方法的差異性。

1 人工邊界及地震動輸入的程序?qū)崿F(xiàn)

目前，隧道動力響應計算模型中，常用的人工邊界有遠置邊界、無限元邊界、粘彈性邊界和等位移邊界等［8］，常用的地震動輸入方法有波動輸入法和振動輸入法［13］。其中，波動輸入法是基于彈性基巖假設，并結(jié)合動力人工邊界理論，把地震動荷載轉(zhuǎn)換為等效節(jié)點力后逐一施加到邊界節(jié)點上。振動輸入法是基于剛性基巖假設，通過在計算模型的底邊界上施加位移、速度或加速度時程來實現(xiàn)地震動輸入。

波動輸入法在施加人工邊界時，需要逐一處理邊界節(jié)點，前處理工作量巨大。鑒于此，文中結(jié)合粘彈性邊界、無限元邊界、自由度綁定邊界和豎向約束邊界的基本理論，采用Python 編程開發(fā)上述人工邊界的自動施加程序。為方便使用，進一步將該程序內(nèi)嵌至ABAQUS 軟件中，新增了專門的隧道動力響應輔助計算模塊（注：原ABAQUS 軟件中無此模塊，文中采用Python 語言編程，結(jié)合其二次開發(fā)功能實現(xiàn)），實現(xiàn)了程序界面操作功能，如圖1 所示。接下來，具體闡述不同邊界條件下波動輸入法的程序?qū)崿F(xiàn)流程。限于篇幅，人工邊界理論推導詳見文獻［5－6］、［14］。

圖1 人工邊界及地震動輸入自動化程序Fig.1 Automatic program of artificial boundary and ground motion input

1.1 粘彈性邊界下的地震動輸入

粘彈性邊界通過在模型邊界處設置并聯(lián)的阻尼器元件和彈簧元件來模擬無限域地基的輻射阻尼和彈性恢復力，避免了高頻失穩(wěn)及低頻漂移的問題。針對粘彈性邊界，常采用波動輸入法將地震波轉(zhuǎn)化成等效節(jié)點荷載輸入。其等效節(jié)點力計算如式（1）［10］：

式中：Fb為邊界節(jié)點上的等效荷載；K、C分別為施加到邊界上的彈簧—阻尼元件的剛度系數(shù)、阻尼系數(shù)，其計算公式與取值參考馬笙杰等［10］研究；Ab為節(jié)點等效面積，分別為邊界節(jié)點的自由場位移、速度及應力。

考慮到等效節(jié)點荷載是基于邊界節(jié)點逐個施加的，且每個節(jié)點處均需要計算對應的等效節(jié)點力時程，文中采用Python語言編制粘彈性邊界及地震動輸入程序，程序流程如圖2所示。

圖2 粘彈性邊界實現(xiàn)及地震動輸入流程Fig.2 Implementation of viscoelastic boundary and seismic input process

1.2 無限元邊界下的地震動輸入

無限元邊界是基于無限單元實現(xiàn)的靜動力人工邊界，在靜力分析時通過有限長的無限單元來模擬半無限地基，能反映地基無窮遠處位移為零的條件，在動力分析時通過在無限單元中內(nèi)置阻尼器來實現(xiàn)粘性邊界的吸收效果［14］。

無限元邊界與粘彈性邊界的地震動施加原理相似，同樣采用Python 語言編程實現(xiàn)無限元邊界及地震動等效節(jié)點力時程的輸入，程序流程如圖3所示。等效節(jié)點荷載計算式如下：

圖3 無限元邊界的地震動輸入流程Fig.3 Flow chart of ground motion input with infinite element boundary

式中：Fb為邊界節(jié)點上的等效荷載；Ab為節(jié)點的等效面積；C為阻尼系數(shù)分別為邊界節(jié)點的自由場速度及應力。

1.3 等位移邊界下的地震動輸入

等位移邊界分為自由度綁定邊界和豎向約束邊界。其中，自由度綁定邊界通過在模型兩側(cè)設置運動綁定約束，約束兩側(cè)邊界上節(jié)點的自由度防止其旋轉(zhuǎn)，迫使土層側(cè)邊同一高度位置的節(jié)點同時運動，從而避免側(cè)邊界對波的反射［5］。豎向約束邊界通過約束兩側(cè)截斷邊界的豎向自由度，以模擬地基模型在剪切波作用下的剪切動力響應。

等位移邊界下的地震動輸入可以采用波動輸入法，即在兩側(cè)邊設置等位移邊界，底邊設置無限元或粘彈性邊界來實現(xiàn)，也可以直接采用振動輸入法。結(jié)合上述粘彈性邊界、無限元邊界的Python 程序，同樣開發(fā)了等位移邊界及地震動輸入（圖1）。

1.4 遠置邊界下的地震動輸入

在結(jié)構(gòu)動力響應分析時間內(nèi)，由于遠置邊界距離結(jié)構(gòu)足夠遠，模型邊界反射波不會對結(jié)構(gòu)響應產(chǎn)生影響，保證了動力響應分析的可靠性。模型邊界寬度根據(jù)下式確定［15］：

式中：L為模型邊界距離結(jié)構(gòu)所在位置的距離；T為結(jié)構(gòu)動力響應分析時間；C為地震波在土體介質(zhì)中的傳播速度。

根據(jù)式（3）可知，當?shù)卣鸩ǔ謺r較長、土體介質(zhì)波速較大時，需要巨大的模型尺寸才可避免邊界對結(jié)構(gòu)位置的影響，并不具備實用性。遠置邊界在動力分析時間內(nèi)是準確的，常被作為波動問題的解析解，以驗證其他邊界的有效性。

2 波動輸入法人工邊界有效性驗證

本節(jié)從自由場模型、土層－隧道相互作用模型2 個方面，分別建立二維地震響應計算模型，對無限元邊界、粘彈性邊界、自由度綁定邊界、豎向約束邊界及相應的地震動輸入方式進行系統(tǒng)性驗證。針對波動輸入方式下的模型驗證時，各人工邊界的選擇見表1。

表1 波動輸入方式下的邊界選擇Table 1 Boundary selection under wave input methods

為方便對比各邊界條件下的計算結(jié)果，從模型底部垂直輸入SV波，其表達式如下［16］：

2.1 自由場模型驗證

考慮不同人工邊界的二維土體自由場模型如圖4所示。模型尺寸取200 m×50 m，單元網(wǎng)格尺寸為1 m×1 m。土體單元采用線彈性的4節(jié)點平面應變單元（CPE4R）模擬，計算參數(shù)見表2。計算中，以遠置邊界自由場模型計算結(jié)果及理論值來驗證人工邊界模型的有效性。為了滿足遠置邊界條件，根據(jù)土體參數(shù)，由式（3）計算的遠置邊界自由場模型寬度至少為589 m，故取該邊界下模型寬度為1 000 m。按上述參數(shù)建立有限元模型后，采用文中開發(fā)的動力計算輔助程序施加各人工邊界及輸入地震動，開展自由場模型動力響應分析。

表2 土體物理力學參數(shù)Table 2 Soil physical mechanics parameters

圖4 不同人工邊界的自由場模型Fig.4 Free-field model considering different artificial boundaries

為分析動力響應結(jié)果，在模型底部及頂部中點設置A、B兩個監(jiān)測點。圖5 為各邊界下監(jiān)測點的水平位移時程結(jié)果。圖5中的理論解顯示，不考慮土體介質(zhì)阻尼的情況下，垂直入射的SV波經(jīng)過0.51 s后到達自由表面和反射波疊加，使自由表面監(jiān)測點A的位移振幅相比入射波放大2倍，而波形保持不變。反射波經(jīng)0.51 s后傳至模型底部，引起模型底部振動，之后被底部人工邊界吸收而不再引發(fā)有限域的振動。各人工邊界條件下，模型監(jiān)測點處的水平位移時程曲線均與理論曲線及遠置邊界模型結(jié)果曲線吻合較好。其中，粘彈性邊界在波動結(jié)束后有一個小幅震蕩，與文獻［10］結(jié)果類似，在實際計算中對結(jié)果產(chǎn)生的誤差可忽略不計?？梢?，上述人工邊界均能較好地吸收地表反射波，具有良好的計算精度，也說明所開發(fā)的人工邊界及地震動輸入程序是有效性的。

圖5 各人工邊界模型監(jiān)測點A、B處的水平位移時程曲線（波動輸入）Fig.5 Time-history curve of horizontal displacement at monitoring points A and B of each artificial boundary models（wave input）

2.2 隧道地震響應模型驗證

為驗證各人工邊界對隧道地震響應計算的有效性，在自由場模型基礎上，進一步建立考慮土層－隧道動力相互作用的橫向地震響應模型見圖6。其中，將遠置邊界模型結(jié)果作為解析解。模型尺寸和土層條件同自由場模型，隧道埋深取為10 m，隧道外半徑取為3 m、襯砌厚度為0.35 m。襯砌采用線彈性梁單元（B21）模擬，計算參數(shù)見表3。

表3 襯砌物理力學參數(shù)Table 3 Lining physical mechanical parameters

圖6 不同人工邊界的自由場模型Fig.6 Free-field model considering different artificial boundaries

同時，引入圓形隧道在橫向剪切地震作用下的4種擬靜力解析解作為襯砌內(nèi)力數(shù)值解的對比（表4）［17－20］。剪切波作用下的圓形隧道解析解，通常是基于土－隧道系統(tǒng)周圍的遠場地層處于純剪狀態(tài)的假設［12］，將剪應力等效為沿45o方向產(chǎn)生的拉、壓主應力，此時襯砌發(fā)生橢圓化的剪切變形，將土體的變形施加到襯砌結(jié)構(gòu)上，即可確定襯砌結(jié)構(gòu)在橫向地震動作用下的內(nèi)力。其中，土體變形由土體所受剪應力的大小以及土－隧道結(jié)構(gòu)的相對剛度來決定。對于半徑為R 的圓形隧道，隧道周邊的遠場土層受到的剪應力通過式（5）計算得到，而土－隧道結(jié)構(gòu)的相對剛度則采用式（6）中的壓縮系數(shù)C和柔度系數(shù)F量化表示［20］：

表4 橫向地震作用下圓形隧道結(jié)構(gòu)內(nèi)力解析解Table 4 Analytical solution of internal force of circular tunnel structure under transverse earthquake

式中：R為隧道外半徑；Es、νs分別為土體彈模、泊松比；γc為隧道周圍2R范圍內(nèi)土體平均剪切應變峰值；τ為周圍土層在剪切地震作用下產(chǎn)生的剪應力；C、F分別為壓縮系數(shù)和柔度系數(shù)；El、νl、Al、Il分別為襯砌的彈模、泊松比、截面積和慣性矩。

需要指出的是，隧道結(jié)構(gòu)與周圍土體介質(zhì)的接觸面是比較粗糙的，不易出現(xiàn)完全滑移的極端條件［12］，因此，本文驗證計算僅考慮土與隧道結(jié)構(gòu)無相對滑移的情況。

圖7 為各邊界條件下襯砌拱肩（θ= 45o）處的內(nèi)力時程曲線。從圖7 可看出，各人工邊界模型的計算結(jié)果在0～1.5 s 之間與遠置邊界模型的結(jié)果基本吻合，其中，粘彈性邊界的計算內(nèi)力結(jié)果在1.5～2 s 之間有較小波動；自由度綁定邊界和豎向約束邊界的計算內(nèi)力結(jié)果在2.5～3 s之間有較小波動?？傮w而言，各人工邊界計算結(jié)果的整體誤差較小，在實際分析中可忽略不計。

圖7 不同人工邊界典型截面內(nèi)力時程曲線（波動輸入）Fig.7 Time-history curve of internal force of typical section with different artificial boundaries（wave input）

結(jié)合土－隧道相互作用模型的參數(shù)，通過式（4）可計算得到土－隧道系統(tǒng)的壓縮系數(shù)C和柔度系數(shù)F分別為0.002 5、1.476。其次，利用自由場模型進行地震響應分析，得到隧道所處地層范圍內(nèi)的平均剪切應變峰值γc為0.002 18。據(jù)此，可計算出地震作用下隧道所處地層平均應變峰值時刻襯砌截面內(nèi)力的解析解，并與動力時程分析中的峰值內(nèi)力結(jié)果對比，如圖8所示。

由圖8中可知，Wang解、Bobet解和Park 解的全截面動軸力結(jié)果比較接近，而Penzien 解嚴重低估了截面動軸力。3種典型人工邊界下的襯砌截面動軸力分布完全一致。與3種典型人工邊界下的動軸力相比，解析解在45°及135°截面及其附近約45°范圍內(nèi)偏小，最大相差8.5%，其它截面位置的動軸力分布與人工邊界下的基本一致。

圖8 波動輸入的峰值時刻襯砌截面動內(nèi)力分布（各人工邊界對比）Fig.8 Distribution of dynamic internal force of lining section at peak time of wave input（comparison of artificial boundaries）

襯砌截面動彎矩的4 種解析解略有差異，其中，在±45°軸線的截面上差別最明顯，彎矩值最大相差約16.02%，Wang解最大，Penzien解次之，Bobet及Park解最小且完全一致。各人工邊界下的襯砌截面動彎矩分布完全一致，不同邊界下的動彎矩在0°～180°范圍內(nèi)與Penzien解基本一致，而在180°～360°范圍內(nèi)則與Bobet解及Park解吻合。對比各邊界下的數(shù)值解與解析解可知，各人工邊界下的內(nèi)力分布及內(nèi)力值均與解析解接近，表明本文開發(fā)的人工邊界及地震動施加程序可用于隧道地震響應分析。對比上述邊界還可看出，相較于無限元邊界和粘彈性邊界，等位移邊界有相當?shù)挠嬎憔?，而且邊界施加簡單，更實用便捷?/p>

3 振動輸入法人工邊界有效性驗證

3.1 自由場模型驗證

同樣采用式（4）所示的剪切波進行驗證。考慮到振動輸入會使地震波在有限域中來回反射，為了反映這種特點，將動力分析時長調(diào)整為6 s。并以遠置邊界計算結(jié)果為解析解，按式（3）計算得到遠置邊界模型寬度需大于1 176.7 m，因此，取自由場模型寬度為2 000 m。結(jié)合上述分析結(jié)果可知，等位移邊界（如自由度綁定邊界）與其他人工邊界相比，施加更簡便。為節(jié)約篇幅，本節(jié)采用自由度綁定邊界開展振動輸入法的有效性驗證，其自由場模型和隧道地震響應模型的尺寸及參數(shù)設置分別沿用2.1節(jié)和2.2節(jié)。

自由場模型計算結(jié)果如圖9 所示，與波動輸入法不同，振動輸入法的模型底部由輸入的位移時程控制，不因自由表面反射波的到來而改變，但下行反射波會經(jīng)過底部邊界再次向自由表面反射，引起模型自由表面位移出現(xiàn)反復的波動。

圖9 不同人工邊界模型監(jiān)測點A、B處的水平位移時程曲線（振動輸入）Fig.9 Time-history curves of horizontal displacement（vibration input）at monitoring points A and B of different artificial boundary models

3.2 隧道地震響應模型驗證

圖10為2種邊界條件下的襯砌拱肩（θ=45o）部位內(nèi)力時程曲線。與波動輸入法明顯不同，由于地震波在模型中的來回反射，土層－隧道體系受到了循環(huán)的地震動荷載作用，導致隧道結(jié)構(gòu)內(nèi)力也出現(xiàn)了循環(huán)響應特征。

圖10 不同人工邊界下的襯砌典型截面內(nèi)力時程曲線（振動輸入）Fig.10 Time-history curve of internal force of typical lining section under different artificial boundaries（vibration input）

圖11為解析解與振動輸入法數(shù)值解的對比結(jié)果。從圖中可以看出，對于襯砌截面動軸力，除Penzien解以外的3種解析解與數(shù)值解都比較接近。對于襯砌截面動彎矩，數(shù)值解與解析解均相差較小。由此說明，自由度綁定邊界對振動輸入法是有效性的。

圖11 峰值時刻襯砌截面內(nèi)力分布曲線（振動輸入）Fig.11 The internal force distribution curve of the lining section at the peak moment（vibration input）

4 波動輸入與振動輸入差異性

在地下結(jié)構(gòu)抗震分析中，由于輸入的地震波持時較長，導致結(jié)構(gòu)的地震響應可能由多個波場共同決定。因此，為探討這種多波場條件下振動輸入與波動輸入的區(qū)別，將土體彈性模量調(diào)整為1 GPa，其余模型參數(shù)不變，則土體介質(zhì)剪切波速cs=438.5 m/s，由此可計算出入射波從底邊界入射后，僅需0.228 s便可經(jīng)過自由表面反射而傳回底部，此時原先的入射波場振動還未結(jié)束，隨即產(chǎn)生了多波場疊加效應。本節(jié)同樣以自由度綁定邊界為例，如圖12（a）所示，波動輸入模型底部A點在0.228～0.5 s間的位移便是由入射波與反射波疊加形成的。而振動輸入模型底部受位移時程控制，不因自由表面反射波的到來而改變。如圖12（b）所示，振動輸入的模型底部反射波還未再次到達自由表面時（0.342 s之前），自由表面僅發(fā)生一次反射，其位移較入射位移放大2倍；當自由表面再次發(fā)生反射（0.342 s之后），此時多波場疊加，位移場放大倍數(shù)超過2倍，土體地震響應加劇。由于底部人工邊界對反射波的具有吸收功能，采用波動輸入的模型頂部B點僅發(fā)生一次反射作用，位移場放大2倍。

圖12 2種地震動輸入下監(jiān)測點A、B的位移時程曲線Fig.12 Displacement time-history curves of monitoring points A and B under two types of ground motion inputs

如圖13、圖14所示，由于地層中多波場疊加效應使得位移場增強，振動輸入時的隧道襯砌的峰值內(nèi)力大于波動輸入模型的計算結(jié)果?？梢?，振動輸入與波動輸入對結(jié)構(gòu)的地震響應有明顯的影響，振動輸入會使地層結(jié)構(gòu)產(chǎn)生更劇烈的循環(huán)加卸載行為，會放大土層及結(jié)構(gòu)的動力響應。因此，在進行地下結(jié)構(gòu)地震響應計算時應根據(jù)具體的場地條件選擇合理的地震動輸入方法，當下臥基巖滿足剛性條件假定可考慮使用振動輸入法，而絕大多數(shù)地下結(jié)構(gòu)的場地條件宜采用波動輸入法。

圖13 2種地震動輸入下襯砌典型截面內(nèi)力時程曲線Fig.13 Time-history curve of internal force of typical lining section under two types of ground motion inputs

圖14 2種地震動輸入下襯砌截面的動內(nèi)力分布（峰值時刻）Fig.14 Dynamic internal force distribution of lining section under two types of ground motion inputs（peak moment）

5 結(jié)論

（1）基于Python 語言，編制了無限元邊界、粘彈性邊界、等位移邊界等常見動力人工邊界施加及其地震動輸入程序。進一步結(jié)合自由場模型及隧道地震響應模型，驗證了波動輸入法和振動輸入法及與之相適應的人工邊界模型的有效性。此基礎上，開發(fā)了隧道地震響應計算輔助程序，大幅度減少了地震響應計算前處理的工作量。

（2）采用等效節(jié)點力的波動輸入時，在模型底邊及兩側(cè)邊設置無限元或粘彈性邊界，以及在兩側(cè)邊設置等位移邊界、底邊設置無限元或粘彈性邊界，計算結(jié)果均與解析解基本重合，具有良好的計算精度，可用于隧道結(jié)構(gòu)地震響應的計算分析。而采用位移、速度或加速度時程的振動輸入時，只適用于遠置邊界或等位移邊界。

（3）自由度綁定邊界簡化了波動輸入法實現(xiàn)的繁瑣程度，且具有與其他人工邊界相當?shù)挠嬎憔?，在計算中更為實用便捷?/p>

（4）與波動輸入相比，振動輸入使得有限域成為半封閉系統(tǒng)，入射波會在底部重復反射，當入射波持時較長或土層波速較大時，土層中發(fā)生多波場疊加，會放大土層及結(jié)構(gòu)的動力響應。因此，下臥基巖滿足剛性條件假定可考慮使用振動輸入法，而絕大多數(shù)地下結(jié)構(gòu)的場地條件宜采用波動輸入法。

（5）文中研究成果適用于二維地震響應分析，擴展到三維條件的人工邊界模型及地震動輸入方法還需進一步完善。