短切玻璃纤维复合材料有限元模拟及其随机分布快速生成算法

Finite Element Simulation of Short-chopped Glass Fiber Composites and a Fast Generation Algorithm for Their Random Distribution

摘要

复合材料的剪切、横向力学性能以及耐久性在很大程度上取决于纤维与基体的界面强度。复合材料界面强度通常依赖纤维和树脂基体的大小、形状、性质和空间分布。文章提出一种适用于复合材料的数值模拟计算方法，选择短切玻璃纤维复合材料为研究对象进行有限元模拟计算，对纤维的随机分布提供一种快速生成算法，分析短切玻璃纤维长度对片状模塑料力学性能的影响。结果表明，短切玻璃纤维片状模塑料弹性模量随纤维长度的增大出现先增大后变缓的趋势，通过比较随机结构和规则结构说明计算方法合理性，纤维分布的随机性对复合材料宏观弹性常数的影响较小。

关键词

短切纤维增强 / 随机分布模型 / 数值模拟 / 纤维长度 / 弹性常数

Key words

Short-chopped fiber reinforcement / Random distribution model / Numerical simulation / Fiber length / Elastic constant

现代工业的快速发展对材料的物理、化学等各方面性能的要求越来越高。传统单一特性材料已经很难满足一些特定场合对材料使用性能的要求。因此，将性能互补的单一材料以各种方式组合形成的复合材料应运而生。复合材料设计灵活，协同效应明显，可以实现单一组分材料无法达到的性能。纤维增强树脂基复合材料是研发力度投入大、应用面广、市场前景广阔的一种复合材料^[1-2]。纤维增强树脂基复合材料凭借其优良的力学性能、较低的密度以及良好的设计性而被应用于化工、航空航天等高端领域^[3]。随着各种类型的树脂以及增强纤维的相继出现及大规模量产，纤维增强树脂基复合材料的原材料价格大幅度下降，因而得到快速发展，应用领域也不断扩大，在生产生活领域占据重要地位^[4-5]。短纤维增强复合材料是纤维增强复合材料中发展最早的一类复合材料，其力学性能远胜传统材料，成型工艺简单，易加工为复杂的结构件。但是，短纤维增强复合材料中增强纤维的含量及其取向等工艺参数对复合材料的形貌结构和力学性能影响巨大，玻璃纤维的横纵数量比及其含量对最终制品也有较大影响^[6]。目前，绝大多数情况下采用实验方法测试短纤维增强复合材料的力学性能，需要大量的时间和成本，因此建立有限元分析模型尤为重要，可以为短纤维增强复合材料用于力学性能要求较高的结构件、支撑件^[7-8]的生产应用提供借鉴。

复合材料界面强度的表征方法众多，如何选择适当的方法来准确反映复合材料的界面性能一直是复合材料力学领域的热门话题。FUKUDA等^[9]采用概率计算的方法研究纤维长度和取向分布对短纤维复合材料强度的影响。ZHU等^[10]考虑纤维和基体位移形变的不均匀性、残余应力、纤维分散强化和基体中错位密度等因素，提出一种分析随机取向短纤维增强复合材料强度的模型，用于预测复合材料强度。MORTAZAVIAN等^[11]研究了纤维的长度、半径和含量对短纤维增强复合材料断裂行为的影响。KOUZNETSOVA等^[12]通过实验及仿真预测研究两种短玻璃纤维增强热塑性塑料的各向异性对拉伸性能的影响。在不同纤维长度、直径和取向分布的数学模型也与分析方法一起用于预测组成材料的应力应变强度和弹性模量。CHEN等^[13]开发一个使用三维代表性体积单元(3D RVE)对连续纤维增强复合材料进行计算均匀化的框架，然而该分析是基于单纤维方阵RVE。MELRO等^[14-15]应用他们所描述的“体积均质化”来评估不同几何参数对单向复合材料弹性性能的影响，基于3D RVE，在选定的RVE窗口内随机排列纤维增强物。考虑的一些几何参数是代表性体积元素的尺寸、纤维直径和相邻纤维之间的间隔。这种利用纤维随机分布的3D RVE预测有效性能的方法必须与周期性边界条件系统相结合，以确定准确的有效性能。OKEREKE等^[16]研究发现，在许多情况下三维球形RVE比三维立方体RVE更能收敛于预期的有效性能。XIA等^[17]使用他们所描述的菱形或平行六面体三维RVE来研究单向(UD)复合材料的响应。其统一周期性边界条件，以预测单向、交叉层和角度层合板的性能。但该预测是基于纤维内含物的有序排列。截至目前，对纤维增强复合材料的强度理论预测仍然处于探索阶段^[18-19]，建立分析临界区域内纤维取向、长度和体积分数等因素的短纤维增强复合材料力学强度模型尤为重要。

本文利用Python语言编程建模，结合平均力学响应计算方法，模拟对短切玻璃纤维复合材料的力学性能进行研究，为复合材料在工业化上的应用提供支撑。

1 实验部分

1.1 材料成分

