考虑温度效应的干纤维预制体压缩蠕变模型

司衍鹏，孙立帅，闫恩玮，2，李玉军，蒋建军，

1.西北工业大学 机电学院，西安 710072

2.中航西安飞机工业集团股份有限公司，西安 710089

干纤维自动铺放-液体成型是一种低成本、高效率的树脂基纤维增强复合材料成型技术体系，在近10 年发展极其迅速，应用该成型工艺制造的零件性能与采用热压罐成型制造的零件性能相当，且适用于大型复杂复合材料结构的整体制造［1-4］。该成型工艺将自动铺放和液体成型的优势综合起来，提高了复合材料结构的减重效率、产品精度和表面质量。

越来越多的国家和组织启动研究项目开展相关研究，例如，俄罗斯已将该成型技术应用于MC-21 单通道客机的中央翼盒、翼壁蒙皮和翼梁结构制造［5］；空客公司在推出的“明日之翼”研究项目中使用了该成型工艺，为下一代支线客机机翼全面采用液体成型工艺做技术验证，目前已实现17 m 验证件测试［6］。

干纤维自动铺放-液体成型主要包括干纤维预成型、树脂填充、固化3 个方面，其中干纤维预成型过程是自动铺放成型的基本过程，干纤维预成型决定了结构成型后的纤维构型，且预制体的质量对制造的复合材料零件的力学性能有重大影响［7］。在预成型过程中，纤维束间的相互挤压引起纤维束截面形状和纤维取向发生变化，进而改变预成型体的纤维束路径及孔隙率，从而影响后续工艺中树脂的渗透性能和成型后结构的力学性能。因此，研究干纤维预制体在预成型过程中的预制体厚度控制具有重要意义。

在干纤维预成型过程中，纤维丝束被挤压在一起，起到弹性弯曲梁的作用，纤维丝束间相互滑动，重新排列，产生摩擦力。上述综合效应导致了干纤维预制体在压实过程中复杂的黏弹性行为［8］。干纤维预制体的黏弹性响应包括蠕变、回复和应力松弛。当移除产生蠕变的载荷后，材料会出现瞬时弹性回复和随时间发展的滞后回复现象，这个过程被称为蠕变回复。

对于干纤维预制体而言，目前相关研究多数集中于应力松弛［9-11］，且对干纤维预成型体在预成型过程中的蠕变/回复行为研究较少，而在实际的干纤维预成型中，保压、释放过程会发生蠕变现象，这使得其在实际应用中仍然是一个悬而未决的问题，特别是在干纤维在预成型过程中，预制体厚度方向表现出非线性、时间相关性和温度相关性，需要一种能够描述干纤维预制件蠕变/回复行为的材料本构模型，该模型应可以准确预测在不同蠕变应力和预成型温度下预制体厚向应变随时间的变化。Celauro 等［12］使用双参数模型（TP）分数阶函数表达式可以准确地反映沥青混合料的蠕变/回复行为，但在较长加载时间下，预测的回复应变与测量数据有较大偏差。Findley 模型（FD）在预测聚合物基复合材料的黏弹性能方面表现出良好的适用性，但它没有考虑材料形状的变化，并且模型的参数对曲线的形状十分敏感［13-14］。由N个（N＞3）Kelvin 模型（GK）并联或串联组成的广义Kelvin 模型可以表示黏弹性材料的黏度和弹性，但由于它包含的参数较多，在实际求解过程相对复杂［15］。

迄今为止，用于分析蠕变/回复非线性黏弹性行为最广泛的方程之一是Burgers 模型（BG）［16］。先前研究表明，干纤维预制体在不同工艺条件下表现出显著的非线性黏弹性响应［17］。Burgers 模型相对于其他模型能更好反映线性黏弹性，具有模型参数少、模拟精度高的特点［13，18-20］。但该模型仅可在特定实验条件下较为准确地描述材料的非线性黏弹性行为，无法用单一模型预测不同蠕变应力和预成型温度下的蠕变/回复行为。综上所述，现有模型主要关注蠕变或回复行为，通常这些模型针对蠕变/回复过程提供2 个单独公式，且准确描述蠕变阶段到回复阶段的过渡较为困难。此外，上述模型均未在不同蠕变应力和预成型温度下，对卸载后材料厚向应变随时间的变化进行分析和预测。由此可见，研究使用单一方程中对材料行为的全面描述对于预测干纤维预制体在不同蠕变应力和预成型温度下的蠕变/回复行为至关重要。

