基于扩散界面法的热管毛细芯气液界面蒸发过程模拟研究

常 建，蔡杰进，谭 冰

(华南理工大学 电力学院，广东 广州 510640)

空间热管堆采用热管作为堆芯向外界热传导的媒介，由于缺少重力的作用，毛细芯对冷却剂的泵送就更重要[1-2]。为更有效地进行传热，对毛细芯的相变研究越来越有必要。由于毛细压力的限制，在毛细芯中所能提供的最大供液流量限制了热管所能承受的最大热输入，过高则会导致毛细芯的烧干。为改善热管性能，毛细结构需优化其泵送能力和渗透率[3-4]。在大多数情况下，毛细芯结构被认为是一个连续的多孔介质，但越来越多的研究发现薄膜蒸发和局部孔隙的马拉格尼对流等微观尺度现象的影响在相变过程中起重要作用。一般来说，薄膜蒸发发生在几个微米长的区域附近的固-液-汽交界处，长期以来一直被认为是在这样的系统中传热的主要模式。薄膜传热的高效率是由大的分离压力梯度以及极低的热阻共同影响的[5-6]。三相接触点附近的强烈蒸发导致了沿半月板界面方向的温度梯度，进而导致表面张力梯度，从而产生热毛细对流[7-8]。

本文通过对单个毛细芯孔隙内的界面蒸发过程进行模拟，以对整个毛细芯的相变提供更准确的模型。首先在相场法的基础上，建立一个从毛细弯液面蒸发到空气中的三维模型，考虑界面处的马拉格尼效应、热浮力效应、界面处的蒸发相变过程以及蒸汽扩散产生的反作用力。然后将模拟结果与实验进行对比验证。最后对界面处的蒸发进行敏感性分析，讨论壁面过热度、接触角和管径对相变过程的影响。

1 薄膜蒸发传热传质模型

扩散界面法不仅是一个简单的数学计算，而且包含着一定的物理意义。它是通过相场变量φ来得到界面层的信息，而不是去直接追踪两种流体界面的变化。它把表面张力等效为场变量的梯度与化学势的乘积，并将它作为一个体积力加入到Navier-Stokes方程中。相场法不但可计算流体界面的对流，而且还保证系统总能量合理减少[9]。

扩散界面法认为界面是一个薄的、物理性质连续变化的区域，且两相的物性参数也在界面连续变化。本文用相场变量φ区分两相，类似于水平集方法。计算时，第1相的值取-1，第2相的值取1，界面处的值从-1到1变化[10]。计算区域内每相所占的份额V可写成相场变量的函数：

(1)

式中，下标v、l分别表示气体和液体。

在模型计算中涉及的物性参数，包括导热系数、黏度、密度以及定压热容等，可定义为：

B=BlVf,l+BvVf,v

(2)

式中，B为相界面对应的物性参数。

λ为混合能量密度，其与表面张力和扩散界面厚度关系如下所示：

(3)

式中，εpf为界面厚度参数，取决于界面处网格的大小。为提高界面变化的准确度，εpf应尽可能小，以保证相场变量φ及两相物理性质在相界面的均匀变化，但同时为保证数值计算的收敛性，相界面的厚度也要与在相界面处的网格密度相匹配。通常取界面处最大网格的1/2作为εpf的值，具有较好的模拟结果。

(4)

式中：ρ为密度；n为界面法向量。对n向量求偏导，并整理可得：

(5)

(6)

式中：u为速度场；ξ为相界面平滑函数：

(7)

在核态沸腾的过程中，通过相界面所发生的传热传质过程是最为重要和特殊的过程，在相变过程中，热流密度和汽化潜热对相变过程的影响十分重要。相变过程的质量转换率可写为：

(8)

式中：q为热流密度；ΔHvl为汽化潜热；Ml为液体摩尔质量。在相界面处，考虑热流仅通过导热方式传导。将热流密度通过傅里叶导热定律展开，式(8)变为：

(9)

式中：T为温度；k为导热系数。上述推导中，主要考虑了流体蒸发造成的界面移动，由于有出口存在，蒸汽的产生未对界面的移动产生直接影响，但蒸发时相变带来的密度变化大，在微小的空间内，汽体从出口扩散时对界面的反作用力不容忽视。本文考虑由于相变导致的气液移动速度，汽体对界面的反作用力fvl可表示为：

(10)

