CN103871059A - 纤维增强复合材料等效弹性参数的计算方法 - Google Patents

Abstract

本发明提供一种纤维增强复合材料等效弹性参数的计算方法，其属于复合材料静力学性能参数计算技术领域。本发明采用XCT技术建立复合材料细观结构模型，模型更加接近真实情况；且纤维束形状、纤维体积比、孔隙率和孔隙的分布等影响弹性模量的重要参数能够自动获得，而不需要进行假设。本发明更为合理，计算准确度更高，流程更为通用，且免去了各种似是而非的假设和复杂的数学推导，更容易让工程人员接受和掌握。

Description

纤维增强复合材料等效弹性参数的计算方法

技术领域：

本发明涉及一种纤维增强复合材料等效弹性参数的计算方法，其属于复合材料静力学性能参数计算技术领域。

背景技术：

纤维增强复合材料具有高比刚度、比强度的特点，还可在极高温度环境下作为结构材料，是航空航天、军事、新能源、汽车等领域的重要材料。纤维增强复合材料分为单向、二维编织、2.5维、三维正交、三维四向等复合材料。不管何种构型的复合材料，当施加低速静态载荷时，应力应变曲线一般包含线性段和非线性段。在线弹性变形阶段应力应变曲线的斜率为复合材料的等效弹性模量。快速准确的计算复合材料的弹性参数，对于评价复合材料的静力学性能以及复合材料结构设计具有重要的实际应用价值。

目前，确定纤维增强复合材料的等效弹性参数的方法有实验法、解析法和数值法。实验法是按照ASTM（American Society of Testing Materials）等试验标准中的相关标准进行静态加载试验，从得到的试验曲线上直接计算出相应的等效弹性参数。实验法需要按照相关标准进行试验，记录响应曲线，操作过程比较繁琐。此外，纤维增强复合材料一般为非各向同性材料，根据复合材料类型的不同，其独立的弹性参数可能为5个或9个，很难通过实验的方法测得所有的弹性参数。

解析法只能计算种类非常有限的结构非常规则的复合材料弹性模量，不能应用于结构复杂的复合材料等效弹性模量计算。目前普遍采用数值法来计算复杂结构复合材料的等效弹性模量。该方法基于均匀化理论，假设复合材料细观结构具有周期性，采用单胞模型来模拟复合材料细观结构，对模型施加周期性边界条件，计算出单胞的应力应变响应，结合单胞的原始尺寸最后计算出等效弹性模量（见博士学位论文：陶瓷基复合材料损伤耦合的宏细观统一本构模型研究，作者：高希光，导师：宋迎东，南京航空航天大学，2007年）。但是，由于工艺上的限制，基体内部通常存在分布不均匀的孔洞，且增强纤维也不可能是均匀分布的。而孔洞和纤维的分布对材料等效弹性模量的影响较大。基于周期性单胞模型很难考虑这些影响因素，因此上述所说的数值法在弹性参数预测时与实验值存在着很大的偏差。

因此，确有必要对现有技术进行改进以解决现有技术之不足。

发明内容：

本发明提供一种纤维增强复合材料等效弹性参数的计算方法，适用于各种纤维增强复合材料的等效弹性参数计算，充分考虑了纤维和孔隙分布不均匀的影响，计算精度高。

本发明采用如下技术方案：一种纤维增强复合材料等效弹性参数的计算方法，其包括如下步骤：

步骤1：将试件样品进行XCT扫描，获得逐层的灰度图像序列；

步骤2：假设试件断层扫描图片序列一共有M幅图片，每幅图片的像素尺寸为宽W像素、高H像素，每个像素的灰度范围是0-255，其中采用(j,i,k)（k=0，1，2～M-1；i=0,1,2～H-1；j=0,1,2～W-1）来表示第k+1幅图片，第i+1行，第j+1列像素，将试件断层扫描图片序列读入内存，建立像素的灰度数组color[]，数组中的元素color[W*H*k+W*i+j]表示像素(j,i,k)的灰度值；

