陈素芳， 谭志勇， 姜 东， 董萼良， 费庆国(.东南大学 工程力学系,南京 0096;.北京临近空间飞行器系统工程研究所,北京 00076;3.南京林业大学 机械电子工程学院,南京 0037;.东南大学 机械工程学院,南京 89)

陶瓷基纤维增强复合材料具有热防护性能好、比强度/比刚度高等优点，广泛应用于航空、航天领域，进而面临极端严酷的高温服役环境。高温环境下结构动态特性分析至关重要，是复合材料热结构动响应计算、振动控制和优化设计的前提。

在热结构动态特性分析方面，杨和振等[1]以复合材料层合板为研究对象，研究了在不同温度条件下复合材料层合板的动力响应，表明复合材料层合板结构的模态参数受到温度变化的显著影响。Spivey[2]对C/SiC热结构操纵面进行了热载荷状态下结构的试验、仿真振动特性分析。Cheng等[3]以加筋板为研究对象，结合材料的热物理性能和力学性能进行热载荷状态下结构的振动特性分析和非均匀温度场对热结构模态参数的影响研究。Wu和Wang等[4]分别利用试验和仿真的方法研究了板在高温环境下热模态参数。通过热模态分析可深入了解温度梯度、材料参数对结构固有频率的影响，为热结构设计提供参考。

为了提高计算效率，复合材料结构的动态分析往往采用等效建模：在常温条件下，仅需要考虑材料的等效弹性性能[5-6]，但高温环境下，材料的弹性模量随温度变化，温度梯度将导致结构内部产生热应力。因此，热结构动态分析时，应着重关注复合材料的弹性和热膨胀性能。热膨胀系数作为表征材料热稳定性的重要参数，对高温工作环境下复合材料的动力学性能影响较大，热应力和热模态分析等都依赖于材料的热膨胀性能。

复合材料的热膨胀性能可以根据其材料组分性能，纤维含量等参数，进行理论分析[7]。20世纪七八十年代，Schapery等[8-9]基于等效弹性矩阵得出了各向同性、多相材料热膨胀系数的上下限。刘书田等[10]利用均匀化的方法，研究了含有均匀分布球形空洞的空心铝材料和单向纤维增强复合材料的热膨胀行为，建立了组分材料性能、体积分数与热膨胀系数之间的关系。梁军等[11]针对定向分布微裂纹的纤维增强复合材料，利用Eshelby理论和Mori-Tanaka方法，预报了复合材料的有效热膨胀系数。尽管人们一直不断探索材料热物理性能预测的理论模型和计算方法，但是由于热物理性能对材料的组成、晶体结构及制备工艺都很敏感，加之复合材料本身结构复杂，所以试验测定仍是最为直接有效的方法。Tognana等[12]利用光学热膨胀仪，测量了不同体积含量的颗粒增强复合材料的热膨胀系数。Pan等[13]采用分离式霍普金斯压杆，结合加热装置的系统，测试3D编织复合材料23 ℃～210 ℃范围内材料的热物理性能。

随着数值计算被引入复合材料热膨胀系数等效研究，有限元法逐步成为继理论分析、试验测定之后的另一种有效手段。燕瑛等[14]基于细观力学有限元法预报了复合材料热膨胀系数，并与理论方法结果进行对比。Karadeniz等[15]基于有限元的方法预报了纤维增强复合材料的横向、纵向热膨胀系数。有限元法的有效性依赖于计算模型的准确建立，且模型要尽可能地模拟材料内部的真实结构。为此，聂荣华等[16]依据二维编织复合材料的细观结构特点，提出一种预测该材料的面内热膨胀系数的单胞模型，充分考虑到编织结构复合材料中的纤维束弯曲和CVI工艺制备陶瓷基复合材料产生的孔洞对热膨胀系数的影响。卢子兴等[17]在三维全五向编织复合材料细观结构模型的基础上，建立了单胞参数化有限元模型，计算编织复合材料的热膨胀系数。Zhang等[18]建立了三维编织复合材料结构不同空间位置的三种单胞模型，施加边界条件，计算三种情况下的等效弹性参数。Lu等[19]在研究纤维、基体连接界面对编织复合材料热物理性能的影响时，提出一种三维有限元模型，将纱线考虑成单向纤维增强复合材料并采用代表性体积单元计算其热物理性能。

