一种微重力下伴有空穴的高温相变蓄热数值模拟方法

摘要

本发明提供一种微重力条件下伴有空穴的高温相变蓄热数值模拟方法，在同时考虑相变、空穴及辐射的条件下，对复杂的相变传热问题进行数值计算，使得技术人员利用计算机便可获取现场高温蓄热容器内固/液相变材料熔化率的变化及相变时流场、温度场、液相分布，从而为优化蓄热容器设计、抑制空穴、提高蓄热效率及减少“热斑”和“热松脱”现象提供重要的参考依据。本发明方法通过基于焓法、有限控制容积法的计算程序，使得微重力条件下伴有空穴的高温蓄热容器内的熔化/凝固这一复杂的相变问题得到了简便、高效的解决，具有重要的实用价值。

说明书

技术领域

本发明涉及一种计算流体力学领域的方法，具体是一种在微重力 下伴有空穴的高温相变蓄热数值模拟方法，即通过对焓法能量控制方 程组进行计算。

背景技术

长期以来，人们在对相变问题的研究过程中忽略了固液密度变化 在容器内产生空穴所带来的影响。通常情况下，物质在固态和液态时 的密度是不相等的。对于相变材料(PCM，Phase Change Material)容器 内质量一定的PCM来说，相变过程中PCM密度的变化相应地引起 PCM体积变化，导致PCM容器内PCM占据的容积以外的剩余容积 增大或缩小。剩余容积的变化引起热阻的变化，进而影响相变过程的 进行。对于LiF而言，凝固时的体积收缩率约为23％，如此高的体积 收缩率，对相变过程的影响是很大的，不能简单地加以忽略，在相变 过程中必须考虑固液密度的变化。

PCM(如LiF、CaF₂等)凝固时体积收缩的后果是在PCM容器内 形成空穴。空穴的存在增大了热阻，并对相变过程进行的方向产生影 响。空穴的存在还可能使PCM容器产生“热斑”和“热松脱”现象， 这两种现象都会造成PCM容器的破坏。

Kerslake(Kerslake T W，Ibrahim M B.Analysis of thermal energy  storage material with change-of phase volumetric effects.Trans.ASME， Journal of Solar Energy Engineering，1993，115(2)：22～31)的一维分析 结果表明，因PCM凝固收缩产生的空穴使PCM容器壁面峰值温度上 升了200K。但是，PCM容器的两侧壁在传热过程中将发挥重要作用， 大量的热会经由侧壁传到PCM容器内部，从而削弱空穴对传热过程 的影响。为量化空穴影响，Kerslake(Kerslake T W，Ibrahim M  B.Two-dimensional model of a space station freedom thermal energy  storage canister.Trans.ASME，Journal of Solar Energy Engineering，1992， 114(5)：114～121)以SSF吸热器中的PCM容器为例，编写了二维(r， z)计算程序，对微重力状态下周期性相变传热过程进行了数值求解。 空穴模型给定了PCM容器内空穴的形状和位置，空穴内热传递过程 包括蒸汽导热和各表面间的辐射换热。Kerslake假定空穴位于外壁处， 且空穴体积保持恒定不变。但是，实际运行过程中，随着PCM的熔 化和凝固的交替进行，容器内空穴的体积也随之减小或增加，并不保 持为定值。

国内除蒋大鹏(“空间太阳能动力装置吸热-储热器相变材料容腔 的热分析计算”.航天工业总公司第五零一设计部硕士学位论文， 1991.7)曾对含空穴的固液相变蓄热过程进行了二维热分析以外， 未开展更深入的研究，尚无相关的文献可供参考。而蒋大鹏所做的热 分析也只是在计算结果中可以看出考虑了空穴，但并未说明如何计算 和处理相变过程中空穴的体积变化。国内在微重力下伴随有空穴生成 和发展的固液相变过程方面的研究还存在空白。

发明内容

本发明的目的在于针对现有技术的不足，提供一种微重力条件下 伴有空穴的高温相变蓄热数值模拟方法，在同时考虑相变、空穴及辐 射的条件下，对复杂的相变传热问题进行数值计算，使得技术人员利 用计算机便可获取现场高温蓄热容器内固/液相变材料(PCM，Phase  Change Material)熔化率的变化及相变时流场、温度场、液相分布， 从而为优化蓄热容器设计、抑制空穴、提高蓄热效率及减少“热斑” 和“热松脱”现象提供重要的参考依据。

为实现上述目标，本发明采用有限容积法对控制体进行积分，采 用焓法来处理相变过程。

本发明为解决其技术问题所采取的技术方案具体为：

一种微重力条件下伴有空穴的高温蓄热容器内固/液相变数值模 拟方法，其特征在于，所述数值模拟方法包括如下具体步骤：

(1)依据高温蓄热容器的几何模型划分液体及固体计算区域， 使用三角形非结构化网格对上述计算区域进行网格划分；

(2)定义操作条件：R向受重力作用，参考压力为一个大气压； 定义基于BOUSSINESQ假设的自然对流模型；

(3)定义非稳态熔化/凝固模型，定义离散坐标辐射模型；

(4)定义相变材料在气、固、液三相下的材料属性，所述材料 属性具体包括密度、比热容、导热率、粘度、熔化热、液相线及固相 线温度、吸收率和发射率；

(5)建立如下的偏微分能量控制方程组：

