CN102819678A - 基于应变场的短纤维增强复合材料纤维取向预测方法 - Google Patents

Abstract

本发明涉及一种基于应变场的短纤维增强复合材料纤维取向预测方法。首先，建立复合材料成形过程的有限元模型；第二，进行成形过程有限元分析，获得坯料全局的应变场，提取有限元分析结果中各个单元的第一、第二、第三主应变大小及方向；第三，根据有限元模型中各单元的主应变方向确定复合材料中纤维定向程度最强及最弱的方向，根据各单元第一、第二、第三主应变大小计算纤维取向定向程度参量。本方法能够预测复合材料成形过程中的纤维取向变化及分布，可应用于材料基体保持为固态且初始纤维取向近似随机分布的短纤维增强复合材料成形工艺，具有预测模型简单、实施难度小的优点，克服了现有方法不适用于成形过程中基体保持为固态情形的不足。

Description

基于应变场的短纤维增强复合材料纤维取向预测方法

技术领域

本发明涉及短纤维增强复合材料的取向预测，尤其涉及一种基于应变场的短纤维增强复合材料纤维取向预测方法。

背景技术

短纤维增强复合材料是一种重要的工程材料，在工程中应用广泛，发展前景广阔。对于纤维取向在三维空间随机或近似随机分布的短纤维增强复合材料，其内部的纤维在成形过程中将不再保持随机取向，即纤维取向将趋于某些方向，同时纤维取向也会相应地偏离另一些方向。短纤维增强复合材料的性能与其内部纤维取向分布状态有很密切的关系，所以短纤维增强复合材料成形中纤维取向的预测工作具有重要意义。

现有的纤维取向预测方法及技术主要是针对纤维悬浮液及聚合物熔体，即基体为液态的情形，而对基体为固态的复合材料纤维取向预测方法仍未见报道。国际上最早由Jeffery分析解决了浸没在粘性流体中单根纤维的运动问题，Jeffery针对牛顿流体中的椭球形单纤维,在忽略纤维粒子之间的相互作用力的基础上,建立了由流体运动引起的纤维动力学演化方程，即Jeffery运动方程。此后，诸多学者发展了纤维取向的理论，建立了相应的纤维取向预测方法，并且一些方法已经得到实际应用。如文献“林兰芬.短纤维增强塑料注射成型中三维纤维取向的数值预测[J].材料科学与工艺,1998,(2).”报道了一种短纤维增强复合材料注塑成形的纤维取向预测方法，该方法是在流动分析的基础上采用Folgar-Tucker取向模型进行纤维取向预测，再如商业软件MoldFlow目前已提供了基于熔体动力学的纤维取向预测分析功能。但上述方法仅适用于基体为液态的情形，不能应用于复合材料成形过程中基体为固态的情形，如短纤维增强复合材料的模锻、挤压、拉拔等。由于目前短纤维增强复合材料的制备及成形技术发展很快，故成形过程中基体保持为固态的短纤维增强复合材料纤维取向预测问题亟待解决，有必要提出相应的纤维取向预测方法。

发明内容

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种基于应变场的短纤维增强复合材料纤维取向预测方法，该方法可应用于基体保持为固态且初始纤维取向近似随机分布的短纤维增强复合材料成形工艺。

技术方案

一种基于应变场的短纤维增强复合材料纤维取向预测方法，其特征在于步骤如下：

步骤1：根据短纤维增强复合材料成形工艺参数及材料属性，建立有限元分析模型；

步骤2：进行成形过程有限元分析，获得坯料全局的应变场，提取有限元分析结果中各个单元的第一、第二、第三主应变大小及方向；以i表示有限元模型中的单元编号，

表示第i个单元第一、第二、第三主应变的大小，

表示第i个单元第一、第二、第三主应变的方向向量；所述应变均指工程应变；

步骤3：计算各单元处纤维与第一、第二及第三主应变的夹角的平均值，

为第i个单元处纤维与第一主应变方向所夹锐角平均值，

计算式为：

为第i个单元处纤维与第二主应变方向所夹锐角平均值，

计算式为：

为第i个单元处纤维与第三主应变方向所夹锐角平均值，

计算式为：

所述

的单位均为弧度；

由于

的解析表达式不可积，下面以给出一种

值的近似计算方法。以

为例，采用二元函数数值积分公式求解，首先将计算式中的积分上下限换算至[-1,1]之间，令

得到

的计算式如下：

