二维三轴编织复合材料预压单胞模型建立及其弹性规律数值预测张芳芳; 段永川【期刊名称】《《纺织学报》》【年(卷),期】2019(040)010【总页数】7页(P85-91)【关键词】体素网格; 刚度预报; 参数化单胞; 二维三轴编织; 复合材料【作者】张芳芳; 段永川【作者单位】燕山大学机械工程学院河北秦皇岛066004; 燕山大学先进锻压成形技术与科学教育部重点实验室河北秦皇岛066004【正文语种】中文【中图分类】TB332二维三轴编织复合材料由3个方向的纤维束在平面内相互编织而成，这样在不同角度下可兼顾材料的横向力学行为，有效克服传统层合板力学性能差的缺点。
由于编织复合材料中纤维束相互交织，因此，编织纤维束的截面形状在不同位置会发生畸变，构造实际的纤维束几何模型已十分困难，如果再对纤维束几何模型离散网格构造有限元模型更加困难。
目前多数学者都是基于一定的模型简化来构造编织复合材料的几何和有限元网格模型。
Tsai等[1]采用弹簧单元建立了简化的单胞网格模型，并对弹性力学参数进行了预测。
Miravete等[2]通过沿编织方向增加4条边界线，建立了简化的“米”字形单胞模型，虽降低了构造模型的难度，但也降低了复合材料弹性、损伤等力学性能的预测精度。
Benzley等[3]通过模拟对比得出，六面体单元具有更多的自由度和抗畸变能力，在力学性能预测时能得到更高的精度[4-5]。
Kim等[6]采用一种体素单元法建立了复合材料模型，该方法在表征复合材料两相材料的交界面时精度降低，但这种方法随着分辨率的提高，可消除这个问题。
Gao等[7]假设纱线横截面为矩形，建立出二维三轴编织预制件的动态几何模型。
严雪等[8]基于传统有限元建模方法建立了考虑纤维束弯曲扭转和空间交错特征的有限元模型，并预测了其弹性性能规律。
Kier等[9]假设纤维束横截面为椭圆形，编织纤维束路径为正弦曲线，利用CAD软件生动再现了二维三轴编织复合材料的几何模型。
张芳芳等[10-11]基于直线的纤维束路径建立了三维编织复合材料模型，并对其进行了损伤分析，然而当纤维束路径波动较大时预测误差会逐渐增大。
目前学者大都利用简化模型对材料的微观形态进行表征，但复合材料纤维束在固化时会受到一定的挤压变形，空间形态改变会对材料的弹性及损伤等性能造成影响。
基于此，本文提出基于力学原理考虑空间挤压的模型快速建立方法。
利用变形后纤维束空间分布及形态再形成体素单胞模型，基于该模型预测了二维三轴编织复合材料的弹性规律。
1 细观几何模型1.1 编织结构二维三轴编织复合材料的编织结构如图1所示。
复合材料由3个方向的纤维束组成，其中与y轴平行的为直纤维束，另外2个方向为编织纤维束，他们沿着轴向纤维束相互交织，在z轴方向上看，轴向纤维束位于2个编织纤维束中间，在厚度方向上轴向纤维束位于材料中间。
图中白色框为一个周期单元，称为一个单胞，其宽度为W，高度为L，编织角指轴向纤维束与编织向纤维束所夹的锐角(见图1中α)，单胞内编织向纤维束长Lb。
在图中所示坐标系下，y轴与复合材料轴向纤维束的轴线平行，x轴与复合材料轴向纤维束的轴向垂直。
z轴与复合材料的厚度方向平行。
图1 二维三轴编织复合材料的结构示意图Fig.1 Structural representation of 2-D triaxial braided composite将纤维束中心线定义为空间函数，选用三次样条曲线描述纤维束路径，其中第i段样条曲线pi(t)的参数方程一般形式为pi(t)=Ai+Bit+Cit2+Dit3(1)式中：t为该样条曲线的参数，0≤t≤1, i=0,1,…,n-2，Ai、Bi、Ci、Di分别为第i 段样条曲线的幂次项系数向量。
该方程是在单胞内定义，在单胞两端面上要满足几何连续性，在样条曲线的起始点处应满足的连续条件为p″0(0)=p″n-2(1)(2)为验证样条曲线建立的正确性，在空间取了5个坐标点，每个坐标点的坐标分别为(0,0,0)，(0.5,0.7,0.4)，(1,0.6,0.0)，(1.5,0.7,-0.4)和(2,1.2,0)。
在边界约束条件下生成的周期样条曲线如图2所示。
可以看出，样条线光滑连续，可以满足实际要求。
图2 周期性空间样条曲线Fig.2 Periodic spline curves参考文献[12]中的纤维束横截面拓扑形状，将纤维束的横截面假设为透镜形，如图3所示。
透镜形纤维束横截面由半径分别为r1和r2的2个圆相交而得，2个圆各自偏离图中坐标系横轴的距离是o1和o2。
参数r1、r2、o1和o2可根据透镜的宽度w、高度h和偏移距离d计算得到。
图3 透镜形纤维束截面Fig.3 Lenticular fiber bundle section透镜形横截面的参数方程为：(3)(4)(5)(6)(7)1.2 编织结构参数根据图1所示几何关系，单胞内编织向纤维长Lb为(8)式中：L为单胞高度，mm；α为编织角，(°)。
定义轴向纤维束横截面面积为Aa及圆弧段所对圆心角为αa，其计算公式分别为：(9)(10)式中：ra为圆弧段的半径，mm；w、h分别为透镜形截面的宽度和高度，mm。
同理，可计算出编织向纤维束横截面面积Ab和圆心角αb。
从图1可以看出，在一个矩形单胞模型中含有2根轴向纤维束和4根编织向纤维束，则编织向纤维束和轴向纤维束的体积分别为：Va=2AaLa(11)Vb=4AbLb(12)(13)式中：Va、Vb分别为单胞中轴线纤维束和编织向纤维束的体积，mm3；Aa、Ab 分别为轴线纤维束和编织向纤维束的横截面面积，mm2；La、Lb分别为单胞中轴线纤维束和编织向纤维束的长度，mm；Vf为单胞中的纤维体积含量，%；L、W、H分别为单胞模型的长度、宽度和高度，mm；φ为纱线填充因子，本文认为编织向纤维束与轴向纤维束的纱线填充因子相同。
2 体素法建立有限元模型2.1 单元材料主方向计算纤维束在空间的走向可以用其中心线在空间的走向进行描述，由于纤维束路径的波动导致纤维束各个部位有着不同的材料主方向。
应用有限元方法对实体模型进行离散后，为准确地描述各单元的材料主方向，需要计算单元形心处的材料主方向，如图4所示。
令第i段曲线中某单元的形心点为P，过该点做垂直于纤维束中心线的平面，计算出该平面与纤维束中心线的交点，交点处的参数坐标为tp，该交点沿纤维束走向的切线向量P′i(tp)计算方程为：P′i(tp)·(P-Pi(tp))=0(14)图4 纤维方向求解Fig.4 Fiber orientation solution依据右手法则，确定单元材料坐标系与总体坐标系的变换方向和角度。
本文采用ANSYS软件对有限元模型进行建立和计算，在该软件内采用LOCAL命令为每个单元建立局部坐标系，该坐标系采用2个相对转动角度θxy和θxz定义，如图5所示。
假设一个单元的纤维方向为L(x,y,z)，L(x,y,z)与图中坐标系的z轴构成一个平面，其法线方向为D(x,y,z)。
图5 材料主方向转角定义Fig.5 Definition of material main direction rotation angle图5中x、y和z轴与全局笛卡尔坐标系的三轴平行，其中θ为y轴与向量D所形成的锐角，β为向量L与z轴所成锐角，2个相对转角θxy和θxz的计算公式为：(15)(16)2.2 挤压力学模型的建立为考虑纤维束间的挤压扭曲影响，首先建立二维三轴编织预制件模型如图6(a)所示。
因预制件未与基质黏接固化，无法构造固化后预制件的周期单胞，对预制件进行预压时，采用材料的真实边界进行约束，一个中心单胞利用其领域内8个单胞进行自然边界的约束。
在预制件有限元模型中纤维上下表面分别加入接触单元，在上下2层分别利用平面进行一定的预压缩，将预制件模型厚度压缩至目标单胞模型厚度，从而得到纤维束呈挤压接触状态下的几何形态。
开发周期单胞提取程序，该程序基于FORTRAN语言编写，具有与ANSYS软件兼容的进出接口，首先将变形后的单元节点信息导出，给出切取单胞的边界盒坐标，在边界面上程序可自动细分单元，完成周期单胞单元提取。
为便于后期处理，提取前对不同方位的纤维束进行集合划分，单元切割程序可保持单元的集合和单元坐标系的继承能力，多个单元可能由一个父单元切割再分而来，这几个单元就要保持与父单元同样的集合和单元坐标系方位。
切割出的增强相单胞网格模型，如图6(b)所示。
切割出的增强相单胞模型单元坐标系方位显示如图6(c)所示。
图6 二维三轴编织预制件网格提取Fig.6 Precast mesh extraction of 2-D triaxial braided. (a) Precast model of 2-D triaxial braided composite；(b) Enhanced phase unit cell mesh；(c) Enhanced phase cell mesh with unit coordinate2.3 体素网格模型的建立应用上述切割算法已经获取了变形后纤维束单胞模型，本文利用体素模型法建立变形后的单胞模型，在该模型区域内填充正方体素单元。
根据已提取的预制件单胞单元来对所填充的正方体素单元进行集合创建与单元坐标系的继承。
利用预制件单胞单元在每根纤维束表面生成表面单元，然后判断体素单元与表面单元的相对位置，当体素单元的形心坐标在表面单元的包络区域内时，将该体素单元纳入该集合内，其他纤维束计算方法相似，最后生成的体素模型如图7所示。
图7 采用体素法建立的复合材料单胞网格模型Fig.7 Single cell mesh model of composite materials by voxel method.(a)Enhanced phase mesh；(b)Matrix mesh；(c) Mesh after superimposed该模型采用全正方六面体网格建立，边界非常容易满足边界节点一一对应的关系。
3 周期性边界条件复合材料内部细观结构周期重复，在外载荷作用下其应力应变场同样满足周期重复性，为降低计算量，选取一个代表单元进行力学行为预测，在单胞2个对应的边界面建立周期对应的耦合约束方程[13]，约束表达式如下:(i,j,k=x,y,z)(17)式中：j+表示沿正方向；j-表示沿负方向；ui为边界面上节点的i方向位移，为全场平均应变；xk为边界面上节点的k方向坐标；为沿xj方向两边界面上对应节点在k方向的坐标差。
图8 编织角和纤维体积含量对二维三轴编织复合材料横向和轴向模量的影响Fig.8 Influence of volume content and angle of weaving on transverse and axial modulus of 2-D triaxial braided composite materials. (a) Elastic modulus of transverse; (b) Elastic modulus of axial; (c) Elastic modulus of thick; (d) Shear modulus of transverse; (e) Shear modulus of axial; (f) Shear modulus of thick; (g) Poisson ratio of transverse; (h) Poissio raton of axial; (i) Poisson ratio of thick4 结果分析以文献[12]中所提供的实验数据进行对比验证，模型中所采用的材料组分力学性能参数如表1所示。
