长纤维增强热塑性复合材料多尺度仿真方法

摘要

本发明公开了一种长纤维增强热塑性复合材料多尺度仿真方法,为解决长纤维增强热塑性复合材料由于应变率效应和各向异性的力学特性在仿真中很难准确模拟的问题，步骤：1.对LFRP平板进行不同方向不同应变率下的拉伸试验；2.对LFRP平板测试区域进行x射线CT扫描；3)对LFRP平板进行多尺度仿真验证：1)粘弹粘塑性材料模型及失效；2)仿真验证：使用ABAQUS软件与Digimat软件耦合建立仿真有限元模型，仿真计算三个应变率下的拉伸试验，其中0.001s‑1和1s‑1采用隐式算法，100s‑1采用显示算法；令有限元样件模型一端固定，另一端延长度方向施加载荷，仿真计算得的各力学特征曲线与试验数据进行对比。

权利要求

1.一种长纤维增强热塑性复合材料LFRP多尺度仿真方法，其特征在于，所述的长纤维增强热塑性复合材料LFRP多尺度仿真方法的步骤如下：

1)对LFRP平板进行不同方向不同应变率下的拉伸试验：

(1)截取试验样件；

(2)力学性能试验；

2)对LFRP平板测试区域进行X射线CT扫描；

3)对LFRP平板进行多尺度仿真验证：

(1)粘弹粘塑性材料模型及失效；

(2)仿真验证：

使用ABAQUS工程模拟的有限元软件与Digimat软件耦合建立仿真有限元模型，仿真计算三个应变率下的拉伸试验，其中0.001s-1和1s-1采用隐式算法，100s-1采用显式算法；创建部件、定义材料属性、定义截面属性、设置分析步、定义载荷和边界条件，进行计算分析，最后对结果进行分析，仿真计算得到的各力学特征曲线与试验数据进行对比；

所述的截取试验样件是指：

(1)LFRP平板是采用注射成型工艺，将热塑性基体和长玻璃纤维混合置于356mm×

100mm×3mm的模具中加工而成，长玻璃纤维质量分数为40％；参照标准ISO527，试验样件是分别沿着与基体流动方向呈0度、45度和90度的夹角从三个LFRP平板上截取下来的；

(2)试验样件的中心宽度为5mm；为防止试验样件上夹持部分圆弧过渡区域提前破坏，修改了试验样件的几何尺寸，将过渡圆弧的半径从30mm减小至5mm；

所述的力学性能试验是指：

(1)根据不同应变率选取不同试验机

应变率范围从0.001s-1到0.1s-1的拉伸试验在万能试验机上进行的，而应变率范围从

1s-1到100s-1的试验是在型号为HTM5020的高速拉伸试验机上完成的，材料拉伸试验的应变率范围涵盖了汽车碰撞过程中车身结构的变形速率；

(2)试验样件应变测量

应变通过散斑制作、标定相机、变形量的采集以及图像处理步骤由DIC数字图像相关系统获得，三维相机通过追踪物体表面散斑变化，采用优化的三维数字图像相关性运算法则，为试验提供二维、三维空间内全视野的形貌、位移及应变数据测量，得到多个应变率下三个方向的应力应变曲线；

(3)分析不同力学特性参数与应变率的关系

长纤维增强复合材料的力学特性既是各向异性又受应变率变化的影响，各向异性与纤维随基体流动的分布情况相关，每种应变率下都是0度方向上的刚度和强度最大，90度最小，另外，材料每个方向上的弹性模量和极限应力随应变率的增大而增大，极限应变则随之减小，通过非线性曲线拟合发现这3个特性参数的变化与应变率呈幂函数关系，其中回归方程为 其中，y为弹性模量、极限应力或极限应变，a、b为拟合参数， 为应变率；

所述的对LFRP平板测试区域进行X射线CT扫描的步骤如下：

1)为了准确了解纤维在LFRP平板中的分布情况，对LFRP平板的中心区域E及其周围区域截取样块进行X射线扫描，考虑纤维在LFRP平板中轴线左右两侧分布是对称的，选择一侧A、B、D、G、H区域，样块尺寸为8mm×6mm×3mm；

2)扫描使用的设备为 的nanoVoxel-2000X射线三维显微镜，能实现低于500nm的成像分辨率；

3)对中心区域E截取样块及周围区域截取样块进行X射线扫描，利用VGStudio形象视觉化软件将扫描的切片图像生成为三维图像，通过对比各厚度上的切片图，能清晰地看到LFRP平板沿厚度方向上是分层的，中间纤维多与流动方向呈明显夹角，且纤维排布更加密集，其余地方纤维多数沿着基体的流动方向排布；

4)量化纤维方向的分布情况

为了进一步量化纤维方向的分布情况，采用二阶方向张量表示，见公式(1)其中：x为基体流动方向、z为LFRP平板厚度方向，其中的矩阵为对称阵，对角线上每项值的范围都是0到1，其和为1，非对角线上每项的范围是-0.5到0.5；二阶张量能进一步转化为特征值(λ1，λ2，λ3)与特征向量 几何形式为椭球，每个特征值分别代表了纤维沿其对应特征向量分布的概率，所以大多数纤维沿着第一特征向量排列，且λ1≥λ2≥λ3；

