CN102495968A - 重力条件下高温蓄热容器内固/液相变的数值模拟方法 - Google Patents

Abstract

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

Description

重力条件下高温蓄热容器内固/液相变的数值模拟方法

技术领域

本发明涉及一种计算流体力学领域的方法，具体是一种在重力条件下对高温蓄热容器内固/液相变的数值模拟方法，即通过对附加自然对流条件下的焓法能量控制方程组进行计算。

技术背景

固液相变问题，即有熔化或凝固峰面运动的传热问题，为纪念德国科学家J.Stefan研究极地冰层的厚度问题又被称为斯蒂芬(Stefen)问题。从1860年Franz Neumann研究半无限大物体的相变问题至今，人们对固液相变问题的研究已经有一百多年的历史。由于该类问题的非线性特征，给问题的求解带来了很大的困难。

目前只有极少数简单情形下的相变问题能够得到解析解。1860年Neumann给出了半无限大区域凝固问题的精确解。其它许多情况下的精确解都是基于Neumann解得到的。存在精确解析解的情况主要集中在常物性、具有简单边界条件和初始条件的一维无限大和半无限大区域，对于一维有限区域及多维情况很难得到解析解。

对于有限区域内的一维相变问题，人们发展了近似解法。近似解法主要有积分法、准稳态法、摄动法、热阻法和级数展开法等。近似解法的应用主要局限于一维单一单调界面相变问题，在极少数特别简单的边界条件下也可推广到二维情况，但对于类似本发明具有复杂边界条件以及存在非单调、多个界面的多维相变问题就显得无能为力了。

对于复杂条件、多维相变问题，数值解法几乎是唯一可行的强有力手段。

相变导热问题的数值解法可以分为两大类。一类称为界面跟踪法或强数值解法，另一类称为固定网格法或弱数值解法。

界面跟踪法包括固定步长法、变空间步长法、变时间步长法、自变量变换法、贴体坐标法和等温面移动法等。界面跟踪法在每一个时间步长都要确定固液两相界面的位置和温度分布，但是通常固液两相界面的形状是不规则的，而且其位置不论在空间还是时间上都是未知的，因此在进行离散化数值求解时需要用特殊的插值方法或采用坐标变换把不规则两相界面变为固定边界。采用坐标变换的方法虽然可以把不规则的移动界面变成一个形状简单的固定界面，但使得原始方程复杂化，这类解法适用于一维问题的求解，但推广到二维和三维情况时方程形式将变得极为复杂，而且对存在非单调、多个界面的情形也是不适用的。

固定网格法或弱数值解法，不需要跟踪固液两相界面的位置，把包含不同相态的求解区域作为整体求解，因而具有很大的灵活性，很容易推广到多维、多界面情况，这类解法有显热容法和焓法。

显热容法把物质的相变潜热看作是在一个很小的温度范围内有一个很大的显热容，从而把分区描述的相变问题转变为单一区域上的非线性导热问题，达到整体求解的目的。显热容法的缺点是当相变温度范围很窄时，如果时间步长稍大，计算过程中就会越过相变区，导致忽略了相变潜热，造成计算结果失真。而对于在单一温度下发生的相变过程，其缺点就更加突出。

焓法是将热焓和温度一起作为待求函数，在整个区域(包括液相、固相和两相界面)建立一个统一的能量方程，利用数值方法求出热焓分布，然后确定两相界面。焓法没有显热容法的缺点，具有方法简单、灵活方便、容易扩展到多维情况等优点，能够求解具有复杂边界条件以及非单调、多个界面的相变问题，但焓法也有不足之处。当网格划分较粗时，相变界面的位置无法清晰地给出，另外，糊态区的物性参数只能按固液相的比例近似计算得到，存在一定的不确定性。当然，前一个缺点可以通过加密网格来克服，后一个缺点将来随着人们对相变过程认识的加深也会逐步得到解决。这两个缺点与焓法的巨大优点相比是微不足道的。因此，目前焓法得到了广泛的应用。对于具有复杂边界条件、相变界面形状复杂以及相变过程中可能存在多个相变界面的相变材料(PCM)容器内的换热过程而言，采用弱数值解法是唯一可行的方法。

经对现有技术的文献检索发现，徐伟强等在《太阳能学报》2005年第3期上发表的论文“热管式吸热器单元热管传热的数值模拟分析”，对微重力条件下固/液相变换热过程进行了分析，但只考虑了导热，不涉及对流及辐射问题。该方法只适用于进行简单科学研究，不能应用于工程实际。

发明内容

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

为实现上述目标，本发明采用有限容积法对控制体进行积分，采用焓法来处理相变过程，并运用BOUSSINESQ假设(即：忽略粘性耗散；除密度外其他物性为常数；对密度仅考虑动量方程中与体积力有关的项，其余各项中的密度仍为常数)以模拟附加有自然对流的相变传热，包括：(1)定义操作条件为重力条件下；(2)定义第二类传热边界条件即热流边界；(3)定义熔化和凝固模型。

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

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

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

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

