内聚力单元在Voronoi元胞的晶间断裂模拟中的应用*

张 帆，赵建平

(南京工业大学机械与动力工程学院，南京 211816)

内聚力单元在Voronoi元胞的晶间断裂模拟中的应用*

张 帆，赵建平

(南京工业大学机械与动力工程学院，南京 211816)

应用二维Voronoi算法进行晶粒的几何建模，采用粘聚力模型模拟晶界的分离现象。以多晶铝为研究对象，研究了无缺陷条件下，施加单向拉伸载荷时，晶界的损伤和失效，以及整个模型的应力分布，并与有缺陷、单向拉伸时，晶界的损伤、应力的分布进行了对比研究，发现在单向拉伸加载的情况下，垂直于拉伸方向，最靠近加载边界的一条晶界带上出现了连续的失效单元，形成晶间断裂“裂纹”。模拟结果表明，随着单向拉伸载荷的增大，多晶铝晶界的损伤会加剧;晶界带上缺陷的数量与分布影响了整个模型的应力分布，因此不同内聚力单元的应力值随着载荷的增加呈现出不同的变化规律。

多晶铝;晶间断裂;晶界损伤

0 引言

为了分析工程结构的耐用性和寿命周期，通常采用唯象的损伤断裂模型［1］，对其宏观破坏行为进行数值模拟，然而，这种宏观同质化［2］的处理方法忽略了导致裂纹萌生与扩展的真正的物理机制。为了综合考虑这些物理机制的影响，必须进行介观或微观尺度下的数值模拟。

在介观尺度下，多晶体金属材料的裂纹主要沿晶界扩展［1］。Luther等［1，3］研究了单向静态加载、不定义初始损伤的条件下，多晶铝介观尺度的代表性体积单元［4］(RVE)的裂纹萌生与扩展;司良英等［5］采用晶体塑性有限元方法(CPFEM)对多晶铝在拉伸试验中的变形行为进行了数值模拟;Iesulauro等［6］应用FRANC2D/L软件，研究了多晶铝晶粒构型、材料模型、晶粒与晶界参数对失效路径、应力分布、失效应变的影响;文献［4］中运用AutoRVE软件，对各向异性多晶铝在对称边界条件下作静力拉伸时的应力状况进行了数值模拟。然而，在晶界上设置“缺陷”对应力分布与晶界损伤的影响研究还鲜见报道。

鉴于此，本文将分别在晶界有、无缺陷的情况下，对多晶铝单向拉伸进行不同载荷大小下的数值模拟，探究载荷与缺陷对应力分布与晶界损伤的影响。

1 基本理论

1.1 模拟晶界的内聚力模型

内聚力模型的核心是用局部的应力——位移之间的关系来描述损伤或断裂过程［7］，即内聚区的本构关系，称为 Traction Separation Law(TSL)。在ABAQUS软件中，内聚区采用一层内聚力单元表示［8］。本文在计算中采用具有线性软化的混合损伤模式［9］下的TSL，其一维形式如图1所示，其中T和δ分别代表界面的内聚力和相对位移。从图中可以看到，在内聚区开始承载时，内聚力单元的初始刚度值［10-11］为Kn，当内聚力达到内聚力强度Tl时界面相对位移达到δm0，此时内聚力单元开始损伤，(1-D)Kn表示材料损伤后的实际刚度;随着外力的进一步增大，界面相对位移增大，内聚力则线性降低直至减小为0，相应地界面相对位移达到临界值δfm，此时内聚力单元完全失效，该处的断裂能达到临界值Gc，Gc的值为 TSL 曲线所包围的面积［8，12］。

图1 具有线性软化的混合损伤模式下的TSL

对于线性软化，ABAQUS采用Camanho和Davila提出的损伤演化变量D定义失效退化［13-14］:

其中δmaxm为加载过程中界面达到的相对位移的最大值。

损伤变量D的初始值为0，一旦损伤开始，在后续载荷的作用下，该变量由0逐渐增长到1。当D=1时，表示内聚力单元完全失效，此时单元刚度完全从模型总体刚度矩阵中删除，意味着此处发生了断裂［10-11，15-16］。

1.2 界面元的混合型失效准则

1.2.1 降阶准则

复杂结构中常见的裂纹形式是混合型裂纹，相应的有初始失效判断准则及失效扩展判断准则［12］。本文采用二次准则［10-11，15-16］作为判断界面元的刚度降阶的失效准则:

