多孔介质是一种有效的传热强化介质^［1-3］，它被广泛应用于建筑、食品、化学品、煤炭、石油开采油气贮藏、热管虹吸材料、化工固定床、航空航天复合材料、人体和动物的组织器官、植物体内的水分和营养输送系统等诸多领域，因此研究多孔介质的热传递在许多工程和科学领域都具有重要意义^［4-7］。

王志国等^［8］基于表征体元建立了“三箱”模型，考虑含湿多孔介质的孔隙率、含湿率、孔道通道分布和迂曲度与导热系数的关系，而实际上，选取局部容积平均范围即表征体元(REV)是很困难的，粗宏观水平的选取只是将多孔介质的特征长度限定在特定范围内，不够全面且不具有代表性。马永亭^［9］采用欧姆定律模型推导多孔介质的等效导热率，且假设多孔介质是由两相（固−液或固−气）组成的带有接触热阻的自相似分布颗粒，并采用谢尔宾斯基地毯模型进行模拟，采用模拟电阻法推导随机多孔介质的等效热导率。该模型只考虑了多孔介质的两相，而未考虑多孔介质加热时发生的相变，所以并不能准确地反映等效导热率，且不能用于三相共存的情形。Datta^［10］通过构造多组分多相建模框架，结合质量守恒方程和能量守恒方程，考虑物质体内质量迁移和物质发生相变时的情况，验证了汉堡肉饼在油炸和接触加热时温度随时间的变化，以及不同位置时的温度分布，并与实验值进行比较。但该框架只是验证了食物内不同位置的温度分布随时间的变化，且采用了许多经验公式，因此不能说明温度变化与孔隙率及温度变化与含水率的关系。胡爱娟等^［11］采用堆积颗粒模型，根据控制方程在考虑相变时得出不同含湿率下加热液体在考虑毛细力与不考虑毛细力下液体分布随时间的变化。但是该模型只是验证了堆积模型在不同含湿率下的受热情况，并不能说明有效导热率的变化。Lu等^［12］通过建立非饱和多孔介质多相流干燥模型，研究了温度、湿度、热空气的速度对干燥过程的影响，在干燥流速较低时利用该模型得到的计算结果与实验结果吻合，但并未说明孔隙率对干燥过程的影响。

本文基于分形理论，考虑多孔介质在含湿情况下，且考虑加热过程中相变的影响，结合加热过程中的质量守恒方程、能量守恒方程、相变公式和傅里叶导热定律，求解有效导热率，且无需求解导热热阻即可直接求出含湿物质的有效导热率。

1 分形多孔介质的结构特点 1.1 孔隙率

由于多孔介质具有自相似的特性，所以能通过分形理论定义多孔介质的孔隙率，定义多孔介质中具有分形结构的团聚体（孔相和固相的细微混合体）所占比例为z，孔相所占比例为x，固相所占比例为y，多孔介质分形模型如图1所示，其中白色表示孔隙，黑色表示固相，深灰色表示团聚体，孔相内含有液相，在图中用浅灰色表示，则有

$ \qquad x + y + {\textit{z}}{\rm { = }}1 $

(1)

孔隙比例在第i步的总和为

$ \qquad \varepsilon {\rm { = }}x+{\textit{z}}x + {{\textit{z}}^2}x + \cdots {\textit { + }}{{\textit{z}}^{i - 1}}x=x\sum\limits_{j = 0}^{i - 1} {{{\textit{z}}^j} = x\dfrac{{{{\textit{z}}^i}- 1}}{{{\textit{z}} - 1}}} $

(2)

式中：ε为^{孔隙率；i为迭代次数；j为最后一次迭代次数。}

根据式(1)得

$\qquad \varepsilon =\dfrac{x}{x + y}\cdot(1-{{\textit{z}}}^{i}) $

(3)

当i取无穷大时，则zⁱ→0，孔相面积比例为

$\qquad {\varepsilon _x} = \dfrac{x}{{x + y}} $

(4)

下文均用ε代替ε_x表示孔隙率。

固相面积比例ε_y为

$ \qquad {\varepsilon _y} = \dfrac{y}{{x + y}} $

(5)

1.2 孔相和固相的分形特性

二维空间的面积分形维数D_f的一般表达式为^［13］