(3)定义操作条件为重力条件，定义基于BOUSSINESQ假设的自然对流模型；

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

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

a)连续性方程

$\frac{u\∂}{\∂z}+\frac{1}{r}\frac{\∂\left(\mathrm{rv}\right)}{\∂r}=0---\left(1\right)$

式中：u，v分别为z方向、r方向的速度分量；

b)动量方程

$\frac{\∂u}{\∂t}+\frac{{\∂u}^2}{\∂z}+\frac{1}{r}\frac{\∂\left(\mathrm{ruv}\right)}{\∂r}=-\frac{1}{\mathrm{\ρ}}\frac{\∂P}{\∂z}+\mathrm{g\β}(T-T_m)+v{\▿}^{2}u---\left(2\right)$

$\frac{\∂v}{\∂t}+\frac{\∂\left(uv\right)}{\∂z}+\frac{1}{r}\frac{\∂\left({\mathrm{rv}}^2\right)}{\∂r}=-\frac{1}{\mathrm{\ρ}}\frac{\∂P}{\∂r}+v({\▿}^{2}v-\left(\frac{v}{r^2}\right))$

式中：P为压力，g为重力加速度，β为相变材料(PCM)液态体积膨胀系数，v为运动粘度。

c)能量方程

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

式中：e为比焓，以比焓的形式把相变的影响考虑进去，适用于整个计算区域。其中焓与温度的关系式可表示成如下的形式：

$T=\left\{\begin{array}{cc}{T}_m+e/c,& e\&le;0\\ T_m,& 0<e<{\mathrm{\&Delta;H}}_m\\ T_m+(e-{\mathrm{\&Delta;H}}_m)/c,& e\&GreaterEqual;{\mathrm{\&Delta;H}}_m\end{array}\right.$

式中：e为比焓，ρ为密度，t为时间，k为导热系数，c为比热容，T为温度，T_m为相变温度，ΔH_m为物质单位质量的相变潜热；

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

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

q_sun＝15424-24.88(T_wall-1020)

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

q_shadow＝-12314-24.88(T_wall-1020)

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

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

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

(9)采用有限容积法、交错网格和乘方格式差分格式对步骤(5)中的偏微分控制方程进行离散化得到代数方程组，并利用步骤(6)的边界条件和步骤(7)的初始条件进行封闭和求解；

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

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

进一步地，在上述步骤(9)中，计算每个单元与相邻单元的换热量时，需要对该单元及相邻单元内是否存在空穴作出判断，以正确计算该单元内相变材料(PCM)的质量以及与相邻单元的换热量，当有空穴存在时，相变材料(PCM)单元内相变材料(PCM)的质量为

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

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

式中i、j、k代表第(i，j，k)个单元，II、JJ、KK代表相变材料(PCM)容器各坐标方向划分的单元数，m代表相变材料(PCM)的质量，ρ代表相变材料(PCM)的密度，V为单元的体积，f_v为空穴体积分数，f_v＝0，无空穴；0＜f_v＜1，单元内存在部分空穴；f_v＝1，单元内全部为空穴。

进一步地，还需要考虑相变材料(PCM)体积变化引起的某些单元内空穴体积的改变，由前后两个时刻固态相变材料(PCM)的质量之差Δm求出每一个时间步长内相变材料(PCM)的体积变化量ΔV，根据体积变化量ΔV对相变材料(PCM)容器的空穴体积进行调整：如果ΔV＞0，说明相变材料(PCM)的体积增加，就减小空穴所占的体积；如果ΔV＜0，说明相变材料(PCM)的体积缩小，相应地增大空穴体积。

进一步地，在步骤(10)中，在求解离散化后的能量方程时采用显式格式，在求解离散化后的动量方程时采用交替方向块迭代法(ADI)与块追赶法(TDMA)相结合的计算方法。

进一步地，采用求解压力耦合方程的半隐方法(SIMPLE算法)求解流场，为防止迭代过程发散，对速度u、v以及压力修正进行亚松弛。

以上是本发明为解决其技术问题所采取的技术方案，本发明的数值模拟方法通过基于焓法和有限控制体积法的计算程序，使得重力条件下高温蓄热容器内的熔化/凝固这一复杂的相变问题得到了简便、高效的解决，具有重要的实用价值。

下面对本发明的上述技术方案中的步骤(9)，即对偏微分控制方程的离散化方法以及代数方程的求解方法进行介绍。

偏微分控制方程的离散方法有泰勒级数展开法、多项式拟合法、控制容积积分法、平衡法以及变分法等。控制容积积分法又称有限容积法，是传热数值计算中广泛采用的方法，该方法的推导过程物理概念清晰，推导结果具有明确的物理意义，离散方程的守恒特性可以得到保证。本发明采用控制容积积分法导出离散方程。

在求解流体流动问题时，为了检测不合理的压力场，广泛采用交错网格。所谓交错网格就是指把速度u、v及压力p(包括其它所有标量和物性参数)分别存储于三套不同网格上的网格系统。其中速度u存在于压力控制容积的东西界面上，速度v存在于压力控制容积的南北界面上，u、v各自的控制容积则以速度所在位置为中心。如图1所示，u控制容积与主控制容积(即压力的控制容积)之间在x方向有半个网格的错位，v控制容积与主控制容积(即压力的控制容积)之间在y方向有半个网格的错位。在交错网格系统中，相邻两点间的压力差恰好构成了两网格点之间速度分量的自然驱动力，从而可以检测出不合理的压力场。但是，采用交错网格也要付出一定的代价。首先增加了计算工作量，在求解u、v方程时必须通过插值计算控制容积各界面上的物性参数和流量，在求解位于主节点上的变量方程时控制容积各界面上的速度u、v也要进行插值计算；其次，三套网格系统(二维情况)中节点编号必须仔细处理才能协调一致，增加了程序编制工作量。

a)通用控制方程的离散

