耿長(zhǎng)春，汪基偉

(1.中交第四航務(wù)工程勘察設(shè)計(jì)院有限公司,廣東廣州 510000; 2.河海大學(xué),江蘇南京 210000)

非協(xié)調(diào)網(wǎng)格在結(jié)構(gòu)動(dòng)力分析中的應(yīng)用

耿長(zhǎng)春1，汪基偉2

(1.中交第四航務(wù)工程勘察設(shè)計(jì)院有限公司,廣東廣州 510000; 2.河海大學(xué),江蘇南京 210000)

水利工程中的結(jié)構(gòu)往往體形巨大且復(fù)雜，采用非協(xié)調(diào)網(wǎng)格技術(shù)可以減少結(jié)點(diǎn)數(shù)量，降低網(wǎng)格剖分難度。文章采用形函數(shù)插值法推導(dǎo)了非協(xié)調(diào)網(wǎng)格過(guò)渡單元的動(dòng)力平衡方程，給出了非協(xié)調(diào)網(wǎng)格過(guò)渡單元的剛度矩陣、質(zhì)量矩陣的計(jì)算公式。對(duì)某一出水池建立了一致網(wǎng)格和非協(xié)調(diào)網(wǎng)格兩種模型，用自編程序?qū)煞N模型進(jìn)行了動(dòng)力分析。通過(guò)對(duì)比分析兩種模型的計(jì)算結(jié)果，證實(shí)非協(xié)調(diào)網(wǎng)格技術(shù)應(yīng)用于此類(lèi)結(jié)構(gòu)是可行的，高效的。

有限元； 非協(xié)調(diào)網(wǎng)格； 自振特性； 動(dòng)力響應(yīng)； 剛度矩陣； 質(zhì)量矩陣

1 非協(xié)調(diào)網(wǎng)格的應(yīng)用

建立合理的有限元模型是混凝土結(jié)構(gòu)有限元分析的首要任務(wù)。在建立有限元模型時(shí)常會(huì)遇到下列問(wèn)題：(1)水工結(jié)構(gòu)(如壩、水閘等)有限元計(jì)算時(shí)，通常取地基在各個(gè)方向上的尺寸為建筑物高度的一倍以上。若采用一致網(wǎng)格，地基在水平方向的網(wǎng)格與上部結(jié)構(gòu)相同，這樣有可能使地基的單元數(shù)超過(guò)上部結(jié)構(gòu)。由于模型總結(jié)點(diǎn)數(shù)受到計(jì)算機(jī)性能的限制，對(duì)于大型結(jié)構(gòu)，不得不采用較粗的網(wǎng)格而降低計(jì)算精度。(2)對(duì)于體型比較復(fù)雜的結(jié)構(gòu)，采用一致網(wǎng)格剖分會(huì)有困難，如水電站蝸殼中的導(dǎo)葉與坐環(huán)。

采用非協(xié)調(diào)網(wǎng)格技術(shù)可以在一定程度上解決上述問(wèn)題。所謂非協(xié)調(diào)網(wǎng)格技術(shù)就是對(duì)結(jié)構(gòu)進(jìn)行網(wǎng)格剖分時(shí)采用了多于一種的網(wǎng)格尺度的單元。在不協(xié)調(diào)網(wǎng)格界面上，以網(wǎng)格尺度較大單元所包含的界面上的結(jié)點(diǎn)作為基本結(jié)點(diǎn)，而界面上與網(wǎng)格尺度較小單元相應(yīng)的結(jié)點(diǎn)作為從結(jié)點(diǎn)，建立從結(jié)點(diǎn)與基本結(jié)點(diǎn)的線性插值關(guān)系，從而導(dǎo)出以基本結(jié)點(diǎn)位移作為求解變量的非協(xié)調(diào)網(wǎng)格協(xié)調(diào)位移解法[1]，這種連接兩種網(wǎng)格尺度單元的區(qū)域即為過(guò)渡單元。如此，可方便網(wǎng)格剖分，節(jié)略單元與結(jié)點(diǎn)，減小計(jì)算規(guī)模。

