方腔内电场强化固液相变传热

全文HTML

--> --> -->

2.问题描述

2.1.物理模型及控制方程

-->

2.1.物理模型及控制方程

如图1所示, 固体PCM材料填充整个正方形腔体, 腔体的左右壁面分别放置高温($ \theta_{\rm h} $

)、高电势($ \phi_1 $

)和低温($ \theta_{\rm c} $

)、低电势($ \phi_0 $

) 的平板电极, 上下壁面绝缘绝热, 自由电荷从左电极板均匀注入. 初始时刻PCM 固定为熔化温度($ \theta_{\rm m} = \theta_{\rm c} $

). 此外, 基于电荷是否可以在固体中进行输运, 介电有机PCM一般可以从理论上分为两类^[36], 即欧姆固体和非欧姆固体. 对于欧姆固体来说, 电流在固体输运过程中遵循欧姆定律, 而对于非欧姆固体来说, 电荷在固体中的输运机制主要通过迁移和扩散来实现. 本文拟采用文献[32]中的假设, 即这里所涉及的介电PCM是一类所谓的非欧姆固体. 此外, 方腔内的流体是不可压缩的牛顿流体, 且流体密度随温度的变化遵循标准的Boussinesq假设. 另一方面, 为了简化问题, 忽略介电PCM在熔化过程中固相和液相物性参数的差异及变化. 如此, 方腔内电场强化固液相变问题可以通过连续方程、动量方程、高斯定律、电场定义方程、电荷守恒方程以及能量方程共同进行描述^[36,37]:

图 1 物理模型示意图
Figure1. Schematic diagram of the physical model.

$ \nabla \cdot {{u}} = 0, $

$\begin{split}\frac{{\partial {{u}}}}{{\partial t}} + \nabla \cdot {{uu}} =\;& - \nabla {p} + \nabla \cdot \left[ {{\upsilon}\left( {\nabla {{u}} + {{\left( {\nabla {{u}}}\right)}^{{\rm T}}}} \right)} \right] \\&+q{{E}} - {{g}}\left[ {1 - {\beta} \left( {{\theta} - {{\theta _{\rm r}}}} \right)} \right], \\[-12pt]\end{split}$

$ \nabla \cdot \left( {{\varepsilon} {\nabla} {\phi} } \right) = - {q}, $

$ {{E}} = - {\nabla {\phi}}, $

$ \frac {\partial {q}}{\partial t}+\nabla \cdot \left[{q \left( {K{{E}} } + {{u}} \right)}\right] = \nabla \cdot \left(D\nabla q \right) , $

$ \frac{\partial H}{\partial t}+\nabla \cdot (c_{\rm p}\theta {{u}}) = \frac{1}{\rho _0} \nabla \cdot (\lambda \nabla \theta), $

式中, 矢量u, g, E分别为流体速度、重力加速度和电场强度; 标量q, $ \theta $

, $ \phi $

, H则分别表示电荷密度、温度、电势和总焓; p, $ \upsilon $

, $ \beta $

, $ \varepsilon $

, K, D, $ c_{\rm p} $

, $ \rho_0 $

, $ \lambda $

分别代表修正压力、液体运动学黏性、液体热膨胀系数、介电常数、离子迁移率、电荷扩散系数、定压比热容、密度以及热导率. 需要指出的是, 总焓H为潜焓与显焓的和, 即

$ H = c_{\rm p}\theta +f_{\rm l} L, $

其中, 液率$ f_{\rm l} $

为PCM熔化过程中, 已经熔化部分占PCM总体积的比例, L为PCM 的潜热.
根据文献[32, 38]该系统实质上可由以下无量纲参数控制:

$ \begin{split}&Ra = \frac{g \beta \left({\theta_{\rm h}}-{\theta_{\rm c}} \right) {\delta}^{3}}{{\upsilon}{\chi}},\;\; T = \frac{\varepsilon\left(\phi_{1}-\phi_{0}\right)}{\mu K}, \\&C = \frac{q_{1} {\delta}^{2}}{\varepsilon\left(\phi_{1}-\phi_{0}\right)},\;\;M = \frac{1}{K}\left(\frac{\varepsilon}{\rho}\right)^{1 / 2},\\&\alpha = \frac{D}{K\left(\phi_{1}-\phi_{0}\right)},\\&Pr = {\upsilon\left(\rho c_{\rm p}\right)}/{\lambda},\;\; Ste = {c_{\rm p}\Delta \theta}/{L}, \end{split}$

这里, $ Ra $

和T分别为热瑞利数和电瑞利数, 用以表示浮升力和库仑力与黏性力的比值; C 是电荷注入强度, C的值越大, 注入到方腔内的电荷就越多; M表示无量纲离子迁移率, 其值由离子和流体的物性参数决定; $ \alpha $

为无量纲电荷扩散系数; $ Pr $

表征流体流动中动量交换与热交换相对比重; $ Ste $

反映PCM显热和潜热的比值. 另外, (8)式中涉及的$ \delta $

, $\,\chi$

以及$ \mu $

