近場(chǎng)動(dòng)力學(xué)與有限元混合模型研究

錢 偉，范存新，沈 峰，夏益兵

(蘇州科技大學(xué) 土木工程學(xué)院，江蘇 蘇州 215000)

以有限單元法為代表的基于偏微分方程的相關(guān)數(shù)值方法已經(jīng)得到了廣泛的運(yùn)用，但在模擬裂紋擴(kuò)展等不連續(xù)問題時(shí)，基于連續(xù)性假設(shè)的有限單元法無(wú)法自發(fā)模擬裂紋擴(kuò)展，國(guó)內(nèi)外學(xué)者通過(guò)設(shè)置界面單元或采用網(wǎng)格重劃分技術(shù)來(lái)處理，但仍存在網(wǎng)格依賴性的問題[1]。

2000年，美國(guó)Sandia國(guó)家實(shí)驗(yàn)室的Silling博士提出了基于非局部作用思想的近場(chǎng)動(dòng)力學(xué)方法[2-3]，在國(guó)際上引起了廣泛關(guān)注，從根本上解決了連續(xù)介質(zhì)力學(xué)在模擬裂紋路徑等不連續(xù)問題時(shí)偏導(dǎo)數(shù)不存在的問題，逐漸成為計(jì)算力學(xué)與工程仿真及相關(guān)領(lǐng)域研究熱點(diǎn)[4]，在斷裂破壞問題中得到了廣泛應(yīng)用[5-8]。然而，近場(chǎng)動(dòng)力學(xué)計(jì)算效率相較有限元而言過(guò)低。為此，國(guó)內(nèi)外學(xué)者嘗試將近場(chǎng)動(dòng)力學(xué)與有限元結(jié)合起來(lái)，充分發(fā)揮兩者優(yōu)勢(shì)。Macek和Silling[9]將研究對(duì)象劃分為PD子域與FE子域以及重疊域，重疊域采用有限元實(shí)體模型，將重疊域中的PD模型以及PD子域視作桿單元，通過(guò)鑲嵌單元技術(shù)，實(shí)現(xiàn)了對(duì)PD與FEM的混合建模；Liu和Hong[10]在近場(chǎng)動(dòng)力學(xué)子域與有限元子域設(shè)置界面單元，在界面單元內(nèi)加入一定的物質(zhì)點(diǎn)，計(jì)算耦合力并將耦合力通過(guò)兩種方式分配到界面單元的結(jié)點(diǎn)上，以此來(lái)實(shí)現(xiàn)兩種方法的混合建模；Seleson等[11-12]通過(guò)在[0，1]區(qū)間內(nèi)變化的混合函數(shù)將局部作用區(qū)域的應(yīng)變能密度與PD應(yīng)變能密度混合實(shí)現(xiàn)了兩區(qū)域的平滑過(guò)渡；Ren等[13]利用不連續(xù)的Galerkin法建立了經(jīng)典近場(chǎng)動(dòng)力學(xué)控制方程，并給出了相關(guān)算例。上述混合模型均在顯式體系下求解，處理靜力問題時(shí)需引入人工阻尼，阻尼的大小影響迭代收斂速度，進(jìn)而導(dǎo)致計(jì)算效率不足。此外，上述文章均局限于通過(guò)低階單元與近場(chǎng)動(dòng)力學(xué)方法實(shí)現(xiàn)混合建模，缺乏對(duì)高階單元與近場(chǎng)動(dòng)力學(xué)混合建模研究。

在此基礎(chǔ)上，本文所建立的模型將研究對(duì)象劃分為PD子域與FE子域，在裂紋出現(xiàn)的區(qū)域，采用近場(chǎng)動(dòng)力學(xué)模型，其他區(qū)域采用等參元模型。通過(guò)桿單元連接PD物質(zhì)點(diǎn)與等參元結(jié)點(diǎn)來(lái)實(shí)現(xiàn)混合建模。該模型無(wú)需引入人工阻尼，提高了計(jì)算效率。同時(shí)，該混合模型在FE子域采用八結(jié)點(diǎn)等參元，提高了計(jì)算精度，最后，通過(guò)對(duì)懸臂梁彈性變形及含裂紋正方形板破壞過(guò)程的模擬，證明了該混合模型的有效性。

1 近場(chǎng)動(dòng)力學(xué)簡(jiǎn)介

1.1 近場(chǎng)動(dòng)力學(xué)方法