把通用控制方程写成圆柱轴对称坐标下的通用形式

$\frac{\&PartialD;\left(\mathrm{\&rho;\&phi;}\right)}{\&PartialD;t}+\frac{\&PartialD;}{\&PartialD;z}\left(\mathrm{\&rho;u\&phi;}\right)+\frac{1}{r}\frac{\&PartialD;}{\&PartialD;r}\left(\mathrm{r\&rho;v\&phi;}\right)=\frac{\&PartialD;}{\&PartialD;z}\left({\mathrm{\&Gamma;}}_{\mathrm{\&phi;}}\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;z}\right)+\frac{1}{r}\frac{\&PartialD;}{\&PartialD;r}\left({\mathrm{r\&Gamma;}}_{\mathrm{\&phi;}}\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;r}\right)+S_{\mathrm{\&phi;}}---\left(4\right)$

式中，

是通用变量，Γ_φ与

是与

相对应的广义扩散系数及广义源项。

引入在z及r方向的对流-扩散总通量J_z、J_r，即

$J_z=\mathrm{\&rho;u\&phi;}-\mathrm{\&Gamma;}\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;z}$ $J_r=\mathrm{r\&rho;v\&phi;}-\mathrm{r\&Gamma;}\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;r}$

(5)

则式(4)可改写为

$\frac{\&PartialD;\left(\mathrm{\&rho;\&phi;}\right)}{\&PartialD;t}+\frac{{\&PartialD;J}_z}{\&PartialD;z}+\frac{1}{r}\frac{{\&PartialD;J}_r}{\&PartialD;r}=S---\left(6\right)$

对控制容积在时间和空间上作积分，方程左边第一项为

${\&Integral;}_t^{t+\mathrm{\&Delta;t}}{\&Integral;}_s^n{\&Integral;}_w^e\frac{\&PartialD;\left(\mathrm{\&rho;\&phi;}\right)}{\&PartialD;t}\mathrm{rdtdzdr}=[{\left(\mathrm{\&rho;\&phi;}\right)}_P-{\left(\mathrm{\&rho;\&phi;}\right)}_P^0]r_p\mathrm{\&Delta;z\&Delta;r}---\left(7\right)$

假设z、r方向的总通量J_z、J_r在各自的界面w、e与s、n上是均匀的，则有

${\&Integral;}_t^{t+\mathrm{\&Delta;t}}{\&Integral;}_s^n{\&Integral;}_w^e\frac{{\&PartialD;J}_z}{\&PartialD;z}\mathrm{rdtdzdr}={\&Integral;}_s^n(J_z^e-J_z^w){\mathrm{\&Delta;tr}}_p\mathrm{dr}=(J_z^e-J_z^w){\mathrm{\&Delta;tr}}_p\mathrm{\&Delta;r}=(J_e-J_w)\mathrm{\&Delta;t}---\left(8\right)$

式中

代表z方向上界面e、w处单位面积上的通量，J_e、J_w代表总面积上的通量。同理可得

${\&Integral;}_t^{t+\mathrm{\&Delta;t}}{\&Integral;}_s^n{\&Integral;}_w^e\frac{1}{r}\frac{{\&PartialD;J}_r}{\&PartialD;r}\mathrm{rdtdzdr}=(J_n-J_s)\mathrm{\&Delta;t}---\left(9\right)$

取S＝S_C+S_Pφ_P(S_P≤0)，将式(6)积分后所得各项整理得

$\frac{{\left(\mathrm{\&rho;\&phi;}\right)}_P-{\left(\mathrm{\&rho;\&phi;}\right)}_P^0}{\mathrm{\&Delta;t}}{\mathrm{\&Delta;zr}}_p\mathrm{\&Delta;r}+(J_e-J_w)+(J_e-J_w)+(J_n-J_s)=(S_C+S_P{\mathrm{\&phi;}}_P){\mathrm{\&Delta;zr}}_p\mathrm{\&Delta;r}---\left(10\right)$

J_e＝(a_E+F_e)φ_P-a_Eφ_E

(11)

J_e＝(a_E+F_e)φ_P-a_Eφ_E

