擴(kuò)展有限元中邊界條件的施加

冷 飛 張 然

(河海大學(xué) 土木與交通學(xué)院， 南京 210098)

?

擴(kuò)展有限元中邊界條件的施加

冷 飛 張 然

(河海大學(xué) 土木與交通學(xué)院， 南京 210098)

擴(kuò)展有限元施加邊界條件的方法與有限元不同，但尚無(wú)文獻(xiàn)進(jìn)行詳細(xì)的說(shuō)明．本文由有限元的控制方程和邊界條件的基本格式出發(fā)，說(shuō)明了擴(kuò)展有限元施加應(yīng)力和位移邊界的原理，以單個(gè)帶裂縫單元為例給出了Standard-XFEM和Shifted-XFEM模式下擴(kuò)展有限元施加邊界條件的具體實(shí)現(xiàn)過(guò)程．最終，通過(guò)一個(gè)受均布荷載的簡(jiǎn)支梁驗(yàn)證了本文介紹方法的正確性．結(jié)果表明：對(duì)于應(yīng)力邊界，應(yīng)由最小位能原理導(dǎo)出的公式計(jì)算等效結(jié)點(diǎn)荷載，不同位移模式下，由外部荷載計(jì)算等效結(jié)點(diǎn)荷載的公式相同，但得到的結(jié)點(diǎn)荷載向量不同．位移邊界則應(yīng)將不同的位移模式在約束處所滿足的方程帶入控制方程求解．雖然均采用上述原理施加邊界條件，但采用不同位移模式時(shí)，實(shí)現(xiàn)邊界條件的具體方法不盡相同．

擴(kuò)展有限單元法； 位移邊界條件； 應(yīng)力邊界條件； 非線性計(jì)算

0 引 言

擴(kuò)展有限單元法最初應(yīng)用于預(yù)設(shè)有裂縫、孔洞、夾雜的線彈性材料的計(jì)算，屬于線性計(jì)算范疇，隨著理論研究的發(fā)展，在非線性材料與結(jié)構(gòu)的計(jì)算中也有越來(lái)越多的應(yīng)用[1-2]．對(duì)于非線性計(jì)算，增量法是一種非常常見的迭代計(jì)算方法．這種方法要求在每一步迭代中均計(jì)算不平衡力，并將其作為結(jié)點(diǎn)荷載繼續(xù)計(jì)算．因此，在非線性計(jì)算中會(huì)遇到給帶裂縫單元施加荷載的問(wèn)題．當(dāng)加載區(qū)域存在不連續(xù)問(wèn)題時(shí)，即使是線性計(jì)算也會(huì)遇到上述問(wèn)題．此外，對(duì)于線性以及非線性材料的實(shí)際結(jié)構(gòu)，裂縫有可能會(huì)出現(xiàn)在約束附近，因而如何給帶裂縫單元施加約束條件也是實(shí)際結(jié)構(gòu)擴(kuò)展有限元分析會(huì)遇到的問(wèn)題之一[3-4]．

在有限元計(jì)算中，無(wú)論是施加荷載還是施加約束，均屬于控制方程的邊界問(wèn)題，即，對(duì)于計(jì)算域Ω，控制方程為

(1)

在其邊界Γ上，應(yīng)滿足

(2)

在擴(kuò)展有限元中，也應(yīng)按以上邊界條件處理帶裂縫單元的位移邊界或應(yīng)力邊界[5]．在實(shí)際計(jì)算中，對(duì)于位移邊界，應(yīng)根據(jù)不同的位移模式，對(duì)邊界處位移采取拉格朗日乘子法[6-8]或罰函數(shù)法[9-10]等具體方法進(jìn)行處理；對(duì)于應(yīng)力邊界，可利用最小位能原理得到等效結(jié)點(diǎn)荷載的計(jì)算公式．

以上方法與常規(guī)有限元處理位移與應(yīng)力邊界條件的方法完全相同，但實(shí)際結(jié)果卻不完全相同．由于沒(méi)有文獻(xiàn)或資料專門對(duì)此問(wèn)題進(jìn)行詳細(xì)解釋，部分研究者在涉及此類問(wèn)題時(shí)按常規(guī)有限元直接規(guī)定或假設(shè)應(yīng)如何施加邊界條件，如文獻(xiàn)[11]則規(guī)定外荷載只作用在常規(guī)自由度卻并未說(shuō)明原因，部分研究人員認(rèn)為某結(jié)點(diǎn)受到約束時(shí)該點(diǎn)常規(guī)自由度與富集自由度的位移均應(yīng)為0等．因而，需要對(duì)擴(kuò)展有限元如果施加邊界條件進(jìn)行明確的說(shuō)明．