所述的粘弹粘塑性材料模型及失效是指：

(1)根据LFRP纤维方向分布和体积含量在厚度上分层的特点进行简化，将代表性体积单元RVE分成两层，分别代表核心层区和表层区；考虑到核心层区由于纤维密集的影响使得纤维长度值更大，故将核心层的长径比设置的略大一些；

(2)在各相材料的定义中，玻璃纤维采用弹性模型来模拟力学特性，热塑性基体采用粘弹粘塑性模型来模拟其在弹性阶段和塑性阶段所表现出的随应变率变化的力学特性，粘弹粘塑性模型分为粘弹性段和粘塑性段两个部分，其中粘弹性段采用Prony级数来表示，见公式(2)，Prony级数常用于描述聚合物基高分子材料的应力松弛和蠕变行为；

且

其中：t为材料松弛时间；Gt为剪切模量；G0代表应变率无穷大的初始剪切模量；Kt为体积模量；K0代表应变率无穷大的初始体积模量； 和 为各分量的松弛时间， 和 为给定分量松弛时间下剪切模量和体积模量的权重系数；假设材料的泊松比μ为常数0.25；由于应力松弛针对的是同一材料，因此，用于描述材料应力松弛行为的各分量的松弛时间相同，即 对应的剪切和体积模量的权重系数也相同，即基体材料粘塑性阶段采用current yield Norton law的Viscoplastic模型，如下所示其中， 是塑性变形应变率，σY是初始屈服应力，η是粘塑性系数，f是粘塑性阶段的应力，R(p)是硬化应力，m是粘塑性指数；

(3)在复合材料级别设置基于应变的Tsai-Hill各向同性失效准则来预测材料失效，表达式见公式(5)，将每个应变率下三个方向的应力-应变曲线输入MX模块，逆向得到参数，最终在MF模块中将上述参数定义给Tsai-Hill失效因子；

其中，X为最大轴向拉伸应变，Y为最大面内拉伸应变，S为最大横观剪切应变；

确定好纤维和基体材料模型参数后，利用软件中MAP模块将测试区域的纤维方向张量和体积分数映射到试验样件的有限元网格单元上。

说明书

长纤维增强热塑性复合材料多尺度仿真方法

技术领域

[0001] 本发明涉及属于汽车用复合材料研究领域的一种仿真方法，更确切地说，本发明涉及一种长纤维增强热塑性复合材料多尺度仿真方法。

背景技术

[0002] 注射成型工艺作为多数LFRP(长纤维增强复合材料)汽车零部件的加工方法，不仅适用于复杂的几何形状并能将多个部件集成。在加工过程中，由于纤维方向在不同位置及厚度上的分布是变化的，材料表现出各向异性的力学特性。现有的CAE仿真中通常使用由缩放因子修正的各向同性材料，为了提高LFRP的设计效率和满足结构抗撞性要求，有必要研究高精度且考虑应变率效应的LFRP 模拟方法来预测其力学行为和失效。

[0003] LFRP的多尺度模拟的目的在于依据复合材料细观力学理论，构建一个包含各相材料本构和细观参数的具有等效力学特性的材料模型。相较于宏观各向异性的本构模型，多尺度模拟能充分考虑材料细观结构的影响，模拟地更加细致准确。另外，相较于体现材料详细纤维基体结构的细观力学模型，多尺度模拟可以进行复杂形状和大尺寸结构的仿真计算，显著提高计算效率，节约建模成本。

[0004] Digimat是一款多尺度材料模拟软件，适用于预测弹性和非线性异质材料的本构关系。该软件使用基于Eshelby夹杂理论的均一化方法和Mori-Tanaka模型进行细观力学计算，生成考虑注射成型纤维方向的细观材料模型给后续的有限元分析。

[0005] Zeng利用Digimat-ABAQUS耦合分析预测30％LFRP注射平板的材料特性，并且评价了复杂仪表板的NVH性能。此外，他还建立了考虑厚度上变化的纤维分布和纤维成束影响的材料模型。Wedgewood构建了能够体现非线性应变率相关、各向异性和失效的LFRP模型，使用弹性-粘塑性材料模型(EVP)模拟基体材料，然而，CAE验证的材料响应比试验数据偏软，失效预测需要进一步改进。

发明内容

[0006] 本发明为了解决长纤维增强热塑性复合材料由于应变率效应和各向异性的力学特性在仿真中很难准确模拟的问题，提供了一种考虑应变率效应长纤维增强热塑性复合材料多尺度模拟方法。

[0007] 为解决上述技术问题，本发明是采用如下技术方案实现的：所述的长纤维增强热塑性复合材料多尺度仿真方法的步骤如下：

[0008] 1)对LFRP平板进行不同方向不同应变率下的拉伸试验：

[0009] (1)截取试验样件；

[0010] (2)力学性能试验；

[0011] 2)对LFRP平板测试区域进行x射线CT扫描；