$ \qquad {D_{\rm {f}}} = {d_{\rm {E}}} - \dfrac{{\ln \varepsilon }}{{\ln {\left(\dfrac{{r}_{\mathrm{min}}}{{r}_{\mathrm{max}}}\right)}}} $

(6)

式中：r_max、r_min 分别为最大、最小孔隙直径，m；d_E为欧氏空间维度数，在二维空间内d_E=2；通过相关文献中的数据计算与验证，在二维空间维度内D_f取值在1～2范围内，不包括1和2。

迂曲度分形维数D_t的一般表达式为^［13］

$ \qquad {D_{\rm{t}}}{\rm{ = }}1 + \frac{{\ln \left[ {\dfrac{{{\tau _{{\rm{av}}}}}}{{{D_{\rm{f}}}}}({D_{\rm{f}}} + {D_{\rm{t}}} - 1)} \right]}}{{\ln \left( {\dfrac{{{L_y}}}{{{r_{\min }}}}} \right)}} $

(7)

式中：τ_av为平均迂曲度；经前人验算，D_t的取值在1.08～1.10范围之间；L_y为孔道两端的特征距离，即物料的厚度，m。

由孔相分布函数得^［13］

$\qquad -\text{d}{{N}_{x}}=\frac{x}{z}{{D}_{\text{f}}}{{(\sqrt{n}{{r}_{\max }})}^{{{D}_{\text{f}}}}}{{r}^{-({{D}_{\text{f}}}+1)}}\text{d}r $

(8)

式中：N_x为孔相单元块数目；r为孔道一般直径，m；负号表示孔隙数目与直径成反比；n为单元块填充个数。

固相分布函数为

$\qquad - {\rm {d}}{N_y} = \dfrac{y}{z}{D_{\rm{f}}}{(\sqrt n {r_{\max }})^{{D_{\rm{f}}}}}{r^{ - ({D_{\rm{f}}} + 1)}}{\rm {d}}r $

(9)

式中，N_y为固相单元块数目。

1.3 孔道迂曲特性

由于多孔介质孔道是迂曲分布的，所以气体通过多孔介质的实际距离为

$\qquad {L_x}(r) = L_y^{{D_{\rm {t}}}}{r^{1 - {D_{\rm {t}}}}} $

(10)

式中，L_x(r)为气体在孔道内通过的距离，m。

1.4 孔相横截面积、固相横截面积、液相横截面积和总面积

由微元法得到多孔介质横截面上孔隙总面积A_x为

$ \begin{split} \quad {A_x} = & - \int_{{r_{\min }}}^{{r_{\max }}} {\dfrac{1}{4}} {\text{π}}{r^2}{\rm {d}}N = \int_{{r_{\min }}}^{{r_{\max }}} {\dfrac{1}{4}} {{\text{{π}}}} {r^2}{D_{\rm{f}}}r_{\max }^{{D_{\rm{f}}}}{r^{ - ({D_{\rm{f}}} + 1)}}{\rm {d}}r =\\ & \dfrac{{{\text{π}} {D_{\rm{f}}}r_{\max }^2}}{{4(2 - {D_{\rm{f}}})}}\left[1 - {\left(\dfrac{{{r_{\min }}}}{{{r_{\max }}}}\right)^{2 - {D_{\rm{f}}}}}\right] \end{split} $

(11)

根据式(5) 可得固相总面积A_y为

$ \qquad {A}_{y}=\dfrac{1-{\varepsilon }_{x}}{{\varepsilon }_{x}} \cdot \dfrac{{\text{π}}{D}_{\rm {f}}{\textit{r}}_{\mathrm{max}}^{2}}{4\left(2-{D}_{\rm {f}}\right)}\left[1-{\left(\dfrac{{r}_{\mathrm{min}}}{{r}_{\mathrm{max}}}\right)}^{2-{D}_{\rm {f}}}\right]$

(12)

图2为含湿多孔介质孔道横截面示意图。

对于液相来说，可以先根据多孔介质内的含湿率求出液体的质量。假设液体均匀分布在多孔介质孔道壁面，则可以根据液体密度求出液体的体积，然后再除以孔道端面的特征距离即可得到空隙内液相横截面积A_x1，不过还要考虑蒸发部分减少的面积^［14］A_x，即