文獻(xiàn)[2]～文獻(xiàn)[8]采用構(gòu)造界面過(guò)渡單元、最小勢(shì)能原理等方法來(lái)解決非協(xié)調(diào)截面的協(xié)調(diào)性，驗(yàn)證了非協(xié)調(diào)網(wǎng)格應(yīng)用于結(jié)構(gòu)靜力分析的合理性。文獻(xiàn)[9]、文獻(xiàn)[10]利用非協(xié)調(diào)網(wǎng)格對(duì)二維無(wú)限地基和二維重力壩進(jìn)行了動(dòng)力特性計(jì)算。本文首先根據(jù)結(jié)點(diǎn)位移協(xié)調(diào)性推導(dǎo)出三維非協(xié)調(diào)網(wǎng)格過(guò)渡單元的剛度矩陣與質(zhì)量矩陣公式，并根據(jù)公式編制了可應(yīng)用于實(shí)際工程中的三維動(dòng)力分析程序。

2 非協(xié)調(diào)網(wǎng)格過(guò)渡單元的剛度矩陣與質(zhì)量矩陣公式推導(dǎo)及程序?qū)崿F(xiàn)

為推導(dǎo)出非協(xié)調(diào)網(wǎng)格過(guò)渡單元的剛度矩陣與質(zhì)量矩陣，本文先列出空間8～20結(jié)點(diǎn)等參單元的剛度矩陣與質(zhì)量矩陣計(jì)算公式。圖1所示為空間8～20結(jié)點(diǎn)等參單元，該單元1～8結(jié)點(diǎn)為角結(jié)點(diǎn)，9～20結(jié)點(diǎn)為中間結(jié)點(diǎn)，任一個(gè)中間結(jié)點(diǎn)均可刪去。

圖1 8～20結(jié)點(diǎn)空間等參單元

空間8～20結(jié)點(diǎn)等參單元的剛度矩陣與單元質(zhì)量矩陣計(jì)算公式為：

(1)

下面以圖2所示過(guò)渡單元為例，來(lái)推導(dǎo)過(guò)渡單元的剛度矩陣和質(zhì)量矩陣。在圖2中，塊體單元a、p和s均為8結(jié)點(diǎn)等參單元，單元a網(wǎng)格尺寸小，稱(chēng)為細(xì)單元；單元p和s網(wǎng)格尺寸大，稱(chēng)為粗單元。細(xì)單元a的有些結(jié)點(diǎn)(最多4個(gè))不與粗單元p和s等單元的結(jié)點(diǎn)相連，而是交與粗單元p和s等單元的表面E上，定義這些結(jié)點(diǎn)為虛結(jié)點(diǎn)(其中虛結(jié)點(diǎn)在單元p上時(shí)記為p*，在單元s上時(shí)記為s*，以此類(lèi)推)，其余結(jié)點(diǎn)為實(shí)結(jié)點(diǎn)。虛結(jié)點(diǎn)的物理量(位移、速度與加速度)不能作為運(yùn)動(dòng)方程中的未知量，只能由與其相交的粗單元的結(jié)點(diǎn)物理量通過(guò)形函數(shù)插值得到。為了簡(jiǎn)單明了地表達(dá)公式，下文的公式中只具體寫(xiě)出有關(guān)虛結(jié)點(diǎn)p*和s*的部分，有關(guān)虛結(jié)點(diǎn)q*和r*的部分由“…”代替。

圖2 非協(xié)調(diào)網(wǎng)格過(guò)渡單元示意

細(xì)單元a的結(jié)點(diǎn)位移為：

ae=[u1v1w1…up*vp*wp*…us*vs*ws*…wn]

(2)

其中：下標(biāo)為1～n的變量為細(xì)單元上實(shí)結(jié)點(diǎn)的位移；下標(biāo)為p*、s*的變量為虛結(jié)點(diǎn)的位移。細(xì)單元中的虛結(jié)點(diǎn)位

移可由與其相交的粗單元結(jié)點(diǎn)位移通過(guò)形函數(shù)插值得到：

(3)

(4)

因此細(xì)單元結(jié)點(diǎn)位移可表示為：

(5)

式(5)可簡(jiǎn)寫(xiě)成：

ae=Ta′e

(6)

其中：a′e和T分別為：

(7)

(8)

細(xì)單元內(nèi)任一點(diǎn)的位移為：

(9)

式中:N為細(xì)單元的形函數(shù)矩陣，表達(dá)式如式(4)所示。

細(xì)單元的應(yīng)變可表示為：

ε=Bae=BTa′e

(10)

式中:B為應(yīng)變轉(zhuǎn)換矩陣。

