基于無(wú)厚度Goodman單元的粘彈性人工邊界模擬

羅 佩， 劉國(guó)明， 王文君, 2, 徐辰奎

(1. 福州大學(xué)土木工程學(xué)院, 福建 福州 350116； 2. 漳州職業(yè)技術(shù)學(xué)院建筑工程系, 福建 漳州 363000)

?

基于無(wú)厚度Goodman單元的粘彈性人工邊界模擬

羅 佩1， 劉國(guó)明1， 王文君1, 2, 徐辰奎1

(1. 福州大學(xué)土木工程學(xué)院, 福建 福州 350116； 2. 漳州職業(yè)技術(shù)學(xué)院建筑工程系, 福建 漳州 363000)

基于粘彈性人工邊界, 推導(dǎo)了三維一致粘彈性人工邊界單元的剛度矩陣和阻尼矩陣. 利用邊界單元厚度對(duì)計(jì)算結(jié)果影響不大， 引入了一種新的一致粘彈性單元的模擬方法， 即無(wú)厚度Goodman單元. 基于無(wú)厚度Goodman單元對(duì)一致粘彈性邊界單元進(jìn)行模擬， 并通過(guò)半無(wú)限空間和均勻半空間的算例進(jìn)行驗(yàn)證. 結(jié)果表明， 使用無(wú)厚度Goodman單元實(shí)現(xiàn)的一致粘彈性邊界單元簡(jiǎn)單、 實(shí)用性強(qiáng)， 并且可以得到足夠的精度. 相比于有厚度的邊界單元， 無(wú)厚度Goodman單元符合邊界單元無(wú)厚度的實(shí)際情況， 理論完善， 且對(duì)于可能存在的邊界不規(guī)則情況適用性強(qiáng).

粘彈性人工邊界； 無(wú)厚度Goodman單元； 人工邊界單元； 動(dòng)力計(jì)算

0　引言

在壩體-地基系統(tǒng)動(dòng)力響應(yīng)分析中， 地基的輻射阻尼的影響是一個(gè)備受關(guān)注的問(wèn)題. 用有限化方法對(duì)無(wú)限地基進(jìn)行處理的時(shí)候， 通過(guò)在劃定的有限范圍地基介質(zhì)中設(shè)置虛擬的邊界， 并確定這些虛擬邊界所需滿足的相應(yīng)條件， 該邊界稱之為人工邊界. 人工邊界條件(artificial boundary condition)的概念是由Alterman[1]最早在1968年提出的. 目前對(duì)人工邊界的研究主要分為全局人工邊界條件和局部人工邊界條件兩大方面. 局部人工邊界條件適應(yīng)性強(qiáng)， 是時(shí)域離散形式的邊界條件， 具備時(shí)空解耦的特性， 形式簡(jiǎn)單、 計(jì)算量小、 計(jì)算耗時(shí)少， 便于通過(guò)編程實(shí)現(xiàn)， 因此， 被廣泛應(yīng)用于有限元方法中.

粘彈性人工邊界是一種連續(xù)的局部人工邊界. 1994年， Deeks等[2]基于粘性邊界的基礎(chǔ)上提出了粘彈性人工邊界. 粘彈性人工邊界等效于將連續(xù)分布的彈簧元件和阻尼器元件設(shè)置于人工截?cái)噙吔缟希?與粘性邊界相比， 模擬出了散射波輻射和地基彈性恢復(fù)性能. 其概念清晰， 便于有限元實(shí)現(xiàn)， 克服了粘性邊界的低頻失穩(wěn)問(wèn)題， 能夠模擬地基的彈性恢復(fù)能力， 具有較高的計(jì)算效率和準(zhǔn)確度， 以及具備良好的魯棒性等優(yōu)點(diǎn)， 已廣泛應(yīng)用于土-結(jié)構(gòu)動(dòng)力相互作用問(wèn)題的分析中. 引入無(wú)厚度的Goodman單元模擬二維和三維一致粘彈性人工邊界單元， 取得了較好的效果. 無(wú)厚度Goodman單元與有厚度的邊界單元相比滿足了邊界單元無(wú)厚度的實(shí)際情況， 理論完善， 且對(duì)于可能存在的邊界不規(guī)則情況適用性強(qiáng)， 施加簡(jiǎn)單.