分别为参考长度、热扩散系数、动力学黏性.
此外, 为了在数据表征方面更好地解释电场强化固液相变过程的内在机理, 引入$ V_{\rm max} $

, $ Fo $

和$ \varPhi $

这三个无量纲参数, 它们的具体定义如下:

$ \begin{split} &V_{\rm max} = {\rm max} \big(\sqrt{u^2+v^2}\,\big), \quad \\&Fo = {\lambda t}{(\rho c_{\rm p})/{\delta}^2}, \quad \varPhi = [{Fo_{\rm b}-Fo_{\rm c}}]/{Fo_{\rm b}}, \end{split}$

其中, $ V_{\rm max} $

是方腔内的流体最大速度, 主要通过水平方向速度u和竖直方向速度v计算而来; $ Fo $

表示介电PCM的无量纲熔化时间, $ \varPhi $

表示电场强化固液相变的实际效果, 其中$ Fo_{\rm b} $

为只有浮升力存在时PCM熔化结束所需时间, $ Fo_{\rm c} $

为施加电场后熔化结束所需时间.
系统的边界条件主要由下面几个部分共同构成, 对于流场而言, 流体运动只发生在液体区域, 固相、固液界面以及腔体四周壁面均为无滑移无渗透边界, 即当$ x \!=\! 0 $

或$ \delta $

, $ y \!=\! 0 $

或$ \delta $

, PCM为固体时:

$ {{u}} = 0. $

对于温度场而言, 主要包括温度边界条件和焓边界条件两个部分, 其数学描述如下^[32]:

$ \begin{split} & x = 0: \quad \theta = \theta _{\rm h},\ H = c_{\rm p}\theta _{\rm h} + L,\\ &x = \delta: \quad \theta = \theta _{\rm c}, \ H = c_{\rm p}\theta _{\rm c},\\ &y = 0,\delta: \quad {\partial \theta}/{\partial{y}} = 0,\ {\partial H}/{\partial{y}} = 0. \end{split}$

此外, 电势边界条件采用与温度类似的边界条件, 即:

$ \begin{split}&x = 0 :\quad \phi = \phi _1, \\ &x = \delta :\quad \phi = \phi _0,\\ &y = 0,\delta :\quad {\partial \phi}/{\partial{y}} = 0. \end{split}$

最后, 电荷密度在左壁面为迪利克雷边界条件, 其他三个壁面为纽曼边界条件:

$ \begin{split}& x = 0 :\quad q = q _1, \\&x = \delta :\quad {\partial q}/{\partial{x}} = 0, \\&y = 0,\delta :\quad{\partial q}/{\partial{y}} = 0. \end{split}$

3.数值方法

在方法选择上采用所谓的介观格子Boltzmann方法(lattice Boltzmann method, LBM), 不同于传统的计算流体力学方法, LBM在处理多物理场耦合以及复杂边界方面具有显著的优势^[39], 并且该方法具有的天然并行性也为开展大规模数值模拟提供了可能^[40,41]. 近些年, LB方法已被广泛应用于诸如多相流^[42–44]、多孔介质流^[45,46] 以及电动流体动力学^{[36,38,47–49]}的相关研究中, 并已成为当下较为流行的计算流体力学方法. 在本文中, 为实现方腔内的电场强化固液相变传热模拟, 采用四个演化方程以求解Navier-Stokes方程、电势方程、电荷守恒方程以及能量方程.
2

3.1.流场格子Boltzmann模型

-->

3.1.流场格子Boltzmann模型

为模拟不可压流体流动, 流场采用Guo等^[50]提出的不可压LB模型进行求解, 该模型的演化方程可以表示如下:

$ \begin{split} &{f_i}\left( {{{x}} + {{{c}}_i}\Delta t,t + \Delta t} \right) - {f_i}\left( {{{x}},t} \right) \\=\;& - \frac{1}{{{\tau _{\rm l}}}}\left[ {{f_i}\left( {{{x}},t} \right) - f_i^{(\rm eq)}\left( {{{x}},t} \right)} \right] + \Delta t{F_i},\end{split}$

其中$ {\tau _{\rm l}} $

为无量纲松弛时间, 它与$ \upsilon $

之间的关系为$ \upsilon = \rho_0 c_{\rm s}^2 (\tau_ {\rm l}-0.5)\Delta t $

, $ {f_i}\left( {{{x}}, t} \right) $

为t时刻x处方向为i的粒子分布函数, 平衡态分布函数$ f_i^{\rm eq} $

由下式给出:

$ f_i^{(\rm eq)}({{x}},t) = {\eta _i}p + {\omega _i}\left[ {\frac{{{{{c}}_i} \cdot {{u}}}}{{c_{\rm s}^2}} + \frac{{{{uu}}:({{{c}}_i}{{{c}}_i} - c_{\rm s}^2{{I}})}}{{2c_{\rm s}^4}}} \right], $

其中$ {\eta _i} $

为模型参数, 它满足$ {\eta _0} = (\omega _0-1)/{c_{\rm s}^2}+\rho_0 $

, $ {\eta _i} = \omega _i/{c_{\rm s}^2}\ (i\neq 0) $

, $ \rho_0 $