式中 σn、τs、τt分别为内聚力单元法向、切向一、切向二的应力;T、S分别为拉伸、剪切强度。

1.2.2 失效准则

在混合模式加载下，失效退化进程判断以能量准则作为判断混合模式失效的标准［15］。本文采用幂指数准则［10-11，15-16］作为判断混合断裂模式中界面元最终破坏的失效准则:

式中，Gic(i=Ⅰ，Ⅱ，Ⅲ)分别代表了在纯i型载荷下，材料单位长度的临界应变能释放率。

2 计算模型与材料参数

2.1 模型分析与假设

多晶铝晶间断裂的实际过程非常复杂，其断裂模式不仅有晶间断裂，还有穿晶断裂。另外，实际的晶粒结构与材料参数很难确切定义，这些因素就给建立与材料介观结构完全一致的模型并考虑所有可能的失效模式，从而进行精确的有限元计算带来相当大的困难［6，9］。

考虑到本研究的主要目的是分析载荷与缺陷对多晶铝代表性体积单元的应力分布与晶界损伤的影响，在有限元建模与计算过程中作了如下假设:

(1)晶粒内部各向同性，各个晶粒具有相同的材料属性。

(2)所有晶界带具有相同的材料参数。

(3)晶粒内部无缺陷，缺陷只存在于晶界带上。

(4)晶间断裂是多晶铝唯一的断裂失效模式。

2.2 有限元模型

本文利用Voronoi tessellation方法［17-19］生成了具有随机晶粒几何形状的49晶粒Voronoi元胞模型。计算模型选用的材料为AA7075多晶铝，平均晶粒度［6］为50 μm，整个Voronoi元胞模型的长和宽均为350μm。

在进行有限元网格划分时，晶粒内部采用CPS4R单元，晶界带采用一层COH2D4单元，共划分26529个单元。

为了考察载荷大小对多晶铝代表性体积单元的应力分布与晶界损伤的影响，在边界BC-1上施加单向拉伸载荷U1，在边界BC-2上设置边界条件UX=0及UY=0如图2所示，其中UX、UY分别表示X轴、Y轴方向的位移。

图2 Voronoi元胞模型和边界条件示意图

2.3 材料参数

对于 AA7075 多晶铝，化学成分［6，20］(质量分数，%)为:Zn 5.6，Cu 1.6，Mg 2.5，Cr 0.25，Al余量。晶粒内部弹性模量E=72000Mpa，泊松比μ=0.33，晶粒内部的塑性硬化则根据图3所示的硬化曲线［6］，通过GetData软件捕捉曲线上的数据点，定义屈服应力与塑性应变。晶界的计算参数［1，6］为Tl=500Mpa，Gc=0.15N/mm，Kn=4e7Mpa。

图3 晶粒内部的塑性硬化曲线

2.4 有限元计算流程

有限元分析一般包括以下三个步骤:前处理、分析计算和后处理。这些步骤在ABAQUS软件中是通过一系列功能模块完成的［21］。本文采用的有限元计算过程如图4所示。

图4 有限元计算流程图

3 计算结果与讨论

3.1 缺陷率对应力值的影响

在不含缺陷、载荷U1=0.005 mm的情况下，晶界带损伤模拟结果如图5所示。

由图5可知，垂直于拉伸方向的一条晶界带出现了连续的失效单元，形成晶间断裂“裂纹”。因此，在该晶界带上分别随机选取50%、40%、30%、20%、10%的单元作为“缺陷”，选取晶界带上的24502、26509号单元作为研究对象，考察含缺陷情况下应力值、损伤值的变化规律，并与无缺陷情况下的计算结果进行对比。

单个内聚力单元的应力值随载荷变化规律:在载荷 U1分别为 0.001 mm、0.002 mm、0.003 mm、0.004 mm、0.005 mm的情况下，单元 24502、26509的应力值随U1的变化曲线如图6所示。从图6a中可以发现，对于单元24502来说，当缺陷率为0%、10%、30%时，随着载荷的增加，单元积分点的应力值逐渐下降，之后突降为0然后保持恒定，这是由于载荷的增加降低了内聚力单元的承载能力，直至单元完全失效被删除，从而失去了承载能力;当缺陷率为20%、40%、50%时，随着载荷的增加，单元积分点的应力值相应地降低至某个极限值，之后基本保持恒定，这是由于晶界带上缺陷的数量与分布影响了整个模型的应力分布，使得该内聚力单元在载荷增加的情况下并未完全失效，而是具有一定承载能力。

