近場動力學框架下的鍵基對應模型

陳 壯，萬 冀，楚錫華*,2

(1.武漢大學 土木建筑工程學院，武漢 430072；2.大連理工大學 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，大連 116023)

1 引 言

近場動力學PD(Peridynamics)是一種基于非局部思想對傳統(tǒng)連續(xù)介質(zhì)理論進行重構(gòu)的物理模型，最早由Silling提出[1]。該方法將微分形式的平衡方程用積分形式進行改寫，使得平衡方程在裂紋等空間不連續(xù)處也可以成立。根據(jù)PD理論建立的數(shù)值模型可以有效地模擬動態(tài)裂紋的擴展現(xiàn)象。裂紋的產(chǎn)生、發(fā)展和分叉等現(xiàn)象在模擬時都不需要外部判斷條件。目前，PD由于其模擬動態(tài)破壞現(xiàn)象的能力，受到研究者們的關(guān)注[2-4]。

PD理論目前主要分為鍵基PD(bond-based PD)[1]和態(tài)基PD(state-based PD)[5]兩種，鍵基PD的局部相互作用機制比較簡單，可以模擬裂紋的動態(tài)擴展現(xiàn)象，但是鍵基PD存在泊松比限制問題(三維及平面應變問題，泊松比μ=1/4；平面應力問題，泊松比μ=1/3)，無法區(qū)分體積變形和偏變形，這給鍵基PD的應用帶來了一定的限制。態(tài)基PD是一種更加一般的PD模型，分為常規(guī)態(tài)(ordinary state-based PD)和非常規(guī)態(tài)(non-ordinary state-based PD)兩種，鍵基PD是態(tài)基模型的一種特例。態(tài)基PD雖然沒有泊松比取值的問題，但是基于態(tài)基PD發(fā)展的數(shù)值模型比鍵基PD要復雜，計算規(guī)模較大。此外，非常規(guī)態(tài)基PD中物質(zhì)點的變形也存在零能模式的問題[6]。許多學者采用修正鍵基PD的方式解決泊松比限制的問題。Gerstle等[7]提出了微極PD模型，在鍵基PD中引入新的微觀彈性常數(shù)，使得鍵基PD獨立的微觀材料參數(shù)在各向同性線彈性體中變成兩個；Zhu等[8]提出了考慮鍵轉(zhuǎn)角的鍵基PD模型，該模型可以很好地模擬泊松比取其他值時材料的力學行為；Zhou等[9]提出了考慮鍵間轉(zhuǎn)角的共軛鍵線彈性模型，并且修正了鍵斷裂的模型；Hu等[10]考慮鍵的法向和切向作用，從而解決了鍵的泊松比限制。

本文提出了一種鍵基對應PD模型，通過定義單根鍵的變形梯度來計算應變，進而計算應力和物質(zhì)點間的近場力。鍵變形梯度的定義在小變形與鍵長趨近無窮小時可以退化為經(jīng)典的變形梯度張量。在鍵對應PD模型中不需要對PD微觀彈性常數(shù)進行計算，直接采用物理方程求應力以及近場力。數(shù)值算例模擬了準靜態(tài)平板受拉和Kalthoff-Winkler板沖擊試驗。平板受拉算例對比鍵基對應PD和傳統(tǒng)鍵基PD的結(jié)果，驗證了鍵基對應PD在解決泊松比限制的有效性；對Kalthoff-Winkler板沖擊試驗的數(shù)值模擬，驗證了鍵基對應PD在模擬動態(tài)裂紋擴展上的有效性。

2 近場動力學簡介

考慮占據(jù)一定空間Ω的物體B，B中某一點x受到周圍一定區(qū)域Hx內(nèi)點x′的作用，如圖1所示，根據(jù)牛頓第二定律，得到PD的動力學方程為

(1)

鍵基PD中定義鍵的伸長率s為

s=(‖ξ+η‖-‖ξ‖)/‖ξ‖

(2)

式中ξ=x′-x為相對位置，η=u′-u表示相對位移。

在鍵基PD中，通過鍵傳遞的近場力f的大小和伸長率成正比，方向沿著鍵的方向:

(3)

c為PD微觀彈性參數(shù)，通過對比鍵基PD的應變能密度和傳統(tǒng)連續(xù)介質(zhì)力學的應變能密度得到[11]

c=18κ/πδ4

(4)

式中κ為體積模量，δ為近場范圍尺寸，關(guān)于c取值的相關(guān)討論可以參考文獻[12]。

鍵基PD模擬各向同性線彈性的材料時，獨立的材料參數(shù)只有一個，這和傳統(tǒng)理論中各向同性線彈性體應該有兩個獨立的材料參數(shù)有差別，而且鍵基PD無法區(qū)分體積變形和偏變形，這些問題導致了鍵基PD只能正確模擬泊松比取限制值的材料。

