陽 棟，王志亮，2

（1.同濟大學土木工程學院，上海200092；2.同濟大學巖土與地下工程教育部重點實驗室，上海200092）

為了對無限地基進行模擬，必須截取有限范圍的地基，然后設置合理的無反射人工邊界模擬無限地基輻射阻尼效應.局部人工邊界由于其時空解耦性、計算耗時少，在地基輻射阻尼效應模擬中獲得廣泛應用.Lysmer等［1］基于一維波動理論提出黏性邊界，在許多波動問題中得到了廣泛應用，也是我國核電站抗震設計規(guī)范建議使用的人工邊界之一［2］.廖振鵬［3］基于透射理論提出了透射邊界，但是計算中可能引起高頻失穩(wěn)，而且編程比較復雜.Deeks等［4］、劉晶波等［2］基于柱面波動方程建立了可以考慮半無限介質彈性恢復能力的二維黏彈性人工邊界，能模擬人工邊界外無限介質彈性恢復性能，有良好的頻率穩(wěn)定性，但是黏彈性邊界不方便考慮初始重力作用等.

無限元概念最早由Ungless［5］提出，后經(jīng)過Bettess［6］，Beer等［7］和Zienkiewicz［8］的改進和發(fā)展，國內學者ZHANG等［9］、葛修潤等［10］也對無限元的研究開展了不少工作.燕柳斌［11］利用無限元法模擬半空間彈性地基和重力壩地基；黃勝等［12］用無限元對隧道地震響應進行模擬，但是沒有考慮自由波場的作用.相對于黏性邊界，無限元不僅能模擬遠場地基對地震波能量的吸收作用，還能夠準確模擬無窮遠處位移為零的邊界條件，采用無限元能夠大大減少單元數(shù)量，節(jié)省計算費用［13］.

本文采用有限元與無限元耦合法模擬地基動力響應，為了考慮人工截斷邊界對波傳播影響，應用波場分離技術，在有限元和無限元分界處施加等效節(jié)點荷載，以實現(xiàn)地震輸入，并對邊界和地震輸入的有效性進行驗證.同時，把該法應用于核電站地震響應計算，對自由場地和考慮土與結構相互作用（SSI）時地震響應進行對比分析.

1 無限元原理與控制方程

ABAQUS中提供了一階和二階無限單元，其在連續(xù)實體靜力分析中能提供剛度，在動態(tài)分析中提供靜止邊界［14］.在進行數(shù)值模擬時，由于主要分析近場區(qū)域結果，故需用有限元模擬，而遠場地基可假定為線性的，用無限元模擬.

假定波沿z軸傳播，則一維波動的平衡方程

對于均勻彈性介質平面波

波動方程uz＝f1（z－ct）代表沿z軸正向傳播的波，而uz＝f2（z＋ct）代表沿z軸負向傳播的波，假設有限元與無限元的分界面為z＝L，因為問題是線性的，所以位移可以疊加，即

當波傳播到界面上時，界面上的應力為

式中：f′1，f′2分別為函數(shù)f1，f2對z的導數(shù).

如果在有限元邊界上設置阻尼器平衡波動力，其阻尼系數(shù)為C，并使下式成立：

為了達到不反射的效果，則反射波引起的位移和應力為零，即f2＝0，f′2＝0，由此解得

通過設置參數(shù)，就可以避免有限元與無限元邊界處波的能量反射.

2 波場分離方法和地震輸入

對于內源波動問題，如動力機器、地鐵等引起的振動問題，用無限元模擬遠場時，可以直接在近場有限元區(qū)域內施加振動作用，當波傳播到有限元邊界時，自動透射進入無限元區(qū)域.而對于外源波動問題，如地震作用，存在地震波散射現(xiàn)象，需要在人工邊界處進行波場分解.根據(jù)不同人工邊界上的波場特點，各邊界采用不同的分解方案，即：底邊界上將總場分解為邊界入射場和邊界外行場，側邊界上將總場分解為自由場和散射場，邊界入射場或自由場可由連續(xù)介質力學模型解析計算得到，邊界外行場或散射場由人工邊界條件模擬并由離散模型依據(jù)數(shù)值分析方法獲得.這種不同邊界分區(qū)分解的方法，在盡量避免計算自由場的前提下較好地解決了人工邊界的平行波模擬問題［15］.

外源作用下，有限元邊界任意節(jié)點運動方程為

