ABAQUS巖土動(dòng)力分析的靜-動(dòng)力統(tǒng)一人工邊界研究

左得奇，蔣良濰，駱 達(dá)，杜美玲，葛學(xué)軍

（1.西南交通大學(xué)土木工程學(xué)院，四川成都 610031；2.高速鐵路線路工程教育部重點(diǎn)實(shí)驗(yàn)室，四川成都 610031；3.中鐵二院工程集團(tuán)有限責(zé)任公司，四川成都 610031）

引言

數(shù)值分析是巖土靜、動(dòng)力問(wèn)題不可或缺的研究手段之一，人工邊界的設(shè)置及動(dòng)力作用輸入方法將直接影響數(shù)值計(jì)算結(jié)果的合理性。因此，各類人工邊界和地震動(dòng)輸入方法相繼被提出。人工邊界可分為全局人工邊界和局部人工邊界，Lysmer等［1］首先提出了以黏性理論為基礎(chǔ)的局部人工邊界，其概念清楚、應(yīng)用方便，但精度較低，容易出現(xiàn)低頻失穩(wěn)問(wèn)題；Deeks等［2］進(jìn)一步發(fā)展提出了黏彈性人工邊界，其引入彈簧單元以模擬地基恢復(fù)力，克服了黏性邊界低頻失穩(wěn)的問(wèn)題；劉晶波等［3］推導(dǎo)了三維形式的黏彈性人工邊界并將其應(yīng)用于地下土-結(jié)構(gòu)的相互作用的計(jì)算中；廖振鵬等［4］提出了一種對(duì)一般無(wú)限域模型具有普遍適應(yīng)性的人工透射邊界，入射波由一維單側(cè)波動(dòng)的疊加構(gòu)成，其具有實(shí)用性廣、與有限元方法結(jié)合方便的優(yōu)點(diǎn)，但也存在低頻失穩(wěn)問(wèn)題；在全局人工邊界方面，Ungless［5］首先提出了無(wú)限元思想；Kim［6］和Yun［7］利用無(wú)限元方法在頻域和時(shí)域內(nèi)研究了二維、三維層狀土-結(jié)構(gòu)的相互作用，其具有使用方便、前處理工作量小的優(yōu)勢(shì)；此外，動(dòng)力作用輸入方式與人工邊界條件的設(shè)置密切相關(guān)，劉晶波等［8］提出了基于黏彈性人工邊界的地震動(dòng)輸入方法；杜修力等［9］考慮場(chǎng)地自由場(chǎng)效應(yīng)推導(dǎo)了下臥剛性基巖條件下的地震動(dòng)輸入方法；馬笙杰等［10］則對(duì)黏彈性人工邊界設(shè)置及地震動(dòng)輸入方法進(jìn)行了詳細(xì)對(duì)比分析；趙武勝等［11］通過(guò)分析人工邊界相互作用推導(dǎo)了地震波垂直及斜入射條件下地震動(dòng)輸入方法，并以無(wú)限元為例進(jìn)行了正確性驗(yàn)證，上述地震動(dòng)輸入方法多為基于動(dòng)力人工邊界推導(dǎo)得出，適用于動(dòng)靜分算的動(dòng)力分量計(jì)算。

對(duì)于線彈性、小變形問(wèn)題，可分別利用靜、動(dòng)力邊界計(jì)算靜、動(dòng)力響應(yīng)分量，再疊加得到完整解。由于土體是一種非線性材料，在力學(xué)作用下可能發(fā)生非線性響應(yīng)甚至塑性變形，疊加原理將不再適用，為此研究者提出了靜-動(dòng)力統(tǒng)一人工邊界。劉晶波等［12］基于彈性半無(wú)限空間中靜力問(wèn)題基本解對(duì)黏彈性人工邊界進(jìn)行節(jié)點(diǎn)約束剛度修正，提出了黏彈性靜-動(dòng)力統(tǒng)一人工邊界，使得靜、動(dòng)力分析能夠連續(xù)的進(jìn)行；高峰等［13］通過(guò)對(duì)比研究發(fā)現(xiàn)，該靜-動(dòng)力統(tǒng)一人工邊界應(yīng)用于初始地應(yīng)力生成誤差較大，并提出了一種基于有限元靜力邊界節(jié)點(diǎn)支反力的黏彈性人工邊界設(shè)置方法。