(12)

J_e＝(a_E+F_e)φ_P-a_Eφ_E

(13)

J_e＝(a_E+F_e)φ_P-a_Eφ_E

(14)

将以上各式代入式(10)，整理得

a_Pφ_P＝a_Eφ_E+a_Wφ_W+a_Nφ_N+a_Sφ_S+b

(15)

其中

a_E＝D_eA(P_Δe)＝D_eA(|P_Δe|)+[|-F_e，0|]

(16)

a_W＝D_wA(P_Δw)＝D_wA(|P_Δw|)+[|F_w，0|]

(17)

a_N＝D_nA(P_Δn)＝D_nA(|P_Δn|)+[|-F_n，0|]

(18)

a_S＝D_sA(P_Δs)＝D_sA(|P_Δs|)+[|F_s，0|]

(19)

$b=S_C{r}_p\mathrm{\&Delta;z\&Delta;r}+a_P^0{\mathrm{\&phi;}}_P^0---\left(20\right)$

$a_P=a_E+a_W+a_N+a_S+a_P^0-S_P{\mathrm{\&Delta;zr}}_p\mathrm{\&Delta;r}---\left(21\right)$

$a_P^0=\frac{{\mathrm{\&rho;}}_P{\mathrm{\&Delta;zr}}_p\mathrm{\&Delta;r}}{\mathrm{\&Delta;t}}---\left(22\right)$

F_e＝(ρu)_er_pΔz ，F_w＝(ρu)_wr_pΔz，F_n＝(ρv)_nr_nΔr，F_s＝(ρv)_sr_sΔr

(23)

$D_e=\frac{{\mathrm{\&Gamma;}}_e{r}_p\mathrm{\&Delta;z}}{{\left(\mathrm{\&delta;z}\right)}_e},$ $D_w=\frac{{\mathrm{\&Gamma;}}_w{r}_p\mathrm{\&Delta;z}}{{\left(\mathrm{\&delta;z}\right)}_w},$ $D_n=\frac{{\mathrm{\&Gamma;}}_n{r}_n\mathrm{\&Delta;r}}{{\left(\mathrm{\&delta;r}\right)}_n},$ $D_s=\frac{{\mathrm{\&Gamma;}}_s{r}_s\mathrm{\&Delta;r}}{{\left(\mathrm{\&delta;r}\right)}_s}$

(24)

其中F表示对流(流动)的强度，D表示扩散的强度，其比值P_Δ＝F/D称为贝克利数，对流动来说，

对于换热，

式(15)可简写为

a_Pφ_P＝∑a_nbφ_nb+b(25)

下标nb代表相邻节点。

b)动量方程的离散

在交错网格中，一般变量

的离散过程与上一节中相同。但对动量方程而言，积分用的控制容积是速度分量u、v各自的控制容积，压力梯度项从源项中分离出来，对于速度u_e的控制容积，该项积分为

${\&Integral;}_s^n{\&Integral;}_P^E(-\frac{\&PartialD;p}{\&PartialD;z})\mathrm{rdzdr}=-{\&Integral;}_s^n\left(p{|}_P^E\right)r_P\mathrm{dr}\&ap;(p_P-p_E)r_P\mathrm{\&Delta;r}---\left(26\right)$

u_e的控制容积的东、西界面上压力是各自均匀的，分别为p_E和p_P。关于u_e的离散方程便具有以下形式

a_eu_e＝∑a_nbu_nb+b+(p_P-p_E)A_e

(27)

式中：u_nb是u_e的邻点速度(u_ee，u_n，u_w以及u_s)；b为不包括压力在内的源项中的常数部分，对非稳态问题为

A_e＝r_PΔr，是压力差的作用面积；系数a_nb的计算公式取决于所采用的格式，见上一节所述。

同理可得速度v的离散方程

a_nv_n＝∑a_nbv_nb+b+(p_P-p_N)A_n

(28)

c)对流-扩散方程的差分格式

对流-扩散方程的差分格式有五种三点格式可以选择，即中心差分、迎风差分、混合格式、指数格式、乘方格式(见表1)。其中中心差分在P_Δ大于一定的数值后会出现解的振荡现象。迎风差分又称为上风差分、施主格子差分，它克服了中心差分的缺点，保证了动量离散方程中各系数永远大于或等于零，避免了解的振荡现象的产生。混合格式综合吸取了中心差分和迎风差分的优点，而且克服了中心差分当P_Δ＞2时出现解的振荡、迎风差分对扩散项的处理不考虑P_Δ的影响等缺点。指数格式在五种三点格式中最精确，但指数的计算比较费时。1979年Patankar提出了乘方格式，在很大的范围内，乘方格式与指数格式的结果非常接近，而计算工作量减少了。本发明中对流-扩散方程的离散采用乘方格式。

表1五种三点格式的A(|P_Δ|)

格式	A(|PΔ|)
		中心差分	1-0.5|PΔ|
迎风差分	1
		混合格式	[|0，1-0.5|PΔ||]
