應力波與缺陷相互作用的宏觀微觀數值模擬＊

郭昭亮，任國武，湯鐵鋼，劉倉理

（中國工程物理研究院流體物理研究所沖擊波物理與爆轟物理重點實驗室，四川 綿陽621999）

在材料斷裂機理研究以及工程應用中，裂紋誘致的斷裂行為，一直都是力學家和工程師關注的焦點。尤其對金屬及其合金等延性材料，裂紋尖端或刻槽尖端附近的塑性變形區，在裂紋的起裂、擴展、止裂過程中起著重要的作用［1－4］。對塑性區形狀大小的估計歷來已久，傳統的做法從宏觀斷裂力學角度按照各種屈服準則確定其形狀。G.R.Irwin估計了理想裂紋尖端塑性區的大小、形狀，并作為彈塑性等效裂紋長度的修正［5］。張亞等［6］針對復合型裂紋，基于小范圍屈服準則，采用俞茂宏統一強度理論，給出了裂尖塑性區的統一解析解。X.Gao等［4］基于小范圍屈服條件，采用Hill屈服準則計算了不同加載情況下的裂尖塑性區狀，發現塑性區呈現蝴蝶狀，且彈塑性區邊界光滑。Y.Huang等［7］利用最大裂紋張開位移估計了塑性區尺寸大小。然而對于真實材料，由于不存在理想的數學裂紋，則可以使用一定幾何形狀的缺陷近似裂紋，缺陷局域的塑性區大小、形狀可以通過計算缺陷處的應力集中，并依據強度理論進行估計。S.Q.Shi等［8］使用簡單的模型計算了刻槽尖端的最大正應力與塑性區的大小。張培源等［9］探討了裂紋尖端為圓弧形的鈍化模型，并依據Tresca屈服條件，得到了鈍化裂紋前緣塑性區的線場分析解，對于含徑向裂紋和圓弧鈍化區的圓盤外圍，得到其線彈性解，形成了一套估計裂紋前端塑性區尺寸的方法。現有的研究指出，對于平面應變與平面應力情況，塑性區的大小形狀以及尖端的正應力分布不同［10］。

在處理裂紋局域的塑性區時，理論大都基于準靜態假定，并且認為塑性區首先出現在裂尖（或缺陷）頂端的邊界。錢才富等［11－12］從微觀斷裂力學的角度更精細地計算了裂紋前緣塑性區和無位錯區，發現應力最大的位置出現在遠離裂尖一定距離處，并指出它們的存在及性質決定著裂紋的擴展行為。

然而，對于由動態載荷誘致的缺陷局域的塑性區形成與演化，以及由此誘發的裂紋萌生、擴展等現象的研究相對較少，研究主要集中在孔洞的長大與演化過程。R.E.Rudd等［13－14］、E.T.Sepp?l?等［15］采用分子動力學程序，模擬了韌性金屬中的孔洞成核、長大，以及孔洞之間的相互作用、貫通。祁美蘭等［16］、王永剛等［17］采用二維LS－DYNA程序，研究了平板撞擊加載下含初始雜質的純鋁樣品中微孔洞的成核、長大與閉合現象。在快速載荷作用下，由于應力波與材料缺陷的相互作用，將出現與準靜態加載時完全不同的物理現象。本文中，利用LS－DYNA 3D有限元程序以及分子動力學方法，從宏觀微觀兩個層次分別模擬含有預置圓形（橢圓形）缺陷的薄板，在動態拉伸載荷作用下，塑性區的形成、演化過程以及隨之而來的裂紋動態擴展過程。

1 計算模型

1.1 宏觀模型

宏觀尺度上，采用LS－DYNA有限元程序，模擬含有預置缺陷的薄板（見圖1）在動態拉伸載荷（見圖2）作用下，缺陷與應力波相互作用對塑性區形成演化的影響，以及由此誘發的裂紋萌生、裂紋動態擴展行為，并通過調節預置缺陷的長軸與短軸比例實現不同的缺陷預置。假定材料動態應力應變曲線如圖3所示。

圖1 建模示意圖Fig.1 Scheme of model

圖2 拉伸載荷Fig.2 History of load

圖3 應力應變曲線Fig.3 Stress－strain curve

