吳 凡， 段慶林

(大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116023)

1 引 言

近場(chǎng)動(dòng)力學(xué)PD(Peridynamics)模型是由Silling[1]提出的一種新興的非局部連續(xù)體理論，其以積分方程取代微分方程，允許位移場(chǎng)出現(xiàn)間斷，能夠?qū)B續(xù)和非連續(xù)問(wèn)題作統(tǒng)一處理。因此，PD對(duì)斷裂破壞問(wèn)題的分析模擬十分方便有效[2]。

PD分為鍵型[1]和態(tài)型[3]。其中，鍵型PD模型簡(jiǎn)單且計(jì)算速度快，不僅應(yīng)用于各種工程問(wèn)題，如水下爆破破冰[4]和邊坡失穩(wěn)[5]等，還廣泛應(yīng)用于各種新型材料的分析，如功能梯度材料[6]、復(fù)合材料[7]以及熱電材料[8]等。然而，PD的計(jì)算精度還有待進(jìn)一步提高[9]，主要是由于積分方程中的區(qū)域積分不能得到精確計(jì)算以及邊界附近鍵缺失導(dǎo)致的邊界效應(yīng)。

針對(duì)積分誤差問(wèn)題，已發(fā)展了多種修正方案。體積修正法[10]對(duì)具有不完整影響域的物質(zhì)點(diǎn)積分給出了以鍵長(zhǎng)為自變量的修正公式。HHB修正方案[11]根據(jù)積分區(qū)域與網(wǎng)格單元重疊的不同情況分別計(jì)算，進(jìn)一步提高了積分的準(zhǔn)確性。Seleson等[12]提出了PA-AC修正方案，能更精確地積分圓形影響域，但計(jì)算代價(jià)較高。Liu等[13]通過(guò)構(gòu)建PD與傳統(tǒng)連續(xù)介質(zhì)力學(xué)之間關(guān)于應(yīng)力與彈性模量的等價(jià)關(guān)系，為PD建立了新的微模量，并稱(chēng)之為數(shù)值微模量，發(fā)現(xiàn)采用數(shù)值微模量能夠顯著提高計(jì)算精度。

邊界效應(yīng)問(wèn)題是由于邊界粒子影響域的殘缺導(dǎo)致鍵缺失，虛假地降低了邊界附近的剛度[14]。針對(duì)該問(wèn)題，Bobaru等[15]根據(jù)邊界處物質(zhì)點(diǎn)的體積修正了鍵剛度。Madenci等[16]采用力密度法和能量法重構(gòu)邊界附近的應(yīng)變能密度，得到了相應(yīng)的修正系數(shù)。Macek等[17]基于虛位移原理引入了剛度修正系數(shù)。章青等[18]通過(guò)修正邊界處的微模量修正邊界效應(yīng)。虛擬節(jié)點(diǎn)法[19]是目前應(yīng)對(duì)該問(wèn)題的主要方法，通過(guò)在邊界外側(cè)布置虛擬節(jié)點(diǎn)，彌補(bǔ)了邊界點(diǎn)影響域的殘缺，有效改善了邊界效應(yīng)問(wèn)題。然而，邊界條件需要施加于虛擬節(jié)點(diǎn)，且需要根據(jù)不同的邊界條件對(duì)虛擬節(jié)點(diǎn)的位移進(jìn)行相應(yīng)的特殊操作，實(shí)現(xiàn)起來(lái)較為繁復(fù)。

本文以一維鍵型近場(chǎng)動(dòng)力學(xué)為研究對(duì)象，對(duì)上述兩個(gè)導(dǎo)致PD精度不高的問(wèn)題分別發(fā)展了改進(jìn)方法。針對(duì)區(qū)域積分問(wèn)題，本文遵循了數(shù)值微模量方法[13]的思想，但提出了根據(jù)應(yīng)變能密度而非應(yīng)力的等價(jià)關(guān)系對(duì)微模量實(shí)施修正。而且，本文方法無(wú)需進(jìn)行體積修正，更為簡(jiǎn)單直接。針對(duì)邊界效應(yīng)問(wèn)題，提出了建立位移和力邊界虛擬鍵的方法，無(wú)需布置虛擬節(jié)點(diǎn)，有效簡(jiǎn)化了方法的實(shí)現(xiàn)。