$\frac{\partial \left(\mathrm{\rho e}\right)}{\partial t}+\frac{\partial \left(\mathrm{\rho ue}\right)}{\partial z}+\frac{1}{r}\frac{\partial \left(\mathrm{r\rho ve}\right)}{\partial r}=k{▿}^{2}T---\left(1\right)$

式中：e为比焓，J/kg；ρ为密度，kg/m³；t为时间，s； T为温度，K；u，v分别为z方向、r方向的速度分量；式中 以比焓e的形式把相变材料相变的影响考虑进去，因而，该表 达式适用于整个求解区域；比焓e与温度T的关系式可表示成如 下的形式：

$T=\left(\begin{array}{cc}{T}_m+e/c,& e\le 0\\ T_m,& 0<e<\Delta H_m\\ T_m+(e-\Delta H_m)/c,e& \ge \Delta H_m\end{array}\right)---\left(2\right)$

式中：T_m为相变温度，K；ΔH_m为物质单位质量的相变潜热，

J/kg；

(6)定义边界条件：定义高温蓄热容器的内壁为第二类传热边 界即热流边界，外壁和侧壁为绝热边界，其它壁面为流/固耦合边界； 其中，所述热流边界为周期性热流边界，设定日照期和阴影期的时 间，写入周期性热流自定义函数，其中：

a)日照期热流密度q_sun具体的函数关系如下：

q_sun＝15424-24.88(T_wall-1020)，W/m²；

b)阴影期热流密度q_shadow具体的函数关系如下：

q_shadow＝-12314-24.88(T_wall-1020)，W/m²；

其中，T_wall为高温蓄热容器的内壁温度；

(7)定义初始条件：对整个计算区域进行初始化，设定初始温 度及初始速度；

(8)设定监视相变材料温度场及液相分数分布，设定监视相变 材料熔化率变化；

(9)对步骤(5)中的偏微分能量控制方程进行离散化，并利用 步骤(6)、(7)定义的边界条件和初始条件进行封闭和求解；

(10)对整个计算区域初始化，设定时间步长及迭代次数，对计 算区域内的代数方程组反复进行迭代计算，直到满足所设定的迭代精 度为止，完成微重力下伴有空穴的高温蓄热容器内固/液相变数值模 拟；

(11)对计算结果进行后处理，绘制出云图及相关曲线。

进一步地，本发明采用内节点法进行区域离散化。

进一步地，本发明采用控制容积积分法导出离散方程。

进一步地，本发明采用完全显式格式进行离散化求解。

进一步地，步骤(6)中，设定日照期为54min，阴影期为36min。

进一步地，步骤(9)中，任一非壁面的单元节点P(i，j，k) 处，上述偏微分能量控制方程的离散化形式为

$E_P^{t+\mathrm{\Delta t}}-E_P^t=[{\left(\mathrm{kr}\right)}_n(T_N^t-T_P^t)-{\left(\mathrm{kr}\right)}_s(T_P^t-T_S^t)]\mathrm{\Delta \theta \Delta z\Delta t}+$

$(k_t\frac{T_T^t-T_P^t}{r_T\Delta {\theta }_T}-k_b\frac{T_P^t-T_B^t}{r_B\Delta {\theta }_B})\mathrm{\Delta r\Delta z\Delta t}+(k_e\frac{T_E^t-T_P^t}{\Delta z_E}-k_w\frac{T_P^t-T_W^t}{\Delta z_W})r_P\mathrm{\Delta r\Delta \theta \Delta t};$

对于内壁节点，其外部接受太阳辐射热流q(θ，t)，内壁节点离散控 制方程化为

$E_P^{t+\mathrm{\Delta t}}-E_P^t={\left(\mathrm{kr}\right)}_s(T_S^t-T_P^t)\mathrm{\Delta \theta \Delta z\Delta t}+q(\theta ,t)r_i\mathrm{\Delta \theta \Delta z\Delta t}+$

$(k_t\frac{T_T^t-T_P^t}{r_T\Delta {\theta }_T}-k_b\frac{T_P^t-T_B^t}{r_B\Delta {\theta }_B})\mathrm{\Delta r\Delta z\Delta t}+(k_e\frac{T_E^t-T_P^t}{\Delta z_E}-k_w\frac{T_P^t-T_W^t}{\Delta z_W})r_P\mathrm{\Delta r\Delta \theta \Delta t},$

其中：

r_i代表PCM容器内壁的外径；

E、W、N、S、T、B分别代表节点P的六个相邻节点：(i+1， j，k)、(i-1，j，k)、(i，j+1，k)、(i，j-1，k)、(i，j，k+1) 和(i，j，k-1)；

(kr)_n、(kr)_s、k_e、k_w、k_t、k_b为：

${\left(\mathrm{kr}\right)}_n=\frac{1}{\frac{\Delta r_N^{-}}{k_N(r_N-\frac{\Delta r_N^{-}}{2})}+\frac{\Delta r_P^{+}}{k_P(r_P+\frac{\Delta r_P^{+}}{2})}},$${\left(\mathrm{kr}\right)}_s=\frac{1}{\frac{\Delta r_S^{+}}{k_S(r_S+\frac{\Delta r_S^{+}}{2})}+\frac{\Delta r_P^{-}}{k_P(r_P-\frac{\Delta r_P^{-}}{2})}},$