在确定了相变质量通量后，可对其他控制方程同样通过添加源项和界面平滑参数的方式，修正在相界面处的表达式。

相场控制方程(C-H方程)：

(11)

(12)

式中：t为时间；γ为迁移率。由于各方程所引入的参数较多，有必要对各参数进行更详细的解释，以帮助理解方程意义。在对气泡运动的研究中，取常温下水的表面张力σ为0.076 N/m。

(13)

在微尺度传热下，这些温差变化会大到足以对流场造成实质性的影响。此外，由于流体会传递热量，因此温度场也会受到流场变化的影响。

动量方程(N-S方程)：

(14)

其中：

(15)

式中：Fg为重力项；I为应力张量。

能量方程：

(16)

式中，cp为比定压热容。

2 计算模型与网格验证

通道内最小的网格单元边长为10 μm，其长度远大于滑移长度，计算中视为不可滑移边界条件。首先，壁面处没有热通量，流体在毛细力的作用下上升，弯液面在重力的作用下重构，并达到稳态，然后添加壁面初始温度，液体在弯液面通过蒸发被转化为蒸汽并通过管中扩散，同时下方的入口处由于毛细作用对流体进行补充，弯液面随着蒸发产生波动。表1列出了计算中所需的物性参数。

表1 数值模拟计算采用物性参数Table 1 Physical property parameter used in numerical simulation calculation

图1为微通道弯液面蒸发计算模型及自适应网格划分。数值模型为采用轴对称的简化的二维通道(图1a)，蓝色部分为液体，灰色部分为气体，初始弯液面构造为圆弧，在相界面初始化过程中，会根据重力和接触角对界面进行重构，入口为压力入口，其值为液面产生的压力，出口压力为0。本研究的重点是弯液面处的蒸发速率，由于弯液面会随着蒸发不断下移，故本文在模拟该模型时，力图在界面梯度较大的区域中创建密集网格。但有时很难预测梯度骤变出现的位置。还有一些瞬态情况，骤变的梯度会移动。这可通过在出现梯度骤变的所有区域创建细化网格来解决，但通常计算成本很高，较好的解决方法是采用自适应网格划分。基础粗化网格用于在一定时间间隔内推进求解，随后使用求解结果基于某个指示函数来细化网格。最后，使用自适应网格再次模拟时间间隔。

本文采用了三角形网格进行区域划分。由于这是一个两相模型，且主要考虑界面处的蒸发，因此自适应基于进行相场变量的梯度。图1d中黑色线为初始界面，红色线为移动后的界面，显示了自适应时间间隔结束时的解。由图1可看出，随着界面的移动，对界面区域，即相场变量的梯度变化大的范围进行了网格细化，达到了所需的计算精度。

本文使用水作为工作液进行了网格无关性研究，如表2所列，生成了4组网格，其对计算结果的影响列于表3。经过对比发现，网格2相比网格1网格数量减少了1/2，计算结果的误差在1%左右，满足网格独立性的验证，对下文的计算中均采用网格3作为网格划分的密度。

a——微通道弯液面蒸发计算模型示意图；b——微通道两相及弯液面处的基础网格； c——根据相界面的移动和梯度变化后经过细化的界面附近网格；d——根据加密后的网格计算得到的移动过后两相界面图1 微通道弯液面蒸发计算模型及自适应网格划分Fig.1 Microchannel meniscus evaporation calculation model and adaptive meshing

表2 网格独立研究中使用的不同网格Table 2 Different grids used in grid independent research

表3 网格独立验证结果Table 3 Grid independent verification result

3 计算结果验证与讨论

Buffone等[11]通过粒子成像系统拍摄不同管径的石英毛细孔隙内的流场和温度场分布，实验工质为乙醇和丙酮，采用ITO镀膜进行电加热，得到了在弯液面附近看到环形涡旋流场。本文采用实验中的相似的物性和条件，模拟了竖直放置的微通道弯液面蒸发，结果如图2所示。模拟所得弯液面附近的流场与Buffone等实验结果进行了对比，验证了数值模型的正确性。在弯液面附近液面中心的流体沿着相界面向三相接触线附近流动，推动壁面附近的流体回流，最终在入口补偿流的作用下，发生偏转，补充位于弯液面中心的流体。这种现象称之为马拉格尼流动。液-汽界面处的马拉格尼对流驱动液体从管中心线的较热区域流向壁的较冷区域，加速了热对流，促进了蒸发。入口处的液体流入速度约为1 mm/s，涡旋附近的流速远大于入口的流速，因此界面附近的流动主要是由温度梯度导致的马拉格尼流引起的。刚流入时，流体流线平行于管轴，直到接近弯液面在涡旋的影响下发生偏转，最后从弯液面的中心流向接触线附近。从Buffone等实验结果也可清晰看到两个涡旋产生在界面附近，同时补给的流体明显发生了偏转，与模拟所得的流场图几乎一致。该结果验证了模型的正确性，表明建立的基于相场法的界面蒸发模型可用于蒸发速率的计算。