通过研究分析干纤维预制体随时间变化的蠕变/回复行为，对比分析4 种蠕变模型在不同条件下的蠕变/回复行为的预测效果，选用参数较少且预测精度较高的Burgers 模型；在此基础上提出1种单一方程形式，且考虑不同蠕变应力和预成型温度的改进的Burgers 模型，然后通过实验数据分析计算及最小二乘法原理获取模型参数；最终结合新设计的实验结果验证该本构模型的准确性。

1 实验材料及测试方法

1.1 实验材料

如图1 所示，使用的纤维增强复合材料是威海光威公司生产的碳纤维机织物CF3031。材料的规格列于表1 中。该织物为3×1 斜纹织物，丝束尺寸为3 K，面积重量为（220±7）g/m2。

表1 机织物CF3031 规格Table 1 Specification of woven fabric CF3031

图1 碳纤维机织物CF3031Fig.1 Carbon fiber woven fabric CF3031

1.2 蠕变-回复实验方法

使用计算机数字控制机床切割成单层边长为50 mm 的方形样品，以0°的方向铺设，形成24 层的预成型体叠层。在干燥环境下研究纤维织物的蠕变/回复行为。图2 为放置在两钢质压板之间的干纤维预制体样品，这两个压板为圆形，直径为100 mm。

图2 蠕变/回复实验装置Fig.2 Creep/recovery experiment apparatus

使用ETM104C 型微机控制电子万能试验机对实验件进行蠕变/回复实验。由于干纤维预制体在实际预成型过程中所受铺放压力一般小于0.8 MPa，CF3031 织物在温度高于60 ℃时，表面的定型剂易发生固化，因此设计蠕变/回复实验方案如表2 所示，包括在室温（25 ℃）不同蠕变应力（0.3、0.5、0.8 MPa）和不同预成型温度（40、50、60 ℃）下同一蠕变应力（0.3 MPa）时的实验，以1 mm/min 的加载速率对试样施加位移控制，用量程为10 kN 的传感器记录试样应变随时间的变化。蠕变/回复实验可分为4 个阶段：加载阶段、蠕变阶段、卸载阶段和蠕变回复阶段。

表2 蠕变/回复实验的测试方案Table 2 Testing programme for creep/recovery experiments

1）加载阶段

在应力加载到预定载荷的过程中，样品主要发生不超过6 min 的弹性变形和非机械变形。

2）蠕变阶段

为使试样充分地蠕变，试样主要承受持续约3 h的恒定载荷，伴随着黏弹性变形和黏塑性变形。

3）卸载阶段

加载3 h 后开始卸载，约150 s，样品主要发生弹性回复。

4）蠕变回复阶段

样品进入自由状态，施加约10 N 的载荷1 h，样品主要发生黏弹性变形。

材料的蠕变/回复过程如图3 所示。

图3 恒定施加应力下蠕变/回复实验的典型应变曲线Fig.3 Typical strain curve from creep/recovery experiment under constant applied stress

2 蠕变/回复计算

2.1 蠕变模型

为了实现复合材料蠕变行为的计算分析，将目前文献中采用的计算模型进行总结，选用Findley 模型、双参数模型、广义Kelvin 模型和Burgers模型4种模型。其中Findley模型应变计算公式为

双参数模型应变计算公式为

式中：a为瞬时应变；t为时间；b为瞬态蠕变的振幅；c为材料常数；ε（t）为t时刻的总材料应变。

广义Kelvin 模型和Burgers 模型均由弹簧、Kelvin 单元和黏壶组成，其模型示意图分别如图4 和图5 所示。

图4 广义 Kelvin 模型示意图Fig.4 Diagram of generalized Kelvin model

图5 Burgers 模型示意图Fig.5 Diagram of Burgers model

广义Kelvin 模型应变计算公式为

式中：n′为广义Kelvin 模型中Kelvin 单元的数量。

Burgers 模型应变计算公式为