由上可知，目前人工邊界的研究主要以黏彈性人工邊界為主，因其在低頻問(wèn)題中計(jì)算穩(wěn)定性好、效率較高，通過(guò)邊界條件轉(zhuǎn)化［13］、引入邊界等效節(jié)點(diǎn)力［9-11］等方法，使得其能廣泛應(yīng)用于各類巖土動(dòng)力問(wèn)題，而其直接應(yīng)用于靜力計(jì)算時(shí)精度較低、邊界設(shè)置較為繁瑣，前處理工作量較大。無(wú)限元在理論上滿足在無(wú)窮遠(yuǎn)處位移為零，波傳至無(wú)窮遠(yuǎn)處衰減為零的客觀性條件，部分耦合了靜、動(dòng)力計(jì)算，且邊界設(shè)置較為簡(jiǎn)單、應(yīng)用方便，大量研究將其直接應(yīng)用于靜力計(jì)算（如初始應(yīng)力狀態(tài)生成［14-16］）和動(dòng)力計(jì)算（內(nèi)、外源振動(dòng)問(wèn)題［14-16］），王飛等［17］基于ABAQUS無(wú)限元的靜力計(jì)算結(jié)果具有可信性的前提，將文獻(xiàn)［11］中的地震動(dòng)輸入方法引入無(wú)限元人工邊界。但實(shí)際上，巖土動(dòng)力計(jì)算中，第一步為進(jìn)行初始應(yīng)力狀態(tài)生成即地應(yīng)力平衡，平衡后位移近似為零，而應(yīng)力狀態(tài)的合理性往往被忽略，無(wú)限元應(yīng)用于靜力步驟的計(jì)算精度會(huì)直接影響初始應(yīng)力狀態(tài)，由于在體力荷載作用下，靜力無(wú)限元是通過(guò)引入邊界節(jié)點(diǎn)力以模擬無(wú)限域?qū)τ邢抻虻募s束支承作用，但目前關(guān)于該節(jié)點(diǎn)力與靜力邊界是否等效的研究較少；在動(dòng)力計(jì)算中，無(wú)限元是一種黏性吸收邊界，其無(wú)法模擬地基的彈性恢復(fù)力及內(nèi)行波輸入，因此，有必要研究基于無(wú)限元進(jìn)行改進(jìn)的靜-動(dòng)力統(tǒng)一人工邊界。

基于此，簡(jiǎn)要介紹了ABAQUS靜、動(dòng)力無(wú)限元人工邊界的基本理論，結(jié)合二維水平地基模型算例，分析無(wú)限元直接應(yīng)用于初始地應(yīng)力計(jì)算單元節(jié)點(diǎn)支承剛度的合理性及計(jì)算精度，提出引入靜力等效節(jié)點(diǎn)力以提高無(wú)限元靜力及內(nèi)源計(jì)算精度的方法；為使得無(wú)限元人工邊界能應(yīng)用于需要考慮模型靜、動(dòng)力綜合效應(yīng)的內(nèi)、外源振動(dòng)問(wèn)題，推導(dǎo)內(nèi)、外源振動(dòng)問(wèn)題中需施加的邊界等效節(jié)點(diǎn)力，編制相關(guān)計(jì)算程序，提出基于無(wú)限元的靜-動(dòng)力統(tǒng)一人工邊界，并對(duì)該邊界的正確性進(jìn)行模型算例驗(yàn)證。

1 靜、動(dòng)力無(wú)限元人工邊界作用原理

根據(jù)振動(dòng)源與研究區(qū)域的位置關(guān)系，可將巖土動(dòng)力分析問(wèn)題大致分為內(nèi)源振動(dòng)和外源振動(dòng)，根據(jù)波場(chǎng)相對(duì)于模型邊界的傳播方向，可將模型總波場(chǎng)分解為內(nèi)行波場(chǎng)和外行波場(chǎng)［19］。內(nèi)源振動(dòng)問(wèn)題的振動(dòng)源位于研究區(qū)域內(nèi)部，如基礎(chǔ)振動(dòng)引起的土中動(dòng)應(yīng)力分析等，波源向外部的無(wú)限、半無(wú)限區(qū)域輻射能量，這個(gè)過(guò)程屬于外行波的散射問(wèn)題；外源振動(dòng)的振動(dòng)源位于研究區(qū)域之外，如地震引起地基和上部結(jié)構(gòu)的振動(dòng)分析等，該類問(wèn)題可分解為從外部無(wú)限域的內(nèi)行波輸入和研究區(qū)域內(nèi)部動(dòng)力響應(yīng)產(chǎn)生的外行波散射。由上可知，計(jì)算模型進(jìn)行人工邊界截?cái)鄷r(shí)，內(nèi)源振動(dòng)只需考慮外行波的吸收，而外源振動(dòng)條件下，外行波吸收和內(nèi)行波輸入均需考慮。