如圖1所示，設(shè)在任意時(shí)刻t,任意空間R內(nèi)任一物質(zhì)點(diǎn)x與其鄰域一定范圍δ內(nèi)的其他物質(zhì)點(diǎn)(x′∈R:‖x′-x‖≤δ)存在單位相互作用力f，根據(jù)牛頓第二定律，可得[14]：

圖1 物質(zhì)點(diǎn)間的相互作用Fig. 1 Interaction of material points

(1)

令式(1)中x′-x=ξ，u′-u=η，ξ為物質(zhì)點(diǎn)間相對(duì)位置，η為物質(zhì)點(diǎn)間相對(duì)位移。

式中：Hx為物質(zhì)點(diǎn)x的鄰域范圍，δ為鄰域范圍尺寸，ρ為物質(zhì)點(diǎn)材料密度，u為物質(zhì)點(diǎn)的位移，b為外力密度。

基于保守場(chǎng)的定義和性質(zhì)，一定存在一個(gè)標(biāo)量函數(shù)w(η,ξ)(物質(zhì)點(diǎn)對(duì)勢(shì)能密度)，使得：

(2)

物質(zhì)點(diǎn)間的作用類似一個(gè)中心彈簧，則

(3)

(4)

μ為一標(biāo)量函數(shù)，用來(lái)表征鍵的斷裂：

(5)

式中:s0為臨界伸長(zhǎng)率。當(dāng)伸長(zhǎng)率小于臨界伸長(zhǎng)率s0時(shí)，μ(t,ξ)=1，鍵未發(fā)生斷裂，否則，μ(t,ξ)=0，鍵斷裂。這里考慮材料的拉壓異性，可以通過(guò)材料的抗拉強(qiáng)度、抗壓強(qiáng)度和彈性模量E去表征：

(6)

在近場(chǎng)動(dòng)力學(xué)理論中，統(tǒng)一定義局部損傷

(7)

式中：0≤φ≤1，φ=0表示材料未損傷，φ=1表示該點(diǎn)完全損傷。

1.2 改進(jìn)的PMB模型

式(3)中的微觀模量函數(shù)c(ξ)可表示為：

c(ξ)=c(0,δ)g(ξ,δ)，

(8)

式中：c(0,δ)為集中函數(shù)，g(ξ,δ)為核函數(shù)，表示遠(yuǎn)程力大小隨物質(zhì)點(diǎn)間距變化的規(guī)律。在改進(jìn)的PMB模型中，考慮遠(yuǎn)程力對(duì)微觀模量的影響，取核函數(shù)[16-17]為：

(9)

根據(jù)近場(chǎng)動(dòng)力學(xué)應(yīng)變能密度與連續(xù)介質(zhì)力學(xué)應(yīng)變能密度相等的原則，可以得到平面應(yīng)力狀態(tài)下的集中函數(shù)為：

(10)

2 有限元簡(jiǎn)介

2.1 有限元方程

有限單元法靜力問題求解方程為

KU=R，

(11)

式中：K為整體剛度矩陣，U為整體位移列陣，R為整體等效結(jié)點(diǎn)荷載列陣。K和R可由式(12)得出

(12)

式中：Ce為選擇矩陣，k為單元?jiǎng)偠染仃嚕琑e為等效荷載列陣，B為應(yīng)變轉(zhuǎn)換矩陣，D為彈性矩陣。

2.2 八結(jié)點(diǎn)等參元

相較于矩形單元，等參元精度高，且適用于復(fù)雜的曲線邊界與曲面邊界，因此得到了廣泛的應(yīng)用。這里在有限元子域采用八結(jié)點(diǎn)等參元。

圖2為四邊形單元在整體坐標(biāo)系xy下的單元結(jié)點(diǎn)分布，在每個(gè)四邊形單元上建立局部坐標(biāo)系ξη，如圖3所示，通過(guò)坐標(biāo)變換，將每個(gè)四邊形單元映射到標(biāo)準(zhǔn)正方形單元，建立了兩種單元的對(duì)應(yīng)關(guān)系。

圖2 四邊形單元Fig. 2 Quadrilateral element

圖3 正方形單元Fig. 3 Square element

設(shè)單元中任意一點(diǎn)的位移是u，v，單元結(jié)點(diǎn)位移為ui，vi(i=1，8)，其位移模式為

(13)

坐標(biāo)變換為：

(14),

式中，Ni=(i=1,8)為八結(jié)點(diǎn)等參元形函數(shù)，其表達(dá)式為：