[0012] 3)对LFRP平板进行多尺度仿真验证：

[0013] (1))粘弹粘塑性材料模型及失效；

[0014] (2)仿真验证：

[0015] 使用ABAQUS工程模拟的有限元软件与Digimat软件耦合建立仿真有限元模型，仿真计算三个应变率下的拉伸试验，其中0.001s-1和1s-1采用隐式算法，100s-1采用显示算法；创建部件、定义材料属性、定义截面属性、设置分析步、定义载荷和边界条件，进行计算分析，最后对结果分析，仿真计算得到的各力学特征曲线与试验数据进行对比。

[0016] 技术方案中所述的截取试验样件是指：

[0017] (1)LFRP平板是采用注射成型工艺，将热塑性基体和长玻璃纤维混合置于 356mm×100mm×3mm的模具中加工而成，长玻璃纤维质量分数为40％；参照标准 ISO527，试验样件是分别沿着与基体流动方向呈0度、45度和90度的夹角从三个LFRP平板上截取下来；

[0018] (2)试验样件的中心宽度为5mm；为防止试验样件上夹持部分圆弧过渡区域提前破坏，修改了试验样件的几何尺寸，将过渡圆弧的半径从30mm减小至5mm。

[0019] 技术方案中所述的力学性能试验是指：

[0020] (1)根据不同应变率选取不同试验机

[0021] 应变率范围从0.001s-1到0.1s-1的拉伸实验在万能试验机上进行的，而应变率范围从1s-1到100s-1的试验是在型号为HTM5020的高速拉伸试验机上完成的，材料拉伸试验的应变率范围涵盖了汽车碰撞过程中车身结构的变形速率；

[0022] (2)试验样件应变测量

[0023] 应变通过散斑制作、标定相机、变形量的采集以及图像处理步骤由DIC数字图像相关系统获得，三维相机通过追踪物体表面散斑变化，采用优化的三维数字图像相关性运算法则，为试验提供二维、三维空间内全视野的形貌、位移及应变数据测量，得到多个应变率下三个方向的应力应变曲线；

[0024] (3)分析不同力学特性参数与应变率的关系

[0025] 长纤维增强复合材料的力学特性既是各向异性又受应变率变化的影响，各向异性与纤维随基体流动的分布情况相关，每种应变率下都是0度方向上的刚度和强度最大，90度最小，另外，材料每个方向上的弹性模量和极限应力随应变率的增大而增大，极限应变则随之减小，通过非线性曲线拟合发现这3个特性参数的变化与应变率呈幂函数关系，其中回归方程为

[0026] 技术方案中所述的对LFRP平板测试区域进行x射线CT扫描的步骤如下：

[0027] 1)为了准确了解纤维在LFRP平板中的分布情况，对LFRP平板的中心区域 E及其周围区域截取样块进行X射线扫描，考虑纤维在LFRP平板中轴线左右两侧分布是对称的，选择一侧A、B、D、G、H区域，样块尺寸为8mm×6mm×3mm；

[0028] 2)扫描使用的设备为 的nanoVoxel-2000X射线三维显微镜，可以实现低于500nm的成像分辨率；

[0029] 3)对中心区域E截取样块及周围截取样块进行X射线扫描，利用VGStudio (一款形象视觉化软件)将扫描的切片图像生成为三维图像，通过对比各厚度上的切片图，可以清晰地看到LFRP平板沿厚度方向上是分层的，中间纤维多与流动方向呈明显夹角，且纤维排布更加密集，其余地方纤维多数沿着基体的流动方向排布；

[0030] 4)量化纤维方向的分布情况

[0031] 为了进一步量化纤维方向的分布情况，采用二阶方向张量表示，见公式(1)[0032]

[0033] 其中：x为基体流动方向、z为LFRP平板厚度方向，该矩阵为对称阵，对角线上每项值的范围都是0到1，其和为1，非对角线上每项的范围是-0.5到 0.5；二阶张量可进一步转化为特征值(λ1，λ2，λ3)(λ1≥λ2≥λ3)与特征向量 几何形式为椭球，每个特征值分别代表了纤维沿其对应特征向量分布的概率，所以大多数纤维沿着第一特征向量排列。

[0034] 技术方案中所述的粘弹粘塑性材料模型及失效是指：

[0035] (1)根据LFRP纤维方向分布和体积含量在厚度上分层的特点进行简化，将RVE分成两层，分别代表核心层区和表层区。考虑到核心层区由于纤维密集的影响使得纤维长度值更大，故将核心层层的长径比设置的略大一些；

[0036] (2)在各相材料的定义中，玻璃纤维采用弹性模型来模拟力学特性,热塑性基体采用粘弹粘塑性模型来模拟其在弹性阶段和塑性阶段所表现出的随应变率变化的力学特性,粘弹粘塑性模型分为粘弹性段和粘塑性段两个部分，其中粘弹性段采用Prony级数来表示，见公式(2)，Prony级数常用于描述聚合物基高分子材料的应力松弛和蠕变行为；

[0037]

[0038]

