基于泰勒展式的混合階次流形方法

，

(三峽大學 土木與建筑學院，湖北 宜昌 443002)

1 研究背景

數(shù)值流形方法(NMM)是石根華博士1991年以拓撲流形和微分流形為基礎，吸收傳統(tǒng)有限元和非連續(xù)體運動學的優(yōu)點，采用有限覆蓋方法把連續(xù)和不連續(xù)變形力學統(tǒng)一于完整的數(shù)學表達式當中[1]。該方法的顯著特點在于單元形狀及其單元上的位移函數(shù)構(gòu)造比較靈活，在處理各種復雜邊界、復雜材料邊界和提高數(shù)值解精度方面優(yōu)勢尤為明顯。因此，該方法在巖土工程中得到了廣泛應用。例如：王水林等[2]利用NMM實現(xiàn)了裂紋擴展過程的數(shù)值模擬，王書法等[3]利用參變量變分原理求解了彈塑性問題，劉紅巖等[4]用數(shù)值流形法模擬了3種節(jié)理巖體在爆炸荷載下的沖擊破壞，等等。

在單元選擇上，筆者放棄傳統(tǒng)三角形單元而采用矩形單元。因為三角形單元具有網(wǎng)格生成復雜，小塊物理覆蓋形成數(shù)量多等缺點[5]，這對前后處理和數(shù)值計算帶來了許多不便。矩形單元不僅能很好地避免這些問題，而且具有更高的數(shù)值解精度[6]。在構(gòu)造高階位移函數(shù)時，國內(nèi)大部分學者采用完全多項式函數(shù)[7]，但是此方法使得升階后的各個廣義自由度完全喪失物理意義。為解決這個問題，筆者嘗試使用泰勒展開式(將高階位移函數(shù)看作是某點的泰勒展開)，再結(jié)合小變形假設下的位移與應變以及轉(zhuǎn)角之間的物理關系，建立了位移函數(shù)與位移、應變、轉(zhuǎn)角之間的函數(shù)關系，使得升階后的各個廣義自由度仍具有物理意義，為引入邊界條件帶來了便利。

在算例分析實現(xiàn)過程中，筆者在結(jié)構(gòu)求解區(qū)域的不同地方使用不同級數(shù)的泰勒展開式( 即不同階次的覆蓋位移場函數(shù)),實現(xiàn)了有限元法與解析法的完美結(jié)合。此種方法不僅大大減少了計算量，還提高了計算精度。

2 矩形流形單元的構(gòu)造及其位移函數(shù)

圖1 重疊網(wǎng)格圖

對于無縫連續(xù)體的數(shù)值計算，流形方法可以采用重疊的數(shù)學網(wǎng)格和物理網(wǎng)格[7]，如圖1所示。圖1中，矩形流形單元1452是由4個片Ω1,Ω4,Ω5,Ω2相互重疊形成的，其中，Ω1是對應于節(jié)點1的片(patch)，由矩形1452構(gòu)成；Ω4是對應于節(jié)點4的片，由共點于節(jié)點4的2個矩形1452和4785構(gòu)成；Ω5是對應于節(jié)點5的片，由共點于節(jié)點5的4個矩形1452，4785，5896和2563所構(gòu)成；Ω2是對應于節(jié)點2的片，由共點于節(jié)點2的2個矩形1452和2563構(gòu)成。在各個片上定義獨立的覆蓋位移場函數(shù)，再通過權(quán)函數(shù)將各個局部場函數(shù)連接在一起，這樣就構(gòu)造出求解域上的總體位移函數(shù)。為了便于形成系統(tǒng)總綱，需引入流形單元，以流形單元1452為例，它是Ω1,Ω4,Ω5,Ω2的公共部分。令w1,w4,w5,w2分別為Ω1,Ω4,Ω5,Ω2上的場函數(shù)，φ1,φ4,φ5,φ2為連接各個片的形函數(shù)，則總體逼近函數(shù)為

(1)

2.1 覆蓋位移場函數(shù)

目前，常采用Lagrange分片插值函數(shù)來構(gòu)造流形單元的局部位移場函數(shù)。由于場函數(shù)定義在獨立的物理覆蓋區(qū)域上，這就導致了不同節(jié)點可選取不同階次的插值函數(shù)來進行插值。這里主要討論零階(場函數(shù)為常數(shù))和一階節(jié)點(場函數(shù)為一次多項式)。對于零階節(jié)點，其數(shù)值求解精度與傳統(tǒng)有限元相同。故此，本文多采用一階節(jié)點。對于一階節(jié)點，水平位移ui(x,y)可以看作是水平位移函數(shù)u(x,y)在(xi,yi)處的一階泰勒展開，即