(15)

式中，ξ,η是定義在標(biāo)準(zhǔn)單元上的局部坐標(biāo)，ξi，ηi(i=1,2,3,4)分別代表標(biāo)準(zhǔn)單元4個(gè)結(jié)點(diǎn)的局部坐標(biāo)值。

3 近場(chǎng)動(dòng)力學(xué)與有限元耦合方案

如圖4所示，將幾何模型劃分為PD子域與FE子域，在PD子域采用近場(chǎng)動(dòng)力學(xué)建模，在FE子域采用八結(jié)點(diǎn)等參元建模。在兩類子域的交界面上，如圖5所示，交界面上的每個(gè)有限元結(jié)點(diǎn)通過(guò)桿單元與其近場(chǎng)范圍內(nèi)的物質(zhì)點(diǎn)相連接[18-19]，其中，PD子域物質(zhì)點(diǎn)對(duì)間相互作用可視為桿單元。為求解近場(chǎng)動(dòng)力學(xué)方程，將研究對(duì)象離散為一系列帶有物性信息的物質(zhì)點(diǎn)，近場(chǎng)動(dòng)力學(xué)積分方程轉(zhuǎn)化為對(duì)有限個(gè)物質(zhì)點(diǎn)的求和，即：

圖4 有限元(FE)與近場(chǎng)動(dòng)力學(xué)(PD)子域Fig. 4 FE and PD region

圖5 FE結(jié)點(diǎn)與PD物質(zhì)點(diǎn)連接示意圖Fig. 5 Interactions between FE nodes and PD nodes through trusses

(16)

對(duì)于靜力問題，令加速度為0，可得近場(chǎng)動(dòng)力學(xué)的平衡方程：

(17)

根據(jù)公式(1)，可將PD方程改寫成矩陣形式

(18)

(19)

桿單元的剛度貢獻(xiàn)矩陣為：

(20)

通過(guò)求解式(12)、式(19)、式(20)，分別得到了等參元、物質(zhì)點(diǎn)對(duì)間、桿單元的剛度貢獻(xiàn)矩陣，其中，物質(zhì)點(diǎn)對(duì)間相互作用視為桿單元，通過(guò)對(duì)剛度貢獻(xiàn)矩陣的集成，形成整體剛度矩陣。最后根據(jù)有限元靜力方程求解位移，實(shí)現(xiàn)了近場(chǎng)動(dòng)力學(xué)與有限元的混合建模。

4 數(shù)值算例

4.1 懸臂梁的彈性變形

懸臂梁幾何尺寸及子域劃分如圖6所示[19]，跨長(zhǎng)1 000 mm，截面高為200 mm，彈性模量E為100 GPa，泊松比ν=1/3，物質(zhì)點(diǎn)間距取為2.5 mm，近場(chǎng)范圍尺寸取δ=4Δx，有限元網(wǎng)格為10 mm×10 mm的八結(jié)點(diǎn)等參單元，右端受大小為1 050 kN/m的均布荷載，探討不同數(shù)值方法對(duì)精度的影響。

圖6 懸臂梁幾何模型Fig. 6 Geometric model of cantilever beam

表1給出了采用不同數(shù)值方法得到的最大水平位移計(jì)算結(jié)果。由表1可得，采用近場(chǎng)動(dòng)力學(xué)方法計(jì)算的最大水平位移與解析解的相對(duì)誤差為2.19%，采用四結(jié)點(diǎn)混合模型時(shí)相對(duì)誤差為0.67%，采用八結(jié)點(diǎn)混合模型的誤差為0.28%，混合模型的精度比近場(chǎng)動(dòng)力學(xué)的計(jì)算精度高，且本文中建立的八結(jié)點(diǎn)混合模型高于四結(jié)點(diǎn)混合模型的計(jì)算精度。證明了本文中提出的混合模型的精確性。其中，當(dāng)不設(shè)置FE子域時(shí)，可根據(jù)有限元靜力求解方程，得到近場(chǎng)動(dòng)力學(xué)計(jì)算結(jié)果。

表1 最大水平位移計(jì)算結(jié)果

4.2 含I型裂紋板受拉破壞分析