$ \qquad {A}_{x1}=\dfrac{\varphi m-\dot{m}}{{\rho }_{\rm {l}}{L}_{y}}=\dfrac{\varphi m-\dfrac{{M}_{\rm {w}}}{RT}\cdot \dfrac{({P}_{\rm {vs}}-{P}_{\rm {vi}})}{t}}{{\rho }_{\rm {l}}{L}_{y}} $

(13)

式中：ṁ为蒸发率，kg·（m³·s）⁻¹；M_w为液相摩尔质量，kg·mol⁻¹；R为通用气体常数，kJ·（mol·K）⁻¹；T为最后加热的温度，K；P_vs为饱和蒸汽分压力，kPa；P_vi为初始蒸汽压力，kPa； $\varphi $ 为含湿率；m为物料总质量，kg；ρ_l为液相密度，kg·m⁻³；t为加热时间。

通过式(9)和(11)可知，在有热风对多孔介质进行加热时，热风通过多孔介质时的横截面积A_x2就相当于孔隙总横截面积A_x减去液相占据的横截面积A_x1，即

$\begin{split} \qquad {A}_{x2}=&{A}_{x}-{A}_{x1}=\dfrac{{\text{π}}\textit{{D}}_{ {{\rm{f}}}}{\textit{r}}_{\mathrm{max}}^{2}}{4\left(2-{D}_{\rm {f}}\right)}\left[1-{\left(\dfrac{{r}_{\mathrm{min}}}{{r}_{\mathrm{max}}}\right)}^{2-{D}_{\rm {f}}}\right]-\\ & \dfrac{\varphi m-\dfrac{{M}_{\rm {w}}}{RT}\cdot\dfrac{\left({P}_{\rm {vs}}-{P}_{\rm {vi}}\right)}{t}}{{\rho }_{\rm {l}}{L}_{y}} \end{split} $

(14)

$\qquad \dot{m}=\dfrac{{M}_{\rm {w}}}{RT}\cdot \dfrac{({P}_{\rm {vs}}-{P}_{\rm {vi}})}{t} $

(15)

特征距离L_y可以认为是多孔介质总横截面积的开平方^［15］，由式(9)、(10)可得多孔介质总截面积A为

$\quad A={A}_{x} + {A}_{y}=\dfrac{1}{{\varepsilon }_{x}}\cdot\dfrac{{{D}}_{\rm {f}}{{r}}_{\mathrm{max}}^{2}}{4\left(2-{D}_{\rm {f}}\right)}\left[1-{\left(\dfrac{{r}_{\mathrm{min}}}{{r}_{\mathrm{max}}}\right)}^{2-{D}_{\rm {f}}}\right] $

(16)

所以特征距离L_y可表示为

$ \quad {L}_{y}=\sqrt{A}={\left \{\dfrac{1}{{\varepsilon }_{x}}\cdot\dfrac{{D}_{\rm {f}}{{r}}_{\mathrm{max}}^{2}}{4\left(2-{D}_{\rm {f}}\right)}\left[1-{\left(\dfrac{{r}_{\mathrm{min}}}{{r}_{\mathrm{max}}}\right)}^{2-{D}_{\rm {f}}}\right]\right\}}^{\frac{1}{2}} $

(17)

2 热量计算 2.1 加热过程中固相吸收的热量

在加热过程中物料温度的变化会影响固相吸收的热量。根据热量计算公式可得

$ \qquad {Q_{\rm {s}}} = \int_{{T_1}}^{{T_2}} {{c_{\rm {s}}}} (1 - \varphi )m{\rm {d}}T = {c_{\rm {s}}}(1 - \varphi )m\Delta T $

(18)

式中：∆T为物料加热前、后的温差，K；Q_s为固相吸收的热量，kJ；c_s为固相的比热容，kJ·（kg·K）⁻¹；T₁、T₂分别为物料加热前、后的温度，℃。

2.2 加热热风具有的热量

在考虑孔道内吸收的热量时，需要考虑液相部分蒸发潜热。孔道内吸收的热量Q_h包含热空气吸收的热量Q_a、发生相变部分的液体吸收的热量Q_q和液相部分吸收的热量Q₁^［16］，则孔道吸收的总热量Q_h表示为

$\qquad {Q_{\rm {h}}} = {Q_{\rm {a}}} + {Q_{\rm {l}}} + {Q_{\rm {q}}} $

