复合材料广泛运用于能源工业领域，通常由两相或者多相材料混合而成。由于复合材料几何结构复杂，在对其进行传热分析时存在一定的困难，因此通常采用等效均匀化方法将复合材料看作具有等效导热系数的均匀材料进行分析。常规的复合材料不含内热源，通过在边界上施加热流，通过保证等效前后的热流一致来研究等效导热系数，比较经典的模型包括热阻串并联模型^[1]、Maxwell模型^[2]等。在电容器、航天器、太阳能热存储、核反应堆、铀矿的开采、煤的堆放等^[3-6]特殊场景中，也存在大量复杂的复合材料传热问题。由于含有内热源，复合材料内的热流分布规律不同于仅有外加热流的情况，传统基于外加热流得到的等效导热系数模型可能不再适用于含内热源复合材料的传热分析^[7]，需要深入研究内热源对热流分布的影响，重新选取等效前后的守恒量来研究等效导热系数^[8]。文献[9-11]通过数值的手段计算了含内热源复合材料的等效导热系数，将复合材料的最高温度或平均温度作为等效前后的守恒量，计算得到了预测最高温度或平均温度的等效导热系数结果，说明了含内热源时的等效导热系数不同于无内热源时的传统等效导热系数模型。采用等效导热系数预测最高温度可用于研究核反应堆中的热点，与安全性相关；采用等效导热系数预测平均温度可用于计算核反应截面，用于反应堆物理热工耦合计算和设计分析。数值分析方法可针对特定工况开展含内热源复合材料的等效导热系数研究，对于影响等效导热系数的规律研究还比较少，理论模型研究还不充分。

为建立复杂的核燃料^[12-13]等效导热系数模型，本文从几何结构较为简单的无限大复合平板入手，采用格林函数法^[14]的思想，推导得到了发热材料均匀分布时的复合平板等效导热系数理论模型，可用于预测含内热源复合平板的平均温度。通过与不含内热源时的热阻串联模型进行对比，说明内热源存在时等效导热系数的特殊性及影响因素。

1 发热板均匀分布的复合平板

本文研究的无限大复合平板如图 1所示黑色部分是发热弥散板，白色部分是基体板。

无限大复合平板由两相材料组成，其中，弥散相为同一种材料，分布在图 1中的黑色区域，宽度一致、热源大小相同，另一相材料为基体相。无限大复合平板可在厚度方向建立一维坐标，坐标原点取为平板中心处，板厚为δ，左端和右端边界分别为-δ/2和δ/2。弥散在平板中的发热板有n块，宽度均为x₀，发热板中心的坐标从左至右依次为x₁，x₂，…，x_n，每个发热板的热源功率线密度均为ϕ。为简化起见，发热板的导热系数λ₁和基体的导热系数λ₂均认为与温度无关。

发热板均匀分布在复合平板内，即：相邻发热板的中心间距均相等，且所有发热板间基体材料的宽度相等，任意选取2块发热板为例，如图 2所示。

相邻2块发热板之间基体材料的宽度Δ为：

$ \varDelta=\frac{\delta-n x_{0}}{n+1} $

(1)

相邻2块发热板的中心间距为：

$ \varDelta+x_{0}=\frac{\delta-n x_{0}}{n+1}+x_{0}=\frac{\delta+x_{0}}{n+1} $

(2)

2 等效导热系数的推导

本节主要推导复合平板的等效导热系数，关注不含内热源情况下的复合平板在外加热流条件下的等效导热系数模型，并讨论含内热源复合平板的等效导热系数。

2.1 无内热源

针对图 1所示的无限大复合平板模型，当弥散板不发热时即为常规的复合平板，将复合平板中的不同材料等效成一种材料，同时保证等效前后的边界条件不变，取复合平板的等效导热系数为λ_e，等效过程如图 3所示阴影部分是不发热的弥散板。

在复合平板上施加热流q，边界温度分别为t₁和t₂(t₂>t₁)，由于热流q在截面上处处相等，基于傅里叶导热定律以及热阻的概念，q为：

$ q=\frac{t_{2}-t_{1}}{\frac{n x_{0}}{\lambda_{1}}+\frac{\delta-n x_{0}}{\lambda_{2}}}=\frac{t_{2}-t_{1}}{\frac{\delta}{\lambda_{e}}} $

