随机多孔碳纤维纸的非线性面外压缩本构模型

石姗姗，吕航宇，吕超雨，苏智博，孙 直

(1. 大连交通大学机车车辆工程学院，辽宁，大连 116028；2. 大连理工大学工业装备结构分析国家重点实验室，辽宁，大连 116023)

碳纤维纸在先进能源以及先进环境材料领域中有广泛应用前景，如功能材料、电化学材料和新能源材料等[1−6]，其主要服役状态为面外压缩。大量的实验和分析研究[7−12]表明，碳纤维纸的面外压缩变形对其力学、热学、电学性能以及流体性能具有显著影响[7−15]。碳纤维纸不同于分布规则的碳纤维编织结构[16−18]，具有随机散乱的结构特征，其微观结构由杂乱的短切碳纤维形成低密度多孔结构，空间几何分布类似于无纺布。这种散乱的多孔结构与现有模型有很大差别，其性能随微孔半径大小而改变[19]，导致碳纤维纸的应力-应变关系不同于常见的均匀材料或连续材料，并呈现出明显的非线性特征。

目前，针对碳纤维纸的研究主要集中于化学处理和微纳观改性领域[20−24]，相关的力学研究较少。对于碳纤维纸的力学性能与耦合性能，多采用实验测量，并根据实验结果进行初步的数值分析。高源等[25]使用随机重构技术建立三维微孔隙结构模型，通过改变孔隙率、纤维直径和燃料电池的气体扩散层厚度等参数观测碳纤维纸的渗透率的变化趋势。GARCÍA-SALABERRI 等[26]充分考虑了碳纤维纸的非线性特点，通过大量的实验数据，拟合了碳纤维纸的非线性应力-应变关系。该模型高度依赖于实验数据的数量和有效性，不利于实际使用，尤其是碳纤维纸型号改变的时候，其应力-应变关系会随着孔隙率、有机质含量、单胞尺寸等参数剧烈变化。BAHRAMI 等[6]基于梁理论，预测压缩载荷作用下碳纤维纸变形的力学解析模型。BAHRAMI 模型考虑了碳纤维纸的微观纤维分布与纤维弯曲变形，但该模型是线性的，只能在很小的范围内与实验数据吻合。ZHANG 等[27]建立了碳纤维纸的随机微观结构，研究发现接触对和孔隙率是决定压缩非线性的2 个关键因素。实际上，碳纤维纸在压缩载荷下，内部孔隙会不断减小，导致结构孔隙率下降，实体率上升，微观结构逐渐密实，从而引发碳纤维纸的弹性模量随压缩应变的增加而快速增加。针对碳纤维纸的非线性应力-应变关系，只有建立一种可以反映其特征的非线性本构模型，才能准确预测碳纤维纸的性能，进而更准确地分析碳纤维纸的力-电性能、力-热性能、力-渗透性能等多场耦合性能。然而，目前尚缺少有效的非线性碳纤维纸本构模型。

为建立有效的非线性本构模型，本文针对碳纤维纸微结构的随机多孔特点，提出了两种本构模型，分别是全面考虑微结构孔隙率、单胞尺寸变化及微结构弯曲与接触变形的对数型模型，以及忽略赫兹接触变形的幂函数型简化模型。首先选取典型六面体微结构作为碳纤维纸的代表性单胞结构，将碳纤维纸的压缩变形分解为接触变形和弯曲变形两部分。分别使用赫兹接触理论和弯曲梁理论分析典型纤维的接触变形和弯曲变形，并考虑压缩变形过程中，结构的孔隙率变化，及其对纤维直径、单胞长度、微结构单元面积等相关几何参数的影响，得到碳纤维纸的非线性本构关系。基于三种商品化的碳纤维纸的应力-应变实验数据，分析对比本文提出的两种本构模型的适用范围。

1 非线性面外压缩模型

1.1 碳纤维纸的代表性六面体单胞

碳纤维纸的多孔结构由杂乱且低密度的短切碳纤维组成，如图1 所示。观察俯视视角和横截面视角的扫描电镜图片，可以发现碳纤维在平面内的分布表现为随机长度和随机角度的四边形(图1(a))；在厚度方向的分布表现为较为规则的堆叠(图1(b))。根据碳纤维纸中碳纤维的实际分布规律，本文选用六面体作为碳纤维纸的代表性单胞，如图2 中虚线所示。六面体单胞四周为4 个矩形，高度是碳纤维直径的4 倍；上、下表面为随机四边形，边长、单胞面积等参数的统计值详见表1。每个代表性六面体单胞包含2 根横向纤维和2 根纵向纤维(第1 层含1 根完整的横向碳纤维，第3 层含有2 个半根的横向碳纤维。同理，第2 层、4 层各含有1 根纵向碳纤维)。