式中：σ为施加在材料上的恒定应力；E0为材料弹性模量；η0为材料的黏度系数；Ei和ηi为每个Kelvin 单元的弹性系数和黏度系数；τi为每个Kelvin单元的松弛时间，且τi=ηi/Ei，i=0，1，…。

2.2 改进的Burgers 模型

应力σ（t）和应变ε（t）之间内在关系的相应数学公式包含具有指数解的整数阶时间导数，即其具体公式通常为

式中：n为应力σ（t）方程中整数阶导数的总数，可被理解为应力方程中包含的导数的最高阶数；m为应变ε（t）方程中整数阶导数的总数，可被理解为应变方程中包含的导数的最高阶数；i和j为阶数，i=1，2，…，n，j=1，2，…，m。

式中：Ki和Si是常数值。其中一个例子是图4 所示的Burgers 模型，它包括一个Kelvin-Voigt 模型（一个与黏壶平行的弹簧）和一个Maxwell 元件（与黏壶串联的弹簧）。Burgers 模型中σ和ε之间的关系可以表示为

Burgers 模型的蠕变应变计算公式为

式中：τ1=η1/E1为松弛时间。

现有研究将蠕变/回复变形划分为可回复和不可回复的部分，一般采用黏弹性模型描述材料的可回复的变形，将不可回复的变形用其他部件表示。在本实验过程中，不同实验条件下材料特性、纤维间孔隙等所产生的非机械应变较为一致，因此将非机械应变元件与原黏弹性模型进行串联，使得改进的Burgers 模型能够考虑非机械应变。在卸载点t0之前，材料在压实过程中由于恒定的应力而发生蠕变应变ε分解为可逆部分和不可逆部分2部分：

式中：εr为可逆部分，是黏弹性应变；εN为不可逆部分，是非机械应变。如图6 所示，由于在压实过程中且未达到蠕变应力σ0之前，干纤维预制体的应变为弹性应变εe与非机械应变εN之和，在卸载过程中产生的应变均为弹性应变εe。

图6 蠕变/回复应变曲线示意图Fig.6 Schematic diagram of creep/recovery strain curve

如图7 所示，在具有Kelvin-Voigt 元件、Maxwell 元件和非机械应变元件的改进的Burgers 模型中，应变被描述为时间t和蠕变应力σ的函数：

图7 改进的Burgers 模型示意图Fig.7 Schematic representation of modified Burgers model

CF3031 干纤维预制体的不稳定性主要体现在2 个方面：一方面材料表面上存在定型剂，这种热敏黏弹性材料使得干纤维预制体的蠕变/回复行为对预成型温度变化十分敏感；另一方面，不同载荷下的预制体纤维丝束间嵌套程度不同。上述问题均会影响材料的蠕变速率。此外，干纤维预制体在实际预成型过程中的黏弹性力学性能不仅与时间有关，还与压辊力和预成型温度有关，而传统的黏弹性本构方程不能直接反映蠕变应力和预成型温度对材料黏弹性特性的影响，例如从Burgers 模型的本构关系中很容易看出，该模型只能表征材料在给定蠕变应力、预成型温度下随时间变化的黏弹性力学特性。由于预成型过程中预制体所受载荷和预成型温度效应的复杂性，研究干纤维预制件的本构模型非常困难。

为解决这一问题，通过进行相同预成型温度不同蠕变应力和同一蠕变应力不同预成型温度下的干纤维预制体蠕变/回复实验，并根据实验数据拟合获取以蠕变应力和预成型温度为变量的材料本构方程系数。最终建立考虑不同蠕变应力和预成型温度的改进的Burgers本构方程形式为

式中：E1和η1为常数；E2（σ，T）和η2（σ，T）为载荷σ和预成型温度T的函数。

2.3 干纤维预制体蠕变回复计算方法

目前，还没有专门面向干纤维预制体的蠕变回复计算方法。参考混凝土等黏弹性材料的蠕变响应［21］，使用Boltzmann 叠加原理建立干纤维预制体蠕变回复方程。Boltzmann 叠加原理认为，当多个载荷共同作用时，每一个负载的增加或减少步骤都是独立的，且能够相互叠加［22］。基于Boltzmann 叠加的应力、应变表示为