對于同種材料，初始載荷σ0及缺陷的形狀位置，決定了薄板中塑性區的形成位置與變形失效行為。

（1）當σ0＜Y／2時，薄板加載處發射彈性波，在遠離缺陷處的應力波相互作用不會形成超過屈服應力的區域，只有在缺陷邊界處的反射會形成塑性區。

（2）當Y／2＜σ0＜Y時，薄板加載處發射彈性波，在遠離缺陷處的應力波相互作用區由于應力超過屈服應力將形成帶狀塑性區，此外由于缺陷邊界處的反射也會形成塑性區。

（3）當σ0＞Y時，加載邊界首先發射彈性前驅波，其后跟隨塑性波，塑性波掃過的區域全部進入塑性，彈性前驅波由于相互作用及在缺陷附近的反射也會在短期內形成塑性區，當塑性波掃過后，將整體進入塑性。

本文中的計算主要側重于第2種載荷情況。薄板材料為20鋼，密度ρ0＝7.85g／cm3，剪切模量G＝82GPa，彈性模量E＝210GPa，泊松比ν＝0.286［18］。采用流體彈塑性本構模型，Grüneisen狀態方程［19］；屈服應力取490MPa，為準靜態屈服應力的2倍，塑性硬化模量取2GPa，為彈性模量的約1%［20］。

1.2 微觀模型

微觀尺度上，采用分子動力學方法模擬橢圓微缺陷在拉伸下的動態裂紋擴展。通過數值求解牛頓方程，獲得整個體系粒子的運動軌跡。采用Verlet速度算法，多粒子相互作用勢為Lennard－Jones勢，并采用多項式截斷［21］：

式中：η、ξ分別是能量尺度和長度尺度，參數a2、a3及rmax可通過式（1）中兩個函數在rspl處的連續、包括其一階和二階導數連續，來共同決定，分別是a2＝0.542 449 4、a3＝0.093 505 27、rspl＝1.244 455和rmax＝1.711 238。由此勢函數，可知其能量最低點對應的平衡位置a0＝21／6ξ。

計算中，采用恒定位移加載方法，即在每個時間步，都對最外層上下邊界的兩個原子層分別施加一個固定的位移增量Δu，這樣可避免邊界復雜的應力狀態，其中Δu＝0.000 6。同時為避免外邊界帶來的波反射，在垂直于裂紋傳播方向考慮加入線性速度梯度，減少模擬的時間。整個模擬研究將限制于拉伸型裂紋，其加載為沿x方向。

依據維里定理，在原子尺度上對i原子，其局域應力狀態可定義為：

式中：α和β是方向，V是原子的體積，對于二維系統而言，是原子所占據的面積。

2 LS－DYNA計算結果

2.1 含偏心圓孔缺陷非對稱薄板塑性區演化過程

圖4為含有偏心圓孔缺陷的薄板，薄板尺寸為50mm×75mm×1mm，缺陷中心位于（25mm，25mm），圓孔直徑10mm。在薄板兩側邊界同時施加如圖1～2所示的拉伸型階躍載荷400MPa，并持續到計算結束。算例中，所施加的載荷滿足Y／2＜σ0＜Y。

圖4 含圓孔缺陷非對稱薄板網格劃分Fig.4 Mesh of an asymmetry sheet with a circularity hole

關注物理量等效塑性應變，可以得到應力波與缺陷相互作用時，塑性區隨時間的演化過程，如圖5所示。在4.5μs時，由于左側應力波與缺陷作用，而右側應力波尚未到達缺陷邊界，所以只形成2個小的塑性變形區。由此可見，缺陷邊界處的塑性區形成與整體應力是否平衡無關，它們是應力波與缺陷相互作用、反射塑性波導致的結果。隨著時間演化，在7.3μs時，左右兩側的應力波在空間遠離缺陷的位置發生相互作用，形成塑性帶，這是應力波在空間相互作用導致的結果。隨后，此塑性帶逐步演化，最終與圓孔附近的塑性區匯合，在空間形成整體的塑性區域。

圖5 含偏心圓孔缺陷非對稱薄板塑性區演化Fig.5 Plastic zone evolution in an asymmetry sheet with a circularity hole

