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

耿長春，汪基偉

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

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

耿長春1，汪基偉2

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

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

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

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

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

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

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

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

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

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

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

(1)

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

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

細單元a的結(jié)點位移為：

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

(2)

其中：下標為1～n的變量為細單元上實結(jié)點的位移；下標為p*、s*的變量為虛結(jié)點的位移。細單元中的虛結(jié)點位

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

(3)

(4)

因此細單元結(jié)點位移可表示為：

(5)

式(5)可簡寫成：

ae=Ta′e

(6)

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

(7)

(8)

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

(9)

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

細單元的應變可表示為：

ε=Bae=BTa′e

(10)

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

細單元的應力可表示為：

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

(11)

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

作用在細單元體積上的作用力為：

(12)

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

(13)

(14)

細單元的等效結(jié)點荷載為：

(15)

(16)

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

Re=∫veBTDBdvae=∫veBTDBTdva′e

(17)

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

(18)

(19)

式中：k′e為細單元剛度矩陣；m′e為細單元質(zhì)量矩陣；c′e為細單元阻尼矩陣。

(20)

式(20)中的細單元剛度矩陣與細單元質(zhì)量矩陣的計算公式可轉(zhuǎn)化為：

(21)

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

在等參單元中任一點的坐標可表示為：

(22)

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

(23)

式中：J為Jacobian矩陣，對整體坐標中的一點P(x,y,z)，其局部坐標(ξ,η,ζ)應符合下式：

(24)

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

(25)

(26)

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

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

3 應用實例

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

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

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

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

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

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

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

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

表1 自振頻率及振型對比

Y向(Y- direction)

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

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

4 結(jié)束語

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

[1] 李同春,李淼,溫召旺,等.局部非協(xié)調(diào)網(wǎng)格在高拱壩應力分析中的應用[J].河海大學學報:自然科學版,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] 王愛民,王勖成.有限元計算中疏密網(wǎng)格間過渡單元的構(gòu)造[J].清華大學學報:自然科學版,1999,39(8):101-103

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

[5] 朱以文,蔡元奇,李偉.等參元逆變換算法在滲流位移耦合場分析中的應用[J].計算力學學報,2002,(2) :233-235

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

[7] 韋未,李同春,牛志偉,等.局部網(wǎng)格加密技術(shù)在混凝土裂縫擴展模擬中的應用[J].華南農(nóng)業(yè)大學學報,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] 賀向麗,李同春.基于逐步擴大網(wǎng)格法的飽和無限地基動力分析[J].工程力學,2010,27(4):149-184

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

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

TV313

A

[定稿日期]2014-08-27