付 浩，仝 睿，宋二祥

(清華大學土木工程系土木工程安全與耐久教育部重點實驗室，北京 100084)

有限元數(shù)值模擬方法是解決土動力問題的重要手段。采用有限元方法進行動力問題計算，一般只能選取有限范圍的計算網(wǎng)格，同時在網(wǎng)格的人為截斷邊界處設定特殊的傳輸邊界，以實現(xiàn)對真實無限域上波動的模擬。近年來，國內(nèi)外對人工傳輸邊界問題已有了廣泛的研究，文獻[1?2]分別對各類傳輸邊界的成果與進展進行了綜述。

Lysmer 和Kuhlemeyer[3]最先提出的粘性邊界通過對一維平面波的分析導出，其做法是在人為截斷邊界上布設阻尼器。因其概念清晰，應用方便，而被廣泛采用。但粘性邊界存在一定缺陷：一方面，粘性邊界的推導中未考慮波的振動幅值隨波的擴散發(fā)生幾何衰減；另一方面，粘性邊界處不能考慮外部介質(zhì)的靜剛度，這就使其存在低頻穩(wěn)定性問題。

Deeks 和Randolph[4]通過分析水平傳播的柱面波給出徑向及豎向波動分析的粘彈性傳輸邊界，相應邊界條件是在邊界上布設阻尼器的同時，也布設彈簧。相比于粘性邊界，粘彈性邊界能模擬人工邊界外介質(zhì)的彈性恢復性能，具有良好的高頻和低頻穩(wěn)定性。劉晶波等[5]對Deeks 和Randolph提出的人工邊界進一步發(fā)展完善，并通過對球面波的分析提出三維粘彈性邊界。

但Deeks 和Randolph 沒有給出針對質(zhì)點在柱面環(huán)向振動波，亦即SH 波的傳輸邊界條件。而這種邊界條件對動力機器基礎扭轉振動以及地震作用下不規(guī)則建筑發(fā)生的扭轉振動進行分析時是需要的。實際上，此類扭轉振動在分層地基中引起的Love 波正是以豎直圓柱面為波陣面沿水平方向向外傳播的。

劉光磊、宋二祥[6]將Deeks 和Randolph 針對單相介質(zhì)的人工邊界拓展到飽和多孔兩相介質(zhì)，建立了飽和兩相介質(zhì)時域動力固結分析的粘彈性人工邊界，計算表明其總體效果較好。但其中針對SH 波的人工邊界條件，針對微分方程近似解的假設偏于簡化，尚有進一步優(yōu)化的空間。

Du 和Zhao[7]基于無限域動力剛度有理近似的連分式展開技術， 建立了針對一般彈性介質(zhì)的高階時域粘彈性傳輸邊界。Zhao 等[8]建立了針對多層介質(zhì)的高階時域粘彈性傳輸邊界。Li 和Song[9]則針對飽和多孔兩相介質(zhì)建立了高階時域粘彈性傳輸邊界，并通過數(shù)值模擬驗證了此類邊界的計算精度。但對一般彈性介質(zhì)中SH 波分析的粘彈性邊界的深入討論還鮮見報道。

本文假設波動以柱面波的形式傳播，推導針對一般彈性介質(zhì)中SH 波進行數(shù)值分析的粘彈性傳輸邊界。由于采用近似解，依推導過程的不同，給出兩種不同的人工邊界，并對這兩種人工邊界的合理有效性進行理論分析和數(shù)值模擬驗證。由于在彈性假設下，剪切波并不會引起超靜孔壓，因此本文的邊界既適用于單相體，也可推廣到飽和土地基等兩相體中。將此人工邊界條件與已有針對P 波、SV 波的人工邊界條件結合可以實現(xiàn)三維波動問題的有限元分析。

1 柱面波平面內(nèi)剪切粘彈性傳輸邊界構造

