關(guān)于材料破壞的PD電算模擬

彭嘯秋

（河海大學力學與材料學院，江蘇南京 210098）

0 引言

固體材料的變形和破壞模擬一直是工程界關(guān)注的經(jīng)典難題，傳統(tǒng)理論使用以拉格朗日描述為基礎(chǔ)的偏微分方程推演本構(gòu)模型，其數(shù)值電算的典型代表即為有限單元法（FEM）。該方法在解決傳統(tǒng)的固體力學問題時非常出色，但在模擬開裂和空洞等大變形問題時往往因為網(wǎng)格變形過大或裂紋尖端奇異等原因?qū)е掠嬎憔认陆瞪踔辆W(wǎng)格失真，解決的方法之一是在計算前重新劃分網(wǎng)格，但這又導致運算程序繁瑣和計算效率低下。隨著斷裂力學和損傷力學的發(fā)展，人們試圖通過添加附加函數(shù)來改善這一困難，但在模擬材料破壞的全過程和微缺陷萌生、開裂延伸以及裂紋間的相互作用等復雜問題方面仍舊面臨挑戰(zhàn)。這一困境的根源在于基于傳統(tǒng)理論的物質(zhì)連續(xù)性假定與開裂衍生的不連續(xù)在本質(zhì)上不相容，偏微分方程所需要的空間導數(shù)在裂紋表面和尖端并不存在。因此，任何源于這些方程的數(shù)值方法在模擬開裂時都將存在這一困難。

作為改變這一現(xiàn)狀的嘗試，美國Sandia國家實驗室的S.A.Silling教授于2000年首先提出了一種不再需要對物體求空間導數(shù)值的新的固體力學理論——PD方法。其初衷是對傳統(tǒng)理論中最基本的數(shù)學描述進行重構(gòu)，將物體的連續(xù)性問題與開裂等不連續(xù)問題統(tǒng)一在相同的方程形式之下。這一基于非局部作用思想的理論使用空間積分方程而非傳統(tǒng)的偏微分方程對一系列離散的空間物質(zhì)點（material point）進行運動推演，從而避免了重新劃分網(wǎng)格和額外添加失效準則等麻煩?？梢园l(fā)現(xiàn)其在面對結(jié)構(gòu)畸變、材料破壞和介質(zhì)不連續(xù)等問題方面具有特殊的優(yōu)勢。

1 PD方法基本理論

PD理論可以視為宏觀尺度下的連續(xù)體分子動力學，根據(jù)牛頓第二定律，可以得到參考構(gòu)型中t時刻任意坐標x處物質(zhì)點的加速度公式：

該公式即為PD理論的基本運動方程，其中：Hx為x坐標的近場區(qū)域，u為位移矢量場，b為體力密度場，ρ為參考構(gòu)型中的密度函數(shù)，f為對偶力函數(shù)，其數(shù)值為物質(zhì)點x'作用于x上的力矢量。定義符號ξ和η分別為參考構(gòu)型中兩物質(zhì)點間的相對位置矢量和相對位移矢量：

則ξ＋η為顆粒間的現(xiàn)時相對位置矢量。x'和x直接的物理相互作用稱為“對偶鍵”（pairwise bond），其作用類似于一根“空間彈簧”。PD理論規(guī)定：物質(zhì)點僅在某個相對作用范圍“δ”（horizon）內(nèi)與臨近物質(zhì)點相互作用，該范圍稱為“近場范圍”，也就是說，當顆粒x'在“近場范圍”之外時將不能夠被“看見”：

公式（1）中的Hx即為以x為中心以δ為半徑的臨近空間鄰域，稱為“近場臨域”：

在“對偶鍵”級別上引入破壞的優(yōu)勢就在于它可以對一點處的局部破壞進行精確描述，而將破壞引入連續(xù)體模型最簡單的方法就是當“對偶鍵”超過預先設(shè)定的拉伸長度時允許它們發(fā)生不可逆的斷裂（受損的“鍵”不可“修復”，也不會生成新“鍵”），其中的張力也將不可持續(xù)?！版I”的拉伸破壞定義如下：

（圖1）

其中，x為包含在μ內(nèi)的函數(shù)自變量，以暗示μ為物體的位置函數(shù)。注意到0≤φ≤1，數(shù)值0代表材料的初始狀態(tài)，數(shù)值1代表一個物質(zhì)點已經(jīng)與所有曾與它關(guān)聯(lián)的點斷開聯(lián)系。因為斷裂的“對偶鍵”無法再承受張力荷載，物體中某些“鍵”的斷裂可能導致其他“對偶鍵”斷裂的連鎖反應，而這又將導致材料的軟化響應，直到多條裂紋接合在一起并形成斷裂面。為了簡便起見，設(shè)定上式中的μ為具有歷史依賴性的標量函數(shù)：

此處，s0為“對偶鍵”破壞的臨界拉伸長度，定義如下：

其中s00和α均為常數(shù)，smin為現(xiàn)時與某物質(zhì)點相關(guān)聯(lián)的所有“對偶鍵”的最小拉伸長度。

2 電算數(shù)值建模

將參考構(gòu)型中已知體積的區(qū)域離散為空間節(jié)點，而所有的節(jié)點粒子在空間中依次排列，形成陣列。為了簡便起見，以下討論使用均質(zhì)材料和線性PD形式。離散公式(1)即可得到運動積分方程的有限和形式：

n為時間步數(shù)，下標為各節(jié)點編號，Vp為節(jié)點p的體積。其中f由下式提供：

其中f為標量形式的對偶力狀態(tài)函數(shù)。加速度值采用中心差分形式：