式中：n′′为独立应力增量的数量；τ为弛豫时间，决定材料响应速度，与整体时间t有所区别，τ只是一个特定瞬时时间点；Δσ（τi）为在弛豫时间τi处引起的应力增量；φ（t，τ）为单位应力下的蠕变。如图8所示，根据叠加原理，蠕变模型可用于计算和分析蠕变回复，即卸载后变形回复相当于反向力蠕变叠加。

图8 Boltzmann 叠加原理在蠕变回复行为分析中的应用Fig.8 Application of Boltzmann superposition principle in analysis of creep-recovery behavior

假定在t0时刻撤除施加在材料表面的恒定应力，基于叠加原理方法，Findley 模型的蠕变回复方程为

双参数模型的蠕变回复方程为

广义Kelvin 模型的蠕变回复方程为

Burgers 模型的蠕变回复方程为

3 计算结果分析与讨论

3.1 干纤维预制体蠕变计算结果

根据最小二乘法原理，编制了非线性拟合程序，在程序中输入合适的模型参数初始值，利用实验数据拟合得到模型各参数值及相关系数，并采用相关系数评判模型拟合精度。图9、图10 中的曲线为分别采用Findley 模型、Burgers 模型、广义Kelvin 模型、双参数模型，在室温（25 ℃）下不同应力（0.3、0.5、0.8 MPa）和同一应力（0.3 MPa）不同预成型温度（25、40、50、60 ℃）下的蠕变曲线的非线性拟合结果，拟合得到的模型参数如表3～表6 所示，表中R2为拟合优度，RMSE 为均方根误差。

表3 Findley 模型蠕变曲线拟合参数Table 3 Fitting parameters of creep curve of Findley model

表4 广义Kelvin 模型蠕变曲线拟合参数Table 4 Fitting parameters of creep curve of generalized Kelvin model

表5 Burgers 模型蠕变曲线拟合参数Table 5 Fitting parameters of creep curve of Burgers model

表6 双参数模型蠕变曲线拟合参数Table 6 Fitting parameters of creep curve of two-parameters model

图9 模型在室温、不同蠕变应力下蠕变计算结果Fig.9 Creep calculation results of model under room temperature and different creep stresses

图10 模型在0.3 MPa、不同预成型温度下蠕变计算结果Fig.10 Creep calculation results of model at 0.3 MPa and different preform temperatures

结合表3～表6 中模型拟合参数和图10 中各模型拟合结果可得，Burgers 模型在不同预成型温度下都表现出了良好的拟合效果，随着预成型温度的增加，Findley 模型、广义Kelvin 模型和双参数模型拟合效果受预成型温度影响较大。

通过不同模型蠕变曲线拟合结果，计算拟合优度R2和均方根误差RMSE。由图11 中的柱状图中R2和RMSE 计算结果可以看出，在室温不同蠕变应力和相同蠕变应力不同预成型温度条件下，Findley 模型和双参数模型拟合效果较差，R2在0.8左右，且RMSE 较大，广义Kelvin 模型拟合效果次之，Burgers 模型的拟合效果最好，R2＞0.95，且RMSE 最小。

图11 不同模型拟合蠕变曲线的R2和RMSE 计算结果Fig.11 Calculation results of R2 and RMSE of different models fitting creep curves

3.2 干纤维预制体蠕变/回复计算结果

基于Boltzmann 叠加原理，采用两参数指数模型、Findley 模型、广义Kelvin 模型和 Burgers模型对蠕变/回复实验数据进行计算，分析得到的蠕变/回复全曲线如图12 所示。这些曲线表明，CF3031 干纤维织物预制体的应变与时间、载荷和预成型温度有关。随着时间、载荷和预成型温度的增加，由于预制体的黏弹性特征，应变在蠕变阶段逐渐增加。由于非机械应变的存在，在不同蠕变应力和不同预成型温度下的蠕变/回复曲线末端，永久变形始终存在。此外，蠕变/回复行为对预成型温度变化较为敏感，随着预成型温度的升高，由于单层干纤维织物表面的定型剂由固态向熔融态转变，材料在压实过程中层间间距逐渐变小，从而导致应变增大。

图12 不同模型在不同蠕变应力和不同预成型温度下蠕变/回复计算结果Fig.12 Creep/recovery calculation results of different models under different creep stresses and preform temperatures

