胡　凱，尹碩輝

(河海大學(xué) 力學(xué)與材料學(xué)院，江蘇 南京 211100)

?

平板斷裂問(wèn)題的多尺度擴(kuò)展有限元法

胡凱，尹碩輝

(河海大學(xué) 力學(xué)與材料學(xué)院，江蘇 南京 211100)

摘要：采用漸進(jìn)彎曲奇異函數(shù)和跨過(guò)裂紋面的不連續(xù)函數(shù)，加強(qiáng)常規(guī)的位移逼近空間，從而使計(jì)算網(wǎng)格獨(dú)立于裂紋，建立了貫穿裂紋Reissner-Mindlin板的多尺度擴(kuò)展有限元法。在裂紋附近區(qū)域采用小尺度網(wǎng)格，其他區(qū)域采用大尺度網(wǎng)格。在計(jì)算代價(jià)不大的情況下，考慮大型結(jié)構(gòu)中小裂紋的存在或者提高裂紋附近的精度。所有尺度單元都采用四結(jié)點(diǎn)四邊形板單元，四邊形任意結(jié)點(diǎn)板單元連接不同尺度單元。用互作用積分法計(jì)算裂尖應(yīng)力強(qiáng)度因子，算例分析檢驗(yàn)了本文方法的精度和有效性。

關(guān)鍵詞：Reissner-Mindlin板；貫穿裂紋；多尺度；擴(kuò)展有限元法；四邊形任意結(jié)點(diǎn)板單元；應(yīng)力強(qiáng)度因子

0引言

板結(jié)構(gòu)在船舶與海洋工程領(lǐng)域有著重要的作用，在服役期不可避免會(huì)出現(xiàn)裂紋，致使結(jié)構(gòu)在滿足強(qiáng)度要求下發(fā)生破壞并可能造成嚴(yán)重的災(zāi)難，因此，有必要研究板結(jié)構(gòu)的力學(xué)行為[1]。鑒于問(wèn)題的復(fù)雜性，數(shù)值模擬是分析開(kāi)裂板力學(xué)性能最有效的方法。擴(kuò)展有限元法(extended finite element method，XFEM)[2]是目前分析不連續(xù)問(wèn)題最有效的數(shù)值方法，已成功用于求解裂紋[3-6]、孔洞夾雜[7]、剪切帶演化[8]、接觸[9]和多相流[10]等問(wèn)題。

結(jié)構(gòu)分析中考慮到裂紋附近(特別是裂尖附近)存在高梯度場(chǎng)量[11]，因此，裂紋附近必須采用小尺度單元。若整個(gè)結(jié)構(gòu)采用小尺度單元，則計(jì)算量巨大。為此，通常在遠(yuǎn)離裂紋區(qū)采用大尺度單元，這樣就出現(xiàn)了不同尺度問(wèn)題，需要采用多尺度計(jì)算方法進(jìn)行求解。多尺度計(jì)算方法可分為基于均勻化的方法和基于層次分解技術(shù)的整體-局部法兩種。基于均勻化的多尺度法存在尺度分離和周期性的局限[12]；基于層次分解技術(shù)的整體-局部法的思想是在關(guān)注的子域采用小尺度網(wǎng)格，在其他子域使用大尺度網(wǎng)格，能克服均勻化法的局限性。文獻(xiàn)[13-16]將這兩種多尺度計(jì)算方法和 XFEM 相結(jié)合進(jìn)行斷裂問(wèn)題的多尺度分析。

基于層次分解技術(shù)的整體-局部法的多尺度計(jì)算，在不同尺度單元間存在一層任意結(jié)點(diǎn)單元。若能構(gòu)造任意結(jié)點(diǎn)單元，則可方便地處理不同尺度單元的連接問(wèn)題。文獻(xiàn)[17]基于Reissner-Mindlin板理論構(gòu)造了四邊形任意結(jié)點(diǎn)板單元。本文建立了Reissner-Mindlin板斷裂分析的多尺度XFEM，所有尺度的單元均采用四結(jié)點(diǎn)四邊形單元，四邊形任意結(jié)點(diǎn)板單元連接不同尺度單元，用互作用積分法計(jì)算應(yīng)力強(qiáng)度因子。并給出算例分析，檢驗(yàn)本文建立的多尺度擴(kuò)展有限元法求解板斷裂問(wèn)題的效果。

1Reissner-Mindlin板理論

圖1是厚度為t的貫穿裂紋板，其中，Γd為裂紋面。根據(jù)Reissner-Mindlin板理論，板內(nèi)任一點(diǎn)的位移可以表示為：

