CFRP三维铣削仿真模型的建立及层间损伤分析

秦旭达 唐心凯 葛恩德 李士鹏 朱圣富

（1 天津大学机械工程学院，天津 300072）

（2 上海飞机制造有限公司航空制造技术研究所，上海 201324）

0 引言

碳纤维增强树脂基复合材料（CFRP）具有比强度和比模量高、抗腐蚀性和耐疲劳性好等诸多优点，在航空航天领域应用广泛［1］。虽然大多数复合材料采用近净成形加工，但通常为了满足零件尺寸及装配要求等，需要对复合材料零件进行铣削加工［2］。但CFRP 是一种典型的难加工材料，在加工过程中易产生层间损伤，严重影响加工质量以及工件的使用寿命［3］。因此，研究CFRP 层间损伤产生机理进而减少加工损伤是当前亟需解决的问题。

由于复合材料试件和刀具价格昂贵以及层间损伤不易检测，仅通过实验方法对加工损伤产生机理进行研究需要很高的实验成本。因此有限元仿真与实验相结合的方式逐渐成为一种较优的分析方法。针对复合材料切削的有限元仿真，大部分国内外学者［4-7］关注二维切削仿真，但二维切削模型无法分析面外失效以及刀具的某些具体角度（螺旋角等）对切削过程的影响，三维切削仿真结果更接近实际，已逐渐成为CFRP 仿真建模的关注重点。CFRP 三维建模需要自定义CFRP 的力学模型，相应理论还未成熟，部分学者已经开始对复合材料三维切削有限元仿真展开研究，例如，A. Vaibhav 等人［8］建立了复合材料三维钻削有限元模型，用Cohesive单元模拟复合材料中的层间分层，分析了钻孔工艺参数对分层损伤的影响。贾振元等人［9］开发了CFRP/Ti 叠层材料三维斜角切削模型，探究了切削不同纤维方向与不同铺层顺序时面下损伤的变化。宿友亮［10］建立了CFRP的三维正交切削有限元模型，利用该模型对微切削过程进行了仿真，对纤维和树脂的断裂演化过程进行了表征。A.Abena 等人［11］利用光滑粒子流体动力学与有限元相结合的方法建立了CFRP 细观三维正交切削模型，仿真精度较好。齐振超［12］建立了三维钻削CFRP 有限元模型，直观表现了分层产生过程和应力场分布规律。综上可知，对于复合材料三维有限元建模的层间损伤研究，三维正交切削模型与钻削模型占绝大多数，且模型中大多未考虑CFRP 的刚度退化部分。以往的研究中还没有铣削加工中纤维方向角对层间损伤的影响规律研究。

本文基于ABAQUS 有限元仿真技术，建立可预测层间损伤和刚度退化的CFRP 三维铣削有限元模型，分析了切削力、层间应力、层间损伤随纤维方向角的变化规律，对于铣削刀具和工艺优化具有一定的技术参考价值。

1 CFRP三维铣削有限元模型

1.1 CFRP层合板损伤演化理论

建立三维铣削模型为宏观模型，单层板设置为等效均质的模型，不区分纤维与基体［13-14］。通过自定义子程序（VUMAT）的开发，构建出三维线弹性本构、三维失效准则及损伤演化在内的复合材料模型，复合材料板层与层之间通过Cohesive单元连接，以模拟层间损伤。

1.1.1 层内损伤演化理论

CFRP 损伤起始判据采用Hashin 失效准则［15-16］，Hashin准则包含四种不同的失效形式。

纤维拉伸失效（σ11≥0）：

纤维压缩失效（σ11＜0）：

垂直于纤维方向拉伸失效（σ22+σ33≥0）：

垂直于纤维方向压缩失效（σ22+σ33＜0）：

式中，Fft、Ffc、Fmt、Fmc分别是四种失效模式的损伤变量；σ11、σ22、σ33、τ12、τ13、τ23是单元积分点处应力张量的分量分别是纤维方向的拉伸、压缩破坏应力，垂直于纤维方向的拉伸、压缩破坏应力；S12、S13、S23分别是1-2 平面、1-3 平面、2-3 平面的剪切破坏应力。

一旦积分点达到损伤起始条件，该点就会进行刚度退化，刚度退化通过刚度退化系数进行控制，而刚度退化系数的演化规律尚未有统一的形式，其中指数损伤演化规律已经被证明是有效的，具有良好的收敛性［17］。由于Hashin 准则包含四种失效条件，与之对应有四种损伤因子di。

纤维拉伸失效（σ11≥0）：

纤维压缩失效（σ11＜0）：

垂直于纤维方向拉伸失效（σ22+σ33≥0）：

