航天器热防护系统普遍采用复合材料，迄今再入大气层速度最高的人造航天器Stardust的热防护系统就采用了酚醛浸渍碳烧蚀(Phenolic impregnated carbon ablator，PICA)，一种纤维增强复合材料作为其隔热材料，其以低密度、耐高温、抗热震、高隔热性能等特点成功保障了采样返回器在再入大气层的严酷气动加热环境下舱内温度维持在设计范围内^[1]。

对于这类以隔热为目的功能性复合材料，准确掌握复合材料在宏观上表现的等效性能，尤其是作为衡量材料隔热性能的导热系数对于准确预测航天器热防护系统再入过程的烧蚀热响应、降低设计冗余度及指导材料研制过程性能设计都具有重要意义。

试验测定复合材料的导热系数是最直接可信的方法。但实际操作过程中，对高温隔热材料导热系数的测定周期较长、成本高，尤其在高温环境下，材料内部形态发生变化甚至发生化学反应，致使其难度较大。目前通过材料在电弧风洞中的温度测量结果进行反向迭代^[2]也可给出等效导热系数，但温度响应本身是依赖于多变量的，且迭代目标函数的选取也缺乏充分的理论依据，因此采用该导热系数计算的烧蚀量与飞行试验之间往往存在较大的差异。

鉴于此，从理论上研究复合材料的导热特性，具有重要的现实及理论价值。国内外开展了大量研究，建立了Maxwell-Eucken、Pilling等诸多模型^[3]，Eucken公式反映了复合物中各组分的导热系数和体积分数对复合物等效导热系数的影响。陈则韶等^[4]根据最小热阻力法则和比等效导热系数相等法则，研究了形状、排布的影响。程远贵^[5]采用等效热阻法建立了高温条件下耐火纤维的传热模型，并引入分形维数的概念对原等效导热系数进行了修正。F.Gori等^[6]给出烧蚀条件下规则编织体复合材料等效导热系数的预测方法。这些方法往往基于某种规则几何体，在某种极限假设下对传热问题进行简化，从而得到等效导热系数依赖于某种因素的解析表达。此外，将有限元、有限体积等数值传热技术拓展应用于纤维、孔隙尺度揭示材料微观结构与宏观材料参数之间的关系的研究也是近年的研究热点^[7-8]。张海峰等^[9]、谢涛等^[10]分别采用热阻网络法和有限体积法对随机结构的颗粒弥散型复合材料的导热系数进行了理论预测。无论是热阻网络方法还是数值模拟方法，首先都依赖于代表性单元(Representative volume element，RVE)的选取以简化分析，对于封闭式掺入或者纤维规则排布的情况下，代表性单元的选取确定较明确，但针对PICA类材料，其纤维排布随机分布、且纤维长度L远大于纤维直径的情况下，如何选取单元体才可以确保其具有代表性是值得讨论的问题。同时，数值模拟微观导热问题时，网格剖分需要同时兼顾纤维半径尺度与纤维径向尺度的刻画，针对L$\gg$r的纤维体这一要求会导致计算量过大而无法接受。

针对现有研究在解决以PICA材料为代表的纤维随机分布复合材料时所面临上述难点，本文在分析PICA类材料的微观形貌特点的基础上，提出了采用随机结构生成方法模拟纤维的随机分布，并以纤维间相交自然形成的连接点为网格点，通过迭代方法确定稳态条件下的温度场分布，从而根据热流通量对等效导热系数进行预测，并以此方法研究揭示影响纤维增强复合材料等效导热系数的关键几何参数。

1.   随机结构复合材料等效导热系数理论模型

1.1   基本物理假设