本文以单向纤维增强复合材料为对象，通过细观力学理论推导，建立一种同时预测复合材料弹性参数和热膨胀系数的方法。并将等效后的均质材料属性和等效前各组分材料属性分别建立相应的等效模型和精细化有限元模型，通过计算结构的热动态特性分析，验证结构热参数等效的可行性。在此基础上，进一步研究不同纤维排布方式、不同单胞模型对参数等效预测结果的影响。

1 理论基础

1.1 热弹性力学理论

热弹性力学中，假设材料的温度和变形可分别计算，不考虑耦合关系；材料始终处于弹性范围内，不考虑因温度产生的非弹性变化；材料变形很小不考虑他们的二阶微分项。因此，胡克定律在均匀常值温度下的本构方程

σij=cijklεkl

(1)

如果考虑温度变化的影响，而不考虑温度和变形之间的耦合关系，则材料的本构方程为

σij=cijklεkl+βij(T-T0)=cijklεkl+βijΔT

(2)

式中：βij是在应变为零时测量的热模量，是一个对称张量；T0为参考温度，即弹性模量cijkl测定时的温度，cijkl也称为刚度系数。

由式(2)求应变，可得

εij=sijklσkl+αijΔT

(3)

式中：αij为各向异性材料的热膨胀系数；sijkl为柔度系数。将式(2)展开为矩阵形式则为

(4)

上式21个弹性系数Cij和6个热系数βij，表示各向异性体全不对称材料。按照简化张量标记，见表1。

表1 简化张量标记Tab.1 Simplified marker of tensor

这样，式(2)、式(3)简化为

σi=cijεj+βi(T-T0)=cijεj+βiΔT

(5)

εi=sijσj+αiΔT

(6)

由上述应力应变式(5)、(6)运算，可得刚度系数、柔度系数、热模量和热膨胀系数之间的关系

cijskl=δik

(7)

βi=-cijαj

(8)

αi=-sijβj

(9)

式中:δ为克罗尼克符号；β为热模量；α为热膨胀系数；c,s分别表示刚度矩阵和柔度矩阵。因此，按照方程式(9)的关系，开展有限元法等效预测过程。

1.2 参数预报步骤

步骤1等效刚度预测

在基准温度T=T0下，温差ΔT=0。本构关系式满足：σ=cijεj，即

(10)

当ε1=1；ε2=ε3=γ4=γ5=γ6=0时，应力矩阵等于刚度矩阵的第一列[c11c21c31c41c51c61]T；依次每施加一个单位正(切)应变所得的应力列阵等于刚度矩阵中相应的一列。

表2和表3分别给出了等效分析的六种工况及六种工况下相应的应变条件。从表3可以看出，单位正(切)应变转化为相应的位移边界条件，可以实现单位应变载荷的施加。如工况一，加载x方向单位正应变(εx=1)，转化为位移边界条件只需在x=0的面上施加固定约束，x=a的面上施加特定位移a，其余四个面x方向自由，y,z方向约束。

表2 等效分析的六种工况Tab.2 Six cases for equivalent analysis

表3 边界条件Tab.3 Boundary conditions

由两种材料构成的细观不均匀结构，宏观上将其等效为均质等效体。等效体的应力等于细观结构的平均应力；等效体的应变等于细观结构的平均应变

(11)

(12)

不考虑温度变化时，则平均应力与平均应变之间满足本构关系

(13)

(14)

(15)

对于正交各向异性材料而言，其刚度系数矩阵D满足

(16)

对应的柔度系数矩阵S

(17)

因此，等效弹性参数及等效主泊松比(v12、v23)和副泊松比(v31)为

(18)

步骤2等效热模量预测

在T0+ΔT温度下(ΔT≠0)，使单胞整体结构不发生任何变形(即ε=0)，则式(2)满足：σ=Dα·ΔT=β·ΔT，即

(19)

采用周期性位移边界条件，施加单位温度差值，经稳态热分析得非均匀温度场。并将此非均匀场作为结构分析的热载荷条件进行静力分析计算单胞的热应力场。因此，每一个单位温度差值下，计算得到的应力列矩阵等于热模量β。