ABAQUS動(dòng)力分析中，將無(wú)限單元與有限單元結(jié)合使用，所關(guān)心的研究區(qū)域（近場(chǎng)有限域）以有限元模擬，而采用無(wú)限元模擬遠(yuǎn)場(chǎng)無(wú)限域。當(dāng)非體力荷載（如集中力等）作用時(shí)，靜力無(wú)限元邊界基于Zienkiewicz等［18］的靜力計(jì)算分析理論，設(shè)定距模型計(jì)算區(qū)域無(wú)窮遠(yuǎn)處位移為零，且假設(shè)從有限元-無(wú)限元交界面至無(wú)窮遠(yuǎn)處，位移呈線性分布。通過(guò)引入滿足上述假定條件的映射坐標(biāo)系以及具有二次收斂速度的位移插值函數(shù)，使得無(wú)限元成為能夠應(yīng)用于靜力計(jì)算的邊界［20］。當(dāng)體力荷載（如重力等）作用時(shí)，ABAQUS程序因無(wú)法在無(wú)限元中定義體力，以在有限元-無(wú)限元邊界自動(dòng)插入節(jié)點(diǎn)力的方式體現(xiàn)無(wú)限域?qū)τ邢抻虻募s束支承作用，但插入的節(jié)點(diǎn)力只為使模型處于某一力學(xué)平衡狀態(tài)，而程序并不判斷插入的節(jié)點(diǎn)力是否代表真實(shí)狀態(tài)［20］，即不檢驗(yàn)其是否對(duì)應(yīng)于變形體真實(shí)位移，故體力荷載作用下直接用無(wú)限元進(jìn)行初始應(yīng)力分析有可能獲得不正確的計(jì)算結(jié)果。

ABAQUS動(dòng)力無(wú)限元通過(guò)在有限元-無(wú)限元邊界節(jié)點(diǎn)設(shè)置一系列法向及切向的阻尼器，模擬半無(wú)限空間介質(zhì)的輻射阻尼［20］，其參考了Lysmer等［1］的黏性邊界理論，依據(jù)彈性介質(zhì)中S波、P波的波動(dòng)方程，推導(dǎo)出反射波為零條件下的邊界阻尼系數(shù)。由此可見，動(dòng)力計(jì)算中無(wú)限元人工邊界是一種無(wú)靜剛度的黏性吸收邊界，不能模擬地基的彈性恢復(fù)力以及內(nèi)行波的能量輸入。

針對(duì)ABAQUS靜、動(dòng)力無(wú)限元邊界所存在的問(wèn)題，論文開展了靜、動(dòng)力統(tǒng)一人工邊界改進(jìn)方法研究。靜力問(wèn)題從二維水平地基內(nèi)初始地應(yīng)力計(jì)算精度出發(fā)，探討無(wú)限元對(duì)體力荷載作用問(wèn)題的適用性及改進(jìn)方法，進(jìn)一步在此基礎(chǔ)上推導(dǎo)內(nèi)、外源振動(dòng)問(wèn)題相應(yīng)的邊界等效節(jié)點(diǎn)力公式，對(duì)應(yīng)利用半無(wú)限地基表面施加動(dòng)力荷載和底部輸入地震動(dòng)的兩種典型情況進(jìn)行驗(yàn)證。

2 無(wú)限元邊界初始地應(yīng)力生成精度探討

2.1 分析模型

由上節(jié)分析可知，將ABAQUS無(wú)限元直接應(yīng)用于體力荷載作用下的初始應(yīng)力狀態(tài)生成，可能會(huì)帶來(lái)不正確的結(jié)果。為此，建立二維彈性水平地基模型（圖1）進(jìn)行算例驗(yàn)證，模型側(cè)面及底面設(shè)置厚度L2=50 m的無(wú)限元人工邊界，有限元主模型區(qū)域50 m×50 m，網(wǎng)格尺寸0.5 m，有限元、無(wú)限元區(qū)域分別采用CPE4、CINPE4單元，地基土密度ρ=2 000 kg/m3、彈性模量E=200 MPa、泊松比v=0.25［10］。

圖1 計(jì)算模型示意圖Fig.1 Calculation model（unit：m）

對(duì)模型施加重力荷載，利用ABAQUS的Geostatic分析步進(jìn)行地應(yīng)力平衡，計(jì)算模型尺寸與計(jì)算精度呈正相關(guān)，采用模型尺寸調(diào)整系數(shù)ξ探討邊界位置對(duì)計(jì)算結(jié)果的影響，ξ定義為主模型調(diào)整后尺寸（L+2L1）與原尺寸（L0）之比即ξ=（L+2L1）/L0，分別按照ξ=1，2，3，4，5，10，15確定有限元計(jì)算區(qū)域，按0.5 m網(wǎng)格尺寸對(duì)模型進(jìn)行剖分。考察主模型中心沿路徑BE的豎向應(yīng)力，并與理論解答σyy=γ（50-y）作對(duì)比（γ為地基土重度，y為計(jì)算點(diǎn)縱坐標(biāo)），記錄計(jì)算耗時(shí)以研究模型尺寸對(duì)計(jì)算結(jié)果及效率的影響，具體方案如表1所示，其中K0為ξ=1并施加靜力邊界條件（兩側(cè)和底部分別采用滾軸、固定邊界）的計(jì)算結(jié)果。

表1 無(wú)限元邊界條件下模型尺寸對(duì)計(jì)算結(jié)果影響工況表Table 1 Effect size model calculation condition table membered infinite boundary conditions m

2.2 計(jì)算結(jié)果