1　粘彈性人工邊界

根據(jù)球坐標(biāo)系內(nèi)的球面波波動(dòng)方程的推導(dǎo)， 人工邊界可以等效為連續(xù)分布的并聯(lián)彈簧-阻尼器系統(tǒng). 因此， 粘彈性邊界的首要問(wèn)題是選取適當(dāng)?shù)膹椈蓜偠群妥枘嵯禂?shù). 切向和法向的彈簧剛度和阻尼系數(shù)分別按照式(1)、 (2)來(lái)取值.

(1)

(2)

其中： KBT、 KBN分別是切向和法向彈簧剛度； CBT、 CBN分別是切向和法向的阻尼系數(shù)； R是人工邊界點(diǎn)至波源的距離； cs、 cP分別是S波和P波的波速； G是介質(zhì)的剪切模量； ρ是介質(zhì)的質(zhì)量密度； αT、 αN分別是粘彈性人工邊界切向和法向的修正系數(shù)， 推薦取0.67、 1.33[3-4].

2　一致粘彈性邊界單元

在進(jìn)行有限元離散時(shí)， 可用有限元形函數(shù)將連續(xù)分布的物理元件轉(zhuǎn)化為耦聯(lián)的邊界單元， 該方法稱為一致粘彈性人工邊界； 也可以通過(guò)對(duì)邊界單元進(jìn)行集中處理， 形成解耦的阻尼器和彈簧元件， 該方法稱為集中粘彈性人工邊界. 劉晶波等[5]將一致粘彈性人工邊界與集中粘彈性人工邊界進(jìn)行對(duì)比， 得出了采用一致粘彈性人工邊界時(shí)的計(jì)算精度較高的結(jié)論.

相對(duì)于采用彈簧-阻尼器元件， 等效實(shí)體邊界單元對(duì)粘彈性人工邊界的模擬更容易實(shí)現(xiàn)， 而且對(duì)人為劃分地基范圍時(shí)可能存在的邊界不規(guī)則情況， 這種邊界單元的適應(yīng)性更強(qiáng). 通過(guò)比較不同等效人工邊界單元厚度h對(duì)計(jì)算結(jié)果的影響， 發(fā)現(xiàn)邊界單元厚度對(duì)計(jì)算結(jié)果的影響不大[3, 5]， 厚度h可以在一個(gè)較大的范圍內(nèi)靈活取值. 當(dāng)然， 由于人工邊界是沒(méi)有厚度的， 采用有厚度的實(shí)體邊界單元模擬畢竟不夠完善， 因此在進(jìn)行有限元?jiǎng)恿τ?jì)算程序的編制時(shí)， 引入無(wú)厚度的Goodman單元模擬一致粘彈性人工邊界單元， 取得了較好的效果. 下面對(duì)無(wú)厚度Goodman單元對(duì)一致粘彈性邊界單元的模擬進(jìn)行介紹.

2.1三維一致粘彈性邊界模型單元?jiǎng)哦扰c阻尼矩陣

采用無(wú)厚度Goodman單元進(jìn)行模擬二維和三維一致粘彈性邊界單元. 三維模型單元采用兩種無(wú)厚度Goodman單元， 即無(wú)厚度六面體8結(jié)點(diǎn)單元和無(wú)厚度五面體6結(jié)點(diǎn)單元. 采用等參單元形式表示這兩種單元， 如圖1所示：

無(wú)厚度六面體8結(jié)點(diǎn)單元形函數(shù)為：

(3)

無(wú)厚度五面體6結(jié)點(diǎn)單元形函數(shù)為：

(4)

位移模式：

(5)

坐標(biāo)變換式：

(6)

(7)

其中：