通常复合材料主要由树脂基体、增强纤维以及其他添加物组成，树脂基体作为黏合剂将纤维黏合在一起。本实验使用的材料为玻璃纤维增强塑料(GFRP)^[20]，玻璃纤维是复合材料的增强体，成型后的零件强度受玻璃纤维影响较大。GFRP为美国欧文斯科宁公司生产的典型的热固性纤维增强复合材料DERAKANE™470/HT-400。表1为GFRP材料组成。

1.2 仪器与设备

平板硫化仪，TZ112，0~60 MPa，上海天质实业有限公司；万能测试机，TCS-2000，20 000 N，台湾高铁检测仪器。

1.3 性能测试与表征

力学性能测试：为便于得到实验所需力学性能参数及保证实验数据准确性，利用平板硫化机压制GFRP板材，按GB/T 1447—2005^[21]进行测试。

1.4 有限元建模基本假设

通过有限元软件分析复合材料，首先必须设计合适的RVE^[22-23]。对于本研究中的RVE模型，假设纤维增强材料与基体之间完美结合，因此不需要指定纤维与基体之间的相间区。为了在定义的RVE窗口内创建纤维的空间随机分布，采用Python开发脚本算法，在ABAQUS内生成随机定位的复合材料内含物。假定材料结构在宏观上是均匀的，但在微观上是不均匀的(形貌由可区分的成分组成，如晶粒、基体、纤维)。因此，本文建立短切玻璃纤维随机分布应力应变模型，该模型包含其几何分布和数值模型。建立简化模型进行应力应变分析原因如下：一方面，标准力学测试中复合材料通常是薄板结构，其长、宽超过厚度的两个数量级，且在宏观上具有对称性，因此可以更有效地简化模型，减少计算量，提高平均刚度计算准确性，在宏观上模型的纤维分布随机性更为明显。根据这种原理，目前国内很多研究人员采用了仿真模型。李凌岩等^[24]在研究基于3D打印短纤维分布对单轴拉伸载荷作用的影响时建立了三维模型，设定颗粒为矩形，得到与实验对应的规律。江真^[25]通过ABAQUS软件，生成了短切碳纤维在乙烯基树脂中随机分布的应力应变模型。另一方面，三维模型确实能够在直观上更加贴近实际情况，但三维实体模型想要完全复制材料的实际情况也不现实^[26-27]。对于微米级别的纤维，网格划分的精确性严重影响计算结果的准确性，而模型计算量将会非常大，难以实现。因此，本文用简化模型可以在不影响模拟精度的同时有效缩短分析时间。

1.5 纤维随机分布的数值实现

不同于一般的均质材料，复合材料在微观上并不均匀，但纤维在树脂内的分布具有一定的规律性，即呈现统计学上的周期性随机分布，由RVE扩展排布即可构成复合材料整体，当外力场作用时RVE表现相似的应力应变场，因此可以用RVE表征整体材料的受力情况。通过对RVE力学响应进行计算可以获得整体材料的弹性模量。为了实现纤维随机分布，本文采用通过Python语言编程，建立纤维随机分布模型，该方法可通过给定几何信息(包括基体大小、纤维数量、纤维长径比分布和体积分数)来建立随机分布的纤维增强复合材料细观结构模型，并将细观结构模型数据存储在xml数据文件中。该方法可对生成的2D细观模型中的纤维数量、体积分数和长径比进行统计分析，利用ABAQUS软件的脚本语言Python读取xml数据文件，在ABAQUS完成几何模型的建立。首先，建立模型，创建基体平面，在基体平面中创建第一个增强纤维Fiber1。接下来调用Random函数，按任意角度对Fiber1进行旋转平移得到Fiber2，把超出边界的部分放到对边基体中。随后进入判断语句，当没有达到所需的纤维数量N(体积分数V)时，继续循环生成纤维。以此类推，得到每一根纤维，当纤维的根数达到所需的体积含量时输出模型。图1为随机分布实现流程。

随后将编写的Python程序导入ABAQUS中生成纤维随机分布模型。利用ABAQUS软件对模型进行网格划分，图2为短切随机分布增强纤维复合材料网格图。

1.6 周期性边界条件及计算方法

计算采用均质化方法对微观结构的全局周期性进行假设，表明整个宏观试样由空间重复的单位网格组成。即通过平均微观结构中属于该点的小代表性部分的应力和应变来估计某点P周围的宏观应力和应变(如图3所示)。通过将局部宏观应变张量应用于点P的微观结构RVE，并对所得的非均匀RVE应力场进行平均，可以建立若干参考点P点的宏观应力±应变关系，从而通过等效RVE总刚度推导出宏观点处的一致刚度。这样，RVE模型上的问题就可以表述为准静态连续体分析的标准问题，从而计算出整个模型的力学响应。

1.6.1 平衡方程

RVE处于平衡状态。在没有力作用的情况下，RVE域 {\Omega }_{\mathrm{R}\mathrm{V}\mathrm{E}}中柯西应力张量 \mathit{\sigma }应该满足：