图1 碳纤维纸扫描电镜图片Fig. 1 SEM images of commercial carbon fiber papers[8]

图2 碳纤维纸的代表性六面体单胞示意图Fig. 2 Schematics of the proposed typical hexahedral unit cell model for the carbon fiber papers

1.2 接触变形

碳纤维纸最重要的力学行为是厚度方向上的压缩变形。当六面体单胞承受压缩应力时，水平方向纤维与相邻的垂直纤维发生紧密接触而产生非线性接触变形，如图3 所示。它是碳纤维纸的主要变形方式之一，也是碳纤维纸发生非线性力学响应的重要原因。

图3 碳纤维纸接触变形模式示意图Fig. 3 Schematics of contact deformation with use of Hertz's contact model

表1 给出了几种商用碳纤维纸的力学性能和几何参数。根据这些数据，可以发现碳纤维纸六面体单胞的平均长度(119 µm～136 µm)远大于平均纤维直径(6.95 µm)。因此，六面体单胞内纤维可以近似认为无限长圆柱体，碳纤维间的接触为无限长圆柱的正交(垂直)赫兹接触。

表1 商用碳纤维纸的力学性能和几何参数Table 1 Mechanical properties and geometric parameters of commercial carbon fiber papers

根据正交赫兹接触[30]圆柱的载荷-位移关系，得到：

式中：P和u分别为单胞所承受的压缩载荷和压缩位移；E、D和ν分别代表碳纤维的弹性模量、单根纤维的直径和碳纤维单丝的泊松比。由式(1)可以分别得到单胞平均压缩应力和赫兹接触变形引起的平均压缩应变：

1.3 弯曲变形

如图4 所示，碳纤维纸代表性单胞在受到压缩载荷时，具有一定长度的纤维与纤维接触点会构成三点弯曲。因此，单胞内除了会产生接触变形外，还会产生弯曲变形。弯曲变形是碳纤维纸的另一种主要变形方式。

图4 碳纤维纸弯曲变形模式示意图Fig. 4 Schematics of bending deformation with use of beam bending model

对碳纤维纸结构的弯曲变形，BAHRAMI 等[6]利用欧拉-伯努利梁理论成功地建立了碳纤维纸六面体单胞的弹性本构模型。由于弯曲纤维所受弯矩和变形与纤维间的夹角无关，欧拉-伯努利梁理论同样适用于六面体单胞中的交叉纤维。因此，压缩载荷下六面体单胞的纤维弯曲变形引起的压缩应变[6]可以表示为：

1.4 考虑赫兹接触影响的对数型模型

碳纤维纸代表性单胞在受到压缩载荷时，主要发生2 种变形，即接触变形和弯曲变形。其中，接触变形主要发生在纤维之间的接触点，如图3 所示；弯曲变形主要体现为纤维轴线在厚度方向上发生挠度，如图4 所示。在模型中，两种变形同时发生，单胞的压缩位移等于两种变形位移之和。根据式(5)和式(8)，可以得到 dε 和 dσ之间的微分关系如下：

由式(10)计算得到的碳纤维占比均大于碳纤维纸的实测实体率，这说明碳纤维纸内含有大量空隙。考虑纤维分布的随机性，假设空隙均匀分布，可以得到单胞的充盈率 η。若 η=1，则说明碳纤维纸被六面体单胞充满；若 η＜0，则碳纤维纸内含有 η的六面体单胞与1−η的空隙。单胞的充盈率 η可由碳纤维纸的实测孔隙率与单胞理论孔隙率之比得到：

随着压缩位移与应变的增加，碳纤维纸会变得越来越密实(碳纤维的体积不变，碳纤维纸的体积减小)，即碳纤维纸的单胞充盈率会随着压缩应变增加而增加。忽略泊松效应对碳纤维纸单胞几何尺寸的影响，单胞名义充盈率 η∗可表示为：

式(12)显示，在压缩载荷下，随着名义应变( ε＜0)的减少，单胞名义充盈率η*可能会大于1，这表示厚度方向上原本未接触的碳纤维发生了更多的接触，导致单胞尺寸变小、数量增加。本文近似地认为其对弹性模量的影响是线性的。

将式(12)代入式(9)，得：

通过对式(13)进行两边积分，可得到碳纤维纸的对数型本构关系：

1.5 同时考虑孔隙率和单胞尺寸的幂函数型模型

若进一步考虑压缩变形时，结构堆叠将对单胞尺寸产生影响。随着名义应变 ε的变化，例如，当名义应变达到−0.5 时，单胞尺寸缩小为l· (1−0.5)，则单胞长度l可写为：