[0039] 其中：t为材料松弛时间；Gt为剪切模量；G0代表应变率无穷大的初始剪切模量；Kt为体积模量；K0代表应变率无穷大的初始体积模量； 和 为各分量的松弛时间， 和为给定分量松弛时间下剪切和体积模量的权重系数；本技术方案中假设材料的泊松比为常数0.25，由于公式4中Gt、Kt和Et的相应系数比相同可以得出

[0040]

[0041] 基体材料粘塑性阶段采用的current yield Norton law的VP模型，如下所示[0042]

[0043] 其中， 是塑性变形应变率，σY是初始屈服应力，η是粘塑性系数，f是粘塑性阶段的应力，R(p)是硬化应力，m是粘塑性指数；

[0044] (3)本专利在复合材料级别设置基于应变的Tsai-Hill各向同性失效准则来预测材料失效，表达式见公式(5)，将每个应变率下三个方向的应力-应变曲线输入MX模块，逆向得到参数，最终在MF中将该函数定义给Tsai-Hill失效因子；

[0045]

[0046]

[0047] 其中，X为最大轴向拉伸应变，Y为最大面内拉伸应变，S为最大横观剪切应变；

[0048] 确定好纤维和基体材料模型参数后，利用软件中MAP模块将测试区域的纤维方向张量和体积分数映射到试验样件的有限元网格单元上。

[0049] 与现有技术相比本发明的有益效果是：

[0050] 1.本发明所述的长纤维增强热塑性复合材料多尺度仿真方法通过LFRP不同方向不同应变率下的拉伸试验，研究了材料各向异性和应变率相关的力学行为，利用非线性拟合的方法确定材料特性的应变率敏感范围。

[0051] 2.本发明所述的长纤维增强热塑性复合材料多尺度仿真方法借助X射线CT 扫描技术获得在测试区域的纤维分布，并采用二阶方向张量量化纤维方向的分布情况，分析纤维随基体流动过程中变化趋势。

[0052] 3.本发明所述的长纤维增强热塑性复合材料多尺度仿真方法采用粘弹粘塑性模型仿真包含应变率效应、失效模式和各向异性的复合材料，在试验样件级别利用有限元仿真计算验证模型有效性。所提出的材料模型可以用于预测汽车车身在不同加载应变率工况下响应，进一步实现对材料细观结构的设计，快速高效地满足抗撞性和轻量化要求的车身部件。

附图说明

[0053] 下面结合附图对本发明作进一步的说明：

[0054] 图1为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法的流程图；

[0055] 图2为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中质量分数40％LFRP平板和试验样件裁剪示意图；

[0056] 图3为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中拉伸试验样件几何尺寸修改图；

[0057] 图4为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中不同位置试验样件截取示意图；

[0058] 图5-1为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中区域E 距上表面0.3mm的电子显微镜扫描切片图；

[0059] 图5-2为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中区域E 距上表面0.9mm的电子显微镜扫描切片图；

[0060] 图5-3为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中区域E 距上表面1.5mm的电子显微镜扫描切片图；

[0061] 图5-4为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中区域E 距上表面2.1mm的电子显微镜扫描切片图；

[0062] 图5-5为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中区域E 距上表面2.7mm的电子显微镜扫描切片图；

[0063] 图6为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中二阶张量定义的椭球；

[0064] 图7为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法的纤维方向张量中axx、axy、ayy在厚度上变化图；

[0065] 图8为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中第一特征向量与X轴夹角在厚度上变化图；

[0066] 图9为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中长纤维增强复合材料纤维体积分数在厚度上变化图；

[0067] 图10为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中区域E、 A、B、G、D、H的第一特征向量与X轴夹角在厚度上变化图；

[0068] 图11-1为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中LFRP 在0度方向不同应变率下的应力应变曲线；

[0069] 图11-2为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中LFRP 在45度方向不同应变率下的应力应变曲线；

[0070] 图11-3为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中LFRP 在90度方向不同应变率下的应力应变曲线；

[0071] 图12-1a为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中 LFRP在0度方向弹性模量试验结果与拟合曲线对比图；

[0072] 图12-2a为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中 LFRP在45度方向弹性模量试验结果与拟合曲线对比图；

[0073] 图12-3a为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中 LFRP在90度方向弹性模量试验结果与拟合曲线对比图；

[0074] 图12-4b为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中 LFRP在0度方向极限应力试验结果与拟合曲线对比图；

[0075] 图12-5b为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中 LFRP在45度方向极限应力试验结果与拟合曲线对比图；

[0076] 图12-6b为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中 LFRP在90度方向极限应力试验结果与拟合曲线对比图；

[0077] 图12-7c为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中 LFRP在0度方向极限应变试验结果与拟合曲线对比图；

[0078] 图12-8c为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中 LFRP在45度方向极限应变试验结果与拟合曲线对比图；

[0079] 图12-9c为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中 LFRP在90度方向极限应变试验结果与拟合曲线对比图；

[0080] 图13-1为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中RVE (代表性体积单元)为0度方向不同应变率计算结果与试验曲线对比图；