\nabla \cdot \sigma = \mathrm{0} (1)

式(1)中： \nabla \cdot \sigma表示柯西应力张量 \mathit{\sigma }的离散度。

1.6.2 本构行为

微观结构成分的力学特性由本构定律描述，该本构定律指定依赖时间和历程的应力-应变关系如下：

{\sigma }_k \left(t\right) = F_k \left\{\mathit{F}\left(\tau \right),\tau \in \left[\mathrm{0},t\right]\right\} , k = \mathrm{1,2} , \cdots , N (2)

式(2)中：N为具有不同材料行为的微观结构成分(例如基体matrix和纤维fiber)的数量；t为进程时间； \mathit{F} \left(\tau \right) = {\left(\nabla \mathrm{0}y\right)}^{\mathrm{c}}， \mathit{F} \left(\tau \right)为任意点在t时刻的变形梯度张量，其实际位置向量为 \mathit{y} \left(\tau \right)， \nabla \mathrm{0}为梯度算子对参考构型的趋向，上标c表示参数共轭； {\sigma }_k \left(t\right)为应变；和 F_k为应力，MPa。

1.6.3 边界条件

对RVE施加周期边界条件。周期性边界条件意味着：(1)在变形过程中，两个相对边缘的形状和空间方向是相同的，并且保持相同。(2)为了在边界上具有应力连续性，作用在相对边上的应力矢量方向相反。

当图1中的边界 {\mathit{\Gamma }}_{\mathrm{43}}和 {\mathit{\Gamma }}_{\mathrm{23}}与其他(保留或独立的)边界 {\mathit{\Gamma }}_{\mathrm{12}}和 {\mathit{\Gamma }}_{\mathrm{14}}以及顶点1、2、4相连接时， {\Omega }_{\mathrm{R}\mathrm{V}\mathrm{E}}的空间周期性得到加强，并且定义了运动学边界约束，其表达式为：

{\mathit{y}}_{\mathrm{43}} \left(s_{\mathrm{43}}\right) = {\mathit{y}}_{\mathrm{12}} \left(s_{\mathrm{12}}\right) - {\mathit{y}}_{\mathrm{1}} + {\mathit{y}}_{\mathrm{4}}(式中s_{\mathrm{43}} = s_{\mathrm{12}})(3)

{\mathit{y}}_{\mathrm{23}} \left(s_{\mathrm{23}}\right) = {\mathit{y}}_{\mathrm{14}} \left(s_{\mathrm{14}}\right) - {\mathit{y}}_{\mathrm{1}} + {\mathit{y}}_{\mathrm{2}}(式中s_{\mathrm{23}} = s_{\mathrm{14}})(4)

式(3)、式(4)中： {\mathit{y}}_i为参考点 i相对某固定参照系的位置向量； {\mathit{y}}_{ij}和 s_{ij}分别为边界上的点的位置向量和局部坐标，其次边界上相反应力向量的要求可以写成 {\mathit{\Gamma }}_{ij}。

相反边界上相反应力向量的表达式如下：

\mathit{\sigma } \cdot n_{\mathrm{43}} \left(s_{\mathrm{43}}\right) = - \mathit{\sigma } \cdot n_{\mathrm{12}} \left(s_{\mathrm{12}}\right)(式中s_{\mathrm{43}} = s_{\mathrm{12}})(5)

\mathit{\sigma } \cdot n_{\mathrm{23}} \left(s_{\mathrm{23}}\right) = - \mathit{\sigma } \cdot n_{\mathrm{14}} \left(s_{\mathrm{14}}\right)(式中s_{\mathrm{23}} = s_{\mathrm{14}})(6)

式(5)、式(6)中： n_{ij}是沿 {\mathit{\Gamma }}_{ij}的单位向外法线。

{\mathit{u}}_i = {\mathit{u}}_i^{\mathrm{0}} , i = \mathrm{1,2} , \mathrm{4} (7)

式(7)中： {\mathit{u}}_i为初始变形梯度张量； {\mathit{u}}_i^{\mathrm{0}}为参考点初始位置。参考点1、2、4的位移是通过对这些顶点的初始位置应用局部宏观变形梯度张量来规定的，这实现了宏观变形场向微观变形场的实际过渡。

1.6.4 形变

假设 F _macro为某宏观点的变形梯度张量，则根据如下公式规定周期RVE点的新位置，可以将该变形梯度张量直接应用于与该点相关RVE模型中。

{\mathit{y}}_i = {\mathit{F}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}} · {\mathit{y}}_{\mathrm{0}i}^{} , i = \mathrm{1,2} , \mathrm{4} (8)

式(8)中： {\mathit{y}}_i为顶点的初始位置，且通过上述公式可以证明RVE模型的点位移等效于 {\mathit{F}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}} = {\overline{\mathit{F}}}_{\mathrm{R}\mathrm{V}\mathrm{E}}。其中 {\overline{\mathit{F}}}_{\mathrm{R}\mathrm{V}\mathrm{E}}为体积平均时效变形梯度张量，这是由RVE变形梯度张量场 F ( y ₀)对RVE初始体积V ₀的直接积分得出的，并按该体积缩放，如式(9)所示。