式中，l0为压缩变形发生前单胞长度。

由于接触等效模量Ec和弯曲等效模量Eb与单胞尺寸的关系不一样，若同时考虑两种等效模量随单胞尺寸的变化，微分方程将无法分离变量，会给求解带来不便。而且随着压缩变形持续增加，接触变形的影响会变小，压缩变形导致的单胞尺寸变化对碳纤维纸应力-应变关系的作用变得明显。因此，本节只分析弯曲等效模量随单胞尺寸的变化，并建立相应的应力-应变关系。

对弯曲等效模量，根据BAHRAMI 等[6]模型，可得：

同时考虑孔隙率的变化，由式(12)、式(15)和式(16)，可以得到：

即：

通过对式(18)进行两边积分，可得到碳纤维纸的幂函数型本构关系：

2 结果与讨论

本节中，采用表1 给出的SGL24AA、SGL25AA和SGL10BA 三种商用碳纤维纸的力学性能和几何参数，将1.4 节和1.5 节中提出的全面考虑微结构孔隙率、单胞尺寸变化及微结构弯曲与接触变形的对数型模型和忽略赫兹接触变形的幂函数型简化模型，与这三种商用碳纤维纸实验所测应力-应变关系曲线进行对比分析。由于SGL10BA 这类碳纤维纸中含有聚四氟乙烯或热塑性树脂等有机物，这将对其力学本构关系产生影响，因此本节将分别讨论有、无额外有机物的碳纤维纸。

2.1 无有机物碳纤维纸的本构模型验证

如图5 所示，三角曲线分别代表了SGL 24AA 和 SGL 25AA 商业碳纤维纸应力-应变曲线的实验测试值；圆点曲线代表对数型本构模型；方形曲线代表幂函数型本构模型；直线代表线性本构模型[6]。可以看出，针对SGL 24AA 和 SGL 25AA 这两类碳纤维纸，幂函数型本构模型与实验测量值呈现出高度吻合，幂函数型模型具有良好的非线性特征。对数型本构模型与实验数据具有相似的模量强化趋势，线性本构模型与实验结果相差最远。

图5 对数型本构模型、幂函数型本构模型与线性本构模型和实验应力-应变曲线的对比图Fig. 5 Comparison of experimental stress-strain curves with the logarithmic-type model，the power-function-type model and the linear bending model

考虑到碳纤维纸最优夹紧力为0.77 MPa[31]，为了有一个更全的覆盖面，本文将讨论的压缩应力区间设定在0.4 MPa～1.2 MPa。如图6 所示，对于SGL 24AA 碳纤维纸，在其压缩应力为0.4 MPa～1.2 MPa 时，线性本构模型分析得到的应变值与实验测得的应变平均值绝对误差值区间为41.8%～76.6%，对数型本构模型分析得到的应变值与实验测得的应变平均值绝对误差值区间为2.0%～41.8%，而幂函数型本构模型分析得到的应变值与实验测得的应变平均值绝对误差值区间为0.1%～35.0%。对于SGL 25AA 碳纤维纸，对数型本构模型分析得到的应变值与实验测得的应变平均值绝对误差值区间降为0.1%～37.5%；幂函数型本构模型与实验测得的应变平均值吻合度最高，最大绝对误差值为1.3%～33.1%；与SGL 24AA 碳纤维纸相比，线性本构模型分析得到的应变值与实验测得的应变平均值之间误差更大，绝对误差值区间为41.5%～79.8%。

图6 对数型本构模型、幂函数型本构模型与线性本构模型和实验应力-应变曲线误差绝对值对比图Fig. 6 Comparison of absolute error values with experimental stress-strain curves between the logarithmic-type model，the power-function-type model and the linear bending model

对比SGL24AA 和SGL25AA 碳纤维纸的应力-应变实验数据与线性、对数型本构及幂函数性本构预测结果，可以发现：首先，本文提出的两种本构模型的相对误差均远小于线性本构模型[6]，其主要原因为线性本构未考虑压缩过程中纤维因接触点增加而发生的应变强化。其次，当碳纤维纸孔隙率增加时(SGL25AA 的孔隙率大于SGL24AA)，本文提出的两种非线性本构模型的相对误差也会增加。这可能是因为孔隙率较大时，微观结构随机性更大，造成预测误差的增加。最后，对比本文提出的两种本构模型，可以发现在较高应变下，考虑接触变形的对数型本构模型和幂函数本构模型模拟结果相差较大，这主要是碳纤维纸孔隙率和弹性模量随应变改变造成的。在压缩应变较低时，弯曲变形明显大于接触变形，忽略接触变形的幂函数型简化模型具有形式更简单的优点，且在压缩应变较低时与对数型本构模型相差不大。随着压缩应变的增加，碳纤维纸的孔隙率减少、弹性模量成倍增加。此时，忽略接触变形与接触刚度的影响随着弹性模量的增加而增加，进而导致对数型本构模型和幂函数本构模型模拟结果相差较大。总体而言，幂函数型模型的参数较少，在准确性和便利性上更适合工程应用。