細(xì)單元的應(yīng)力可表示為：

σ=DBae=DBTa′e=[S1S2…Sn]Ta′e

(11)

式中：Si為應(yīng)力轉(zhuǎn)換矩陣，Si=DBi；D為彈性矩陣。

作用在細(xì)單元體積上的作用力為：

(12)

將式(9)代入式(12)，得到：

(13)

(14)

細(xì)單元的等效結(jié)點(diǎn)荷載為：

(15)

(16)

由虛功原理可得平衡方程：

Re=∫veBTDBdvae=∫veBTDBTdva′e

(17)

將式(17)代入式(16)可得：

(18)

(19)

式中：k′e為細(xì)單元?jiǎng)偠染仃?；m′e為細(xì)單元質(zhì)量矩陣；c′e為細(xì)單元阻尼矩陣。

(20)

式(20)中的細(xì)單元?jiǎng)偠染仃嚺c細(xì)單元質(zhì)量矩陣的計(jì)算公式可轉(zhuǎn)化為：

(21)

比較式(1)和式(21)可知，非協(xié)調(diào)網(wǎng)格過(guò)渡單元的單元?jiǎng)偠染仃噆′e是在一般單元的單元?jiǎng)偠染仃噆e的兩邊分別乘以轉(zhuǎn)換矩陣TT和T；非協(xié)調(diào)網(wǎng)格過(guò)渡單元的質(zhì)量矩陣m′e也是在一般單元的單元質(zhì)量矩陣me的兩邊分別乘以轉(zhuǎn)換矩陣TT和T。此時(shí)得到的單元質(zhì)量矩陣也是一致質(zhì)量矩陣。

在等參單元中任一點(diǎn)的坐標(biāo)可表示為：

(22)

式中：N(ξ,η,ζ)為三維等參單元的形函數(shù)矩陣；X，Y，Z為單元結(jié)點(diǎn)坐標(biāo)向量。對(duì)式(22)微分可得：

(23)

式中：J為Jacobian矩陣，對(duì)整體坐標(biāo)中的一點(diǎn)P(x,y,z)，其局部坐標(biāo)(ξ,η,ζ)應(yīng)符合下式：

(24)

可用牛頓迭代法求解式(24)，迭代格式為：

(25)

(26)

式中：Jn=J(ξn,ηn,ζn)；Nn=N(ξn,ηn,ζn)。

用牛頓-拉斐遜法求解非線性方程組時(shí)，在真實(shí)解附近具有二階收斂速度，可較快地將整體坐標(biāo)轉(zhuǎn)換為其相應(yīng)的局部坐標(biāo)。

3 應(yīng)用實(shí)例

本文利用上述過(guò)渡單元的公式編制了空間非協(xié)調(diào)網(wǎng)格動(dòng)力分析程序，對(duì)某一出水池建立了一致網(wǎng)格和非協(xié)調(diào)網(wǎng)格兩種模型，并用自編程序?qū)υ摻Y(jié)構(gòu)進(jìn)行動(dòng)力計(jì)算。一致網(wǎng)格模型如圖3所示，整個(gè)模型共劃分了65 257個(gè)結(jié)點(diǎn)，52 956個(gè)單元。非協(xié)調(diào)網(wǎng)格模型如圖4～圖6所示，該模型共劃分了41 619個(gè)結(jié)點(diǎn)，32 320個(gè)單元，其中虛結(jié)點(diǎn)有1 840個(gè)，比一致網(wǎng)格模型節(jié)省了36%的結(jié)點(diǎn)數(shù)目。

圖3 一致網(wǎng)格計(jì)算模型

圖4 非協(xié)調(diào)網(wǎng)格計(jì)算模型

圖5 非協(xié)調(diào)網(wǎng)格橫向剖面示意

圖6 非協(xié)調(diào)網(wǎng)格縱向剖面示意

表1給出了兩種模型下結(jié)構(gòu)前九階自振頻率及振型描述。

由表1可以看出，兩種有限元模型得到的結(jié)構(gòu)自振頻率值非常接近，最大相差5.03%，振型相同，這說(shuō)明了非協(xié)調(diào)網(wǎng)格技術(shù)應(yīng)用于求解這類(lèi)工程自振特性是合理可行的。