(19)

$ \begin{split} \qquad {Q}_{\rm {a}}=&{\displaystyle {\int }_{0}^{t}{c}_{\rm{a}}}{\rho }_{\rm{a}}{\textit{v}}{{\textit{A}}}_{{\textit{x}}2}\Delta {\textit{T}}{\rm {d}}t={c}_{\rm{a}}{\rho }_{\rm {a}}{\textit{vt}}{{\textit{A}}}_{{\textit{x}}2}\Delta {\textit{T}}=\\ &{c}_{\rm {a}}{\rho }_{\rm {a}}vt \left\{\dfrac{{\text{π}} {D}_{\rm {f}}{r}_{\mathrm{max}}^{2}}{4\left(2-{{\textit{D}}}_{\rm {f}}\right)}\left[1-\left({\dfrac{{r}_{\mathrm{min}}}{{r}_{\mathrm{max}}}}\right)^{2-D_{\rm{f}}}\right]-\right.\\ & \left. \dfrac{\varphi m-\dfrac{{M}_{\rm {w}}}{RT}\cdot\dfrac{\left({P}_{\rm {vs}}-{P}_{\rm {vi}}\right)}{t}}{{\rho }_{\rm {l}}{L}_{y}}\right\}\Delta T \\[-30pt] \end{split} $

(20)

$\qquad {Q_{\rm {l}}} = \int_{{T_1}}^{{T_2}} {{c_{\rm {l}}}} (\varphi m - \dot m){\rm {d}}T = {c_{\rm {l}}}(\varphi m - \dot m)\Delta T $

(21)

式中：c_a为空气的比热容，kJ·（kg·K）⁻¹；c_l为液体的比热容，kJ·（kg·K）⁻¹；v为热风速率，m·s⁻¹；ρ_a为空气的密度，kg·m⁻³ 。

由蒸发率及汽化潜热所需热量的计算公式可得

$\qquad {Q}_{\rm {q}}=\dot{m}{L}_{\rm {t}}=\dfrac{{M}_{\rm {w}}}{RT}\cdot\dfrac{({P}_{\rm {vs}}-{P}_{\rm {vi}}){L}_{\rm {t}}}{t} $

(22)

式中，L_t为液相蒸发潜热，kJ·kg⁻¹。

2.3 加热物质吸收的总热量

多孔介质吸收的总热量Q为不同相吸收热量的总和，即

$\begin{split} \qquad Q=& {Q}_{\rm {s}} + {Q}_{\rm {h}}={Q}_{\rm {s}} + {Q}_{\rm {a}} + {Q}_{\rm {l}} + {Q}_{\rm {q}} ={c}_{\rm {s}}(1-\varphi )m\Delta T +\\ & {c}_{\rm {a}}{\rho }_{\rm {a}}vt\Delta T\left\{\dfrac{{D}_{\rm {f}}{{\textit{r}}}_{\mathrm{max}}^{2}}{4(2-{D}_{\rm {f}})}\left[1-\left({\dfrac{{r}_{\mathrm{min}}}{{r}_{\mathrm{max}}}}\right)^{2-D_{\rm{f}}}\right] -\right.\\ &\left. \dfrac{\varphi m-\dfrac{{M}_{\rm {w}}}{RT}\cdot\dfrac{({P}_{\rm {vs}}-{P}_{\rm {vi}})}{t}}{{\rho }_{\rm {l}}{L}_{y}}\right\} + {c}_{\rm {l}}(\varphi m-\dot{m})\Delta T + \\ &\dfrac{{M}_{\rm {w}}}{RT}\cdot\dfrac{({P}_{\rm {vs}}-{P}_{\rm {vi}}){L}_{\rm {t}}}{t} \\[-20pt] \end{split}$

(23)

3 有效导热率的相关计算

傅里叶导热定律可用于多孔介质的导热，即

$\qquad {\lambda _{{\rm {eff}}}} = \dfrac{{Q{L_y}}}{{\Delta TA}} $

(24)

式中，λ_eff为多孔介质的有效导热率，W·(m·K)⁻¹。

因为多孔介质结构复杂，所以热量通过多孔介质的不同相时，导热路径不同，由此得出改进后的傅里叶导热定律^［17］为

