张芳芳 刘 才

燕山大学国家冷轧板带装备及工艺工程技术研究中心，秦皇岛，066004

0 引言

三维编织复合材料以纤维束空间交织成立体网状为主要结构特征，具有损伤容限大、比模量大、比强度高等优点，在航空航天等领域得到广泛应用。

由于编织结构复杂，材料的损伤模式也比较复杂，因此，通过实验很难观察到复合材料内部的损伤演变过程。而有限元法可以弥补上述不足，因此越来越多的学者采用有限元法模拟其损伤演变过程。通过有限元法研究复合材料损伤时的损伤演化方法包括刚度折减法、连续介质损伤力学方法和基于断裂力学的损伤演化方法。其中，刚度折减法是对材料的刚度进行直接折减，操作过程方便，很容易在有限元中实现，但折减系数的选取依赖经验［1］。连续介质损伤力学方法中的损伤演化方程是基于热力学框架下导出的，损伤演化规律与耗散势相关。Kachanov［2］首次引入了“连续性因子”的概念来描述低应力脆性蠕变损伤。之后很多学者提出了损伤分析模型，Ladeveze等［3］将应力与应变张量分解为正负两部分，弹性模量的退化由损伤参数来描述，演化规律由热动态力控制。基于断裂力学的损伤演化方法认为复合材料在最终破坏时，材料的耗散能与材料的断裂能相等，损伤演化与材料不同裂纹形式的断裂能相关。Camanho等［4］首先提出了一个与材料断裂能相关的损伤演化模型，并引入到内聚力单元中对复合材料的分层过程进行了分析。Lapczyk等［5］提出了正交各向异性的损伤本构模型，预测了弹脆性材料的渐进损伤过程。Fang等［6］基于Murakami－Ohno损伤理论建立了正交各向异性损伤本构模型，预测了三维四向编织复合材料的渐进损伤过程。文献［7－8］基于区域叠合技术预测了复合材料的弹性性能。应用有限元方法预测复合材料性能时，大多采用共节点方法建立复合材料单胞模型，为施加周期性边界条件，要求相对边界面上节点一一对应，这样即使采用四面体进行网格划分也十分困难，区域叠合技术在模型建立上避免了上述共节点等网格划分的困难，但不便于基体中非叠合区域的结果提取与显示。综上所述，刚度折减法中折减系数大小影响结果预测，结合区域叠合技术与断裂能量法建立描述复合材料渐进损伤的模型，提高了建模效率的同时也弥补了上述不足。

本文 基 于 区 域 叠 合 技 术［9－10］，利 用 ANSYS APDL语言建立了参数化编织件网格模型，结合单胞增强相网格提取算法，实现三维编织复合材料单胞模型的参数化建立，提出了含损伤刚度匹配方法，使区域叠合技术中两相模型重合区域的材料刚度与实际增强相材料刚度相匹配。基于Murakami损伤理论建立了正交各向异性损伤本构模型，通过三维Hashin与Mises准则判断增强相和基体的初始损伤，通过等价位移控制损伤变量的演变。利用ANSYS用户子程序接口（Usermat）开发材料模型子程序，应用该模型分别对典型大小编织角三维四向编织复合材料在单向拉伸载荷作用下的渐进损伤过程和拉伸强度进行数值预测。通过基体损伤结果映射方法，实现了基于区域叠合技术所建基体模型渐进损伤演变过程的云图显示。

1 编织复合材料细观有限元模型

计算机图像分析技术发现［11］，编织复合材料中由于纤维束的相互挤压，其横截面形状接近菱形的4个角经倒圆所形成的图形。文献［12］考虑了内部纤维束的真实形态，假设内部纤维束横截面为一个内切椭圆的八边形，如图1所示，建立了三维四向编织复合材料内部细观结构单胞模型，其空间拓扑几何关系如图2所示。本文采用此单胞模型进行分析。

图1 内部纤维束横截面形状

图2 纤维束空间拓扑几何关系

图1和图2中纤维束内切椭圆的长短半轴分别为a和b，h为单胞模型的高度，其中

式中，γ为纤维束内部编织角。

根据单胞模型空间拓扑几何关系，本文采用ANSYS APDL语言应用六面体单元建立三维四向编织预制件网格模型，结合基于Fortran语言编写的单胞增强相网格提取算法，提取出了增强相单胞网格模型，如图3a所示。建立了复合材料单胞整体区域（包括单胞中所有增强相和基体相所占几何空间）网格模型，如图3b所示。将增强相网格模型与整体区域网格模型在空间叠合，组成用于区域叠合有限元技术分析的复合材料单胞网格模型，如图3c所示。

图3 采用区域叠合技术建立单胞网格模型