{\overline{\mathit{F}}}_{\mathrm{R}\mathrm{V}\mathrm{E}} = \frac{\mathrm{1}}{V_{\mathrm{0}}} \underset{{\mathit{y}}_{\mathrm{0}}\in v_{\mathrm{0}}}{\int } \mathit{F}({\mathit{y}}_{\mathrm{0}})\mathrm{d}{V}_{\mathrm{0}} (9)

1.6.5 应力

平均RVE应力张量可以根据RVE分析来定义。由RVE平均应力张量 {\mathit{\sigma }}_{\mathrm{R}\mathrm{V}\mathrm{E}}对RVE模型的任意做功位移 \delta \mathit{u}所做的内部功 \delta W_{\mathrm{i}\mathrm{n}\mathrm{t}}的变化必须等于外部载荷对RVE所做的功 \delta W_{\mathrm{e}\mathrm{x}\mathrm{t}}。由RVE平均应力完成的内部功可以写成式(10)：

\delta W_{\mathrm{i}\mathrm{n}\mathrm{t}} = \underset{V}{\int } {\overline{\mathit{\sigma }}}_{\mathrm{R}\mathrm{V}\mathrm{E}} : \nabla \delta \mathit{u} \mathrm{d} V = {\overline{\mathit{\sigma }}}_{\mathrm{R}\mathrm{V}\mathrm{E}} : \underset{\Gamma }{\int } n\delta \mathit{u}d\mathit{\Gamma } (10)

式(10)中：V为RVE的当前体积；n为变形后的RVE表面 \mathit{\Gamma }上的单位向外法向量；对当前构型取梯度算子 \nabla。结合周期性条件式(3)和式(4)，将内功关系式(10)转化为如下方程：

\delta W_{\mathrm{i}\mathrm{n}\mathrm{t}} = {\overline{\mathbf{\sigma }}}_{\mathrm{R}\mathrm{V}\mathrm{E}} : \left(n_{\mathrm{12}}^{\mathrm{*}}\left(\delta {\mathit{u}}_{\mathrm{1}}-\delta {\mathit{u}}_{\mathrm{2}}\right)‖{\mathit{y}}_{\mathrm{1}}-{\mathit{y}}_{\mathrm{4}}‖+n_{\mathrm{12}}^{\mathrm{*}}\left(\delta {\mathit{u}}_{\mathrm{1}}-\delta {\mathit{u}}_{\mathrm{4}}\right)‖{\mathit{y}}_{\mathrm{1}}-{\mathit{y}}_{\mathrm{2}}‖\right) (11)

式(11)中： n_{ij}^{\mathrm{*}}为当前构型中连接顶点 i和 j的直线的单位向外法线； ‖m‖为任意向量的长度 m。

通过约束条件来执行应用的边界条件，不仅意味着约束点的位置由保留点规定，而且还意味着约束点上的载荷必须转移到保留点上。因此，外部负载的虚拟功仅在保留的顶点和受限RVE边界的保留部分上执行(顶点1、2、4以及边界条件 {\mathit{\Gamma }}_{\mathrm{12}}， {\mathit{\Gamma }}_{\mathrm{14}}如图1所示)。由于应力在由式(5)和式(6)表示的RVE边界上的连续性，边界条件保留部分上分布的外部载荷 P^{{}^{\mathrm{e}\mathrm{x}\mathrm{t}}}也从公式中抵消，如下式所示：

P_{\mathrm{12}}^{\mathrm{e}\mathrm{x}\mathrm{t}} \left(S_{\mathrm{12}}\right) = \mathrm{0} , P_{\mathrm{14}}^{\mathrm{e}\mathrm{x}\mathrm{t}} \left(S_{\mathrm{14}}\right) = \mathrm{0} (12)

其次外部载荷的功可以写为式(13)：

\delta W_{\mathrm{e}\mathrm{x}\mathrm{t}} = f_{\mathrm{1}} \cdot \delta {\mathit{u}}_{\mathrm{1}} + f_{\mathrm{2}} \cdot \delta {\mathit{u}}_{\mathrm{2}} + f_{\mathrm{4}} \cdot \delta {\mathit{u}}_{\mathrm{4}} (13)

式(13)中： f_i为施加在顶点 i每单位密度的外力。

式(11)和式(13)中： \delta {\mathit{u}}_i为保留节点的位移变化，可以用变形梯度的变化量表示，也可以用式(8)表示。

\delta {\mathit{u}}_i = \delta {\mathit{F}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}} \cdot {\mathit{y}}_{\mathrm{0}i} , i = \mathrm{1,2} , \mathrm{4} (14)

将式(11)和式(13)计算在内，考虑到变形梯度张量 {\mathit{F}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}的任意变化都应满足质量变化，最终得到平均应力张量的表达式如下：