$ \qquad {\lambda _{{\rm {eff}}}} = \dfrac{{Q{L_y}}}{{\Delta T{A_{x1}}}} + \dfrac{{Q{L_y}}}{{\Delta T{A_y}}} + \dfrac{{Q{L_x}}}{{\Delta T{A_{x2}}}} $

(25)

联立式(8)、(10)、(11)、(12)、(14)和(21)，得到含湿多孔介质的有效热导率为

$ \begin{split} &{{\lambda }_{\rm {eff}}}=\dfrac{Q}{\Delta T}\centerdot \left(\dfrac{{{L}_{y}}}{{{A}_{x1}}} + \dfrac{{{L}_{y}}}{{{A}_{y}}} + \dfrac{{{L}_{x}}}{{{A}_{x2}}}\right) = \Biggr\{{{c}_{\rm {s}}}\left(1 - \varphi \right)m + \Biggr.\\ & {{c}_{\rm {a}}}{{\rho }_{\rm {a}}}vt\left\{\dfrac{{\text{π}}{{D}_{\rm {f}}}r_{\max }^{2}}{4\left(2 - {{D}_{\rm {f}}}\right)}\left[1 - {{\left(\dfrac{{{r}_{\min }}}{{{r}_{\max }}}\right)}^{2 - {{D}_{\rm {f}}}}}\right] - \right. \\ &\left.\dfrac{\varphi m - \dfrac{{{M}_{\rm {w}}}}{RT}\centerdot \dfrac{\left({{P}_{\rm {vs}}} - {{P}_{\rm {vi}}}\right)}{t}}{{{\rho }_{\rm {l}}}{{L}_{y}}}\right\} + {{c}_{\rm {l}}}\varphi m + \dfrac{{{M}_{\rm {w}}}}{RT}\centerdot \\ &\left.\dfrac{\left({{P}_{\rm {vs}}} - {{P}_{\rm {vi}}}\right){{L}_{\rm {t}}}}{t\Delta T}\right\}\centerdot \left\{\dfrac{\dfrac{{{\rho }_{\rm {l}}}}{\varepsilon }\centerdot \dfrac{{\text{π}}{{D}_{\rm {f}}}r_{\max }^{2}}{4\left(2 - {{D}_{\rm {f}}}\right)}\left[1 - {{\left(\dfrac{{{r}_{\min }}}{{{r}_{\max }}}\right)}^{2 - {{D}_{\rm {f}}}}}\right]}{\varphi m - \dfrac{{{M}_{\rm {w}}}}{RT}\centerdot \dfrac{\left({{P}_{\rm {vs}}} - {{P}_{\rm {vi}}}\right)}{t}} + \right. \\ & \dfrac{1}{\left(1 - {{\varepsilon }_{x}}\right) \left\{\dfrac{1}{\varepsilon }\centerdot \dfrac{{\text{π}}{{D}_{\rm {f}}}r_{\max }^{2}}{4\left(2 - {{D}_{\rm {f}}}\right)}\left[1 - {{\left(\dfrac{{{r}_{\min }}}{{{r}_{\max }}}\right)}^{2 - {{D}_{\rm {f}}}}}\right]\right\}^{\frac{1}{2}}} + \\ & \left. \dfrac{{{\left\{\dfrac{1}{\varepsilon }\centerdot \dfrac{{\text{π}}{{D}_{\rm {f}}}r_{\max }^{2}}{4\left(2 - {{D}_{\rm {f}}}\right)}\left[1 - {{\left(\dfrac{{{r}_{\min }}}{{{r}_{\max }}}\right)}^{2 - {{D}_{\rm {f}}}}}\right]\right\}}^{\frac{{{D}_{\rm {t}}}}{2}}}{{r}^{1 - {{D}_{\rm {t}}}}}}{\dfrac{{\text{π}}{{D}_{\rm {f}}}r_{\max }^{2}}{4\left(2 - {{D}_{\rm {f}}}\right)}\left[1 - {{\left(\dfrac{{{r}_{\min }}}{{{r}_{\max }}}\right)}^{2 - {{D}_{\rm {f}}}}}\right] - \dfrac{\varphi m - \dfrac{{{M}_{\rm {w}}}}{RT}\centerdot \dfrac{\left({{P}_{\rm {vs}}} - {{P}_{\rm {vi}}}\right)}{t}}{{{\rho }_{\rm {l}}}{{L}_{y}}}} \right \} \end{split}$