表4、表5分别给出了热传导分析时加载的三种工况和对应的边界条件。求解的热膨胀系数是在T=T0时的值，如果材料组分(纤维、基体)的性能在温度变化区间[T0，T1]内不是常量，可将原区间细分为若干子区间，认为每个子区间内材料的性能参数为常数。考虑上述因素，因此在施加边界条件时，选取的温差值a、b、c在整个温差ΔT范围内，材料参数都保持不变，可以不必细分温度区间。

编程获取各个单元应力、单元体积，按式(20)体积平均，求解平均热应力

(20)

表4 热传导分析的三种工况Tab.4 Three cases for heat transfer analysis

表5 给定的边界条件Tab.5 The given boundary conditions

由步骤1可知，单位应变场的加载不易实现，转化为位移场加载。同理，单位长度的温差在未经热分析之前更不易确定和加载。因而，首先对结构两对面施加任意的温度差值，经稳态热分析，可求得单元的温度梯度，然后将平均热应力归一化，即得单位温度梯度下结构的热应力，即等效热模量。

(21)

式中：βe表示等效热模量；Δ表示单位长度的温度梯度的平均值。

步骤3等效热膨胀系数

由第1步得到的等效刚度矩阵，求逆，得等效柔度矩阵；由第2步求得等效热模量，按式(9)求得到等效热膨胀系数。其中，等效模型的密度根据纤维、基体各自的体积分数求解计算。

图1给出了单向纤维增强复合材料热膨胀系数预测流程图。首先对单胞模型进行刚度预测，得到复合材料的等效刚度矩阵，然后对该单胞模型进行热传导分析，获得非均匀分布温度场，将其作为温度载荷，加载于结构进行结构分析，获得结构的平均热应力，由式(21)进而得等效热模量。最后，由式(9)求得结构的等效热膨胀系数。

图1 热膨胀系数等效预测流程图Fig.1 Prediction flow chart for equivalent thermal expansion coefficient

2 单胞模型

通过对组成材料的性能的了解，再由给定纤维增强结构的几何参数，应用解析的方法来预报复合材料的宏观力学性能，就是细观力学的研究方法。细观力学理论假设了纤维理想黏结，均匀分布在整个基体中，且材料中的孔隙和气泡很小。复合材料的纱线是由许多纤维与基体合成的单向纤维增强复合材料柱。按照纤维在基体中排布方式的不同，取正方形排列、六边形排列，如图2(a)、(b)所示，图2(c)正方形排列纤维束对应的代表性体积单元，图2(d)六边排列纤维束对应的代表性体积单元,只要对图2(c)、(d)进行平移变换、镜面对称就能形成整个宏观结构材料。

(a)纤维束正方形排列(b)纤维束六边形排列(c)正方形排列对应的胞元(d)六边形排列对应的胞元

图2 纤维排列分布及选取的单胞

Fig.2 Fiber distribution and selection of the unit cell

图3给出了单胞有限元模型。单胞II和IV分别为单胞Ⅰ和Ⅲ代表性体积单元(RVE)。缩减的RVE细观结构模型，一般会极大增加周期性边界条件的施加难度。但由于等效预测给定边界条件在单胞和代表性体积单元上的作用面有所不同，因此，取单胞I、III作为研究对象的同时，代表性体积单元(单胞II、IV)也作为研究对象。网格划分时保证网格在空间上满足周期性的要求。设单胞中基体为各项同性，纤维为横观各向同性材料。纤维、基体的弹性性能、热物理性能的材料参数如表6、表7所示。

(a)单胞I(b)单胞Ⅱ(c)单胞Ⅲ(d)单胞Ⅳ

图3 不同单胞的有限元模型

Fig.3 Different finite element model of unit cell

