徐 磊 葉志才 任青文 杜小凱

(1.河海大學(xué)水利水電工程學(xué)院, 南京 210098;2.河海大學(xué)土木工程學(xué)院, 南京 210098;3.宿遷市水務(wù)局, 江蘇宿遷 223800;4.中國水電工程顧問集團(tuán)公司, 北京 100120)

在對(duì)結(jié)構(gòu)-地基動(dòng)力相互作用問題的計(jì)算分析中,如何考慮無限地基的輻射阻尼效應(yīng)是核心問題之一[1].若人工截取有限區(qū)域地基來模擬無限地基,則會(huì)使得散射波在人工截取的邊界上產(chǎn)生反射而導(dǎo)致模擬失真.目前,在截取的有限域邊界上設(shè)置合理的動(dòng)力人工邊界條件是解決這一問題的有效方法[2].近年來,相關(guān)學(xué)者已提出了多種動(dòng)力人工邊界條件,如粘性邊界、旁軸近似邊界、透射邊界和粘彈性邊界等[3].其中,粘彈性動(dòng)力人工邊界[4]以其精度較高、能模擬人工邊界外半無限介質(zhì)彈性恢復(fù)性能且有良好的頻率穩(wěn)定性而在結(jié)構(gòu)-地基動(dòng)力相互作用問題的有限元分析中得到了廣泛的應(yīng)用.

對(duì)于實(shí)際工程中的結(jié)構(gòu)-地基動(dòng)力相互作用有限元分析,一般均需采用三維計(jì)算模型,隨著工程仿真建模的日益精細(xì)化,有限元分析模型存在大量的邊界結(jié)點(diǎn),數(shù)量通?？梢赃_(dá)到數(shù)千甚至更多.當(dāng)采用粘彈性動(dòng)力人工邊界時(shí),由于各邊界結(jié)點(diǎn)與散射波源(結(jié)構(gòu))之間距離不等、各結(jié)點(diǎn)所控制的邊界面積不等以及各點(diǎn)所在區(qū)域的地基材料也會(huì)有所差別,因而,對(duì)于各邊界結(jié)點(diǎn),需逐點(diǎn)計(jì)算其彈簧系數(shù)及阻尼系數(shù)并施加彈簧單元及阻尼器單元,有限元前處理工作量浩繁.為了便于處理,通常采用一些近似方法來進(jìn)行粘彈性邊界的施加,如各邊界結(jié)點(diǎn)取散射波源至其所在人工邊界面的平均距離作為其與散射波源實(shí)際距離的近似等[5],這些近似勢必會(huì)對(duì)邊界條件的精度造成一定程度的影響.

為此,基于粘彈性動(dòng)力人工邊界理論大型商用有限元軟件ABAQUS,開發(fā)相關(guān)程序,實(shí)現(xiàn)了對(duì)粘彈性動(dòng)力人工邊界的精確自動(dòng)施加.算例分析驗(yàn)證了程序的合理性和有效性.此外,還將本文方法和通常的近似施加方法進(jìn)行了對(duì)比分析,結(jié)果表明,近似施加方法會(huì)引起較大計(jì)算誤差,本文所提出方法的施加簡便性、施加精度均優(yōu)于通常的近似施加方法.

1 粘彈性動(dòng)力人工邊界

粘彈性動(dòng)力人工邊界可以等效為在計(jì)算模型截?cái)噙吔缟线B續(xù)分布的并聯(lián)彈簧-阻尼器系統(tǒng),其中,彈簧系數(shù)Kb及阻尼系數(shù)Cb的計(jì)算公式如下:

式中,G、ρ分別為人工邊界處介質(zhì)的剪切模量、密度;R 為散射波源至人工邊界的距離;c 為人工邊界處介質(zhì)的波速,法向阻尼系數(shù)取為縱波速cp,切向阻尼系數(shù)取為剪切波速cs;α的取值取決于粘彈性動(dòng)力人工邊界的類型及設(shè)置方向,見表1.

表1 α的取值

粘彈性動(dòng)力人工邊界是一種連續(xù)分布的人工邊界條件.當(dāng)采用有限元法將計(jì)算區(qū)域離散化后,連續(xù)的人工邊界面也隨之離散化,此時(shí),在有限元計(jì)算模型中,為了實(shí)現(xiàn)粘彈性動(dòng)力人工邊界的施加,需要在邊界結(jié)點(diǎn)的法向和切向分別設(shè)置并聯(lián)的彈簧單元和阻尼器單元,相應(yīng)的彈簧系數(shù)K 和阻尼系數(shù)C 除了與式(1)有關(guān)外,還與邊界結(jié)點(diǎn)的等效控制面積相關(guān),計(jì)算表達(dá)式如下

式中,Gi、ρi、ci及Ai分別為與邊界結(jié)點(diǎn)相關(guān)聯(lián)單元的材料物理力學(xué)參數(shù)及面積.

