陽 棟，王志亮，2

（1.同濟(jì)大學(xué)土木工程學(xué)院，上海200092；2.同濟(jì)大學(xué)巖土與地下工程教育部重點(diǎn)實(shí)驗(yàn)室，上海200092）

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

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

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

1 無限元原理與控制方程

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

假定波沿z軸傳播，則一維波動(dòng)的平衡方程

對(duì)于均勻彈性介質(zhì)平面波

波動(dòng)方程uz＝f1（z－ct）代表沿z軸正向傳播的波，而uz＝f2（z＋ct）代表沿z軸負(fù)向傳播的波，假設(shè)有限元與無限元的分界面為z＝L，因?yàn)閱栴}是線性的，所以位移可以疊加，即

當(dāng)波傳播到界面上時(shí)，界面上的應(yīng)力為

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

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

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

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

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

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

外源作用下，有限元邊界任意節(jié)點(diǎn)運(yùn)動(dòng)方程為

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

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

在側(cè)邊界處

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

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

3 邊界及地震輸入的可行性驗(yàn)證

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

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

圖2 A，B和C點(diǎn)水平位移時(shí)程曲線Fig.2 Time histories of horizontal displacement at nodes A，Band C

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

4 工程應(yīng)用實(shí)例

4.1 計(jì)算模型和地震輸入

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

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

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

4.2 自由場分析

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

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

4.3 SSI效應(yīng)分析

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

5 結(jié)論

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

圖7 自由場C點(diǎn)及考慮SSI效應(yīng)時(shí)C1點(diǎn)的主應(yīng)力時(shí)程曲線Fig.7 Time histories of principal stress at point Cin free－field and point C1with a consideration of SSI effect

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

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

［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］.巖土工程學(xué)報(bào)，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］ 廖振鵬.工程波動(dòng)理論導(dǎo)論［M］.北京：科學(xué)出版社，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］.巖土工程學(xué)報(bào)，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］.水利學(xué)報(bào)，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)忠，楊建平，等.地下工程地震動(dòng)力響應(yīng)及抗震研究［J］.巖石力學(xué)與工程學(xué)報(bào)，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］ 姜袁，陳燈紅，姚艷華.有限元－無限元法在重力壩應(yīng)力分析中的應(yīng)用［J］.武漢大學(xué)學(xué)報(bào)：工學(xué)版，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］ 趙建鋒，杜修力，韓強(qiáng)，等.外源波動(dòng)問題數(shù)值模擬的一種實(shí)現(xiàn)方式［J］.工程力學(xué)，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］ 杜修力，趙密.基于黏彈性邊界的拱壩地震反應(yīng)分析方法［J］.水利學(xué)報(bào)，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.