从图6b中可以发现，对于单元26509来说，当缺陷率≤10%时，随着载荷的增加，单元积分点的应力值大体上逐渐下降;当缺陷率在20% ～50%之间时，随着载荷的增加，单元积分点的应力值相应地降低至某个极限值，之后基本保持恒定。

图6 不同缺陷率下单个内聚力单元的应力值随载荷的变化情况

3.2 缺陷率对损伤值的影响

3.2.1 单个内聚力单元的损伤值随载荷变化规律选取单元号为13159、13733、1295、13162 的内聚力单元，在载荷 U1分别为0.001 mm、0.002 mm、0.003 mm、0.004 mm、0.005 mm的情况下，分别考察这些单元积分点的SDEG值随U1的变化如图7a所示。从图中可以发现，随着载荷的增加，这些单元的SDEG值基本保持恒定。

定义SDEG值在0.9～1之间，内聚力单元失效。在载荷U1分别为0.001 mm、0.002 mm、0.003 mm、0.004 mm、0.005 mm 的情况下，单元 24502、26509的损伤值随载荷的变化曲线分别如图7b、7c所示。从图7b中可以发现，对于单元24502来说，随着载荷的增加，单元积分点的SDEG值均保持在0.9～1之间，意味着该单元已经失效。

图7 不同缺陷率下单个内聚力单元的损伤值随载荷的变化情况

从图7c中可以发现，对于单元26509来说，随着载荷的增加，单元积分点的SDEG值出现波动，之后基本保持恒定。

图8 不同缺陷率下单个内聚力单元的损伤值随增量步的变化情况

3.2.2 单个内聚力单元的损伤值随增量步变化规律

选取24502号单元，在不同的缺陷率下，考察该单元积分点的SDEG值随增量步的变化如图8所示。从图中可以发现，随着增量步的增加，该单元积分点的SDEG值均经过一个阶跃式上升。当缺陷率为0时，经过阶跃式上升之后SDEG值以缓慢的速率上升;当缺陷率在0% ～50%之间时，经过阶跃式上升之后SDEG值基本保持恒定。

4 结论

(1)在不同的缺陷率下，随着载荷的增加，不同内聚力单元的应力值呈现出不同的变化规律，这是由于晶界带上缺陷的数量与分布影响了整个模型的应力分布，各个内聚力单元因而处于不同的受力状态。

(2)在不同的缺陷率下，随着载荷的增大，已经失效的内聚力单元的SDEG值经过波动之后，均保持在0.9～1之间，这是由于内聚力单元的损伤程度随着载荷的增大而不断加深，在达到极限状态之后，SDEG值基本保持恒定。

(3)随着增量步的增加，内聚力单元积分点的SDEG值均经过一个阶跃式上升。当缺陷率为0时，经过阶跃式上升之后SDEG值以缓慢的速率上升;当缺陷率在0% ～50%之间时，经过阶跃式上升之后SDEG值基本保持恒定。

［1］T.Luther，C.Könke.Analysis of crack initiation and propagation in polycrystallinemeso-andmicrostructures ofmetal materials［C］//III European Conference on Computational Mechanics.Lisbon，Portugal: ［s.n.］，2006:1-2.

［2］E.Iesulauro，A.R.Ingraffea，S.Arwade，et al.Simulation of grain boundary decohesion and crack initiation in aluminum microstructuremodels.Fatigue and Fracture Mechanics:33rd Volume，ASTM STP 1417，W.G.Reuter and R.S.Piascik，Eds.，American Society for Testing and Materials，West Conshohocken，PA，715-728，2002.

［3］T.Luther，C.Könke.Investigation of crack growth in polycrystalline mesostructures［C］//17th International Conference on the Application of Computer Science and Mathematics in Architecture and Civil Engineering.Weimar，Germany: ［s.n.］，2006:12-14.

［4］任淮辉，李旭东，李俊琛，等.二维多晶体材料微结构的力学响应计算［J］.武汉科技大学学报，2008，31(4):427-430.