經(jīng)地應(yīng)力平衡后，K0～K7模型位移量值均小于10-7m，僅從位移大小角度來(lái)看，地應(yīng)力平衡效果均較好。然而考察主模型區(qū)域豎向應(yīng)力分布（圖2，圖3）、計(jì)算精度及耗時(shí)（圖4）隨ξ的變化可知：

圖2 不同模型尺寸下主模型區(qū)域豎向應(yīng)力云圖Fig.2 Vertical stress contour of the main model with different model sizes

（1）施加傳統(tǒng)靜力邊界條件的K0模型計(jì)算結(jié)果與理論解完全一致，無(wú)限元邊界條件下計(jì)算精度受模型邊界遠(yuǎn)近影響，若地應(yīng)力平衡結(jié)果要滿足工程要求，邊界需要達(dá)到主模型區(qū)域15～20倍以上的遠(yuǎn)置距離。當(dāng)ξ=1時(shí)，豎向應(yīng)力明顯偏離自重應(yīng)力理論分布，呈現(xiàn)出沿中軸線BE向模型兩側(cè)邊界逐漸減小的形式，BE路徑豎向應(yīng)力與理論值的最大相對(duì)誤差達(dá)-33.6%；隨著ξ增大，邊界逐漸遠(yuǎn)置，豎向應(yīng)力分布趨于橫向均勻，逐漸接近理論（靜力邊界條件下）結(jié)果；當(dāng)ξ=15時(shí)，模型尺寸已達(dá)750 m×750 m，豎向應(yīng)力與理論值最大相對(duì)誤差仍有-6.1%，僅近似滿足工程精度要求。

（2）豎向應(yīng)力計(jì)算誤差與深度基本呈正相關(guān)，即ABAQUS靜力無(wú)限元自動(dòng)插入節(jié)點(diǎn)力的等效剛度與靜力邊界剛度的差異隨深度增加而變大。若需將誤差控制在工程允許誤差（5%以內(nèi)），ξ=1時(shí)，僅2.5 m深度范圍內(nèi)的計(jì)算結(jié)果可信（圖3）；而按傳統(tǒng)經(jīng)驗(yàn)取3～5倍主模型區(qū)域尺寸設(shè)置邊界進(jìn)行建模（圖4中粉紅色區(qū)域ξ=3～5），計(jì)算結(jié)果可信深度范圍為8.5～14 m，仍難以滿足主模型全區(qū)域的精度需要。

圖3 各模型沿路徑BE豎向應(yīng)力分布情況Fig.3 Vertical stress distribution along path BE

用圖4中所示的三次多項(xiàng)式對(duì)相對(duì)誤差曲線進(jìn)行擬合，得到尺寸調(diào)整系數(shù)ξ與誤差百分比ω的擬合函數(shù)（決定系數(shù)0.982）。根據(jù)擬合函數(shù)進(jìn)行外延，若誤差控制在5%，ξ應(yīng)不小于20，計(jì)算量急劇增加，因此有必要提出一種適用于主模型區(qū)域高精度、高效率分析的靜力無(wú)限元邊界改進(jìn)方法并使之能應(yīng)用于動(dòng)力計(jì)算的初始狀態(tài)生成。

圖4 豎向應(yīng)力誤差和計(jì)算耗時(shí)與尺寸調(diào)整系數(shù)關(guān)系Fig.4 The relationship between vertical stress error，calculation time and size adjustment coefficient

3 基于無(wú)限元的靜-動(dòng)力統(tǒng)一人工邊界改進(jìn)方法

3.1 靜力無(wú)限元邊界改進(jìn)

3.1.1 基于靜力邊界等效節(jié)點(diǎn)力的無(wú)限元?jiǎng)偠刃拚?/p>

體力荷載作用下，ABAQUS無(wú)限元靜力計(jì)算產(chǎn)生誤差的根本原因是程序引入的節(jié)點(diǎn)力等效剛度與靜力狀態(tài)的實(shí)際邊界剛度不匹配，因此，若考慮在施加無(wú)限元人工邊界的同時(shí)，對(duì)有限元-無(wú)限元邊界節(jié)點(diǎn)施加傳統(tǒng)靜力邊界（滾軸、固定）計(jì)算得到的邊界支反力，由于邊界支承剛度相等，主模型區(qū)域受力狀態(tài)應(yīng)與靜力邊界條件下完全一致。下面對(duì)其進(jìn)行理論驗(yàn)證：

靜力計(jì)算中，總體平衡方程為

式中：K為剛度矩陣；U為位移向量；F為節(jié)點(diǎn)力向量。

將式（1）寫成如下形式：

式中：KI、KIB、KBI、KBI為剛度矩陣的分塊子矩陣；UI、UB分別表示模型內(nèi)部與邊界節(jié)點(diǎn)的位移；FI、FB分別表示模型內(nèi)部與邊界節(jié)點(diǎn)的節(jié)點(diǎn)力。

由于邊界及內(nèi)部節(jié)點(diǎn)力一部分由外荷載引起，一部分由約束反力引起，式（2）可以進(jìn)一步表示為