在区域叠合技术中，均采用等参单元对模型进行离散。增强相网格模型与整体区域网格模型在空间叠合后，通过建立增强相单元节点与整体区域单元节点自由度间的耦合方程，使两模型单元节点的变形相协调。区域叠合技术中，周期边界条件［13］施加在整体区域单元上，由于整体区域模型为规则的长方体，因此很容易满足周期边界条件中对相对边界面上节点一一对应的要求。

2 渐进损伤演变模型

2.1 初始损伤判断准则

纤维束初始损伤判断采用三维Hashin准则［14］。基体采用Mises准则，形式如下：

式中，σ11、σ22、σ33为材料点正应力分量；σ12、σ23、σ31为材料点切应力分量；σm为基体的破坏强度。

2.2 基于材料断裂能的损伤演变模型

为减小局部损伤的网格依赖性，引入单元特征长度［15］建立有限元网格与组分材料断裂能的联系，即假设组分材料不同破坏模式的断裂能量密度为常数，破坏应变随着有限元网格尺寸的变化而改变。在此假设有限单元的特征长度是单元体积的三次立方根，破坏平面的面积是单元特征长度的平方［5-6］。当组分材料局部破坏时，单元的释放能与单元的弹性应变能相等，即

式中，l为有限单元的特征长度；GI、εIf和σIf分别为I型破坏模式的断裂能量密度、等价峰值应变和等价峰值应力。

定义组分材料破坏点的等价位移为

根据上述形式，可以得到不同破坏模式对应的等价位移和等价应力［6］，如表1所示，其中，L、T、Z分别表示纤维束的三个主轴方向，即L表示纤维束轴向，T、Z均表示纤维束横向；t表示拉伸，c表示压缩；＜x＞＝ （x＋｜x｜）／2。

表1 不同破坏模式对应的等价位移和等价应力

当满足初始损伤判断准则即当等价位移超过初始损伤等价位移时，可以用组分材料损伤演变方程控制损伤演变，不同破坏模式下的损伤演变方程为

基于Murakami损伤模型，可利用二阶损伤张量描述增强相和基体损伤，损伤张量为

式中，D、Di、ni分别为损伤张量、损伤张量主值和主方向单位矢量。

损伤演变过程中，有效应力σ＊与名义应力σ的关系为

式中，I为单位矩阵。

采用应变能等效假设，将损伤张量主值引入到材料刚度矩阵中，即

式中，C为未发生损伤时的材料刚度矩阵；C（D）为包含损伤的材料刚度矩阵。

损伤张量主值取该方向拉压损伤变量中的最大值［6］，其表达形式为

基体为各向同性材料，因此其各损伤主值相同。为提高含损伤本构有限元算法的收敛性，引入黏性规则化算法［16］，在当前t＋Δt时刻的损伤张量主值可变换为

将新的损伤张量主值代入式（9），可获得当前损伤状态下的材料刚度矩阵。

3 含损伤刚度匹配方法

区域叠合技术中将基体材料属性赋予整体区域网格模型，为使增强相模型与基体模型重合区域的材料刚度与实际增强相材料刚度相匹配，需要对赋予增强相模型的材料刚度进行刚度匹配处理［10］，本文针对复合材料损伤问题，对刚度匹配方法做了进一步修正，从而使损伤分析过程中，增强相模型与基体模型重合区域的材料刚度与实际增强相材料刚度相匹配。

含损伤刚度匹配方法是通过对增强相积分点材料刚度矩阵的修正实现的，其形式为

式中，CFF（D）和CFM（D）分别为根据增强相材料属性和基体材料属性组建的材料刚度矩阵；CF（D）为赋予该增强相积分点处的材料刚度矩阵。

在损伤分析过程中，对于变形历史过程中未发生过损伤的增强相积分点，根据当前增强相积分点处的应变状态，分别按CFF（D）和CFM（D）计算应力，将计算出的应力分别代入Hashin和Mises准则中进行初始损伤判断，如果不满足初始损伤判断准则，则不对CFF（D）和CFM（D）进行处理；如果满足损伤判断准则，结合式（10）和式（11）计算损伤张量主值，并分别代入式（9）计算CFF（D）和CFM（D），从而获得当前损伤状态下的材料刚度矩阵；对于之前已经发生过损伤的增强相积分点，按式（10）和式（11）对损伤张量主值进行更新。按照式（12）计算更新后的材料刚度矩阵CF（D），将CF（D）返回主程序完成材料损伤本构计算，从而使两相模型重合区域的材料刚度与实际增强相材料刚度相匹配。

4 损伤分析流程

渐进损伤分析过程主要包括有限元平衡方程的应力求解、损伤模式判断和基于材料断裂能的刚度退化，分析流程如图4所示。

图4 渐进损伤分析流程图

在渐进损伤分析过程中，对于初始载荷增量步，采用初始材料刚度求解非线性平衡方程并计算各积分点的应力状态。对于以后的每一载荷增量步，采用上一载荷增量步结束时的材料刚度求解非线性平衡方程并计算各积分点的应力状态。