[0081] 图13-2为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中RVE (代表性体积单元)为45度方向不同应变率计算结果与试验曲线对比图；

[0082] 图13-3为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中RVE (代表性体积单元)为90度方向不同应变率计算结果与试验曲线对比图；

[0083] 图14-1为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中CAE (计算机辅助工程)为0度方向不同应变率耦合计算结果与试验曲线对比图；

[0084] 图14-2为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中CAE (计算机辅助工程)为45度方向不同应变率耦合计算结果与试验曲线对比图；

[0085] 图14-3为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中CAE (计算机辅助工程)为90度方向不同应变率耦合计算结果与试验曲线对比图；

[0086] 图15为本发明所述的长纤维增强热塑性复合材料多尺度仿真方法中质量分数40％LFRP不同方向不同应变率极限应力的极坐标图。

具体实施方式

[0087] 下面结合附图对本发明作详细的描述：

[0088] 本发明所述的考虑应变率效应的长纤维增强热塑性复合材料多尺度仿真方法流程图如图1所示，首先沿着与基体流动方向呈0度、45度和90度的夹角从多个平板上截取样件，做应变率范围从0.001s-1到100s-1的拉伸试验；同时对平板测试区域截取样块进行X射线扫描准确了解纤维在LFRP中的分布，并采用二阶方向张量量化纤维方向的分布情况；最后在扫描结果的基础上，建立表征纤维方向张量和体积分数的各向异性材料模型实现LFRP的多尺度仿真。

[0089] 所述的长纤维增强热塑性复合材料多尺度仿真方法如下：

[0090] 1.对LFRP平板进行不同方向不同应变率下的拉伸试验：

[0091] 1)截取试验样件

[0092] (1)参阅图2，LFRP平板是采用注射成型工艺，将热塑性基体和长玻璃纤维混合置于356mm×100mm×3mm的模具中加工而成，长玻璃纤维质量分数为40％。参照标准ISO527，试验样件是分别沿着与基体流动方向呈0度、45度和90度的夹角从三个LFRP平板上截取下来；

[0093] (2)参阅图3，试验样件的中心宽度为5mm；为防止试验样件上夹持部分圆弧过渡区域提前破坏，修改了试验样件的几何尺寸，将过渡圆弧的半径从30mm 减小至5mm。

[0094] 2)力学性能试验

[0095] 由于纤维分布的变化，平板各位置的性能有所不同，为保证测量结果的一致性，计算试验样件中心区域的应变来进行比较。

[0096] (1)根据不同应变率选取不同试验机

[0097] 应变率范围从0.001s-1到0.1s-1的拉伸实验在万能试验机上进行的，而应变率范围从1s-1到100s-1的试验是在高速拉伸试验机(HTM5020)上完成的，材料拉伸试验的应变率范围基本涵盖了汽车碰撞过程中车身结构的变形速率。

[0098] (2)试验样件应变测量

[0099] 应变通过散斑制作、标定相机、变形量的采集以及图像处理等步骤由DIC(数字图像相关系统)分析获得，三维相机通过追踪物体表面散斑变化，采用优化的三维数字图像相关性运算法则，为试验提供二维、三维空间内全视野的形貌、位移及应变数据测量，得到多个应变率下三个方向的应力应变曲线。

[0100] (3)分析不同力学特性参数与应变率的关系

[0101] 分析不同应变率下的应力应变曲线，曲线的终点代表材料的极限应力和应变，此时试验样件完全断裂，应力随后骤降为0。长纤维增强复合材料的力学特性既是各向异性又受应变率变化的影响。各向异性与纤维随基体流动的分布情况相关，每种应变率下都是0度方向上的刚度和强度最大，90度最小。另外，材料每个方向上的弹性模量和极限应力随应变率的增大而增大，极限应变则随之减小，通过非线性曲线拟合发现这3个特性参数的变化与应变率呈幂函数关系，其中回归方程为

[0102] 2.对LFRP平板测试区域进行x射线CT扫描：

[0103] 1)参阅图4，为了准确了解纤维在LFRP平板中的分布情况，对LFRP平板的中心区域E及其周围区域(考虑纤维在LFRP平板中轴线左右两侧分布是对称的，选择一侧A、B、D、G、H区域)截取样块进行X射线扫描，样块尺寸为8mm ×6mm×3mm。

[0104] 2)扫描使用的设备为三英精密仪器 的nanoVoxel-2000X射线三维显微镜，可以实现低于500nm的成像分辨率。

[0105] 3)对中心区域E截取样块及周围截取样块进行X射线扫描，利用VGStudio (一款形象视觉化软件)将扫描的切片图像生成为三维图像，通过对比各厚度上的切片图，如图5，可以清晰地看到LFRP平板沿厚度方向上是分层的，中间纤维多与流动方向呈明显夹角，且纤维排布更加密集，其余地方纤维多数沿着基体的流动方向排布。

[0106] 4)量化纤维方向的分布情况

[0107] 为了进一步量化纤维方向的分布情况，采用二阶方向张量表示，见公式(1)[0108]