本文將嚴(yán)格由邊界條件方程出發(fā)，給出擴(kuò)展有限元中standard XFEM模式與Shifted-XFEM模式施加位移邊界條件和應(yīng)力邊界條件的具體方法，并通過(guò)相應(yīng)的算例驗(yàn)證其正確性．

1 Standard XFEM模式下邊界條件的施加

1.1 Standard XFEM模式下的位移邊界條件

Standard XFEM的位移模式為

(3)

(4)

在某一點(diǎn)施加位移約束時(shí)，應(yīng)由式(4)計(jì)算得到該點(diǎn)滿足式(2)的邊界條件，并將其作為控制方程的約束條件．一般而言，該約束條件為關(guān)于u和a的表達(dá)式．

圖1為帶裂縫的富集單元，邊長(zhǎng)為2l，豎直方向裂縫經(jīng)過(guò)單元形心．結(jié)點(diǎn)1、結(jié)點(diǎn)2和結(jié)點(diǎn)3分別受到x向、雙向和y向約束．u1～u8分別表示各結(jié)點(diǎn)x方向和y方向的常規(guī)自由度，a1～a8分別表示1到4節(jié)各結(jié)點(diǎn)x方向和y方向的富集自由度．因此，單元的位移向量可寫為

(5)

圖1 開裂單元受結(jié)點(diǎn)集中荷載

對(duì)于結(jié)點(diǎn)1，由式(3)，其x向位移為u1+a1，因其受到支座的約束作用，應(yīng)有u1+a1=0，這就是控制方程應(yīng)滿足的位移邊界條件．同理也可得到其他結(jié)點(diǎn)的位移邊界條件．若采用拉格朗日乘子法處理該邊界條件，控制方程還應(yīng)滿足以下條件：

(6)

其中λ1，λ2，λ3和λ4為拉格朗日乘子．將上述4個(gè)位移邊界條件的表達(dá)式引入整體剛度矩陣K中后即可得到施加位移邊界條件后的新整體剛度矩陣K*．受版面所限，本文只列出將式(6)中第一個(gè)子式引入整體剛度矩陣后的新K*，即

如此，就可以求出單元的位移向量u．

由該例可看出，對(duì)于標(biāo)準(zhǔn)模式的擴(kuò)展有限元，雖然其控制方程與有限元相同，但當(dāng)帶裂縫單元受到位移約束時(shí)，不能像有限元一樣簡(jiǎn)單認(rèn)為被約束結(jié)點(diǎn)的常規(guī)自由度與富集自由位移均為0，而應(yīng)由位移插值模式得出對(duì)應(yīng)自由度u和a的關(guān)系，再耦合進(jìn)剛度矩陣．

1.2 Standard XFEM模式下的應(yīng)力邊界條件

按照最小位能原理[12]，易知擴(kuò)展有限元常規(guī)自由度上的結(jié)點(diǎn)荷載與有限元完全相同，富集自由度所對(duì)應(yīng)的結(jié)點(diǎn)荷載為

(7)

其中，b為體力，t為面力，F(xiàn)為集中力．

由上式可見，富集自由度結(jié)點(diǎn)荷載的計(jì)算方法與常規(guī)有限元是一致的．仍以圖1所示單元為例，在其1、4結(jié)點(diǎn)施加豎直向上的集中荷載，集中荷載大小為F．由有限元中結(jié)點(diǎn)荷載的計(jì)算方法及式(7)，可得該單元的等效節(jié)點(diǎn)荷載向量為

(8)

現(xiàn)假設(shè)泊松比μ=0，根據(jù)式(6)施加約束條件，由式(8)所示荷載向量可計(jì)算得到該單元內(nèi)的應(yīng)力分布σ為

(9)

如果用有限元法計(jì)算圖1所示單元，當(dāng)采用片狀裂縫模型時(shí)和分離裂縫模型時(shí)，計(jì)算得到的應(yīng)力分別為