2.2 含中心橢圓孔缺陷對稱薄板塑性區演化及失效過程

圖6為含有橢圓孔缺陷的薄板，薄板尺寸為50mm×50mm×1mm，橢圓孔長軸10mm，短軸4mm，中心位置在（25mm，25mm），加載形式如圖1～2所示。

圖6 含橢圓孔缺陷對稱薄板網格劃分Fig.6 Mesh of a symmetry sheet with an elliptic hole

當塑性應變累積到一定程度，局部將會出現失效。為此，在程序中設置了失效應變0.04。應力波與缺陷相互作用后，塑性區域隨時間的演化過程以及由應變累積導致的局部失效行為，如圖7所示。

在4.4μs時，橢圓邊緣同時形成4個小的塑性變形區，隨著應力波的相互作用，塑性區迅速形成，在7.0μs時，塑性應變最大值的位置已經上移到橢圓邊界外。由此誘致的裂紋萌生首先出現在遠離缺陷邊界，然后裂紋與缺陷貫通并向外擴展。有趣的是，隨著裂紋的擴展，由于卸載波的相互作用，在與加載平行的方向，將出現裂紋的萌生、擴展，最終形成十字架形的斷裂形式。塑性區的形成演化、裂紋的萌生擴展都與應力波和缺陷的相互作用密切相關，它展現了動態載荷作用下獨特的變形與失效模式。

圖7 含中心橢圓孔缺陷對稱薄板塑性區演化及失效過程Fig.7 Plastic zone evolution and failure process in a symmetry sheet with an elliptic hole

3 分子動力學計算結果

對于微觀情況，所選用模擬體系的幾何形狀是一個包含橢圓裂紋的二維密堆積體系，其長度沿x方向是lx＝600a0，沿y方向是ly＝693a0，而橢圓裂紋尺寸是沿x方向的長軸為30ξ，沿y方向的短軸為10ξ。總模擬原子數目為約480000。模擬采用微正則系綜（粒子數、體積和總能守恒），此系綜適合研究動態裂紋的非平衡過程。初始的模擬體系溫度趨近于零。在實際模擬中，將長度參數ξ、能量參數η和粒子質量m都設置為參考量，那么所有量的量綱均為1。由德拜頻率，時間步長選擇為Δt＝0.003t0，t0＝ξ－1（η／m）1／2，每隔1000Δt獲得一幅圖像。

計算結果如圖8所示。在計算至99t0時，遠離缺陷上下邊緣的位置出現位錯堆積，在隨后的6t0內，局部空洞化，由此誘致裂紋的萌生。裂紋沿著縱向擴展，與缺陷邊緣貫穿，最終導致整體的斷裂失效。在微觀模擬中，即使不考慮無位錯區，含有預置缺陷的薄板在動態拉伸下應力最大的位置也將出現在遠離缺陷邊緣，這與宏觀模擬（見圖7）類似。分子動力學模擬中，缺陷橫向邊緣附近已經存在了相當數量的位錯，并出現了應變局域化現象。但由于系統提供的能量不足以使這些位錯發展為裂紋，因而沒有得到圖7中所示十字架形的斷裂形式。

圖8 含有橢圓缺陷的二維密堆積體系的分子動力學模擬Fig.8 Molecular dynamic simulation in 2Dclosed－pack system with an elliptic hole

4 討 論

宏觀和微觀兩個尺度的計算表明，在動態載荷作用下探討含缺陷結構的變形與失效行為時，必須考慮應力波與缺陷相互作用所帶來的特有現象。它與含缺陷結構在準靜態載荷作用下的變形與失效行為有著本質的區別。

在準靜態加載下，不考慮結構中的應力波傳播所帶來的影響，認為材料內部任何時刻都處在應力平衡態。

圖9 材料內部缺陷分布示意圖Fig.9 Flaw distribution within material

如圖9所示，Ω為關注的材料區域，?Ω為該區域的外邊界，內部存在分布缺陷，可用內邊界?Ωi、?Ωi＋1、?Ωi＋2等效。準靜態加載下，應力將實時在區域Ω內達到平衡，并形成相應的應力分布。考慮缺陷的幾何形狀，那么就不可能出現所謂理想的尖端，也就不會出現應力奇點，此時應力分布為解析函數。根據最大模原理［22］，應力絕對值的最大值只可能出現在邊界上，也即是說，如果結構內部存在應力極大值，則必然分布于缺陷的邊界上。