[0109] 通常情况下，x为基体流动方向、z为LFRP平板厚度方向，该矩阵为对称阵，对角线上每项值的范围都是0到1，其和为1，非对角线上每项的范围是-0.5 到0.5。二阶张量可进一步转化为特征值(λ1，λ2，λ3)(λ1≥λ2≥λ3)与特征向量 几何形式为椭球，见图6。每个特征值分别代表了纤维沿其对应特征向量分布的概率，所以大多数纤维沿着第一特征向量排列。

[0110] a.对于采用注射成型工艺的LFRP平板，绝大多数纤维分布在xy平面内，方向向量主要由axx、axy(ayx)、ayy来确定，其它项(azz、axz、azx和ayz)近似为 0，当axx＝1，axy＝ayy＝0时，纤维沿着x轴分布；当ayy＝1，axy＝axx＝0时，纤维沿着 y轴分布；当axx＝ayy＝0.5，axy＝0时，纤维随机分布；当axx＝ayy＝axy＝0.5时，纤维沿着45度方向分布。将区域E的切片图沿厚度方向均匀分成40组，先计算出每一组的平均纤维方向张量，如图7所示，纤维分布在厚度方向上可沿中间面对称划分成的5个区域，将其命名为表层区、过渡层区和核心层区。

[0111] b.求得每组 与x轴夹角，如图8所示在表层区第一特征向量与x轴的夹角约为10度，说明大多数纤维基本上沿着基体流动方向分布，过渡层区的夹角从 0度到60度都有且λ1和λ2约为0.5，说明该区域内纤维呈随机分布，核心层区的夹角约为-70度，与图5显示的情况一致，研究发现纤维接近横向的排列是由于纤维受到基体低剪切流动和高拉伸流动的影响而造成的。

[0112] c.参阅图9，图中列出LFRP平板纤维体积分数在厚度方向上的变化，可以看出在核心层区和模具上下表面体积分数会明显增加。这说明如果在材料仿真中各厚度材料都设置同一纤维含量值是无法反映纤维的真实分布情况的。

[0113] d.为了研究纤维在测试区域随基体流动的分布情况，本技术方案对比了所有扫描区域纤维第一特征向量与x轴的偏角的变化趋势，如图10。其中区域A 和B接近模具内基体流动的末端，区域G和D接近模具一侧的边缘，区域H和E 都在平板中间位置且基体从H流向E。分析曲线可以发现A、B、D和G区域的纤维方向与E差别很大，尤其是核心层区纤维与x轴的偏角明显减小，是由于模具腔壁影响基体剪切和拉伸流动造成的。相反，代表了平板中间区域的H和E，其偏角基本未发生变化。

[0114] 综上可以看出，测试区域所在LFRP平板中间位置的纤维分布情况很稳定，区域E的方向张量即可代表；而在靠近平板边缘的地方，核心层区的纤维的分布变化较大，表层区比较稳定。

[0115] 3.对LFRP平板进行多尺度仿真验证

[0116] 1)粘弹粘塑性材料模型及失效

[0117] 本发明利用Digimat(一款多尺度非线性材料仿真软件)中MF模块(一种用于多相材料非线性材料本构预测工具)来计算包含材料细观结构参数及各相材料模型的RVE(代表性体积单元)在不同应变率下不同方向材料的应力应变曲线。所谓的RVE是复合材料细观力学中的一个概念，它的平均性质能够描述复合材料宏观有效性能，因此要求RVE满足一定的尺度关系，即它相对于结构尺度要充分小，在复合材料结构中可以看作一个点；另外相对夹杂(也就是长纤维) 的尺度又要充分大，能包含足够多的夹杂。

[0118] (1)根据LFRP纤维方向分布和体积含量在厚度上分层的特点进行简化，将RVE分成两层，分别代表核心层区和表层区。考虑到核心层区由于纤维密集的影响使得纤维长度值更大，故将核心层层的长径比设置的略大一些。

[0119] (2)在各相材料的定义中，玻璃纤维采用弹性模型来模拟力学特性。热塑性基体采用粘弹粘塑性模型来模拟其在弹性阶段和塑性阶段所表现出的随应变率变化的力学特性。粘弹粘塑性模型分为粘弹性段和粘塑性段两个部分，其中粘弹性段采用Prony级数来表示，见公式2，Prony级数常用于描述聚合物基高分子材料的应力松弛和蠕变行为。

[0120]

[0121]

[0122] 其中，t为材料松弛时间；Gt为剪切模量；G0代表应变率无穷大的初始剪切模量；Kt为体积模量；K0代表应变率无穷大的初始体积模量； 和 为各分量的松弛时间， 和为给定分量松弛时间下剪切和体积模量的权重系数。本技术方案中假设材料的泊松比为常数0.25，由于公式4中Gt、Kt和Et的相应系数比相同可以得出

[0123]

[0124] 基体材料粘塑性阶段采用的current yield Norton law的VP模型，如下所示[0125]

[0126] 其中， 是塑性变形应变率，σY是初始屈服应力，η是粘塑性系数，f是粘塑性阶段的应力，R(p)是硬化应力，m是粘塑性指数。