$k_e=\frac{\Delta z_E}{\frac{{\left(\Delta z_E\right)}^{+}}{k_E}+\frac{{\left(\Delta z_E\right)}^{-}}{k_P}},$$k_w=\frac{\Delta z_W}{\frac{{\left(\Delta z_W\right)}^{+}}{k_P}+\frac{{\left(\Delta z_W\right)}^{-}}{k_W}}$

$k_t=\frac{r_T\Delta {\theta }_T}{\frac{{\left(r_T\Delta {\theta }_T\right)}^{+}}{k_T}+\frac{{\left(r_T\Delta {\theta }_T\right)}^{-}}{k_P}},$$k_b=\frac{r_B\Delta {\theta }_B}{\frac{{\left(r_B\Delta {\theta }_B\right)}^{+}}{k_P}+\frac{{\left(r_B\Delta {\theta }_B\right)}^{-}}{k_B}}.$

进一步地，步骤(10)中，时间步长的稳定性限制条件为：

$\mathrm{\Delta t}\le \frac{(\rho c_P/k)r_P\mathrm{\Delta r\Delta \theta \Delta z}}{(\frac{1}{\frac{\mathrm{\Delta r}}{r_N-\frac{\mathrm{\Delta r}}{4}}+\frac{\mathrm{\Delta r}}{r_P+\frac{\mathrm{\Delta r}}{4}}}+\frac{1}{\frac{\mathrm{\Delta r}}{r_P-\frac{\mathrm{\Delta r}}{4}}+\frac{\mathrm{\Delta r}}{r_S+\frac{\mathrm{\Delta r}}{4}}})2\mathrm{\Delta \theta \Delta z}+\frac{2r_P\mathrm{\Delta r\Delta \theta }}{\mathrm{\Delta z}}+\frac{2\mathrm{\Delta r\Delta z}}{r_P\mathrm{\Delta \theta }}}.$

进一步地，本发明采用计算PCM容器内每一个时间步长的前后 两个时刻t和t+Δt时刻的固态PCM的质量之差Δm作为发生相变的 PCM的质量。

进一步地，每一个时间步长内都要检查当前时刻容器内PCM的 总质量，若存在偏差，则对Δm进行偏差修正，从而保证质量守恒； 如果Δm不为0，说明PCM发生了相变。

进一步地，如果修正后的Δm＞0，说明固态PCM的质量增加， PCM发生了凝固；如果修正后的Δm＜0，说明固态PCM的质量减少， PCM发生了熔化；如果修正后的Δm＝0，说明固态PCM的质量没有 发生变化，PCM正处于显式吸热或放热状态。

进一步地，如果Δm不为0，说明PCM发生了相变，由Δm求出 每一个时间步长内PCM的体积变化量ΔV，根据计算得到的PCM体 积变化量ΔV，对PCM容器的空穴体积进行调整：如果ΔV＞0，说明 PCM的体积增加，就减小空穴所占的体积；如果ΔV＜0，说明PCM 的体积缩小，相应地增大空穴体积。

进一步地，PCM空穴内的辐射换热计算基于以下假设：(1)所 有空穴表面均为漫反射灰体表面；(2)PCM表面吸收所有波长的辐 射；(3)空穴内蒸气不参与辐射换热；空穴界面间的辐射换热主要集 中在径向，轴向和周向的辐射换热相比弱得多，为简化计算，忽略轴 向和周向辐射换热。

进一步地，空穴有效传热系数k_veff按下式确定

$k_{\mathrm{veff}}=k_{\mathrm{PCM}}+k_{\mathrm{reff}}=k_{\mathrm{PCM}}+\ln \left(\frac{r_w}{r_v}\right)·\frac{r_v\sigma (T_w\left(t\right)+T_v\left(t\right))(T_w^2\left(t\right)+T_v^2\left(t\right))}{\frac{1}{{ϵ}_{\mathrm{PCM}}}+\frac{r_v}{r_w}(\frac{1}{{ϵ}_w}-1)},$

式中k_reff为空穴界面间辐射换热的折合当量导热系数，σ为斯蒂 芬玻尔兹曼常数(σ＝5.67×10^-8W·m^-2·K^-4)，T_W为外壁温度，T_V为PCM在空穴界面处的温度，ε_PCM为PCM发射率，ε_W为容器壁发射 率，v代表空穴界面。

本发明方法通过基于焓法和有限控制体积法的计算程序，使得微 重力下伴有空穴的高温蓄热容器内的熔化/凝固这一复杂的相变问题 得到了简便、高效的解决，具有重要的实用价值。

下面对本发明中偏微分控制方程的离散化方法以及代数方程的 求解方法进行介绍。

为进行数值计算，首先要进行区域离散化，即把求解区域划分成 许多个互不重叠的子区域，由子区域内有代表性的点即节点上的待求 变量值来表示区域内连续变化的待求变量场。根据节点在子区域内位 置的不同，区域离散化方法可分为两类：外节点法和内节点法。由于 内节点法处理物性变化的情况比较方便，故本发明优选采用内节点 法。

把求解区域离散化后，就可以将控制微分方程在各个节点上离散， 针对每个节点建立相应的离散方程。控制微分方程的离散方法有泰勒 级数展开法、多项式拟合法、控制容积积分法、平衡法以及变分法等。 控制容积积分法又称有限容积法，是传热数值计算中广泛采用的方 法。该方法推导过程的物理概念清晰，推导结果具有明确的物理意义， 离散方程的守恒特性可以得到保证。优选地，本发明采用控制容积积 分法导出离散方程。