PICA材料由高热导率的纤维增强体和低热导率的不定型基质复合而成，其中碳纤维直径约为14~16 um，长度约为1600 um^[11]，其中纤维体积分数10vol%，随机分布。图1为PICA材料的微观形貌图^[1]。研究复合材料等效导热系数，关键在于确定给定条件下的内部能量通量。在温度梯度的作用下，能量以总体热阻力最小的方式进行传递，可以将其视为两部分：即以增强相纤维作为传递通道的能量输运、以聚合物基质相为传递通道的能量输运。通常情况下，聚合物基质在复合材料内部以连续相的形式存在，影响能量传递的几何因素主要是截面积的变化。基质材料的导热系数较纤维增强体一般低两个数量级，保证了材料良好的隔热性能；其各向同性使其能量传递并不具有方向选择性。与此相对，作为增强相的纤维体具有较好的导热能力，为能量输运提供了短路通道，但随机分布的纤维是介于连续相与离散相之间的，限制其导热能力的主要因素在于纤维空间搭接关系所形成的能量通道的有效性。另一方面，恰是纤维的空间排布及纤维自身在径向、横向的导热系数差异使复合材料整体的导热系数表现出各向异性的特点。因此复合材料导热系数的复杂性主要由增强相的纤维引起。

图  1  酚醛浸渍碳烧蚀(PICA)材料微观形貌^[1]

Figure  1.  Micrograph of phenolic impregnated carbon ablator (PICA) material^[1]

下载: 全尺寸图片 幻灯片

根据该型材料的微观图像^[1]，且考虑到基质材料导热比纤维导热能力低两个数量级^[3]，假设：(a) 纤维为给定长度的线段，所有纤维半径为均一值r；(b) 能量主要通过纤维编织体形成的路径进行传递；(c) 在纤维编织体的基础上引入基质材料对导热系数的贡献修正。

1.2   随机生成结构方法

在二维空间中，单一某根纤维的空间位置可以通过起点坐标$({x_0},{y_0})$、空间方向角$\theta$及线段长度$l$唯一指定。相应地在三维空间中可以通过纤维起始点坐标$({x_0},{y_0},{z_0})$，线段在空间的方向角$(\theta ,\varphi )$及线段长度$l$予以确定，如图2所示。当纤维增强体是规则编织体时，这些参数及纤维之间的交点坐标是可以预先指定的，但对于类PICA材料，采用随机函数生成这些参数，并采用随机结构生成法(Random generate-growth method，RGGM)^[12]逐根添加纤维模拟材料内纤维的随机排布。由于纤维间的位置关系是未知的，添加新纤维时需要依次求解判断其与当前其余纤维的位置关系。若存在交点，则记录下该点坐标、形成该交点的两根纤维编号及该节点在纤维上的相临节点编号。其中考虑到相临节点可能会由于其后添加的纤维发生变化，引入链表结构记录节点间的相互关系。如此循环添加下一根纤维直至纤维的体积分数$\left({\text{π}} r\displaystyle\sum\limits_i {{l_i}}/V \right)$或单位体积的纤维根数n/V等统计参数符合试验观测结果，以生成的线段簇描述随机分布的纤维体，从而建立一整套的节点信息及连接关系。

图  2  纤维生长模型示意图

Figure  2.  Sketch of fibers generation model

l—Fiber length; θ, φ—Directional angle

下载: 全尺寸图片 幻灯片

1.3   空间节点及相互联接关系的确定

二维空间，可以通过求解两根纤维的参数化方程确定交点及其是否在纤维长度范围内。三维空间中，两根求空间线段的位置关系要复杂些，采用数值方法很难严格定义交点的存在。因此采用空间直线的距离予以确定，假定$d < 2r$，且公垂线在位于两根纤维长度范围内，则认为存在交点，设置其两根纤维上公垂线的中点为用于计算的空间节点。

具体求解过程如下：

设有空间线段AB，其起点A点坐标为$({x_1},{y_1},$${z_1})$，终点B坐标为$({x_2},{y_2},{z_2})$，则该线段上任意的一点P的坐标$(X,Y,Z)$可以表示为

另一空间线段CD，其起点C点坐标为$({x_3},{y_3},$${z_3})$，终点D坐标为$({x_4},{y_4},{z_4})$，该线段段任意一点Q的坐标$(U,V,W)$可以表示为

当参数$0 \leqslant s \leqslant 1$时，P为线段AB上的点；参数$0 \leqslant t \leqslant 1$时，Q为线段CD上的点。

P、Q两点之间的距离为

求线段AB、CD最短距离，即求目标函数$f\left( {s,t} \right)$的最小值：

即令其满足：

若参数$s,t$的值满足$0 \leqslant s \leqslant {\rm{1}}$及$0 \leqslant t \leqslant {\rm{1}}$，说明P、Q落在两条线段内，两条线段的距离为PQ的长度。若两线段间距离小于半径和，且P、Q落在两条线段内，则生成该点坐标，并记录下该点在两条线段内的相邻节点。从而建立整个纤维体的空间连接关系。