指数格式	|PΔ|/(exp(|PΔ|)-1)
		乘方格式	[|0，(1-0.1|PΔ|)5|]

d)能量方程的离散

采用容积平衡法，得内外节点的离散差分方程与前面无对流情况时的相似，但增加了与相邻单元的对流换热项。为保证计算稳定，对流项采用迎风差分格式，得到通过节点P所在控制容积界面的对流换热量如下

w界面

e界面

s界面

n界面

其中，ρ_w、ρ_e、ρ_s、ρ_n分别为界面w、e、s和n上的流体密度，由界面两侧单元的密度插值计算得到。F′_w、F′_e、F′_s、F′_n分别为界面w、e、s和n的流通面积，当考虑整个圆周时分别为2πr_iΔr、2πr_iΔr、2π(r_i-Δr/2)Δz、2π(r_i+Δr/2)Δz。

参考式(13)，得相变材料(PCM)容器内部能量方程的离散化方程为：

${\mathrm{\&rho;}}_P(e_P^{t+\mathrm{\&Delta;t}}-e_P^t)r_P\mathrm{\&Delta;r\&Delta;z}=(k_e\frac{T_E^t-T_P^t}{{\mathrm{\&Delta;z}}_e}-k_w\frac{T_P^t-T_W^t}{{\mathrm{\&Delta;z}}_w})\mathrm{\&Delta;r\&Delta;t}$

$+[{\left(\frac{\mathrm{kr}}{\mathrm{\&delta;}}\right)}_n(T_N^t-T_P^t)-{\left(\frac{\mathrm{kr}}{\mathrm{\&delta;}}\right)}_s(T_P^t-T_S^t)]\mathrm{\&Delta;z\&Delta;t}$

$+{\mathrm{\&rho;}}_w{F}_w^{\&prime;}[e_W^t\max (u_P^t,0)-e_P^t\max ({-u}_P^t,0)]+{\mathrm{\&rho;}}_e{F}_e^{\&prime;}[e_E^t\max (-u_E^t,0)-e_P^t\max (u_P^t,0)]$

$+{\mathrm{\&rho;}}_s{F}_s^{\&prime;}[e_S^t\max (v_P^t,0)-e_P^t\max ({-v}_P^t,0)]+{\mathrm{\&rho;}}_n{F}_n^{\&prime;}[e_N^t\max ({-v}_N^t,0)-e_P^t\max (v_N^t,0)]---\left(33\right)$

对于容器壁上的节点，由于采用壁面无滑移边界条件，壁面处流通量为0，故容器壁上节点离散控制方程不包括对流换热项。

e)空穴体积变化的计算与处理

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

当有空穴存在时，相变材料(PCM)单元内相变材料(PCM)的质量应为

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

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

(34)

式中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)的质量之差作为发生相变的相变材料(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/ρ₁)

(36)

根据计算得到的相变材料(PCM)体积变化量ΔV，对相变材料(PCM)容器的空穴体积进行调整。如果ΔV＞0，说明相变材料(PCM)的体积增加，就减小空穴所占的体积；如果ΔV＜0，说明相变材料(PCM)的体积缩小，相应地增大空穴体积。因假定空穴位于外壁处，故减小空穴所占的体积，按从内向外的顺序进行，而增大空穴体积按相反的顺序进行。

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

f)代数方程的求解方法

最初的微分方程经离散化后得到上述一系列的代数方程，问题的求解转化为对上述代数方程的求解。代数方程的求解可分为直接解法和迭代法两大类。最基本的直接解法有克莱姆法则、Gauss消去法等，但只适于求解未知数个数极少的情形，对于通常划分为大量单元具有很多个节点的离散方程来说，计算量非常大，实际上是不可行的。

描写传热和流体流动问题的控制方程大多数是非线性的。对非线性问题，离散方程中的系数可能都是未知量的函数。因此，通常采用迭代法进行传热和流体流动问题的求解。

比较常用的迭代解法有点迭代法、块迭代法、交替方向块迭代法(ADI-Alternating Direction Implicit Method)等。点迭代法又叫显式迭代法，如雅可比迭代法、高斯-塞德尔迭代法、超/欠松弛迭代法(SOR/SUR)等。块迭代法把求解区域分成若干块，每块由一条或数条网格线组成。在同一块中各节点的值用代数方程的直接解法求得。采用块迭代法可使迭代次数显著减少，缩短计算时间。交替方向块迭代法在迭代计算过程中扫描的方向交替进行，即先逐行进行一次扫描，再逐列进行一次扫描，两次扫描组成一轮迭代。由于边界条件传播的速度加快，加快了迭代的收敛。实际求解过程中可能采用上述几种方法的组合。

本发明把能量方程和动量方程分开求解，在求解能量方程时采用显式格式，在求解动量方程时采用ADI与TDMA(求解一维问题的追赶法)相结合的方法。例如，对通用方程的离散化方程(15)，采用ADI与TDMA相结合的方法，把方程右边的a_Wφ_W和a_Eφ_E两项归并到常数项b中，用b′表示，则方程(15)化为