假設波動以柱面波的形式傳播，采取如圖1所示無限大單位厚度圓盤進行分析。選用柱坐標系，坐標軸分別為r、θ、z，并記各方向的位移分別為ur、uθ、uz。震源位于圓盤中心。

圖1 軸對稱無限平面圓盤Fig.1 Axisymmetric infinite plane disk

柱面波可解耦為P 波、SH 波和SV 波三類波，分別對應上述圓盤平面內(nèi)徑向、平面內(nèi)剪切和平面外剪切振動的運動形式[10]。對于P 波和SV波，Deeks 和Randolph 已給出了在邊界處位移與應力的關系式，并構造了合理有效的粘彈性傳輸邊界，分別如式(1)、式(2)：

以下采用兩種略有不同的兩種推導方法，分別給出圓盤平面內(nèi)剪切波動分析的兩種粘彈性傳輸邊界。

1.1 粘彈性傳輸邊界一

按給定柱坐標系，設震源在 θ向振動，振動將以SH 波的形式向四周傳播。在柱坐標系中，列出振動的平衡方程、物理方程及幾何方程如式(3)。由于振動具有繞軸中心對稱的性質(zhì)，因此各物理量對 θ的導數(shù)為0。

式中，τrθ為圓盤平面內(nèi)切向的剪應力。

由此可得到SH 波的波動方程：

為求解以上方程，引入如式(5)定義的位移勢函數(shù) ?：

圖2 第一種粘彈性傳輸邊界模型Fig.2 The first viscoelastic transport boundary model

1.2 粘彈性傳輸邊界二

除了前文提到的推導方法外，使用另一種轉化方法，可以得到不同的粘彈性傳輸邊界。注意到：

注意到此時邊界處的剪切應力與切向位移的關系與Deeks 和Randolph 文章中平面內(nèi)法向的邊界條件表達式接近，故可采用如圖3 所示的元件模型以模擬該邊界條件，選取適當?shù)脑?shù)，與元件系統(tǒng)相連的邊界微元面的應力與位移可滿足式(17)的關系。對于距震源R處的平面內(nèi)切向邊界，其單位面積上模擬遠場介質(zhì)作用的等效物理元件的力學參數(shù)分別為k=2G/R，c=ρVS，m=2ρR。

圖3 第二種粘彈性傳輸邊界模型Fig.3 The second viscoelastic transport boundary model

1.3 上列兩人工邊界的合理性論證

以上對同一問題的分析得出兩種不同的人工邊界條件。那么哪一種更為合理呢？本節(jié)試圖回答此問題。這里提出的一個觀點是，當振動頻率很低時，人工邊界條件應退化為靜力邊界條件，也就是人工邊界的剛度應等于靜力剛度。符合此要求的人工邊界條件，至少其低頻精度較好。

由于一般彈性力學教科書中沒有給出此問題的靜剛度，這里給出相應的推導。選取有孔圓盤，內(nèi)徑為r0，圓盤在孔內(nèi)側受平面內(nèi)切應力q0。由微元體的平衡方程可得到：

式中，C1、C2為與邊界條件有關的待定系數(shù)。

考慮無限大有孔圓盤，其邊界條件為：

與前面給出的粘彈性邊界條件對比可見，第二種傳輸邊界其彈簧剛度與靜剛度相等，其低頻精度會較好，在理論上也應該更為精確。隨后的數(shù)值模擬計算也驗證了這一推斷。

事實上，第二種傳輸邊界同樣通過勢函數(shù)的近似解推出，但卻得到了與第一種傳輸邊界不同的結果。這是由于式(7)并非勢函數(shù)的精確解，隨后利用此近似解進行推導的不同代入過程即會得到不同的結果。從數(shù)學的角度來說，二者均為近似解，但第二種方法恰可以滿足靜剛度的要求，且具有更精簡的表達式與元件模型。

2 數(shù)值計算驗證

2.1 模型及參數(shù)