图4分别给出了四种单胞施加6种基于单位应变载荷(εx，εy，εz，γxy，γyz，γzx)的边界条件得到的应力分布云图。由图4(a)、(b)可以看出，由于单胞Ⅱ是单胞Ⅰ的1/4代表性体积单元，在施加周期性位移边界条件后静力分析过程中，两者的应力云图分布基本一致：基体的应力明显高于纱线所受应力，由于基体的弹性性能高于纱线，按照纤维、基体刚度大小分配原则，导致纤维所受应力小于基体应力。同时，从工况一对比来看，单胞Ⅰ、Ⅱ的应力细微有所区别，这主要由于单胞Ⅱ受边界约束的六个表面相比于单胞Ⅰ的六个边界约束面同时分布有纤维和基体，进而使两种材料所受应力重新分配。

表6 纤维和基体的弹性性能(C纤维/SiC基体)Tab.6 Elastic properties for the fiber and matrix

表7 纤维和基体的热物理性能Tab.7 Thermophysical properties for the fiber and matrix

(a) 单胞I应力云图

(b) 单胞II应力云图

(c) 单胞III应力云图

(d) 单胞IV应力云图图4 不同单胞在六种单位应变工况下的应力云图Fig.4 Stress distribution of unit cells in six units strain conditions

按式(11)～(16)，计算得等效刚度矩阵。表8、表 9分别给出了四种单胞的等效刚度矩阵。其中，“/”前为单胞I、III的等效刚度矩阵，“/”后为代表性体积单元(单胞II、IV)的等效刚度矩阵。对比表8、表 9的计算结果可以看出，同一类型纤维排列方式下，两种单胞取法的计算结果非常接近。仅在C44、C55、C66中存在微小差异。由等效刚度矩阵求逆，按式(18)可求得等效弹性参数。

表8纤维束正方形排列对应单胞Ⅰ/Ⅱ的等效刚度矩阵

Tab.8EquivalentstiffnessmatrixofunitcellI/IIunderfiberbundlesquarearrangement

GPa

表9纤维束六边形排列对应单胞III/IV的等效刚度矩阵

Tab.9EquivalentstiffnessmatrixofunitcellI/IIunderfiberbundlehexagonalarray

GPa

由工程弹性常数互等关系，可求得副泊松比。从表10中可以看出，两种排列方式对应的单胞等效后都表现横观各向同性，但由于纤维束六边形排列的不对称性，导致其等效结果稍稍偏离横观各向同性的特点；纤维排列方式对横向弹性模量影响较大，正方形排列的横向弹性模量大于六边形的横向弹性模量；而轴向弹性模量和纤维排列方式无关；同一种纤维排布方式的单胞，其单胞尺寸的大小对剪切模量影响很大。这主要由代表性体积单元(单胞Ⅱ、Ⅳ)被施加周期性边界条件的作用面和单胞(I、III)有所不同造成的。

由平均热应力计算得等效热模量，按式(20)、(21)及式(9)求得细观结构等效热膨胀系数,如表11所示。

3 结果比较

为验证等效参数的准确性、合理性，将纤维增强复合材料进行精细化有限元建模和等效建模。计算结构非线性热模态，边界条件两边简支，整体结构施加100 ℃均匀温度场，对比精细化模型和等效模型的热模态频率及振型。如图5所示，首先建立复合材料的梁模型(尺寸:800 mm×10 mm×10 mm)，通过对比梁的精细化有限元模型(正方形排布精细化模型180万个单元，六边形排布精细化模型208万个单元)和均质等效模型(8万个单元)，验证x方向等效的准确性。然后，建立复合材料板(尺寸:2.5 mm×100 mm×200 mm)，通过对比板的精细化有限元模型(正方形排布精细化模型11.25万个单元，六边形排布精细化模型为130万个单元)与其均质等效模型(6万个单元)，验证y,z方向参数等效预测的准确性。梁、板结构精细模型，在随着单元网格数目逐渐加密的同时，有限元的解基本保持不变，且单元网格数目已达两百万个单元，解的收敛性保持较好；等效模型网格数目虽大大下降，但由于等效模型材料属性已预测为均质的横观各向同性材料，目前的单元个数已经可以满足精度要求。而且，从等效前、后动态特性对比结果可以看出，两者有限元解的误差很小，满足精度要求。

表10不同单胞对应的等效弹性参数

Tab.10EquivalentelasticparametersofunitcellsGPa

表11不同单胞的等效热膨胀系数(×10-6)