a_Pφ_P＝a_Nφ_N+a_Sφ_S+b′

(37)

式中b′＝a_Eφ_E+a_Wφ_W+b。由上式，采用TDMA算法即可求出同一条线上的各节点待求变量值φ_P、φ_N以及φ_S，完成该条线上的求解后，再求解下一条线上的变量值，如此扫描完整个求解区域，即完成了一轮迭代。再结合ADI方法，让相邻两轮迭代的扫描方向交替变化，先按列(或行)扫描完成一轮计算，再按行(或列)扫描进行计算，两次扫描完成一轮迭代，从而加快迭代收敛的速度。

为防止迭代求解过程发散，在求解过程中一般采用欠松弛处理，即

${\mathrm{\&phi;}}_P^n={\mathrm{\&phi;}}_P^{n-1}+\mathrm{\&alpha;}({\mathrm{\&phi;}}_P^{*}-{\mathrm{\&phi;}}_P^{n-1})---\left(38\right)$

上标n代表迭代轮数，n-1为上一轮迭代。

为本轮计算得到的未经欠松弛处理的值，

α为松弛因子，把松弛因子组合到迭代计算式中，得

$\left(\frac{a_P}{\mathrm{\&alpha;}}\right){\mathrm{\&phi;}}_P^n=a_N{\mathrm{\&phi;}}_N^n+a_S{\mathrm{\&phi;}}_S^n+b^{\&prime;}+\left(\frac{1-\mathrm{\&alpha;}}{\mathrm{\&alpha;}}\right)a_P{\mathrm{\&phi;}}_P^{n-1}---\left(39\right)$

在进行另一个方向的扫描求解时，只需把上式与b′中的下标互换即可。

g)求解N-S方程的SIMPLE算法

以速度和压力为原始变量的N-S方程求解的关键，是如何求解压力场，或者在假定了一个压力场之后如何改进这个压力场。目前广泛采用的是SIMPLE算法及其衍化派生出来的各种算法。

SIMPLE算法是由Patankar与Spalding在1972年提出的，是Semi-Implicit Method for Pressurelinked Equations的缩写，即求解压力耦合方程的半隐方法。其基本思想是按给定的压力场，求出速度u和v，这样按给定的压力场求得的速度场可能不满足质量守恒的要求，要对给定的压力场进行修正。将求得的速度值代入连续方程，根据速度与压力的关系，得到压力修正方程，由此方程求得压力修正值，用压力修正值来修正压力和速度，以修正后的速度重新计算动量离散方程的系数，开始下一层次的计算。如此反复，直到获得收敛的解。下面首先构造出压力修正方程。

假定原来的压力为P^*，相应的速度为u^*、v^*。记压力修正值为P′，相应的速度修正值为u′、v′，则修正后的压力和速度分别为

P＝P^*+P′，u＝u^*+u′，v＝v^*+v′。

将修正后的压力和速度代入到动量离散方程中，可得

其中u^*、v^*是根据P^*从动量离散方程求解得到的，因此它们满足动量方程，即

$a_e{u}_e^{*}=\mathrm{\&Sigma;}{a}_{\mathrm{nb}}{u}_{\mathrm{nb}}^{*}+b+(p_P^{*}-p_E^{*})A_e---\left(41\right)$

式(40)与式(41)两式相减，得

a_eu′_e＝∑a_nbu′_nb+b+(p′_P-p′_E)A_e

(42)

上式右边第一项为邻点速度修正值对该节点速度修正值的影响，第二项为相邻两节点间的压力修正值之差对该节点速度修正值的影响。忽略邻点速度修正值的影响，只考虑压力修正值的影响，得速度修正方程

a_eu′_e＝(p′_P-p′_E)A_e

(43)

即

$u_e^{\&prime;}=\frac{A_e}{a_e}(p_P^{\&prime;}-p_E^{\&prime;})=d_e(p_P^{\&prime;}-p_E^{\&prime;})---\left(44\right)$

同理可得速度v的修正方程

v′_n＝d_n(p′_P-p′_N)

(45)

式中 $d_e=\frac{A_e}{a_e},$ $d_n=\frac{A_n}{a_n}$

利用速度修正值来修正速度，即

$u_e=u_e^{*}+d_e(p_P^{\&prime;}-p_E^{\&prime;})---\left(46\right)$

$v_n=v_n^{*}+d_n(p_P^{\&prime;}-p_N^{\&prime;})---\left(47\right)$

修正后的速度场应能够满足连续性方程。首先对连续性方程在时间间隔Δt内对主控制容积作积分，得

$\frac{{\mathrm{\&rho;}}_P-{\mathrm{\&rho;}}_P^0}{\mathrm{\&Delta;t}}{\mathrm{\&Delta;zr}}_P\mathrm{\&Delta;r}+[{\left(\mathrm{\&rho;u}\right)}_e-{\left(\mathrm{\&rho;u}\right)}_w]r_P\mathrm{\&Delta;r}+[{\left(\mathrm{\&rho;v}\right)}_n{r}_n-{\left(\mathrm{\&rho;v}\right)}_s{r}_s]\mathrm{\&Delta;z}=0---\left(48\right)$