［5］SIL Y，K.Tie.Simulation of polycrystalline aluminum tensile test with crystal plasticity finite element method［J］.Trans.Nonferrous Met.Soc.China，2007(17):1412-l416.

［6］E.Iesulauro.Decohesion of grain boundaries in statistical representations of aluminium polycrystals［D］.Cornell University，2002.

［7］吴治辉.韧性材料损伤断裂过程的数值模拟［D］.西安:西北工业大学，2003.

［8］王艳飞，巩建鸣，蒋文春，等.基于内聚力模型的 AISI4135高强钢氢致滞后断裂数值模拟［J］.金属学报，2011，47(5):594-596.

［9］吴波，魏悦广，谭建松，等.纳米晶Ni晶间断裂的数值模拟［J］. 金属学报，2009，45(9):1077-1081.

［10］Carlos G.D?vila，Pedro P.Camanho，Marcelo F.de Moura.Mixed-mode decohesion elements for analyses of progressive delamination［C］//42nd AIAA/ASME/ASCE/AHS/ASC Structures，Structural Dynamics and Materials Conference.Seattle，Washington: ［s.n.］，2001:2-6.

［11］崔浩，李玉龙，刘元镛，等.基于粘聚区模型的含填充区复合材料接头失效数值模拟［J］.复合材料学报，2010，27(2):162-164.

［12］黄刘刚.内聚力模型的分析及有限元子程序开发［D］.郑州:郑州大学，2010.

［13］ABAQUSVersion6.9 Analysis User’s Manual.2009.

［14］陈国华，陈灵杰.解相关路径的裂纹扩展模拟及其应用［J］. 压力容器，2012，29(10):37-38.

［15］赵宁，欧阳海彬，戴建京，等.内聚力模型在结构胶接强度分析中的应用［J］.现代制造工程，2009(11):128-129.

［16］刘红霞.复合材料分层损伤的数值模拟［D］.西安:西北工业大学，2006.

［17］司良英，邓关宇，吕程，等.基于Voronoi图的晶体塑性有限元多晶几何建模［J］.材料与冶金学报，2009，8(3):193-197.

［18］苏金明，阮沈勇.MATLAB实用教程(第2版)［M］.北京:电子工业出版社，2008.

［19］倪黎，李旭东.应用TransMesh进行微观组织结构的有限元网格划分［J］.兰州理工大学学报，2005，31(2):24-28.

［20］T.Venugopal，K.Srinivasa Rao，K.Prasad Rao.Studies on friction stir welded AA 7075 aluminum alloy［J］.Trans.Indian Inst.Met，2004，57(6):660.

［21］石亦平.ABAQUS有限元分析实例详解［M］.北京:机械工业出版社，2006.

(编辑 李秀敏)

The Application of Cohesive Elements in Numerical Simulations of the Intergranular Fracture in Voronoi Cellular Models

ZHANG Fan，ZHAO Jian-ping
(School of Mechanical and Power Engineering，Nanjing University of Technology，Nanjing 211816，China)

The 2D Voronoi algorithms were used for geometric modeling of crystalline grains.A cohesive interfacemodel was used to simulate the process of grain-boundary separation.Polycrystalline aluminum was chosen as the research object.Under the condition of no defectswhile the uniaxial stretching load was applied，damage and failure of the grain boundaries as well as the stress distribution of the whole model were investigated.Based on this，a comparative study on grain-boundary damage and stress distribution was carried out in the case of defective while uniaxially stretched.It has been found that in the case of uniaxial stretch continuous failure elements appeared on the grain boundary perpendicular to the stretching direction，closest to the loading surface.This led to the formation of the so-called intergranular fracture‘crack’.The simulation results showed that as the uniaxial tensile load increased the damage to the grain boundaries of polycrystalline aluminum would exacerbate，and due to the effect of quantity and distribution of the defects in grain boundaries on stress distribution of the entire model different cohesive elements showed a different variation in stress value with the increase of the load.

polycrystalline aluminum;intergranular fracture;grain-boundary damage

张帆(1987—)，男，南京人，南京工业大学机械与动力工程学院硕士，主要研究方向为失效分析与寿命预测，(E-mail)gzlza@163.com。

TH122;TG111.91

A

1001-2265(2013)11-0026-04

2013-03-20

“十二五”国家科技支撑计划项目(2011BAK06B03和2012BAK13B03)