图2 竖直放置的微通道中流场图Fig.2 Flow field diagram in vertically placed microchannel

不同壁面过热度下弯液面蒸发时的温度场如图3所示。可看出初始时壁面附近温度分层分布，因为蒸发的作用，界面处于饱和温度，相对于过热的液体，界面温度最低，且从界面向下温度逐渐增高。随着液体在马拉格尼流涡旋作用下，界面处的高温流体被带入管道的中心，最终在弯液面附近产生较大的温度梯度，且越靠近壁面，温度梯度越大，说明壁面附近的蒸发越强烈[12-14]。对比图2中的流场图可看出，温度场的分布和流场的分布具有一致性。不同于毛细孔隙在层流下的壁面分层传热，马拉格尼涡旋的作用下，壁面处的温度可更快地传递到蒸发界面，增强了换热。对于不同的过热度，具有相似的温度场分布，过热度的影响将在接下来的结果讨论中分析。

图3 不同壁面过热度下弯液面蒸发时的温度场Fig.3 Temperature field of meniscus evaporation under different wall superheat degrees

图4为不同尺寸的微通道水平放置时在重力作用下涡旋发生偏移的模拟结果与Hemanth等的实验结果对比。在模拟中，毛细孔隙的轴线是水平的，因此重力效应会影响弯液面附近的流动方式[15]。图中观察到的不对称流动是浮力和马拉格尼效应的综合结果。如图中实验结果所示，图中每个部分均叠加有多个图像，以提高对比度。通过沿垂直中心平面获得的长曝光颗粒条纹图像表征了浮力和热毛细作用对流动的相对影响，涡旋的偏移与模拟结果一致。图中的箭头指向弯液面上停滞点的位置，模拟中该位置处的界面由于补偿流的作用温度相对最低，热毛细作用最小，与两侧所产生的切应力处于平衡状态，当在界面处发生蒸发时，它使弯液面区域液体较离开弯液面的回流液体相对温度低。温度的这种差异同时会产生浮力效应，其强度随管径的增加而增加。对于研究的最小直径为75 μm的毛细孔隙，浮力的影响可忽略不计，并且观察到与竖直的毛细孔隙相同的对称流型(图4a)。对于200 μm的毛细孔隙，涡旋的不对称性开始产生，发现逆时针底部涡旋大于顺时针上部涡旋，这表明存在不对称的环形涡旋，前面提到的停滞点不再位于液面的中心，开始向上偏移。浮力效应随着毛细孔隙直径的增加而增大，如图4c、d所示，逆时针涡旋的大小随着毛细孔隙直径的增加而变大，停滞点沿弯液面向上移动，浮力增加。可看出，该流动被分为两个区域：由表面张力引起运动的上部区域和由密度差驱动的下部区域。同时在模拟中还可看到，直径的增大导致重力在界面的影响也越明显，图4a、b中，上下两侧的三相接触点大致在同一垂直平面上，图4c、d中，则由于重力的作用，下方界面向前伸展，在加热时，这种不平衡增大了下方的蒸发界面，加剧了涡旋的偏移。因此，浮力驱动涡旋的尺寸随毛细直径的增加而增大。在界面处产生的相对温差驱动了浮力引起的流动，其大小取决于蒸发的质量通量。

图5示出了不同直径毛细孔隙弯液面处的蒸发速率。x0为管径，x为界面与对称轴的轴向距离。从图中可看出，对称轴附近的蒸发速率均较小，且在约60%的长度区域内，蒸发速率均变化不大，在靠近壁面的区域内，蒸发速率快速增长，在接近壁面处，不同管径的界面具有相同的蒸发速率。在对称轴附近，直径较小的弯液面蒸发速率较高，其值最多仅有壁面附近蒸发率的1/10，这是由于较小的弯液面从壁面处由涡旋带来的热量更多，其温度梯度更大，因而小管径的弯液面局部蒸发速率更大。从图5还可看出，蒸发主要集中在壁面附近，壁面附近的蒸发速率与管径尺寸无关。

