付 浩，仝 睿，宋二祥

(清華大學(xué)土木工程系土木工程安全與耐久教育部重點(diǎn)實(shí)驗(yàn)室，北京 100084)

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

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

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

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

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

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

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

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

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

圖1 軸對(duì)稱無(wú)限平面圓盤Fig.1 Axisymmetric infinite plane disk

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

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

1.1 粘彈性傳輸邊界一

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

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

由此可得到SH 波的波動(dòng)方程：

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

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

1.2 粘彈性傳輸邊界二

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

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

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

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

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

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

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

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

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

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

2 數(shù)值計(jì)算驗(yàn)證

2.1 模型及參數(shù)

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

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

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

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

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

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

2.2 突加荷載下的計(jì)算結(jié)果分析

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

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

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

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

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

2.3 脈沖荷載下的結(jié)果分析

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

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

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

3 結(jié)論

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