式中：

同理，對于垂直位移vi(x,y)，有

式中：

基于小變形假設，

這里ω表示過點(xi,yi)的任意矢量繞z軸的轉(zhuǎn)角。因此位移矢量Ui(x,y)=(ui(xi,yi),vi(x,y))T，可以寫成

Ui=Tidi, Ui∈Ωi。

(2)

這樣節(jié)點自由度便有了明確的物理意義，便于引入邊界條件。

對于高于一階的節(jié)點，以二階為例，其水平位移u(x,y)在(xi,yi)的二階泰勒展式為

由上述推導可知：當覆蓋位移函數(shù)階次大于一階時，使用泰勒展式并不能使升階后的各個廣義自由度都具有明確的物理意義。

2.2 權(quán)函數(shù)

流形方法里的權(quán)函數(shù)與有限元法里的插值形函數(shù)定義原則一致。圖2為四節(jié)點矩形單元。

圖2 四節(jié)點矩形單元

矩形的邊長分別為2a,2b。矩形的兩邊分別與x軸，y軸平行。為方便運算，引入局部坐標系x′y′局部坐標系的原點取在矩形形心，x′軸和y′軸分別和x,y軸平行。其坐標變換關系為

x=x0+ax′ ,y=y0+by′ 。

(3)

式中：

其中(xi,yi)(i=1,2,3,4)是節(jié)點i的整體坐標。

(4)

將式(3)代入式(4)中可得整體坐標系中的形函數(shù)為：

(5)

則水平位移函數(shù)u、垂直位移函數(shù)v和局部逼近位移函數(shù)U分別為：

(6)

式中，單元整體形函數(shù)S和單元整體自由度D分別為:

S=[φ1T1φ2T2φ3T3φ4T4];

D=[d1d2d3d4]T。

2.3 單元應變矩陣

將式(6)的單元位移函數(shù)代入單元的物理和幾何方程,可得到單元應變矩陣

式中：B稱為應變矩陣；Ld是平面問題的微分算子，且

對于全一階節(jié)點

3 流形單元分析

為計算和編程的方便，本文采用改進的罰函數(shù)法與廣義節(jié)點法相結(jié)合的方式(參看本文第3.2節(jié))來施加邊界約束條件。有約束的修正泛函為

根據(jù)變分原理，令δπ(u)=0得到平衡方程Kp=q。其中：p是所有節(jié)點總的全體自由度；K是總體剛度矩陣，由各個單元剛度矩陣Ke組裝形成，即

3.1 數(shù)值積分

計算單元剛度矩陣和等效荷載列陣時，必然會涉及到積分運算。隨著單元覆蓋位移函數(shù)階次的升高，單元剛度矩陣和等效荷載列陣的被積函數(shù)的方次也隨著增加[8]，數(shù)值積分的階次必將相應升高，這就要求具有很高精度的積分方案。本文采取高斯積分的方式來進行數(shù)值運算。

3.2 邊界處理

在施加位移約束時，有的只需要約束邊界上的幾個固定點，有的則需要對整條邊施加約束。對固定點的約束，筆者采用改進的罰函數(shù)法[9]，即將傳統(tǒng)的單一雙向約束改為更加靈活的雙向約束(在2個約束方向上分別設置互不相干的剛性彈簧)。而對固定邊的約束，筆者采用文獻[10]中提到的方法：將固定邊界的約束處理轉(zhuǎn)化為對廣義節(jié)點的約束處理，通過對廣義節(jié)點的處理來實現(xiàn)對固定邊界約束條件的滿足。對約束處理一般化的敘述為：當廣義節(jié)點受雙向固定約束時，取消該廣義節(jié)點覆蓋函數(shù)在單元位移函數(shù)中的組成部分；當廣義節(jié)點只受x或y向固定約束時，在該節(jié)點上使用只包含自由方向位移的覆蓋函數(shù)。由此得到流形單元的位移函數(shù)來進行其他公式(剛度矩陣、荷矩陣、慣性力矩陣等)的推導，則可準確地對固定約束問題進行處理。