[0127] (3)本专利在复合材料级别设置基于应变的Tsai-Hill各向同性失效准则来预测材料失效，表达式见公式5，将每个应变率下三个方向的应力-应变曲线输入MX模块(Digmat的材料数据库)，逆向得到参数，最终在MF中将该函数定义给Tsai-Hill失效因子。

[0128]

[0129] 其中，X为最大轴向拉伸应变，Y为最大面内拉伸应变，S为最大横观剪切应变。

[0130] 确定好纤维和基体材料模型参数后，利用软件中MAP模块(一种将工艺分析映射到结构分析网格上的工具)将测试区域的纤维方向张量和体积分数映射到试验样件的有限元网格单元上。

[0131] 2)仿真验证

[0132] 使用ABAQUS(一套工程模拟的有限元软件)与Digimat软件耦合建立仿真有限元模型，仿真计算三个应变率下的拉伸试验，其中0.001s-1和1s-1采用隐式算法，100s-1采用显示算法。创建部件、定义材料属性、定义截面属性、设置分析步、定义载荷和边界条件，进行计算分析，最后对结果分析。仿真计算得到的各力学特征曲线与试验数据进行对比，具有良好的一致性。

[0133] 实施例：

[0134] 本发明接下来结合实例介绍利用本发明提出的长纤维增强热塑性复合材料多尺度仿真方法的过程。

[0135] 1.对LFRP平板进行不同方向不同应变率下的拉伸试验

[0136] LFRP平板由 采用注射成型工艺将热塑性基体(PP)和长玻璃纤维混合置于356mm*100mm*3mm的模具中加工，纤维质量分数为40％。试验样件(参照标准ISO527)分别沿着与基体流动方向呈0度、45度和90度的夹角从多个平板上截取下来，如图2所示，中心宽度为5mm。为防止样件上夹持部分圆弧过渡区域提前破坏，修改了样件的几何尺寸，将过渡圆弧的半径从30mm 减小至5mm，如图3所示。

[0137] 试验得到的6个应变率下3个方向样件的应力应变曲线如图11-1至图11-3所示，曲线的终点代表材料的极限工程应力和极限工程应变，此时样件完全断裂，应力随后骤降为0。从曲线上可以看出，该材料的力学特性既是各向异性又受应变率变化的影响。各向异性与纤维随基体流动的分布情况相关，每种应变率下都是0度方向上的刚度和强度最大，90度最小。另外，材料每个方向上的弹性模量和极限应力随应变率的增大而增大，极限应变则随之减小，通过非线性曲线拟合发现这3个特性参数的变化与应变率呈幂函数关系，回归方程见公式7，拟合系数a、b和平方差见表1。

[0138] 12-1a至12-9c显示了试验数据与拟合曲线的关系，通过曲线的发展趋势可以看出LFRP在0.001s-1到100s-1的应变率范围内的特性变化最为敏感。

[0139]

[0140] 其中：y代表材料弹塑性力学性能参数，为应变率。

[0141] 表1

[0142]

[0143] 2.对LFRP平板测试区域进行x射线CT扫描

[0144] 对平板的区域E及其周围(考虑纤维在平板中轴线左右两侧分布是对称的，选择一侧A、B、D、G、H区域)截取样块进行X射线扫描，样块尺寸为8mm*6mm*3mm，如图4所示。扫描使用的设备为 的nanoVoxel-2000X射线三维显微镜，可以实现低于500nm的成像分辨率。

[0145] 利用VGStudio(一款用于工业CT数据分析和可视化的软件平台)将扫描的切片图像生成为3D图像，通过对比各厚度上的切片图，如图5，可以清晰地看到LFRP沿厚度方向上是分层的，中间纤维多与流动方向呈明显夹角，且纤维排布更加密集，其余地方纤维多数沿着基体的流动方向排布。

[0146] 采用二阶方向张量的概念来表示，见公式1。

[0147] 通常情况下，x为基体流动方向、z为平板厚度方向，该矩阵为对称阵，对角线上每项值的范围都是0到1，其和为1，非对角线上每项的范围是-0.5到 0.5。二阶张量可进一步转化为特征值(λ1，λ2，λ3)(λ1≥λ2≥λ3)与特征向量 几何形式为椭球，见图8。每个特征值分别代表了纤维沿其对应特征向量分布的概率，所以大多数纤维沿着第一特征向量排列。

[0148] 对于采用注射成型工艺的LFRP，绝大多数纤维分布在xy平面内，方向向量主要由axx、axy(ayx)、ayy来确定，其它项(azz、axz、azx和ayz)近似为0，当 axx＝1，axy＝ayy＝0时，纤维沿着x轴分布；当ayy＝1，axy＝axx＝0时，纤维沿着y轴分布；当axx＝ayy＝0.5，axy＝0时，纤维随机分布；当axx＝ayy＝axy＝0.5时，纤维沿着45度方向分布。本文将区域E的切片图沿厚度方向均匀分成40组，先计算出每一组的平均纤维方向张量，如图7，纤维分布在厚度方向上可沿中间面对称划分成的5个区域，将其命名为表层区、过渡层区和核心层区。