垂直于纤维方向压缩失效（σ22+σ33＜0）：

式中，LC为单元的特征长度，Gft、Gfc、Gmt、Gmc为纤维拉伸、纤维压缩、基体拉伸、基体压缩的断裂能为单元积分点处应变张量的分量。

损伤因子di控制着刚度退化，如公式（9）所示［17］，即当满足起始失效准则后，损伤产生并累积，损伤因子从0 开始增加，当达到1 时，刚度退化为最小值。当dft达到1 时，即沿纤维方向拉伸完全失效，单元被删除。

1.1.2 层间损伤演化理论

层间通过Cohesive 单元连接并传递应力，Cohesive 单元应力应变相应如图1 所示，在切削载荷作用下，Cohesive 单元发生线弹性变形，达到临界等效应力后产生损伤发生刚度退化，最终单元删除，层间开裂。

图1 Cohesive单元应力应变关系Fig.1 Stress-strain of cohesive elements

采用二次名义应力准则作为层间单元损伤起始判据，此准则用公式（10）表示。当各方向的名义应力比的平方和达到1时，损伤产生。

式中，tn、ts、tt为界面法向和剪切方向的瞬时应力、为界面法向和剪切方向的最大名义应力。

关于Cohesive单元的损伤演化，基于能量的损伤演化规律得到广泛应用，采用基于断裂能的Powerlow判据确定断裂能在混合模式下的关系如下：

式中，Gn、Gs、Gt为界面法向及剪切方向的瞬时断裂能为法向及剪切方向破坏所需的临界断裂能，指数常数β为1.6［18］。

仿真中使用的Cohesive 单元的性能参数见表1［18］。

表1 Cohesive单元性能参数1）Tab.1 Property parameters for Cohesive element

1.2 有限元模型设置

三维铣削CFRP 有限元模型主要由三大部分组成：刀具、层合板、层间部分，由于尺寸及层数的增加会导致计算效率急剧降低，而层间应力及损伤主要集中在切削区，故将模型简化为如图2 所示，工件的尺寸为5 mm×3 mm×2 mm，分为4 层，每层厚0.5 mm；层与层之间设置Cohesive 单元。远离切削区的一侧1 mm 宽的工件区域设置为完全约束。由于刀具的刚度远大于工件，为节省计算时间，将刀具定义为刚体。刀具的具体结构参数见表2。通过控制刀具上参考点的移动和转动以模拟铣削运动过程，铣削的工艺参数见表3。

图2 铣削CFRP有限元模型Fig.2 Finite element model for milling CFRP

表2 刀具结构参数Tab.2 Tool geometry parameters

表3 切削工艺参数Tab.3 Cutting process parameters

网格的设置在涉及高变形问题的模拟中十分重要，CFRP 单层板的单元网格类型设置为C3D8R（8节点三维线性实体，减缩积分沙漏控制单元），层间Cohesive 单元网格类型设置为COH3D8（8 节点三维粘聚力单元），刀具网格类型为C3D4（4 节点线性实体）。网格尺寸大小的设置要考虑计算精度及时间效率的平衡，经过仿真测试，宏观切削模型网格在0.1 mm以内精度较好，大于0.1 mm仿真精度急剧变差，故单层板单元网格大小设置为0.1 mm×0.1 mm×0.1 mm；层间单元厚度为0.01 mm，单元网格大小为0.1 mm×0.1 mm×0.01 mm；CFRP 层合板的单元数量为34 500 个。刀具为刚体，不参与计算，网格大小设置为0.3 mm，单元数量为76 000个。

刀具与CFRP 的接触设置为面与节点（刀具的表面与工件的所有节点）的接触，摩擦因数为0.5［19］。与此同时，为避免切削过程中工件内部单元相互嵌入，对整体模型采用自接触。

模型中时间总长为0.29 s，整个模型需在16 GB运行内存的计算机上运行28 h。

2 实验设计

为了研究铣削过程中加工损伤类型及分布，采用德玛吉五轴加工中心DMU80T 对四种典型的纤维方向角的CFRP 层合板进行铣边加工。实验中所用的工件材料为恒神股份有限公司的CFRP 单向层合板，碳纤维型号为T700，质量密度为1.21 g/cm3，材料的其他性能参数见表4。

表4 CFRP材料性能参数Tab.4 Property parameters for CFRP

工件尺寸大小为250 mm×120 mm×2.5 mm。如图3 所示，CFRP 层合板通过两块压板固定以降低切削过程中振动，刀具采用Φ10 mm 硬质合金铣刀，铣刀结构参数见表2。铣刀通过逆铣的方式沿着工件的长边缘不断向前进给，工件厚度为2.5 mm，因而切削深度ap为2.5 mm，其他工艺参数见表3。采用Kistler9257B 三向测力仪对切削力进行测量，采用VHX-1500 超景深电子显微镜对加工表面形貌进行观测。