這里使用COMSOL 計算軟件對無限體中的柱面波進行模擬，以驗證上述兩種粘彈性傳輸邊界的有效性，并分析其精度。算例考慮二維軸對稱平面應變問題，使用圓環(huán)網(wǎng)格，內(nèi)徑為10 m，外徑為50 m。計算模型如圖4 所示。

材料參數(shù)參考中硬砂土層選取，G=79.6 MPa，材料密度ρ=2650 kg/m3，由此可知波速V=173 m/s。

對人為截斷邊界分別設為粘彈性傳輸邊界一、粘彈性傳輸邊界二和固定邊界進行對比分析。在內(nèi)邊界分別施加如圖5 所示兩種扭矩荷載。計算時長為1 s。同時還給出網(wǎng)格足夠大(外徑500 m)，計算時間內(nèi)反應不受邊界影響的結果作為真實解的參考值，這與文獻[12?13]中參考解的處理方法一致。

圖4 模型網(wǎng)格示意圖Fig.4 Model grid diagram

圖5 荷載示意圖Fig.5 Load diagram

選取模型內(nèi)徑處(r=10 m)、A點(r=20 m)和B點(r=40 m)處數(shù)據(jù)進行分析。

2.2 突加荷載下的計算結果分析

圖6 為模型在突加荷載作用下，模型內(nèi)徑處(r=10 m)的環(huán)向位移計算結果。將位移值除以F0R/G，進行無量綱化處理。其中，F(xiàn)0為荷載峰值，R為模型外邊徑，此處為50 m。比較了固定邊界與本文提到的兩種傳輸邊界同參考解之間的差別。結果表明，在突加荷載作用下，模型內(nèi)徑處變形量不斷發(fā)展。在反射波回傳之前，固定邊界模型的變形量相差不大。但是當反射波傳回之后，固定邊界模型與參考解之間可以看到明顯的差別。而兩類傳輸邊界模型都和參考解較為接近，且第二種傳輸邊界模型更為精確。

圖7 和圖8 顯示了在突加荷載作用下，三種模型中A點和B點的位移大小及其與參考解的比較。可以看出，距荷載點越遠，固定邊界與參考解之間的偏差會越大。同時固定邊界偏差開始明顯的時間點會受到反射波到達時間的影響。而兩種傳輸邊界條件依然可以較好地反映外部土作用，得到的結果和參考解較為接近。第二種傳輸邊界模型同樣更為精確。

圖6 突加荷載下內(nèi)徑處環(huán)向位移Fig.6 Circular displacement at inner diameter under sudden load

圖7 突加荷載下A 點環(huán)向位移Fig.7 Circular displacement of point A under sudden load

圖8 突加荷載下B 點環(huán)向位移Fig.8 Circular displacement of point B under sudden load

2.3 脈沖荷載下的結果分析

圖9 和圖10 給出了沖擊荷載作用下A點和B點的環(huán)向位移計算結果。在沖擊荷載作用下，位移量更小，固定邊界模型會出現(xiàn)明顯的震蕩情況，而兩種傳輸邊界則可以較好地模擬實際情況。同樣，第二種傳輸邊界得到的結果會更為精確。

圖9 沖擊荷載下A 點環(huán)向位移Fig.9 Circular displacement of point A under impact load

圖10 沖擊荷載下B 點環(huán)向位移Fig.10 Circular displacement of point B under impact load

3 結論

本文通過分析水平傳播的柱面波，推導給出了針對一般彈性無限介質(zhì)內(nèi)SH 波傳播問題數(shù)值分析的粘彈性邊界，對所給人工邊界的合理有效性進行了理論分析和數(shù)值模擬檢驗。由于在推導過程采用相應微分方程的近似解，不同的推導過程可給出不同的人工邊界條件，兩種人工邊界的計算結果均可達到滿意的精度，但滿足靜力剛度條件的第二種傳輸邊界精度更高。盡管本文只分析了單相體的情況，但由于在彈性假設下，剪切波并不會引起超靜孔壓，因此本文的邊界既適用于單相體，也可推廣到飽和土地基等兩相體中。