随机多尺度短切碳纤维复合结构模型中RVE 尺寸效应和方向模量的均一化响应

邱伊健, 郑萍, 程香平, 郭玉松

(1. 江西省科学院 应用物理研究所，江西 南昌 330001; 2. 南昌师范学院，江西 南昌 330032)

0 引言

手工铺设、长丝缠绕、蒸压和真空袋处理技 术[1]，通常与连续纤维增强热固性复合材料相关，它通常应用于短期或一次性使用的产品，并且对产品质量和性能的要求较高。然而，当需要大量的此类产品时，这些技术是有限制的。对于当前许多应用来说，强度或刚度峰值的提高并不是这种复合材料的主要需求。而是在满足特定强度和刚度需求情况下，人们往往希望批量生产这种产品，使得用于连续纤维的加工技术在机器时间和人工时间方面显得极其冗长和昂贵。通常聚合物基复合材料中的基体不具有承受很大载荷的能力，其主要作用之一是将负载从一根纤维转移到另一根纤维。它影响复合材料的抗压、抗剪强度，也有助于提高复合材料的抗断裂性能和能量吸收能力。

短纤维增强树脂基复合材料(SFERC)[2-4]由于其良好的经济性和相对可靠的机械性能已经被广泛关注。近年来，由于复合材料制造工艺的进步，短切纤维增强弹性复合材料在中低端工程结构中的力学性能能够满足大多数使用需求，在某种程度上能够替代连续纤维增强复合材料(CFRC)[5-6]。由于短切纤维在基体中复合状态的无序性，并呈现卷曲 状[7]，使得这种材料在细观层面存在极其复杂的结构。如何将SFERC 中复杂的细观结构建立准确的均一化方法，是一项具有挑战性的工作。

已有大量研究工作提出SFERC 的均质化方 法[8-13]，传统的方法是依靠实验技术来表征材料的有效力学性能。这些性能用适当的现象学本构模型来表达时，可以作为传统工程材料力学分析的基础。然而，现代方法论倾向于用均一化方案和微观力学模型结合来取代这种现象学的方法，它是通过在宏观和微观两个尺度下同时求解，以降低求解非线性和时变响应的两个尺度相关计算代价的解耦方法。复合材料均一化理论最早由Hill[9]和Ogden[10]提出。但是由于复合材料均一化场是复杂结构的非线性场，在受到外部载荷时场量的非线性变化限制了其实际应用。随后Castaneda[11-12]提出了2 阶均一化方案并由Kouznetsova 等[13]成熟应用。由于平均场理论组分均一化响应的限制，Chen 等[14]采用全场法得到了复合结构的局部应力和应变场响应。全场法比平均场法具有更高的微观场分辨率。最近Caylak 等[15]研究了结合平均场和全场均一化方法的框架下建立多态不确定模型，用模糊随机变量发展多态不确定平均场和全场方法，在纤维体积分数不确定几何形状和材料参数不确定的情况下，确定纤维增强复合材料的随机有效性能。

除此之外，人们还将各种其他领域的研究方法引入复合材料的细观力学研究[16-20]，例如机器学习和人工智能记忆法[18,21-24]，即便它们目前与传统的直接数值计算法相比并没有带来切实的计算效益，但是可以作为未来本研究领域的一个新方向。

与直立纤维相比，卷曲纤维的随机离散化计算很难采用直接的有限元分析(FEA)方法来处理。直立纤维可以直接离散到代表体积元(RVE)结构中，而卷曲纤维由于其模量在各部位是不一致的，需要先将卷曲纤维的模量进行局部均一化，然后对均一化完成后的局部样本进行离散。传统均一化方法通常采用一个尺度，即RVE 尺度，直接把短切纤维离散到RVE 结构中，然后对RVE 进行均一化。这种离散算法最典型的是2 阶或4 阶随机纤维取向张量的统计离散法。而本文要解决的问题在于卷曲纤维模量是随着局部坐标位置不断变化的，它不能像直立纤维那样直接通过3 个方向的正向模量和3 个方向的剪切模量来表征。