考慮如圖7所示含I型裂紋正方形板，模型尺寸及材料參數(shù)如下，邊長(zhǎng)為50 mm，在板中間預(yù)制一條長(zhǎng)為10 mm的裂紋，彈性模量為30 GPa，ft為2.01 MPa，fc為20.1 MPa，泊松比為1/3。將正方形板劃分為兩個(gè)有限元子域與一個(gè)近場(chǎng)動(dòng)力學(xué)子域，物質(zhì)點(diǎn)間距取為0.5 mm，近場(chǎng)范圍尺寸取δ=4Δx，有限元網(wǎng)格為2 mm×2 mm八結(jié)點(diǎn)等參元，采用位移加載，每一步位移增量為1.0×10-8m。

圖7 正方形板幾何模型Fig. 7 Geometric model of square plate

混合模型計(jì)算得到的正方形板裂紋擴(kuò)展過(guò)程如圖8所示，當(dāng)加載到57步時(shí)(此時(shí)位移荷載為5.7×10-7m),預(yù)制裂紋的裂尖出現(xiàn)損傷；隨著荷載的進(jìn)一步增大，裂紋發(fā)生擴(kuò)展，如圖8(b)(c)所示；當(dāng)加載至92步時(shí)，(此時(shí)位移荷載為9.2×10-7m)，裂紋貫穿整個(gè)正方形板，構(gòu)件發(fā)生破壞。

圖8 裂紋擴(kuò)展示意圖Fig. 8 Crack propagation process

4.3 多裂紋板受拉破壞分析

考慮如圖9所示的含兩條斜裂紋的正方形板，模型尺寸及材料參數(shù)如下，板的邊長(zhǎng)為50 mm，初始裂紋長(zhǎng)度為10 mm，b為10 mm，彈性模量為30 GPa，ft為2.01 MPa，fc為20.1 MPa，泊松比為1/3。將正方形板劃分為兩個(gè)有限元子域與一個(gè)近場(chǎng)動(dòng)力學(xué)子域，物質(zhì)點(diǎn)間距取為0.5 mm,近場(chǎng)范圍尺寸取δ=4Δx，有限元網(wǎng)格為2 mm×2 mm的八結(jié)點(diǎn)等參單元，采用位移加載，每一步位移增量為1.5×10-8m。

圖9 正方形板幾何模型Fig. 9 Geometric model of square plate

混合模型計(jì)算得到的正方形板裂紋擴(kuò)展過(guò)程如圖10所示，加載到39步時(shí)(此時(shí)位移荷載為5.85×10-7m)，裂紋裂尖出現(xiàn)損傷；隨著荷載的進(jìn)一步增大，裂紋發(fā)生擴(kuò)展并交匯，如圖10(b)(c)所示；當(dāng)加載至63步時(shí)，(此時(shí)位移荷載為9.45×10-7m)，裂紋貫穿整個(gè)正方形板，構(gòu)件發(fā)生破壞。圖10(e)為多維虛內(nèi)鍵計(jì)算結(jié)果[20]。對(duì)比可得，混合模型與虛內(nèi)鍵計(jì)算結(jié)果基本吻合。

圖10 裂紋擴(kuò)展示意圖Fig. 10 Crack propagation process

5 結(jié) 論

PD理論通過(guò)求解積分方程模擬材料斷裂破壞行為，在分析裂紋擴(kuò)展等不連續(xù)問題時(shí)具有顯著優(yōu)勢(shì)，然而，PD理論計(jì)算效率相較FEM而言過(guò)低。為兼顧兩者優(yōu)勢(shì)，采用近場(chǎng)動(dòng)力學(xué)與有限元混合建模的方法，建立了新的混合模型，該模型將研究對(duì)象劃分為PD子域與FE子域，在PD子域采用近場(chǎng)動(dòng)力學(xué)建模型，F(xiàn)E子域采用等參元模型，通過(guò)桿單元連接PD物質(zhì)點(diǎn)與等參單元結(jié)點(diǎn)來(lái)實(shí)現(xiàn)混合建模。此混合模型在求解處理靜力問題時(shí)無(wú)需引入人工阻尼，采取類似有限元的方法，對(duì)剛度集成，形成整體剛度矩陣，然后根據(jù)有限元支配方程求解靜力問題，提高了計(jì)算效率。同時(shí)在FE子域采用八結(jié)點(diǎn)等參元建模，提高了混合模型的計(jì)算精度。最后，采用所建立的混合模型計(jì)算分析了懸臂梁的彈性變形和含裂紋正方形板的破壞過(guò)程，取得了較好的結(jié)果，為斷裂破壞問題的解決提供了一種新思路。