圖1 PD中物質(zhì)點的非局部相互作用

Fig.1 Interaction of material points

(5)

在鍵基PD中，s0的取值和極限能量釋放率Gc、鍵基PD微觀彈性參數(shù)c以及近場范圍尺寸δ有關(guān)

(6)

(7)

損傷值越大，表明連接點x的鍵斷裂的數(shù)目越多，點x附近的損傷越嚴重。

3 鍵變形定義

fx x′=RU

(8)

在二維問題中

(9)

(10)

在三維問題中，為簡化表述只畫出Δα的變化情況，其余兩個平面類似，如圖3所示。

圖2 二維情況下鍵變形的定義

Fig.2 Deformation of bond in 2D problems

圖3 三維情況下鍵角度變化在x-y平面的投影

Fig.3 Bond’s projection onx,yplane

(11,12)

R各分量具體是

R11=cos(Δβ)cos(Δθ)

(13a)

R12=cos(Δβ)sin(Δθ)

(13b)

R13=-sin(Δβ)

(13c)

R21=-cos(Δα)sin(Δθ)+

sin(Δα)sin(Δβ)cos(Δθ)

(13d)

R22=-cos(Δα)cos(Δθ)+

sin(Δα)sin(Δβ)sin(Δθ)

(13e)

R23=sin(Δα)cos(Δβ)

(13f)

R31=sin(Δα)sin(Δθ)+

cos(Δα)sin(Δβ)cos(Δθ)

(13g)

R32=-sin(Δα)cos(Δθ)+

cos(Δα)sin(Δβ)sin(Δθ)

(13h)

R33=cos(Δα)cos(Δβ)

(13i)

(14)

式中R為鍵的剛體轉(zhuǎn)動，U為鍵的伸長，r1為變形前的鍵長，r2是變形后的鍵長，假設垂直于鍵方向沒有發(fā)生變形。fx x′反映了鍵在局部坐標系下的變形情況，這種變形定義的合理性證明如下。

設物體受到外力作用后產(chǎn)生的變形張量場為F0(x,y,z)，對于該物體內(nèi)部任意一條鍵lM N,點的坐標M(x1,y1,z1)和N(x2,y2,z2)，在小變形假設下有R≈I，I為單位矩陣，鍵lM N的變形根據(jù)式(8)計算：

(15)

假設鍵的長度極小，點M和N無限接近，鍵lM N的局部變形張量F0(x1,y1,z1)≈F0(x2,y2,z2)，實際上鍵lM N在小變形以及鍵長無限小的情況下就等同于過點M或N的線元。根據(jù)變形張量的定義，在局部坐標下，鍵lM N的變形張量為

(16)

可以發(fā)現(xiàn)，在鍵長趨近無窮小以及小變形的情況下，鍵的變形張量定義和經(jīng)典連續(xù)介質(zhì)力學中變形張量的定義是一致的。

接下來通過鍵的變形張量計算應變和應力，在應力的基礎(chǔ)上計算近場力。首先進行坐標轉(zhuǎn)換，將局部坐標系轉(zhuǎn)換到全局坐標系中，T為坐標轉(zhuǎn)換矩陣，

(17)

計算位移梯度張量

(18)

在小應變下的應變定義為

(19)

線彈性條件下的應力

σx x′=Dεx x′

(20)

定義局部作用面積Sx x′來描述通過鍵傳遞應力的作用面積，Sx x′的取值和網(wǎng)格劃分的尺寸Δx有關(guān)，本文取Sx x′=(Δx)2為物質(zhì)點的截面積。根據(jù)鍵傳遞應力來計算近場力，考慮到鍵長和物質(zhì)點體積對近場力作用的影響，引入影響函數(shù)ω=ω(ξ,Vx′)，其中ξ為點x和x′的相對位置矢量，Vx′為點x′的體積，本文忽略距離的影響，只考慮物質(zhì)點體積的影響

ω(ξ,Vx′)=Vx′/VH

(21)

式中VH為近場域的體積。

通過鍵傳遞應力計算近場力

(22)

式中nx x′為鍵的方向，

nx x′=ξ/‖ξ‖

(23)

將式(22)代入式(1)

(24)

采用了這種近場力計算方法的PD模型稱為鍵基對應PD，在鍵基對應PD中，不需要計算傳統(tǒng)鍵基PD的微觀材料常數(shù)，直接利用物理方程求應力，進而計算近場力。

4 數(shù)值算例

4.1 準靜態(tài)平板拉伸