(10)

(11)

觀察式(9)至式(11)所示的應(yīng)力分布可發(fā)現(xiàn)：采用片狀裂縫模型時(shí)，該單元在荷載方向?yàn)檩S拉構(gòu)件，因而單元應(yīng)力呈均勻分布，且只在荷載方向存在正應(yīng)力；采用分離模型時(shí)，裂縫兩側(cè)分別為兩個(gè)角點(diǎn)受力的矩形單元，類似兩個(gè)偏拉構(gòu)件，因而兩個(gè)單元內(nèi)不僅在荷載方向有正應(yīng)力，還有剪應(yīng)力存在，正應(yīng)力在單元的左右兩側(cè)最大，中間為零，剪應(yīng)力在單元的上下兩側(cè)最大，中間為零；采用擴(kuò)展有限元方法計(jì)算時(shí)，所得到的應(yīng)力結(jié)果與分離式裂縫模型的結(jié)果相同，從而可認(rèn)為本文說(shuō)明的施加荷載的方法是正確的．考慮到采用分離裂縫模型計(jì)算過(guò)程中重新剖分網(wǎng)格的工作量及難度，擴(kuò)展有限元在分析裂縫問(wèn)題時(shí)具有明顯的優(yōu)點(diǎn)．

對(duì)于更為一般的情況，討論集中荷載不作用于結(jié)點(diǎn)時(shí)的等效結(jié)點(diǎn)荷載．以圖2所示單元為例，在單元上邊緣±l/2處施加大小為F、方向向上的集中荷載．該情況下，無(wú)論單元視為一個(gè)整體，還是將裂縫兩側(cè)視為獨(dú)立的兩部分，單元均應(yīng)為單軸受拉受力狀態(tài)．

由式(7)可計(jì)算得到圖2單元的等效結(jié)點(diǎn)荷載為

(12)

圖2 開裂單元受非結(jié)點(diǎn)集中荷載

這與單向軸拉應(yīng)力狀態(tài)(即式(10))的等效結(jié)點(diǎn)荷載完全一致，因而可再次驗(yàn)證本文施加結(jié)點(diǎn)荷載方法的正確性．對(duì)不在節(jié)點(diǎn)上的集中荷載進(jìn)行積分即可得到面荷載或體荷載的等效結(jié)點(diǎn)荷載，受于篇幅的限制，本文不再進(jìn)行累述．

2 Shifted-XFEM模式下邊界條件的施加

2.1 Shifted-XFEM模式下的位移邊界條件

加強(qiáng)函數(shù)同樣選為Heaviside函數(shù)時(shí)，Shifted-XFEM的位移模式為

(13)

其中，H(xi)為節(jié)點(diǎn)處的Heaviside函數(shù)在i結(jié)點(diǎn)處的數(shù)值．

2.2 Shifted-XFEM模式下的應(yīng)力邊界條件

根據(jù)最小位能原理，富集自由度上所對(duì)應(yīng)的荷載為

(14)

仍采用圖2所示單元進(jìn)行說(shuō)明．根據(jù)式(14)計(jì)算不難發(fā)現(xiàn)，富集自由度上對(duì)應(yīng)的等效結(jié)點(diǎn)荷載均為0，因此，等效結(jié)點(diǎn)荷載向量可以表示為

(15)

同樣假設(shè)該單元的泊松比μ=0．根據(jù)2.1節(jié)方法施加位移邊界，由式(15)確定的結(jié)點(diǎn)荷載可求出節(jié)點(diǎn)位移向量，計(jì)算得到該單元內(nèi)的應(yīng)力分布σ與式(9)所示結(jié)果相同，即可得到正確的應(yīng)力解答．

對(duì)圖2所示的集中力不作用于單元結(jié)點(diǎn)的一般情況，由式(14)可以得到單元的等效結(jié)點(diǎn)荷載為

(16)

式(16)所給出的結(jié)點(diǎn)荷載與式(10)所示單軸受拉應(yīng)力狀態(tài)在Shifted-XFEM模式下的結(jié)點(diǎn)力相互平衡．

限于篇幅，本文不再討論由集中力等效結(jié)點(diǎn)荷載積分得到分布力等效結(jié)點(diǎn)荷載的過(guò)程．