图12 中的曲线图为分别采用Findley 模型、Burgers 模型、广义Kelvin 模型和双参数模型，在室温（25 ℃）、不同蠕变应力（0.3、0.5、0.8 MPa）和相同蠕变应力（0.3 MPa）、不同预成型温度（25、40、50、60 ℃）条件下的蠕变/回复曲线的非线性预测结果，拟合得到的模型参数如表7～表10所示。

表7 广义Kelvin 模型蠕变/回复曲线拟合参数Table 7 Fitting parameters of creep/recovery curve of generalized Kelvin model

表8 Burgers 模型蠕变/回复曲线拟合参数Table 8 Fitting parameters of creep/recovery curve of Burgers model

表9 Findley 模型蠕变/回复曲线拟合参数Table 9 Fitting parameters of creep/recovery curve of Findley model

表10 双参数模型蠕变/回复曲线拟合参数Table 10 Fitting parameters of creep/recovery curve of tow-parameter model

通过不同模型蠕变/回复曲线拟合结果，计算得到各模型的R2和RMSE 如图13 所示。结合表7～表10 中模型拟合参数和图12 中各模型拟合结果可得，在室温不同蠕变应力条件下，4 种模型预测值与实验数据吻合较好，R2均较高，且Burgers 模型的RMSE 最小。在相同蠕变应力不同预成型温度条件下，Findley 模型、双参数模型和广义Kelvin 模型随着预成型温度的升高拟合效果逐渐变差，其中广义Kelvin 模型拟合效果最差，R2均值在0.8 左右，且RMSE 较大。

图13 不同模型拟合蠕变/回复曲线的R2和RMSE 计算结果Fig.13 Calculation results of R2 and RMSE of different models fitting creep/recovery curves

可以清晰看出，在不同实验条件下采用Burgers 模型均能够良好分析CF3031 干纤维织物预制体的蠕变/回复行为，且整体拟合效果优于其他模型，R2均值大于0.95，且RMSE 最小。

4 本构方程建立与验证

4.1 参数识别

基于对不同模型拟合效果对比，选用拟合效果最好的Burgers 模型，对传统Burgers 模型进行改进。与传统Burgers 模型相比，改进的Burgers模型引入了蠕变应力σ0和预成型温度T作为经验本构方程系数的变量，并引入由于材料特性、纤维间孔隙等所产生的非机械应变εN。其蠕变/回复全过程应变描述为

式中：ε1（t）为从压实过程到卸载结束时σ0增加的应变；根据叠加原理将t0时刻卸载产生应变等效为施加反向力σ0产生的应变ε2（t-t0），则ε1（t）-ε2（t-t0）是卸载响应阶段的总应变；H（x）为阶跃函数，其函数表达式为

现有研究通过实验数据与拟合结果分析，基于叠加原理方法的材料模型均会高估材料的回复能力［23］，一方面是因为CF3031 干纤维织物表面的定型剂为热敏黏弹性材料，定型剂随着预成型温度升高逐渐由固态转变为熔融态，导致各干纤维织物层间相互粘连；另一方面蠕变应力的增加会加剧纤维丝束间的嵌套行为，因此在回复过程中无法将t0时刻卸载产生的应变，等效为施加与蠕变应力大小相同方向相反力σ0所产生的应变。为考虑不同应力和预成型温度对回复行为产生的影响，引入系数A，定义系数A为蠕变应力σ0和预成型温度T的函数。基于上述分析，单一方程形式的CF3031 干纤维预制体蠕变/回复材料模型为

图14、图15 分别为室温不同蠕变应力下蠕变柔量随时间变化的曲线和相同蠕变应力不同预成型温度下蠕变柔量随时间变化曲线。

图14 不同蠕变应力下蠕变柔量随时间变化曲线Fig.14 Creep compliance versus time under different creep stresses

图15 不同预成型温度下蠕变柔量随时间变化曲线Fig.15 Creep compliance versus time at different preform temperatures

蠕变柔量J（t）计算公式为

如图14、图15 所示，蠕变柔量J（t）在加载过程中均随着时间的延长逐渐增大，在卸载阶段蠕变柔量均随着时间的增加逐渐减少。加载阶段蠕变应力越大、预成型温度越高，蠕变柔量变化越小；卸载阶段蠕变应力越小、预成型温度越低，蠕变柔量变化越大。