式中：uli，，分別為l節(jié)點在i方向的位移、速度和加速度；Cl為節(jié)點l的阻尼系數(shù)；Kl為節(jié)點l的剛度系數(shù)；和分別為模擬邊界自由場和散射場在l節(jié)點沿i方向施加的等效荷載；ml為節(jié)點l質量.

在模擬散射場時，有限元邊界節(jié)點上施加的應力由無限元提供，無需在邊界節(jié)點上附加額外的荷載；模擬自由場時，在有限元邊界節(jié)點上需要施加包括克服無限單元阻尼所需的應力和自由場在有限元邊界處的應力.利用式（7）進行波動輸入的前提是首先確定自由場，而自由場一般采用解析法或數(shù)值法確定.杜修力等［16］給出了避免計算自由場前提下的黏彈性邊界等效荷載計算方法，在此基礎上去掉彈簧力作用，即可得到無限元邊界等效荷載的計算公式.下面以平面問題S波從有限元底部邊界垂直入射為例，其等效荷載計算式為在底邊界上

在側邊界處

式中：Δt1，Δt2分別為l節(jié)點處入射S波和地表反射S波的時間延遲，即Δt1＝h／cS，Δt2＝（2 H－h(huán)）／cS，其中H為地表面到有限元區(qū)域底部邊界的距離，h為l節(jié)點到有限元底邊界的距離，cS為S波波速；Al為有限元邊界節(jié)點l的影響面積；ClN，ClT為法向和切向阻尼系數(shù)，其值為ClN＝ρcS，ClT＝ρcP，cP為縱波波速.等效地震荷載下標代表節(jié)點號和分量方向，上標代表節(jié)點所在有限元邊界界面的外法線方向，與坐標軸方向一致為正.

在ABAQUS計算時，先用FORTRAN90編制程序計算等效荷載，然后選取邊界節(jié)點，通過調用子程序在邊界節(jié)點上施加等效荷載.

3 邊界及地震輸入的可行性驗證

為了驗證無限元和等效荷載輸入模擬自由場地震的有效性，本文通過一個簡單的模型計算，驗證半無限地基中波動傳播規(guī)律.對于二維彈性半空間模型，從底部垂直入射一S波，入射波方程如式（13）所示，求二維彈性半平面內各點的位移反應

已知巖石地基的剪切模量G＝4.233 6GPa，泊松比ν＝0.25，密度ρ＝2 100kg·m－3，波速cS＝1 400m·s－1，cP＝2 425m·s－1.選取有限元的計算區(qū)域寬為762m，高381m，用平面應變四邊形減縮積分單元剖分網(wǎng)格，單元尺寸為Δx＝Δy＝19.05 m，有限單元數(shù)為800個，節(jié)點數(shù)為861個.對于無限元區(qū)域，選取頂部中點（C點）作為輻射中心，延長C點與有限元邊界節(jié)點的連線，并使無限單元的長度等于有限元邊界節(jié)點到輻射中心的長度，建立無限元網(wǎng)格.總的計算時間為2.0s，取固定增量步長為0.005s，計算網(wǎng)格模型如圖1a所示.

圖2 A，B和C點水平位移時程曲線Fig.2 Time histories of horizontal displacement at nodes A，Band C

由一維波動理論可知，波分別在t＝0.136s和t＝0.272s時到達模型中部（B點）和頂部（C點）.由于頂部是自由表面，故頂部的最大振幅相比入射波放大2倍，波在頂部經(jīng)自由表面反射后向下傳播，在t＝0.408s時反射波到達中部，在t＝0.545s時反射波到達底部（A點），在t＝0.945s時波完全經(jīng)過底部有限元邊界并透射進入無限元區(qū)域，此時有限元區(qū)域中所有節(jié)點的位移均變?yōu)榱?從圖2可看出，模擬結果與理論值吻合較好，這也說明上述采用無限元結合等效荷載輸入來模擬波在半無限地基中傳播的方法是正確的.

4 工程應用實例

4.1 計算模型和地震輸入

設有一核電站建在上述場地上（見圖1b），把核電站簡化為一矩形結構，尺寸114.3m×38.1m，單元尺寸取為19.05m.混凝土結構密度為2 643 kg·m－3，彈性模量為31.027GPa，泊松比為0.15.在ABAQUS中采用隱式積分方法求解.