2 一維鍵型近場(chǎng)動(dòng)力學(xué)模型

近場(chǎng)動(dòng)力學(xué)的控制方程為

(1)

(2)

PD一般采用規(guī)則分布的節(jié)點(diǎn)(node)進(jìn)行空間離散[21](一維情況如圖1所示)，優(yōu)點(diǎn)是計(jì)算速度快，但計(jì)算精度還有待提高。離散后，式(2)可寫(xiě)為

(3)

式中Δx為格子間距(grid spacing)。顯然，距邊界不小于δ的節(jié)點(diǎn)(圖1的實(shí)心圓點(diǎn))影響域完整，而邊界附近的節(jié)點(diǎn)(圖1的實(shí)心方點(diǎn))影響域不完整，發(fā)生鍵缺失，導(dǎo)致邊界效應(yīng)。

圖1 一維近場(chǎng)動(dòng)力學(xué)空間離散

3 高精度近場(chǎng)動(dòng)力學(xué)方法

為提高近場(chǎng)動(dòng)力學(xué)的計(jì)算精度，本文針對(duì)一維鍵型PD的微模量和邊界效應(yīng)分別發(fā)展了改進(jìn)方法。

3.1 基于離散應(yīng)變能等效的微模量修正

PMB材料的一維鍵型PD在均勻拉伸下粒子的應(yīng)變能線密度為

(4)

式中w=0.5cs2ξ為鍵的微勢(shì)能。相同情況下，經(jīng)典連續(xù)介質(zhì)力學(xué)理論中一維桿的應(yīng)變能線密度為

(5)

式中E為彈性模量，ε為線應(yīng)變，A為桿的橫截面積。注意到此種情況下有s=ε，則通過(guò)式(4，5)的等效可得到PMB材料的PD微模量[20]

c=2EA/δ2

(6)

應(yīng)指出的是，上述等效的前提是式(4)的積分得到精確計(jì)算。然而，在離散情況下，由于數(shù)值積分的不準(zhǔn)確，直接采用式(6)將導(dǎo)致較大誤差。

針對(duì)該問(wèn)題，本文提出的解決辦法是在上述等效中采用離散而非連續(xù)情況下的PD應(yīng)變能密度?？紤]受均勻拉伸的一維PD模型，任意粒子的應(yīng)變能密度為該粒子所有鍵的微勢(shì)能之和，即

(7)

(8)

式中m=[δ/Δx],[]表示向下取整。將式(8)代入式(7)可得

(9)

再通過(guò)式(9)與式(5)的等價(jià)即可得到修正的微模量

(10)

后面的數(shù)值算例結(jié)果將表明，采用式(10)的修正微模量將能顯著提高計(jì)算精度。

3.2 邊界虛擬鍵

邊界附近粒子會(huì)由于影響域殘缺造成鍵缺失，進(jìn)而導(dǎo)致邊界效應(yīng)。本文處理該問(wèn)題的總體思路是通過(guò)建立邊界附近粒子與邊界之間的虛擬鍵來(lái)彌補(bǔ)這些粒子的鍵缺失，從而改善邊界效應(yīng)，分為固定邊界和力邊界。

3.2.1 固定邊界

固定邊界指的是將一維桿的端點(diǎn)位移固定為0。在PD方法中，該條件的實(shí)施通常是將離端點(diǎn)最近的粒子位移固定為0。由于邊界效應(yīng)的存在，這會(huì)導(dǎo)致端點(diǎn)附近區(qū)域位移和應(yīng)變的較大誤差。本文將通過(guò)增加邊界虛擬鍵來(lái)改善邊界效應(yīng)，然而該虛擬鍵的大小如何確定仍然未知，這關(guān)系到固定邊界條件是否能夠得到精確施加。注意到，結(jié)構(gòu)和載荷均對(duì)稱(chēng)時(shí)一維桿的中心點(diǎn)處位移恒為0(如 圖2 的桿B)，這為確定虛擬鍵的大小提供了思路。