式中：上標(biāo)f表示外荷載引起的節(jié)點(diǎn)力分量；而上標(biāo)r表示約束引起的節(jié)點(diǎn)力分量，在施加靜力邊界時(shí)即為靜力邊界約束引起的節(jié)點(diǎn)力，而施加無(wú)限元邊界時(shí)則為無(wú)限元邊界引起的節(jié)點(diǎn)力。

施加靜力邊界UB=0時(shí)，由于邊界節(jié)點(diǎn)在約束或指定位移方向不會(huì)對(duì)邊界以內(nèi)節(jié)點(diǎn)產(chǎn)生節(jié)點(diǎn)力［13］，故F rI=0，于是

可解得

將F rB代回式（3），仍可解得UB=0，且原UI的解也不變，上述推導(dǎo)從物理角度可解釋為對(duì)模型施加約束條件和對(duì)模型施加與約束條件UB=0相對(duì)應(yīng)的等效節(jié)點(diǎn)力F rB是等效的，因此可在需要去除靜力約束的情況下，施加相應(yīng)的約束反力以達(dá)到不改變?cè)o力平衡狀態(tài)的目的，下文將F rB稱為靜力等效節(jié)點(diǎn)力。

而當(dāng)施加無(wú)限元人工邊界時(shí)，程序自動(dòng)插入的節(jié)點(diǎn)力FINFB并不一定滿足FINFB=F rB，故將FINFB代回式（4）不一定得到真實(shí)位移邊界條件UB=0，計(jì)算所得的應(yīng)力狀態(tài)亦不一定與靜力邊界一致。若在施加無(wú)限元人工邊界的同時(shí)，對(duì)主模型區(qū)域與無(wú)限元交界面處一并施加靜力等效節(jié)點(diǎn)力F rB，則主模型區(qū)域的應(yīng)力場(chǎng)將保持與靜力邊界條件下一致。

3.1.2 模型驗(yàn)證

由上節(jié)推導(dǎo)可知，可引入靜力等效節(jié)點(diǎn)力以提高無(wú)限元的初始地應(yīng)力場(chǎng)計(jì)算精度，具體實(shí)現(xiàn)步驟為：（1）建立有限元模型，對(duì)模型施加靜力邊界條件，計(jì)算完成后，提取模型左、右及底部支座反力即靜力等效節(jié)點(diǎn)力。（2）去除靜力邊界，施加無(wú)限元人工邊界、導(dǎo)入步驟（1）中計(jì)算所得的應(yīng)力場(chǎng)并施加有限元-無(wú)限元靜力等效節(jié)點(diǎn)力即可完成初始地應(yīng)力生成，靜力等效節(jié)點(diǎn)力在后續(xù)分析過(guò)程中保持不變。

為方便對(duì)比，取第2.1節(jié)中尺寸調(diào)整系數(shù)ξ=1的模型進(jìn)行驗(yàn)證。將施加靜力等效節(jié)點(diǎn)力前后的計(jì)算結(jié)果、誤差連同理論解繪入圖5，由圖可知，僅設(shè)置無(wú)限元邊界的模型豎向應(yīng)力誤差可達(dá)-33.6%，而采用文中引入靜力等效節(jié)點(diǎn)力方法與理論解完全一致，說(shuō)明引入靜力等效節(jié)點(diǎn)力的方法是合理的。

圖5 引入靜力等效節(jié)點(diǎn)力方法與僅設(shè)置無(wú)限元邊界計(jì)算所得豎向應(yīng)力Fig.5 Vertical stress difference between the static equivalent nodal force and the infinite boundary

3.2 靜-動(dòng)力統(tǒng)一人工邊界等效節(jié)點(diǎn)力分析

3.2.1 內(nèi)源振動(dòng)問(wèn)題

對(duì)于需要考慮靜、動(dòng)力綜合效應(yīng)（即靜動(dòng)力合算）的內(nèi)源振動(dòng)問(wèn)題，動(dòng)力計(jì)算只需考慮外行波吸收，無(wú)需考慮內(nèi)行波輸入，而無(wú)限元本身可較好地模擬半無(wú)限空間的輻射阻尼以吸收外行波，因此僅需保證靜力計(jì)算精度，故內(nèi)源振動(dòng)問(wèn)題的等效節(jié)點(diǎn)力等于靜力等效節(jié)點(diǎn)力：式中：下標(biāo)α為主模型區(qū)域邊界位置；β為等效節(jié)點(diǎn)力方向。

3.2.2 外源振動(dòng)問(wèn)題