本例以Koyna地震的S波［14］作為輸入波，對自由場及考慮SSI效應下的地震響應進行分析.由于實際記錄的是地表地震波，其中已包含地基的放大效應，在應用到底部基巖時，可取其幅值的一半.通常記錄的地震波數(shù)據(jù)中混雜了各種干擾，直接對加速度時程進行二次積分得到的位移時程會產(chǎn)生漂移失真，因此要對加速度時程進行校準和濾波，達到消除噪聲的目的.此處用地震分析軟件Seismosignal對原始的加速度時程數(shù)據(jù)進行處理，圖3是經(jīng)校準和濾波后的加速度時程曲線、位移時程曲線及加速度反應譜.

圖3 輸入波的加速度時程、位移時程曲線及加速度反應譜Fig.3 Acceleration and displacement time histories and acceleration response spectrum of input wave

4.2 自由場分析

圖4 無限元自由場地震加速度反應譜Fig.4 Acceleration response spectrum of free－field modeled with infinite element

此處以圖1a的自由場為例，分別進行施加重力與不施加重力的地震模擬，前者在地震輸入前單獨建立一個時間步長為1s的重力加載步，總的分析時間為12.0s.圖4a，b分別為無重力與考慮重力作用下，自由場中A，B，C三點的加速度反應譜，圖5a為自由場頂部C點的位移時程曲線.從圖4和5a可看出，重力對地震加速度和位移響應都無影響，這是因為重力不影響水平方向力的平衡.在周期為0.25s和0.70s時，場地頂部的加速度反應譜放大超過2倍，而周期為1.00s時，場地對加速度反應譜的放大作用最不明顯，由于I類核電場地的特征周期約為0.25s，故0.25s處的加速度反應譜值最大.圖4a中A點曲線為有限元底部地震波入射點處的加速度反應譜，由于其中包括反射波的影響，故其與圖3c中入射波的加速度反應譜有所區(qū)別.C點位移是入射波和反射波共同作用的結果，對比圖5a和圖3b，可以看出自由場位移響應值為入射位移波的2倍.通過與一維地震分析程序SHAKE91［17］的計算結果作對比（如圖4a、圖5a所示），也可以看出本文采用的邊界和地震動輸入是有效的.

4.3 SSI效應分析

對圖1b中模型輸入上述Koyna波進行計算，并選取節(jié)點A1～G1對結果進行分析.圖6為考慮SSI效應后的地震加速度反應譜，對比圖6a和4a可以看出，考慮SSI后，結構底部處C點的加速度反應譜最大值從1.92g下降到0.94g（其中g為重力加速度），且在周期為0.25s附近，加速度反應譜下降，而在周期為0.30s到4.00s范圍內，加速度反應譜增大，即SSI效應使短周期段的加速度反應譜值降低，而長周期段的加速度反應譜值增加.對于結構底部標高的地基節(jié)點C1，E1，F(xiàn)1，G1，從圖6b可看出，在結構與有限元邊界的中間點F1處，加速度反應譜最大值為2.36g，大于自由場中的加速度反應譜最大值1.92g，即在S波的激勵下，SSI效應使結構附近地基的加速度反應譜降低，而遠離結構處的加速度反應譜峰值相比自由場反而略有升高.圖5b為自由場頂部C點和考慮SSI效應后相同位置處C1點的位移時程曲線圖，可以看出S波作用下，SSI效應對地基的水平位移無影響.圖7為地表C點在自由場和考慮SSI效應時最大主應力和最小主應力的時程曲線，由于結構底部和地基之間將產(chǎn)生較大的摩擦力，從圖中容易發(fā)現(xiàn)SSI效應使得地基中C點的最大主應力和最小主應力的幅值均增大2倍之多.

5 結論

（1）采用近場地基由有限元模擬，遠場土體由無限元模擬，并使用本文所提出的等效荷載計算公式計算等效節(jié)點力，在有限元和無限元交界處輸入地震作用，此法能夠有效仿真地震波在半無限地基中的傳播過程.

圖7 自由場C點及考慮SSI效應時C1點的主應力時程曲線Fig.7 Time histories of principal stress at point Cin free－field and point C1with a consideration of SSI effect

（2）對于巖石地基，當?shù)卣鸩◤牡撞炕鶐r傳播到地表時，加速度反應譜峰值和位移振幅都放大2倍，通過自由場地施加重力和無重力的對比分析，可以看出S波作用下，重力對地震作用時地基的加速度反應譜和水平位移響應無明顯影響.

（3）在S波作用下，SSI效應使短周期段的加速度反應譜值降低，而長周期段的加速度反應譜值增加；結構附近地基的加速度反應譜峰值相比自由場降低，而遠離結構處加速度反應譜峰值相比自由場略有提高；SSI作用對地基的水平位移無影響，但使結構底部地基的最大主應力和最小主應力幅值都有增大趨勢.