u=zβy(x,y)；v=-zβx(x,y)；w=w(x,y)，

(1)

圖1　貫穿裂紋的Reissner-Mindlin板

其中：βx和βy為板中面法線分別繞x軸和y軸的轉(zhuǎn)角；w為板的撓度。

板中面應(yīng)變可以表示為：

(2)

橫向剪應(yīng)變?yōu)椋?/p>

(3)

對(duì)于各向同性彈性板，應(yīng)力-應(yīng)變本構(gòu)關(guān)系可表示為：

(4)

其中：E和ν分別為彈性模量和泊松比；k為橫向剪應(yīng)力修正因數(shù)，一般取5/6。

彎矩M和橫向剪力Q定義為：

(5)

將式(4)代入式(5)得到：

M=Dbψ，Q=Dsγ，

(6)

板的勢(shì)能為：

(7)

圖2　多尺度網(wǎng)格(粗黑線為裂紋)

2多尺度擴(kuò)展有限元法

2.1　四邊形任意結(jié)點(diǎn)板單元

裂紋附近采用小尺度單元，其他區(qū)域采用大尺度單元。所有尺度均采用四邊形四結(jié)點(diǎn)單元，則在不同尺度單元間存在一層四邊形任意結(jié)點(diǎn)單元，如圖2所示。下面給出四邊形任意結(jié)點(diǎn)單元的形函數(shù)。

基于Np個(gè)多項(xiàng)式的點(diǎn)插值，用uh(ξ)逼近位移場(chǎng)u(ξ)，則uh(ξ)可表示為：

(8)

對(duì)于四結(jié)點(diǎn)四邊形板單元，多項(xiàng)式基函數(shù)為：

p(ξ)=[1,ξ,η,ξη]T，

(9)

其中：ξ、η為等參元局部坐標(biāo)。

由點(diǎn)插值可得：

uh(ξ)=aTp(ξ)=UTq-1p(ξ),

(10)

其中：

圖3　四邊形任意結(jié)點(diǎn)板單元

由式(8)～ 式(10)可得四結(jié)點(diǎn)四邊形板單元形函數(shù)為：

(11)

類似地，通過(guò)增加一些特殊的基函數(shù)以滿足點(diǎn)插值，就可以得到四邊形任意結(jié)點(diǎn)板單元的形函數(shù)。通過(guò)單元邊上需要滿足的插值類型來(lái)選取特殊的基函數(shù)。考慮一個(gè)線性的四邊形任意結(jié)點(diǎn)板單元，稱之為(4+k+m)結(jié)點(diǎn)單元，見(jiàn)圖3。單元上的這些結(jié)點(diǎn)可以分為3類：第1類結(jié)點(diǎn)為4個(gè)角點(diǎn)；第2類結(jié)點(diǎn)有k個(gè)，位于四邊形單元的上下兩條邊上；第3類結(jié)點(diǎn)有m個(gè)，位于四邊形單元的左右兩條邊上。

多項(xiàng)式基函數(shù)可表示為[17]：

(12)

q=p(ξi)。

(13)

由式(10)可得，(4+k+m)結(jié)點(diǎn)單元的形函數(shù)為：

[N1,…,N4+k+m]T=q-1p(ξ)。

(14)

四邊形任意結(jié)點(diǎn)板單元根據(jù)結(jié)點(diǎn)的不同會(huì)有很多種情況，而四結(jié)點(diǎn)四邊形板單元只是眾多情況中的一種，因此,單元的插值函數(shù)能夠構(gòu)造成同一種形式，這樣有利于編程。并且四邊形任意結(jié)點(diǎn)板單元上所有的結(jié)點(diǎn),一定要是周圍四結(jié)點(diǎn)四邊形板單元的結(jié)點(diǎn)。

2.2　積分方案

為了保證充分積分應(yīng)變場(chǎng)，采用如下的積分方案：(Ⅰ)四結(jié)點(diǎn)四邊形板單元。該部分單元和常規(guī)擴(kuò)展有限元單元的積分方案一樣。(Ⅱ)四邊形任意結(jié)點(diǎn)板單元。本文構(gòu)造的四邊形任意結(jié)點(diǎn)板單元的形函數(shù)在內(nèi)部上出現(xiàn)不連續(xù)，而數(shù)值積分中是不允許出現(xiàn)不連續(xù)的，所以將單元?jiǎng)澐殖梢恍┧慕Y(jié)點(diǎn)四邊形子單元，這些子單元的邊是形函數(shù)斜坡不連續(xù)處。采用2×2個(gè)高斯積分點(diǎn)使每個(gè)子單元內(nèi)形函數(shù)都是線性插值的。