$Q_1^j=\frac{\mathrm{\π}}{8}{\∫}_{-1}^1{\∫}_{-1}^1\sin \frac{\mathrm{\π}(x+1)}{4}\arctan \left(\tan \frac{\mathrm{\π}(x+1)}{4}\sqrt{{\left(\frac{1+e_2^i}{1+e_1^i}\right)}^2{\cos }^2\frac{\mathrm{\π}(y+1)}{4}{\left(\frac{1+e_3^i}{1+e_1^i}\right)}^2{\sin }^2\frac{\mathrm{\π}(y+1)}{4}}\right)\mathrm{dxdy}$

令：

$F(x,y)=\frac{\mathrm{\π}}{8}\sin \frac{\mathrm{\π}(x+1)}{4}\arctan \left(\tan \frac{\mathrm{\π}(x+1)}{4}\sqrt{{\left(\frac{{1+e}_2^i}{1+e_1^1}\right)}^2{\cos }^2\frac{\mathrm{\π}(y+1)}{4}{\left(\frac{{1+e}_3^i}{{1+e}_1^i}\right)}^2{\sin }^2\frac{\mathrm{\π}(y+1)}{4}}\right)$

则：

$Q_1^i=\underset{-1}{\overset{1}{\∫}}\underset{-1}{\overset{1}{\∫}}F(x,y)\mathrm{dxdy}=\underset{k=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{H}_k{H}_{j}F(x_k,x_j)$

H_k,H_j,x_k,x_j根据n值和积分公式类型进行选取，当选用五点高斯积分时，取值如下：

x1	x2	x3	x4	x5
					-0.90618	-0.53847	0	0.53847	-0.90618
H1	H2	H3	H4	H5
					0.23693	0.47863	0.56889	0.47863	0.23693

将H_k,H_j,x_k,x_j取值带入计算式即可求出的值，按同样方法即可求出。

步骤4：第i个单元处的纤维取向预测如下：

①当

值在0至1之间时，

值与纤维沿

方向定向程度呈反比；

②当

值在1至π/2之间时，值与纤维偏离方向的程度呈正比；

③当纤维沿

方向的定向程度介于

与

之间时，

的值表示了纤维与

所夹锐角的平均值，当

值小于1时，

值与纤维偏向

方向的程度呈反比；当

值大于1时，

值与纤维偏离方向的程度呈正比。

有益效果

本发明提出的一种基于应变场的短纤维增强复合材料纤维取向预测方法，该方法可应用于基体保持为固态且初始纤维取向近似随机分布的短纤维增强复合材料成形工艺，如短纤维增强复合材料的模锻、挤压、拉拔等，解决了现有纤维取向预测方法不适用于基体保持为固态的复合材料成形工艺的问题；该方法可预测制件内部纤维沿哪个方向定向程度最强，沿哪个方向定向程度最弱，同时能得出纤维沿上述方向的定向程度参量值；此外，该方法还具有预测模型简单、实现难度小的优点。

附图说明

图1是Csf/Mg复合材料挤压成形过程示意图；

1-凸模，2-Csf/Mg复合材料，3-凹模；

图2是有限元模型的几何模型机初始网格；

图3是模拟结果中工程应变的第一主应变方向向量图；

图4是第一主应变方向向量图的局部放大图；

图5是模拟结果中工程应变的第一主应变大小云图；

图6是模拟结果中工程应变的第二主应变方向向量图；

图7是第二主应变方向向量图的局部放大图；

图8是模拟结果中工程应变的第二主应变大小云图；

图9是模拟结果中工程应变的第三主应变方向向量图；

图10是第三主应变方向向量图的局部放大图；

图11是模拟结果中工程应变的第三主应变大小云图；

图12是复合材料成形后各个单元处的Q₁值云图；

图13是复合材料成形后各个单元处的Q₂值云图；

图14是复合材料成形后各个单元处的Q₃值云图；

图15是有限元模型中单元1747的位置示意图；

图16是有限元模型中单元909的位置示意图；

图17是复合材料成形后纤维照片拍照位置；

图18是复合材料A点显微照片；

图19是复合材料B点显微照片。

具体实施方式

现结合实施例、附图对本发明作进一步描述：

以真空吸渗挤压制备Csf/Mg复合材料（短切碳纤维增强镁基复合材料）棒材为例，其工艺过程主要是利用渗铸和液态挤压成形原理，使注入挤压筒中的液态镁合金在冲头压力作用下渗入增强纤维预制体中，发生压力下结晶凝固，并随之从挤压凹模口挤出，成形复合材料棒材，附图1为复合材料挤出的过程示意图。在复合材料从挤压凹模口挤出之前，复合材料中的纤维取向近似为三维空间随机取向，而在复合材料在结晶凝固的状态下被挤出之后，复合材料中的纤维取向发生了剧烈改变。在Abaqus平台中基于本发明所提出的方法针对该工艺的纤维取向预测步骤如下：

步骤1：建立有限元分析模型。在Abaqus平台中，采用轴对称有限元分析模型，其几何模型及初始网格划分如附图2所示，模具定义为刚性体，模具与坯料之间的摩擦系数定义为1.2，在坯料的顶端定义恒定的竖直向下速度1mm/s。为避免材料变形过程中网格过度畸变，将坯料的网格定义为ALE网格(任意的欧拉-拉格朗日网格)。

步骤2：将步骤1中建立的有限元模型提交计算，获得包含工程应变张量(Abaqus平台中以“NE”表示)的输出数据库文件（“.odb”文件）。挤压终了时工程应变张量的第一、第二、第三主应变大小及方向如附图3至附图11所示。令i表示有限元模型中的单元编号，

表示各单元工程应变张量的第一、第二、第三主应变的大小， 表示各单元工程应变张量的第一、第二、第三主应变的方向向量。

均可从输出数据库文件中获得。

步骤3：计算各单元处纤维与第一、第二及第三主应变的夹角的平均值，

为第i个单元处纤维与第一主应变方向所夹锐角平均值，计算式为：

为第i个单元处纤维与第二主应变方向所夹锐角平均值，

计算式为：

为第i个单元处纤维与第三主应变方向所夹锐角平均值，

计算式为：

所述

的单位均为弧度。

采用数值积分的方法计算

以

为例，采用二元函数数值积分公式求解，首先将计算式中的积分上下限换算至[-1,1]之间，令

得到

的计算式如下：

$Q_1^j=\frac{\mathrm{\π}}{8}{\∫}_{-1}^1{\∫}_{-1}^1\sin \frac{\mathrm{\π}(x+1)}{4}\arctan \left(\tan \frac{\mathrm{\π}(x+1)}{4}\sqrt{{\left(\frac{1+e_2^i}{1+e_1^i}\right)}^2{\cos }^2\frac{\mathrm{\π}(y+1)}{4}{\left(\frac{1+e_3^i}{1+e_1^i}\right)}^2{\sin }^2\frac{\mathrm{\π}(y+1)}{4}}\right)\mathrm{dxdy}$

令：

$F(x,y)=\frac{\mathrm{\π}}{8}\sin \frac{\mathrm{\π}(x+1)}{4}\arctan \left(\tan \frac{\mathrm{\π}(x+1)}{4}\sqrt{{\left(\frac{{1+e}_2^i}{1+e_1^1}\right)}^2{\cos }^2\frac{\mathrm{\π}(y+1)}{4}{\left(\frac{{1+e}_3^i}{{1+e}_1^i}\right)}^2{\sin }^2\frac{\mathrm{\π}(y+1)}{4}}\right)$

则：

$Q_1^i=\underset{-1}{\overset{1}{\∫}}\underset{-1}{\overset{1}{\∫}}F(x,y)\mathrm{dxdy}=\underset{k=1}{\overset{n}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{H}_k{H}_{j}F(x_k,x_j)$

H_k,H_j,x_k,x_j根据n值和积分公式类型进行选取，当选用五点高斯积分时，取值如下：

x1	x2	x3	x4	x5
					-0.90618	-0.53847	0	0.53847	-0.90618
H1	H2	H3	H4	H5
					0.23693	0.47863	0.56889	0.47863	0.23693

将H_k,H_j,x_k,x_j取值带入计算式即可求出

的值，

按同样方法即可求出。

本步骤所述的工作可通过Abaqus平台中的python脚本语言完成，采用编写命令脚本的方法对输出数据库中所有单元的进行计算，并将计算结果写入输出数据库以便进行结果的可视化显示。

的计算结果如附图12—附图14所示。

步骤4：第i个单元处的纤维取向预测结果如下：

①

方向为纤维定向程度最强的方向，即材料变形后纤维偏向该方向的程度最强，纤维与

方向所夹锐角的平均值最小，其值为的取值范围在0至1之间；

取值越小，说明纤维沿

方向定向程度越强；

取值近似为0时，表示第i个单元处的纤维与

方向近似平行；

②

方向为纤维定向程度最弱的方向，即材料变形后纤维偏离该方向的程度最强，纤维与

方向所夹锐角的平均值最大，其值为

的取值范围在1至π/2之间；

取值越大，说明纤维偏离

方向的程度越强；

取值近似为π/2时，表示第i个单元处的纤维与

方向近似垂直；

③纤维沿

方向的定向程度介于

与

之间，纤维与所夹锐角的平均值为

两两正交且取值大于

小于

当

值小于1时，说明材料变形后纤维偏向了

方向，

取值越小，说明纤维偏向

方向的程度越强；当

值大于1时，说明材料变形后纤维偏离了

方向，

取值越大，说明纤维偏离

方向的程度越强；

下面以编号为909的单元为例分析纤维取向预测结果。909单元的位置如附图15所示。从输出数据库中可以读出单元909的纤维取向数据： ${\stackrel{\→}{n}}_3^{909}=\left(\mathrm{1,0,0}\right),$ $Q_1^{909}=0.391,$ $Q_2^{909}=1.335,$ $Q_3^{909}=1.336.$ 可见

两两正交， $Q_1^{909}\≤Q_2^{909}\≤Q_3^{909}.$

①

方向为单元909处纤维定向程度最强的方向，即材料发生塑性变形后纤维偏向该方向的程度最强，纤维与的夹角的平均值最小，其值为

②

方向为单元909处纤维定向程度最弱的方向，即材料发生塑性变形后纤维偏离该方向的程度最强，纤维与

的夹角平均值最大，其值为

③单元909处纤维沿

方向的定向程度介于与之间。纤维与

的夹角平均值为 $Q_2^{909}=1.335.$

再以编号为1747的单元为例分析纤维取向预测结果。1747单元的位置如附图16所示。从输出数据库中可以读出单元1747的纤维取向数据：

${\stackrel{\→}{n}}_2^{1747}=\left(\mathrm{0,0,1}\right),$ ${\stackrel{\→}{n}}_3^{1747}=(0.514,-\mathrm{0.858,0}),$ $Q_1^{1747}=0.314,$ $Q_2^{1747}=1.326,$ $Q_3^{1747}=1.437.$ 可见

两两正交， $Q_1^{1747}\≤Q_2^{1747}\≤Q_3^{1747}.$

①

方向为单元1747处纤维定向程度最强的方向，即材料发生塑性变形后纤维偏向该方向的程度最强，纤维与

的夹角的平均值最小，其值为 $Q_1^{1747}=0.314.$

②

方向为单元1747处纤维定向程度最弱的方向，即材料发生塑性变形后纤维偏离该方向的程度最强，纤维与

的夹角平均值最大，其值为 $Q_3^{1747}=1.437.$

③单元1747处纤维沿

方向的定向程度介于

与

之间。纤维与

的夹角平均值为 $Q_2^{1747}=1.326.$

下面是对本实施实例预测结果的实验验证：

分别在有限元模型中编号909及1747单元处对应的位置切取复合材料试样，即附图17中的A、B两点，通过扫描电镜获得A、B两点的显微照片，A点显微照片见附图18，B点显微照片见附图19。通过显微照片中的纤维截面进行纤维取向的计算统计，得出A点处纤维沿

的夹角平均值为0.632(预测值为

沿

的夹角平均值为1.146(预测值为

沿

的夹角平均值为1.218(预测值为

通过显微照片中的纤维截面进行纤维取向的计算统计，得出B点处纤维沿的夹角平均值为0.478(预测值为沿

的夹角平均值为1.167(预测值为

沿

的夹角平均值为1.397(预测值为

的取值范围在0至π/2之间，上述预测值与实验值最大相差均在0.25以内，可见预测结果与实验测定结果吻合，说明通过本方法预测的复合材料纤维取向状态可信。

Claims (1)

1.一种基于应变场的短纤维增强复合材料纤维取向预测方法，其特征在于步骤如下：

步骤1：根据短纤维增强复合材料成形工艺参数及材料属性，建立有限元分析模型；

步骤2：进行成形过程有限元分析，获得坯料全局的应变场，提取有限元分析结果中各个单元的第一、第二、第三主应变大小及方向；以i表示有限元模型中的单元编号，

表示第i个单元第一、第二、第三主应变的大小，

表示第i个单元第一、第二、第三主应变的方向向量；所述应变均指工程应变；

步骤3：计算各单元处纤维与第一、第二及第三主应变的夹角的平均值，

为第i个单元处纤维与第一主应变方向所夹锐角平均值，

计算式为：

为第i个单元处纤维与第二主应变方向所夹锐角平均值，

计算式为：

为第i个单元处纤维与第三主应变方向所夹锐角平均值，

计算式为：

所述的单位均为弧度；

步骤4：第i个单元处的纤维取向预测如下：

①当值在0至1之间时，值与纤维沿

方向定向程度呈反比；

②当

值在1至π/2之间时，

值与纤维偏离

方向的程度呈正比；

③当纤维沿

方向的定向程度介于

与

之间时，

的值表示了纤维与

所夹锐角的平均值，当

值小于1时，

值与纤维偏向

方向的程度呈反比；当

值大于1时，

值与纤维偏离方向的程度呈正比。

Images