Tab.11Theequivalentcoefficientofthermalexpansionofunitcells(×10-6) GPa

梁

板图5 梁、板几何模型Fig.5 Geometry model of beam and plate8

图6和图7分别给出了梁、板等效前、后精细有限元模型和等效模型。(a)为正方形排列纤维束精细化梁模型，(b)是六边形排布的精细化梁模型，纤维体积分数均为40%。图(c)为单胞等效后的均质等效模型，材料参数值如表10和表11所示。

1⁃1(a)正方形排布精细化模型1⁃1(b)六边形排布精细化模型1⁃1(c)等效模型

图6 梁等效前、后有限元模型

Fig.6 Finite element model of beam between equivalent before and after

2⁃2(a)正方形排布精细化模型2⁃2(b)六边形排布精细化模型2⁃2(c)等效模型

图7 板等效前、后有限元模型

Fig.7 Finite element model of plate between equivalent before and after

表12和表13给出了梁模型等效前、后结构的热模态频率和对应的误差。其中，正方形排布纤维束精细模型的热模态频率与均质等效模型的热模态频率最大误差为0.4%；六边形排布纤维束最大误差0.7%，均在合理范围内。可见，由本文提出的预测复合材料热膨胀系数的方法是合理有效的，从而验证了本文预测复合材料热膨胀系数及工程弹性参数方法的正确性。

表12 梁等效前后热模态频率对比Tab.12 Comparison of beam thermal modal frequencybetween equivalent before and after

表13 梁等效前后模态频率对比Tab.13 Comparison of beam thermal modal frequencybetween equivalent before and after

表14和表15给出了复合材料板在等效前、后结构的热模态及其对应误差。由表中结果可以看出，相比于梁等效前后对比结果，板在等效前、后存在一定的误差。说明在整个等效预测过程中，x方向预测精度较高；y,z方向的预测虽有一定误差，但也在合理范围之内。

由上述四个表对比分析可知，不同类型单胞之间，单胞I、II比单胞III、IV预测准确，说明纤维排布方式对参数等效预测存在影响；同一纤维排布方式下，代表性体积单元越大，预测数据误差更小，即单胞I比单胞II预测结果合理，单胞III比单胞IV 预测结果精准。这主要是由于缩减的代表性体积单元极大增加了周期性边界条件的施加难度。除非在极其复杂的非线性问题中，否则，代表性体积单元并不作为等效预测分析的首选。

表14 板等效前后模态频率对比Tab.14 Comparison of plate thermal modal frequencybetween equivalent before and after

表15 板等效前后模态频率对比Tab.15 Comparison of plate thermal modal frequencybetween equivalent before and after

4 结 论

(1)依据细观力学基础，推导了复合材料热膨胀系数的公式，建立一种简洁的等效热、力学性能的有限元预测方法。针对预测结果，分别进行精细化建模和等效建模，对比等效前、后结构的热模态及振型，验证了本文方法的可行性及有效性。

(2)本文所提出的等效热膨胀系数公式和有限元预测分析方法，不但预测结构的热膨胀系数，而且可以一并预测复合材料弹性参数。

(3)不同纤维排布方式下，正方形排布纤维束预测结果比六边形排布纤维束预测精度高；同一纤维排布方式下，单胞的预测结果较缩减后的RVE预测结果精度高。

(4)该预测热膨胀系数的方法，可以预测复合材料结构在室温下热膨胀系数的值，还可以预测不同高温下的复合材料的热膨胀系数。为其他新型构型的复合材料弹性、热膨胀性能的预报提供参考。

参 考 文 献

[1] 杨和振, 李华军. 温度变化下复合材料层合板的试验模态分析[J]. 复合材料学报, 2008, 25(2): 149-155.

YANG Hezhen, LI Huajun. Experimental modal analysis of the composite laminates with temperature variation[J].Acta Materiae Compositae Sinica, 2008, 25(2): 149-155.

[2] SPIVEY N D. High-temperature modal survey of a hot-structure control surface[C]∥ICAS 27th International Congress of the Aeronautical Sciences. Nice, France: ICAS, 2010.

[3] CHENG H, LI H, ZHANG W, et al. Effects of radiation heating on modal characteristics of panel structures[J]. Journal of Spacecraft and Rockets, 2015, 52(4): 1228-1235.