由图12 中蠕变/回复实验数据以及图15 参数识别示意图获取E0、εN、η0的数值。

根据不同蠕变应力和不同预成型温度下蠕变/回复实验数据中非机械应变εN和卸载纯弹性应变εe计算E1、εN，发现不同实验条件下E0、εN计算结果趋于一致，且在拟合过程中发现η1的大小与图16 中红色椭圆形标注的曲线倾斜程度有关，随着η1的增大，此部分拟合曲线与x轴的夹角越小，随着时间的增加此部分曲线逐渐趋于与x轴平行，根据拟合结果取η0为一较大值（大于100）。因此本构模型系数E0、εN的计算结果和η0的取值为

图16 参数E0、εN、η1识别示意图Fig.16 Schematic diagram of parameters E0，εN，η1 identification

为获取本构模型中E1、τ1、A的表达式，首先将蠕变应力σ0和T进行归一化处理，其次选用二次多项式作为模型系数E1、τ1、A的表达式形式，最后根据不同实验条件下的实验结果对二次多项式的系数进行求解，求得E1、τ1、A与蠕变应力σ0和预成型温度T相关的函数表达式为

式中：σ0∈（-1，1），T∈（-1，1）。将式（22）、式（23）代入式（20），得到能够预测CF3031 干纤维预制体在不同蠕变应力和预成型温度下蠕变/回复行为的单一方程形式本构模型。

将T=25 ℃和σ0=0.3 MPa 归一化后，分别代入得到的本构方程式（20）中，采用最小二乘法原理分别对室温下不同蠕变应力和相同应力在不同预成型温度下的蠕变柔量进行拟合，得到的蠕变柔量三维拟合图如图17、图18 所示，拟合优度R2均大于0.97，说明不论从数据的大小还是数据的变化趋势，实验结果与拟合曲面均表现出良好的一致性。

图17 室温下不同蠕变应力的蠕变柔量三维拟合图Fig.17 Three dimensional fitting diagram of creep compliance of different creep stresses at room temperature

图18 相同蠕变应力在不同预成型温度下的蠕变柔量三维拟合图Fig.18 Three dimensional fitting diagram of creep compliance of the same creep stress at different preform temperatures

4.2 模型验证

为了验证建立的CF3031 干纤维预制体的非线性黏弹性本构模型的准确性，采用与实验相同的干纤维预制体样品，以蠕变应力σ0、回复保载应力σ1、预成型温度T、卸载时刻t0为变量进行实验设计，模型验证实验方案如表11 所示。

表11 本构方程验证实验的测试程序Table 11 Test procedure of constitutive equation validation experiment

图19 为本构模型预测结果与实验数据的比较结果，验证实验与模型预测曲线基本吻合。结果表明，该单一方程形式的材料模型能够实现在不同蠕变应力和预成型温度下干纤维预制体蠕变/回复的良好预测，这对于实际预成型工艺过程中材料的非线性黏弹性行为的预测较为便利。

图19 预测曲线与实验曲线对比Fig.19 Comparison between prediction curves and experimental curves

5 结论

采用了一系列蠕变/回复实验来研究CF3031 干纤维预制件在不同蠕变应力和预成型温度下的黏弹性行为，对于干纤维预制件的蠕变/回复行为进行计算分析。特别是，蠕变应力和预成型温度对CF3031 干纤维预制体的黏弹性行为的影响已经通过使用改进的Burgers 模型进行了预测与量化。主要研究结果列举如下。

1）采用4种模型进行了CF3031干纤维预制体的蠕变计算分析，发现Burgers 模型在不同的蠕变应力和预成型温度下的预制体蠕变的计算效果略优于广义Kelvin模型、Findely模型和双参数模型。

2）基于Burgers 模型和Boltzmann 叠加原理，提出了考虑蠕变应力和预成型温度的经验本构模型，仅使用单一方程即可全面描述CF3031干纤维预制体的在不同的蠕变应力和预成型温度下的黏弹性行为。

3）改进的Burgers 模型能够被用来分析预测CF3031 干纤维预制体蠕变/回复的黏弹性行为，且实验结果与预测结果吻合较好，证明了该本构模型的可行性。