1. 引言
近年来，随着工业的迅速发展，纤维增强复合材料以其具有良好的强度与刚度等特性广泛应用于航空航天工业与新能源发电中[1]-[3]。但传统的纤维增强复合材料存在着韧性差，易出现裂纹等缺点[4]，使得纤维增强复合材料的适用具有较大的局限性，目前纤维杂化是提高复合材料力学性能的有效途径之一。通过使用两种或更多种纤维作为增强相，可以将多种纤维的优点结合在一起以弥补传统纤维增强复合材料的不足[5] [6]。然而由于杂化复合材料强度和刚度的变化与多种变量有关，如果通过实验方法对材料属性进行验证，往往需要大量的实验次数，成本昂贵的同时还会造成资源的浪费，因此非常必要提出一种高效的对复合材料进行多变量数据分析的方法。
目前有诸多学者对纤维增强复合材料的多个方面进行了研究。例如在模型的生成上，Chen [7]等人通过生成周期性空间排列模型对纤维素微原纤维增强生物复合材料的载荷传递特性进行研究。虽然周期性分布保留了纤维的确定性和稳定性等原有特性，能够极大简化分析过程，但是在实际计算中这种稳定性可能会对强度预测和疲劳寿命的结果产生极大地影响。在算法的选择上，目前主要有初始周期性振动模型、硬核摇动模型和随机顺序扩展三种方法，Buryachenko [8]等人对初始周期性振动模型和硬核摇动模型两种方法进行了详细的阐述，Wang [9]等人提出了一种单向长纤维增强复合材料硬核摇动模型算法，研究了纤维分布的随机性对横向模量的影响。然而，每个算法也有其一定的局限性，例如随机顺序扩展算法能够有效地提升模型的总体积分数，但是对于杂化材料会出现较大的限制；硬核模型虽然自然，但是却可能出现大小纤维束集中的情况，因此对目前现有算法进行改进显得十分有必要。
基于此，本文对碳纤维/玻璃纤维杂化复合材料进行研究，提出一种将随机顺序扩展算法与硬核模型算法相结合的方法进行生成随机性模型，利用控制变量法研究了总体积分数和大小径纤维数量比两个参数对材料力学性能的影响，同时对提出的将随机顺序扩展算法与硬核模型算法相结合的算法进行可靠性分析，进一步验证了算法的可靠性。
2. 仿真前处理
2.1. 材料属性的确定
在本研究中，主要需要树脂填充层、大直径纤维与小直径纤维三种材料的属性进行确定。Cao [10]等人给出了一些常用材料及其工程常数，通过对比分析，本文选取了树脂填充层材料为环氧树脂8552，大径纤维选用材料为玻璃纤维，小径纤维选用碳纤维IM7，其工程常数如表1所示。
Table 1. Material engineering constants
表1. 所需材料工程常数
弹性常数 |
碳纤维 |
玻璃纤维 |
环氧树脂 |
轴向的弹性模量E1/GPa |
276 |
74 |
4.08 |
横向的弹性模量E2/GPa |
19 |
74 |
4.08 |
面内剪切模量G12/GPa |
27 |
30.8 |
1.478 |
主泊松比ν12 |
0.2 |
0.2 |
0.38 |
主泊松比ν23 |
0.2 |
0.2 |
0.38 |
纵横剪切模量G12/GPa |
7 |
30.8 |
1.478 |
其中ν23的值需要额外计算，计算公式如式(1)所示。
(1)
除此之外，由于玻璃纤维与环氧树脂是各向同性的关系，因此在计算时不需要计算两者的ν23。
2.2. 仿真参数的选择
由于本次研究中选取大径纤维与小径纤维的材料为玻璃纤维和碳纤维，两者的直径分别为11 μm和5 μm，因此在微观模型中模拟的直径分别定义为10个单元和5个单元。在释放掉变量大小纤维直径比之后，在进行模型分析时主要将模型分为两个参数控制的模型：其一为总体积分数，其二为大小纤维束的数量比。
2.3. 边界条件的确定
在本研究中，边界条件设置为在三个零位移面上施加位移约束，即面位移为0，并分别对零位移面的对面进行位移的施加，分别为0.1 μm，10 μm，10 μm，位移量根据模型的长宽与厚度的1/10进行控制。同时，体积单元所有平面在变形后仍保持为平面。
2.4. 仿真算法的确定
本文采用一种将顺序扩展算法和硬核模型算法相结合的全新算法生成模型，在生成大圆直径时主要是采用顺序扩展算法，该算法起始点的设立会直接影响总体积分数，目前顺序扩展算法中常见的初始圆心有四种分布方法：模型几何中心、四个坐标极点、底边中心点与随机起始点。由于所研究的材料为杂化材料，几何中心与底边中心点的分布方法更佳，本文采用以矩形的几何中心作为起始圆心的方法生成模型，相较于底边中心点的方法能够更好的展示随机性。在小圆圆心生成则采用了硬核模型算法，具体算法流程图如图1所示。
Figure 1. Algorithm flowchart
图1. 算法流程图
对于随机生成圆心是否达标的判定条件主要通过三个因素：① 是否满足在随机生成的范围内；② 是否与已经生成的大圆相交；③ 返检次数是否已经超过了设计要求。
3. 数值模拟结果与分析
本节通过控制变量法对由于参数变化带来的模型力学属性的改变进行分析，主要参数为大小径纤维数量比和总体积分数。
3.1. 总体积分数对模型力学性能影响结果
当大径纤维与小径纤维的数量比为2:1时，分别改变总体积分数的值(0.4、0.5和0.6)，建立的三种模型图、变形云图应力观测图如图2，图3所示。
(a) 总体积分数0.4 (b) 总体积分数0.5 (c) 总体积分数0.6
Figure 2. Different population integral number models generate graphs
图2. 不同总体积分数模型生成图
(a) 总体积分数0.4 (b) 总体积分数0.5 (c) 总体积分数0.6
Figure 3. Deformation nephogram stress observation map of different models
图3. 不同模型变形云图应力观测图
在相同条件下对每个总体积分数随机生成三个模型，对每个模型所得的xy平面、yz平面、xz平面的反作用力进行求和处理，分别记作轴向、径向和垂向并进行求和计算，再将求和值带入到弹性模量的计算公式中进行计算，对此操作进行三次反复可得表格数据如表2所示。
Table 2. Macroscopic elastic constants of hybrid composites with three overall integral number models
表2. 三种总体积分数模型杂化复合材料宏观弹性常数
总体积分数 |
方向 |
弹性常数 |
模型一 |
模型二 |
模型三 |
0.4 |
轴向 |
弹性模量E/GPa |
53.900 |
53.954 |
53.397 |
泊松比υ |
0.3063 |
0.3065 |
0.3055 |
径向 |
弹性模量E/GPa |
8.2131 |
8.1326 |
8.1989 |
泊松比υ |
0.4323 |
0.4251 |
0.4325 |
垂向 |
弹性模量E/GPa |
8.2070 |
8.1773 |
8.2006 |
续表
|
|
泊松比υ |
0.4206 |
0.4359 |
0.4048 |
0.5 |
轴向 |
弹性模量E/GPa |
53.900 |
53.954 |
53.397 |
泊松比υ |
0.3063 |
0.3065 |
0.3055 |
径向 |
弹性模量E/GPa |
8.2131 |
8.1326 |
8.1989 |
泊松比υ |
0.4323 |
0.4251 |
0.4325 |
垂向 |
弹性模量E/GPa |
8.2070 |
8.1773 |
8.2006 |
泊松比υ |
0.4206 |
0.4359 |
0.4048 |
0.6 |
轴向 |
弹性模量E/GPa |
53.900 |
53.954 |
53.397 |
泊松比υ |
0.3063 |
0.3065 |
0.3055 |
径向 |
弹性模量E/GPa |
8.2131 |
8.1326 |
8.1989 |
泊松比υ |
0.4323 |
0.4251 |
0.4325 |
垂向 |
弹性模量E/GPa |
8.2070 |
8.1773 |
8.2006 |
泊松比υ |
0.4206 |
0.4359 |
0.4048 |
由表2可知，在总体积分数为0.4的杂化纤维模型中，材料的轴向平均弹性模量为53.75 Gpa，最大偏差值仅为千分之六，通过算法生成的模型力学性能较为平均，并且Trias [10]给出这个偏差值是较小的。并且在总体积分数为0.5和0.6的杂化纤维模型中，材料的平均轴向弹性模量64.142 GPa和73.35 GPa，最大偏差值分别为0.15%和0.11%，这也说明了当纤维排布越紧密时，此算法的稳定性越高，极其适合应用于高比例总体积分数的模型生成。此外还不难发现，模型的轴向弹性模量大小和纤维的总体积分数呈正相关，在特定大径纤维与小径纤维数量比例下的杂化纤维模型，其抵抗外界作用力所产生的形变就越小，并且在此算法下，这个数值变化是偏向线性的。
通过观察总体积分数为0.4，0.5和0.6的三种杂化纤维模型，可以看出径向与垂向的弹性模量与泊松比的平均值相近，差值都小于0.5%。但是两个数值有极性分布的特点，即当两个方向中有一个方向上的泊松比或是弹性模量远大于或者远小于平均值时，与其垂直的方向上的弹性模量与泊松比往往是更贴近于平均值。这就意味着虽然两个相互垂直方向上的纤维对力学参数的抵抗作用和排列方式有关，但是往往会有趋同性，即两者同时趋大或者趋小。而这种趋势是与轴向的力学性能直接挂钩的：当轴向力学性能良好时，其径向与垂向的力学性能就会偏差，反之亦然。因此虽然大小径纤维数量比例相同且模型的总体积分数相同，但是模型的力学性能仍然和算法与模型有直接的关联，而这种关系与体积分数变化并没有直接关系。因此即使是相同参数下的模型，也会出现模型适用环境的差别。
3.2. 大小径纤维数量比对模型力学性能影响结果
当总体积分数为0.5时，分别改变大小径纤维数量比的值(1:2、1:2.5和1:3)，建立的三种模型图、变形前后应力对比图如图4，图5所示。
由于在总体积分数确定，改变大小径纤维数量比情况下主要需要进行分析比较的数据是轴向的弹性模量，因此在进行相同的总体体积分数不同的纤维数量比的数据比较中，对轴向弹性模量进行计算，不对周向和垂向弹性模量进行计算。表3为各模型轴向弹性模量计算表。
当总体积分数为0.55时，同样的，分别改变大小径纤维数量比的值(1:2、1:2.5和1:3)，建立的三种模型图、变形前后应力对比图如图6，图7所示。
(a) 大小径纤维数量比1:2 (b) 大小径纤维数量比1:2.5 (c) 大小径纤维数量比1:3
Figure 4. The model of fiber quantity ratio of different large small path was generated
图4. 不同大小径纤维数量比模型生成图
(a) 大小径纤维数量比1:2 (b) 大小径纤维数量比1:2.5 (c) 大小径纤维数量比1:3
Figure 5. Stress comparison diagram before and after deformation
图5. 变形前后应力对比图
Table 3. Macroscopic elastic constant of hybrid composite materials
表3. 杂化复合材料宏观弹性常数
|
模型一 |
模型二 |
模型三 |
大小径纤维数量比 |
1:2 |
1:2.5 |
1:3 |
弹性模量E(GPa) |
64.090 |
68.078 |
71.118 |
泊松比υ |
0.2932 |
0.2919 |
0.2883 |
(a) 大小径纤维数量比1:2 (b) 大小径纤维数量比1:2.5 (c) 大小径纤维数量比1:3
Figure 6. The model of fiber quantity ratio of different large small path was generated
图6. 不同大小径纤维数量比模型生成图
(a) 大小径纤维数量比1:2 (b) 大小径纤维数量比1:2.5 (c) 大小径纤维数量比1:3
Figure 7. Stress comparison diagram before and after deformation
图7. 变形前后应力对比图
表4为在总体积分数为0.55时，各模型轴向弹性模量计算表。
Table 4. Macroscopic elastic constant of hybrid composite materials
表4. 杂化复合材料宏观弹性常数
|
模型一 |
模型二 |
模型三 |
大小径纤维数量比 |
1:2 |
1:2.5 |
1:3 |
弹性模量E(GPa) |
68.906 |
72.76 |
75.567 |
泊松比υ |
0.2862 |
0.2838 |
0.2832 |
通过比较相同总体积分数下不同纤维数量比的模型后可以得出，由于模型中小径纤维数量增加了，即模型中有了更多的碳纤维材料，其轴向的力学性能也就相应地获得了提升。因此在进行模型生成时可以根据该模型使用的条件不同来对纤维数量比进行灵活的调整，即当纤维杂化材料需要抵抗拉压变形时可以提高小纤维的比例以提升其轴向的力学性能，而当应用场景要求材料拥有良好的剪切性能时可以增大大径(玻璃纤维)的数量比，来满足使用需求。
3.3. 模型可靠性分析
由于本文所使用的模型是基于全新的算法(随机顺序扩展算法与硬核模型算法相结合)下生成，所以需要对模型的可靠性进行分析。本文采用对不同面积的树脂模型进行分析，选取了窗口面积50 × 50、60 × 60和70 × 70的模型进行比较分析，变量参数总体积分数设定为0.55，可以得到三个模型及其变形前后应力对比图如图8，图9所示。
(a) 窗口面积50 × 50 (b) 窗口面积60 × 60 (c) 窗口面积70 × 70
Figure 8. Different window area models generate graphs
图8. 不同窗口面积模型生成图
(a) 窗口面积50 × 50 (b) 窗口面积60 × 60 (c) 窗口面积70 × 70
Figure 9. Stress comparison diagram before and after deformation
图9. 变形前后应力对比图
与相同总体积分数下不同纤维数量比相同，不同窗口面积需要进行分析比较的数据也是轴向的弹性模量，因此在计算中仅对轴向弹性模量进行计算。表5为各模型轴向弹性模量计算表。
Table 5. Axial elastic modulus table of different window area models
表5. 不同窗口面积模型轴向弹性模量表
|
模型一 |
模型二 |
模型三 |
窗口面积 |
50 × 50 |
60 × 60 |
70 × 70 |
弹性模量E(GPa) |
66.898 |
68.906 |
64.087 |
泊松比υ |
0.2856 |
0.2862 |
0.2899 |
由表5可以得出，即使窗口的面积不同，并不会直接影响到模型的力学性能，虽然窗口面积越小在同样参数的条件下的模型参数值会有偏差，但是泊松比仍然是趋于稳定的，这就意味着模型仍然保持着较好的一致性。因此工程实际中在条件允许的情况下仍然优先选择足够大的窗口面积。
4. 结论
本文主要研究碳纤维/玻璃纤维杂化复合材料不同的材料参数对其力学性能的影响，主要有如下几个结论：
(1) 建立了一种新的模型算法，即将随机顺序扩展算法与硬核模型算法相结合的方法，该算法不仅满足了可靠性的要求，同时还提供了良好的随机性属性，可以应用于现实中对杂化纤维材料进行模拟。
(2) 在材料大小径纤维数量比确定情况下，材料的轴向弹性模量大小和纤维的总体积分数呈正相关；在材料总体积分数确定情况下，提高小径纤维数量可以提高轴向的力学性能，提高大径纤维数量可以提高材料的剪切性能。
NOTES
*通讯作者。