對(duì)于靜動(dòng)力合算的外源振動(dòng)問(wèn)題，則需滿足靜力計(jì)算精度并同時(shí)考慮外行波吸收和內(nèi)行波輸入。其中，外行波吸收可由無(wú)限元本身解決，對(duì)于內(nèi)行波輸入，需克服人工邊界的能量吸收才能由域外透射至域內(nèi)，文獻(xiàn)［11］考慮無(wú)限域與有限域之間的相互作用力和克服人工邊界的阻尼力，導(dǎo)出了外源振動(dòng)條件下用于地震動(dòng)輸入的等效節(jié)點(diǎn)力公式（下文稱之為動(dòng)力等效節(jié)點(diǎn)力）；靜力計(jì)算部分可采用上文第3.2.1節(jié)靜力等效節(jié)點(diǎn)力方法以保證精度，故外源振動(dòng)條件下，靜-動(dòng)力統(tǒng)一人工邊界等效節(jié)點(diǎn)力應(yīng)為模擬內(nèi)行波輸入的動(dòng)力等效節(jié)點(diǎn)力和保證靜力計(jì)算精度的靜力等效節(jié)點(diǎn)力共同組成。

將動(dòng)力等效節(jié)點(diǎn)力［11］與靜力等效節(jié)點(diǎn)力（式（6））疊加即得到外源振動(dòng)條件下靜-動(dòng)力統(tǒng)一人工邊界的等效節(jié)點(diǎn)力，下面以（以圖1為例）給出平面剪切S波和壓縮P波沿y方向從模型底部垂直向上傳播時(shí)，有限元-無(wú)限元邊界等效節(jié)點(diǎn)力計(jì)算公式。

平面P波垂直向上傳播時(shí)，考慮底部入射波f1和地表反射波f2，則無(wú)限元靜-動(dòng)力統(tǒng)一人工邊界等效節(jié)點(diǎn)力，在主模型左側(cè)AD邊界為：

主模型右側(cè)C F邊界：

主模型底部DF邊界：

同理可得剪切S波垂直向上傳播時(shí)等效節(jié)點(diǎn)力為：

主模型左側(cè)A D邊界：

主模型右側(cè)CF邊界：

主模型底部BE邊界：

式中：下標(biāo)l、r、b分別表示模型左側(cè)、右側(cè)和底面；x、y分別表示x，y方向；A i為對(duì)應(yīng)節(jié)點(diǎn)的控制面積；λ為第一Lamé常數(shù)。

3.3 實(shí)現(xiàn)流程及適用性驗(yàn)證

3.3.1 實(shí)現(xiàn)流程

靜、動(dòng)力等效節(jié)點(diǎn)力可分別基于式（6）和式（11）～式（22）進(jìn)行編程計(jì)算。對(duì)于內(nèi)源振動(dòng)問(wèn)題，等效節(jié)點(diǎn)力即為靜力等效節(jié)點(diǎn)力，對(duì)于外源振動(dòng)問(wèn)題，等效節(jié)點(diǎn)力則由靜力等效節(jié)點(diǎn)力和動(dòng)力等效節(jié)點(diǎn)力共同組成。完整的靜、動(dòng)力統(tǒng)一人工邊界實(shí)現(xiàn)流程如圖6所示，文中研究以Python語(yǔ)言編制等效節(jié)點(diǎn)力輸入程序。

圖6 靜-動(dòng)力統(tǒng)一人工邊界計(jì)算實(shí)現(xiàn)流程Fig.6 The implementation method of static-dynamic unified artificial boundary

3.3.2 內(nèi)源振動(dòng)問(wèn)題驗(yàn)證

仍考慮第2.1中的二維水平地基模型，靜力荷載僅為重力，在模型地表中心B點(diǎn)施加如圖7所示的動(dòng)荷載，峰值200 kN，作用時(shí)間0.5 s，計(jì)算時(shí)長(zhǎng)4 s，對(duì)模型施加內(nèi)源振動(dòng)問(wèn)題對(duì)應(yīng)的等效節(jié)點(diǎn)力（式6），求解時(shí)間步長(zhǎng)為0.001 s。為驗(yàn)證計(jì)算結(jié)果準(zhǔn)確性，另建立一遠(yuǎn)置邊界模型作為近似理論解答，根據(jù)文獻(xiàn)［19］的建議，監(jiān)測(cè)點(diǎn)到人工邊界的最近距離應(yīng)大于求解時(shí)間與波速乘積的一半，設(shè)置模型尺寸為1 280 m×640 m，采用傳統(tǒng)靜力邊界條件，材料參數(shù)、單元尺寸及類型均與主模型區(qū)域一致。首先進(jìn)行地應(yīng)力平衡，隨后進(jìn)行動(dòng)力計(jì)算。

圖7 輸入動(dòng)荷載時(shí)程Fig.7 Dynamic load time history

考慮到地應(yīng)力平衡后模型位移近似為零，位移響應(yīng)主要取決于動(dòng)力荷載，應(yīng)同時(shí)考察模型應(yīng)力與位移響應(yīng)，監(jiān)測(cè)模型B點(diǎn)的豎向位移及中心G點(diǎn)（深度24.75 m）的豎向應(yīng)力時(shí)程，分別如圖8（a）和（b）所示。由圖可知，采用靜-動(dòng)力統(tǒng)一人工邊界方法計(jì)算結(jié)果與遠(yuǎn)置邊界基本一致，精度良好；除此之外，采用文中方法耗時(shí)305 s，僅為遠(yuǎn)置邊界模型耗時(shí)（86 631 s）的3.52%，因此該方法應(yīng)用于靜、動(dòng)力合算問(wèn)題可極大提高計(jì)算效率。