本文对于卷曲纤维的方向性模量需要在纤维局部进行模量路径上微分转置处理，将偏轴模量转置到正轴模量的分量上，然后对卷曲纤维进行离散。与传统算法不同的是，本文的创新点在于在有限元模型的构建中采用双尺度，即单胞(UC)和RVE 尺度同时进行分析。首先将卷曲纤维的模量均一化到UC 尺度，获取UC 结构的均一化模量，然后采用本文提出的算法将UC 尺度模型离散到RVE 结构中，进行二次均一化。为了使得双尺度均一化的实现，本文采用一种新型的随机场，称为空间方向敏感(HE)算法[25]，将空间随机场分为两部分，其中一部分为确定场量，另一部分为随机场量。

1 多尺度均一化实验

1.1 卷曲短切纤维弹性场微分的坐标变换

为描述单根卷曲短切纤维结构模量，把短切纤维假设为 Hsiao-Daniel 正弦波形式的卷曲状 态[8]。其卷曲周期为ω，振幅为A，如图1 所示。图1 中，θ为轴与Ox轴的旋转角度，轴表示局部纤维长度方向，轴、轴分别表示与长度方向垂直的纤维截面方向。

图1 卷曲短纤维长度微分的坐标变换Fig. 1 Coordinate transformation of length differential of CCF

纤维的波动性方程在Oxyz坐标系(后文简称主坐标系)下可以定义为

式中：u为纤维波动的纵坐标值。

在卷曲纤维的任意一端截取一段长度为dx的波动段，并且在平行于这段纤维主方向的切线上建立转置纤维局部坐标系(后文简称转置坐标系)。其中坐标轴平行于切线方向，坐标轴垂直于切线方向，而坐标轴与Oyz轴保持一致。因此主坐标系与转置坐标系下的应力转换可以表示为

式中：σxx、σyy、σzz分别为主坐标系3 个方向的正应力；τyz、τxz、τxy分别为主坐标系3 个方向的剪应力；Tij为转置矩阵(i、j为矩阵分量的下角标)；分别为转置标系的3 个正应力；分别为转置坐标系的3 个方向剪应力。

这里的Tij-1为矩阵Tij的逆矩阵，并且

同理，相应的应变转换方程可以表示为

式中：εx、εy、εz和分别为主坐标系和转置坐标系3 个方向的正应变；γx、γy、γz和分别为主坐标系和转置坐标系下 3 个方向的剪应变。

通过在式(4)两边分别左乘对角矩阵Rij，

将左边正应变和剪应变化简成一致形式，右边左乘一个单位矩阵E=Rij-1Rij，将右边也简化成一致形式，得到

因此主坐标系下的应力应变关系可以表示为

式中：S和Sij分别为主坐标系和转置坐标系下的柔量矩阵，二者存在如下转换关系：

在转置坐标系下，纤维可以看成主轴方向平行于Ofxf轴的横观各向异性材料，对于横观各向异性的柔量矩阵采用Hill 矩阵可以写成如下形式：

1.2 短切纤维的随机场

1.2.1 Digimat 原生算法

Digimat 软件采用2 阶张量取向填充算法，它是一种以概率随机取向形式统计的一种纤维离散方法，这种算法一般用于直立纤维，并且纤维结构中不附带基体相的情况。对于短纤维复合材料，定义微观结构的最重要的随机变量是纤维体积分数vf、纤维取向和纤维长度lf。由图2(a)可知，可以根据两个纤维角φf和θf来定义单根纤维取向。

图2 纤维取向Fig. 2 Fiber orientation

图2 中，Ω为复合材料的结构域，P(q)为 第q根纤维的方向矢量，x1′、x2′、x3′分别为纤维方向矢量分解的3 个正交向量，e1′、e2′、e3′分别为纤维方向矢量分解的3 个正交基向量。

式中：P1f、P2f、P3f分别为纤维取向的3 个正交方向矢量分量；α1f、α2f、α3f分别为纤维与总体坐标系所呈的3 个Euler 角。