{\overline{\mathit{\sigma }}}_{\mathrm{R}\mathrm{V}\mathrm{E}} = \frac{\mathrm{1}}{A} \left[f_{\mathrm{1}}{\mathit{y}}_{\mathrm{1}}+f_{\mathrm{2}}{\mathit{y}}_{\mathrm{2}}+f_{\mathrm{4}}{\mathit{y}}_{\mathrm{4}}\right] = \frac{\mathrm{1}}{A} \sum _{i=\mathrm{1,2},\mathrm{4}} f_i{\mathit{y}}_i (15)

式(15)中：A为RVE的当前面积。

如果力 f_i处于旋转平衡状态(这是RVE水平收敛解的情况)，则 {\overline{\mathit{\sigma }}}_{\mathrm{R}\mathrm{V}\mathrm{E}}的对称性得到保证。基于有限元虚拟场载荷平衡的上述平均应力推导等价于RVE应力场 \mathit{\sigma } \left(\mathit{y}\right)对RVE当前体积V的简单体积平均值，如下式所示：

{\overline{\mathit{\sigma }}}_{\mathrm{R}\mathrm{V}\mathrm{E}} = \frac{\mathrm{1}}{V} \underset{\mathit{y}\in V}{\int } \mathit{\sigma }\left(\mathit{y}\right)\mathrm{d}V (16)

1.6.6 刚度的一致性

当微观方法在非线性有限元程序框架内实现时，需要在每个宏观积分点上建立一致的刚度矩阵。由于该方法在宏观层面上不存在特定形式的系统行为，因此刚度矩阵必须直接由宏观应力与宏观变形的数值关系来确定。在宏观积分点处的一致的宏观刚度矩阵是由RVE模型导出的。通过将后者简化为作用在RVE保留顶点上的力与这些顶点的位移之间的关系，来计算刚度矩阵。因此，将整个RVE-SYSTEM方程组重新排列成如下形式：

\left[\begin{array}{cc}{\mathit{K}}_{\mathrm{r}\mathrm{r}}& {\mathit{K}}_{\mathrm{r}\mathrm{p}}\\ {\mathit{K}}_{\mathrm{p}\mathrm{r}}& {\mathit{K}}_{\mathrm{p}\mathrm{p}}\end{array}\right] \left[\begin{array}{c}\delta {\mathit{u}}_{\mathrm{r}}\\ \delta {\mathit{u}}_{\mathrm{p}}\end{array}\right] = \left[\begin{array}{c}\mathrm{0}\\ \delta f_{\mathrm{p}}\end{array}\right] (17)

式(17)中：用 {\mathit{K}}_{\mathrm{r}\mathrm{r}}、 {\mathit{K}}_{\mathrm{r}\mathrm{p}}、 {\mathit{K}}_{\mathrm{p}\mathrm{r}}、 {\mathit{K}}_{\mathrm{p}\mathrm{p}}划分刚度矩阵 \delta {\mathit{u}}_{\mathrm{p}}和 \delta f为规定保留顶点的迭代位移和外力列阵； \delta {\mathit{u}}_{\mathrm{r}}为相关节点的迭代位移列阵。

在式(17)中，假设对于收敛解，相关节点的(残余)反作用力可以忽略。可以重新计算式(17)，得到将位移变化与顶点力变化关系，从而推导刚度矩阵 {\mathit{S}}_{\mathrm{R}\mathrm{V}\mathrm{E}}。

{\mathit{S}}_{\mathrm{R}\mathrm{V}\mathrm{E}} \delta {\mathit{u}}_{\mathrm{p}} = \delta f_p (18)

{\mathit{S}}_{\mathrm{R}\mathrm{V}\mathrm{E}} = {\mathit{K}}_{\mathrm{p}\mathrm{p}} - {\mathit{K}}_{\mathrm{p}\mathrm{r}} {\mathit{K}}_{\mathrm{r}\mathrm{r}}^{-\mathrm{1}} {\mathit{K}}_{\mathrm{r}\mathrm{p}} (19)

接下来，需要对有关位移变化和力变化的式(18)进行变换，得到有关宏观应力和变形变化的表达式如下：

{{}_{}{}^{\mathrm{4}}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}} : \delta {\mathit{F}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}^{\mathrm{c}} = \delta {\sigma }_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}} (20)

式(20)中： {\mathit{F}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}和 {\mathit{\sigma }}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}分别可由 {\overline{\mathit{F}}}_{\mathrm{R}\mathrm{V}\mathrm{E}}和 {\overline{\mathit{\sigma }}}_{\mathrm{R}\mathrm{V}\mathrm{E}}认定，其中四阶张量 {{}_{}{}^{\mathrm{4}}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}表示宏观积分点水平所需的一致刚度。对于指定顶点为1、2、4的周期性RVE模型通过式(18)定义的简化后的RVE模型的刚度矩阵为 {\mathit{S}}_{\mathrm{R}\mathrm{V}\mathrm{E}}，可划分为式(21)：