步骤3：确定基体、碳纤维和孔隙的参考灰度值colorM、colorF和colorV，理论上，孔隙的灰度值为0，基体和碳纤维的灰度值分别开展XCT实验，获得纯基体和纯碳纤维的灰度值，且实验参数与试件扫描时的参数相同，纯基体和纯碳纤维试件采用反应烧结SiC板和碳纤维布，碳纤维布用树脂固定；

步骤4：创建单元阵列，其通过一个像素代表一个三维8节点正方体单元，按照步骤2的设定，创建的单元阵列的长、宽和高分别包含W、H和M个单元；

步骤5：确定单元的弹性参数，已知碳纤维和基体的弹性张量分别为

和

其中f和m分别表示纤维和基体，i,j,k,l=1,2,3，孔洞的弹性张量

设定为一个接近零而不为零的数值；

采用线性插值函数来计算单元(j,i,k)的弹性张量，假设co=color[W*H*k+W*i+j]，在纤维束区域，如果colorV<co<=colorF，则单元(j,i,k)的弹性张量

$E_{\mathrm{ijkl}}=\frac{(\mathrm{co}-\mathrm{colorV})}{(\mathrm{colorF}-\mathrm{colorV})}{E}_{\mathrm{ijkl}}^f+\frac{(\mathrm{co}-\mathrm{colorF})}{(\mathrm{colorV}-\mathrm{colorF})}{E}_{\mathrm{ijkl}}^v---\left(1\right)$

如果colorF<co<=colorM，则单元的弹性张量

$E_{\mathrm{ijkl}}=\frac{(\mathrm{co}-\mathrm{colorM})}{(\mathrm{colorF}-\mathrm{colorM})}{E}_{\mathrm{ijkl}}^f+\frac{(\mathrm{co}-\mathrm{colorF})}{(\mathrm{colorM}-\mathrm{colorF})}{E}_{\mathrm{ijkl}}^m---\left(2\right)$

在基体区域，则单元的弹性张量可用下式计算：

$E_{\mathrm{ijkl}}=\frac{(\mathrm{co}-\mathrm{colorM})}{(\mathrm{colorV}-\mathrm{colorM})}{E}_{\mathrm{ijkl}}^v+\frac{(\mathrm{co}-\mathrm{colorV})}{(\mathrm{colorM}-\mathrm{colorV})}{E}_{\mathrm{ijkl}}^m---\left(3\right)$

步骤6：施加边界条件，对步骤4和步骤5建立的单元阵列依次施加如下位移边界条件：

第1组：x＝0,u_x＝0;y＝0,u_y＝0;z＝0,u_z＝0;x＝w,u_x＝w*0.001  (4a)

第2组：x＝0,u_x＝0;y＝0,u_y＝0;z＝0,u_z＝0;y＝l,u_y＝l*0.001  (4b)

第3组：x＝0,u_x＝0;y＝0,u_y＝0;z＝0,u_z＝0;z＝h,u_z＝h*0.001  (4c)

第4组：x＝0,u_x＝u_y＝u_z＝0;x＝w,u_z＝w*0.001  (4d)

第5组：x＝0,u_x＝u_y＝u_z＝0;x＝w,u_y＝w*0.001  (4e)

第6组：y＝0,u_x＝u_y＝u_z＝0;y＝l,u_y＝l*0.001  (4f)

步骤7：计算出复合材料等效弹性参数，依次计算出步骤6所述边界条件下的单元阵列的应力分布和应变分布，然后采用方程（5）计算出平均应力和平均应变

并写成公式（6）所示的向量形式：

${\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{\mathrm{ij}}=\frac{1}{V}{\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&sigma;}}_{\mathrm{ij}}\mathrm{dv},{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{\mathrm{ij}}=\frac{1}{V}{\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&epsiv;}}_{\mathrm{ij}}\mathrm{dv}---\left(5\right)$