需要说明的是，使用欧拉角或方向单位向量时，3 个变量中只有2 个是独立的。所有4 个独立的随机变量Pf、φf、θf和lf(或者分别是欧拉角或方向矢量分量中的2 个)需要提供适当的概率分布，来表示纤维体积分数、取向和长度的不确定性。

对于2 阶取向张量填充算法(2 阶纤维张量取向描述法)来填充聚集无序短纤维的取向问题，因为它们可以与纤维取向的变换矩阵和概率分布直接相关。这里简要介绍一下张量取向aij的定义的推导，首先考虑包含根据图2(b)中有限数量的纤维的任意参考体积Ω，其中包含了有限数目nf根短纤维，如图2(b)所示。任意纤维号为q的取向是根据其单独的纤维取向向量P(q)定义的，其向量组成P(q)根据式(9)来定义。因此，对于参考体积为Ω，纤维总数为nf的2 阶纤维张量取向由式(10)定义：

式(10)用于域Ω内的所有纤维。如果Ω内的纤维数目接近无穷大，即nf→∞，则用连续定义 式(11)来代替2 阶纤维张量取向的离散定义(见 式(10))：

式中：f(Pf)为纤维取向Pf出现的概率密度分布函数。考虑到概率分布中期望值的定义：

根据式(12)，2 阶纤维张量取向的分量aij定义了纤维(归一化)方向矢量分量乘积的期望值，因此参考式(9)，即乘积的期望值，分别为离散和连续纤维取向分布的欧拉角。由于式(10)和式(11)中乘积的交换性，2 阶张量取向是对称的，即aij=aji。因此，对于完全排列的纤维，全部朝向x1方向，2 阶纤维张量取向的形式为

而

则表示纤维取向完全无序。

1.2.2 HE 算法

在获取短切纤维卷曲弹性应力应变场变换之后，还需要考虑卷曲纤维在基体相中的离散分布。因此需要在体积元中建立合适的纤维随机场来描述离散纤维强化基体复合结构的弹性本构。本文通过RVE 结构来研究离散短纤维在基体中的离散行为，其坐标系为Orxryrzr，而单根短切纤维通过UC结构来表征，其坐标系用Ouxuyuzu表示，它内部包含了单根短切纤维，与附着于纤维壁呈一定厚度的界面层，以及部分UC 结构中包裹纤维及其界面层的基体材料(见图3)。

图3 单胞纤维在RVE 中的离散角度Fig. 3 Discrete angle of single cell fiber in RVE

图3 中，α为UC 坐标系旋转矢量与主坐标系Ox轴的夹角，β为旋转矢量在主坐标系Oyz平面上的投影与Oy轴的夹角。

本文中介绍一种新型随机场模型，这种模型称为HE 算法，本文将该算法应用于针对卷曲短纤维随机场的离散化模型。HE 算法首先将一个空间随机场H(x，y，z，δ)分解为一个确定域和一个随机域：