對(duì)于實(shí)際工程問題,由于有限元計(jì)算分析模型的邊界結(jié)點(diǎn)數(shù)量眾多,且需逐點(diǎn)按式(2)計(jì)算其彈簧系數(shù)及阻尼系數(shù)并施加彈簧單元及阻尼器單元,有限元前處理工作十分困難,若采用近似處理又難免引入誤差.為此,在上述基礎(chǔ)上,基于ABAQ US 編制了用于粘彈性動(dòng)力人工邊界精確自動(dòng)施加的fortran 程序.

2 粘彈性動(dòng)力人工邊界精確自動(dòng)施加

近年來,作為國際上最先進(jìn)的大型通用商業(yè)非線性有限元分析軟件之一,A BAQ US 在結(jié)構(gòu)-地基動(dòng)力相互作用等問題中的應(yīng)用日益廣泛[6].ABAQUS 中的彈簧阻尼器單元也為施加粘彈性動(dòng)力人工邊界提供了方便.此外,用戶還可以通過二次開發(fā)接口[7]或編程讀寫input 計(jì)算文件等方式實(shí)現(xiàn)更為復(fù)雜的計(jì)算功能.為此,基于粘彈性動(dòng)力人工邊界理論以及ABAQ US 的input 計(jì)算文件的相關(guān)格式,編制了用以實(shí)現(xiàn)粘彈性動(dòng)力人工邊界精確自動(dòng)施加的fortran程序.該程序在獲取有限元計(jì)算模型相關(guān)信息的基礎(chǔ)上,直接生成用以施加粘彈性動(dòng)力人工邊界的命令文本,該文本完全符合input 計(jì)算文件的相關(guān)格式,將其放入相應(yīng)input 文件的指定位置即可完成基于ABAQ US 的粘彈性動(dòng)力人工邊界精確自動(dòng)施加.

為了實(shí)現(xiàn)粘彈性動(dòng)力人工邊界精確自動(dòng)施加,首先需要完成有限元計(jì)算模型的網(wǎng)格剖分、材料分區(qū)等,在此基礎(chǔ)上,獲取所需的有限元網(wǎng)格的結(jié)點(diǎn)坐標(biāo)、各單元的結(jié)點(diǎn)編號(hào)及其材料號(hào)、需要施加粘彈性動(dòng)力人工邊界的邊界結(jié)點(diǎn)編號(hào)、各邊界結(jié)點(diǎn)所在邊界面的法線方向以及散射波源位置等信息.

在上述基礎(chǔ)上,編制了相應(yīng)的fortran 程序,實(shí)現(xiàn)了基于ABAQUS 的粘彈性動(dòng)力人工邊界條件精確自動(dòng)施加,圖1 給出了相應(yīng)的程序流程圖.

圖1 程序流程圖

下面對(duì)程序中較為關(guān)鍵的問題進(jìn)行扼要的說明,首先是如何計(jì)算關(guān)聯(lián)單元邊界面(即位于截?cái)噙吔缟系膯卧?的面積,一般而言,與某一邊界結(jié)點(diǎn)相關(guān)聯(lián)的單元有多個(gè),對(duì)任一關(guān)聯(lián)單元,可首先依據(jù)單元各結(jié)點(diǎn)坐標(biāo)與邊界結(jié)點(diǎn)坐標(biāo)之間的關(guān)系,并結(jié)合邊界結(jié)點(diǎn)所在面的法向來確定在單元邊界面的所有結(jié)點(diǎn).其次是如何計(jì)算單元邊界面的面積,對(duì)于二維問題,邊界面的面積計(jì)算實(shí)質(zhì)上是計(jì)算兩點(diǎn)間線段的長度,易于解決;對(duì)于三維問題,單元邊界面形狀一般為三角形或四邊形,對(duì)于三角形,依其頂點(diǎn)坐標(biāo)容易計(jì)算其面積,而對(duì)于四邊形,當(dāng)只知其4 個(gè)頂點(diǎn)坐標(biāo)而不知其頂點(diǎn)順序時(shí),無法應(yīng)用通常的四邊形面積公式進(jìn)行求解,此時(shí),可先應(yīng)用海倫公式計(jì)算由此4 個(gè)頂點(diǎn)中任意3 個(gè)頂點(diǎn)所確定的4 個(gè)不同三角形的面積,則四邊形面積即為這4 個(gè)三角形面積之和的一半.需要說明的是,以上關(guān)于面積計(jì)算的說明中僅考慮不含中結(jié)點(diǎn)的單元,對(duì)于非此類單元,一般均可將其單元形態(tài)近似為不含中結(jié)點(diǎn)的單元來加以處理.

3 算例驗(yàn)證及對(duì)比分析

為了驗(yàn)證2 中所編程序的正確性和有效性,進(jìn)行了如下算例分析.

考慮半無限空間z ≥0,分析模型如圖2 所示,人工截取計(jì)算區(qū)域取為|x|≤10 m 、|y|≤10 m、|z|≤10 m.模型介質(zhì)的物理力學(xué)參數(shù)如下:彈性模量E =50 MPa;泊松比v=0.3;剪切波速cs=83.62 m/s;縱波速cp=156.48 m/s;密度ρ=2 750 kg/m3.為了便于對(duì)計(jì)算成果進(jìn)行分析,在分析模型上設(shè)置了2 個(gè)觀察點(diǎn)A、B,其坐標(biāo)分別為(0,0,10),(5,0,10).在自由面上沿z 方向加一集中面源(作用范圍為-5 m ≤x ≤5 m,-5 m ≤y ≤5 m,z=10 m)時(shí)變荷載f(t), f(t)的函數(shù)表達(dá)式如式(3)所示.