3應(yīng)力強(qiáng)度因子計(jì)算

采用互作用積分法[18]計(jì)算Reissner-Mindlin板的彎矩和剪力強(qiáng)度因子。

對(duì)于Ⅰ型和Ⅱ型應(yīng)力強(qiáng)度因子，互作用積分為：

(15)

互作用積分與應(yīng)力強(qiáng)度因子KI、KⅡ和KⅢ的關(guān)系為：

(16)

選取輔助狀態(tài)為I型、II型或III型，根據(jù)式(16)就可分別求出 I型、II型或III型的應(yīng)力強(qiáng)度因子。

4數(shù)值算例

圖4　中心裂紋平板

算例1中心裂紋問(wèn)題

含一中心裂紋的正方形平板(見(jiàn)圖4)，板的兩端受到彎矩M0的作用。板材料的彈性模量E=210 GPa，泊松比ν=0.3。板的尺寸L=20 m，裂紋長(zhǎng)度為2a，板的厚度為t=2a，裂紋傾角為α。

無(wú)限大板含一中心斜裂紋的解析解為[19]：

(17)

其中：φ(1)和ψ(1)從積分方程中計(jì)算得來(lái)，當(dāng)板的厚度t=2a時(shí)，φ(1)=0.82，ψ(1)=0.68。

圖5為擴(kuò)展有限元計(jì)算網(wǎng)格。首先，對(duì)整個(gè)區(qū)域采用大尺度單元進(jìn)行離散，見(jiàn)圖5a，稱之為大尺度網(wǎng)格。然后，對(duì)裂紋附近單元進(jìn)一步細(xì)化，本文將一個(gè)大單元細(xì)化為3×3的小單元，如圖5b所示，稱之為多尺度網(wǎng)格。為了便于比較，將所有單元都細(xì)化為3×3的小單元，如圖5c所示，稱之為小尺度網(wǎng)格。

圖5　擴(kuò)展有限元計(jì)算網(wǎng)格

表1列出了不同類型網(wǎng)格下模型的自由度數(shù)目。由表1可知：多尺度問(wèn)題的自由度要遠(yuǎn)小于小尺度問(wèn)題的自由度。對(duì)于圖5a，由于裂紋附近單元尺度太大，以致無(wú)法用互作用積分法計(jì)算應(yīng)力強(qiáng)度因子。

裂紋附近存在高梯度應(yīng)力，多尺度擴(kuò)展有限元網(wǎng)格在裂紋附近采用了小尺度單元，因此，采用較少自由度網(wǎng)格的多尺度擴(kuò)展有限元法就能獲得較高精度的結(jié)果。多尺度擴(kuò)展有限元法能大大節(jié)省計(jì)算量。

圖6　應(yīng)力強(qiáng)度因子隨裂紋傾角的變化

裂紋傾角/(°)自由度數(shù)目多尺度小尺度0276012432152772124443032401246845286812444603240124687527721244490276012432

算例2邊裂紋問(wèn)題

圖7為含一邊裂紋的長(zhǎng)方形板，在板的兩端受到彎矩M0的作用，板的高寬比c/b=2，寬厚比b/h=2，裂紋長(zhǎng)度為a。板材料的彈性模量E=210 GPa，泊松比ν=0.3。

采用11×22個(gè)四結(jié)點(diǎn)四邊形單元離散板，形成的網(wǎng)格稱為大尺度網(wǎng)格。然后,對(duì)裂紋所在單元及外2層單元細(xì)化為3×3個(gè)單元，形成多尺度網(wǎng)格，如圖8所示。將大尺度網(wǎng)格的每個(gè)單元細(xì)化為3×3個(gè)單元，形成的網(wǎng)格稱為小尺度網(wǎng)格。

圖7 邊裂紋板圖8 多尺度擴(kuò)展有限元網(wǎng)格

圖9　歸一化的KⅠ隨裂紋長(zhǎng)度的變化

5結(jié)論

(1)建立了Reissner-Mindlin板斷裂問(wèn)題的多尺度擴(kuò)展有限元法，所有尺度單元都采用四結(jié)點(diǎn)四邊形板單元，用四邊形任意結(jié)點(diǎn)板單元來(lái)連接不同尺度單元。

(2)裂紋附近用小尺度單元進(jìn)行劃分，其他地方采用大尺度單元?jiǎng)澐?，這種多尺度擴(kuò)展有限元法既能降低計(jì)算成本，又能得到較好的精度。

以后的工作將根據(jù)誤差估計(jì)來(lái)科學(xué)地選取需要采用小尺度單元的區(qū)域，研究關(guān)于板斷裂問(wèn)題的自適應(yīng)多尺度擴(kuò)展有限元法。