[0149] 其次求得每组 与x轴夹角，如图8，可以看出在表层区第一特征向量与x 轴的夹角约为10度，说明大多数纤维基本上沿着基体流动方向分布，过渡层区的夹角从0度到60度都有且λ1和λ2约为0.5，说明该区域内纤维呈随机分布， Core区的夹角约为-70度，与图5显示的情况一致，研究发现纤维接近横向的排列是由于纤维受到基体低剪切流动和高拉伸流动的影响而造成的。

[0150] 图9列出材料纤维体积分数在厚度方向上的变化，可以看出在核心层区和模具上下表面体积分数会明显增加。这说明如果在材料仿真中各厚度材料都设置同一纤维含量值是无法反映纤维的真实分布情况的。

[0151] 为了研究纤维在测试区域随基体流动的分布情况，本文对比了所有扫描区域纤维第一特征向量与x轴的偏角的变化趋势，如图10。其中区域A和B接近模具内基体流动的末端，区域G和D接近模具一侧的边缘，区域H和E都在平板中间位置且基体从H流向E。分析曲线可以发现A、B、D和G区域的纤维方向与E差别很大，尤其是核心层区纤维与x轴的偏角明显减小，是由于模具腔壁影响基体剪切和拉伸流动造成的。相反，代表了平板中间区域的H和E，其偏角基本未发生变化。综上可以看出，测试区域所在平板中间位置的纤维分布情况很稳定，区域E的方向张量即可代表；而在靠近平板边缘的地方，核心层区的纤维的分布变化较大，表层区比较稳定。

[0152] 3.对LFRP平板进行多尺度仿真验证

[0153] 1)VEVP材料模型及失效

[0154] 用Digimat中MF模块选取应变率为0.001s-1、1s-1和100s-1的试验曲线作为研究对象进行参数的定义，RVE相关参数设置见表2。定义玻璃纤维模型参数，包括拉伸弹性模量E＝72GPa、密度ρ＝2.54g/cm3和泊松比γ＝0.22等。

[0155] 表2

[0156]

[0157]

[0158] 本专利使用软件中MX模块中逆向优化算法来确定VEVP模型的参数，将0 度方向三个应变率下的曲线作为输入条件，计算得到Prony级数的参数为 G0＝3020MPa，K0＝9060MPa， 利用准静态(应变率0.001s-1)三个方
向的试验曲线在MX模块中逆向出模型中σY＝5MPa，R(p) 相关的硬化模量为3.55MPa，硬化指数为179。用0度方向三个应变率下的曲线逆向得到η＝2.3，m＝4.0。通过逆向得到的参数见表3，定义给Tsai-Hill失效因子。

[0159] 表3

[0160]-1
应变率(s ) 轴向应变(mm/mm) 横向应变(mm/mm) 剪切应变(mm/mm)
0.001 0.0597 0.0430 0.0534
1 0.0472 0.0221 0.0366
100 0.0368 0.0250 0.0188

[0161] 在软件MF模型块中计算不同应变率不同角度下材料RVE的响应曲线，计算结果和试验数据的对比如图13-1至图13-3，其中0度100s-1的曲线在塑性段应力的最大误差为18.3％，其余的曲线在弹塑性阶段的特性变化和失效判断上都有很好的一致性。将材料参数映射到拉伸样件有限元网格上，样件有限元模型在厚度方向上划分为10层。再用CAE模块生成两个材料文件，除了长径比其余参数都相同。

[0162] 2)仿真验证

[0163] 使用ABAQUS(一套工程模拟的有限元软件)与Digimat软件耦合建立仿真有限元模型，仿真计算三个应变率下的拉伸试验。创建部件，定义材料属性，其中材料模型为digimat中所创建模型，并采用实体单元，定义截面属性、设置分析步、定义载荷和边界条件，令有限元样件模型一端固定，另一端延长度方向施加载荷，划分网格在厚度方向上分为10层进行分析，最后对计算结果分析。仿真计算得到的各力学特征曲线与试验数据进行对比，具有良好的一致性。

[0164] 图15给出了仿真计算的极限应力随加载角度和加载应变率的变化图，可以清楚看出LFRP各向异性和应变率相关的力学特性，为进一步进行精确LFRP碰撞模拟提供了仿真基础。

[0165] 综上所述，本发明提出的考虑应变率效应的长纤维增强热塑性复合材料多尺度模拟方法通过对试验曲线中特性参数的非线性曲线拟合确定了材料的应变率敏感区间为0.001s-1到100s-1。利用X-rayCT扫描观察到测试区域内的纤维方向分布比较稳定，与区域E的纤维方向张量相同。在材料模型建立中，定义基体材料为粘弹粘塑性模型来预测LFRP的应变率效应，最终用有限元耦合计算验证了模拟的有效性，仿真与试验结果的一致性很好。
本发明选取长玻纤中的玻璃纤维为例，对于碳纤维、玄武岩纤维等长纤维增强热塑性复合材料同样适用。