为流体密度, $ {\omega _i} $

与${{c}}_i$

分别为权系数以及离散速度, 值由所选取的DnQq模型(n维q速模型)确定, 针对当前的二维问题, 速度场使用经典的D2Q9模型, 对应的$ {\omega _i} $

和${{c}}_i$

定义如下:

$ {{\omega }_{i}}=\left\{ \begin{aligned} & 4/9, \;\;\;i=0, \\ & 1/9,\;\;\; i=1-4, \\ & 1/36, \;\; i=5-8, \end{aligned} \right.$

$ {{{c}}_i} = \left\{ {\begin{aligned}&{(0,0),}\;\;\;\qquad\qquad\qquad\qquad\qquad\qquad{i = 0,}\\&{c\left( {\cos \left[ {(i - 1)\frac{\pi }{2}} \right],\sin \left[ {(i - 1)\frac{\pi }{2}} \right]} \right),}\;\;{i = 1—4,}\\&\sqrt 2 c\left( \cos \left[ {(2i - 1)\frac{\pi }{4}} \right],\right.\left.\sin \left[ {(2i - 1)\frac{\pi }{4}} \right] \right),\\& \quad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\;\;{i = 5—8.}\\\end{aligned}} \right.$

此外, $ c_{\rm s} = c/\sqrt 3 $

为格子声速, 其中$ c $

为离散速度, 可以通过网格步长$ \Delta x $

和时间步长$ \Delta t $

计算得到$ c = \Delta x/\Delta t $

. 另一方面, 源项$ F_i $

采用Guo等^[51]提出的外力格式:

$ \begin{split}\;& {F_i}({{x}},t) ={\omega _i}\left( {1 - \frac{1}{{2{\tau _{\rm{l}}}}}} \right) \bigg[ \frac{{{{{c}}_i} \cdot {{F}}}}{{c_{\rm s}^2}} \\& \qquad \qquad + \frac{ {\left({{Fu}} + {{uF}}\right):\left( {{{{c}}_i}{{{c}}_i} - c_{\rm s}^2{{I}}} \right)}}{{2c_{\rm s}^4}} \bigg], \end{split}$

其中$ {{F}} = q{{E}}+{{g}}\beta \left(\theta-\theta_{\rm c} \right) $

. 最后, 宏观速度$ {{u}} $

以及压力项p由以下公式计算得到:

$ {{u}} = \sum\limits_i {{{{c}}_i}{f_i}} + \frac{{\Delta t}}{2}{{F}}, $

$ p = \frac{{c_{\rm s}^2}}{{1 - {\omega _0}}}\left( {\sum\limits_{i \ne 0} {{f_i}} - {\omega _0}\frac{{{{\left| {{u}} \right|}^2}}}{{2c_{\rm s}^2}}} \right). $

3.2.电荷密度的格子Boltzmann模型

-->

3.2.电荷密度的格子Boltzmann模型

电荷密度的演化方程主要是借鉴Luo等^[48]于2016年提出的模型, 其演化方程由下式给出:

$\begin{split} &{h_i}\left( {{{x}} + {{{c}}_i}\Delta t,t + \Delta t} \right) - {h_i}\left( {{{x}},t} \right) \\=\;& - \frac{1}{{{\tau _{\rm q}}}}\left[ {{h_i}\left( {{{x}},t} \right) - h_i^{(\rm eq)}\left( {{{x}},t} \right)} \right], \end{split}$

其中电荷密度的平衡态分布函数为

$\begin{split} \;& h_i^{(\rm eq)}({{x}},t) = {q{\omega _i}} \left\{ 1+\frac{ {{{{c}}_i} \cdot \left( {K{E}} +{{u}} \right) } } {{c_{\rm s}^2}}\right.\\& \qquad \qquad \left.+ \frac{ \left[ {{c}}_i \cdot \left(K{{E}}+{{u}} \right) \right]^{2} - c_{\rm s}^2 \left| K{{E}}+{{u}} \right|^{2} } {{2c_{\rm s}^4}} \right\}. \end{split}$

松弛时间$ {\tau _{\rm q}} $

与电荷扩散系数D的关系如下:

$ {D} = {c_{\rm s}^{2}}\left({\tau _{\rm q}} -0.5 \right) \Delta t. $

最后, 电荷密度可通过下面的表达式计算得到:

$ q = \sum\limits_{i} h_{i}. $

电荷密度计算所使用的网格模型与速度场相同, 即D2Q9.
2

3.3.电势的格子Boltzmann模型

-->

3.3.电势的格子Boltzmann模型

从电势控制方程的形式上看, 其本质是一类典型的Poisson方程. 为此, 这里采用Chai和Shi^[52]提出的求解Poisson方程的LB 模型, 其演化方程可以表述如下:

$\begin{split}&{g_i}\left( {{{x}} + {{{c}}_i}\Delta t,t^\prime + \Delta t} \right) - {g_i}\left( {{{x}},t^\prime} \right) \\=\;& - \frac{1}{{{\tau _{\phi}}}}\left[ {{g_i}\left( {{{x}},t^\prime} \right) - g_i^{(\rm eq)}\left( {{{x}},t^\prime} \right)} \right] + \Delta t {\hat{\omega}}_{i} \xi R, \end{split}$

这里$ t^\prime $

为迭代伪时间, $ \xi $

为自定义的扩散系数. 为保持数值稳定, 在接下来的数值模拟过程中, 取$ \xi $

等于0.5^[48]. 此外, 源项R为$ R = q/\varepsilon $

. 平衡态分布函数定义如下:

$ g_{i}^{(\rm eq)}=\left\{ \begin{aligned} &\left( {{{\hat{\omega }}}_{0}}-1 \right)\phi , \;\; i=0, \\ & {{{\hat{\omega }}}_{i}}\phi ,\;\;\;\;\;\qquad i=1-4, \end{aligned} \right. $

由于上述平衡态是线性平衡态形式, 为加快运算效率, 这里的离散速度采用D2Q5模型^[52], 该模型的权系数取值如下:

$ {{\hat{\omega }}_{i}}=\left\{ \begin{aligned} & 0, \;\;\quad i=0, \\ & 1/4, \;\; i=1-4. \end{aligned} \right. $

另一方面, 格子声速$ {\hat {c}}_{\rm s}^2 = c^2/2 $

, 松弛时间通过关系式$ \xi = {\hat {c}}_{\rm s}^2\left(\tau_{\phi}-0.5 \right)\Delta t $

计算得到. 而电势则通过下式计算得到:

$ {\phi} = \frac {1}{1-{\hat{\omega} _0}} { \sum\limits_{i} g_{i}}. $

最后, 电场E可以借助电势的一阶分布函数计算得到^[38], 即

$ {{E}} = \frac{1}{{\tau _\phi}{\hat {c}}_{\rm s}^2}\sum\limits_{i}{{c}}_{i}{g_i}. $

3.4.基于焓的格子Boltzmann模型

-->

3.4.基于焓的格子Boltzmann模型

能量方程采用Lu等^[53]提出的双松弛时间LB模型, 不同于现有的部分固液相变LB模型^[37,53], 该模型能够有效抑制相变模拟过程中固液界面处的非物理扩散现象, 其演化方程为

$ \begin{split} &l_{i}\left({{x}}+{{c}}_{i} \Delta t, t+\Delta t\right)= \\\;& l_{i}^{(\rm eq)}({{x}}, t)+\left(1-\frac{1}{2 \tau_{\rm s}}-\frac{1}{2 \tau_{\rm a}}\right) \left[l_{i}({{x}}, t)-l_{i}^{(\rm eq)}({{x}}, t)\right] \\& -\left(\frac{1}{2 \tau_{\rm s}}-\frac{1}{2 \tau_{\rm a}}\right) \left[{l}_{\bar{i}}({{x}}, t)-{l}_{\bar{i}}^{(\rm eq)}({{x}}, t)\right],\\[-10pt] \end{split}$

其中$ \tau_{\rm s} $

和$ \tau_{\rm a} $

分别为对称松弛时间和反对称松弛时间, $ {l}_{\bar{i}}({{x}}, t) $

和$ {l}_{\bar{i}}^{(\rm eq)}({{x}}, t) $

分别是$ l_{i}({{x}}, t) $

和$ l_{i}^{(\rm eq)}({{x}}, t) $

相反方向的分布函数和平衡态分布函数. $ l_{i}^{(\rm eq)}({{x}}, t) $

由下式给出:

$ l_i^{\rm eq} = \left\{ {\begin{aligned}&{H - {c_{\rm p}}\theta + {\omega _i}{c_{\rm p}}\theta \left( {1 - \frac{{|{{u}}{|^2}}}{{2c_{\rm s}^2}}} \right),}&{i = 0,}\\&{{\omega _i}{c_{\rm p}}\theta \left[ {1 + \frac{{{{{c}}_i} \cdot {{u}}}}{{c_{\rm s}^2}} + \frac{{{{({{{c}}_i} \cdot {{u}})}^2}}}{{2c_{\rm s}^4}} - \frac{{|{{u}}{|^2}}}{{2c_{\rm s}^2}}} \right],}&{i \ne 0,}\end{aligned}} \right.$

$ \tau_{\rm a} $

与$ \lambda $

的关系为

$ \frac{\lambda}{\rho _0 c_{\rm p}} = c_{\rm s}^2(\tau_{\rm a}-0.5)\Delta t, $

而$ \tau_{\rm s} $

可通过以下等式得到:

$ {1}/{\tau_{\rm s}}+ {1}/{\tau_{\rm a}}= 2. $

焓的计算方式为

$ H = \sum\nolimits_{i} l_{i}. $

最后, $ f_{\rm l} $

以及$ \theta $

可由焓计算得到:

$ {f_{\rm l}} = \left\{ {\begin{aligned}&0,\qquad\qquad{H \leqslant {H_{\rm s}}},\\&{\frac{{H - {H_{\rm s}}}}{{{H_{\rm l}} - {H_{\rm s}}}}},\;\;{{H_{\rm s}} < H < {H_{\rm l}}},\\&1,\qquad\qquad{H \geqslant {H_{\rm l}}},\end{aligned}} \right. $

$ \theta \!=\! \left\{ {\begin{aligned}&{ {H}/{{{c_{\rm p}}}}},\qquad\qquad\qquad\qquad\quad {H \leqslant {H_{\rm s}}},\\&{\frac{{{H_{\rm l}} - H}}{{{H_{\rm l}} - {H_{\rm s}}}}{\theta _{\rm s}} \!+\! \frac{{H - {H_{\rm s}}}}{{{H_{\rm l}} - {H_{\rm s}}}}{\theta _{\rm l}}},\;\;\;{{H_{\rm s}} \!<\! H \!<\! {H_{\rm l}}},\\&{{\theta _{\rm l}} + \frac{{H - {H_{\rm l}}}}{{{c_{\rm p}}}}},\qquad\qquad\quad{H \geqslant {H_{\rm l}}},\end{aligned}} \right. $

其中$ H_{\rm s} = c_{\rm p}\theta _{\rm s} $

和$ H_{\rm l} = c_{\rm p}\theta_{\rm s}+L $

分别代表固相线和液相线总焓, $ \theta _{\rm s} $

和$ \theta _{\rm l} $

分别为固相线温度和液相线温度, 本文中考虑单区域熔化问题, 即$ \theta _{\rm s} = \theta _{\rm l} = \theta _{\rm m} $

. 网格模型与速度场及电荷密度相同.

5.结果与讨论

5.1.纯电场力作用下的固液相变问题探究

-->

5.1.纯电场力作用下的固液相变问题探究

为了更好地理解电场的存在对方腔中PCM熔化的影响, 首先模拟了只有电场力($ T = 600 $

)存在时介电PCM的熔化过程. 此外, 作为对比, 也同样模拟了只有浮升力($ Ra = 10000 $

)存在时的熔化过程. 图4 为PCM熔化过程中液体体积分数($ f_{\rm l} $

)以及液相流体流动最大速度($ V_{\rm max} $

)随无量纲时间($ Fo $

)的变化曲线. 如图4(a)所示, 虽然在熔化初期($Fo\leqslant 1.1$

)和熔化中期($1.1 < Fo \leqslant 2.8$

), 电场驱动下PCM的熔化速率没有浮升力驱动下PCM 的熔化速率快, 但是在熔化后期($Fo > 2.8$

), 浮升力驱动下PCM的熔化速率逐渐减慢, 而电场力驱动下的PCM始终以较快的速率进行熔化, 熔化结束所消耗的总时间为$ 4.62 $

, 与浮升力驱动下PCM熔化结束所需时间$ 6.15 $

相比, 减少了大约$ 24\% $

, 可以看到, 电场力的存在对PCM熔化特别是在熔化后期有一定的加速效果.

图 4 PCM熔化过程中(a)液体体积分数$f_{\rm l}$

和(b)液相最大速度$V_{\rm max}$

随时间的变化趋势. 图中点A, B, C和D分别对应时刻$Fo=0.63$

, $1.90$

, $2.90$

和$4.16$

Figure4. Time evolutions of (a) the liquid fraction and (b) the maximum fluid velocity. Points A, B, C and D correspond to the dimensionless time $Fo=0.36$

, $1.90$

, $2.90$

and $4.16$

, respectively.

图5上栏给出电场力驱动下PCM熔化到$ Fo = 0.63 $

, $ 1.90 $

, $ 2.90 $

, $ 4.16 $

四个时刻的电荷密度和流场分布图, 同样时刻浮升力驱动下的温度场以及流场的分布作为对比如下栏所示. 其中, 有流场的一侧为已熔化PCM, 右侧为未熔化部分, 二者交界处的曲线段为固液界面. 从图5可以看出, 熔化刚开始时, 电场力驱动下PCM液相区域对流尚未形成, 因此图4(b)中对应的流场最大速度基本为0. 此时导热熔化占优, 而浮升力驱动时, 虽然在$ Fo = 0.63 $

时, PCM的熔化主要也由导热传热贡献, 熔化速率同只有电场力存在时保持一致, 但是前者液相区流动涡结构已经形成, 并且不断发展. 从图4(b)可以看到, 流场最大速度不断增大, 对流传热逐渐增强, 结果是在浮升力的作用下, 高温壁面附近的热流体沿着高温壁面向上流动, 在上壁面的限制下, 液体携带大量的热量流向固液界面, 加速上半部分PCM的熔化. $ Fo = 1.90 $

时, 如图5(b) 所示, PCM 上半区已熔化部分比下半区多. 而此时对于只有电场力存在的情况, 随着液相区域不断扩大, 电场力逐渐克服流体黏性, 在液相区域形成沿着中心水平线对称的一对反螺旋一级漩涡, 同时由于受到上下壁面的限制, 两个角涡在固液界面的上下端附近形成, 流动强度迅速增大(如图4(b)所示). $ Fo = 2.90 $

时, 可以很明显地看到, PCM 的上半部分与下半部分熔化较快, 而中间部分熔化较慢, 这是由于电场作用下流场形成的反螺旋涡结构使得热壁面附近的热流分为两股, 分别沿着上下壁面流向固液界面的上下两端, 加速了这两部分PCM的熔化. 值得注意的是, 在液体区域中心出现一个“电荷空白区域”, 这是对称电极下电对流的主要特征, 可以简单解释为撞击固液界面后冷却的流体沿着固液界面在中心汇聚后, 又沿水平中心线流向热壁面, 此时在中心区域电荷的迁移速度小于流体流动的速度, 由热电极板注入的电荷无法沿中心水平线到达冷壁面, 从而形成了电荷空白区域. 此时对于浮升力驱动的情况, PCM 的上半部分已经熔化完毕, 流场最大速度不再增加, 这就意味着对流传热对熔化的贡献逐渐减少, 热传导成为热量传递的主要方式, 但是PCM 的热导率较低, 同时随着液相体积的扩大, 热阻也不断增大, 最终导致熔化速率不断减慢, 由图4(a) 可以看到, 在最后阶段, PCM熔化速率趋于0. 而对于电场力驱动的情况, PCM的上下部分熔化完毕以后, 熔化速率虽然有一定程度的放缓, 但是由于电对流强烈的非线性不稳定性(图4(b)中最大速度呈现不规则波动), 液相区域漩涡的数量以及形状不断地发生变化, 使得冷热流体得到充分混合, 提高了热壁面的换热效率, 因此PCM未熔化部分也可以及时熔化完毕.

图 5 电场力驱动下PCM熔化到不同时刻电荷密度、流场以及相界面分布(上)和浮升力驱动下PCM熔化到不同时刻温度、流场以及相界面分布(下)　(a)$Fo=0.63$

; (b)$Fo=1.90$

; (c)$Fo=2.90$

; (d)$Fo=4.16$

Figure5. Distributions of charge density (temperature), fluid field and liquid-solid interface of PCM at four representative time of (a)$Fo=0.63$

, (b)$Fo=1.90$

, (c)$Fo=2.90$

and (d)$Fo=4.16$

in presence of electric field force (top) and buoyancy force (bottom)

2

5.2.电场力与浮升力共同作用下的固液相变问题探究

-->

5.2.电场力与浮升力共同作用下的固液相变问题探究

通过5.1节的分析可以知道, 电场的存在对PCM熔化有一定的加速效果, 值得注意的是, 只有电场力驱动下PCM的中间部分最后熔化完, 而只有浮升力驱动时PCM的下半部分最后熔化完(如图5(d)所示). 本节探究当两种驱动力共同存在时, 电瑞利数T对方腔内PCM熔化速率的影响, 模拟中$ Ra $

, M, $ Pr $

, $ Ste $

分别固定为$ 10000 $

, $ 10 $

, $ 1.0 $

, $ 0.1 $

. 不同电瑞利数T下, PCM的液体体积分数$ f_{\rm l} $

随时间$ Fo $

的变化如图6所示. 可以看出, T较小时($T\leqslant 800$

), 外加电场的存在对PCM的熔化速率影响不大, 特别地, 当T取$ 400 $

时, PCM的熔化时间反而延长; 当T 较大时($T\geqslant 800$

), 随着T 的增大, PCM 熔化速率不断增大, 总的熔化时间不断减小.

图 6 不同电瑞利数T下液体体积分数$f_{\rm l}$

随时间的变化曲线
Figure6. Time evolutions of the liquid fraction $f_{\rm l}$

for different electric Rayleigh numbers T.

图7为电瑞利数T分别取$ 0 $

, $ 400 $

, $ 800 $

, $ 1600 $

, $ 2400 $

, PCM熔化到同一时刻($ Fo = 1.9 $

)时, 电荷密度、温度以及流场的分布图. 从图中可以看到, 当$ T = 400 $

时, 库仑力较弱, 浮升力起着主导性作用, 液体区域只有一个沿着顺时针旋转的一级漩涡, 与只存在浮升力时相似. 电荷从热电极板注入后, 大量的电荷受到热流的影响, 传输至上壁面附近, 相应地, 电场力的存在增大了此处的对流速度, 使得更多的热量被传输至固液界面, 加速PCM上半部分的熔化. 对于PCM 下半部分, 根据5.1节的讨论, 电场力的存在会使一股热流沿着下壁面朝固液界面流去, 而浮升力的存在又使得冷流体下沉后沿着下壁面流向热电极板, 两股流体流动方向相反, 但是由于浮升力强于库仑力, 所以流动方向与只存在浮升力时保持一致, 此时电场力的存在对下半区域的对流起阻碍的作用. PCM熔化到$ Fo = 1.9 $

时刻, T取$ 0 $

和$ 400 $

时, 直线$ x = 1/(4\delta) $

上速度水平分量u的分布如图8(a)所示, 而T 取$ 0 $

, $ 400 $

, $ 800 $

以及$ 1600 $

时PCM固液界面的位置如图8(b)所示. 从图8(a)可以看到, 施加电场以后, 靠近上壁面流体的速度水平分量较大, 这是电场力驱动下引起电对流的结果, 而在靠近下壁面的区域, 流体因受到电场力的阻碍作用流速较小, 该结果与以上分析一致. 所以相应地, 从图8(b)可以看到, 电场的存在加速了PCM上半区域的熔化速率, 而对下半区而言, 熔化速率有所减慢, 但是因为熔化由导热主导, 减慢并不明显. 这也很好地解释了为什么T较小时施加电场反而使得后期的熔化速率减慢. 当T增大到$ 800 $

时, 如图7(c)所示, 在PCM液相下半区域形成一个小漩涡, 根据5.1节的分析, 这是电场力影响的结果, 而小漩涡的形成以及不断发展, 导致下半区冷热流体充分混合, 从图8(b)可以看到, 此时PCM上半区和下半区的熔化速率都得到提高, 固液界面移动较快, 熔化完毕所需时间有所减少(见图6). 电场力继续增大, 从图7(d)可以明显地看到, 当T 达到$ 1600 $

时, PCM液相区域有多个涡形成, 需要指出的是, 下半区所形成的漩涡大小仅次于中心一级漩涡, 对流传热使得下半区固液界面移动速率得到提高(见图8(b)). 最后, PCM在$ T = 2400 $

时, 除了中间少部分区域以外, 其余大部分已经熔化完毕, 流动结构变得极其混乱, 此时库仑力起绝对主导作用, 热电极板的热量在强烈电对流的作用下被快速输送至固液界面, 加速了PCM熔化过程, 总的来看, 熔化结束花费的时间$ Fo = 2.1 $

, 与未加电场($ Fo = 6.2 $

)相比, 大约缩短66.1%, 电场的存在起到了加速熔化的作用.

图 7 不同电瑞利数下$Fo=1.9$

时电荷密度(上)、温度(中)以及流场(下)分布:　(a) $T=0$

; (b) $T=400$

; (c) $T=800$

; (d) $T=1600$

; (e) $T=2400$

Figure7. Distribution of the charge density (top), temperature (middle) and streamlines (bottom) at $Fo=1.9$

under different electrical Rayleigh numbers: (a) $T=0$

; (b) $T=400$

; (c) $T=800$

; (d) $T=1600$

; (e) $T=2400$

图 8 $Fo=1.9$

时　(a)直线$x=1/(4\delta)$

上速度水平分量u的分布; (b)固液界面位置
Figure8. (a) Distributions of the horizontal velocity along line of $x=1/(4\delta)$

and (b) the liquid-solid interface at $Fo=1.9$

除上述电瑞利数影响外, 接下来研究斯蒂芬数$0.05 \leqslant Ste \leqslant 1.0$

、离子迁移率$3\leqslant M\leqslant $

$ 40$

以及普朗特数$1\leqslant Pr \leqslant 20$

对电场强化固液相变的影响, 模拟过程中热瑞利数$ Ra $

固定为$ 10000 $

.
3

5.2.1.斯蒂芬数

$ Ste $

对电场强化相变传热过程的影响

-->

5.2.1.斯蒂芬数$ Ste $对电场强化相变传热过程的影响

图9(a)为不同电瑞利数下斯蒂芬数$ Ste $

对电场加速熔化效果的影响, 图9(b)为$ T = 1000 $

时, $ Ste = 0.05 $

和$ Ste = 1.0 $

的PCM熔化过程中, 液相最大速度随液率$ f_{\rm l} $

变化趋势. 可以看到, 对于给定的T, 随着$ Ste $

的增大, 加速百分比$ \varPhi $

呈现下降趋势, 这表示, PCM的潜热值越大, 电场对熔化的加速效果越好. 此外, 还可以得到当电瑞利数T较大时, $ Ste $

的影响逐渐减弱的结论. 实际上, 潜热值越小表示PCM发生固液相变所需热量越少, 熔化速率越快, 熔化完毕所需时间越短, 这样电场驱动引起的电对流发生作用时间便减少.

图 9 (a)斯蒂芬数$Ste$

对电场加速熔化效果的影响; (b)液相最大速度$V_{\rm max}$

随液体体积分数$f_{\rm l}$

的变化趋势
Figure9. (a) The effect of Stefan number $Ste$

on the melting time saving; (b) evolutions of the maximum fluid velocity $V_{\rm max}$

versus the liquid fraction $f_{\rm l}$

.

图10和图11分别为T取$ 1000 $

和$ 3200 $

时, 不同的斯蒂芬数$ Ste = 0.05 $

, $ Ste = 1.0 $

下PCM 熔化过程中液体体积分数$ f_{\rm l} $

分别为0.3, 0.5, 0.9 时电荷密度、温度以及流场分布图. 可以看到, 当$ T = $

$ 1000 $

时, $ Ste = 0.05 $

的PCM 在熔化到$ f_l = 0.5 $

时(图10(a2)), 液相中心已经形成明显的电荷空白区并且热流体在下半部分形成的反向漩涡作用下, 撞向固液界面, 加速下板附近PCM的熔化, 固液界面发生明显扭曲, 而$ Ste = 1.0 $

的PCM熔化到$ f_{\rm l} = 0.5 $

时(图10(b2)), 电对流并未充分发展. 此外, 从图9(b)可以明显地看到, $ Ste = 0.05 $

时, PCM熔化到$ f_{\rm l} = 0.45 $

时, 流场涡结构充分发展, 且由于强烈的不稳定性最大速度开始发生不规则波动, 这有利于冷热流体的充分混合, 而对于$ Ste = $

$ 1.0 $

的情况, 流场涡结构发展相对滞后, 在熔化到$ f_{\rm l} = 0.7 $

以后, 涡结构才基本发展充分, 电场强化相变传热作用时间相对缩减, 导致最终加速效果受到影响. 与$T\! =\! 1000$

时不同, 当$T\! = \! 3200$

时(图11), 较强的电场力能够使电对流在熔化早期便充分发展, 使得PCM在熔化初期的导热主导传热机制快速过渡到对流传热占优阶段, 如图11(a1)和图11(b1)所示, $ Ste = 0.05 $

以及$ Ste = 1.0 $

的PCM在熔化到$ f_{\rm l} = 0.3 $

时, 固液界面均已经发生了不同程度的扭曲, PCM在电场力的主导下以较高的速率进行熔化, $ Ste $

的影响不明显.

图 10 $T=1000$

时, 不同的斯蒂芬数下, PCM熔化到$f_{\rm l}=0.3$

, $f_{\rm l}=0.5$

, $f_{\rm l}=0.9$

时电荷密度、温度以及流场分布　(a)$Ste= $

$ 0.05$

; (b)$Ste=0.1$

Figure10. With $T=1000$

, distributions of charge density, temperature field and streamlines of PCM melting at $f_{\rm l}=0.3$

, $f_{\rm l}=0.5$

, $f_{\rm l}=0.9$

for different Stefan number: (a) $Ste=0.05$

; (b) $Ste=0.1$

图 11 $T=3200$

时, 不同的斯蒂芬数下, PCM熔化到$f_{\rm l}=0.3$

, $f_{\rm l}=0.5$

, $f_{\rm l}=0.9$

时电荷密度、温度以及流场分布　(a)$Ste= $

$ 0.05$

; (b)$Ste=0.1$

Figure11. With $T=3200$

, distributions of charge density, temperature field and streamlines of PCM melting at $f_{\rm l}=0.3$

, $f_{\rm l}=0.5$

, $f_{\rm l}=0.9$

for different Stefan number: (a) $Ste=0.05$

; (b) $Ste=0.1$

.

3

5.2.2.离子迁移率M与普朗特数

$ Pr $

对电场强化相变传热过程的影响

-->

5.2.2.离子迁移率M与普朗特数$ Pr $对电场强化相变传热过程的影响

电场对PCM熔化的加速效果与无量纲迁移率M以及普朗特数$ Pr $

的关系如图12 所示, T固定为$ 1000 $

, 虚线表示没有外加电场的情况. 可以看到, 对于给定的$ Pr $

, 随着M的增大, 加速百分比$ \varPhi $

呈下降趋势, 这表示电场力对M较小的PCM的相变传热强化效果较好. 特别地, 对于$ Pr = 1.0 $

的情况, 当$M\leqslant 10$

时, 随着M的增大, 电场对PCM熔化的加速效果急剧下降, 而当$M\geqslant 10$

以后, 电场对PCM熔化过程的影响逐渐减弱甚至消失($ M = 40 $

). 为了解释此现象, 以黏性速度$u_{\upsilon} = \upsilon / \delta$

作为特征速度将动量守恒方程写成无量纲形式:

图 12 无量纲离子迁移率M与普朗特数$Pr$

对电场强化相变传热效果的影响
Figure12. The effect of M and $Pr$

on the melting time saving.

$ \begin{split}\frac{{\partial {{u}}}}{{\partial t}} + \nabla \cdot {{uu}} =\;& - \nabla {p} + \nabla \cdot \left[ {{\upsilon}\left( {\nabla {{u}} + {{\left( {\nabla {{u}}}\right)}^{{\rm T}}}} \right)} \right] \\&+\left(\frac{T}{M}\right)^2qC{{E}} +\frac{Ra}{Pr} \theta {{e}} .\\[-15pt] \end{split}$

可以看到, 外力项包含库仑力驱动项和浮升力驱动项, 随着M的增大, 库仑驱动力项逐渐减小, 驱动力的减小使得液相流动速度减小, 对流传热受到影响, 熔化速率便逐渐放缓. 另外, T可以看成是$ 10^3 $

量级, 当$M\leqslant 10$

时, 库仑力在流场外力中占优, M的改变直接关系整个流场外力的大小, 在这种情况下, M的增大使得对流换热效率急剧下降. 而当$M\geqslant 10$

时, 虽然M的增大仍然会导致库仑力项逐渐减弱, 但是由于库仑力项的量级发生改变, 浮升力的影响逐渐显现, 总的流场外力减小并不明显. 当M足够大时, 电场力对流场外力的贡献可以忽略, 电场的存在对PCM的熔化也就没有实质性影响. 实际上, 无量纲离子迁移率表征静电方程组与Navier-Stokes方程的耦合紧密程度^[38], M 越大, 耦合越不紧密. $ Pr = 1.0 $