式中,0 ≤t ≤2 s.

圖2 算例分析模型

采用8 結(jié)點(diǎn)六面體單元進(jìn)行網(wǎng)格剖分,有限元模型如圖3 所示,網(wǎng)格尺寸Δx=Δy =Δz=1 m,散射波源取為荷載作用面的中心.在截?cái)噙吔缟?采用節(jié)2中編制的fortran 程序?qū)崿F(xiàn)粘彈性動(dòng)力人工邊界的精確自動(dòng)施加,記為C.B.;為了對(duì)比分析本文方法與通常近似方法之間的差別,還采用了后一方法對(duì)上述算例進(jìn)行了粘彈性邊界的施加,在施加過程中,各邊界結(jié)點(diǎn)取散射波源至其所在人工邊界面的平均距離作為其與散射波源實(shí)際距離的近似,記為A.B.;此外,還給出了遠(yuǎn)置粘彈性動(dòng)力人工邊界(本文采取在模型各個(gè)方向均擴(kuò)大20 倍范圍后設(shè)置人工邊界的方法)的有限元解作為精確解,記為E.B..

圖3 有限元計(jì)算網(wǎng)格圖

動(dòng)力時(shí)程分析的時(shí)間步長取為0.01 s,計(jì)算總時(shí)長為3 s.圖4 ～5 給出了觀測點(diǎn)A、B 在不同邊界條件下的z 向位移反應(yīng)時(shí)程曲線.從中可以看出,與精確解相比,本文所提出的粘彈性動(dòng)力人工邊界精確自動(dòng)施加方法(C.B.)與近似施加方法(A.B.)均可以給出與精確解(E.B.)總體趨勢一致的計(jì)算結(jié)果,但也可明顯發(fā)現(xiàn),C.B.方法所給出的計(jì)算結(jié)果與A.B.方法所給出的計(jì)算結(jié)果相比,與精確解更為一致,而A.B.方法所給出的計(jì)算結(jié)果與精確解相比有較大誤差,這是由于其在設(shè)置粘彈性動(dòng)力人工邊界時(shí)所采用的簡化近似處理而造成的.以上分析表明,本文所提出的粘彈性動(dòng)力人工邊界精確自動(dòng)施加方法在計(jì)算精度上相比通常的近似施加方法有明顯的優(yōu)勢,可以給出更為精確的計(jì)算結(jié)果,此外,由于本文所提出方法為自動(dòng)施加,其在施加速度上較通常近似施加方法也有明顯的優(yōu)勢,尤其是對(duì)于計(jì)算規(guī)模較大、邊界結(jié)點(diǎn)較多的結(jié)構(gòu)與地基動(dòng)力相互作用分析問題.

4 結(jié) 語

在結(jié)構(gòu)-地基動(dòng)力相互作用問題的有限元分析中,粘彈性動(dòng)力人工邊界得到了廣泛的應(yīng)用.但對(duì)于規(guī)模較大的復(fù)雜計(jì)算模型而言,粘彈性動(dòng)力人工邊界的施加存在著前處理工作量浩繁、施加精度不高等問題.為此, 本文基于粘彈性動(dòng)力人工邊界理論及ABAQ US 軟件平臺(tái),提出了粘彈性動(dòng)力人工邊界精確自動(dòng)施加方法,編制了相關(guān)程序并進(jìn)行了驗(yàn)證.算例將本文方法和通常的近似施加方法進(jìn)行了對(duì)比分析,結(jié)果表明,近似施加方法會(huì)引起較大計(jì)算誤差,本文所提出方法的施加簡便性、施加精度均優(yōu)于通常的近似施加方法.

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

[2] 劉晶波,王振宇,杜修力等.波動(dòng)問題中的三維時(shí)域粘彈性動(dòng)力人工邊界[ J] .工程力學(xué),2005,22(6):46-51.

[3] 廖振鵬.工程波動(dòng)理論導(dǎo)論[ M] .2 版.北京:科學(xué)出版社, 2002.

[4] Deeks A J, Randolph M F.Axisymmetric Time-domain Transmitting Boundaries[ J] .Journal of Engineering Mechanics, 1994,120(1):25-42.

[5] 張燎軍,張慧星, 王大勝等.黏彈性人工邊界在ADINA中的應(yīng)用[ J] .世界地震工程,2008,24(1):12-16.

[6] 莊海洋,陳國興,胡曉明.兩層雙柱島式地鐵車站結(jié)構(gòu)水平向非線性地震反應(yīng)分析[ J] .巖石力學(xué)與工程學(xué)報(bào),2006, 25(S1):3074-3079.

[7] 葉志才, 徐 磊,王 超.基于ABAQUS 的三維錨桿單元開發(fā)[ J] .三峽大學(xué)學(xué)報(bào):自然科學(xué)版, 2008, 30(5):29-32.