式中：为空间场平均值(它又可以分为α和β角度方向的空间场(见图3)；δ为初始随机性；ξ为常变量；λi和φi分别为解析协方差内核的本征值和本征函数；n为HE 的截断阶。通过求解下面的Fredholm 积分方程[26-30]解出λi和φi的值：

式中：x1(x1,y1,z1)和x2(x2,y2,z2)分别为Ofxfyfzf坐标系内随机选择的两个点；VRVE为RVE 的体积；C(x1,x2)为假设协方差函数的指数形式，

UC 纤维的离散通过两个旋转张量角度α和β通过旋转变换将场量由主坐标系旋转至离散UC 坐标系，如图3 所示。本文采用Galerkin 有限元法[31]对特征函数进行离散。

1.3 短切卷曲纤维强化树脂基三维随机排列的概率模型

在卷曲纤维外层一定厚度范围内存在一层界面相，其属性是树脂基体材料，并且界面相承担了更密集的应力传递，它比远离纤维的基体形变量更大，因此在直径为dCCF的纤维外部增加一层外径为di的界面层。考虑含基体相的独立短切纤维复合结构UC，如图4所示。图4 中，LCCF为CCF的长度。

图4 包含单根波动碳纤维的单胞结构 Fig. 4 Unit cell with single wavy CCF

CCF 的随机卷曲度[6]可以定义为

式中：x为微尺度域上的位置向量；LCCF为CCF 的长度。CCF UC 的无限小片段看成是离轴横断各向同性层(见图4)。通过把Chamis 的简化微观力学模型[24]应用于短切纤维随机场，可以得到5 个独立的弹性属性E11、E22、G12、G23、V12。

对UC 的E1和υ12使用混合等效应变准则时[30]，可以得到

式中：VCCF(x,δ)为碳纤维的体积；ECCF1为碳纤维长度方向模量；Vm为基体体积；Em为基体弹性模量；Vi(x,δ)为随机界面的体积；Ei(x,δ)为随机界面弹性模量；v12(x,δ)为基体泊松比；vCCF12为碳纤维泊松比；νi为随机界面的泊松比。

界面材料的体积分数可以表示为

式中：di包含了纤维和界面层；dCCF为不包含界面层的单根纤维外经。横向有效模量可以通过复合材料基础力学推导而来[32-37]：

同理，随机剪切模量(G23(x，δ)，G12(x，δ)=G13(x，δ))也可以通过将式(19)中的Em替换为Gm，ECCF2替换为GCCF12得到。随后的计算可以通过第2 节所述短切纤维卷曲度坐标变换公式计算。同时，短切碳纤维增强相的长径比影响可以间接通过局部体积分数和纤维卷曲度来考虑。

本文提出的构建空间随机CCF 排列模型的方法可以通过两个取向角来表征(αCCF和βCCF)，分别表示UC 在RVE 中的α和β的离散角。这里将CCF 看成是随空间变化的随机场，并且通过HE 方法将其表示为

式中： ＜·＞表示平均值。x轴、y轴、z轴方向的角度均值和相关长度是决定CCF 随机取向的主要统计参数因素。在实际的聚合物基体中，每根CCF不可能完全对齐，如图5 所示，图中〈θCCF〉、θird分别为平均波动和随机波动。因此，取向的CCF复合材料具有空间变化的各向异性。

图5 RVE 空间变化的随机排列角度Fig. 5 Random arrangement angle of RVE space variation

应力张量和应变张量在RVE 坐标系Orxryrzr和UC 坐标系Ouxuyuzu之间来回变换，如式(21)所示：

随机转动张量根据式(20)中的角度可以表 示为

式中：Rkl表示随机旋转张量变换。

1.4 平均应力，应变及其后处理方法

在三维RVE 结构中，每个离散的UC 可以看成是包含在其中的一个子域，如果要从子域中提取应力和应变，则在对结果进行后处理时不能使用关键点自由度(Kdofs)，因为 Kdofs 是针对完整的RVE，而不是子域。应力和应变必须以某种方式平均，涉及到应力和应变在子域上的积分。由于将应力和应变的积分过程作为计算机代码来实现对于大多数用户是一项要求很高的工作，下面的简化可以保持数学上的严谨性。高斯定理(散度定理)给出了下列积分公式：

式中：[fx,fy,fz]为三维空间中的连续可微向量场；∂Ω为三维空间中感兴趣域的边界；[nx,ny,nz]为∂Ω的单位向外法线；S为∂Ω的表面积。特别地，如果在三维空间中引入[fx,fy,fz]=[gxhx,gyhy,gzhz]，这里g和h是两个单独的连续可微向量场，类似于f，高斯定理给出了下列分段积分：

利用这些数学公式，可以导出三维情况下的平均应变和应力，为此使用高斯定理得出以下平均应变的表达式：

1.5 数值模拟算法实现和材料参数

1.5.1 算法实现

本文提出的算法实现过程如图6 所示。首先对DIGIMAT 脚本进行RVE 参数，材料模型等相关数据输入，包括HE 算法、RVE 设计尺寸、边界条件、纤维目标体积分数、纤维半径、纵横比、纤维填充数目、纤维最小间距、界面相直径与纤维直径之比(di/dCCF)等参数。然后程序根据HE 算法的周期边界条件对RVE 进行填充，得到相应的RVE FE 结构模型。在MATLAB 软件脚本对RVE FE 结构模型进行Halpin-Tsai 方程坐标变换，将宏观RVE 尺度随机场量离散到细观的UC 尺度下，并得到以UC 为最小单位的离散化的空间随机分布模量场，将离散化数据在Abaqus CAE 软件中基于空间随机分布模量场进行相应的有限元计算，从而得到相应的数值解。

图6 HE 算法的流程图Fig. 6 Flowchart of space direction-sensitive hoteling expansion algorithm

1.5.2 材料属性和实验材料制备

表1 为RVE 中3 种材料结构的材料参数和统计参数。需要说明的是，基体材料和界面材料是同种物质，但是界面材料由于紧邻CCF 相，受到其强化作用的影响，其材料力学行为与基体材料在微力学模型中存在一些区别，但是这种细微的调整在微力学模型中是必要的，也是可以接受的。虽然CCF 填充树脂基复合材料复合强度随着CF 含量的增加递增，但是短纤维强化与长纤维不同，短纤维增强树脂基体随着纤维含量的增加，复合结构脆性断裂的风险呈指数递增。因此为兼顾韧性和强度，短纤维增强环氧树脂基复合材料的纤维体积分数应该低于3%。本文选取最大体积分数为1.24%。

基体材料考虑塑性变形，采用J2 塑性模型[30]，并且考虑的是指数和线性强化准则，CCF 只考虑弹性效应，相应参数如表1 所示。

表1 RVE 结构各相材料属性Table 1 Material properties of each phase in RVE structure

RVE 结构微力学实验模型如图7 所示。选取平均长度为1 mm 的短切碳纤维作为增强相，均匀分散于环氧树脂基体中，并用聚四氟乙烯模具将其定型固化为正方体形状，如图7(a)所示。力学实验选取尺寸为20 mm×20 mm×20 mm 的立方体，填充纤维体积分数分别选取0.2%、0.41%、0.61%、0.82%、1.03%和1.24% 6 组进行实验。其中，图7(b)所示为纤维体积分数为1.24%的样本中短切碳纤维在基体内随机分布的细观显微图像。然后将固化的立方体试样再次置于模具中，加入无增强相的环氧树脂固化后，成为拉伸力学试样，如图7(c)所示。

图7 复合结构实验样本及其细观显微图像Fig. 7 Composite structure test samples and their mesoscopic images

2 实验结果与讨论

2.1 实验样本的纤维长度及2 阶张量取向角的概率分布测定

为确保数值RVE 模型参数与实验样本的一致性，采用样品基体的化学焚化法提取纤维，通过图像处理系统对它们进行长度扫描，统计图7 所示的样本中短纤维的长度概率分布情况。图8(a)所示为测试样本纤维长度的统计数据。图8 中，f(LCCF)为短纤维长度的概率密度分布函数，而F(LCCF)为纤维长度归一化概率分布函数，它们之间的关系为f(LCCF)=dF(LCCF)/dx。

测试样本中，短纤维长度分布为0~4.8 mm，纤维长度分布概率整体呈线性递增，同时存在一些概率跳跃分布点，这些跳跃点在较短纤维区间 (< 1 mm)分布较明显。

同时还确定了纤维取向分布，它需要有关微观结构几何形状的完整的三维信息。如果纤维和基体的密度存在较大不同，例如对于大多数玻璃纤维，碳纤维增强的聚合物，可以通过X 射线计算机断层扫描以无损方式研究微观结构。为此，在大量连续的旋转角度下，通过X 射线源对样本进行射线照相。对拍摄的图像进行计算，可以重建微观结构；可获得基于体素的显微组织三维图像结果并跟踪每根纤维，从而提供必要的纤维方向数据，以根据 2 阶张量取向法确定纤维取向张量。通过X 射线计算机断层扫描获取纤维取向矢量p，并得出样本中每根纤维的纤维角度。图8(b)、图8(c)给出了实验样本中根据图3 所示角度确定的纤维取向分布曲线，该样本取向纤维以仰角α和投影角β确定。其中投影角α和β的概率分布都是呈高斯正态分布，β的概率分布统计的周期为π，即[-90°，90°]，而α的概率分布统计的周期为π/2，即[-45°，45°]，从而导致β和α分别呈现单峰和双峰。

图8 实验样本中随机取向纤维的长度和方向角的统计分布Fig. 8 Statistical distribution of length and orientation angle of randomly oriented fibers in experimental samples

2.2 RVE 尺寸效应及边界影响区长度

2.2.1 RVE 尺寸效应

RVE 通常在下部长度尺内引入，来描述其上部长度尺度的性质，其定义基于代表性这个关键词。如果RVE 能够再现从有关长度尺度的无限大体积材料中获得的特性，则认为RVE 具有代表性，这应该作为确定它的尺寸的一个标准。代表性通常是根据感兴趣材料上部长尺度的有效性质来定义的。只要材料在其上部长度尺度上是均匀的，就可以在其下部长度尺度上引入RVE 作为材料的适当体积。在多尺度模型中，在上部尺度中的任何有限维在下部尺度上都被认为是无限的。换言之，RVE 的维数必须在下部长度尺度上是有限的，而它在其上部长度尺度上可以被视为数学上的无穷小。在较低长度尺度上的无限维体积总是具有代表性，但不是理想的RVE 候选体。

根据Hill[9]的理论，RVE 尺寸对数值模拟的精度影响很大程度上取决于边界条件影响的区域深度。它是指从曲线或曲面到边界效应减小点的距离，是一个特征度量。虽然它沿RVE 的表面逐点变化，但它是有界的，其上限称为衰减长度t，t是给定材料的确定值。受边界条件影响的区域的体积Vs影响。Vs可以等效地看作是衰减长度厚度的壳，是RVE 的最外层。由于对于给定的材料，t通常是一个固定值，通过衰减长度到RVE 的距离，应力和应变场不再受强加边界条件的影响，例如，是否通过施加均匀位移或均匀牵引来完成。在将具有衰减长度厚度的壳体从RVE 上剥离后，留下的芯中的应力和应变场应与通过分析无限大体积材料获得的应力和应变场相同。例如，考虑具有边长a、衰减长度t和体积为VRVE(RVE 的体积)的立方RVE。随着RVE 的大小增加，a增加，而t保持不变。当t相对于a变得可以忽略时，两个体积的比值r趋于消失(见式(31))，此时RVE 预测的有效性无限接近其上部尺度。

目前实验测定衰减长度采用数字散斑全场应变仪测定试样在加载前后的局部应变量的不同来计算衰减长度，而数值方法可以直接通过FEA 计算 获得。

在此模拟中，采用空间方向敏感性Hotelling膨胀算法对CCF 镶嵌环氧树脂的RVE 进行填充。为研究RVE 尺寸对力学精度的影响，考虑的RVE尺 寸 分 别 为 0.25 mm×0.25 mm×0.25 mm 、 0.5 mm×0.5 mm×0.5mm 和1 mm×1 mm×1 mm，并对3 种尺寸的RVE 立方体基于HE 算法进行填充。力学实验则选择图7(b)所示实验样本。界面刚度假设为0 Pa，并且在RVE 中实现了弹性刚度随机离散分布。RVE 结构的网格划分选取的是8 节点砖块单元，单元尺寸均为5 μm×5 μm×5 μm。模型采用的是不考虑孔隙问题的理想化情况。

2.2.2 RVE 边界影响区长度

图9 显示了3 种尺寸的RVE 结构中，采用周期 边界条件HE 填充碳纤维的具体过程，体积分数由0.2%递增至1.24%，图中每列所对应的体积分数是一致的。其中VRVE=13mm 的最大纤维填充数量为 330根，VRVE=0.53mm的最大纤维填充数量为44根，而VRVE=0.253mm的最大纤维填充数量为8根。

图9 3 种尺寸RVE 的纤维填充过程(纤维横纵比：145DCCF=6.9 μm，纤维扭曲度=0.3)Fig. 9 Fiber filling process of CCF for RVE of three sizes(Fiber aspect ratio: 145DCCF=6.9 μm, Torsion of fiber= 0.3)

结构边界条件的约束方法是通过有限元节点施加线性约束方程，并将周期边界条件施加到带有周期性的RVE 结构中。对所属面(见图10)施加约束方程，对于面(也可以是边或者顶点)的每个节点上的任意自由度i，边界约束方程如下：

图10 RVE 面边界约束方程Fig. 10 Surface boundary conditions of RVE

段落 true="1">式中：ui为第i个分量的位移，i=1, 2, 3；Lx、Ly、Lz分别为x轴、y轴、z轴方向的初始长度。

同时，对3 种不同尺寸的RVE 在x′轴方向上施加0.03 mm 的位移边界条件，通过数值计算，获取了对应的边界影响区长度(见图11)。

由于纤维与基体材料存在非常大的平均应力差，图11 中把短纤维和界面材料从基体中去除，从而可以准确获取边界影响区长度的细节信息。从图11 中可以看到：基体材料在此边界条件的影响下，应力分布为20~55 MPa 之间；通过对比不同x轴方向深度的切面应力分布图可知，由于存在边界影响区，RVE 结构由外到内存在不同的应力分布梯度；VCCF=1.24%(短切碳纤维体积分数)的RVE结构，应力影响区长度t取值为15~20 μm 之间，它受到RVE 尺寸的影响并不明显。在边界影响区范围内，外层平均应力值明显大于内层。当t值超过20 μm 时，平均应力值不随x深度变化而改变。

图11 衰减长度t对RVE 计算精度的影响，VCCF=1.24%Fig. 11 Influence of attenuation lengthton the accuracy of RVE calculation,VCCF=1.24%

为比较随机纤维体积分数，UC尺寸(见图3)，以及纤维扭曲度在3种尺寸RVE中对t的影响，对RVE结构模型进行了系统的研究，如图12所示。首先研究不同VCCF下3种尺寸RVE的衰减长度与实验测定的值进行对比(见图11(a))。由于短切纤维的不断填充，RVE结构内两相界面的指数式增加以及相界面无序性的增大，使得应力影响区效应变得明显，整体的t值随Vccf呈指数变大，其对t值的影响范围约为17 μm。可是通过图12(a)对比，在其他条件(各相组成、VCCF、扭曲度等)相同前提下，t值与RVE尺寸无关，而实验测定的t值比数值模拟的结果略大。

此外，填充纤维的弯曲度对t值也存在一定的影响(见图12(b))，扭曲度的影响与t值基本呈线性关系，但是弱于VCCF的影响，其影响范围约为4 μm。图12(c)所示为UC的尺寸对t值的影响。对于UC结构(见图3)，考虑它是由纤维长度LCCF和基体长度Lm共同组成，并且LCCF为确定的值，通常是通过控制Lm来调节UC的尺寸，因此这里采用Lm/LCCF的值(UC中基体长度与纤维长度的比值)来确定UC的尺寸。结果显示，UC中基体相越大，t值越小，其对t值的影响范围约为1.6 μm。

图12 RVE 体积，纤维扭曲度和UC 结构与RVE 应力影响区长度关系Fig. 12 Relationship among RVE volume, fiber tortuosity, UC structure and the length of RVE

对比图3 可知，对RVE 衰减长度t的影响因素有很多，但是影响程度为VCCF>扭曲度 >Lm/LCCF>VRVE。

2.3 HE 算法的模量均一化响应

由于CCF 在基体中的随机分散，复合结构中细观层面的结构模量存在差异性，这主要体现在纤维和基体界面层的模量扰动。过大的扰动也给结构模量均一化过程带来了相应的误差。本文研究6 组方向模量，包括3 个方向杨氏模量、3 个剪切模量和3 个泊松比。

图13 所示为尺寸VRVE=0.5 mm3，VCCF=1.24 以及旋转角度α=45°、β=90° RVE 中的6 个方向模量和3 个泊松比的计算结果。从图13 中可以看出：E1′、E2′的均值在270~330 GPa 之间，由于受到旋转角度的影响，E3′均值模量略小，其范围在230~300 GPa之间；G1′2′、G2′3′和G3′1′均值在500 GPa 以上；泊松比ν1′2′、ν2′3′和ν3′1′的均值在0.22~0.27 之间，并且基体相分布略大于纤维相，与实际情况相符。

图13 VCCF=1.24%以及旋转张量角度α=45°，β=90°时的模量空间随机分布Fig. 13 Spatial random distribution of modulus withVCCF=1.24% and rotation tensor angleα=45° andβ=90°

RVE结构中方向模量的平均响应主要受到旋转角度(见图3)的影响，因此通过固定一个旋转角度(α或β)来研究方向模量受到另一个角度旋转的变化，旋转角度α和β在RVE中的取值范围均为0°~90°之间。图14和图15分别显示了将旋转角度α和β固定(45°)，方向模量随另一个角度的变化情况。结果显示，复合结构中平均剪切模量变化范围为550~600 GPa之间，平均杨氏模量变化范围为280~320 GPa之间，泊松比的平均变化范围为0.231~0.257。

图14 随角β改变的平均方向模量的变化Fig. 14 Average directional modulus varying with rotation angleβ

图15 随角α改变的平均方向模量的变化Fig. 15 Average directional modulus varying with rotation angleα

由图14和图15可见：G1′2′、G2′3′和E2’受角度β的变化较大，G2′3′、E1′和E2′受角度α的影较大；而泊松比均一化结果是随机分布的(见图14(b))，与旋转角度变化没有明显关系；当α=45°时，6个方向模量随β的0-90°变 化 值 的 大 小 关 系 为G2'3'>G1′2′>G3′1′>E2′>E1′>E3'；β=45°时，6个方向模量随α的0°~90°变化值的大小关系为G2′3′>E2′>E1′>G1′2′>G3′1′>E3'；受两个旋转角度影响最大的方向模量均为G2′3′，它随两个角度的变化值分别为50 GPa和72 GPa。表明G2′3′对于旋转角度更敏感，无论是哪个角度的旋转，在2′3′平面上，有更多的随机短纤维排布由纬向转变为径向，或者由径向转变成纬向。

除此之外，为验证本文中方向模量均一化的准确性，与拉伸实验样本(见图7)测定的结果进行了对比，对比结果如表2 所示。

表2 当α=β=45°时，均一化模量与实验对比Table 2 Comparison of the homogenization modulus with the experiment when α=β=45°

根据表2 中所示，实验测定的和有限元模量均一化计算出的结果差异在2~7 GPa 之间，并且实验测定的数值比有限元计算出的数值偏大，可以根据图11(a)中的结果来说明，因为实验测定的应力影响区厚度较数值分析的大，而应力影响区存在较大的集中应力，因此会对模量的均一化评估产生偏大的影响。

3 结论

本文提出一种采用Hotelling expansion 模型进行CCF 的三维随机填充的均一化方法，并对该方法生成的RVE 结构进行了研究。首先通过数字散斑全场应变测量仪测定了CCF 增强环氧树脂基复合结构的衰减长度，并利用数值模拟的方法研究衰减长度的影响因素。然后利用所提出的均一化模型研究了短纤维空间分布引起的方向模量各向异性响应行为。根据数值实验结果，讨论了旋转角度对短纤维的空间随机分布对复合材料方向模量的影响。得出主要结论如下：

1)提出一种有效预测随机分布短纤维复合材料方向模量均一化的解决方案。在研究复合材料均一化方向模量时，由于不同方向模量与旋转角度之间存在的敏感度各异，导致各方向模量均一化结果随旋转角度变化存在大小差异。

2)短纤维在空间上的随机分布状态的不同导致了RVE 衰减长度尺寸的变化。受实验所需成本和时间的限制，仅对RVE 体积VRVE、纤维扭曲度，UC 结构与衰减长度t的测试，并得出其相应的变化情况。数值计算结果对比实验数据的结果较小主要是由于除去内嵌的纤维，真实的材料中存在气泡，夹杂物等微小孔洞，而数值分析并没有将其考虑进去。这一干扰效应的屏蔽及其对有效数值计算精度的修正还有待于进一步的研究。