3 算例驗(yàn)證

本文采用的算例與文獻(xiàn)[13]中所計(jì)算的受均布荷載的簡(jiǎn)支梁算例相同．采用自編的鋼筋混凝土擴(kuò)展有限元程序HohaiRCFE-P進(jìn)行計(jì)算，將計(jì)算結(jié)果與規(guī)范結(jié)果進(jìn)行比較．

如圖3所示簡(jiǎn)支梁，截面尺寸為250mm×600mm，長(zhǎng)度為7.2m，承受標(biāo)準(zhǔn)值為10.0kN/m的均布荷載．采用C25混凝土，配有兩根直徑為18mm的HRB335縱向受拉鋼筋，混凝土保護(hù)層厚度為35mm．采用增量法進(jìn)行計(jì)算，分20個(gè)荷載步施加，每一荷載步的增量為0.5kN/m，每一步迭代中將不平衡力視為結(jié)點(diǎn)荷載按本文方法進(jìn)行施加．共劃分876個(gè)混凝土單元，單元尺寸為100mm×50mm；鋼筋采用2結(jié)點(diǎn)桿單元(共73個(gè))，網(wǎng)格剖分如圖3所示．粘結(jié)滑移關(guān)系采用Houde-雙折線形式．

圖3 受均布荷載的簡(jiǎn)支梁

表1列出了計(jì)算得到的裂縫寬度，并與《混凝土結(jié)構(gòu)設(shè)計(jì)規(guī)范》(GB50010-2010)裂縫寬度公式計(jì)算結(jié)果進(jìn)行了對(duì)比．需要說(shuō)明的是，規(guī)范裂縫寬度公式考慮了長(zhǎng)期效應(yīng)(擴(kuò)大系數(shù)為1.5)，而本文擴(kuò)展有限元計(jì)算中沒(méi)有考慮荷載長(zhǎng)期作用的影響，因而表1中按規(guī)范計(jì)算的裂縫寬度除以了1.5的擴(kuò)大系數(shù)．

表1 最大裂縫寬度計(jì)算值比較 (單位：mm)

從表1可以看出，裂縫寬度較小時(shí)，與規(guī)范公式計(jì)算結(jié)果相比，本文方法計(jì)算得到的裂縫寬度誤差稍大，但當(dāng)裂縫寬度接近各類環(huán)境下的裂縫寬度限值時(shí)，本文結(jié)果與與規(guī)范裂縫寬度計(jì)算值吻合良好，從而驗(yàn)證了本文所述方法的正確性．

4 結(jié) 論

本文由有限元的控制方程和邊界條件的基本格式出發(fā)，說(shuō)明了擴(kuò)展有限元施加應(yīng)力和位移邊界的方法，以單個(gè)帶裂縫單元為例給出了Standard-XFEM和Shifted-XFEM模式下擴(kuò)展有限元施加邊界條件的具體實(shí)現(xiàn)過(guò)程．最終，通過(guò)一個(gè)受均布荷載的簡(jiǎn)支梁驗(yàn)證了本文介紹方法的正確性．由本文的工作，可以得到以下結(jié)論：

1)擴(kuò)展有限元施加邊界條件的原理與有限元相同，對(duì)位移邊界，結(jié)點(diǎn)實(shí)際位移應(yīng)與約束條件協(xié)調(diào)；對(duì)應(yīng)力邊界，可由最小位能原理確定等效結(jié)點(diǎn)荷載．

2)位移邊界條件的具體實(shí)現(xiàn)方法因位移模式的不同而不同．Standard-XFEM模式下，應(yīng)將位移邊界條件所滿足的方程帶入基本方程求解；Shifted-XFEM模式下位移邊界條件對(duì)富集自由度的廣義位移沒(méi)有約束作用．

3)不同位移模式下，由外部荷載計(jì)算等效結(jié)點(diǎn)荷載的公式相同，但得到的結(jié)點(diǎn)荷載向量不同．

4)擴(kuò)展有限元中，集中荷載是否作用在結(jié)點(diǎn)上產(chǎn)生不同的等效結(jié)點(diǎn)荷載，其應(yīng)力計(jì)算結(jié)果與采用分離式單元模型的有限元計(jì)算結(jié)果相同．除Standard-XFEM和Shifted-XFEM模式外，擴(kuò)展有限元還有其他常用位移模式．不同位移模式下，施加邊界條件的原理都是相同的，限于篇幅，本文不再重復(fù)．