圖8 沖擊荷載作用下B點(diǎn)位移及G點(diǎn)豎向應(yīng)力響應(yīng)時(shí)程Fig.8 Displacement response of point B and the vertical stress response of point G

3.3.3 外源振動(dòng)問(wèn)題驗(yàn)證

（1）計(jì)算模型及工況

基于第2.1中的二維水平地基模型，靜力荷載僅為重力，底部x方向根據(jù)地震位移時(shí)程曲線（圖9，最大幅值0.01 m），利用式（17）確定的等效節(jié)點(diǎn)力輸入地震波，計(jì)算時(shí)長(zhǎng)4 s，計(jì)算對(duì)比3種邊界形式下模型動(dòng)力響應(yīng)差異，其中，工況Ⅰ為傳統(tǒng)無(wú)限元人工邊界；工況Ⅱ?yàn)橹髂Ｐ蛡?cè)面施加文獻(xiàn)［11］所推得的動(dòng)力等效節(jié)點(diǎn)力，但不考慮靜力等效節(jié)點(diǎn)力；工況Ⅲ為主模型側(cè)面施加本文所提出的靜-動(dòng)力統(tǒng)一人工邊界。

圖9 輸入波位移時(shí)程曲線Fig.9 The displacement time history of seismic wave

首先進(jìn)行地應(yīng)力平衡，隨后進(jìn)行動(dòng)力計(jì)算。監(jiān)測(cè)模型頂部A、B、C點(diǎn)和底部E、F、G點(diǎn)位移時(shí)程以及內(nèi)部G點(diǎn)的豎向應(yīng)力時(shí)程。

（2）計(jì)算結(jié)果分析

由彈性波動(dòng)理論可知，當(dāng)?shù)卣鸩ù怪毕蛏蟼鞑ブ恋乇頃r(shí)，反射波與入射波幅值大小及方向均應(yīng)一致，地表位移幅值為輸入波的2倍；剪切波在傳播過(guò)程中，不會(huì)引起單元正應(yīng)力而僅引起切應(yīng)力，因此單元正應(yīng)力應(yīng)基本保持不變，即等于前期的靜力計(jì)算結(jié)果?？捎缮鲜?項(xiàng)規(guī)律，判斷各邊界條件下位移及應(yīng)力計(jì)算結(jié)果的合理性。

頂部A、B、C點(diǎn)及底部D、E、F點(diǎn)的位移響應(yīng)情況分別如圖10（a）～（f）所示，由圖可知：

（1）工況Ⅰ由于僅設(shè)置了無(wú)限元人工邊界而未考慮靜、動(dòng)力等效節(jié)點(diǎn)力，地表位移（圖10（a））響應(yīng)顯著偏小，同時(shí)模型底部（圖10（d））未出現(xiàn)地表反射波，均與理論響應(yīng)曲線差別較大，地表A、B、C三點(diǎn)間同步性較差，模型底部D、E、F點(diǎn)間也出現(xiàn)類似情況，說(shuō)明該邊界條件會(huì)造成波動(dòng)場(chǎng)的不均勻；

圖10 模型頂部及底部測(cè)點(diǎn)位移響應(yīng)時(shí)程Fig.10 Displacement response of the top and bottom points

（2）分別考慮了動(dòng)力等效節(jié)點(diǎn)力和靜、動(dòng)力等效節(jié)點(diǎn)力的工況Ⅱ、工況Ⅲ，頂部及底部位移響應(yīng)均與理論值基本相同，且地表A、B、C點(diǎn)及底部D、E、F點(diǎn)響應(yīng)時(shí)程曲線同步性較好，說(shuō)明波的傳播是均勻的，僅從位移角度來(lái)看，2種邊界條件下的計(jì)算精度均較好。

3種工況下，G點(diǎn)豎向正應(yīng)力響應(yīng)分別如圖11（a）～（c）所示，由圖可知：

圖11 模型中部G點(diǎn)單元豎向應(yīng)力響應(yīng)時(shí)程Fig.11 Vertical normal stress response of element in position G

（1）工況Ⅰ和工況Ⅱ單元豎向正應(yīng)力響應(yīng)量值（-326.2 kPa）均與理論值（-495.0 kPa）差異較大，其應(yīng)力波動(dòng)基準(zhǔn)值與理論值相比均較小，精度較低；

（2）工況Ⅲ單元豎向正應(yīng)力響應(yīng)量值均穩(wěn)定在理論值附近，波動(dòng)最大幅值僅為理論值的0.011%，可見，采用本文靜-動(dòng)力統(tǒng)一人工邊界獲得的應(yīng)力響應(yīng)是可靠的。