{\mathit{S}}_{\mathrm{R}\mathrm{V}\mathrm{E}} \left[\begin{array}{ccc}\mathit{S}\left(\mathrm{11}\right)& \mathit{S}\left(\mathrm{12}\right)& \mathit{S}\left(\mathrm{14}\right)\\ \mathit{S}\left(\mathrm{21}\right)& \mathit{S}\left(\mathrm{22}\right)& \mathit{S}\left(\mathrm{24}\right)\\ \mathit{S}\left(\mathrm{41}\right)& \mathit{S}\left(\mathrm{42}\right)& \mathit{S}\left(\mathrm{44}\right)\end{array}\right] (21)

广义矩阵如下：

\mathit{S} \left(ij\right) = \left[\begin{array}{cc}{\mathit{S}}_{\left(ij\right)}^{\mathrm{11}}& {\mathit{S}}_{\left(ij\right)}^{\mathrm{12}}\\ {\mathit{S}}_{\left(ij\right)}^{\mathrm{21}}& {\mathit{S}}_{\left(ij\right)}^{\mathrm{22}}\end{array}\right] , i , j = \mathrm{1,2} , \mathrm{4} (22)

式(15)中：每个矩阵 \mathit{S} \left(ij\right)都可以认为是二阶张量 \mathit{S} \left(ij\right)的表示。通过这种解释，式(18)被改写为式(23)：

\sum _{j=\mathrm{1,2},\mathrm{4}} \mathit{S}\left(ij\right)\cdot \delta f_i , i = \mathrm{1,2} , \mathrm{4} (23)

由 {\mathit{\sigma }}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}} = {\overline{\mathit{\sigma }}}_{\mathrm{R}\mathrm{V}\mathrm{E}}的条件在各积分点处成立，则应力张量 \delta {\mathit{\sigma }}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}的变化由式(13)推导为：

\delta {\mathit{\sigma }}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}} = - \frac{{\mathit{\sigma }}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}}{A} \delta A + \frac{\mathrm{1}}{A} \sum _{i=\mathrm{1,2},\mathrm{4}} \left(\delta f_i{\mathit{y}}_i+f_i\delta {\mathit{y}}_i\right) (24)

二维构型下RVE面积的变化量 \delta A为：

\delta A = A {\mathit{F}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}^{-\mathrm{c}} : \delta {\mathit{F}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}^{\mathrm{c}} (25)

综上所述，结合公式(25)、(23)、(14)、(24)，可以推导出：

\delta {\mathit{\sigma }}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}} = \left[\frac{\mathrm{1}}{A}\sum _{i=\mathrm{1,2},\mathrm{4}}\left({\mathit{y}}_i\sum _{i=\mathrm{1,2},\mathrm{4}}s\left(ij\right){\mathit{y}}_{\mathrm{0}j}+f_i\mathit{I}{\mathit{y}}_{\mathrm{0}i}\right)-{\mathit{\sigma }}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}{\mathrm{F}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}^{-\mathrm{c}}\right] : \delta {\mathit{F}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}^{\mathrm{c}} (26)

式(26)中：方括号中的表达式等于宏观积分点处所需刚度张量 {{}_{}{}^{\mathrm{4}}}_{\mathrm{m}\mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}}；I为单位张量。

基于上述推导计算可以确定模型在计算过程中每一分析步的整体力学响应。

2 结果与讨论

2.1 短切玻璃纤维复合材料力学性能有限元分析结果

材料属性采用GFRP/MY750本构参数^[12，28]，通过ABAQUS中的用户自定义子程序(UMAT)输入材料参数，分析整体受力时的应力应变。表2为短切玻璃纤维复合材料参数。

在输入材料本构关系之后，完成网格划分，建立有限元模型。图4为随机分布纤维数量为30个且体积分数相同、短切纤维长度不同的细观模型。纤维长度分别为(8、12、19、26、37 mm)。从图4可以看出，由于模型纤维集中分布的区域且纤维方向分布一致的区域纤维所承受的应力明显大于方向不一致的纤维且单一分布的纤维，这是由于与受力方向一致的纤维会传递更多的载荷。同时，当纤维长度较短时，复合材料纤维端部的应力集中比较明显，随着纤维长度的增加，纤维端部应力集中得到缓解。这是由于基体大小不变的情况下，纤维长度越小，纤维在同一大小基体中分布相对较长纤维不集中，与受力方向一致的纤维根数也越多。相反地，与受力方向一致的长纤维的根数较少，因此短纤维的端部树脂基体更容易出现应力集中的情况。

根据1.6的计算方法，提取设置分析步的应力张量、做功位移和积分张量做功，表3为随机分布不同纤维长度模型计算结果。从表3可以看出，当纤维长度较短时(8~19 mm)，弹性模量较小，增大趋势较缓，说明当纤维长度较小时，纤维的增强作用不明显。当纤维长度大于19 mm时，复合材料的弹性模量明显增大，并且在19~26 mm之间时涨幅最大。纤维长度在26~37 mm时，弹性模量仍有增大，但是涨幅稍有降低。在该模型中纤维长度为26 mm时，其整体响应的有效横向模量大约为理论值的两倍，说明在基体大小不变的情况下，纤维的长度越大，纤维增强的作用更为剧烈，且与基体的相互作用更为明显。