計算時設(shè)定時間步長Δt為常數(shù)。定義“初始微彈脆性”(Prototype Micro-elastic Brittle, PMB)材料，其受力狀態(tài)場函數(shù)為：

彈性常數(shù)c和“對偶鍵”拉伸長度s為：

其中K為體積模量，δ為近場范圍。設(shè)標量函數(shù)f僅依賴于“對偶鍵”的拉伸長度，因為f并不依賴于方向系數(shù)，表明這種材料為“各向同性”。PMB模型僅有一個屬性參數(shù)即體積模量K，這是因為PMB模型預先設(shè)定了泊松比（三維模型中為1/4，二維模型中為1/3）。在本文中，PMB材料在其初始狀態(tài)下為各向同性，而隨著一些特定方向“對偶鍵”的斷裂，將導致隨后材料響應的各向異性。

3 電算程序流程圖

見（表 1）

4 數(shù)值算例

算例1：材料裂紋的擴展

數(shù)據(jù)：材料尺寸為5×5m，材料密度ρ=300kg/m3，體積模量K=19.2GPa，單元顆粒尺度Δx=5×10-4m，近場范圍δ=1.8×10-3m，s0=5×10-2，邊界初始速度狀態(tài)V=15m/s。圖片時間間隔為6.0×10-3s。見（圖 2）

算例2：球體偏心炸裂

數(shù)據(jù)：爆炸沖擊速度V=300m/s，球體半徑為0.02m，密度ρ=3000kg/m3，體積模量 K=15GPa，s00=5×10-3，α=0.25，單元顆粒尺度Δx=5×10-4m，近場范圍δ=1.5×10-3m，圖片時間間隔為 3.0×10-5s。見（圖 3）

以上電算算例表面，在使用PD構(gòu)件的破壞模型中，材料的開裂和破壞自然發(fā)生，顯示了相較與傳統(tǒng)理論，PD方法在處理不連續(xù)問題時的特殊優(yōu)勢。

5 結(jié)語

PD方法對于開裂模型的主要優(yōu)勢在于當某個裂紋開展時不需要額外補充任何有關(guān)它的速度、方向、分叉、不穩(wěn)定性以及延遲性的相關(guān)條件。所有這些現(xiàn)象的出現(xiàn)都是運動方程和包含了“對偶鍵”數(shù)量級破壞的本構(gòu)模型共同作用的結(jié)果。而且它也不需要對開裂尖端的位置延伸保持跟蹤和關(guān)注，另一個方面，PD理論關(guān)于本構(gòu)模型的模擬中f對ξ的依賴性有別于傳統(tǒng)理論模型。因為這些原因，我們可以計算任意數(shù)量的裂紋，包括它們之間的相互作用。

有些無網(wǎng)格方法（如SPH）的提出被視為是專門為了處理破壞模型。而現(xiàn)在提及的PD方法和這些無網(wǎng)格方法之間確實存在著差別：PD模型中裂紋的自發(fā)產(chǎn)生是基于成對節(jié)點之間的相互作用和運動方程運算的自然結(jié)果。因此不需要統(tǒng)計大量的鄰近節(jié)點來確定空間系數(shù)。這種對偶作用的后果之一就是PD方法不會受到類似SPH中“不穩(wěn)定張力”的影響。

作為一種全新的方法，PD理論仍有諸多不足和有待完善的地方，如嚴格的收斂準則和計算誤差估計，模型的自適應性和均質(zhì)化，與傳統(tǒng)的彈、塑、粘性材料模型的關(guān)聯(lián)等，相信隨著研究的深入，PD理論將在工程材料和結(jié)構(gòu)破壞方面得到更廣泛的運用。

（圖2）

（圖3）

[1] 龍馭球. 有限元法概論[M]. 北京：人民教育出版社.1978

[2] 李臥東,王元漢,陳曉波. 無網(wǎng)格法在斷裂力學中的應用[J]. 巖石力學與工程學報, 2001, 20(4)：462~466

[3] Silling S A. Reformulation of elasticity theory for discontinuities and long-range forces[J].Journal of the Mechanics and Physics of Solids,2000, 48(1)：175~209

[4] 黃丹,章青,喬丕忠,沈峰. 近場動力學方法及其應用[J]. 力學進展, 2010, 40(4)：448~459

[5] Gerstle W, Sau N. Peridynamic modeling of plain and reinforced concrete structures. In：18th International Conference on Structural Mechanics in Reactor Technology(SMiRT 18), Bejing, China,2005

[6] Macek R W, Silling S A. Peridynamics via finite element analysis[J]. Finite Elements in Analysis and Design, 2007, 43(15)：1169~1178

[7] Silling S A, Epton M,Weckner O, et al. Peridynamic states and constitutive modeling[J]. Journal of Elasticity, 2007,88(2)：151~184

[8] Silling S A. Linearized theory of peridynamic states. Sandia National Laboratory Report,SAND2009-2458, Albuquerque, New Mexico, 2009

[9] Park M L, Lehoucq R B, Plimpton S J et al.Implementing peridnamics within a molecular dynamic code[J]. Computer Physics Communications,2008,179(11)：777~783

[10] Gerstle W, Sau N, Silling S A. Peridynamic modeling of concrete structures[J]. Nuclear Engineering and Design, 2007, 237(12-13)：1250~1258

[11] Silling S A, Askari E. A meshfree method based on the peridynamic model of solid mechanics[J].Computer and Structures, 2005, 83(17-18)：1526~1535

[12] Killc B, Madenci E. Prediction of crack paths in a quenched glass plate by using peridynamic theory[J]. International Journal of Fracture,2009, 156(2)：165~177

[13] 劉謀斌,宗智,常建忠. 光滑粒子動力學方法的發(fā)展與應用[J]. 力學進展, 2011, 41(2)：217~234