2.2 含有有机物的碳纤维纸本构模型验证

SGL10BA 等碳纤维纸中含有较大比例的聚四氟乙烯或其他热塑性树脂，这将对其力学本构关系产生影响，本节将讨论本文提出的本构模型对这类碳纤维纸的适用性。考虑到聚四氟乙烯或热塑性树脂的影响，需要在现有的对数型本构模型和幂函数型本构模型中引入修正因子 γ，其中 γ代表碳纤维纸中碳纤维含量比值。由于有机物的弹性模量均小于碳纤维模量，因此可认为有机物对碳纤维纸模量贡献较小，碳纤维纸的弹性模量主要由其碳纤维在实体材料的比例 γ决定，式(14)和式(19)可以分别修正为：

图7(a)为SGL10BA 碳纤维纸的应力-应变曲线。其中黑色三角为实验测试值[29]；圆点曲线代表对数型本构模型；方形曲线代表幂函数型本构模型；黑色直线代表线性本构模型[6]。

图7 SGL 10BA 碳纤维纸Fig. 7 SGL 10BA carbon fiber paper

对于SGL10BA 型碳纤维纸，纸的一侧存在5%的聚四氟乙烯薄膜，因此碳纤维的实际含量需要相应减少5%，修正因子γ 为0.95。观察图7(b)所示的应力-应变曲线可以发现，本文所提出的对数型本构模型与幂函数型本构模型依旧具有一定的准确性，两种模型的准确性均优于传统的线性本构模型。

基于以上结果对比，可以发现碳纤维纸中的额外树脂对其应力-应变关系有比较明显的影响。对比三种本构模型的预测结果，全面考虑接触变形、弯曲变形与孔隙率变化的对数型本构模型具有较好的鲁棒性，适用范围较广；简化后的幂函数型本构模型能够较准确的预测树脂含量较低的碳纤维纸的应力-应变关系；传统的线性本构模型[6]对各种碳纤维纸的预测均有较大的误差。本文提出的两种本构模型的相对误差均远小于线性本构模型，其主要原因为线性本构未考虑压缩过程中纤维因接触点增加而发生的应变强化。对比本文提出的两种本构模型，可以发现对数型模型的误差最小，更适合含有机物碳纤维纸的理论与数值分析。

不过，在有额外有机物的条件下，三种本构模型与实验结果的相对误差均高于无有机物碳纤维纸，说明虽然有机物在碳纤维纸中含量较低，仍对其微观变形机制和压缩本构关系有较大影响。有必要根据碳纤维纸内的有机物分布，进行更深入的研究。

3 结论

本文针对碳纤维纸微结构随机多孔的特点，提出了两种非线性面外压缩本构模型，分别是全面考虑微结构孔隙率、单胞尺寸变化及微结构弯曲与接触变形的对数型模型，以及忽略赫兹接触变形的幂函数型简化模型。基于三种商品化碳纤维纸的应力-应变实验数据，将本文提出的两种本构模型与传统的线性本构模型进行分析对比，得到以下结论：

(1) 对于本文涉及的多种典型碳纤维纸，作者提出的幂函数型本构模型与对数型本构模型均能够描述材料的非线性应力-应变关系。两种本构模型与实验数据的相对误差均明显小于传统的线性本构模型，可以为相关结构的仿真分析提供有效的本构模型基础。

(2) 对于含树脂等有机物较少的碳纤维纸，幂函数型本构模型与实验数据的相对误差最小，在最优加紧应力附近的相对误差可低至0.1%。同时，由于幂函数型本构模型具有参数少的特点，仅通过一对应力-应变数据即可建立幂函数型本构，比较便于通过实验测试得到。

(3) 对于含树脂等有机物较多的碳纤维，由于这些随机分布有机物的模量较低，降低了碳纤维纸的等效模量。本文引入了体积分数修正因子 γ，以考虑有机物杂质对碳纤维纸本构模型的影响。与实验数据的对比显示，修正后的两种本构模型与实验数据的相对误差仍明显小于传统的线性本构模型，但其相对误差的离散性较大，仍有待于进一步修正。