图5为随机分布不同纤维长度应力-应变曲线。从图5可以看出，在纤维长度较小时，由于单根纤维与树脂界面面积相对稍长纤维较小，在相同作用力下，较短纤维分担的拉应力较小，整体结构稳定承载能力强。随着纤维长度继续增加，单根纤维上承受的拉应力逐渐变大，整体形变剧烈增大。因此，整体模型的拉伸模量增大，承载能力减弱。随着纤维长度的继续增加，纤维承载作用接近饱和，因此变化逐渐放缓。

2.2 纤维空间分布对比分析

微观结构的空间排列可能对非均质材料的宏观响应有重要影响。为了讨论计算方法对宏观模型结构的随机性行为的影响，将完全规则结构(非均质堆叠)的总体响应与随机结构的总体响应进行比较。选择GFRP材料参数(整体分布最稳定的12 mm样本)。完全规则结构由方形RVE建模，包含30个规则分布的纤维，基体大小为50 mm×50 mm×1 mm。图6为纤维规则分布的模型，图7为算法生成的纤维随机分布的模型，图8为10组模型的变形云图。从图7和图8可以看出，生成了10个不同的随机RVE。假定这10个RVE的平均行为能够代表具有给定非均质体积分数的真实随机结构。首先考虑宏观试样的整体应力应变。在这种情况下，宏观变形场是均匀的，因此不需要进行完整的宏观建模，只需要对具有适当边界条件的独立RVE进行分析即可。因此，分别比较随机结构和规则结构受力情况的模拟的结果。

从图6可以看出，规则结构树脂基体发生巨大形变响应，并且每一根纤维的应力集中位置也相同。由于每一个纤维的排列方向相同(即受力方向相同)，分配给整体上的各个宏观点的所有载荷基本相同。从图7和图8可以看出，随机结构表现的整体响应更为明显，受力条件不变的情况下未出现较大的形变响应。

图9为规则结构和随机结构弹性模量对比。从图9可以看出，在受力相同的情况下，规则结构的加载相对平稳，而随机结构中的一些纤维分布方向集中带已经积累了显著的塑性应变。因此，与随机结构相比，规则结构导致高估了整体刚度。而随机组平均值与理论值的偏差仅为2%左右，因此在整个计算过程中，随机性的影响可以小于纤维长度的影响对整体模型承载能力的影响，因此该模型的合理性也得到了验证。

表4为随机结构有效弹性模量及有效屈服强度。从表4可以出，随机分布模型整体响应的有效屈服强度基本位于125~130 MPa之间，其与理论值误差为13.82%。

选择有效屈服强度最大(随机分布结构3)和最小(随机分布结构9)两组数据的应力进行分析，图10为分析结果。从图10可以看出，由于施加外力均匀递增，随机分布结构3表现出较为稳定的整体响应。随机分布结构9由于纤维分布的随机性，方向不同的纤维在相互作用下会出现纤维冲突处的应力集中，其中最大值约为141 MPa，最小值约为37 MPa。其原因可能为随机生成的纤维在某一方向上出现了大量集中导致纤维分担的应力较小，这情况也符合纤维分布的分散性以及工艺生产中的实际情况。

3 结论

本文利用Python语言对ABAQUS进行二次开发，生成短切纤维随机分布应力模型，建立RVE，并结合宏观刚度矩阵与积分节点矢量组合划分对模型整体的宏观响应进行计算，得到了模型的应力-应变关系以及弹性常数。

在周期性边界条件下，分析不同纤维长度对玻璃纤维增强环氧树脂片状模力学性能的影响，随着纤维长度的增加，玻璃纤维复合材料的弹性模量出现逐渐增大的趋势。

空间分布对复合材料的性能常数的影响。结果表明：纤维随机分布模型计算值更接近理论值，特别是弹性常数的计算。这主要是因为纤维随机分布模型更接近材料的真实结构。

在一定范围内，短切玻璃纤维长度的增加，有利于提高玻璃纤维增强环氧树脂片状模塑料的弹性模量。但随着纤维长度的增加，其在基体内部不均匀分布的情况加剧，应力水平出现差异，从而导致承载集中。因此，需要进一步优化复合材料的制备工艺来提高短切玻璃纤维的分散性和浸润性，从而进一步提高复合材料的力学性能。

参考文献

[1] 王明道,姜宁,赵新涛,纤维增强热塑性树脂基复合材料固相复合成型工艺研究进展[J].材料科学与工艺,2024,32(4):95-110.

[2] 党婧,王凯,肇研.碳纤维增强热塑性树脂基复合材料界面性能研究[J].化工新型材料,2024,52(4):142-145.

[3] 陈平,陆春,于祺,纤维增强热塑性树脂基复合材料界面研究进展[J].材料科学工艺,2007,15(5):665-669.