考慮圖2桿A的PD求解，其左邊界x=xL為固定邊界。參照?qǐng)D2桿B，將桿A以左邊界為對(duì)稱(chēng)面進(jìn)行鏡像映射，得到如圖2桿C所示的虛擬對(duì)稱(chēng)區(qū)域，再基于此分析缺失的鍵(即由邊界截?cái)嗟逆I)，即可確定所需的虛擬鍵的大小。

(11)

其中

(12)

伸長(zhǎng)率可由其定義寫(xiě)為

(13)

(14)

(15)

將式(15)代入式(12)，并將結(jié)果代入式(11)即可得到粒子x的固定邊界虛擬鍵為

(xL

式(16)為連續(xù)情況下固定邊界虛擬鍵的計(jì)算公式。空間離散后，式中連續(xù)積分通過(guò)離散節(jié)點(diǎn)求和計(jì)算，對(duì)于邊界區(qū)域(圖2的陰影區(qū)域)距離固定邊界的第K個(gè)節(jié)點(diǎn)XK，其位移為UK，則其固定邊界虛擬鍵為

(17)

在PD計(jì)算中，為所有處于固定邊界區(qū)域的節(jié)點(diǎn)按式(17)建立邊界虛擬鍵，能有效改善邊界效應(yīng)。

圖2 固定邊界條件施加方案

3.2.2 力邊界

力邊界的邊界效應(yīng)也是由邊界區(qū)域粒子的鍵缺失導(dǎo)致。因此，類(lèi)似于上文固定邊界的處理方式，對(duì)力邊界區(qū)域的粒子也通過(guò)建立邊界虛擬鍵彌補(bǔ)鍵缺失，改善邊界效應(yīng)。與處理固定邊界不同，虛擬鍵的大小不由對(duì)稱(chēng)關(guān)系導(dǎo)出，而是基于遠(yuǎn)離邊界區(qū)域沒(méi)有邊界效應(yīng)這一事實(shí)，詳細(xì)闡述如下。

(xR

(0<γ<δ)(19)

(20)

由式(20)可求得伸長(zhǎng)率

(21)

將式(21)代入式(19)得到距離力邊界γ的粒子的力邊界虛擬鍵為

(0<γ<δ)(22)

同樣，式(22)也為連續(xù)情況下力邊界虛擬鍵的計(jì)算公式。空間離散后，對(duì)于距離力邊界第I個(gè)節(jié)點(diǎn)XI，其虛擬鍵應(yīng)為

(23)

式中伸長(zhǎng)率需由

(24)

計(jì)算得到

(25)

圖3 邊界力施加方案

將式(25)代入式(23)可得到節(jié)點(diǎn)XI的力邊界虛擬鍵

(26)

4 數(shù)值算例

4.1 兩端拉伸桿

如圖4所示的1 m長(zhǎng)桿，橫截面積A=1 m2，兩端受集中力F=5000 kN，材料密度為2400 kg/m3，楊氏模量為10 GPa。PD離散采用100個(gè)均勻分布的節(jié)點(diǎn)，并用動(dòng)態(tài)松弛法[22]計(jì)算至準(zhǔn)靜態(tài)。

圖4 兩端拉伸桿件

經(jīng)典PD方法和本文發(fā)展的高精度PD方法 hp -PD(high-precision PD)計(jì)算得到的應(yīng)變沿桿分布如圖5所示?？梢钥闯?，經(jīng)典PD的結(jié)果邊界效應(yīng)明顯，應(yīng)變?cè)跅U兩端區(qū)域發(fā)生劇烈的虛假震蕩，遠(yuǎn)離邊界區(qū)域的應(yīng)變也嚴(yán)重偏離解析解，誤差很大。本文通過(guò)實(shí)施微模量修正和邊界修正發(fā)展的 hp -PD 方法得到的應(yīng)變分布與解析解εe=5×10-4完全吻合，誤差在10-25量級(jí)，達(dá)機(jī)器精度。為厘清兩種修正方法各自的作用，圖5還比較了對(duì)經(jīng)典PD僅作微模量修正(標(biāo)識(shí)為PD+c)和僅作邊界修正(標(biāo)識(shí)為PD+b)兩種方法的數(shù)值結(jié)果?？梢钥闯?，兩種修正方法的作用相對(duì)獨(dú)立，一方面，微模量修正能顯著提升遠(yuǎn)離邊界處的計(jì)算精度，但對(duì)邊界附近的虛假震蕩無(wú)用；另一方面，邊界修正能消除邊界附近的虛假震蕩，但不能改善內(nèi)部區(qū)域的計(jì)算精度。因此，為全面提升PD方法的性能，需同時(shí)采用微模量和邊界修正。

采用該算例進(jìn)一步考察節(jié)點(diǎn)影響域大小對(duì)計(jì)算精度的影響。選取五種影響域，大小分別為δ=2Δx,3Δx,4Δx,5Δx,6Δx，并取桿中點(diǎn)處的應(yīng)變與解析解作比較，結(jié)果如圖6所示。本文的 hp -PD 方法在不同影響域下的計(jì)算結(jié)果均與解析解完全一致，成為機(jī)器精度下的解析解。標(biāo)準(zhǔn)PD的所有計(jì)算結(jié)果均嚴(yán)重偏離于解析解，誤差很大。一方面，減小影響域，PD計(jì)算精度明顯下降，這主要是由于影響域內(nèi)積分點(diǎn)數(shù)目的減少導(dǎo)致了更大的積分誤差；另一方面，即使對(duì)于最大影響域δ=6Δx的情況，PD的計(jì)算誤差仍然達(dá)到10%左右，與本文的 hp -PD 差距甚遠(yuǎn)。繼續(xù)增大影響域，PD的計(jì)算精度可得到些許改善，但同時(shí)也會(huì)導(dǎo)致計(jì)算量的急劇增大，因而不建議采用。

圖5 兩端拉伸桿應(yīng)變分布比較

4.2 一端固定一端拉伸桿

該算例主要用于考察本文方法施加固定邊界條件的準(zhǔn)確性。仍采用上一算例的長(zhǎng)桿，材料及數(shù)值離散均不變，只是將桿左端由力邊界改為固定邊界，如圖7所示。經(jīng)典PD和本文發(fā)展的 hp -PD 計(jì)算得到的應(yīng)變沿桿分布如圖8所示。與4.1節(jié)結(jié)果一致，經(jīng)典PD導(dǎo)致了很大的計(jì)算誤差并表現(xiàn)出明顯的邊界效應(yīng)。 hp -PD 的結(jié)果與解析解完全一致，不僅說(shuō)明其能精確施加固定邊界條件，而且再次展現(xiàn)了其高精度。

圖6 兩端拉伸算例不同影響域的應(yīng)變計(jì)算精度比較

圖7 一端固定一端拉伸桿

圖8 一端固定一端拉伸桿應(yīng)變分布比較

5 結(jié) 論

本文分別通過(guò)修正微模量和建立邊界虛擬鍵有效改善了一維鍵型PD方法的積分誤差和邊界效應(yīng)，建立了高精度PD方法 hp -PD ，顯著提升了PD的計(jì)算精度，準(zhǔn)靜態(tài)常應(yīng)變問(wèn)題能達(dá)到機(jī)器精度。應(yīng)強(qiáng)調(diào)的是，本文發(fā)展的 hp -PD 方法不破壞PD方法原有的優(yōu)點(diǎn)，不引入待定的計(jì)算參數(shù)，幾乎不增加額外的計(jì)算量，無(wú)需采用繁瑣的體積修正，也不引入需額外處理的虛擬粒子，實(shí)現(xiàn)起來(lái)簡(jiǎn)單方便，值得進(jìn)一步向二維和三維問(wèn)題拓展。