陳 壯,萬 冀,楚錫華*,2
(1.武漢大學(xué) 土木建筑工程學(xué)院,武漢 430072;2.大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室,大連 116023)
近場動力學(xué)PD(Peridynamics)是一種基于非局部思想對傳統(tǒng)連續(xù)介質(zhì)理論進(jìn)行重構(gòu)的物理模型,最早由Silling提出[1]。該方法將微分形式的平衡方程用積分形式進(jìn)行改寫,使得平衡方程在裂紋等空間不連續(xù)處也可以成立。根據(jù)PD理論建立的數(shù)值模型可以有效地模擬動態(tài)裂紋的擴(kuò)展現(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的局部相互作用機(jī)制比較簡單,可以模擬裂紋的動態(tài)擴(kuò)展現(xiàn)象,但是鍵基PD存在泊松比限制問題(三維及平面應(yīng)變問題,泊松比μ=1/4;平面應(yīng)力問題,泊松比μ=1/3),無法區(qū)分體積變形和偏變形,這給鍵基PD的應(yīng)用帶來了一定的限制。態(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要復(fù)雜,計算規(guī)模較大。此外,非常規(guī)態(tài)基PD中物質(zhì)點的變形也存在零能模式的問題[6]。許多學(xué)者采用修正鍵基PD的方式解決泊松比限制的問題。Gerstle等[7]提出了微極PD模型,在鍵基PD中引入新的微觀彈性常數(shù),使得鍵基PD獨(dú)立的微觀材料參數(shù)在各向同性線彈性體中變成兩個;Zhu等[8]提出了考慮鍵轉(zhuǎn)角的鍵基PD模型,該模型可以很好地模擬泊松比取其他值時材料的力學(xué)行為;Zhou等[9]提出了考慮鍵間轉(zhuǎn)角的共軛鍵線彈性模型,并且修正了鍵斷裂的模型;Hu等[10]考慮鍵的法向和切向作用,從而解決了鍵的泊松比限制。
本文提出了一種鍵基對應(yīng)PD模型,通過定義單根鍵的變形梯度來計算應(yīng)變,進(jìn)而計算應(yīng)力和物質(zhì)點間的近場力。鍵變形梯度的定義在小變形與鍵長趨近無窮小時可以退化為經(jīng)典的變形梯度張量。在鍵對應(yīng)PD模型中不需要對PD微觀彈性常數(shù)進(jìn)行計算,直接采用物理方程求應(yīng)力以及近場力。數(shù)值算例模擬了準(zhǔn)靜態(tài)平板受拉和Kalthoff-Winkler板沖擊試驗。平板受拉算例對比鍵基對應(yīng)PD和傳統(tǒng)鍵基PD的結(jié)果,驗證了鍵基對應(yīng)PD在解決泊松比限制的有效性;對Kalthoff-Winkler板沖擊試驗的數(shù)值模擬,驗證了鍵基對應(yīng)PD在模擬動態(tài)裂紋擴(kuò)展上的有效性。
考慮占據(jù)一定空間Ω的物體B,B中某一點x受到周圍一定區(qū)域Hx內(nèi)點x′的作用,如圖1所示,根據(jù)牛頓第二定律,得到PD的動力學(xué)方程為
(1)
鍵基PD中定義鍵的伸長率s為
s=(‖ξ+η‖-‖ξ‖)/‖ξ‖
(2)
式中ξ=x′-x為相對位置,η=u′-u表示相對位移。
在鍵基PD中,通過鍵傳遞的近場力f的大小和伸長率成正比,方向沿著鍵的方向:
(3)
c為PD微觀彈性參數(shù),通過對比鍵基PD的應(yīng)變能密度和傳統(tǒng)連續(xù)介質(zhì)力學(xué)的應(yīng)變能密度得到[11]
c=18κ/πδ4
(4)
式中κ為體積模量,δ為近場范圍尺寸,關(guān)于c取值的相關(guān)討論可以參考文獻(xiàn)[12]。
鍵基PD模擬各向同性線彈性的材料時,獨(dú)立的材料參數(shù)只有一個,這和傳統(tǒng)理論中各向同性線彈性體應(yīng)該有兩個獨(dú)立的材料參數(shù)有差別,而且鍵基PD無法區(qū)分體積變形和偏變形,這些問題導(dǎo)致了鍵基PD只能正確模擬泊松比取限制值的材料。
圖1 PD中物質(zhì)點的非局部相互作用
Fig.1 Interaction of material points
(5)
在鍵基PD中,s0的取值和極限能量釋放率Gc、鍵基PD微觀彈性參數(shù)c以及近場范圍尺寸δ有關(guān)
(6)
(7)
損傷值越大,表明連接點x的鍵斷裂的數(shù)目越多,點x附近的損傷越嚴(yán)重。
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是變形后的鍵長,假設(shè)垂直于鍵方向沒有發(fā)生變形。fx x′反映了鍵在局部坐標(biāo)系下的變形情況,這種變形定義的合理性證明如下。
設(shè)物體受到外力作用后產(chǎn)生的變形張量場為F0(x,y,z),對于該物體內(nèi)部任意一條鍵lM N,點的坐標(biāo)M(x1,y1,z1)和N(x2,y2,z2),在小變形假設(shè)下有R≈I,I為單位矩陣,鍵lM N的變形根據(jù)式(8)計算:
(15)
假設(shè)鍵的長度極小,點M和N無限接近,鍵lM N的局部變形張量F0(x1,y1,z1)≈F0(x2,y2,z2),實際上鍵lM N在小變形以及鍵長無限小的情況下就等同于過點M或N的線元。根據(jù)變形張量的定義,在局部坐標(biāo)下,鍵lM N的變形張量為
(16)
可以發(fā)現(xiàn),在鍵長趨近無窮小以及小變形的情況下,鍵的變形張量定義和經(jīng)典連續(xù)介質(zhì)力學(xué)中變形張量的定義是一致的。
接下來通過鍵的變形張量計算應(yīng)變和應(yīng)力,在應(yīng)力的基礎(chǔ)上計算近場力。首先進(jìn)行坐標(biāo)轉(zhuǎn)換,將局部坐標(biāo)系轉(zhuǎn)換到全局坐標(biāo)系中,T為坐標(biāo)轉(zhuǎn)換矩陣,
(17)
計算位移梯度張量
(18)
在小應(yīng)變下的應(yīng)變定義為
(19)
線彈性條件下的應(yīng)力
σx x′=Dεx x′
(20)
定義局部作用面積Sx x′來描述通過鍵傳遞應(yīng)力的作用面積,Sx x′的取值和網(wǎng)格劃分的尺寸Δx有關(guān),本文取Sx x′=(Δx)2為物質(zhì)點的截面積。根據(jù)鍵傳遞應(yīng)力來計算近場力,考慮到鍵長和物質(zhì)點體積對近場力作用的影響,引入影響函數(shù)ω=ω(ξ,Vx′),其中ξ為點x和x′的相對位置矢量,Vx′為點x′的體積,本文忽略距離的影響,只考慮物質(zhì)點體積的影響
ω(ξ,Vx′)=Vx′/VH
(21)
式中VH為近場域的體積。
通過鍵傳遞應(yīng)力計算近場力
(22)
式中nx x′為鍵的方向,
nx x′=ξ/‖ξ‖
(23)
將式(22)代入式(1)
(24)
采用了這種近場力計算方法的PD模型稱為鍵基對應(yīng)PD,在鍵基對應(yīng)PD中,不需要計算傳統(tǒng)鍵基PD的微觀材料常數(shù),直接利用物理方程求應(yīng)力,進(jìn)而計算近場力。
考慮受到單軸拉伸的正方形板,如圖4所示,板邊長為7 m,厚度為0.175 m,上下邊界加上固定位移ub=0.001 m,材料參數(shù)列入表1。采用顯式積分算法計算方程(24),根據(jù)文獻(xiàn)[11]計算臨界的時間步長為Δtc=8.1×10-5s,考慮加入的人工阻尼項對計算的影響,取小于限制值的時間步長為Δt=1×10-7s。
對比鍵基PD和鍵基對應(yīng)PD的計算結(jié)果,y向位移和x向位移分布如圖5所示。兩種PD模型計算的y方向位移分布均勻,但是x方向位移受到邊界效應(yīng)的影響,存在一定的誤差??紤]到鍵基泊松比限制的問題,將鍵基PD和鍵基對應(yīng)PD在y=-0.4375 m的x方向位移與理論解(u*=-νεyx,εy是施加邊界荷載后y方向產(chǎn)生的應(yīng)變,ν為泊松比)進(jìn)行對比,如圖6所示,相對誤差計算公式為‖uP D-u*‖/‖u*‖×100%??梢钥闯?,鍵基對應(yīng)PD比鍵基PD更接近于理論解,驗證了鍵基對應(yīng)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)
為驗證鍵基對應(yīng)PD模擬動態(tài)裂紋擴(kuò)展的能力,進(jìn)行Kalthoff-Winkler板沖擊試驗的數(shù)值模擬[13],考慮到原試驗沖擊塊形狀對平板動態(tài)裂紋擴(kuò)展影響較小,將問題簡化為二維,計算模型尺寸如 圖7 所示。根據(jù)文獻(xiàn)[14],采用對稱化建模,提高計算效率。材料參數(shù)和計算參數(shù)列入表2,網(wǎng)格劃分尺寸Δx=0.001 m,近場范圍尺寸δ=3.015Δx。PD斷裂準(zhǔn)則采用鍵斷裂準(zhǔn)則,模擬板受沖擊作用發(fā)生脆性破壞的現(xiàn)象。
圖6 鍵基PD與鍵基對應(yīng)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
裂紋擴(kuò)展的模擬結(jié)果如圖8所示,裂紋在初始裂紋尖端附近開始擴(kuò)展,之后裂紋會先沿著與荷載呈41°夾角的方向進(jìn)行擴(kuò)展,隨后沿著與荷載呈67°夾角的方向進(jìn)行擴(kuò)展,同時有裂紋分叉的現(xiàn)象出現(xiàn)。這與已有的數(shù)值模擬結(jié)果[14,16](70°)以及實驗結(jié)果[15](70°)比較接近。隨著荷載向下傳遞,損傷開始在梁的下端產(chǎn)生,之后向著預(yù)制裂紋所在的方向發(fā)展,這和文獻(xiàn)[14]的現(xiàn)象一致。
圖8 Kalthoff-Winkler板沖擊試驗裂紋擴(kuò)展隨時間變化
Fig.8 Crack path of Kalthoff-Winkler plate
本文結(jié)合鍵基PD和態(tài)基PD的優(yōu)點,提出了鍵基PD對應(yīng)模型,解決了鍵基PD泊松比限制的問題。新的鍵張量變形張量在小變形以及鍵長趨近無限小的條件下可以退化為經(jīng)典的變形張量。在此基礎(chǔ)上利用物理方程計算應(yīng)力,再通過應(yīng)力來計算近場力,建立了經(jīng)典連續(xù)介質(zhì)力學(xué)與鍵基PD的聯(lián)系。數(shù)值算例驗證了當(dāng)泊松比不取限制值時,鍵基對應(yīng)PD相對于傳統(tǒng)鍵基PD在預(yù)測平板x方向位移的正確性。對Kalthoff-Winkler板沖擊試驗的模擬驗證了鍵基對應(yīng)PD模擬動態(tài)裂紋擴(kuò)展的有效性。
鍵基對應(yīng)模型相比于傳統(tǒng)鍵基模型,計算成本稍高,對于三維問題中對鍵變形梯度的描述本文也只給出了公式而沒有作進(jìn)一步的討論,同時對于鍵變形的定義與傳統(tǒng)變形梯度定義之間的關(guān)系還有待進(jìn)一步的理論證明。