[4] ATIQAH A, JAWAID M, SAPUAN S M, et al. Thermal properties of sugar palm/glass-fiber reinforced thermoplastic polyurethane hybrid composites[J]. Composite Structures, 2018, 202(5): 954-958.

[5] 蔡建.玻璃钢成型技术[J].工程塑料应用,2003(2):66-70.

[6] ZHANG D H, HE M, QIN S H, et al. Effect of fiber length and dispersion on properties of long glass fiber reinforced thermoplastic composites based on poly(butylene terephthalate)[J]. Rsc Advances, 2017, 7(25): 15439-15454.

[7] 胡婉欣.玻璃纤维增强尼龙6吸波复合材料的制备与性能研究[D].西安:西安建筑科技大学,2024.

[8] UNTERWEGER C, BRÜGGEMANN O, FÜRST C. Effects of different fibers on the properties of short-fiber-reinforced polypropylene composites[J]. Composites Science & Technology, 2014, 103: 49-55.

[9] FUKUDA H, CHOU T W. A probabilistic theory of the strength of short-fibre composites with variable fibre length and orientation[J]. Journal of Materials Science, 1982, 17(4): 1003-1011.

[10] ZHU Y T, ZONG G, MANTHIRAM A, et al. Strength analysis of random short-fibre-reinforced metal matrix composite materials[J]. Journal of Materials Science 1994, 29(23): 6281-6286.

[11] MORTAZAVIAN S, FATEMI A. Effects of fiber orientation and anisotropy on tensile strength and elastic modulus of short fiber reinforced polymer composites[J]. Composites Part B:Engineering, 2015, 72: 116-129.

[12] KOUZNETSOVA V, BREKELMANS W A M, BAAI JENS F P T. An approach to micro-macro modeling of heterogeneous materials[J]. Computational Mechanics, 2001, 27: 37-48.

[13] CHEN Y L, GHOSH S. Micromechanical analysis of strain rate-dependent deformation and failure in composite microstructures under dynamic loading conditions[J]. International Journal of Plasticity, 2012, 32-33: 218-247.

[14] MELRO A R, CAMANHO P P, PINHO S T. Influence of geometrical parameters on the elastic response of unidirectional composite materials[J]. Composite Structures, 2012, 94(11): 3223-3231.

[15] MELRO A R, CAMANHO P P, PINHO S T. Generation of random distribution of fibres in long-fibre reinforced composites[J]. Composites Science and Technology, 2008, 68(9): 2092-2102.

[16] OKEREKE M I, BUCKLEY C P, SIVIOUR C R. Compression of polypropylene across a wide range of strain rates[J]. Mechanics of Time-Dependent Materials, 2012, DOI: 10.1007/s11043-012-9167-z.

[17] XIA Z H, ZHANGY F, ELLYIN F. A unified periodical boundary conditions for representative volume elements of composites and applications[J]. International Journal of Solids and Structures, 2003, 40(8): 1907-1921.

[18] DING Y, CAMANHO P P, SILVA A. A new quantitative method to evaluate the spatial distribution of fibres in composites: The degree of randomness[J]. Journal of Materials Research and Technology, 2023, 26: 7248-7258.

[19] 北京航空航天大学.一种颗粒随机分布的复合材料2D细观结构建模方法:CN202010287490.3[P].2020-08-14.

[20] MORIOKA K, TOMITA Y, TAKIGAWA K. High-temperature fracture properties of GFRP composite for aerospace applications[J]. Materials Science and Engineering: A, 2001, 319: 675-678.

[21] 中华人民共和国国家质量监督检验检疫总局,中国国家标准化管理委员会. 纤维增强塑料拉伸性能试验方法:GB/T1447—2005 [S].北京:中国标准出版社,2005.

[22] 尚士朋.基于商业有限元软件的周期性材料等效性能预测均匀化实施新方法[D].大连:大连理工大学,2017.

[23] OKEREKE M I, AKPOYOMARE A I. A virtual framework for prediction of full-field elastic response of unidirectional composites[J]. Computational Materials Science, 2013, 70:82-99.

[24] 李凌岩,周金宇,马馨远.基于3D打印的三维短纤维增强复合材料强度预测分析[J].工程塑料应用,2022,50(1):109-115.

[25] 江真.短切碳纤维/乙烯基酯树脂片状模塑料拉伸性能分析[D].哈尔滨:哈尔滨工业大学,2018.

[26] VAUGHAN T J, MCCARTHY C T. A combined experimental-numerical approach for generating statistically equivalent fibre distributions for high strength laminated composite materials[J]. Composites Science and Technology, 2010, 70(2): 291-297.

[27] 章继峰,王振清,周健生,基于Python-Abaqus复合材料代表性体积元的数值模型[J].宇航材料工艺,2009,39(3):25-29.

[28] HINTON M J, KADDOUR A S, SODEN P D. Failure criteria in fibre-reinforced-polymer composites[M]. Amsterdam: Elsevier Ltd, 2004.