［1］ Lysmer J，Kulemeyer R L.Finite dynamic model for infinite media［J］.Journal of Engineering Mechanics Division，ASCE，1969，95（4）：759.

［2］ 劉晶波，谷音，杜義欣.一致粘彈性人工邊界及粘彈性邊界單元［J］.巖土工程學報，2006，28（9）：1070.

LIU Jingbo，GU Yin，DU Yixin.Consistent viscous－spring artificial boundaries and viscous－spring boundary elements［J］.Journal of Geotechnical Engineering，2006，28（9）：1070.

［3］ 廖振鵬.工程波動理論導論［M］.北京：科學出版社，2002.

LIAO Zhenpeng.Introduction to wave motion theories in engineering［M］.Beijing：Science Press，2002.

［4］ Deeks A J，Randolph M F.Axisymmetric time－domain transmitting boundaries［J］.Journal of Engineering Mechanics Division，ASCE，1994，120（1）：25.

［5］ Ungless R F.An infinite finite element［D］.Vancouver：University of British Columbia，1973.

［6］ Bettess P.More on infinite element［J］.International Journal for Numerical Methods in Engineering，1980，15（11）：1613.

［7］ Beer G，Meek J L.Infinite domain element［J］.International Journal for Numerical Methods in Engineering，1981，17（3）：43.

［8］ Zienkiewicz O C.A novel boundary infinite element［J］.International Journal for Numerical Methods in Engineering，1983，19（3）：393.

［9］ ZHANG Chuhan，ZHAO Chongbin.Coupling method of finite and infinite elements for strip foundation wave problems［J］.Earthquake Engineering and Structural Dynamics，1987，15（7）：839.

［10］ 葛修潤，谷先榮，豐定祥.三維無限元和節(jié)理無界元［J］.巖土工程學報，1986，8（5）：9.

GE Xiurun，GU Xianrong，F(xiàn)ENG Dingxiang.Threedimensional infinite domain elements and joint infinite domain elements［J］.Journal of Geotechnical Engineering，1986，8（5）：9.

［11］ 燕柳斌.用三維映射無限元模擬重力壩地基［J］.水利學報，1991，2（10）：7.

YAN Liubin.Modeling of gravity dam foundation using three dimensional mapped infinite element［J］.Journal of Hydraulic Engineering，1991，2（10）：7.

［12］ 黃勝，陳衛(wèi)忠，楊建平，等.地下工程地震動力響應及抗震研究［J］.巖石力學與工程學報，2009，28（3）：483.

HUANG Sheng，CHEN Weizhong，YANG Jianping，et al.Research on earthquake－induced dynamic response and aseismic measures for underground engineering［J］.Chinese Journal of Rock Mechanics and Engineering，2009，28（3）：483.

［13］ 姜袁，陳燈紅，姚艷華.有限元－無限元法在重力壩應力分析中的應用［J］.武漢大學學報：工學版，2009，42（3）：322.

JIANG Yuan，CHEN Denghong，YAO Yanhua.Application of finite element and infinite element method to stress analysis of gravity dam［J］.Engineering Journal of Wuhan University，2009，42（3）：322.

［14］ SIMULIA Inc.ABAQUS theory manual［EB／OL］.［2011－04－20］.USA：https：／／www.sharcnet.ca／Software／Abaqus610／Documentation／ddoc／v6.10／books／stm／default.htm.

［15］ 趙建鋒，杜修力，韓強，等.外源波動問題數(shù)值模擬的一種實現(xiàn)方式［J］.工程力學，2007，24（4）：52.

ZHAO Jianfeng，DU Xiu1i，HAN Qiang，et al.An approach to numerical simulation for external source wave motion［J］.Engineering Mechanics，2007，24（4）：52.

［16］ 杜修力，趙密.基于黏彈性邊界的拱壩地震反應分析方法［J］.水利學報，2006，37（9）：1063.

DU Xiuli，ZHAO Mi.Analysis method for seismic response of arch dams in time domain based on viscous－spring artificial boundary condition［J］.Journal of Hydraulic Engineering，2006，37（9）：1063.

［17］ Idriss I M，Sun J I.SHAKE 91：a computer program for conducting equivalent linear seismic response analyses of horizontally layered soils deposits［EB／OL］.［2011－04－20］.USA：http：／／nees.ucdavis.edu／publications／Shak91＿manual.pdf.