图3 铣削实验设置Fig.3 Milling experimental device

3 仿真结果分析

3.1 切削力

图4 为铣削CFRP 单向层合板的三向切削力随纤维方向角的变化规律。由于仿真中模型厚度的增加会导致计算时间的急剧增加，因此实验与仿真均取单位厚度的切削力（N/mm）进行对比验证，其中实验与仿真的铣削力分别取稳定切削区域的切削力峰值的平均值。通过对比切削不同纤维方向角的CFRP 产生的切削力可得，y向切削力（进给切削力）与z向切削力仿真与实验吻合精度较高，误差控制在20%以下，x向切削力仿真精度稍差。由于切削仿真是在刀具完全锋利并没有磨损的理想条件下进行的，因此仿真得到的切削力大部分相比实验切削力偏小。

图4 CFRP单向板切削力随纤维方向角的变化Fig.4 Variation of cutting force for CFRP unidirectional laminates with fiber orientation angle

由图4 可得，x向切削力和y向（进给方向）切削力均随着纤维方向角的增大先增加后减小，但纤维方向角为45°时x向切削力最大，纤维方向角为90°时y向切削力最大。z向切削力随着纤维方向角的增加呈增大趋势，纤维方向角为90°和135°时z向切削力较大。如表5 所示，铣削过程中纤维切削角（纤维切削角为切削速度与纤维方向所成的夹角）是动态变化的，由于径向切宽为1 mm，刀具直径为10 mm，切削纤维方向角为0°的CFRP 时，纤维切削角为大钝角（135°＜β≤180°），此过程破坏形式主要为水平方向的剪切破坏，故所需的三向切削力均较小；当纤维方向角为45°时，纤维切削角为小锐角（0°＜β≤45°），破坏形式主要为沿纤维方向的剪切破坏，x向与y向分力大致相等。而当切削纤维方向角为90°的CFRP时，纤维切削角为大锐角（45°＜β≤90°），纤维束大部分直接被切断，故所需的y向（进给方向）切削力最大；纤维方向角为135°时，纤维切削角为小钝角（90°＜β≤135°），破坏形式主要为弯曲断裂，所需的切削力小于切削纤维方向角为90°的CFRP 时所需的切削力。

表5 不同纤维方向角下纤维切削角的变化Tab.5 Variation of cutting angle for CFRP unidirectional laminates with different fiber orientation angle

3.2 层间应力

在刀具切削CFRP 过程中，层内与层间均产生应力，如图5所示，但由于层间单元的破坏应力低，层间应力先达到极值产生损伤，导致层间产生开裂，形成撕裂、分层等加工损伤。层间应力与加工损伤有密切的关系，但层间应力通过实验很难测得，而利用有限元仿真中的Cohesive 单元可提取CFRP 在铣削加工过程中层间应力大小。

加工过程中层间应力随着刀具的移动不断变化，故本文选取与切削区有一定距离（保证整个切削过程未被切除）的4个单元（16个节点），如图6所示，提取节点上的层间应力随时间的变化曲线，取其16组数据的平均值作为层间应力大小的参考标准，在相同位置点的层间应力的平均值越大，说明层间损伤越严重。

图5 切削过程中CFRP上的应力分布Fig.5 Stress distribution on CFRP during cutting

图6 切削过程中CFRP上的应力分布Fig.6 Extraction point of interlaminar stress

经计算处理，图7 为有限元模型中铣刀铣削CFRP 产生的层间应力随纤维方向角的变化曲线，采用Mises等效应力表达层间应力大小［20］。

图7 不同纤维方向角下层间应力变化Fig.7 Variation of interlaminar stress with different fiber orientation angle

由图7看出，层间应力随着纤维方向角的增加而增大。纤维方向角为0°层间应力最小，纤维方向角为135°时层间应力最大。纤维方向角为0°时，轴向力最小，相应法向正应力小，且切削过程容易从界面处断裂，层间单元所受的剪切应力小，故层间Mises总应力最小。纤维方向角为45°时，轴向力与剪切力均增大，层间Mises 总应力增加。纤维方向角为90°、135°时轴向力差距不大，即法向正应力相近，纤维方向角为90°时，切削刃能及时切断大部分纤维，但纤维方向角为135°时，纤维在刀刃挤压作用下先发生弯曲，这使得层间单元所受的剪切应力最大，故层间Mises总应力最大。

