编织复合材料热物理性能快速预测及实验验证

张芳芳，段永川，高安娜

(1.燕山大学 机械工程学院，河北 秦皇岛 066004；2.先进锻压成形技术与科学教育部重点实验室(燕山大学)，河北 秦皇岛 066004)

0 引言

编织复合材料是由两股或多股纤维束在空间按一定规则相互编织缠绕，形成编织预制件后，经树脂浸渍固化而成。成型过程中纤维束保持连续，成型件的整体性能好，并具有比模量大、比强度高以及性能可设计等突出优点[1]，在航天领域得到了广泛应用[2]。但航天环境昼夜温差大，容易使编织复合材料结构件产生一定的变形和内应力，进而影响结构件间的连接匹配和安装精度，因此编织复合材料的热传导和热膨胀性能的预测及设计研究具有重要的意义。

目前许多学者通过有限元方法建立细观模型从而研究其热物理性能[3-4]。程伟等[5]建立了一种单胞模型，称之为简化“米”字型，应用于复合材料导热性能的研究。刘振国等[6]为了预测材料的热学性能，建立了一种更加贴近实际工程的单胞模型，基于该模型对三维四向编织复合材料在热载荷下的热学性能进行了模拟预测。而Soheil[7]从细观尺度出发，运用有限单元法，创建了细观的单胞模型，应用于三维复合材料的等效热膨胀系数的获取。夏彪等[8]基于有限单元方法，选取了三维四向和三维五向两种复合材料，研究了固定角度下的等效热传导系数和等效热膨胀系数。姜黎黎[3]采用螺旋型纤维等效模型研究了热载荷对三维四向编织复合材料拉伸性能的影响。目前学者们采用的模型有两点要求：一是要求纤维束与基体界面处节点共用；二是要求单胞模型相对面上节点一一对应。同时基体是由流动到纤维束相互挤压缝隙中的树脂固化而成，为空间多边形的复杂体(如图1(b)所示)。这样建立单胞模型时很难自动完成网格划分，需要反复调整网格密度。成功的模型往往计算规模较大，如果进一步对复合材料的热物理性能进行优化将更加困难。迫切需要寻找新的方法预测复合材料热物理性能。

本文提出基于耦合思想的编织复合材料热学性能快速预测方法，在该方法中，利用数学约束方程替代纤维束与基体界面处节点共用的网格约束，放松了原模型中的第一条要求。在边界上通过耦合关系转化了周期性边界条件，放松了原模型中的第二条要求。同时整体区域模型采用六面体为主的网格划分，提高了计算精度、降低了计算规模。利用该方法对三维四向编织复合材料热物理性能进行了快速预测，经与试验和传统方法对比，验证了该方法的精度和可行性，同时为提高复合材料优化设计的效率奠定了基础。

图1 复合材料共节点网格模型建立过程

Fig.1 Process of establishing composite co-node mesh model

1 单胞模型的建立

1.1 网格模型的建立

首先采用六面体网格建立编织复合材料中的纤维束模型，利用编制的网格切割程序提取单胞区域内的纤维束模型，当六面体单元经过单胞模型的边界面切割提取时，余下的部分可以采用六面体、三棱柱或者四面体单元重新进行网格划分，根据切割面与单元的位置关系，一共有14种相对位置关系；当三棱柱单元再次被单胞边界面切割提取时，切割面与单元的相对位置关系共有19种；当四面体单元再次被单胞边界面切割提取时，切割面与单元的相对位置关系共有3种。根据这些单元可能的被切割情况，进行单元结构形式的重新构建，并编写了相应程序[9]，在该程序中，通过读入模型的单元节点坐标、单元节点连接关系和要进行切割提取的单胞边界面坐标，该程序便可进行自动的判断和网格重组划分，并最终输出为有限元软件可快速读入的代码数据文件。本文利用该方法进行单胞模型区域内网格模型的切割提取。为了便于观看，此处只显示一个方向的纤维束模型提取，如图2所示。

图2 纤维束切割提取过程示意图

Fig.2 Schematic diagram of fiber bundle extraction process

将切割提取出的增强相模型与整体区域模型叠合，得到复合材料单胞网格模型。如图3所示。

图3 基于耦合思想的单胞网格模型建立过程

Fig.3 Process of establishing composite unit cell mesh model by coupling method

1.2 自由度协调关系的处理

从图3模型建立过程中可以看出，增强相网格和整体区域网格在空间上会存在区域的重叠，因此在叠合区域内节点的自由度会被重复定义，需要协调自由度。本文基于耦合法，对单胞模型重合区域自由度进行协调处理，处理的对象分别是增强相和整体区域这两个模型，对两者形状函数、位移场以及温度场的协调处理分别如式(1)～(3)所示：