${\left[\begin{array}{cccccc}{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{11}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{22}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{33}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{23}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{13}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{12}\end{array}\right]}^T,{\left[\begin{array}{cccccc}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{11}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{22}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{33}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{23}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{13}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{12}\end{array}\right]}^T---\left(6\right)$

将六组边界条件的计算结果填入矩阵（7），即可得到复合材料的等效刚度矩阵对

求逆后得到复合材料的等效柔度矩阵最后由等式(8)可计算出复合材料的等效弹性参数

$\stackrel{\&OverBar;}{S}=\left[s_{\mathrm{ij}}\right],i,j=1~6$

E₁₁＝1/s₁₁,E₂₂＝1/s₂₂,E₃₃＝1/s₃₃

μ₁₂＝-E₁*s₁₂,μ₁₃＝-E₁*s₁₃,μ₂₃＝-E₂*s₂₃

G₂₃＝1/s₄₄,G₁₃＝1/s₅₅,G₁₂＝1/s₆₆     (8)。

所述步骤3中如果没有开展独立的实验来确定colorM和colorF，则通过在color[]数组中寻找最大灰度值colorMax和最小灰度值colorMin，认为colorM=colorMax，colorV=colorMin，在材料内部没有其他杂质时，打开一个断层图片，锁定一个纤维束区域，然后寻找纤维束区域的最大值和最小值，然后取其平均值就是colorF。

本发明具有如下有益效果：

（1）.采用XCT技术建立复合材料细观结构模型，模型更加接近真实情况；

（2）.纤维束形状、纤维体积比、孔隙率和孔隙的分布等影响弹性模量的重要参数能够自动获得，而不需要进行假设；

（3）.相比于现有的计算方法，本发明更为合理，计算准确度更高，流程更为通用，且免去了各种似是而非的假设和复杂的数学推导，更容易让工程人员接受和掌握。

附图说明：

图1是单元阵列示意图，其中w、l和h分别表示单元阵列在x、y和z方向上的尺寸。

图2是平纹编织碳纤维增强碳化硅复合材料试件图。

图3是平纹编织碳纤维增强碳化硅复合材料XCT灰度图片序列，采用XCT机扫描图2所示试件后得到。

图4平纹编织碳纤维增强碳化硅复合材料单元阵列，其中颜色a表示经纱，颜色b表示纬纱，颜色c表示基体。

图5是平纹穿刺碳纤维增强碳化硅复合材料试件图。

图6是平纹穿刺碳纤维增强碳化硅复合材料XCT灰度图片序列，采用XCT机扫描图5所示试件后得到。

图7平纹穿刺碳纤维增强碳化硅复合材料单元阵列，其中颜色d表示经纱，颜色e表示纬纱，颜色f表示厚度方向的穿刺纤维束，颜色g表示基体。

具体实施方式：

本发明纤维增强复合材料等效弹性参数的计算方法包括如下步骤：

步骤1：将试件样品进行XCT扫描，获得逐层的灰度图像序列；

步骤2：假设试件断层扫描图片序列一共有M幅图片，每幅图片的像素尺寸为宽W像素、高H像素，每个像素的灰度范围是0-255，本发明采用(j,i,k)（k=0，1，2～M-1；i=0,1,2～H-1；j=0,1,2～W-1）来表示第k+1幅图片，第i+1行，第j+1列像素，将试件断层扫描图片序列读入内存，建立像素的灰度数组color[]，数组中的元素color[W*H*k+W*i+j]表示像素(j,i,k)的灰度值；

步骤3：确定基体、碳纤维和孔隙的参考灰度值colorM、colorF和colorV，理论上，孔隙的灰度值为0，基体和碳纤维的灰度值必须分别开展XCT实验，获得纯基体和纯碳纤维的灰度值，且实验参数必须与试件扫描时的参数相同，纯基体和纯碳纤维试件可以考虑采用反应烧结SiC板和碳纤维布，碳纤维布需要用树脂固定；

如果没有开展独立的实验来确定colorM和colorF，也可以采用下面的方法，在color[]数组中寻找最大灰度值colorMax和最小灰度值colorMin，认为colorM=colorMax，colorV=colorMin，前提是材料内部没有其他杂质，打开一个断层图片，锁定一个纤维束区域，然后寻找纤维束区域的最大值和最小值，然后取其平均值就是colorF；

步骤4：创建单元阵列，本方法的基本思路是一个像素代表一个三维8节点正方体单元，按照步骤2的设定，创建的单元阵列的长、宽和高分别包含W、H和M个单元；

步骤5：确定单元的弹性参数，已知碳纤维和基体的弹性张量分别为

和

（f和m分别表示纤维和基体，i,j,k,l=1,2,3），孔洞的弹性张量

可以人为设定为一个接近零而不为零的数值；

采用线性插值函数来计算单元(j,i,k)的弹性张量，假设co=color[W*H*k+W*i+j]，在纤维束区域，如果colorV<co<=colorF，则单元(j,i,k)的弹性张量

$E_{\mathrm{ijkl}}=\frac{(\mathrm{co}-\mathrm{colorV})}{(\mathrm{colorF}-\mathrm{colorV})}{E}_{\mathrm{ijkl}}^f+\frac{(\mathrm{co}-\mathrm{colorF})}{(\mathrm{colorV}-\mathrm{colorF})}{E}_{\mathrm{ijkl}}^v---\left(1\right)$

如果colorF<co<=colorM，则单元的弹性张量

$E_{\mathrm{ijkl}}=\frac{(\mathrm{co}-\mathrm{colorM})}{(\mathrm{colorF}-\mathrm{colorM})}{E}_{\mathrm{ijkl}}^f+\frac{(\mathrm{co}-\mathrm{colorF})}{(\mathrm{colorM}-\mathrm{colorF})}{E}_{\mathrm{ijkl}}^m---\left(2\right)$

在基体区域，则单元的弹性张量可用下式计算：

$E_{\mathrm{ijkl}}=\frac{(\mathrm{co}-\mathrm{colorM})}{(\mathrm{colorV}-\mathrm{colorM})}{E}_{\mathrm{ijkl}}^v+\frac{(\mathrm{co}-\mathrm{colorV})}{(\mathrm{colorM}-\mathrm{colorV})}{E}_{\mathrm{ijkl}}^m---\left(3\right)$

步骤6：施加边界条件，对步骤4和步骤5建立的单元阵列依次施加如下位移边界条件：

第1组：x＝0,u_x＝0;y＝0,u_y＝0;z＝0,u_z＝0;x＝w,u_x＝w*0.001  (4a)

第2组：x＝0,u_x＝0;y＝0,u_y＝0;z＝0,u_z＝0;y＝l,u_y＝l*0.001  (4b)

第3组：x＝0,u_x＝0;y＝0,u_y＝0;z＝0,u_z＝0;z＝h,u_z＝h*0.001  (4c)

第4组：x＝0,u_x＝u_y＝u_z＝0;x＝w,u_z＝w*0.001  (4d)

第5组：x＝0,u_x＝u_y＝u_z＝0;x＝w,u_y＝w*0.001  (4e)

第6组：y＝0,u_x＝u_y＝u_z＝0;y＝l,u_y＝l*0.001  (4f)

步骤7：计算出复合材料等效弹性参数，依次计算出步骤6所述边界条件下的单元阵列的应力分布和应变分布，然后采用方程（5）计算出平均应力

和平均应变

并写成公式（6）所示的向量形式：

${\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{\mathrm{ij}}=\frac{1}{V}{\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&sigma;}}_{\mathrm{ij}}\mathrm{dv},{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{\mathrm{ij}}=\frac{1}{V}{\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&epsiv;}}_{\mathrm{ij}}\mathrm{dv}---\left(5\right)$

${\left[\begin{array}{cccccc}{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{11}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{22}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{33}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{23}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{13}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{12}\end{array}\right]}^T,{\left[\begin{array}{cccccc}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{11}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{22}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{33}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{23}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{13}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{12}\end{array}\right]}^T---\left(6\right)$

将六组边界条件的计算结果填入矩阵（7），即可得到复合材料的等效刚度矩阵对

求逆后得到复合材料的等效柔度矩阵

最后由等式(8)可计算出复合材料的等效弹性参数

$\stackrel{\&OverBar;}{S}=\left[s_{\mathrm{ij}}\right],i,j=1~6$

E₁₁＝1/s₁₁,E₂₂＝1/s₂₂,E₃₃＝1/s₃₃

μ₁₂＝-E₁*s₁₂,μ₁₃＝-E₁*s₁₃,μ₂₃＝-E₂*s₂₃

G₂₃＝1/s₄₄,G₁₃＝1/s₅₅,G₁₂＝1/s₆₆      (8)。

实施例1

请参照图2至图4所示，平纹编织碳纤维增强陶瓷基复合材料等效弹性参数计算：

如图2所示为典型的平纹编织碳纤维增强陶瓷基复合材料试件，增强相为碳纤维，基体为碳化硅基体。碳纤维的弹性参数如下：

E₁₁＝221Gpa,E₂₂＝E₃₃＝13.8Gpa,G₁₂＝G₁₃＝9.0Gpa,G₂₃＝9.2Gpa,μ₁₂＝μ₁₃＝0.2,μ₂₃＝0.25

基体的弹性参数如下：

E＝450Gpa,μ＝0.25

步骤1：将试件样品进行XCT扫描，获得图3所示的逐层的灰度图像序列；

步骤2：断层扫描图片序列一共有16幅图片，每幅图片的像素尺寸为宽154像素、高83像素，每个像素的灰度范围是0-255，将试件断层扫描图片序列读入内存，建立像素的灰度数组color[]，数组中的元素color[12782*k+154*i+j]表示像素(j,i,k)的灰度值；

步骤3：确定基体、碳纤维和孔隙的参考灰度值colorM、colorF和colorV，在color[]数组中寻找最大灰度值colorMax和最小灰度值colorMin分别等于240和0，设定colorM=colorMax，colorV=colorMin，打开图3中的#2断层图片，锁定一个纤维束区域，纤维束区域的灰度最大值和最小值分别为68和54，其平均值61就是colorF的数值；

步骤4：创建如图3所示的长、宽和高分别包含154、83和16个单元的单元阵列；

步骤5：确定单元的弹性参数，假设孔洞的弹性模量E等于10Pa，采用线性插值函数（1）-（3）来计算单元(j,i,k)的弹性张量，并赋予对应单元；

步骤6：对步骤4和步骤5建立的单元阵列施加方程（4）所述的位移边界条件；

步骤7：计算出单元阵列的应力分布。然后根据公式（5）计算出平均应力和平均应变，并写成如（6）所述的向量形式。将六组边界条件的计算结果填入矩阵（7），即可得到复合材料的等效刚度矩阵

对求逆后得到复合材料的等效柔度矩阵

由

根据等式(8)可计算出复合材料的等效弹性参数如下：

E₁₁＝162GPa,E₂₂＝158GPa,E₃₃＝73GPa,μ₁₂＝0.12,μ₁₃＝0.24,μ₂₃＝0.24

G₂₃＝37GPa,G₁₃＝31GPa,G₁₂＝50GPa。

实施例2

请参照图5至图7所示，平纹穿刺碳纤维增强陶瓷基复合材料等效弹性参数计算：

如图5所示为典型的平纹穿刺碳纤维增强陶瓷基复合材料试件。增强相为碳纤维，基体为碳化硅基体。碳纤维的弹性参数如下：

E₁₁＝221Gpa,E₂₂＝E₃₃＝13.8Gpa,G₁₂＝G₁₃＝9.0Gpa,G₂₃＝9.2Gpa,μ₁₂＝μ₁₃＝0.2,μ₂₃＝0.25

基体的弹性参数如下：

E＝450Gpa,μ＝0.25

步骤1：将试件样品进行XCT扫描，获得图6所示的逐层的灰度图像序列；

步骤2：断层扫描图片序列一共有5幅图片，每幅图片的像素尺寸为宽124像素、高159像素，每个像素的灰度范围是0-255，将试件断层扫描图片序列读入内存，建立像素的灰度数组color[]。数组中的元素color[19716*k+124*i+j]表示像素(j,i,k)的灰度值；

步骤3：确定基体、碳纤维和孔隙的参考灰度值colorM、colorF和colorV，在color[]数组中寻找最大灰度值colorMax和最小灰度值colorMin分别等于255和0，设定colorM=colorMax，colorV=colorMin，打开图6中的#2断层图片，锁定一个纤维束区域，纤维束区域的灰度最大值和最小值分别为109和69，其平均值89就是colorF的数值；

步骤4：创建如图6所示的长、宽和高分别包含124、159和5个单元的单元阵列；

步骤5：确定单元的弹性参数。假设孔洞的弹性模量E等于10Pa，采用线性插值函数（1）-（3）来计算单元(j,i,k)的弹性张量，并赋予对应单元；

步骤6：对步骤4和步骤5建立的单元阵列施加方程（4）所述的位移边界条件；

步骤7：计算出单元阵列的应力分布，然后根据公式（5）计算出平均应力和平均应变，并写成如（6）所述的向量形式，将六组边界条件的计算结果填入矩阵（7），即可得到复合材料的等效刚度矩阵

对

求逆后得到复合材料的等效柔度矩阵

由

根据等式(8)可计算出复合材料的等效弹性参数如下：

E₁₁＝135GPa,E₂₂＝139GPa,E₃₃＝115GPa,μ₁₂＝0.15,μ₁₃＝0.16,μ₂₃＝0.16

G₂₃＝37GPa,G₁₃＝29GPa,G₁₂＝35GPa。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下还可以作出若干改进，这些改进也应视为本发明的保护范围。

Claims (2)

1.一种纤维增强复合材料等效弹性参数的计算方法，其特征在于：包括如下步骤

步骤1：将试件样品进行XCT扫描，获得逐层的灰度图像序列；

步骤2：假设试件断层扫描图片序列一共有M幅图片，每幅图片的像素尺寸为宽W像素、高H像素，每个像素的灰度范围是0-255，其中采用(j,i,k)（k=0，1，2～M-1；i=0,1,2～H-1；j=0,1,2～W-1）来表示第k+1幅图片，第i+1行，第j+1列像素，将试件断层扫描图片序列读入内存，建立像素的灰度数组color[]，数组中的元素color[W*H*k+W*i+j]表示像素(j,i,k)的灰度值；

步骤3：确定基体、碳纤维和孔隙的参考灰度值colorM、colorF和colorV，理论上，孔隙的灰度值为0，基体和碳纤维的灰度值分别开展XCT实验，获得纯基体和纯碳纤维的灰度值，且实验参数与试件扫描时的参数相同，纯基体和纯碳纤维试件采用反应烧结SiC板和碳纤维布，碳纤维布用树脂固定；

步骤4：创建单元阵列，其通过一个像素代表一个三维8节点正方体单元，按照步骤2的设定，创建的单元阵列的长、宽和高分别包含W、H和M个单元；

步骤5：确定单元的弹性参数，已知碳纤维和基体的弹性张量分别为

和

其中f和m分别表示纤维和基体，i,j,k,l=1,2,3，孔洞的弹性张量

设定为一个接近零而不为零的数值；

采用线性插值函数来计算单元(j,i,k)的弹性张量，假设co=color[W*H*k+W*i+j]，在纤维束区域，如果colorV<co<=colorF，则单元(j,i,k)的弹性张量

$E_{\mathrm{ijkl}}=\frac{(\mathrm{co}-\mathrm{colorV})}{(\mathrm{colorF}-\mathrm{colorV})}{E}_{\mathrm{ijkl}}^f+\frac{(\mathrm{co}-\mathrm{colorF})}{(\mathrm{colorV}-\mathrm{colorF})}{E}_{\mathrm{ijkl}}^v---\left(1\right)$

如果colorF<co<=colorM，则单元的弹性张量

$E_{\mathrm{ijkl}}=\frac{(\mathrm{co}-\mathrm{colorM})}{(\mathrm{colorF}-\mathrm{colorM})}{E}_{\mathrm{ijkl}}^f+\frac{(\mathrm{co}-\mathrm{colorF})}{(\mathrm{colorM}-\mathrm{colorF})}{E}_{\mathrm{ijkl}}^m---\left(2\right)$

在基体区域，则单元的弹性张量可用下式计算：

$E_{\mathrm{ijkl}}=\frac{(\mathrm{co}-\mathrm{colorM})}{(\mathrm{colorV}-\mathrm{colorM})}{E}_{\mathrm{ijkl}}^v+\frac{(\mathrm{co}-\mathrm{colorV})}{(\mathrm{colorM}-\mathrm{colorV})}{E}_{\mathrm{ijkl}}^m---\left(3\right)$

步骤6：施加边界条件，对步骤4和步骤5建立的单元阵列依次施加如下位移边界条件：

第1组：x＝0,u_x＝0;y＝0,u_y＝0;z＝0,u_z＝0;x＝w,u_x＝w*0.001  (4a)

第2组：x＝0,u_x＝0;y＝0,u_y＝0;z＝0,u_z＝0;y＝l,u_y＝l*0.001  (4b)

第3组：x＝0,u_x＝0;y＝0,u_y＝0;z＝0,u_z＝0;z＝h,u_z＝h*0.001  (4c)

第4组：x＝0,u_x＝u_y＝u_z＝0;x＝w,u_z＝w*0.001  (4d)

第5组：x＝0,u_x＝u_y＝u_z＝0;x＝w,u_y＝w*0.001  (4e)

第6组：y＝0,u_x＝u_y＝u_z＝0;y＝l,u_y＝l*0.001  (4f)

步骤7：计算出复合材料等效弹性参数，依次计算出步骤6所述边界条件下的单元阵列的应力分布和应变分布，然后采用方程（5）计算出平均应力

和平均应变

并写成公式（6）所示的向量形式：

${\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{\mathrm{ij}}=\frac{1}{V}{\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&sigma;}}_{\mathrm{ij}}\mathrm{dv},{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{\mathrm{ij}}=\frac{1}{V}{\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&epsiv;}}_{\mathrm{ij}}\mathrm{dv}---\left(5\right)$

${\left[\begin{array}{cccccc}{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{11}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{22}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{33}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{23}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{13}& {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{12}\end{array}\right]}^T,{\left[\begin{array}{cccccc}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{11}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{22}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{33}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{23}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{13}& {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{12}\end{array}\right]}^T---\left(6\right)$

将六组边界条件的计算结果填入矩阵（7），即可得到复合材料的等效刚度矩阵

对

求逆后得到复合材料的等效柔度矩阵

最后由等式(8)可计算出复合材料的等效弹性参数

$\stackrel{\&OverBar;}{S}=\left[s_{\mathrm{ij}}\right],i,j=1~6$

E₁₁＝1/s₁₁,E₂₂＝1/s₂₂,E₃₃＝1/s₃₃

μ₁₂＝-E₁*s₁₂,μ₁₃＝-E₁*s₁₃,μ₂₃＝-E₂*s₂₃

G₂₃＝1/s₄₄,G₁₃＝1/s₅₅,G₁₂＝1/s₆₆      (8)。

2.如权利要求1所述的纤维增强复合材料等效弹性参数的计算方法，其特征在于：所述步骤3中如果没有开展独立的实验来确定colorM和colorF，则通过在color[]数组中寻找最大灰度值colorMax和最小灰度值colorMin，认为colorM=colorMax，colorV=colorMin，在材料内部没有其他杂质时，打开一个断层图片，锁定一个纤维束区域，然后寻找纤维束区域的最大值和最小值，然后取其平均值就是colorF。

Images