对于未发生过损伤的积分点，将计算出的应力代入初始损伤判断准则中，如果不满足初始损伤判断准则，则进入下一载荷增量步求解，如果满足初始损伤判断准则，则按照基于材料断裂能的损伤演变模型计算相应的损伤张量主值。对于已发生过损伤的积分点，则根据基于材料断裂能的损伤演变模型对损伤张量主值进行更新。按照更新后的损伤张量主值对积分点的材料刚度进行退化，并进入下一载荷增量步求解，如此循环完成分析。

5 数值分析结果与讨论

采用文献［17］提供的实验数据作为验证算例。组分材料的性能参数如表2所示。

表2 组分材料性能参数

利用上述材料参数预测典型大小编织角三维四向编织复合材料的拉伸强度和渐进损伤演变过程。表3给出了两种编织角试件的拉伸强度、断裂应变预测值与文献［17］中提供的实验数据的对比。由于实验所用试件中可能存在孔隙等细观缺陷，而本文所建模型无法考虑这些可能的缺陷，因此预测值略高于实验值。

表3 试件实验参数与数值预测结果

5.1 基体损伤结果映射法

在后处理时，区域叠合技术所建整体区域模型不便于基体渐进损伤过程的结果显示，因此本文应用基体损伤结果映射法。在该方法中，首先建立与传统方法所建基体具有相同几何空间的模型（此模型没有传统方法中对网格的限制和要求），并采用自由网格划分，将该网格模型与整体区域模型在空间叠合，对于该网格模型中的节点，在全局坐标系下判断出包含该节点的整体区域单元，根据整体区域单元节点的全局坐标，计算出该节点在相应整体区域单元中的自然坐标值，并将计算出的形状函数值代入下式：

式中，Ni为单元节点的形状函数；Di为单元在节点i处的损伤张量。实现根据整体区域单元节点的损伤张量插值出该网格模型节点的损伤张量。

当网格模型的所有节点都通过此方法插值完成后，将网格模型中每个单元的所有节点的损伤张量进行平均，用平均后的损伤张量表示网格模型单元的损伤状态。利用此网格模型单元的损伤结果信息实现基体损伤云图的显示。

5.2 损伤演变过程

基体损伤云图的提取采用上述结果映射法。由于增强相模型具有对称性，为了更好地观察其损伤演变过程，提取其中一个方向纤维束的损伤演变云图进行展示。数值模拟中小编织角增强相模型在纤维轴向拉伸破坏模式下的损伤单元最多，大编织角增强相模型在纤维横向拉伸破坏模式下的损伤单元最多，因此分别提取小编织角增强相在纤维轴向拉伸破坏模式下的损伤演变云图和大编织角增强相在纤维横向拉伸破坏模式下的的损伤演变云图，见图5a和图6a。大小编织角复合材料的基体损伤演变云图分别见图5b和图6b。

图5 小编织角损伤演变云图

从图5a可以看出，小编织角复合材料损伤首先发生在纤维束交错面处，随着拉伸位移的增大，损伤区域沿着纤维束横向和表面逐渐扩展，当应变达到0.0073时，损伤的逐渐累积导致增强相最终失去承载能力。从图5b可以看出，基体损伤首先发生在基体与纤维束交错面处，随着拉伸位移的逐渐增大，损伤区域沿着基体棱边扩展并逐渐贯通融合，随着基体损伤的累积和增强相承载能力的丧失，导致基体最终失效。在整个变形历史过程中，纵向载荷主要由纤维束承担，出现损伤后，纤维束拉伸损伤区域快速扩展，直至最终脆性断裂。因此，小编织角复合材料以纤维拉伸破坏为主要破坏模式，与实验破坏断面的扫描结果［18］一致。

图6 大编织角损伤演变云图

从图6a可以看出，当应变达到0.0052时，大编织角复合材料首先在纤维束表面出现横向损伤，随着拉伸位移的增大，损伤区域沿着纤维束横向扩展。在拉伸过程中，纤维束逐渐向拉伸方向转动，与小编织角不同，纤维束受纵向和横向综合作用，损伤以横向为主。从图6b可以看出，当应变达到0.0062时，基体在与纤维束相交错的棱边发生初始损伤，损伤区域随着拉伸位移的增加沿棱边逐渐扩展，直至失效。因此，大编织角复合材料拉伸强度主要由基体的强度和纤维束的横向强度控制，上述结果与实验破坏断面的扫描结果［18］一致。

6 结论

（1）本文基于区域叠合技术，利用ANSYS APDL语言建立了参数化编织件模型，结合单胞增强相网格提取算法，实现三维编织复合材料单胞模型的参数化建立，提高了单胞模型的建模效率。