上式， 8結(jié)點(diǎn)單元時(shí)，d=4； 6結(jié)點(diǎn)單元時(shí)，d=3.

(8)

其中：

三維的勁度矩陣和阻尼矩陣為二維積分， 采用高斯積分確定.

整體坐標(biāo)系下的單元?jiǎng)哦染仃噆和阻尼矩陣c為:

(9)

其中， 轉(zhuǎn)換矩陣：

結(jié)點(diǎn)轉(zhuǎn)換矩陣：

高斯點(diǎn)雅克比行列式為：

局部坐標(biāo)與整體坐標(biāo)的關(guān)系為：

(10)

(11)

對(duì)于二維的模型情況易于采用退化的三維模型得到， 不再進(jìn)行公式推導(dǎo).

3　數(shù)值算例

采用自主開發(fā)的適用于二、 三維混凝土壩的靜動(dòng)力分析的動(dòng)力線性有限元程序DYCDAM.f. 三維實(shí)體單元采用六面體8結(jié)點(diǎn)、 五面體6結(jié)點(diǎn)和四面體4結(jié)點(diǎn)等參單元， 采用無(wú)厚度的六面體8結(jié)點(diǎn)單元和五面體6結(jié)點(diǎn)單元模擬人工邊界. 靜荷載包括自重、 水荷載、 法向面荷載、 集中力的作用， 靜力計(jì)算可以考慮分期施工與分期蓄水過(guò)程. 采用Wilson-θ法進(jìn)行動(dòng)力分析， 動(dòng)荷載包括地震波自由場(chǎng)輸入、 給定結(jié)點(diǎn)的動(dòng)荷載時(shí)程以及人工邊界的等效地震動(dòng)輸入. 與動(dòng)水附加質(zhì)量矩陣結(jié)合， 可考慮壩-水動(dòng)力相互作用. 通過(guò)輸入2的維數(shù)， 自動(dòng)退化為二維計(jì)算程序. 以下通過(guò)兩個(gè)算例， 分別建立二維和三維模型進(jìn)行動(dòng)力計(jì)算分析， 并與解析解進(jìn)行對(duì)比， 以說(shuō)明粘彈性人工邊界對(duì)波動(dòng)問(wèn)題處理是有效的. 同時(shí), 通過(guò)算例驗(yàn)證本程序的可行性.

3.1二維半無(wú)限空間自由邊界受沖擊荷載作用問(wèn)題

建立二維半無(wú)限空間模型， 模型尺寸L=400 m，H=200 m. 根據(jù)對(duì)稱性， 取一半即L=200 m，H=200 m進(jìn)行計(jì)算， 計(jì)算區(qū)域如圖2(a)所示. 材料的彈性模量為E=0.1 GPa， 泊松比ν=0.22， 密度ρ=2 000 kg·m-3， 介質(zhì)的水平向波速取為143 m·s-1. 區(qū)域頂部為自由表面， 在區(qū)域的兩側(cè)邊和底邊處分別施加固定邊界條件和等效粘彈性邊界條件. 網(wǎng)格尺寸為10 m×10 m， 固定邊界條件模型剖分為441個(gè)結(jié)點(diǎn)， 400個(gè)平面四結(jié)點(diǎn)四邊形等參單元； 等效粘彈性邊界模型剖分為482個(gè)結(jié)點(diǎn)， 440個(gè)平面四結(jié)點(diǎn)四邊形等參單元. 為獲得擬精確解， 將截?cái)噙吔缭O(shè)置在離計(jì)算區(qū)域足夠遠(yuǎn)的地方， 即將模型擴(kuò)大為L(zhǎng)=1 600 m，H=800 m(計(jì)算時(shí)取一半大小即L=800 m，H=800 m)， 以保證在計(jì)算時(shí)間內(nèi)計(jì)算關(guān)心的區(qū)域不會(huì)有沖擊荷載傳遞到截?cái)噙吔绾蠓瓷浠貋?lái)的波動(dòng)影響. 網(wǎng)格尺寸為10 m×10 m， 共剖分6 561個(gè)結(jié)點(diǎn)， 6 400個(gè)平面四結(jié)點(diǎn)四邊形等參單元.