建立两个动态数组，分别存储描述纤维几何位置信息及节点位置信息；同时针对每根纤维建立一个链表依次存储该纤维上的节点编号。对于任意节点，首先确定该节点由哪两根纤维相交而成，进而在该纤维对应链表中确定与该节点相邻的上下节点，并给出相互间距离。因此，内部任一节点仅与沿两根纤维与四个节点发生能量交换。

1.4   导热模型及有效导热系数的计算

针对任一节点建立能量平衡方程，通过迭代得到稳定条件下节点温度的分布情况：

式中：$\rho$表示密度；${c_{\rm{p}}}$表示等压比热；$T$表示温度；${Q_{i,{\rm{in}}}}$表示由相邻节点进入的热流量。离散后得到

则下一时刻的温度$T_0'$为

其中：$T_0'$表示代求点下一时刻温度；$\Delta t$表示推进时间步长；$\Delta V$表示控制体体积；${S_i}$纤维横截面积；${\lambda _{\rm{f}}}$表示纤维导热系数；${T_i}$临近第$i$点温度；${l_i}$临近第$i$点与该点的距离。

在稳态下，当纤维排布、体积分数及边界条件给定情况下，材料内部的温度分布仅表现为上述变量的函数，而与纤维的材料特性无关，也即与$\left( {{\lambda _{\rm{f}}}{S_i}\Delta t} \right)/\left( {\rho {c_{\rm{p}}}\Delta V} \right)$无关，为使迭代过程满足基本物理直观、保证最终结果收敛，可以选取$\Delta t$令当前点下一时刻温度介于周围节点温度与当前时刻温度形成的区间内，并取前后时刻的温度相对变化量作为收敛判据^[9]，即

温度相对变化小于预先给定值$\varepsilon$，则判定其收敛，否则进行下一轮迭代，直至达到稳态温度场。为加快收敛可设置节点初始温度与其坐标相关。

根据得到的稳态温度场，可以给出任意截面上的热流通量，进而依据等效导热系数的定义确定等效导热系数。

取上边界固定温度为${T_{{\rm{b1}}}}$，下边界固定温度为${T_{{\rm{b2}}}}$，通过试件的热流通量为

其中：${T_{i,{\rm{next}}}}$表示材料内部临近边界$i$点的节点；$r$表示纤维半径。

上下界面的距离为${d_0}$，$S$为总面积，则热流通量：

结合式(10)、式(11)，等效导热系数${\lambda _{{\rm{eff}}}}$可如下给出：

可见，纤维增强体的等效导热系数是纤维径向导热的线性函数，中间的面积比反映了该纤维增强体中的纤维体积分数，最后一项是无量纲温差与无量纲节点间距给出的温度梯度的平均，可见与试件高度无关、所加载的温度，仅纤维体的几何拓扑结构的函数。

1.5   基质材料对等效导热系数的贡献