考慮受到單軸拉伸的正方形板，如圖4所示，板邊長為7 m，厚度為0.175 m，上下邊界加上固定位移ub=0.001 m，材料參數(shù)列入表1。采用顯式積分算法計算方程(24)，根據(jù)文獻[11]計算臨界的時間步長為Δtc=8.1×10-5s，考慮加入的人工阻尼項對計算的影響，取小于限制值的時間步長為Δt=1×10-7s。

對比鍵基PD和鍵基對應PD的計算結(jié)果，y向位移和x向位移分布如圖5所示。兩種PD模型計算的y方向位移分布均勻，但是x方向位移受到邊界效應的影響，存在一定的誤差?？紤]到鍵基泊松比限制的問題，將鍵基PD和鍵基對應PD在y=-0.4375 m的x方向位移與理論解(u*=-νεyx，εy是施加邊界荷載后y方向產(chǎn)生的應變，ν為泊松比)進行對比，如圖6所示，相對誤差計算公式為‖uP D-u*‖/‖u*‖×100%。可以看出，鍵基對應PD比鍵基PD更接近于理論解，驗證了鍵基對應PD解決泊松比限制問題的有效性。

圖4 平板拉伸

Fig.4 Geometry and boundary conditions of the plate

表1 平板拉伸計算參數(shù)

Tab.1 Mechanical parameters

參數(shù)數(shù)值密度/kg·m-38000彈性模量/GPa192泊松比0.1阻尼系數(shù)/kg·m-3·s-14×108

圖5 平板拉伸位移分布云圖(單位：m)

Fig.5 Displacements of bond-based PD and bond-based corresponding PD(unit：m)

4.2 Kalthoff-Winkler板沖擊試驗

為驗證鍵基對應PD模擬動態(tài)裂紋擴展的能力，進行Kalthoff-Winkler板沖擊試驗的數(shù)值模擬[13]，考慮到原試驗沖擊塊形狀對平板動態(tài)裂紋擴展影響較小，將問題簡化為二維，計算模型尺寸如 圖7 所示。根據(jù)文獻[14]，采用對稱化建模，提高計算效率。材料參數(shù)和計算參數(shù)列入表2，網(wǎng)格劃分尺寸Δx=0.001 m，近場范圍尺寸δ=3.015Δx。PD斷裂準則采用鍵斷裂準則，模擬板受沖擊作用發(fā)生脆性破壞的現(xiàn)象。

圖6 鍵基PD與鍵基對應PD沿著y=-0.4375 m上x方向位移分布與理論解的對比

Fig.6xdisplacements alongy=-0.4375 m and comparison with the theoretical solution

圖7 Kalthoff-Winkler板沖擊試驗建模

Fig.7 Geometric condition of the Kalthoff-Winkler plate

裂紋擴展的模擬結(jié)果如圖8所示，裂紋在初始裂紋尖端附近開始擴展，之后裂紋會先沿著與荷載呈41°夾角的方向進行擴展，隨后沿著與荷載呈67°夾角的方向進行擴展，同時有裂紋分叉的現(xiàn)象出現(xiàn)。這與已有的數(shù)值模擬結(jié)果[14,16](70°)以及實驗結(jié)果[15](70°)比較接近。隨著荷載向下傳遞，損傷開始在梁的下端產(chǎn)生，之后向著預制裂紋所在的方向發(fā)展，這和文獻[14]的現(xiàn)象一致。

圖8 Kalthoff-Winkler板沖擊試驗裂紋擴展隨時間變化

Fig.8 Crack path of Kalthoff-Winkler plate

5 結(jié) 論

本文結(jié)合鍵基PD和態(tài)基PD的優(yōu)點，提出了鍵基PD對應模型，解決了鍵基PD泊松比限制的問題。新的鍵張量變形張量在小變形以及鍵長趨近無限小的條件下可以退化為經(jīng)典的變形張量。在此基礎(chǔ)上利用物理方程計算應力，再通過應力來計算近場力，建立了經(jīng)典連續(xù)介質(zhì)力學與鍵基PD的聯(lián)系。數(shù)值算例驗證了當泊松比不取限制值時，鍵基對應PD相對于傳統(tǒng)鍵基PD在預測平板x方向位移的正確性。對Kalthoff-Winkler板沖擊試驗的模擬驗證了鍵基對應PD模擬動態(tài)裂紋擴展的有效性。

鍵基對應模型相比于傳統(tǒng)鍵基模型，計算成本稍高，對于三維問題中對鍵變形梯度的描述本文也只給出了公式而沒有作進一步的討論，同時對于鍵變形的定義與傳統(tǒng)變形梯度定義之間的關(guān)系還有待進一步的理論證明。