需要指出的是，這里的結論只適合于外加載荷在材料內部尚未形成塑性區的情況，因為一旦塑性區形成，應力極大值就會向材料內部移動。于是在準靜態假設下，塑性區只可能首先產生在邊界，并隨著載荷的增加由邊界向材料內部擴展。

動態載荷作用下的情況則不同，區域內的應力并非實時達到平衡，考慮應力波的傳播，應力波前緣應力值會出現跳躍，這就破壞了區域內的應力可導性，此時應力函數非解析，不能使用最大模原理理解應力極大值的位置。在發生應力波相互作用后，以及應力波與缺陷相互作用后，局域應力非平衡時，材料內部的應力狀態非常復雜，很難用解析方式對應力極大值的位置給出合適的描述，尤其是對于存在不規則幾何的缺陷，因此選擇LS－DYNA程序以及分子動力學方法通過數值模擬給出其發展過程。

應力波與缺陷相互作用時，可以將內部缺陷視為內邊界，應力波在內邊界上進行反射，反射波與加載波的相互作用形成塑性區。此塑性區必然形成于遠離缺陷邊緣的位置，并由此誘導裂紋的萌生，引發裂紋的動態擴展。考慮應力波與缺陷的相互作用對塑性區形成演變以及材料失效行為的影響，雖然目前只進行了初期的數值模擬，卻得到了完全不同于準靜態分析的有趣現象。這可以為動態載荷作用下結構的防護提供相應的支持。

5 結 論

在宏觀上利用LS－DYNA有限元程序，在微觀上使用分子動力學程序，分別模擬了含有預置缺陷的薄板在動態拉伸載荷作用下的塑性區的形成演化以及裂紋的萌生與擴展過程。主要結論如下：

（1）動載作用下，從本質上講塑性區是塑性波掃過的區域，它可以由應力波在缺陷處的反射以及應力波在空間的相互作用形成。

（2）宏觀LS－DYNA程序模擬與微觀分子動力學模擬共同指出：在動態載荷作用下，即使不考慮無位錯區的存在，裂紋萌生的位置也可能出現在遠離缺陷邊緣的位置。從數學的角度看，這個現象的出現是由于動載作用下，應力波與缺陷的相互作用改變了應力函數的解析性。這個改變直接導致應力極大值空間分布不再保持準靜態載荷作用下的邊界分布特征。對于真實材料在動態載荷作用下缺陷附近的裂紋萌生，建議同時考慮無位錯區的存在以及應力波效應帶來的影響。

（3）由于裂紋擴展導致的卸載波相互作用，將在與加載平行的方向形成新的塑性區，誘導新的裂紋萌生，產生類十字架形的失效形式。這體現了應力波相互作用在裂紋萌生過程中所產生的重要影響。

［1］Fields R J，de Wit R.Plastic zone formation around an arresting crack［C］∥Knauss W G，Rosakis A J.Non－linear fracture：Recent advances.Springer，1990：231－238.

［2］Khan S M A，Khraisheh M K.A new criterion for mixed mode fracture initiation based on the crack tip plastic core region［J］.International Journal of Plasticity，2004，20（1）：55－84.

［3］Bian Li－chun，Kim K S.The minimum plastic zone radius criterion for crack initiation direction applied to surface cracks and through－cracks under mixed model loading［J］.International Journal of Fatigue，2004，26（11）：1169－1178.

［4］Gao Xin，Wang Han－gong，Kang Xing－wu，et al.Analysis solutions to crack tip plastic zone under various loading conditions［J］.European Journal of Mechanics A：Solids，2010，29（4）：738－745.

［5］Janssen M，Zuidema J，Wanhill R J H.Fracture Mechanics［M］.Spon Press，2004.

［6］張亞，強洪夫，楊月誠.復合型裂紋小范圍屈服下裂尖塑性區統一解［J］.機械工程學報，2007，43（2）：50－54.Zhang Ya，Qiang Hong－fu，Yang Yue－cheng.Unified solutions to mixed mode crack tip under small scale yielding［J］.Chinese Journal of Mechanical Engineering，2007，43（2）：50－54.