对于复合材料中基质材料的存在状态，假设两种情况：(1) 基质均匀填充纤维体构成的孔隙间。对于这种情况，假设基质与纤维体整体构成并联模型，则复合材料整体的等效导热系数$\lambda _{{\rm{eff}}}^{\rm{'}}$可以由下式给出：

其中，${\alpha _{\rm{m}}}$表示在该方向的基质面积比。

(2) 基质围绕在纤维外形成包裹层。对于这种情况，可以选取圆柱体作为代表性单元，采用并联模型或有限元方法计算径向的等效纤维导热系数$\lambda _{\rm{f}}'$，

其中，$R$表示含包裹层的半径。进而依据等效纤维导热系数求出材料整体的等效导热系数。

2.   算法验证

2.1   与理论结果的对比

以规则的二维纤维编织体为例，在纤维相对材料尺度为无限长且为直线情况下，可知某根纤维上的温度分布为关于y的线性函数，且该线性函数的参数由材料上下表面边界温度确定，与纤维的排布角度无关。因此，不同方向角的纤维之间不会存在温度差异，也就不会有能量传递。在这种情况下，整个温度场就是线性的。

图3为边界处的纤维示意图。可根据第一层的交叉点的温度，给出边界上通过每根纤维的热流：

其中：${T_1}$、${T_2}$表示材料上下面温度；$H$表示材料厚度。可见通过每根纤维的热流量仅与其空间角度有关，而与位置无关。

图  3  边界处的纤维示意图

Figure  3.  Sketch of fibers at boundary

h_i—Dietance to the boundary; θ₁, θ₂—Boundary angle

下载: 全尺寸图片 幻灯片

通过第一簇纤维的热流通量为

其中：${d_1}$表示第一簇纤维间距；$L/{d_1}$表示单位宽度上第一簇纤维的根数。

通过第二簇纤维束的热流通量与此类似，则总的热流通量可表示为

根据等效导热系数定义，可以给出其计算式：

即纤维编织体等效导热系数取决于纤维的本征导热系数、纤维截面积、方向角度及疏密程度。

假设取一平板，令其上下壁面为等温边界，设上面板温度为500 K，下面板温度为300 K，上下间距为0.5 mm，左右壁面为绝热边界。取纤维直径为 7.5 um，纤维为无限长，在边界处截断。其排布为规则排布的两簇，在第一种情况下方向角度分别40°、130°，即正交，同族纤维间的间距0.04 mm；第二种情况下方向角分别为30°、135°，第一族纤维间距为0.04 mm，第二族纤维间距为0.02 mm。设置各节点的初始温度为300 K，取迭代收敛判据为1.0×10⁻³ K。

图4为一般规则排布下纤维编织体内的温度场分布。可见温度在厚度方向均表现为线性分布。依据温度给出的规则正交和一般规则排布下的无量纲导热系数$\dfrac{{{\lambda _{{\rm{eff}}}}}}{{{\lambda _{\rm{f}}}}}$分别为0.1863、0.2343。相对理论解0.1875、0.2344的偏差分别为0.64%和0.006%。迭代得到的温度场与理论解一致，表明该模型用于预测实际碳纤维编织体等效热导率是有效的。其中，等效导热系数与理论解之间的偏差实际体现了在纤维稀疏情况下试件尺寸相对于纤维间距不充分大时纤维编织体表观参数的真实物理波动。

图  4  一般规则排布情况下的温度场

Figure  4.  Temperature distribution in general regular arrangement situation

下载: 全尺寸图片 幻灯片

2.2   与实验结果的对比

将本文方法与碳纤维随机增强环氧树脂复合材料的实验结果^[13]进行对比，其中碳纤维的导热系数取为200 W/(m·K)，基体材料的导热系数取为0.2 W/(m·K)。图5为质量分数5wt%情况下的纤维和温度分布情况。图6为不同纤维质量分数下的理论预测结果与实验结果的对比情况。

图  5  质量分数5%情况下的碳纤维随机增强环氧树脂复合材料中纤维和温度分布情况

Figure  5.  Fibers and temperature distribution for carbon fiber reinforced epoxy resin (5wt%)

下载: 全尺寸图片 幻灯片

图  6  碳纤维随机增强环氧树脂复合材料数值模拟结果与实验结果的比较

Figure  6.  Comparison of simulation and experimental results for carbon fiber reinforced epoxy resin

下载: 全尺寸图片 幻灯片

可见，数值预测结果与实验测量结果所给的等效导热系数的数值大致相当，其随纤维质量分数的变化趋势一致。但在个别状态下与实验结果的偏差明显。分析原因可能是由于本文模型也未考虑接触热阻的影响，且随机生成的结构未能完全刻画材料的真实模型状态。对于如何表征随机结构、选取两种随机结构等效的评判参数是值得深入探讨的另一问题。

3.   结果与讨论

3.1   问题描述

将上述方法应用于增强相纤维随机分布的编织体结构,预测随机结构材料的等效导热系数并分析影响等效导热的相关因素。计算状态设为二维条件，取试件宽度为10 mm，厚度为5 mm，其中内部纤维的直径为15 um，设置上表面温度为500 K，下表面温度为300 K，两侧为绝热边界条件。假设试件内部的纤维排布完全随机，取温度场的收敛判据取为1.0×10⁻³ K。首先给出纤维长度为1.6 mm、3.2 mm、4.8 mm、6.4 mm、8.0 mm、9.6 mm情况下，编织体的温度分布及其对应的等效导热系数、纤维体积分数，梳理纤维特征长度与等效导热系数间的关系，其次给出纤维数目为255、300、350、400、450、500、550情况下纤维编织体的温度分布及对应的等效导热系数、体积分数，分析纤维数目与等效导热系数间的关系，最后研究了固定纤维体积分数条件下，不同纤维直径、长度的影响。此外，将节点间纤维段在垂直方向的平均长度定义为有效长度，即$l = \dfrac{1}{{{n'}}}\displaystyle\sum\limits_i {{l_i}\left| {\sin {\theta _i}} \right|}$，其在一定程度上刻画了纤维体的联通性，进而尝试归纳其与等效导热之间的关系。

3.2   随机结构复合材料等效导热系数与纤维特征长度间关系

图7为相同纤维数目(n=255)不同纤维长度下的温度分布。受纤维随机分布的影响，温度分布呈现高度的不规则性，但整体上随着纤维特征长度的增加，温度场更趋于均匀化。虽然温度场的变化并不明显，但等效导热对纤维特征长度是高度敏感的，可导致成倍的变化。一方面这是由于在纤维根数一定的情况下，纤维特征长度的增加，必然引起纤维体积分数的增加；另一方面纤维特征长度增加，导致纤维间的连接增强，一定程度上引起热阻网络由串联连接向并联连接的转化，更有效地形成能量通路，整体上降低了纤维编织体的热阻，提高了等效导热系数。

图  7  不同纤维特征长度下酚醛浸渍碳烧蚀材料的温度分布

Figure  7.  Temperature distribution of phenolic impregnated carbon ablator under different intrinsic fiber length ((a) 1.6 mm; (b) 3.2 mm; (c) 4.8 mm; (d) 6.4 mm; (e) 8.0 mm; (f) 9.6 mm)

下载: 全尺寸图片 幻灯片

表1为不同纤维长度下的体积分数、有效长度及无量纲等效导热系数。在PICA材料的典型纤维体积分数范围内(≈0.1)，无量纲等效导热系数约为0.02，即等效导热系数是纤维径向本征导热系数的2%左右，而基质的本征导热系数在纤维本征导热系数的1%以下，说明通过纤维的热传导是能量传递中的主要形式，可以用纤维体的等效导热系数进行复合材料的整体导热系数的估算。此外，等效导热系数与特征长度呈非线性关系，当特征长度大于材料尺寸时，特征长度的增加将不再影响体积分数及有效长度，等效导热系数趋于某一定值。即当纤维体的整体几何尺度相对于纤维特征长度较小时，则特征长度不再影响等效导热系数。该结论同样提示在进行材料等效导热系数测量时，必须要求试样尺度选择合理或符合实际应用条件，否则试样尺度会引起测量结果的变化。该定性结论与文献[14]的结果一致。

表  1  不同纤维特征长度下的酚醛浸渍碳烧蚀材料等效导热系数

Table  1.  Effective thermal conductivity of phenolic impregnated carbon ablator under different intrinsic fiber length

Fiber length/mm	Volume fraction	Effective length/mm	Effective thermal conductivity $\dfrac{{{\lambda _{{\rm{eff}}}}}}{{{\lambda _{\rm{f}}}}}$
1.60	1.01×10−1	1.17×10−1	2.10×10−2
3.20	1.71×10−1	7.52×10−2	6.40×10−2
4.80	2.13×10−1	6.16×10−2	8.64×10−2
6.40	2.30×10−1	5.65×10−2	8.92×10−2
8.00	2.39×10−1	5.37×10−2	9.05×10−2
9.60	2.42×10−1	5.34×10−2	9.09×10−2

下载: 导出CSV 

| 显示表格

3.3   随机结构复合材料等效导热系数与纤维根数的关系

图8为相同纤维长度(l=1.6 mm)、不同纤维数目下的温度分布。表2为不同纤维数目下相关结果的对比情况。纤维体积分数、有效长度与纤维数目近似线性关系，前者为正相关，后者为负相关。随着纤维数目的增加，体积分数升高、有效长度降低两种因素共同作用都促使等效导热系数增强。纤维数目由250增加到500时，纤维的无量纲等效导热系数从1.96×10⁻²增加到6.49×10⁻²，表现出非线性关系。

图  8  不同纤维数目下酚醛浸渍碳烧蚀材料的温度分布

Figure  8.  Temperature distribution of phenolic impregnated carbon ablator under different fiber number ((a) 250; (b) 300; (c) 350; (d) 400; (e) 450; (f) 500)

下载: 全尺寸图片 幻灯片

表  2  不同纤维数目下的等效导热系数

Table  2.  Effective thermal conductivity under different fiber number

Fiber number	Volume fraction	Effective length/mm	Effective thermal conductivity $\dfrac{{{\lambda _{{\rm{eff}}}}}}{{{\lambda _{\rm{f}}}}}$
250	9.90×10−1	1.17×10−1	1.96×10−2
300	1.20×10−1	1.02×10−1	2.93×10−2
350	1.40×10−1	9.15×10−2	3.74×10−2
400	1.59×10−1	8.19×10−2	4.80×10−2
450	1.81×10−1	7.22×10−2	5.59×10−2
500	2.01×10−1	6.65×10−2	6.49×10−2

下载: 导出CSV 

| 显示表格

3.4   随机结构复合材料等效导热系数与有效长度的关系

表3为在相同纤维体积分数(20vol%)，n=500下不同纤维直径下的无量纲等效导热系数。可见在相同体积分数下，等效导热系数存在数量级的变化，纤维体积分数并不是决定纤维体隔热性能的单一变量，以往文献拟合给出的等效导热系数与体积分数的关系^[15]应当只适用某一特定构型的纤维体，不存在通用的拟合公式。随着纤维直径减小，长度增加，纤维间形成的连接节点数增多，有效长度逐渐减小，形成很好的能量通路，等效导热系数增加。反之，当纤维直径增加长度减小时，有效长度逐渐增加，等效导热系数降低，当纤维体不足以形成能量的传递通道时，则本文方法失效，纤维体成为理论上的完全隔热体，等效导热系数趋于零。以当前算例为例，当纤维直径为0.035 mm，长度为0.61 mm时，纤维体的无量纲等效导热系数为2.73×10⁻¹³。理论分析表明采用短切纤维有利于材料隔热性能的增强。

表  3  相同体积分数(20vol%)下不同纤维尺寸的等效导热系数

Table  3.  Effective thermal conductivity under different fiber size with same volume fraction (20vol%)

Case	Radius×length/mm	Effective length/mm	Effective thermal conductivity $\dfrac{{{\lambda _{{\rm{eff}}}}}}{{{\lambda _{\rm{f}}}}}$
Case1	0.010×2.65	4.46×10−2	7.83×10−2
Case2	0.015×1.60	6.65×10−2	6.49×10−2
Case3	0.020×1.13	8.58×10−2	4.15×10−2
Case4	0.025×0.88	1.01×10−1	1.93×10−2
Case5	0.030×0.72	1.13×10−1	2.14×10−3
Case6	0.035×0.61	1.20×10−1	2.73×10−13

下载: 导出CSV 

| 显示表格

4.   结 论

针对随机结构的纤维增强复合材料，提出了一种预测等效导热系数的理论方法。该方法物理直观简明，可为随机结构的材料设计及性能预测提供参考。既可将其视为一种基于随机结构的无网格数值方法，也可以将其视为热阻网络法在复杂随机结构的一种拓展。相较于宏观的有限元方法，其回避了代表性单元的选取问题；相较于格子玻尔兹曼方法模拟细观结构中的热量传递时依赖纤维尺度的网格划分，该方法的网格依据纤维排布自然而成，从而降低了对网格尺度的要求。通过与规则情形下的解析结果对比，证明了该方法的有效性，进而将其应用全随机的纤维体，研究了影响等效导热系数的相关因素。

针对随机结构的纤维体，当纤维特征长度超过导热方向的穿透长度时，则不再影响材料内部等效导热系数，也说明复合材料的等效导热系数并非材料固有属性，还依赖于材料试件的尺度与纤维特征长度的相对长短；单位空间的纤维根数与等效导热系数呈非线性正相关关系，但等效导热系数并非体积分数的单变量函数，引入有效长度的概念作为描述纤维体的联通性的补充参数，其有效长度愈短，等效导热系数愈大。