图5 不同直径毛细孔隙界面处的蒸发速率Fig.5 Evaporation rate at interface of different diameter capillarie interstices

图6示出管径为0.4 mm时不同过热度ΔT和接触角下毛细孔隙界面处的蒸发速率。从图6a可看出，蒸发速率是过热度的函数，呈正比关系，过热度的增加并未改变弯液面的相对蒸发速率，即过热度对涡旋的位置未产生影响，壁面处的蒸发速率远大于弯液面中心停滞点。图3中的温度场也同样证明了这一结论，壁面处的界面温度几乎恒定，且温度梯度变化大，蒸发剧烈；中心处由于补偿流的作用，温度梯度小，蒸发慢。

微通道尺寸：a——75 μm；b——200 μm； c——400 μm；d——762 μm图4 微通道水平放置时在重力作用下 涡旋发生偏移的结果Fig.4 Vortex shift results of microchannels placed horizontally under action of gravity

壁面附近的蒸发速度是壁面上液体厚度的函数，而这个厚度是由界面的接触角决定的，为接触角的函数。从图6b可看出壁面处液膜厚度与蒸发速率的关系，其值相对于接触角的增加，变化不大。在接近壁面处，接触角减小导致更好的润湿性，液体更有利于向壁面流动，补偿壁面处的蒸发，且液膜的厚度更小，蒸发相对更剧烈，考虑到轴向液膜的长度，则接触角小的液膜蒸发速率远大于大接触角的情况。本文润湿性在蒸发过程中更主要的是提供更薄的液膜，在厚度一致的情况下，蒸发的速率基本相同，接触角小意味着有更长的薄膜蒸发区域，因此在整个蒸发表面，有更大的蒸发速率。

图6 管径为0.4 mm时不同过热度和接触角下毛细孔隙界面处的蒸发速率Fig.6 Evaporation rates at capillary interstice interface under different superheat degrees and contact angles with pipe diameter of 0.4 mm

4 半月板上蒸汽扩散通量的不均匀性

通过图6可知界面处某点的蒸发速度是液面厚度和温度梯度的函数，对于整个界面的蒸发速率，需对界面作积分以体现在轴向上的液面分布。毛细孔隙界面处蒸发模型示意图如图7所示，d为半径，r为界面上点到对称轴距离，r0为弯液面的等效半径。这里假设重力对液面的影响简化为圆弧，设弯液面最低点的高度为0，则弯液面的某点处的位置表示为：

图7 毛细孔隙界面处蒸发模型示意图Fig.7 Diagram of evaporation model at interface of capillary interstice

r0=d/cosθ

(17)

(18)

实际中，热管毛细芯的流体均是选择润湿性较好的流体，以达到流体的泵送目的，其接触角大多在60°以内，对图6b中60°接触角进行膜厚与蒸发速率的曲线拟合得到：

m″=22.5+(2.1er/0.11+1.3×10-6er/0.023)/cosθ

(19)

式中，m″为蒸发速率，g/(m2·s)。

5 结论

本文通过对单个毛细芯微孔内的毛细蒸发过程进行建模，考虑了马拉格尼效应引起的界面剪切作用和液体区域的热浮力效应，研究了微管内部的液体在弯液面的传热和传质过程，具体结论如下。

1) 在弯液面和管壁的三相接触处较在弯液面中心处的蒸发更强。从弯液面中心到壁的温度梯度导致从管的中心朝向壁的马拉格尼流。

2) 重力的作用导致液体中的涡旋不对称，这种不对称对小管径影响不大，在大于0.4 mm的情况下更显著。

3) 液面中心处的蒸发速率均较小，在靠近壁面的区域内，蒸发速率快速增长，在接近壁面处，不同管径的界面具有相同的蒸发速率。

4) 壁面附近的蒸发速度是壁面上液体厚度和过热度的函数，而该厚度是由界面的接触角决定的，为接触角的函数。

5) 随着孔隙在毛细芯中的液体有更多的比表面积，文中可看出接触线附近蒸发速率最大，因此孔隙减小，蒸发速率增大。但较小的孔隙会导致较差的渗透率，发生烧干，需综合考虑孔隙减小对蒸发的影响。