(1)

其中，x、y、z为增强相单元节点在整体区域单元中的形状函数；xi、yi、zi和Ni分别为整体区域的单元节点坐标和形状函数值；m为整体区域单元节点个数。

(2)

其中，ui，vi，wi为整体区域单元的节点位移，u、v、w为增强相单元节点位移。

(3)

其中，T和Ti分别为增强相和整体区域的单元节点温度。

1.3 材料热学性能匹配

复合材料由增强相和基体相组成，在建模过程中需要对两相的网格模型分别赋予增强相和基体相材料属性。从图3建模过程中可见，由于整体区域网格与增强相网格重叠区域的基体材料属性的存在，使得单胞网格模型中，产生了附加刚度，故需要对增强相单元刚度进行修正，修正公式为

(4)

2 单胞模型温度场周期性边界条件

六面体单胞模型如图4所示，建立坐标系，分别对它的12条棱和8个顶点进行命名，对应单胞模型的各边及节点。为了便于有限元软件处理计算，将施加在节点的温度载荷与单胞网格模型的某一方向的自由度进行对应并逐一进行替换。即三维空间中三个坐标轴上的节点温度载荷与对应的自由度替换：x轴方向与节点4对应，y轴方向与节点5对应，z轴方向与节点2对应。为了便于温度载荷的施加，令节点1上的温度T1为0，三个坐标方向的平面约束方程为

(5)

对于x=x1和x=x2与y=y1和y=y2面的相交棱边l1～l4上，应分别满足

(6)

对于x=x1和x=x2与z=z1和z=z2面的相交棱边l5～l8上，应分别满足

(7)

对于y=y1和y=y2与z=z1和z=z2面的相交棱边上l9～l12，应分别满足

(8)

对角节点约束方程，变换为

(9)

图4 六面体单胞模型边及节点标号

Fig.4 Signs of hexahedral unit cell model edges and nodes

3 传热及热膨胀系数的计算方法

材料的热传导系数Ki随着单胞模型施加的边界条件的不同而变化，建立边界条件与Ki的函数关系，即可求出不同情况下的Ki，公式如下

(10)

其中，qi为热流输出面在i方向的平均热流密度，ai为单胞i方向的边长，温度差ΔT为用户输入参数。由于未涉及非线性问题，因此在计算材料的等效热物理性能时不需要考虑ΔT的取值大小。

进行温度周期边界条件施加时，只需在单胞模型中的4、5和2节点上施加温度载荷即可。将这三个节点的温度载荷分别代入到式(10)，可得单胞模型在x,y、z三个方向的等效热传导系数

(11)

单胞模型在x、y、z三个主轴方向的等效热膨胀系数为

(12)

其中，u41、u51、u21分别表示节点4、5、2在x方向的位移。

纤维与树脂材料相互混合形成的束状的混合物定义为纤维束。纤维束的横向导热性能参数，可通过Hashin的上限估计公式[10]求得，即

(13)

纤维束纵向导热性能参数用混合定律公式[10]来确定，即

KL=Kf 1Vf+Km(1-Vf),

(14)

利用Schapery(SH)公式计算纤维束的热膨胀系数[11]

(15)

α2=(1+μm)αm(1-Vf)+

(1+μf12)αf1Vf-αf1μf12,

(16)

其中，Kf1，Kf2，Km分别为纤维束轴向、横向以及基体的导热系数，Vf为纤维束中的纤维体积百分比，上标+表示上限值，μm,μf12，分别为基体的泊松比和纤维的主泊松比，αf 1,αf 2分别为纤维的轴向和横向热膨胀系数，αm为基体的热膨胀系数。

4 有限元分析结果的验证与对比

为验证本文采用耦合法预测复合材料热物理性能的正确性，本文采用三维四向编织复合材料试件，进行热膨胀系数测试试验。

三维四向编织复合材料的组成成分是T300纤维和环氧树脂，热物理性能参数见表1[5]。

采用全桥方式连接应变片，利用因瓦合金(又名不胀钢)进行温度补偿。贴片位置示意图如图5所示，其中应变片1、2、3、4组成全桥电路，用于测量试件的轴向热应变；应变片5、6、7、8组成全桥电路,用于测量复合材料试件的横向热应变。

表1 组分材料热物理性能参数

Tab.1 Thermo-physical properties of constituents

材料热传导系数W·K-1·m-1热膨胀系数10-6K-1K11K22α11α22T300纤维8.001.00-0.33.1环氧树脂0.180.1831.731.7

图5 应变片粘贴位置示意图

Fig.5 Schematic diagram of the position of strain gauges

复合材料试件与因瓦合金组成的全桥电路如图6所示。

