吳 凡， 段慶林

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

1 引 言

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

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

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

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

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

2 一維鍵型近場動力學(xué)模型

近場動力學(xué)的控制方程為

(1)

(2)

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

(3)

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

圖1 一維近場動力學(xué)空間離散

3 高精度近場動力學(xué)方法

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

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

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

(4)

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

(5)

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

c=2EA/δ2

(6)

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

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

(7)

(8)

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

(9)

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

(10)

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

3.2 邊界虛擬鍵

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

3.2.1 固定邊界

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

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

(11)

其中

(12)

伸長率可由其定義寫為

(13)

(14)

(15)

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

(xL

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

(17)

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

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

3.2.2 力邊界

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

(xR

(0<γ<δ)(19)

(20)

由式(20)可求得伸長率

(21)

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

(0<γ<δ)(22)

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

(23)

式中伸長率需由

(24)

計算得到

(25)

圖3 邊界力施加方案

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

(26)

4 數(shù)值算例

4.1 兩端拉伸桿

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

圖4 兩端拉伸桿件

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

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

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

4.2 一端固定一端拉伸桿

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

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

圖7 一端固定一端拉伸桿

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

5 結(jié) 論

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