(3)

式(3)化简可得：

$ \frac{\lambda_{e}}{\lambda_{2}}=\frac{1}{\frac{n x_{0}}{\delta} \cdot \frac{\lambda_{2}}{\lambda_{1}}+\frac{\delta-n x_{0}}{\delta}} $

(4)

令每个弥散板的无量纲尺寸μ为：

$ \mu=\frac{x_{0}}{\delta} $

(5)

则复合材料中弥散材料的填充率φ为：

$ \varphi=\frac{n x_{0}}{\delta}=n \mu $

(6)

填充率是指弥散相材料的宽度占总复合平板宽度的份额。同时，令无量纲导热系数κ为：

$ \kappa=\frac{\lambda_{2}}{\lambda_{1}} $

(7)

则式(4)可化简为：

$ \frac{\lambda_{e}}{\lambda_{2}}=\frac{1}{1+\varphi(\kappa-1)}=\frac{1}{1+n \mu(\kappa-1)} $

(8)

等效导热系数模型与经典的热阻串联模型一致，等效导热系数的大小与材料的填充率及其材料导热系数有关。

2.2 有内热源

同样对于图 1所示的复合平板，如果弥散板含内热源，几何均匀化的方法一致，将不同材料等效成一种材料，同时，热源也需要进行均匀化，可按照总热源不变的原则将热源均匀分布到复合平板上，同时保证等效前后的边界条件不变，简单起见，边界温度均设为t₀，含内热源复合平板的等效过程如图 4所示。

由于复合平板内部含分布式内热源，不同位置处的热流不再处处相等，与图 3中不含内热源的情况不同，因此，不能采用热阻法直接得到等效导热系数，需要建立能量守恒方程重新推导。等效过程中，选取平均温度作为守恒量。如图 4所示，等效后平板的等效导热系数λ_{e, s}，线功率密度减小为φϕ，边界温度维持不变为t₀，控制方程为：

$ \frac{\mathrm{d}}{\mathrm{d} x}\left(\lambda_{e, s} \cdot \frac{\mathrm{d} t}{\mathrm{~d} x}\right)+\varphi \phi=0,-\frac{\delta}{2}<x<\frac{\delta}{2} $

(9)

容易求得其温度分布为：

$ t=-\frac{\varphi \phi}{2 \lambda_{e, s}} x^{2}+\frac{\varphi \phi \delta^{2}}{8 \lambda_{e, s}}+t_{0} $

(10)

对整个平板上的温度分布进行几何平均，得到其平均温度为：

$ t_{\text {avg }}=\frac{\int_{-\frac{\delta}{2}}^{\frac{\delta}{2}}\left(-\frac{\varphi \phi}{2 \lambda_{e, s}} x^{2}+\frac{\varphi \phi \delta^{2}}{8 \lambda_{e, s}}+t_{0}\right) \mathrm{d} x}{\delta}=\frac{\varphi \phi \delta^{2}}{12 \lambda_{e, s}}+t_{0} $

(11)

化简得到可预测复合平板平均温度的等效导热系数为：

$ \lambda_{e, s}=\frac{\varphi \phi \delta^{2}}{12\left(t_{\text {avg }}-t_{0}\right)}=\frac{n x_{0} \phi \delta}{12\left(t_{\text {avg }}-t_{0}\right)} $

(12)

由式(12)可知，获得复合平板的平均温度，即可得到其等效导热系数。

针对图 1所示复合平板模型，可以分别建立基体材料和弥散发热材料各个部分的控制方程，从左至右依次为：