对于差分格式的选取，Thibault曾经做过详细的研究。通过求解 三维导热问题，比较了三种显式格式、四种ADI格式和两种隐式格 式的计算结果与精确解的误差、易编程性、所需计算时间以及对计算 存储量的要求，指出最好的是两种ADI格式，完全显式格式排在第 三位。完全显式格式的主要缺点是受稳定性的限制，时间步长不能取 得过大。对本发明的研究来说，由于相变过程中PCM的密度发生变 化，容器内存在空穴，随着相变过程的进行，空穴体积增加或缩小， 这就决定了时间步长不能取得太大，以避免计算结果失真。当时间步 长较小时，隐式格式的优势也就不存在了，而且由于需要迭代，计算 时间反而比显示格式要长。因此，本发明优选采用完全显式格式进行 离散化求解。

将PCM容器按r×θ×z(即柱坐标方式，r代表径向、θ代表周向、 z代表轴向)划分网格单元。任一单元节点P(i，j，k)的控制容积 如图1所示。节点P所在的轴向截面如图2所示。

将能量控制方程在节点P(i，j，k)的控制容积(参见图1)上作 积分：

$\int \int \int \int \frac{\partial \left(\mathrm{\rho e}\right)}{\partial t}\mathrm{dVdt}=\int \int \int \int \frac{1}{r}\frac{\partial }{\partial r}\left(\mathrm{kr}\frac{\partial T}{\partial r}\right)\mathrm{dVdt}\int \int \int \int +\frac{1}{r}\frac{\partial }{\partial \theta }\left(\frac{k}{r}\frac{\partial T}{\partial \theta }\right)\mathrm{dVdt}+$

$\int \int \int \int \frac{\partial }{\partial z}\left(k\frac{\partial T}{\partial z}\right)\mathrm{dVdt}---\left(3\right)$

表达式(3)等号左边非稳态项在时间和空间上均采用阶梯形分布， 即

${\int }_r^{r+\mathrm{\Delta r}}{\int }_{\theta }^{\theta +\mathrm{\Delta \theta }}{\int }_z^{z+\mathrm{\Delta z}}{\int }_t^{t+\mathrm{\Delta t}}\frac{\partial \left(\mathrm{\rho e}\right)}{\partial t}\mathrm{rdrd\theta dzdt}=\rho (e^{t+\mathrm{\Delta t}}-e^t)r_P\mathrm{\Delta r\Delta \theta \Delta z}=E_P^{t+\mathrm{\Delta t}}-E_P^t---\left(4\right)$

式中E_P代表节点P(i，j，k)所在控制容积内PCM的焓值，m_P为节点P(i，j，k)所在控制容积内PCM的质量；

表达式(3)等号右边扩散项在时间上采用阶梯形分布，在空间上 采用分段线性分布，则有

${\int }_r^{r+\mathrm{\Delta r}}{\int }_{\theta }^{\theta +\mathrm{\Delta \theta }}{\int }_z^{z+\mathrm{\Delta z}}{\int }_t^{t+\mathrm{\Delta t}}\frac{1}{r}\frac{\partial }{\partial r}\left(\mathrm{kr}\frac{\partial T}{\partial r}\right)\mathrm{rdrd\theta dzdt}=[{\left(\mathrm{kr}\right)}_n(T_N^t-T_P^t)-{\left(\mathrm{kr}\right)}_s(T_P^t-T_S^t)]\mathrm{\Delta \theta \Delta z\Delta t}---\left(5\right)$

${\int }_r^{r+\mathrm{\Delta r}}{\int }_{\theta }^{\theta +\mathrm{\Delta \theta }}{\int }_z^{z+\mathrm{\Delta z}}{\int }_t^{t+\mathrm{\Delta t}}\frac{1}{r}\frac{\partial }{\partial \theta }\left(\frac{k}{r}\frac{\partial T}{\partial \theta }\right)\mathrm{rdrd\theta dzdt}=(k_t\frac{T_T^t-T_P^t}{r_T\Delta {\theta }_T}-k_b\frac{T_P^t-T_B^t}{r_B\Delta {\theta }_B})\mathrm{\Delta r\Delta z\Delta t}---\left(6\right)$

${\int }_r^{r+\mathrm{\Delta r}}{\int }_{\theta }^{\theta +\mathrm{\Delta \theta }}{\int }_z^{z+\mathrm{\Delta z}}{\int }_t^{t+\mathrm{\Delta t}}\frac{\partial }{\partial z}\left(k\frac{\partial T}{\partial z}\right)\mathrm{rdrd\theta dzdt}=(k_e\frac{T_E^t-T_P^t}{\Delta z_E}-k_w\frac{T_P^t-T_W^t}{\Delta z_W})r_P\mathrm{\Delta r\Delta \theta \Delta t}---\left(7\right)$

式中，E、W、N、S、T、B分别代表节点P的六个相邻节点：(i+1， j，k)、(i-1，j，k)、(i，j+1，k)、(i，j-1，k)、(i，j，k+1)和(i，j， k-1)。(kr)_n、(kr)_s、k_e、k_w、k_t、k_b意义如下。

${\left(\mathrm{kr}\right)}_n=\frac{1}{\frac{\Delta r_N^{-}}{k_N(r_N-\frac{\Delta r_N^{-}}{2})}+\frac{\Delta r_P^{+}}{k_P(r_P+\frac{\Delta r_P^{+}}{2})}},$${\left(\mathrm{kr}\right)}_s=\frac{1}{\frac{\Delta r_S^{+}}{k_S(r_S+\frac{\Delta r_S^{+}}{2})}+\frac{\Delta r_P^{-}}{k_P(r_P-\frac{\Delta r_P^{-}}{2})}},$