由以上分析可知，雖然工況Ⅱ邊界條件下的位移響應(yīng)與理論值基本一致，但應(yīng)力響應(yīng)出現(xiàn)較大失真。巖土動(dòng)力計(jì)算中，第一步往往是地應(yīng)力平衡，動(dòng)力加載前模型位移很小（＜10-5m），動(dòng)位移響應(yīng)主要取決于動(dòng)力荷載，然而由于初始地應(yīng)力的誤差影響后續(xù)動(dòng)應(yīng)力響應(yīng)結(jié)果精度。

3.3.4 邊界條件適用性探討

從監(jiān)測(cè)點(diǎn)位移及應(yīng)力響應(yīng)與理論解對(duì)比可知，工況Ⅰ位移、應(yīng)力響應(yīng)量值及規(guī)律均明顯偏離理論解，因此傳統(tǒng)直接設(shè)置無(wú)限元的邊界條件不適用于外源振動(dòng)問(wèn)題；工況Ⅱ位移響應(yīng)與理論解基本一致但初始地應(yīng)力出現(xiàn)較大系統(tǒng)誤差，因此得到的動(dòng)應(yīng)力響應(yīng)并不可靠，工況Ⅱ邊界條件僅適用于靜、動(dòng)力分算的外源振動(dòng)問(wèn)題；采用文中的無(wú)限元靜-動(dòng)力統(tǒng)一人工邊界（工況Ⅲ）的位移及應(yīng)力響應(yīng)均與理論解相符，且可根據(jù)求解問(wèn)題的類型考慮按需引入靜力、動(dòng)力等效節(jié)點(diǎn)力，使得其對(duì)靜力問(wèn)題和靜、動(dòng)力分算（合算）的內(nèi)、外源振動(dòng)問(wèn)題均具有較強(qiáng)適用性。綜上，3種工況對(duì)應(yīng)的邊界條件對(duì)應(yīng)的適用范圍如表2所示。

表2 常用人工邊界設(shè)置方法及其對(duì)巖土靜、動(dòng)力問(wèn)題適用性Table 2 The common artificial boundary′s applicability to rock and soil static and dynamic problems

4 結(jié)論

通過(guò)分析ABAQUS靜、動(dòng)力無(wú)限元理論，利用二維水平地基模型，分析討論了無(wú)限元直接應(yīng)用于體力荷載初始應(yīng)力計(jì)算時(shí)的精度較差且對(duì)模型尺寸要求大的局限性，針對(duì)性的提出了引入靜力等效節(jié)點(diǎn)力的改進(jìn)方法，建立了無(wú)限元靜-動(dòng)力統(tǒng)一人工邊界并進(jìn)行了算例驗(yàn)證，得出以下結(jié)論：

（1）因ABAQUS靜力無(wú)限元引入節(jié)點(diǎn)力等效剛度與靜力邊界的差異，無(wú)限元邊界直接應(yīng)用于體力荷載下初始地應(yīng)力的計(jì)算結(jié)果誤差較大。算例表明，50 m×50 m水平地基的豎向自重應(yīng)力計(jì)算結(jié)果嚴(yán)重偏小，最大相對(duì)誤差可達(dá)-33.6%；邊界遠(yuǎn)置方法雖能一定程度改善精度，但所需的模型尺寸過(guò)大、計(jì)算效率極低。故有必要研究適用于無(wú)限元人工邊界的高精度、高效率的改進(jìn)方法。

（2）以無(wú)限元應(yīng)用于靜力工況及內(nèi)、外源振動(dòng)的普遍適用性為目標(biāo)，基于模型約束條件和對(duì)應(yīng)節(jié)點(diǎn)力的等效性，提出引入靜力等效節(jié)點(diǎn)力提高初始應(yīng)力計(jì)算精度的方法，并通過(guò)模型計(jì)算驗(yàn)證了其具有與傳統(tǒng)靜力邊界一致的效果；進(jìn)一步基于彈性波動(dòng)理論考慮有限域-無(wú)限域相互作用和邊界阻尼效應(yīng)，推導(dǎo)了內(nèi)、外源振動(dòng)條件下的等效節(jié)點(diǎn)力，提出了無(wú)限元靜-動(dòng)力統(tǒng)一人工邊界并編制了相應(yīng)的等效節(jié)點(diǎn)力輸入程序。

（3）針對(duì)半無(wú)限地基表面的動(dòng)力加載和地基深部振動(dòng)向地表傳播的兩種典型內(nèi)、外源振動(dòng)情況，以模型位移與應(yīng)力響應(yīng)為監(jiān)測(cè)量驗(yàn)證了無(wú)限元靜-動(dòng)力統(tǒng)一人工邊界的有效性，分析了傳統(tǒng)無(wú)限元與僅考慮動(dòng)力等效節(jié)點(diǎn)力的邊界條件對(duì)靜、動(dòng)力合算問(wèn)題的局限性，論文提出的無(wú)限元靜-動(dòng)力統(tǒng)一人工邊界對(duì)內(nèi)、外源振動(dòng)的靜、動(dòng)力分（合）算問(wèn)題均可具有較好的適用性。