$ \begin{aligned} &\frac{\mathrm{d}}{\mathrm{d} x}\left(\lambda_{2} \frac{\mathrm{d} t}{\mathrm{~d} x}\right)=0,-\frac{\delta}{2} \leqslant x \leqslant x_{1}-\frac{x_{0}}{2} \\ &\left\{\begin{array}{c} \frac{\mathrm{d}}{\mathrm{d} x}\left(\lambda_{1} \frac{\mathrm{d} t}{\mathrm{~d} x}\right)+\phi=0, x_{i}-\frac{x_{0}}{2} \leqslant x \leqslant x_{i}+\frac{x_{0}}{2} \\ \frac{\mathrm{d}}{\mathrm{d} x}\left(\lambda_{2} \frac{\mathrm{d} t}{\mathrm{~d} x}\right)=0, x_{i}+\frac{x_{0}}{2} \leqslant x \leqslant x_{i+1}-\frac{x_{0}}{2} \end{array}\right. \\ &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \vdots\\ &\frac{\mathrm{d}}{\mathrm{d} x}\left(\lambda_{1} \frac{\mathrm{d} t}{\mathrm{~d} x}\right)+\phi=0, x_{n}-\frac{x_{0}}{2} \leqslant x \leqslant x_{n}+\frac{x_{0}}{2} \\ &\frac{\mathrm{d}}{\mathrm{d} x}\left(\lambda_{2} \frac{\mathrm{d} t}{\mathrm{~d} x}\right)=0, x_{n}+\frac{x_{0}}{2} \leqslant x \leqslant \frac{\delta}{2} \\ &i=1,2,3, \cdots, n-1 \end{aligned} $

(13)

根据边界条件，以及发热弥散相和基体相材料界面处的热流和温度连续条件，可以逐一求解方程并联立解得整个复合平板的温度分布，通过沿x方向上的几何平均，即可得到复合平板的平均温度。

本文从另外一种思路来求解。根据格林函数法，含分布式内热源复合材料的温度场等于每个热源以及边界条件单独作用时温度场的叠加。以含3块发热板的复合平板为例进行说明，如图 5所示，阴影部分是不发热时的弥散板，黑色部分是发热的弥散板。因此，求解整个含内热源复合平板温度分布的问题可以转化为先求解某一块弥散板带内热源时的温度分布，然后再将所有弥散板单独发热时的温度分布进行叠加。

对于含多块弥散材料的复合平板，当仅有一块弥散板发热时，可分解为3个部分，对于该发热板两侧的部分，可以分别看作是无内热源的复合平板，仍然可以利用热阻串联模型分析，对于发热板可单独求解含内热源的控制方程，其等效过程和思路如图 6所示。该复合平板的温度分布示意如图 7所示。

假设发热板中心坐标为x_i，则发热板左侧和右侧的弥散板个数分别为n_l和n_r，共有n-1块，数量多少与左右两侧的厚度成正比，具体为：

$ n_{l}+n_{r}=n-1 $

(14)

$ \frac{n_{l}}{n_{r}}=\frac{\frac{\delta}{2}+x_{i}-\frac{2 \delta+(1-n) x_{0}}{2(1+n)}}{\frac{\delta}{2}-x_{i}-\frac{2 \delta+(1-n) x_{0}}{2(1+n)}} $

(15)

求解得到：

$ n_{l}=(n-1) \cdot \frac{\frac{\delta}{2}+x_{i}-\frac{2 \delta+(1-n) x_{0}}{2(1+n)}}{\delta-\frac{2 \delta+(1-n) x_{0}}{1+n}} $

(16)

$ n_{r}=(n-1) \cdot \frac{\frac{\delta}{2}-x_{i}-\frac{2 \delta+(1-n) x_{0}}{2(1+n)}}{\delta-\frac{2 \delta+(1-n) x_{0}}{1+n}} $

(17)

同时，发热板左右两部分的弥散板填充率为：

$ \varphi_{l}=\frac{n_{l} x_{0}}{\frac{\delta-x_{0}}{2}+x_{i}} $

(18)

$ \varphi_{r}=\frac{n_{r} x_{0}}{\frac{\delta-x_{0}}{2}-x_{i}} $

(19)

参考式(8)的热阻串联模型，结合式(18)和式(19)，即可得到发热板左右两端复合平板等效导热系数λ_{e, l}、λ_{e, r}。此时可以列出拼接板的温度分布控制方程：

$ \left\{\begin{array}{l} \frac{\mathrm{d}}{\mathrm{d} x}\left(\lambda_{e, l} \frac{\mathrm{d} t}{\mathrm{~d} x}\right)=0,-\frac{\delta}{2} \leqslant x \leqslant x_{i}-\frac{x_{0}}{2} \\ \frac{\mathrm{d}}{\mathrm{d} x}\left(\lambda_{1} \frac{\mathrm{d} t}{\mathrm{~d} x}\right)+\phi=0, x_{i}-\frac{x_{0}}{2} \leqslant x \leqslant x_{i}+\frac{x_{0}}{2} \\ \frac{\mathrm{d}}{\mathrm{d} x}\left(\lambda_{e, r} \frac{\mathrm{d} t}{\mathrm{~d} x}\right)=0, x_{i}+\frac{x_{0}}{2} \leqslant x \leqslant \frac{\delta}{2} \end{array}\right. $

(20)

当边界温度为0，根据3部分界面处的温度和热流连续条件，可以解得拼接板的温度分布：

$ t=\left\{\begin{array}{l} \left(\frac{\phi x_{0}}{\lambda_{e, l}}-\frac{2 c \lambda_{e, r}}{\delta \lambda_{e, l}}\right) x+\frac{\phi x_{0} \delta}{2 \lambda_{e, l}}-c \frac{\lambda_{e, r}}{\lambda_{e, l}},-\frac{\delta}{2} \leqslant x \leqslant x_{i}-\frac{x_{0}}{2} \\ \ \ -\frac{\phi}{2 \lambda_{1}} x^{2}+\left[\frac{\phi\left(x_{i}+\frac{x_{0}}{2}\right)}{\lambda_{1}}-\frac{2 \lambda_{e, r} c}{\lambda_{1} \delta}\right] x+ \\ \ \ \frac{\phi x_{0}\left(x_{i}-\frac{x_{0}}{2}+\frac{\delta}{2}\right)}{\lambda_{e, l}}-\frac{\phi}{2 \lambda_{1}}\left(x_{i}-\frac{x_{0}}{2}\right)\left(x_{i}+\frac{3}{2} x_{0}\right)+ \\ \ \ \frac{c\left(2 x_{i}-x_{0}\right) \lambda_{e, r}}{\delta}\left(\frac{1}{\lambda_{1}}-\frac{1}{\lambda_{e, l}}\right)-c \frac{\lambda_{e, r}}{\lambda_{e, l}},\\ x_{i}-\frac{x_{0}}{2} \leqslant x \leqslant x_{i}+\frac{x_{0}}{2} \\ \ \ -\frac{2 c}{\delta} x+c, x_{i}+\frac{x_{0}}{2} \leqslant x \leqslant \frac{\delta}{2} \end{array}\right. $

(21)

式中：

$ c=\frac{\frac{\phi x_{0}^{2}}{2 \lambda_{1}}+\frac{\phi x_{0}\left(x_{i}-\frac{x_{0}}{2}+\frac{\delta}{2}\right)}{\lambda_{e, l}}}{1+\frac{\lambda_{e, r}}{\lambda_{e, l}}+\frac{x_{0}}{\delta}\left(2 \frac{\lambda_{e, r}}{\lambda_{1}}-\frac{\lambda_{e, r}}{\lambda_{e, l}}-1\right)+\frac{2 x_{i}}{\delta}\left(\frac{\lambda_{e, r}}{\lambda_{e, l}}-1\right)} $

(22)

将温度分布在复合平板的厚度上进行积分并取几何平均，化简即得到一块弥散板发热时的平均温度为：

$ \begin{aligned} t_{\text {avg }, i}=& \frac{\phi x_{0}}{\delta}\left[\frac{3(n-1) \delta x_{0}+(3 n-1) x_{0}^{2}}{24 \lambda_{1}}+\right.\\ & \frac{\delta^{2}-\delta x_{0}(n-1)-n x_{0}^{2}}{8 \lambda_{2}}-\\ & \left.\frac{\delta+\left[(n+1) \frac{\lambda_{2}}{\lambda_{1}}-n\right] x_{0}}{2 \lambda_{2}\left(\delta+x_{0}\right)} x_{i}^{2}\right] \end{aligned} $

(23)

由式(23)可以看出，只有一块弥散板发热时的复合平板平均温度与总弥散板的块数、两相材料的导热系数、弥散板和复合平板的宽度以及热源的位置有关。

根据叠加原理，含分布式内热源复合平板的平均温度等于每块弥散板单独发热时的平均温度以及边界条件t₀的叠加，基于式(23)可以得到完整复合平板的平均温度为：

$ \begin{aligned} t_{\mathrm{avg}} &=\frac{\phi x_{0}}{\delta}\left[n \frac{3(n-1) \delta x_{0}+(3 n-1) x_{0}^{2}}{24 \lambda_{1}}+\right.\\ & \ \ \ \ \ \ \ \ n \frac{\delta^{2}-\delta x_{0}(n-1)-n x_{0}^{2}}{8 \lambda_{2}}-\\ &\left.\frac{\delta+\left[(n+1) \frac{\lambda_{2}}{\lambda_{1}}-n\right] x_{0}}{2 \lambda_{2}\left(\delta+x_{0}\right)} \sum\limits_{i=1}^{n} x_{i}^{2}\right]+t_{0} \end{aligned} $

(24)

由于发热弥散板为均匀分布，结合式(1)和式(2)，则所有发热板中心位置坐标平方的求和为：

$ \begin{gathered} \sum\limits_{i=1}^{n} x_{i}^{2}=2\left(\frac{\delta+x_{0}}{2(n+1)}\right)^{2}\left(1^{2}+3^{2}+\cdots+\left(2 \frac{n}{2}-1\right)^{2}\right)= \\ 2\left(\frac{\delta+x_{0}}{2(n+1)}\right)^{2} \frac{n(n+1)(n-1)}{6}= \\ \left(\delta+x_{0}\right)^{2} \frac{n(n-1)}{12(n+1)} \end{gathered} $

(25)

将式(25)代入式(24)化简后得到：

$ \begin{gathered} t_{\text {avg }}=\frac{\phi x_{0}}{\delta}\left[\frac{n x_{0}^{2}}{12 \lambda_{1}}+\frac{n\left(\delta^{2}-x_{0}^{2}\right)}{8 \lambda_{2}}+\right. \\ \left(\frac{1}{\lambda_{1}}-\frac{1}{\lambda_{2}}\right) \frac{x_{0}}{2}\left(\delta+x_{0}\right) \frac{n(n-1)}{6}- \\ \left.\frac{1}{2 \lambda_{2}}\left(\delta+x_{0}\right)^{2} \frac{n(n-1)}{12(n+1)}\right]+t_{0} \end{gathered} $

(26)

进一步将式(26)代入式(12)，化简后可得含内热源复合平板等效导热系数的无量纲形式：

$ \begin{gathered} \frac{\lambda_{e, s}}{\lambda_{2}}=1 /\left\{\kappa \mu^{2}+\frac{3}{2}\left(1-\mu^{2}\right)-\right. \\ \left.(\mu+1)(n-1)\left[\frac{\mu+1}{2(n+1)}-\mu(\kappa-1)\right]\right\} \end{gathered} $

(27)

由式(27)可知，预测含内热源复合平板平均温度的等效导热系数模型，与不含内热源时复合平板的等效导热系数模型式(8)，有较大差异。含内热源时，等效导热系数的表达式更加复杂，不仅与材料导热系数相关，还与热源的几何尺寸μ和数量n分别相关，体现一定的分布效应，主要来自分布式热源的影响。而在式(8)无内热源的等效导热系数模型中，弥散材料的影响则体现出整体性，与弥散材料的填充率(即数量与尺寸的乘积)相关，主要由弥散材料的几何比例决定。

同时，从式(27)可以看出，即使内热源的存在显著改变了平均温度及其等效导热系数规律，但内热源的大小没有直接影响等效导热系数模型中。

3 结论

1) 内热源的存在改变了复合材料中的热流分布，其等效导热系数的研究不应再选择热流作为守恒量，可选择平均温度作为守恒量。

2) 预测含内热源复合平板平均温度的等效导热系数模型与两相材料导热系数大小、发热材料数量、发热材料相对尺寸，即热源的分布情况相关，与热源的绝对大小无关。

本文提出了研究含内热源复合材料等效导热系数的初步思路和方法，为今后研究内热源分布更为复杂的传热问题(如内热源分布非均匀、颗粒弥散型燃料球)奠定了基础。