本文選用El-centro(1940,NS) 地震波，最大加速度峰值為1.962 m/s2，時(shí)間歷時(shí)20 s，共有1 000個(gè)時(shí)刻，對(duì)兩個(gè)模型進(jìn)行了時(shí)程分析，一致網(wǎng)格模型費(fèi)時(shí)50 h，而非協(xié)調(diào)網(wǎng)格模型只用了26 h。在出水池結(jié)構(gòu)上取一個(gè)點(diǎn)來(lái)分析結(jié)構(gòu)在橫向地震波作用下的位移響應(yīng)。為節(jié)約篇幅，圖7列出兩種模型下該點(diǎn)在Y向和Z向的位移時(shí)程曲線。

表1 自振頻率及振型對(duì)比

Y向(Y- direction)

Z向(Z- direction)圖7 兩種模型下該點(diǎn)的位移時(shí)程曲線

由圖7可見(jiàn)，兩種模型所得的A兩點(diǎn)在Y向、Z向的位移的最大誤差為6.26 %，其中最大值的誤差只有3.70 %，說(shuō)明非協(xié)調(diào)網(wǎng)格技術(shù)應(yīng)用于這類(lèi)工程有足夠的精度，滿(mǎn)足工程精度要求。

4 結(jié)束語(yǔ)

在滿(mǎn)足工程精度要求的情況下，在大型工程動(dòng)力分析中應(yīng)用非協(xié)調(diào)網(wǎng)格技術(shù)可大為減少結(jié)點(diǎn)數(shù)目，極大提高計(jì)算效率，節(jié)省計(jì)算時(shí)間，該技術(shù)可廣泛地應(yīng)用到實(shí)際工程中。

[1] 李同春,李淼,溫召旺,等.局部非協(xié)調(diào)網(wǎng)格在高拱壩應(yīng)力分析中的應(yīng)用[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版,2003,31(1):42～45

[2] Zienkiew iczO C. The Finite ElementM Ethod[J].Third edition.London:MCGRAM-HLL Book Company Limited,1977,158-160

[3] 王愛(ài)民,王勖成.有限元計(jì)算中疏密網(wǎng)格間過(guò)渡單元的構(gòu)造[J].清華大學(xué)學(xué)報(bào):自然科學(xué)版,1999,39(8):101-103

[4] 朱以文,徐晗,蔡元奇.交界面非協(xié)調(diào)網(wǎng)格連結(jié)的新約束方程法[J].武漢大學(xué)學(xué)報(bào):工學(xué)版,2003,36(6):47-69

[5] 朱以文,蔡元奇,李偉.等參元逆變換算法在滲流位移耦合場(chǎng)分析中的應(yīng)用[J].計(jì)算力學(xué)學(xué)報(bào),2002,(2) :233-235

[6] 強(qiáng)天馳,寇曉東,周維垣.三維有限元網(wǎng)格加密界面協(xié)調(diào)方法及在大壩開(kāi)裂分析中的應(yīng)用[J].巖石力學(xué)與工程學(xué)報(bào),2000,19(5):562-566

[7] 韋未,李同春,牛志偉,等.局部網(wǎng)格加密技術(shù)在混凝土裂縫擴(kuò)展模擬中的應(yīng)用[J].華南農(nóng)業(yè)大學(xué)學(xué)報(bào),2007,28(4):112～116

[8] Carpnteria.Valentes, Ferrarag. Experimental and Numerical Fracture Modeling of a Gravitydam[J].BAZANT,ZP.Fracture Mechanics of Concrete Structures. New York: Elsevier, 1992: 351-360

[9] 賀向麗,李同春.基于逐步擴(kuò)大網(wǎng)格法的飽和無(wú)限地基動(dòng)力分析[J].工程力學(xué),2010,27(4):149-184

[10] 陳禮平,劉國(guó)明.非協(xié)調(diào)網(wǎng)格在重力壩動(dòng)力特性有限元計(jì)算中的應(yīng)用[J].福州大學(xué)學(xué)報(bào):自然科學(xué)版,2005,33(1):84-88

耿長(zhǎng)春(1988～)，男，碩士，助理工程師，主要從事計(jì)算方法研究和結(jié)構(gòu)設(shè)計(jì)；汪基偉(1962～)，男，教授，博士，博導(dǎo)，主要從事鋼筋混凝土結(jié)構(gòu)限裂配筋研究和結(jié)構(gòu)靜動(dòng)力分析。

TV313

A

[定稿日期]2014-08-27