按照上述方法處理位移邊界，從理論上嚴格符合固定約束的物理意義，同時簡化了處理工作，并經(jīng)筆者大量算例證實是有效和準確的。

4 數(shù)值算例

4.1 懸臂梁

如圖3所示懸臂梁，在末端受到垂直向下的集中力P1=P2=20 N,長度L=40 cm,寬h=10 cm，按平面應力問題進行分析。材料的彈性模量E=3×104kN/cm2,泊松比μ=0.25，不計自重。

圖3 懸臂梁及其網(wǎng)格

懸臂梁平面問題的解析解[11]為

本文采用全零階、零階和一階混合、全一階3種近似覆蓋位移函數(shù)進行計算研究分析。圖4為y=5 cm時的位移曲線。

圖4 y=5 cm時y方向的位移對比

從圖4所示的沿懸臂軸線在y方向位移結(jié)果可以看出，采用零階覆蓋函數(shù)的流形方法，其位移與有限元法矩形單元結(jié)果相同，與解析解偏離甚遠；而采用零階和一階混合近似覆蓋函數(shù)的流形方法，其位移解的精度有所提高。1-0-0-0型NMM幾乎和有限元結(jié)果相同；1-1-1-0型NMM和全一階NMM結(jié)果很接近。荷載附近區(qū)域，1階節(jié)點布置越多，數(shù)值解越精確。且位移存在如下關系：

全一階>1-1-1-0>1-0-1-0>1-0-0-0>零階=有限元。

當采用全一階覆蓋位移函數(shù)時，整個求解域的位移解精度有了明顯提高，基本上落在解析解上。

當采用混合階次流形方法時，其數(shù)值解的精度是由節(jié)點的布置、1階節(jié)點的數(shù)量和荷載附近的一階節(jié)點數(shù)這3方面共同決定的。

圖5 純彎曲梁及其網(wǎng)格

4.2 純彎曲梁

純彎曲梁的位移函數(shù)的解析解[12]為

如圖5所示純彎曲梁，右邊為計算模型及單元網(wǎng)格，M為左右兩端施加的線性分布力矩，為10 MPa，長8 m，寬2 m,，E=200 GPa，μ=0.25。由于對稱性和反對稱性，僅需對梁的1/4進行計算。根據(jù)對稱性和反對稱性，y軸和x軸上，x位移均為0，另外為消除結(jié)構(gòu)在y方向上的剛體位移，可令原點O處y方向上的位移為0(假設O點y方向真實位移為v0，則所有節(jié)點y方向的真實位移應為v+v0。其中，根據(jù)解析公式算得v0=4×10-4m),分別采用全零階、零階和一階混合、全一階3種近似覆蓋位移函數(shù)進行計算研究分析。

表1 y=1 m時，梁的右上部分y方向位移解

根據(jù)對稱性可推算出y=1 m時，上半部分梁在y方向位移對比如圖6所示。

圖6 y=1 m時y方向位移對比

根據(jù)圖表可得出：全1階NMM數(shù)值解幾乎全部落在解析解曲線上，零階NMM，1-0-0-0型NMM，1-1-0-0型NMM，1-1-1-0型NMM也與解析解相當接近。

5 結(jié) 語

數(shù)值流形方法是把連續(xù)和不連續(xù)變形力學統(tǒng)一于完整的數(shù)學表達式中的更一般化的數(shù)值解析方法。它不僅能計算不連續(xù)面的位錯、滑動、開裂和旋轉(zhuǎn)等大位移的靜力和動力學問題，對于連續(xù)體的彈性靜力學分析，也優(yōu)于有限元法，具有編程和前后處理簡單、數(shù)值解精度高等特點。但其理論體系和分析方法從提出到實際應用僅有20多年時間，很多方面還不成熟，仍在不斷發(fā)展之中。

對于無縫連續(xù)體的數(shù)值計算，本文(只使用1套網(wǎng)格，即數(shù)學網(wǎng)格與物理網(wǎng)格重合)直接采用矩形格子作為數(shù)學網(wǎng)格，相比三角形網(wǎng)格，減少了單元數(shù)和節(jié)點數(shù)，簡化了前處理和編程工作；并運用泰勒展式來提高位移階次，使得廣義自由度具有明確的物理意義，方便引入邊界條件，而實際處理邊界時采用改進的罰函數(shù)和廣義節(jié)點法，從理論上嚴格符合固定約束的物理意義。在滿足單元協(xié)調(diào)性的前提下，同時也提高了數(shù)值解的精度。以上算例表明，此種方法能很好地解決部分材料變形分析問題，具有一定的實用和研究價值。