參考文獻(xiàn)：

[1]張振國(guó),張娜,張秀麗,等.復(fù)合材料薄層板常溫固化收縮及翹曲變形[J].河南科技大學(xué)學(xué)報(bào)(自然科學(xué)版),2009,30(6):1-5.

[2]BELYTSCHKOT,BLACKT.Elasticcrackgrowthinfiniteelementswithminimalremeshing[J].Internationaljournalfornumericalmethodsinengineering,1999,45(5):601-620.

[3]SUKUMARN,CHOPPDL,MORANB.Extendedfiniteelementmethodandfastmarchingmethodforthree-dimensionalfatiguecrackpropagation[J].Engineeringfracturemechanics,2003,70(1):29-48.

[4]YUTT,RENQW.Modelingcrackinviscoelasticmediausingtheextendedfiniteelement[J].ScienceChinatechnologicalsciences,2011,54(6):1599-1606.

[5]ASADPOUREA,MOHAMMADIS,VAFAIA.Modelingcrackinorthotropicmediausingacoupledfiniteelementandpartitionofunitymethods[J].Finiteelementsinanalysisanddesign,2006,42(13):1165-1175.

[6]LIUP,YUT,BUITQ,etal.Transientthermalshockfractureanalysisoffunctionallygradedpiezoelectricmaterialsbytheextendedfiniteelementmethod[J].Internationaljournalofsolidsandstructures,2014,51(11):2167-2182.

[7]SUKUMARN,CHOPPDL,MOЁSN,etal.Modelingholesandinclusionsbylevelsetsintheextendedfinite-elementmethod[J].Computermethodsinappliedmechanicsandengineering,2001,190(46):6183-6200.

[8]SAMANIEGOE,BELYTSCHKOT.Continuum-discontinuummodellingofshearbands[J].Internationaljournalfornumericalmethodsinengineering,2005,62(13):1857-1872.

[9]SIAVELISM,GUITONMLE,MASSINP,etal.LargeslidingcontactalongbrancheddiscontinuitieswithX-FEM[J].Computationalmechanics,2013,52(1):201-219.

[10]CHESSA J,BELYTSCHKO T.An enriched finite element method and level sets for axisymmetric two-phase flow with surface tension[J].International journal for numerical methods in engineering,2003,58(13):2041-2064.

[11]HOLL M,LOEHNERT S,WRIGGERS P.An adaptive multiscale method for crack propagation and crack coalescence[J].International journal for numerical methods in engineering,2013,93(1):23-51.

[12]NEMAT-NASSER S,HORI M.Micromechanics:overall properties of heterogeneous materials[M].Amsterdam:North-Holland,1999.

[14]RANNOU J,GRAVOUIL A,BA?ETTO-DUBOURG M C.A local multigrid XFEM strategy for 3D crack propagation[J].International journal for numerical methods in engineering,2009,77(4):581-600.

[15]HOLL M,ROGGE T,LOEHNERT S,et al.3D Multiscale crack propagation using the XFEM applied to a gas turbine blade[J].Computational mechanics,2014,53(1):173-188.

[16]王振,余天堂.模擬三維裂紋問(wèn)題的多尺度擴(kuò)展有限元法[J].巖土力學(xué),2014,35(9):2702-2707.

[17]SOHN D,IM S.Variable-node plate and shell elements with assumed natural strain and smoothed integration methods for nonmatching meshes[J].Computational mechanics,2013,51(6):927-948.

[18]DOLBOW J,MOЁS N,BELYTSCHKO T.Modeling fracture in Mindlin-Reissner plates with the extended finite element method[J].International journal of solids and structures,2000,37(48):7161-7183.

[19]SIH G.Mechanics of fracture 3:plates and shells with cracks[M].The Netherlands:Noordhoff,1977.

[20]DIRGANTARA T,ALIABADI M H.Stress intensity factors for cracks in thin plates[J].Engineering fracture mechanics,2002,69(13):1465-1486.

文獻(xiàn)標(biāo)志碼：A

中圖分類號(hào)：TU311;P315.9

DOI:10.15926/j.cnki.issn1672-6871.2016.02.003

文章編號(hào)：1672-6871(2016)02-0011-06

收稿日期：2015-10-20

作者簡(jiǎn)介：胡凱(1989-),男,湖北黃岡人,碩士生，主要從事計(jì)算力學(xué)與工程仿真方面的研究.

基金項(xiàng)目：國(guó)家自然科學(xué)基金項(xiàng)目(51179063)