在A點(diǎn)處施加如圖2(b)所示的沖擊荷載， 計(jì)算分析時(shí)間步長(zhǎng)為0.01 s， 共計(jì)算600步， 歷時(shí)6 s.A點(diǎn)(0， 0)和B點(diǎn)(190， 0)處的位移時(shí)程曲線圖如圖3所示.

從圖3中可以看出， 將截?cái)噙吔缭O(shè)置在離計(jì)算區(qū)域足夠遠(yuǎn)時(shí)， 在計(jì)算時(shí)間內(nèi)計(jì)算關(guān)心的區(qū)域不會(huì)有沖擊荷載傳遞到截?cái)噙吔绾蠓瓷浠貋?lái)的波動(dòng)影響， 因此觀測(cè)點(diǎn)不會(huì)受到邊界反射波的影響. 采用固定邊界時(shí)， 由于計(jì)算區(qū)域的范圍較小， 在計(jì)算時(shí)間內(nèi)沖擊荷載傳遞到固定邊界時(shí)會(huì)產(chǎn)生反射， 進(jìn)而對(duì)各觀測(cè)點(diǎn)的位移時(shí)程曲線產(chǎn)生較大的影響， 其位移產(chǎn)生了震蕩， 與擬精確解相比誤差較大. 而在截?cái)噙吔缟喜捎玫刃д硰椥匀斯み吔纾?在計(jì)算區(qū)域與固定邊界模型相同的情況下， 各個(gè)觀測(cè)點(diǎn)的位移不會(huì)受到邊界反射波的影響. 通過(guò)對(duì)比計(jì)算得出的位移時(shí)程曲線與擬精確解， 可以看出二者非常接近.

算例證明了本程序的合理性， 同時(shí)也證明了采用粘彈性邊界在處理二維半無(wú)限空間的內(nèi)源波動(dòng)問(wèn)題時(shí)是有效的.

3.2經(jīng)典Lamb問(wèn)題

考慮經(jīng)典Lamb問(wèn)題， 建立三維半無(wú)限空間的計(jì)算模型， 計(jì)算時(shí)各單位為無(wú)量綱. 模型尺寸為1×1×0.5， 根據(jù)對(duì)稱性， 在水平方向上取其四分之一大小即0.5×0.5×0.5， 如圖4(a)所示， 共劃分為1 662個(gè)結(jié)點(diǎn)和1 300個(gè)單元. 在x、y、z方向上劃分為10×10×10的網(wǎng)格， 網(wǎng)格尺寸為0.05. 模型的剪切模量G=16， 泊松比ν=0.25， 介質(zhì)密度ρ=1， 時(shí)間步長(zhǎng)為0.01. 等效一致粘彈性邊界單元通過(guò)在模型的四周和底面設(shè)置厚度為0的三維Goodman單元， 并將人工邊界單元的最外層結(jié)點(diǎn)固定. 輸入的荷載為作用于半空間表面的集中荷載， 如圖4(b)所示.

計(jì)算時(shí)， 人工邊界參數(shù)αT=0.67，αN=1.33. S波和P波波速均取4 m·s-1. 計(jì)算結(jié)果如圖5所示， 取模型上的觀測(cè)點(diǎn)分別距加載中心O點(diǎn)距離為r=0.20和r=0.40， 分別計(jì)算出各個(gè)觀測(cè)點(diǎn)在解析解、 粘彈性人工邊界和固定邊界下的z向位移反應(yīng).

從圖5所示的位移時(shí)程曲線可以看出， 對(duì)應(yīng)不同觀測(cè)點(diǎn)， 當(dāng)采用固定邊界時(shí)， 位移產(chǎn)生了較大幅度的波動(dòng)， 相對(duì)于解析解存在著很大的誤差； 而使用一致粘彈性人工邊界時(shí)， 在接近加載點(diǎn)附近的結(jié)點(diǎn)和靠近邊界的結(jié)點(diǎn)上， 其計(jì)算結(jié)果與解析解非常接近， 計(jì)算的精度較高. 同時(shí)也可以看出， 在對(duì)粘彈性人工邊界的模擬上， 采用等效實(shí)體粘彈性邊界單元代替彈簧-阻尼器， 將使施加更為簡(jiǎn)單， 而且可以在邊界不規(guī)則的情況下使用.

