關于材料破壞的PD電算模擬

彭嘯秋

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

0 引言

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

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

1 PD方法基本理論

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

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

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

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

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

（圖1）

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

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

其中s00和α均為常數，smin為現時與某物質點相關聯的所有“對偶鍵”的最小拉伸長度。

2 電算數值建模

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

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

其中f為標量形式的對偶力狀態函數。加速度值采用中心差分形式：

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

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

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

3 電算程序流程圖

見（表 1）

4 數值算例

算例1：材料裂紋的擴展

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

算例2：球體偏心炸裂

數據：爆炸沖擊速度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構件的破壞模型中，材料的開裂和破壞自然發生，顯示了相較與傳統理論，PD方法在處理不連續問題時的特殊優勢。

5 結語

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

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

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

（圖2）

（圖3）

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

[2] 李臥東,王元漢,陳曉波. 無網格法在斷裂力學中的應用[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] 劉謀斌,宗智,常建忠. 光滑粒子動力學方法的發展與應用[J]. 力學進展, 2011, 41(2)：217~234