[4] WU D, WANG Y, SHANG L, et al. Experimental and computational investigations of thermal modal parameters for a plate-structure under 1 200 °C high temperature environment[J]. Measurement, 2016, 94: 80-91.

[5] JIANG Dong, ZHANG Dahai, FEI Qingguo, et al. An approach on identification of equivalent properties of honeycomb core using experimental modal data[J]. Finite Elements in Analysis and Design, 2014, 90: 84-92.

[6] JIANG Dong, LI Yanbin, FEI Qingguo, et al. Prediction of uncertain elastic parameters of braided composites[J]. Composite Structures, 2015, 126, 123-131.

[7] 黄争鸣. 复合材料细观力学引论[M].北京:科学出版社, 2004: 69-73.

[8] SCHAPERY R A. Thermal expansion coefficients of composite materials based on energy principles[J]. Journal of Composite Materials, 1968, 2(3): 380-404.

[9] ROSEN B W, HASHIN Z. Effective thermal expansion coefficients and specific heats of composite materials[J]. International Journal of Engineering Science, 1970, 8(2): 157-173.

[10] 刘书田, 程耿东. 基于均匀化理论的复合材料热膨胀系数预测方法[J]. 大连理工大学学报, 1995, 35(4): 451-457.

LIU Shutian, CHENG Gengdong. Homogenization-based method for predicting thermal expansion coefficients of composite materials[J]. Journal of Dalian University of Technology, 1995, 35(4): 451-457.

[11] 梁军, 陈晓峰. 含缺陷纤维增强复合材料热膨胀系数预报[J]. 哈尔滨工业大学学报, 1997, 29(3): 36-38.

LIANG Jun, CHEN Xiaofeng. Thermal expansion coefficients of fiber composite materials containing matrix microcracks[J] Journal of Harbin Institute of Technology University, 1997, 29(3): 36-38.

[12] TOGNANA S, SALGUEIRO W, SOMOZA A, et al. Influence of the filler content on the thermal expansion behavior of an epoxy matrix particulate composite[J]. Materials Science and Engineering: B, 2009, 157(1): 26-31.

[13] PAN Z, GU B, SUN B. Longitudinal compressive behavior of 3D braided composite under various temperatures and strain rates[J]. Applied Physics A, 2014, 118(4): 1315-1337.

[14] 燕瑛，李剑峰. 复合材料热膨胀性能的细观分析模型与预报[J]. 北京航空航天大学学报, 2013, 39(8): 1069-1085.

YAN Ying, LI Jianfeng. Prediction of thermal expansion properties for composites by micromechanical analysis model[J] Journal of Beijing University of Aeronautics and Astronautics, 2013, 39(8): 1069-1085.

[15] KARADENIZ Z H，KUMLUTAS D. A numerical study on the coefficients of thermal expansion of fiber reinforced composite materials[J]. Composite Structures, 2007, 78(1): 1-10.

[16] 聂荣华, 矫桂琼, 王波. 二维编织C/Sic复合材料的热膨胀系数预测[J]. 复合材料学报, 2008, 25(2): 109-114.

NIE Ronghua, JIAO Guiqiong, WANG Bo. Prediction on coefficient of thermal conductivity for 2D braided C/SiC composites[J].Acta Materiae Compositae Sinica, 2008, 25(2): 109-114.

[17] 卢子兴, 王成禹, 夏彪. 三维全五向编织复合材料弹性性能及热物理性能的有限元分析[J].复合材料学报, 2013(3): 160-167.

LU Zixing, WANG Chengyu, XIA Biao. Finite element analysis of elastic property and thermo-physical properties of three-dimensional and full five-directional braided composites[J]. Acta Materiae Compositae Sinica, 2013(3): 160-167.

[18] ZHANG C, XU X. Finite element analysis of 3D braided composites based on three unit-cells models[J]. Composite Structures, 2013, 98: 130-142.

[19] LU Z, WANG C, XIA B, et al. Effect of interfacial properties on the thermophysical properties of 3D braided composites: 3D multiscale finite element study[J]. Polymer Composites, 2014, 35(9): 1690-1700.