$k_e=\frac{\Delta z_E}{\frac{{\left(\Delta z_E\right)}^{+}}{k_E}+\frac{{\left(\Delta z_E\right)}^{-}}{k_P}},$$k_w=\frac{\Delta z_W}{\frac{{\left(\Delta z_W\right)}^{+}}{k_P}+\frac{{\left(\Delta z_W\right)}^{-}}{k_W}}$

$k_t=\frac{r_T\Delta {\theta }_T}{\frac{{\left(r_T\Delta {\theta }_T\right)}^{+}}{k_T}+\frac{{\left(r_T\Delta {\theta }_T\right)}^{-}}{k_P}},$$k_b=\frac{r_B\Delta {\theta }_B}{\frac{{\left(r_B\Delta {\theta }_B\right)}^{+}}{k_P}+\frac{{\left(r_B\Delta {\theta }_B\right)}^{-}}{k_B}}$

由式(4)、(5)、(6)、(7)得偏微分能量方程(1)的离散化形式 如下

$E_P^{t+\mathrm{\Delta t}}-E_P^t=[{\left(\mathrm{kr}\right)}_n(T_N^t-T_P^t)-{\left(\mathrm{kr}\right)}_s(T_P^t-T_S^t)]\mathrm{\Delta \theta \Delta z\Delta t}+$

$(k_t\frac{T_T^t-T_P^t}{r_T\Delta {\theta }_T}-k_b\frac{T_P^t-T_B^t}{r_B\Delta {\theta }_B})\mathrm{\Delta r\Delta z\Delta t}+(k_e\frac{T_E^t-T_P^t}{\Delta z_E}-k_w\frac{T_P^t-T_W^t}{\Delta z_W})r_P\mathrm{\Delta r\Delta \theta \Delta t}---\left(8\right)$

对于内壁节点，其外部接受太阳辐射热流q(θ，t)，节点离散控制方 程化为

$E_P^{t+\mathrm{\Delta t}}-E_P^t={\left(\mathrm{kr}\right)}_s(T_S^t-T_P^t)\mathrm{\Delta \theta \Delta z\Delta t}+q(\theta ,t)r_i\mathrm{\Delta \theta \Delta z\Delta t}+$

$(k_t\frac{T_T^t-T_P^t}{r_T\Delta {\theta }_T}-k_b\frac{T_P^t-T_B^t}{r_B\Delta {\theta }_B})\mathrm{\Delta r\Delta z\Delta t}+(k_e\frac{T_E^t-T_P^t}{\Delta z_E}-k_w\frac{T_P^t-T_W^t}{\Delta z_W})r_P\mathrm{\Delta r\Delta \theta \Delta t}---\left(9\right)$

式中r_i代表PCM容器内壁的外径。

如果某单元的相邻单元为空穴单元，则式(8)、(9)中计算与相 邻单元的换热量时，应采用空穴复合传热系数按与空穴单元相邻的非 空穴单元的换热进行计算。

${\left(\mathrm{kr}\right)}_n=\frac{2k}{\frac{\mathrm{\Delta r}}{r_N-\frac{\mathrm{\Delta r}}{4}}+\frac{\mathrm{\Delta r}}{r_P+\frac{\mathrm{\Delta r}}{4}}},$${\left(\mathrm{kr}\right)}_s=\frac{2k}{\frac{\mathrm{\Delta r}}{r_P-\frac{\mathrm{\Delta r}}{4}}+\frac{\mathrm{\Delta r}}{r_S+\frac{\mathrm{\Delta r}}{4}}},$

整理得

$T_P^{t+\mathrm{\Delta t}}=T_P^t+\frac{1}{\rho c_P{r}_P\mathrm{\Delta r\Delta \theta \Delta z}}(\frac{T_N^t-T_P^t}{\frac{\mathrm{\Delta r}}{r_N-\frac{\mathrm{\Delta r}}{4}}+\frac{\mathrm{\Delta r}}{r_P+\frac{\mathrm{\Delta r}}{4}}}+\frac{T_S^t-T_P^t}{\frac{\mathrm{\Delta r}}{r_P-\frac{\mathrm{\Delta r}}{4}}+\frac{\mathrm{\Delta r}}{r_S+\frac{\mathrm{\Delta r}}{4}}})2\mathrm{k\Delta \theta \Delta z\Delta t}+$

$\frac{1}{\rho c_P{r}_P\mathrm{\Delta r\Delta \theta \Delta z}}[(\frac{T_T^t-T_P^t}{r_P\mathrm{\Delta \theta }}-\frac{T_P^t-T_B^t}{r_P\mathrm{\Delta \theta }})\mathrm{k\Delta r\Delta z\Delta t}+\frac{k{r}_P\mathrm{\Delta r\Delta \theta \Delta t}}{\mathrm{\Delta z}}(T_E^t+T_W^t-2T_P^t)]$