參考文獻：

[1] SHI Gen-hua.Discontinuous Deformation Analysis and Numerical Models of Static and Dynamic Block System[D].Berkeley, USA: Department of Civil Engineering, University of California, 1989.

[2] 王水林，葛修潤.流形元方法在模擬裂紋擴展中的應用[J].巖石力學與工程學報,1997,16(5): 405-420.(WANG Shui-lin, GE Xiu-run.Application of Manifold Method in Simulating Crack Propagation [J].Chinese Journal of Rock Mechanics and Engineering, 1997， 16(5)：405-420.(in Chinese))

[3] 王書法, 朱維申, 李術才, 等.巖體彈塑性分析的數(shù)值流形方法[J].巖石力學與工程學報,2002, 21(6):900-904.(WANG Shu-fa, ZHU Wei-shen, LI Shu-cai,etal.Numerical Manifold Method of Elastic-plastic Analysis for Rock Mass[J].Chinese Journal of Rock Mechanics and Engineering, 2002，21(6):900-904.(in Chinese))

[4] 劉紅巖, 王貴和.節(jié)理巖體沖擊破壞的數(shù)值流形方法模擬[J].巖土力學,2009,31(11):3523-3527.(LIU Hong-yan, WANG Gui-he.Simulation of Impact Failure of Jointed Rock Mass by Numerical Manifold Method[J].Rock and Soil Mechanics, 2009,31(11): 3523-3527.(in Chinese))

[5] 宋俊生，大西有三.高階四邊形單元的流形方法[J].巖石力學與工程學報,2003,22(6): 932-936.(SONG Jun-sheng, OHNISHI YUZO.High-order Rectangular Element of Manifold Method[J].Chinese Journal of Rock Mechanics and Engineering, 2003, 22(6): 932-936.(in Chinese))

[6] 蔡永昌，廖林燦，張湘?zhèn)?高精度四節(jié)點四邊形流形單元[J].應用力學學報,2001,(2):75-80.(CAI Yong-chang, LIAO Lin-can, ZHANG Xiang-wei.A High Precision Quadrnodel-quadrangle Manifold Element[J].Journal of Applied Mechanics, 2001, (2):75-80.(in Chinese))

[7] 石根華，裴覺民.數(shù)值流形方法與非連續(xù)變形分析[M].北京：清華大學出版社，1997.(SHI Gen-hua, PEI Jue-min.Numerical Manifold Method and Discontinuous Deformation Analysis[M].Beijing: Tsinghua University Press, 1997.(in Chinese))

[8] 李樹忱, 李術才, 張京偉, 等.數(shù)值流形方法的數(shù)學推導及其應用[J].工程力學,2007, 24(6):36-42.(LI Shu-chen, LI Shu-cai, ZHANG Jing-wei,etal.Numerical Manifold Method from Mathematical Theory and Its Application[J].Journal of Engineering Mechanics, 2007, 24(6): 36-42.(in Chinese))

[9] 王芝銀, 李云鵬.在數(shù)值流形方法中的幾點改進[J].巖土工程學報,1998,20(6):33-36.(WANG Zhi-yin, LI Yun-peng.Some Improvements in the Numerical Manifold Method [J].Chinese Journal of Geotechnical Engineering, 1998, 20(6): 33-36.(in Chinese))

[10] 朱愛軍, 鄧安福, 曾祥勇,等.巖體工程數(shù)值流形方法的固定邊界約束處理方法[J].巖石力學與工程學報, 2005,24(19): 3582-3587.(ZHU Ai-jun, DENG An-fu, ZENG Xiang-yong,etal.Constraints Treatment of Fixed Boundary in Numerical Manifold Method[J].Chinese Journal of Rock Mechanics and Engineering, 2005, 24 (19): 3582-3587.(in Chinese))

[11] 徐芝綸.彈性力學簡明教程(第3版)[M].北京: 高等教育出版社, 2002.(XU Zhi-lun.Simple Tutorial of Elasticity(Third Edition)[M].Beijing: Higher Education Press, 2002.(in Chinese))

[12] 吳家龍．彈性力學[M]．北京：高等教育出版社，2001.(WU Jia-long.Elastic Mechanics [M].Beijing: Higher Education Press, 2001.(in Chinese))