[1] Theiner Y, Hofstetter G． Numerical Prediction of Crack Propagation and Crack Widths in Concrete Structures[J]．Steel Construction, 2009, 31(8)：1832-1840．

[2] 方修君, 金 峰, 王進(jìn)廷．用擴(kuò)展有限元方法模擬混凝土的復(fù)合型開裂過(guò)程[J]．工程力學(xué), 2007, 24(S1)：46-52．

[3] Mo?s N, Béchet E, Tourbier M． Imposing Essential Boundary Conditions in the Extended Finite Element Method[J]． International Journal for Numerical Methods in Engineering, 2006, 67(12)：1641-1669．

[4] 余天堂．摩擦接觸裂紋問(wèn)題的擴(kuò)展有限元法[J]．工程力學(xué), 2010(4)：84-89．

[5] 余天堂．?dāng)U展有限單元法 - 理論、應(yīng)用及程序[M]．南京：科學(xué)出版社, 2014．

[6] éric Béchet, N Mo?s, B Wohlmuth． A Stable Lagrange Multiplier Space for Stiff Interface Conditions within the Extended Finite Element Method[J]． International Journal for Numerical Methods in Engineering, 2008, 78(8)：931-954．

[7] Dolbow J, Mo?s N, Belytschko T． An Extended Finite Element Method for Modeling Crack Growth with Frictional Contact[J]． Computer Methods in Applied Mechanics & Engineering, 2001, 190(s51-52)：6825-6846．

[8] Haziza F． The Extended Finite Element Method and Its Implementation in 2D in the Aster Code[D]．Stockholm：Royal Institute of Technology, 2006．

[9] Liu F, Borja R I． A Contact Algorithm for Frictional Crack Propagation with the Extended Finite Element Method[J]． International Journal for Numerical Methods in Engineering, 2008, 76(10)：1489-1512．

[10] Kim T Y, Dolbow J, Laursen T． A Mortared Finite Element Method for Frictional Contact on Arbitrary Interfaces[J]． Computational Mechanics, 2007, 39(3)：223-235．

[11] 方修君, 金 峰, 王進(jìn)廷．基于擴(kuò)展有限元法的粘聚裂紋模型[J]．清華大學(xué)學(xué)報(bào)(自然科學(xué)版), 2007, 47(3)：344-347．

[12] 王勖成．有限單元法[M]．北京：清華大學(xué)出版社, 2002．

[13] 李 俊．基于鋼筋混凝土擴(kuò)展有限元的鋼筋應(yīng)力求解方法研究[D]．南京：河海大學(xué), 2015．

[責(zé)任編輯 周文凱]

Implementation of Boundary Conditions in Extended Finite Element Method

Leng Fei Zhang Ran

(College of Civil & Transportation Engineering, Hohai Univ., Nanjing 210098, China)

It is different to implement boundary conditions in extended finite element method(XFEM) and finite element method(FEM). But there is no detailed statement in any

so far. Starting with the basic format of control equations and boundary conditions in FEM, the paper introduces how to impose the force boundary conditions and displacement boundary conditions in XFEM. The detailed processes of the implementation of boundary conditions under Standard-XFEM and Shifted-XFEM are given by an example of a single cracked element. Finally, the method in the paper is validated by an example of a simply supported beam under uniform load. The results show that, for the stress boundary conditions, equivalent nodal loads should be calculated by the formula derived from the principle of minimum potential energy, under different displacement modes; although, the formula of equivalent nodal load has no different; the nodal load vectors are different. For the displacement boundary conditions, the different equations satisfying with constraint condition for different displacement mode should be introduced into the control equation. Although the above basic method is the same, how to implement it is different for different displacement mode.

extended finite element method; displacement boundary condition; force boundary condition; nonlinear calculation

10.13393/j.cnki.issn.1672-948X.2016.05.011

2016-03-23

國(guó)家自然科學(xué)基金(51179063)

張 然(1993-)，男， 碩士研究生，研究方向?yàn)榛炷两Y(jié)構(gòu)數(shù)值計(jì)算．E-mail：popeyecross@outlook.com

TU31

A

1672-948X(2016)05-0059-05