$=\frac{\mathrm{k\Delta t}}{\rho c_P{r}_P\mathrm{\Delta r\Delta \theta \Delta z}}[(\frac{T_N^t}{\frac{\mathrm{\Delta r}}{r_N-\frac{\mathrm{\Delta r}}{4}}+\frac{\mathrm{\Delta r}}{r_P+\frac{\mathrm{\Delta r}}{4}}}+\frac{T_S^t}{\frac{\mathrm{\Delta r}}{r_P-\frac{\mathrm{\Delta r}}{4}}+\frac{\mathrm{\Delta r}}{r_S+\frac{\mathrm{\Delta r}}{4}}})2\mathrm{\Delta \theta \Delta z}+\frac{\mathrm{\Delta r\Delta z}}{r_P\mathrm{\Delta \theta }}(T_T^t+T_B^t)$

$+\frac{r_P\mathrm{\Delta r\Delta \theta }}{\mathrm{\Delta z}}(T_E^t+T_W^t)]+\{1-\frac{\mathrm{k\Delta t}}{\rho c_P{r}_P\mathrm{\Delta r\Delta \theta \Delta z}}[(\frac{1}{\frac{\mathrm{\Delta r}}{r_N-\frac{\mathrm{\Delta r}}{4}}+\frac{\mathrm{\Delta r}}{r_P+\frac{\mathrm{\Delta r}}{4}}}+\frac{1}{\frac{\mathrm{\Delta r}}{r_P+\frac{\mathrm{\Delta r}}{4}}+\frac{\mathrm{\Delta r}}{r_S+\frac{\mathrm{\Delta r}}{4}}})2\mathrm{\Delta \theta \Delta z}$

$+\frac{2\mathrm{\Delta r\Delta z}}{r_P\mathrm{\Delta \theta }}+\frac{2r_P\mathrm{\Delta r\Delta \theta }}{\mathrm{\Delta z}}\left]\right\}{T}_p^t---\left(10\right)$

保证稳定性的条件是的系数大于等于0，即

$\mathrm{\Delta t}\le \frac{(\rho c_P/k)r_P\mathrm{\Delta r\Delta \theta \Delta z}}{(\frac{1}{\frac{\mathrm{\Delta r}}{r_N-\frac{\mathrm{\Delta r}}{4}}+\frac{\mathrm{\Delta r}}{r_P+\frac{\mathrm{\Delta r}}{4}}}+\frac{1}{\frac{\mathrm{\Delta r}}{r_P-\frac{\mathrm{\Delta r}}{4}}+\frac{\mathrm{\Delta r}}{r_S+\frac{\mathrm{\Delta r}}{4}}})2\mathrm{\Delta \theta \Delta z}+\frac{2r_P\mathrm{\Delta r\Delta \theta }}{\mathrm{\Delta z}}+\frac{2\mathrm{\Delta r\Delta z}}{r_P\mathrm{\Delta \theta }}}---\left(11\right)$

此式即为步长的稳定性限制条件。

下面对空穴体积变化的计算与处理进行介绍。

本发明为准确模拟容器内伴随着空穴体积变化的相变传热过程， 引入一个变量——空穴体积分数f_v，表示单元内空穴体积的相对大 小。f_v为0时，代表单元内无空穴；f_v为1时，代表单元内被空穴充 满；f_v介于0和1之间时，代表单元有空穴存在。通过引入空穴体积 分数f_v，本发明提出了计算空穴体积变化及空穴调整的算法，建立了 伴随有空穴生成和发展的高温相变蓄热容器的换热过程数学模型，完 成了微重力下伴随着空穴体积变化的相变传热过程的热分析计算。

由于PCM容器内存在空穴，这就要求在对控制方程离散化后， 计算每个单元与相邻单元的换热量时，需要对该单元及相邻单元内是 否存在空穴作出判断，以正确计算该单元内PCM的质量以及与相邻 单元的换热量。

当有空穴存在时，PCM单元内PCM的质量应为

m_i，j，k＝ρ_i，j，kV_i，j，k(1-f_vi，j，k)

i＝1，2，……，II；j＝1，2，……，JJ；k＝1，2，……，KK

                                              (12)

式中i、j、k代表第(i，j，k)个单元，II、JJ、KK代表PCM容器各 坐标方向划分的单元数，f_v为空穴体积分数，具体含义如下：f_v＝0， 无空穴；0＜f_v＜1，单元内存在部分空穴；f_v＝1，单元内全部为空穴。

另外，随着PCM相变的发生，发生相变的PCM的密度随之发生 变化，必须考虑PCM体积变化引起的某些单元内空穴体积的改变， 即f_v随之变化。容器内PCM的体积变化量可由发生相变的PCM的质 量计算得到。每一个时间步长内，发生相变的PCM的质量即为前后 两个时刻PCM容器内固态或液态PCM的质量之差。在相变过程中由 于PCM体积膨胀或收缩，液态PCM会充填容器内的空穴，或者从某 单元中抽出，使得该单元变为含有部分空穴甚至全空，这样每一个时 间步长内计算的液态PCM的质量前后是不一致的，而固态PCM的质 量是不变的。为避免错误的计算结果，本发明采用计算PCM容器内 每一个时间步长的前后两个时刻：t和t+Δt时刻的固态PCM的质量 之差Δm作为发生相变的PCM的质量，即

$\mathrm{\Delta m}=\underset{i,j,k}{\Sigma }[(1-f_{\mathrm{li},j,k}^{t+\mathrm{\Delta t}})·(1-f_{\mathrm{vi},j,k}^{t+\mathrm{\Delta t}})·{\rho }_{i,j,k}^{t+\mathrm{\Delta t}}-(1-f_{\mathrm{li},j,k}^t)·(1-f_{\mathrm{vi},j,k}^t)·{\rho }_{i,j,k}^t]·V_{i,j,k}---\left(13\right)$