3.3 层间损伤

CFRP由于层间结合强度低，在切削过程中易产生层间损伤，图8为铣削CFRP层合板层间损伤模拟，由图可定性的观察出切削过程是否产生了层间损伤，但定量的去表征层间损伤的大小对于对比分析十分重要。

图8 铣削过程分层损伤模拟Fig.8 Simulation of delamination in milling process

切削过程中，层间单元在切削力的作用下应力达到极值后会产生刚度退化，当刚度退化到0.99后，材料已基本失去承载能力，可作为层间损伤区域［8］。故本文将刚度退化系数达到0.99的区域作为层间损伤区域。针对铣削加工的层间损伤提出二维层间损伤因子。二维层间损伤因子是层间损伤区域的面积与被切除材料面积的比值，如图9所示。

图9 二维层间损伤因子表征示意图Fig.9 Schematic diagram of 2D interlaminar damage factor characterization

二维层间损伤因子（FD）的表达式如下所示：

式中，AD表示层间损伤区域的面积，AT为被切除材料的面积，AT=ae×L，其中ae为侧铣加工时的切宽，L为工件被加工一侧的长度。

在右旋铣刀铣削纤维方向角为45°的CFRP有限元结果中，提取上、下二层Cohesive单元的刚度退化如图10所示，层间损伤区域的面积AD为红色虚线内的白色区域及红色区域面积之和。利用Matlab软件对切削稳定区域的层间损伤图片进行二值化处理，计算出层间损伤面积AD占图片面积的比值以求得分层损伤面积AD的实际大小，进而求得二维层间损伤因子FD。对四种纤维方向角的有限元模型中切削稳定区域后处理得到各自的二维层间损伤因子FD见图11。

图10 不同位置的层间损伤Fig.10 Interlaminar damage in different position

图11 二维层间损伤因子Fig.11 2D interlaminar damage factor

图12 纤维方向角对层间损伤的影响Fig.12 Effect of fiber orientation angle on interlaminar damage

由图10 与图11 可以看出，当刀具切削四种纤维方向角的CFRP 时，均呈现上层层间单元损伤大于下层层间单元损伤的规律。这是因为刀具为右螺旋铣刀，切削时切削刃产生一个向上的切削分力CFRP 上表层的上部没有约束，故层间损伤最大，层间损伤后产生毛刺、撕裂等加工损伤。CFRP 下表面在切削应力挤压下也有部分失效，但上侧有材料的约束，故下表面的层间损伤较小。

由于分层、撕裂、毛刺等加工损伤集中在工件的上表面，故针对工件的上层层间单元进行分析，由图11 和图12 得层间损伤区域随纤维方向角的增加而增大。纤维方向角为0°时，层间应力最小，几乎没有层间损伤，加工表面质量最好；纤维方向角为45°时，受沿着纤维方向剪切失效的影响，纤维与基体界面开裂，但由于纤维方向角的影响，未切除的纤维束在刀具进给方向的反方向，纤维束不易被掀起，CFRP内部产生较小的层间损伤，已切削的部分则因纤维与基体界面开裂形成丝状毛刺。纤维方向角为90°时，大部分纤维束被直接切断，未切断的纤维束由于纤维束与刀具进给方向成90°，受剪切应力的作用产生层间损伤，纤维束被掀起，形成撕裂以及片状毛刺，层间损伤较大；纤维方向角为135°时，未切除的纤维束在刀具进给方向的同方向，大部分纤维束不能被切断，受弯矩应力的作用CFRP 内部沿着纤维方向产生层间损伤，纤维束基本都被掀起形成撕裂和片状毛刺。由此可得，纤维方向角为135°时，层间损伤最大。

4 结论

通过建立的包含CFRP 层间、层内单元损伤演化的三维铣削CFRP 有限元模型，探究了切削力、层间应力与层间损伤随纤维方向角的变化关系，同时针对层间损伤基于刚度退化理论提出了二维层间损伤因子，实现了动态切削过程中的层间损伤仿真预测。已验证的三维铣削CFRP 有限元模型可用于对于刀具结构以及加工工艺参数的优化。通过仿真结果分析，得出以下结论：

（1）通过对比宏观铣削模型与实验测得的切削力可知，x向和y向切削力均随着纤维方向角的增大先增加后减小，纤维方向角为90°和135°时z向切削力较大。仿真与实验结果变化规律一致。

（2）层间应力随着纤维方向角的增加而增大，纤维方向角为135°时，层间应力最大。

（3）层间损伤与层间应力随纤维方向角的变化趋势一致，纤维方向角为0°时，层间损伤最小；纤维方向角为135°时，层间损伤最严重。