［7］Huang Yi，Chen Jing－jie，Liu Gang.A new method of plastic zone size determined based on maximum crack opening displacement［J］.Engineering Fracture Mechanics，2010，77（14）：2912－2918.

［8］Shi S Q，Puls M P.A simple method of estimating the maximum normal stress and plastic zone size at a shallow notch［J］.International Journal of Pressure Vessels and Piping，1995，64（1）：67－71.

［9］張培源，張曉敏，嚴波，等.裂尖曲率對裂紋前緣塑性區的影響［J］.應用力學學報，2004，21（4）：93－96.Zhang Pei－yuan，Zhang Xiao－min，Yan Bo，et al.Plastic zone affected by crack tip curvature［J］.Chinese Journal of Applied Mechanics，2004，21（4）：93－96.

［10］布洛克 D.工程斷裂力學［M］.王克仁，譯.北京：科學出版社，1980.

［11］錢才富，姜忠軍，陳平，等，裂紋尖端塑性區和無位錯區的微觀模擬［J］.金屬學報，2004，40（2）：159－162.Qian Cai－fu，Jiang Zhong－jun，Chen Ping，et al.Micro－simulation of crack tip plastic zone and dislocation－free zone［J］.Acta Metallurgica Sinica，2004，40（2）：159－162.

［12］錢才富，李慧芳，崔文勇.Ⅰ型裂紋尖端塑性區和無位錯區及其對裂紋擴展的影響［J］.材料研究學報，2007，21（6）：599－603.Qian Cai－fu，Li Hui－fang，Cui Wen－yong.ModeⅠcrack tip plastic zone，dislocation－free zone and their effects on crack propagation［J］.Chinese Journal of Materials Research，2007，21（6）：599－603.

［13］Rudd R E，Belak J F.Void nucleation and associated plasticity in dynamic fracture of polycrystalline copper：An atomistic simulation［J］.Computational Materials Science，2002，24（1／2）：148－153.

［14］Rudd R E.Void growth in BCC metals simulated with molecular dynamics using the Finnis－Sinclair potential［J］.Philosophical Magazine，2009，89（34／35／36）：3133－3161.

［15］Sepp?l?E T，Belak J，Rudd R E.Onset of void coalescence during dynamic fracture of ductile metals［J］.Physical Review Letters，2004，93（24）：245503.

［16］祁美蘭，賀宏亮，王永剛，等.動態沖擊下純鋁中微損傷演化的仿真研究［J］.振動與沖擊，2007，26（8）：133－135.Qi Mei－lan，He Hong－liang，Wang Yong－gang，et al.Simulation of micro void evolution in the pure aluminum under dynamic loading［J］.Journal of Vibration and Shock，2007，26（8）：133－135.

［17］王永剛，劉宏偉.強沖擊載荷下含雜質的純鋁中微孔洞長大的動力學行為［J］.高壓物理學報，2010，24（4）：248－254.Wang Yong－Gang，Liu Hong－wei.Dynamic behavior of void growth in aluminum with a preexisting flaw under intense impact loading［J］.Chinese Journal of High Pressure Physics，2010，24（4）：248－254.

［18］工程材料實用手冊委員會.工程材料實用手冊：結構鋼、不銹鋼［M］.北京：中國標準出版社，1988.

［19］LS－DYNA Keyword user’s manual－2003［M］.California：Livermore，2003.

［20］郭昭亮，劉倉理，湯鐵鋼.預置圓孔膨脹環動態斷裂行為研究［J］.實驗力學，2010，25（5）：546－552.Guo Zhao－liang，Liu Cang－li，Tang Tie－gang.On the expanding ring dynamic fracture behavior with a preset circular hole［J］.Journal of Experimental Mechanics，2010，25（5）：546－552.

［21］Wagner N J，Holian B L，Voter A F.Molecular－dynamics simulations of two－dimensional materials at high strain rates［J］.Physical Review A，1992，45（12）：8457－8470.

［22］潘永亮，汪琥庭，汪芳庭，等.復變函數［M］.北京：科學出版社，2004.