式中f_l为PCM液相体积分数。

为保证质量守恒，每一个时间步长内都要检查当前时刻容器内 PCM的总质量，若存在偏差，则对Δm进行偏差修正。判断修正后的 Δm的值，如果Δm＞0，说明固态PCM的质量增加，PCM发生了凝固； 如果Δm＜0，说明固态PCM的质量减少，PCM发生了熔化；如果 Δm＝0，说明固态PCM的质量没有发生变化，PCM正处于显式吸热 或放热状态。如果Δm不为0，说明PCM发生了相变，由前后两个时 刻固态PCM的质量之差Δm，就可以求出每一个时间步长内PCM的 体积变化量ΔV，即

ΔV＝Δm·(1/ρ_s-1/ρ_l)          (14)

式中：ρ_s表示固态PCM的密度，ρ_l表示液态PCM的密度。 根据计算得到的PCM体积变化量ΔV，对PCM容器的空穴体积进行 调整。如果ΔV＞0，说明PCM的体积增加，就减小空穴所占的体积； 如果ΔV＜0，说明PCM的体积缩小，相应地增大空穴体积。因假定 空穴位于外壁处，故减小空穴所占的体积，按从内向外的顺序进行， 而增大空穴体积按相反的顺序进行。

调整空穴体积应遵循的原则是：如果f_v＜1，同时mf_l＞0，即单元 非全空，有液态PCM存在，则可以从中削减PCM的体积，即增加空 穴体积；如果f_v＞0，即单元未充满，则可以向单元内充填液态PCM， 即缩小空穴体积。

下面对空穴有效传热系数的计算与处理进行介绍。

由于PCM凝固时体积收缩以及吸热器设计上的考虑，在PCM容 器内总有一定比例的空穴存在。对LiF来说，空穴占PCM容器的体 积百分比处于9％～21％之间，随PCM的熔化和凝固而变动。PCM和 容器壁由热膨胀而引起的体积变化忽略不计。

微重力下空穴的分布尚无成熟的理论可以引用。在STS-62飞行任 务进行了PCM容器的空间飞行试验(Namkoong David，Jacqmin D and  Szaniszlo.Effect of microgravity on material undergoing melting and  freezing-the TES experiment.AIAA 95-0614)。试验后的PCM容器进行 了层析X射线照相法研究，试验用PCM容器及照相结果如图3所示。 从图中可以看出，凝固后的固态LiF集中在PCM容器靠近辐射排热 器后部的一端(图3中靠近位置9处，PCM容器在该处的温度最低)， 而另一端由于润湿性的作用也聚集了一部分LiF。容器内部形成的空 穴集中位于容器内温度稍高的区域，而且倾向于朝向加热热流最大的 方向(图3中容器正上方)。对于本发明研究的PCM容器来说，由于 热量从内壁被循环工质带走，因此当处于阴影期时，内壁处温度较低， PCM首先在内壁处发生凝固，并沿径向逐渐向外扩展，最终在外壁 处形成空穴。故可以假定空穴位于靠近容器外壁处，形成一柱状环形 空间。

假定空穴内充满PCM蒸气，蒸气压力很低，1121K时LiF蒸气压 力只有0.933Pa，可以认为空穴内基本上充满了LiF蒸气。由于蒸气 压力很低，空穴内蒸气质量可以忽略不计。故透过空穴的换热包括空 穴导热和空穴界面间的辐射换热。空穴内轴向温度梯度很小，因此空 穴内蒸气温度分布按径向稳态传热方程确定

$\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T}{\partial r}\right)=0---\left(15\right)$

方程(15)的解如下

T(r)＝Alnr+B       (16)

式中

$A=\frac{T\left(r_0\right)-T\left(r_v\right)}{\ln (r_0/r_v)},$$B=\frac{T\left(r_0\right)-\ln [T\left(r_0\right)-T\left(r_v\right)]}{\ln (r_0/r_v)}$

下标0代表PCM容器外壁的内表面，v代表与PCM交界处的空穴表 面。

空穴内的辐射换热计算基于以下假设：(1)所有空穴表面均为漫 反射灰体表面；(2)PCM表面吸收所有波长的辐射；(3)空穴内蒸 气不参与辐射换热。在二维和三维情况下，相变界面以及空穴界面的 形状通常是扭曲的、不规则的，界面之间还会互相遮挡，计算空穴界 面之间的辐射换热角系数非常困难。虽然采用蒙特卡洛法、霍特尔 (Hottel)区域法以及离散传播法等方法来计算辐射换热可以得到比 较精确的结果，但这些方法本身都需要大量的计算时间和计算机存储 量，将使得整个相变过程的计算时间大大延长，因此是不经济的。此 外，空穴界面间的辐射换热主要集中在径向，轴向和周向的辐射换热 相比弱得多，为简化计算，忽略轴向和周向辐射换热。

为量化空穴内蒸气导热与空穴界面间辐射换热的相对大小，本发 明对PCM容器径向一维情况下空穴界面辐射换热与LiF蒸气导热进 行了比较。

空穴内LiF蒸气的导热系数k_LiF根据气体运动论按如下公式计算

$k_v=1.457\times {10}^{-3}\sqrt{T}---\left(17\right)$