(26)

4 分析与讨论

图3为砂石有效导热率的模拟值与实验值^［18］比较，此时面积分形维数取1.8，迂曲分形维数取1.1，此处选取孔隙率为0.42。从图中可以看出，有效导热率随着含湿率的增加而增大。在孔隙率较小时，有效导热率增加速度稍慢，在含湿率大于0.1时有效导热率随着含湿率的增加基本保持匀速增加。因为随着含湿率的增加，液相占比增加，孔内之前被空气占据的空间变为液体，又因为液体导热率大于气体导热率，所以总体有效导热率会随着含湿率的增加而增大。将模拟值与实验值进行对比发现，当含湿率分别为0.041、0.214、0.525、0.716和0.977时，有效导热率模拟值与实验值的差值分别为0.12、0.36、0.21、0.18和0.03 W·（m·K）^{− 1}。除了在计算和部分数据小数位取舍过程存在计算误差外，实验数据也可能存在测量等其他方面的误差，但可以看出模拟值与实验值偏差很小，所以该模型的准确性得到了很好的验证。

图4为面积分形维数取1.8，迂曲分形维数取1.1，含湿率分别为0.4、0.5和0.6时有效导热系数随孔隙率的变化。从图中可以看出，当含湿率一定时，有效导热率随着孔隙率的增加而减小，且随着孔隙率的增加含湿率对有效导热率的影响逐渐减小。通过比较发现，在含湿率为0.6时有效导热率随孔隙率的变化更为显著；相同孔隙率下，有效导热率随着含湿率的增加而增大；当孔隙率较大时，不同含湿率下的有效导热率更为接近。原因是固相导热率大于液相导热率，而液相导热率又大于气相导热率，所以在孔隙率较大时，气相导热率影响较大，有效导热率下降速度趋缓且趋近于同一稳定值。

图5为含湿率为0.15、迂曲分形维数取1.1时，不同面积分形维数下孔隙率对有效导热率的影响。从图中可以看出，同一孔隙率下，有效导热率随面积分形维数的增加而减小，其原因是面积分形维数越大孔道结构越复杂，热量通过孔道的时间会随之增加，所以会导致总体导热率减小；在面积分形维数一定时，有效导热率随孔隙率的增加而减小。因为气体导热率小于液体和固体导热率，所以孔隙率越大气体占比越大，故导致总体导热率随之减小。

图6为面积分形维数为1.6、含湿率为0.15时，不同迂曲分形维数下孔隙率对有效导热率的影响。从图中可以看出，迂曲分形维数对有效导热率的影响较小，且在孔隙率较大时迂曲分形维数对有效导热率基本无影响，即不同迂曲分形维数下的有效导热率曲线接近重合，在孔隙率较小(如ε=0.1)时，能够清楚地观察到有效导热率随迂曲分形维数的增大而减小。因为在含有液相时，迂曲面将被液体覆盖，且孔道迂曲度的增加使传热距离增加，所以会导致有效导热率下降。而且在孔隙率较大时，多孔介质内气体导热占主导，所以迂曲分形维数在此时对有效导热率的影响更小。

5 结 论

（1）根据完全分形理论，建立了在含湿多孔介质的传热模型。该模型与传统的假设热阻和相关经验公式法相比，能通过具体分形维数和实际测算更好地反映含湿多孔介质的传热情况。

（2）在含湿率较高时，有效导热率较大。因为液相导热率大于气相导热率，所以随着含湿率的增加气相占比减少，多孔介质整体有效导热率增加，同样，在含湿率一定的情况下，孔隙率越大，有效导热率越低。

（3）面积分形维数越大，有效导热率越低。因为面积分形维数越大，孔道表面越粗糙，将会增加传热的路径，导致有效导热率减小。迂曲分形维数对有效导热率的影响较小，并且在孔隙率较大时，分形维数对有效导热率的影响很小。所以孔隙率越大分形维数对有效导热率的影响越小。

（4）通过比较模拟值与实验值发现，当孔隙率一定时，有效导热率随着含湿率的增加呈线性增加，只是在含湿率较低时有效导热率增加得较为缓慢。因为在含湿率较小时，液体会附着在孔道壁上，所以有效导热率增加得缓慢。