（2）结合区域叠合技术与断裂能量法，基于Murakami损伤理论建立了正交各向异性损伤本构模型，应用含损伤刚度匹配方法，实现对典型大小编织角三维编织复合材料渐进损伤过程和拉伸强度的数值预测。

（3）根据基体损伤结果映射法，将整体网格模型的单元损伤状态映射到采用传统方法建立基体几何模型并采用自由网格划分后的网格模型上，实现区域叠合技术所建复合材料中基体渐进损伤演变过程的结果提取与显示。

（4）小编织角复合材料以纤维拉伸破坏为主要破坏模式，材料呈脆性断裂，具有较高的强度。大编织角复合材料中纤维束处于非均匀的受拉状态，拉伸强度主要受基体强度和纤维束的横向强度控制，拉伸强度较低。数值预测结果与实验值均吻合较好。

［1］Blackketter D M，Walrath D E，Hansen A C.Modeling Damage in a Plain Weave Fabric－reinforced Composite Material［J］.Journal of Composites Technology and Research，1993，15（2）：136－142.

［2］Kachanov L M.On the Time to Failure under Creep Conditions Izv［J］.ANSSSR，Otd.Tekhn.Nauk，1958，8：26－31.

［3］Ladeveze P，Lemaitre J.Damage Effective Stress in Quasi Unilateral Conditions［C］／／Proceeding of the 16th International Congress of Theoretical and Applied Mechanics.Lyngby，Denmark，1984.

［4］Camanho P P，Dávila C G.Mixed－mode Decohesion Finite Elements for the Simulation of Delamination in Composite Materials［J］.NASA－Technical Paper，2002，211737（6）：1－37.

［5］Lapczyk I，Hurtado J A.Progressive Damage Modeling in Fiber－reinforced Materials［J］.Composites Part A：Applied Science and Manufacturing，2007，38（11）：2333－2341.

［6］Fang G D，Liang J，Wang B L.Progressive Damage and Nonlinear Analysis of 3DFour－directional Braided Composites under Unidirectional Tension［J］.Composite Structures，2009，89（1）：126－133.

［7］Jiang W G.Implementation of Domain Superposition Technique for the Nonlinear Analysis of Composite Materials［J］.Journal of Composite Materials，2013，47（2）：243－249.

［8］张芳芳，姜文光，于春蕾，等.基于区域叠合有限元技术预测三维编织复合材料弹性性能［J］.燕山大学学报，2012，36（3）：219－223.Zhang Fangfang，Jiang Wenguang，Yu Chunlei，et al.Prediction of Elastic Properties of 3DBraided Composite Using Domain Superposition Technique［J］.Journal of Yanshan University，2012，36（3）：219－223.

［9］Jiang W G.A Computer and a Method of Modelling a Woven Composite Material：International Patent，PCT／GB2008／000713［P］.2008－10－16.

［10］Jiang W G，Hallett S R，Wisnom M R.Development of Domain Superposition Technique for Woven Composites［M］／／Camanho P P，Davila C G，Pinho S T，et al.Mechanical Response of Composites.Berlin：Springer，2008：281－291.

［11］李嘉禄，刘谦.三维编织复合材料中纤维束横截面形状的研究［J］.复合材料学报，2001，18（2）：9－13.Li Jialu，Liu Qian.Study on Fiber Tows Crosssection in 3－D Braided Composites［J］.Acta Materiae Compositae Sinica，2001，18（2）：9－13.

［12］徐焜，许希武.四步法三维矩形编织复合材料的细观结构模型［J］.复合材料学报，2006，23（5）：154－160.Xu Kun，Xu Xiwu.On the Microstructure Model of Four－step 3DRectangular Braided Composites［J］.Acta Material Compositae Sinica，2006，23（5）：154－160.

［13］Jiang W G，Yan L J.Implementation of Stress Loading Repetitive unit Cell Finite Element Model［C］／／Proceeding of the 3rd International Conference on Heterogeneous Material Mechanics.Shanghai，2011：816－819.

［14］Hashin Z.Failure Criteria for Unidirectional Fiber Composite［J］.Journal of Applied Mechanics，1980，47：329－334.

［15］Bažant Z P，Oh B H.Crack Band Theory for Fracture of Concrete［J］.Materials and Structures，1983，16（3）：155－177.

［16］Duvaut G，Lions J L.Inequalities in Mechanics and Physics［M］.Berlin：Springer，1976.

［17］修英姝.四步法三维编织复合材料力学性能的有限元分析［D］.天津：天津工业大学，2001.

［18］卢子兴，冯志海，寇长河，等.编织复合材料拉伸力学性能的研究［J］.复合材料学报，1999，16（3）：129－134.Lu Zixing，Feng Zhihai，Kou Changhe，et al.Studies on Tensile Properties of Braided Structural Composite Materials［J］.Acta Materiae Compositae Sinica，1999，16（3）：129－134.