4　結(jié)論

基于粘彈性人工邊界推導(dǎo)了三維一致粘彈性人工邊界單元的剛度矩陣和阻尼矩陣. 通過(guò)理論推導(dǎo)引入了一種新的一致粘彈性單元的模擬方法， 即采用無(wú)厚度Goodman單元實(shí)現(xiàn)粘彈性邊界單元的模擬. 數(shù)值算例表明， 使用無(wú)厚度Goodman單元實(shí)現(xiàn)的一致粘彈性邊界單元施加簡(jiǎn)單、 實(shí)用性強(qiáng)， 并且具有足夠的精度. 無(wú)厚度Goodman單元與有厚度的邊界單元相比滿足了邊界單元無(wú)厚度的實(shí)際情況， 理論完善， 且對(duì)于可能存在的邊界不規(guī)則情況適用性強(qiáng)， 施加簡(jiǎn)便， 體現(xiàn)出一致粘彈性人工邊界在進(jìn)行動(dòng)力計(jì)算時(shí)的優(yōu)越性.

[1] ALTERMAN Z S, KARAL F C. Propagation of elastic waves in layered media by finite-difference methods[J]. Bull Seism Soc Am, 1968, 58(1): 367-398.

[2] DEEKS A J, RANDOLPH M F. Axisymmetric time-domain transmitting boundaries[J]. Journal of Engineering Mechanics, 1994, 120(1): 25-42.

[3] 谷音, 劉晶波, 杜義欣. 三維一致粘彈性人工邊界及等效粘彈性邊界單元[J]. 工程力學(xué), 2007, 24(12): 31-37.

[4] 趙建峰, 杜修力, 韓強(qiáng), 等. 外源波動(dòng)問(wèn)題數(shù)值模擬的一種實(shí)現(xiàn)方式[J]. 工程力學(xué), 2007, 24(4): 52-58.

[5] 劉晶波, 谷音, 杜義欣. 一致粘彈性人工邊界及粘彈性邊界單元[J]. 巖土工程學(xué)報(bào),2006, 28(9): 1 070-1 075.

(責(zé)任編輯： 洪江星)

Simulation of viscous-spring artificial boundaries on the basis of zero-thickness Goodman element

LUO Pei1, LIU Guoming1, WANG Wenjun1, 2, XU Chenkui1

(1. College of Civil Engineering， Fuzhou University， Fuzhou， Fujian 350116， China；2. Department of Architectural Engineering, Zhangzhou Vocational Technical College, Zhangzhou, Fujian 363000， China)

The stiffness and damping matrixes of 3D consistent artificial boundary are deduced through viscous-spring artificial boundary theory. Using thickness of the boundary element has little effect on the results, the zero-thickness Goodman element is developed as a new simulation method of consistent viscous-spring artificial boundary. The viscous-spring boundary element is simulated by the zero-thickness Goodman element, which is verified through the example of the half infinite space and uniform half space problems. The results show that this method is simple, practical, and can get enough accuracy. Compared with thickness of the boundary element, the zero-thickness Goodman element conforms to the actual situation, which the boundary element is no thick, have perfect theory and well adaptability to the possible irregular boundary conditions.

viscous-spring artificial boundary; zero-thickness Goodman element; artificial boundary element; dynamic calculation

10.7631/issn.1000-2243.2016.01.0104

1000-2243(2016)01-0104-06

2014-07-17

劉國(guó)明(1963-)， 博士， 教授， 主要從事水工結(jié)構(gòu)工程專業(yè)研究， lgm6379@163.com

福建省自然科學(xué)基金資助項(xiàng)目(2011J01309)

TV312

A