图6 应变片组成的全桥电路图

Fig.6 Full-bridge circuit diagram composed of strain gauges

图6中C代表在复合材料试件上粘贴的应变片，I代表在因瓦合金上粘贴的应变片，UI为输入端电压，UO为测量端电压。当环境温度改变时，应变片的阻值会因为温度的变化而发生改变，从而引起测量端电压的改变，若将因瓦合金上的应变片与复合材料试件上的应变片并联，当环境温度改变时，粘贴在因瓦合金上的应变片与粘贴在复合材料试件上的应变片的阻值发生相同的改变，因此测量端电压不变化，消除了由于应变片本身阻值改变而引起的附加应变。

复合材料热物理性能试验现场如图7所示，进行复合材料试件的等效热膨胀系数测量时，采用德国INSPEKT Table100 kN电子万能高温试验机，利用德维创动态数据采集系统对应变进行记录测量。

图7 热物理性能测试平台

Fig.7 Thermal physical performance test platform

首先对因瓦合金的温度补偿性能进行测试，在单独的因瓦合金上采用全桥方式粘贴应变片，将恒温箱中的温度从室温逐渐升高到42 ℃左右，温度保持一段时间后，继续升温到62 ℃左右，记录整个过程中因瓦合金的应变曲线，发现上下波动量小于1×10-6ε，验证了因瓦合金的温度补偿性能。

将贴好应变片的复合材料试件与因瓦合金一同放入测试环境中，通过温度控制箱，将保温箱中温度从室温逐渐升高到42 ℃，保持此温度不变，待试件变形稳定一段时间后，再将温度逐渐升高到62 ℃，记录沿复合材料试件轴向和横向的应变。试验重复3次，数据稳定性较好，选取其中一组试验数据曲线如图8所示。

分别取42 ℃和62 ℃温度保持期时复合材料试件的平均应变，利用式(12)计算复合材料试件的热膨胀系数。测得复合材料试件横向热膨胀系数为5.01×10-6K-1，轴向热膨胀系数为-1.09×10-6K-1。

采用耦合法建立与试验试件相同编织角即编织角为19.2°的三维四向编织复合材料模型，数值预测的横向热膨胀系数为4.76×10-6K-1，与试验结果相差5.25%；数值预测的轴向热膨胀系数为-1.01×10-6K-1，与试验结果相差7.92%。

图8 三维四向编织复合材料热膨胀性能测试曲线

Fig.8 Test curve of thermal expansion performance of three-dimensional four-directional braided composites

将采用耦合法分析预测的材料热传导和热膨胀系数，与采用传统方法预测的结果进行对比，如表2所示。

表2 基于耦合法与传统方法预测的热物理性能结果对比

Tab.2 Comparison of thermo-physical properties predicted using cupling method and traditional FEM

传统有限元法(收敛结果)耦合法粗网格中等密度网格细网格(收敛结果)相对误差节点数132758121112900954348—单元数714407108492710352379—Kx/(W·K-1·m-1)0.820.790.810.820.0%Ky/(W·K-1·m-1)0.850.820.840.850.0%Kz/(W·K-1·m-1)4.094.074.094.100.2%αx/(10-6K-1)4.734.714.744.760.6%αy/(10-6K-1)4.734.714.744.760.6%αz/(10-6K-1)-1.00-0.87-1.00-1.011.0%

其中，采用耦合法建立的单胞模型，随着网格密度的逐渐增大，模型对材料热物理性能的预测结果没有明显差距，可认为获得收敛解。两种方法预测结果之间的相对误差较小，且本文建立的模型中耦合法的计算规模仅为传统方法的40%，进一步验证了耦合法用于复合材料热物理性能分析的快速性。

另外，对表2中的数据进行横观比较，不难发现，所示的热物理性能具有横观各向同性，分析认为材料的截面属性和形状特征都会影响材料物理性能参数的分布特性及数值大小，其横观各向同性特征本质上是由于材料的横截面具有对称性。

4 结论

1)提出了预测编织复合材料热物理性能的耦合方法，模型采用六面体为主的网格划分，建立了增强相单元与整体区域单元网格的温度场自由度协调关系、增强相与整体区域叠合后的热学性能匹配方法。

2)构造了单胞模型的表面、棱边上的节点与角节点之间的耦合关系，实现了仅在角节点处施加温度载荷来完成整个单胞模型边界条件的快速建立。

3)建立了复合材料试件热膨胀性能测试及温度补偿实验系统，在测试温度范围内该系统的应变偏差小于1×10-6ε，验证了该系统的可靠性；经数值预测结果与试验结果的对比，验证了数值预测方法的可行性和正确性，为编织复合材料热物理性能的快速优化奠定了基础。