将速度修正关系式代入上式，得到压力修正值方程

a_Pp′＝a_Ep′_E+a_Wp′_W+a_Np′_N+a_Sp′_S+b

(49)

其中

a_E＝ρ_ed_er_PΔr，a_W＝ρ_wd_wr_PΔr，a_N＝ρ_nd_nr_nΔz，a_S＝ρ_sd_sr_sΔz，

a_P＝a_E+a_W+a_N+a_S，

(50)

$b=\frac{{\mathrm{\&rho;}}_P-{\mathrm{\&rho;}}_P^0}{\mathrm{\&Delta;t}}{\mathrm{\&Delta;zr}}_P\mathrm{\&Delta;r}+[{\left({\mathrm{\&rho;u}}^{*}\right)}_e-{\left({\mathrm{\&rho;u}}^{*}\right)}_w]r_P\mathrm{\&Delta;r}+[{\left({\mathrm{\&rho;v}}^{*}\right)}_n{r}_n-{\left({\mathrm{\&rho;v}}^{*}\right)}_s{r}_s]\mathrm{\&Delta;z}---\left(51\right)$

b的值代表了一个控制容积不满足连续性条件的剩余质量的大小，故可以用各控制容积剩余质量b的绝对值最大值以及各控制容积剩余质量b的代数和作为收敛判据，当这两个值都为很小的值时，就可以认为速度场迭代已经收敛。

采用SIMPLE算法求解N-S方程的步骤如下：

1)先假定一个速度分布，以此计算动量离散方程的系数和常数项；

2)假定一个压力场；

3)依次求解动量方程，得到速度场；

4)求解压力修正值方程；

5)根据压力修正值修正压力和速度；

6)利用改进后的速度重新计算动量离散方程的系数，并用改进后的压力场作为下一层次迭代计算的初值。重复步骤3)～6)，直到获得收敛的解。

在采用SIMPLE算法求解流场时，为防止迭代过程发散，一般对速度u、v以及压力修正要进行亚松弛。不同的具体问题，松弛因子取值不一定相同，应根据具体情况进行适当调整。Patankar推荐的速度松弛因子为0.5，压力松弛因子为0.8。

图2所示为重力条件下的高温相变换热计算程序框图。

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

附图说明

图1为交错网格示意图。

图2为重力条件下的高温相变换热计算程序框图。

图3为相变材料(PCM)高温蓄热容器物理模型示意图。

图4为采用本发明的数值模拟方法得到的计算结果和日本国家航空实验室热管吸热器蓄热单元热循环实验结果比较。

具体实施方式

下面结合附图和具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图3所示，本发明涉及的物理模型系日本国家航空实验室Takeshi等人所研究的空间站太阳能热动力发电系统热管吸热器高温蓄热容器。该高温蓄热容器外环直径89mm，内环直径51mm，内壁、外壁、侧壁厚度1.5mm，壁面材料均采用C/C复合功能梯度材料，其热导率达到400W/(m.K)，轴向长度为25mm，相变材料(PCM)为80.5LiF-19.5CaF₂，轨道周期为日照期40min，阴影期50min。

由于该实验物理模型简单、数据详尽且具有代表性，本发明对该实验进行数值计算，并将计算结果与实验值比较，以验证其准确性，具有很强的说服力。其主要实现步骤如下：

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

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

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

(4)在“定义”一栏中选择“材料”用来设定流场的介质：定义流场内的流质为相变材料(PCM)80.5LiF-19.5CaF₂，壁面材料为C/C复合功能梯度材料，定义空穴里为稀薄80.5LiF-19.5CaF₂蒸汽；

(5)建立包括偏微分方程(1)(3)在内的偏微分控制方程组；

(6)定义蓄热容器内壁为热流边界，外壁和侧壁为绝热边界，其它壁面为流固耦合边界；写入周期性热流自定义函数，日照期40min，阴影期50min；

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

q_sun＝15424-24.88(T_wall-1020)

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

q_shadow＝-12314-24.88(T_wall-1020)

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

(7)对整个计算区域进行初始化，定义初始温度为900K，初始速度为0m/s；

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

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

(10)设定迭代时间步长为60s，迭代精度达到10e-6，对计算区域内的代数方程组反复进行迭代计算，直到满足所设定的迭代精度为止，便可以认为计算的结果已接近实际效果。

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

图4所示为本发明计算稳定后结果与日本国家航空实验室热管吸热器蓄热单元热循环实验结果比较。通过与实验结果比较后可以看出：本发明数值计算所得热管壁温和相变材料(PCM)外壁面温度与日本国家航空实验室对应处各点的实验温度值相差较小，不过本发明计算所得的温度值稍许偏高，但最多也不超过12K，另外其温度变化趋势无论是日照期还是阴影期，与实验曲线变化趋势基本接近，这验证了本发明所建立的重力条件下热管吸热器蓄热单元计算方法、模型及程序的合理性与准确性。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (5)

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

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

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

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

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

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