将空穴两界面间的辐射换热组合到LiF蒸气导热中去，即综合考 虑透过空穴的导热和辐射换热效应，空穴有效传热系数k_veff可按下式 计算：

$k_{\mathrm{veff}}=k_{\mathrm{LiF}}+k_{\mathrm{reff}}=k_{\mathrm{LiF}}+\ln \left(\frac{r_w}{r_v}\right)·\frac{r_v\sigma (T_w\left(t\right)+T_v\left(t\right))(T_w^2\left(t\right)+T_v^2\left(t\right))}{\frac{1}{{ϵ}_{\mathrm{PCM}}}+\frac{r_v}{r_w}(\frac{1}{{ϵ}_w}-1)}---\left(18\right)$

式中k_reff为空穴界面间辐射换热的折合当量导热系数，σ为斯蒂芬 玻尔兹曼常数(σ＝5.67×10^-8W·m^-2·K^-4)，T_W为外壁温度，T_V为PCM 在空穴界面处的温度，ε_PCM为PCM发射率，ε_W为容器壁发射率，v 代表空穴界面。

图4为计算得到的空穴有效传热系数k_veff与PCM蒸气导热系数 k_LiF之比k_veff/k_LiF在整个轨道周期内的变化，可以看出k_veff/k_LiF在 2.29～4.57之间变动，从中可以得到k_veff/k_LiF在整个周期内的算术平均 值为3.407，即空穴辐射换热量平均为空穴导热量的2.407倍。考虑到 PCM容器侧壁的导热作用，在三维计算中按比较保守的估计，取空 穴辐射换热量为LiF蒸气导热量的2倍，即空穴复合传热系数 k_eff＝3.0*k_LiF。

图5示出了本发明的微重力下伴有空穴的高温相变换热计算程序 框图。

附图说明

图1为任一单元节点P(i，j，k)的控制容积示意图。

图2为节点P(i，j，k)所在的轴向截面示意图。

图3为试验用PCM容器及照相结果。

图4为计算得到的空穴有效传热系数k_veff与PCM蒸气导热系数 k_LiF之比k_veff/k_LiF在整个轨道周期内的变化。

图5所示为微重力下伴有空穴的高温相变换热计算程序框图。

图6为PCM容器的几何外形示意图。

图7为PCM容器的轴向截面示意图。

图8为本发明计算稳定后结果与美国NASA方案中高温相变蓄热 单元热循环计算结果比较。

具体实施方式

为了更好地理解本发明的技术方案，以下结合附图以及具体实例 作进一步描述，实施例是对本发明技术特征的支持，而不是限定。

本发明涉及的物理模型系美国宇航局(NASA，National  Aeronautics and Space Administration)所研究的空间站太阳能热动力发 电系统热管吸热器高温相变蓄热单元。PCM容器外形如图6所示， 容器的轴向截面如图7所示，PCM蓄热容器外环直径89mm，内环 直径51mm，内壁、外壁、侧壁厚度1.5mm，壁面材料均采用 Hayness188，轴向长度为25mm，相变材料为LiF，轨道周期为日照 期54min，阴影期36min。

由于该实验物理模型简单、数据详尽且具有代表性，本发明对该 NASA方案进行数值计算，并将计算结果与NASA计算结果进行比较， 以验证其准确性，具有很强的说服力。图5所示为微重力下伴有空穴 的高温相变换热计算程序框图。其主要实现步骤如下：

(1)划分计算区域、连续体、边界及网格；

(2)设定操作条件：参考压力为一个大气压；

(3)定义层流、非稳态熔化/凝固模型，离散坐标辐射模型；

(4)在“定义”一栏中选择“材料”用来设定流场的介质：定义 流场内的流质为相变材料LiF，壁面材料为Hayness188，定义空穴里 为稀薄LiF蒸汽；

(5)建立基于能量方程的偏微分控制方程组(参见表达式(1))；

(6)定义边界条件：定义高温蓄热容器的内壁为热流边界，外 壁和侧壁为绝热边界，其它壁面为流/固耦合边界；其中，所述热流 边界为周期性热流边界，写入周期性热流自定义函数，日照期54min， 阴影期36min；日照期热流密度为q_sun＝15424-24.88(T_wall-1020)，W/m²； 阴影期热流密度为q_shadow＝-12314-24.88(T_wall-1020)，W/m²；

(7)初始化，定义初始温度为950K，初始速度为0m/s；

(8)监视内壁面温度变化，相变材料区温度及液相分数变化；

(9)对以上包含边界条件的已封闭的偏微分控制方程组进行离 散(参见表达式(8)、(9))，并对其迭代求解；

(10)对整个计算区域初始化，迭代时间步长设为60s即可，迭 代精度达到10e-6，便可以认为计算的结果已接近实际效果；

(11)通过上述求解过程得出了各计算节点上的解后，需要通过 线值图、矢量图、等直线图、流线图、云图等方式对计算结果进行表 示。

图8所示为本发明计算稳定后结果与美国NASA方案中高温相变 蓄热单元热循环计算结果比较。通过与NASA方案计算结果比较后可 以看出：本发明数值计算所得一个周期内热管壁面、蓄热容器外壁面 温度变化范围和趋势与NASA方案计算所得到的温度变化范围和趋 势比较接近，对应处各点的温度值和NASA方案计算值最多相差不超 过16K，这验证了本发明所建立的微重力条件下伴有空穴的高温相 变蓄热单元计算方法、模型及程序的合理性与准确性。