a)连续性方程

$\frac{\&PartialD;u}{\&PartialD;z}+\frac{1}{r}\frac{\&PartialD;\left(\mathrm{rv}\right)}{\&PartialD;r}=0$

式中：u，v分别为z方向、r方向的速度分量；

b)动量方程

$\frac{\&PartialD;u}{\&PartialD;t}+\frac{{\&PartialD;u}^2}{\&PartialD;z}+\frac{1}{r}\frac{\&PartialD;\left(\mathrm{ruv}\right)}{\&PartialD;r}=-\frac{1}{\mathrm{\&rho;}}\frac{\&PartialD;P}{\&PartialD;z}+\mathrm{g\&beta;}(T-T_m)+v{\&dtri;}^{2}u$

$\frac{\&PartialD;v}{\&PartialD;t}+\frac{\&PartialD;\left(uv\right)}{\&PartialD;z}+\frac{1}{r}\frac{\&PartialD;\left({\mathrm{rv}}^2\right)}{\&PartialD;r}=-\frac{1}{\mathrm{\&rho;}}\frac{\&PartialD;P}{\&PartialD;r}+v({\&dtri;}^{2}v-\left(\frac{v}{r^2}\right))$

式中：P为压力，g为重力加速度，β为相变材料(PCM)液态体积膨胀系数，v为运动粘度。

c)能量方程

$\frac{\&PartialD;\left(\mathrm{\&rho;e}\right)}{\&PartialD;t}+\frac{\&PartialD;\left(\mathrm{\&rho;ue}\right)}{\&PartialD;z}+\frac{1}{r}\frac{\&PartialD;\left(\mathrm{r\&rho;ve}\right)}{\&PartialD;r}={k\&dtri;}^{2}T$

式中：e为比焓，以比焓的形式把相变的影响考虑进去，适用于整个计算区域，其中焓与温度的关系式可表示成如下的形式：

$T=\left\{\begin{array}{cc}{T}_m+e/c,& e\&le;0\\ T_m,& 0<e<{\mathrm{\&Delta;H}}_m\\ T_m+(e-{\mathrm{\&Delta;H}}_m)/c,& e\&GreaterEqual;{\mathrm{\&Delta;H}}_m\end{array}\right.$

式中：e为比焓，ρ为密度，t为时间，k为导热系数，c为比热容，T为温度，T_m为相变温度，ΔH_m为物质单位质量的相变潜热；

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

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

q_sun＝15424-24.88(T_wall-1020)

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

q_shadow＝-12314-24.88(T_wall-1020)

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

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

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

(9)对步骤(5)中的偏微分控制方程进行离散化得到代数方程组，并利用步骤(6)的边界条件和步骤(7)的初始条件进行封闭和求解；

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

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

2.根据权利要求1所述的重力条件下高温蓄热容器内固/液相变的数值模拟方法，其特征在于，在上述步骤(9)中，计算每个单元与相邻单元的换热量时，需要对该单元及相邻单元内是否存在空穴作出判断，以正确计算该单元内相变材料(PCM)的质量以及与相邻单元的换热量，当有空穴存在时，PCM单元内PCM的质量为

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

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

式中：i、j、k代表第(i，j，k)个单元，II、JJ、KK代表PCM容器各坐标方向划分的单元数，m代表相变材料(PCM)的质量，ρ代表相变材料(PCM)的密度，V为单元的体积，f_v为空穴体积分数，f_v＝0，无空穴；0＜f_v＜1，单元内存在部分空穴；f_v＝1，单元内全部为空穴。

3.根据权利要求2所述的重力条件下高温蓄热容器内固/液相变的数值模拟方法，其特征在于，还需要考虑相变材料(PCM)体积变化引起的某些单元内空穴体积的改变，由前后两个时刻固态相变材料(PCM)的质量之差Δm求出每一个时间步长内相变材料(PCM)的体积变化量ΔV，根据体积变化量ΔV对相变材料(PCM)容器的空穴体积进行调整：如果ΔV＞0，说明相变材料(PCM)的体积增加，就减小空穴所占的体积；如果ΔV＜0，说明相变材料(PCM)的体积缩小，相应地增大空穴体积。

4.根据权利要求1至3任一项所述的重力条件下高温蓄热容器内固/液相变的数值模拟方法，其特征在于，在步骤(10)中，在求解离散化后的能量方程时采用显式格式，在求解离散化后的动量方程时采用交替方向块迭代法(ADI)与块追赶法(TDMA)相结合的计算方法。

5.根据权利要求4所述的重力条件下高温蓄热容器内固/液相变的数值模拟方法，其特征在于，采用求解压力耦合方程的半隐方法(SIMPLE算法)求解流场，为防止迭代过程发散，对速度u、v以及压力修正进行亚松弛。

Images