流动不稳定是制约自然循环系统安全稳定运行的重要因素^[1]，特别是开式自然循环系统，由于其系统压力较低，易于发生沸腾、流动闪蒸等剧烈的相变，同时缺乏能够容纳系统压力波动，抑制空泡份额变化及流量振荡的封闭气空间，因此开式自然循环系统更容易发生流动不稳定现象^[2]。

在开式自然循环流动不稳定中，闪蒸不稳定是重要的不稳定类型之一。闪蒸不稳定发生时，上升段内出现间歇性的闪蒸流动或者闪蒸起始点沿上升通道周期性地上下移动，从而导致自然循环流动出现振荡。目前，针对开式自然循环闪蒸不稳定的研究多为实验研究^[3-6]，理论分析研究也集中于采用数值解析方法进行求解^[7-8]，而采用近似分析方法对闪蒸不稳定进行研究的文献很少^[9]。本文提出了描述开式自然循环系统闪蒸流动不稳定的线性均相流模型，并基于系统控制原理中的状态空间分析方法对闪蒸不稳定边界进行分析。

1 数学模型建立 1.1 开式自然循环闪蒸流动物理过程及基本假设

图 1给出了开式自然循环闪蒸流动时，加热段和上升段内的温度压力的变化情况。加热段出口流体处于过冷状态，随着流体向上流动，其承受的重位压头和相应的饱和温度均不断降低，导致流体运动至上升段某一位置时达到当地饱和温度而快速汽化，形成汽液两相流动，该位置称为闪蒸起始点。

系统出现稳定闪蒸流动时，闪蒸起始点固定在上升段的某一位置上，而系统呈现闪蒸不稳定特性时，闪蒸起始点间歇性地消失或沿上升段周期性上下移动，当闪蒸起始点消失(或下移)时，上升段内空泡份额降低，自然循环驱动降低，使得自然循环流量降低，上升段内流体温度增加，进而导致闪蒸起始点出现(或上移)。随后上升段内空泡份额升高，循环流量上升，上升段内流体温度降低，又导致闪蒸起始点再次消失(或下移)。如此反复形成了自持性的闪蒸流动不稳定。

为了建立自然循环系统闪蒸不稳定的线性均相流模型，本文针对上述物理过程提出如下假设：

1) 流体为一维流动，流型的影响可忽略不计；

2) 上升段内流动过程(包括闪蒸过程)为热力学平衡过程，同时忽略加热段内过冷沸腾的影响；

3) 考虑到系统内温度和压力的变化范围较小，故忽略温度和压力变化对于液相密度ρ_l和气相密度ρ_g的影响。

1.2 基本控制方程及边界条件

针对上升段内的闪蒸流动过程，列出均相流的质量守恒方程和动量守恒方程:

$ {\partial \over {\partial t}}\left[ {{\rho _g}\alpha + {p_1}\left( {1 - \alpha } \right)} \right] + {{\partial G} \over {\partial z}} = 0$

(1)

$ \eqalign{ & - {{\partial P} \over {\partial z}} = {{\partial G} \over {\partial t}} + \left[ {\left( {1 - \alpha } \right){p_1} + \alpha {\rho _g}} \right]g + {{\xi {G^2}\varphi _f^2} \over {2D}} + \cr &{\partial \over {\partial z}}\left\{ {{G^2}\left[ {{{{{\left( {1 - x} \right)}^2}} \over {\left( {1 - \alpha } \right){p_1}}} + {{{x^2}} \over {{\rho _g}\alpha }}} \right]} \right. \cr} $

(2)

式中：t为时间，ρ为密度，α为空泡份额，G为质量流速，z为距加热段入口距离，P为压力，D为当量直径，φ_f为分液相折算系数，ξ为摩擦阻力系数，下标g表示气相，下标l表示液相。同时列出能量守恒方程时，忽略了外界力做功及动能变化引起的能量变化：

$ \eqalign{ &{q \over A} = {\partial \over {\partial t}}\left[ {{\rho _1}\left( {1 - \alpha } \right){h_1} + {\rho _g}\alpha {h_g}} \right] + \cr &{\partial \over {\partial z}}\left\{ {G\left[ {\left( {1 - x} \right){h_1} + x{h_g}} \right]} \right\} \cr} $

(3)

式中：q为线功率密度，A为流道横截面积，h_l和h_g分别为液相焓值和气相焓值，x为干度。

另外，根据漂移流模型得到上升段两相区内流体的质量含气率x和空泡份额α之间的关系，列出辅助方程：

$ \alpha = {x \over {{C_0}\left[ {{{x\left( {{\rho _1} - {\rho _g}} \right)} \over {{\rho _1}}} + {{{\rho _g}} \over {{\rho _1}}}} \right] + {{{\rho _g}{v_{gj}}} \over G}}}$

(4)

式中：C₀为分布参数；υ_gj为气相漂移速度，其具体数值的选取参照Zuber&Findlay公式^[10]。

最后列出了开式自然循环上升段内闪蒸流动的边界方程，即

加热段入口位置

$ {P_{c.in}} = {P_e} - {1 \over 2}\left( {1 + {\xi _{c,in}}} \right){G^2}$

(5)

$ {x_{c.in}} = \left( {{h_{c.in}} - {h_{fa}}} \right)/{h_{gl}}$

(6)

上升段入口位置：

$ {P_{c.ex}} = {P_b} - {1 \over 2}{\xi _b}{G^2}$

(7)

上升段出口位置：

$ {P_a} = {P_e} - {1 \over 2}\left( {{\xi _{r,ex}} - 1} \right)\varphi _f^2{G^2}$

(8)

式中：下标a为上升段出口位置，b为上升段入口位置，c表示加热段区域，e为外界环境压力，in和ex分别为入口和出口位置，f表示饱和状态的参数值，x_c,in为相对于上升段出口处饱和状态的加热段入口相对干度，h_gl为汽化潜热。

1.3 方程的无量纲化

对上述方程中参数进行无量纲化：

$ \eqalign{ &\hat G = {G \over {{G_{co}}}},\hat t = {{{G_{c0}}t} \over {{L_c}{\rho _1}}},\hat z = {z \over {{L_c}}},\hat p = {p \over {G_{co}^2/{\rho _1}}} \cr &{R_{gl}} = {{{\rho _g}} \over {{\rho _1}}},{{\hat v}_{gj}} = {{{v_{gj}}{\rho _1}} \over {{G_{co}}}},{N_{sub}} = {{{c_p}\Delta {T_{in}}} \over {{h_{gl}}}}\left( {{{{\rho _1}} \over {{\rho _g}}} - 1} \right) \cr &Fr = {{{{\left( {{G_{co}}/{\rho _1}} \right)}^2}} \over {g{L_c}}},{N_{pch}} = {{{P_e}q{L_C}} \over {A{G_{co}}{h_{gl}}}}\left( {{{{\rho _1}} \over {{\rho _g}}} - 1} \right) \cr &{{\hat h}_f} = {{{h_f} - {h_{fa}}} \over {{h_{fb}} - {h_{fa}}}},{N_f} = {{{h_{fb}} - {h_{fa}}} \over {{h_{gl}}}}\left( {{{{\rho _1}} \over {{\rho _g}}} - 1} \right) \cr} $

(9)

式中：L为长度，c_p为液相定压比热容，ΔT_in为加热段入口过冷度。将上述无量纲参数代入方程(1)~(3)，得到如下的无量纲控制方程：

$ {\partial \over {\partial \hat t}}\left[ {\left( {1 - \alpha } \right)\alpha {R_{gl}}} \right] + {{\partial \hat G} \over {\partial z}} = 0$

(10)

$ \eqalign{ &{{\partial \hat G} \over {\partial \hat t}} + {\partial \over {\partial z}}\left[ {{{{{\left( {1 - x} \right)}^2}} \over {\left( {1 - \alpha } \right)}} + {{{x^2}} \over {{R_{gl}}\alpha }}} \right]{{\hat G}^2} = {{\partial \hat P} \over {\partial \hat z}} - \cr &\left[ {\left( {1 - \alpha } \right) + \alpha {R_{gl}}} \right]{1 \over {Fr}} - {1 \over 2}\xi {{{{\hat G}^2}} \over {\hat D}}\varphi _f^2 \cr} $

(11)

$ \eqalign{ &{\partial \over {\partial \hat t}} + \left[ {\left( {1 - \alpha } \right){{{N_f}} \over {{N_{pch}}}}{{\hat h}_f} + {{\alpha \left( {1 - {R_{gl}}} \right)} \over {{N_{pch}}}}} \right] + \cr &{\partial \over {\partial \hat z}}\left[ {\left( {1 - x} \right){{{N_f}} \over {{N_{pch}}}}{{\hat h}_f}\hat G + {{x{{\hat G}_g}\left( {1 - {R_{gl}}} \right)} \over {{R_{gl}}{N_{pch}}}}} \right] = {{\hat q} \over {\hat A}} \cr} $

(12)

为了对上述的无量纲控制方程进行简化，对其进行量纲分析。结果发现由于R_gl数值接近0.001，导致αR_gl等项与其余各项相比要小2~3个数量级，可以认为是小量而忽略。此时控制方程简化为

$ {\partial \over {\partial \hat t}}\left( {1 - \alpha } \right) + {{\partial \hat G} \over {\partial z}} = 0$

(13)

$ \eqalign{ &{{\partial \hat G} \over {\partial \hat t}} + {\partial \over {\partial z}}\left\{ {\left[ {{{{{\left( {1 - x} \right)}^2}} \over {\left( {1 - \alpha } \right)}} + {{{x^2}} \over {{R_{gl}}\alpha }}} \right]{{\hat G}^2}} \right\} \cr & = {{\partial \hat P} \over {\partial \hat z}} - {{1 - \alpha } \over {Fr}} - {1 \over 2}\xi {{{{\hat G}^2}} \over {\hat D}}\varphi _f^2 \cr} $

(14)

$ \eqalign{ &{\partial \over {\partial \hat t}}\left[ {\left( {1 - \alpha } \right){{{N_f}} \over {{N_{pch}}}}{{\hat h}_f} + {\alpha \over {{N_{pch}}}}} \right] = {{\hat q} \over {\hat A}} - \cr &{\partial \over {\partial \hat z}}\left[ {\left( {1 - x} \right){{{N_f}} \over {{N_{pch}}}}{{\hat h}_f}\hat G + {{x{{\hat G}_g}} \over {{R_{gl}}{N_{pch}}}} - {{x{{\hat G}_g}} \over {{N_{pch}}}}} \right] \cr} $

(15)

同理对式(4)所示的辅助方程和式(5)~(8)所示的边界方程无量纲化和量级分析，得到无量纲的辅助方程和边界方程，具体过程不再冗述。

2 系统状态空间表达 2.1 线性增量方程的提出

假定系数具有线性特性，通过外加线性小扰动的方法，对系统方程进行线性化，得到描述系统动态特性的线性增量方程。系统变量采用下式表示：

$ \left\{ \matrix{ \hat G = {{\hat G}_0} + \Delta \hat G\exp \left( {\hat s\hat t} \right) \hfill \cr \hat P = {{\hat P}_0}\Delta \hat P\exp \left( {\hat s\hat t} \right) \hfill \cr x = {x_0}\Delta x\exp \left( {\hat s\hat t} \right) \hfill \cr \alpha = {\alpha _0}\Delta \alpha \exp \left( {\hat s\hat t} \right) \hfill \cr} \right.$

(16)

式中：$ {\hat s} $为无量纲化的Laplace变化系数，其定义式如下

$ \hat s = s{L_C}{\rho _1}/{G_{co}}$

(17)

将式(16)代入式(13)~(15)，得到线性增量的控制方程，其中单相区域的表达形式如下

$ {{d\Delta \hat G} \over {dz}} = 0$

(18)

$ \hat s\Delta \hat G = {{d\Delta \hat P} \over {d\hat z}} - \xi {{\hat G} \over {\hat D}}\Delta \hat G$

(19)

$ \eqalign{ &\left[ {{{\Delta x} \over {{R_{gl}}{N_{pch}}}} + {{{N_f}} \over {{N_{pch}}}}{{{{\hat h}_f}} \over {d\hat P}}\left| {_{{{\hat P}_0}}} \right.\Delta \hat P} \right]\hat s + {{1 - {R_{gl}}} \over {{R_{gl}}{N_{pch}}}}{\partial \over {\partial \hat z}} \cr &\left[ {{{\hat G}_0} + \Delta x + {x_0}\Delta \hat G} \right] + {{{N_f}} \over {{N_{pch}}}}{\partial \over {\partial \hat z}}\left[ {{{{{\hat h}_f}} \over {d\hat P}}\left| {_{{{\hat P}_0}}{{\hat G}_0}\Delta \hat P + {{\hat h}_{fo}}\Delta \hat G} \right.} \right] = 0 \cr} $

(20)

而在两相区域内，线性扰动方程的表达形式为

$ d\Delta \hat G - \hat s\Delta \alpha = 0$

(21)

$ \eqalign{ &\hat s\Delta \hat G + {d \over {d\hat z}}\left[ {\left( {{{1 - {x_0}^2} \over {1 - {\alpha _0}}}} \right)\hat G_0^2\left( {{{\Delta \alpha } \over {1 - {\alpha _0}}}} \right) - {{2\Delta x} \over {1 - {x_0}}} + {{2\Delta \hat G} \over {{{\hat G}_0}}} + {{x_0^2\hat G_0^2} \over {{R_{gl}}{\alpha _0}}}\left( {{{2\Delta x} \over {{x_0}}} + {{2\Delta \hat G} \over {{{\hat G}_0}}} - {{\Delta \alpha } \over {{\alpha _0}}}} \right)} \right] \cr & = {{\Delta \alpha } \over {{F_R}}} + {{d\Delta \hat P} \over {d\hat z}} - {{\xi x_0^2\hat G_0^2} \over {2D}}\left[ {{{2\Delta \hat G} \over {{{\hat G}_0}}} + {1 \over {\varphi _{f0}^2}}{{d\Delta \varphi _f^2} \over {dx}}\left| {_{_{{{\hat P}_0}}}\Delta x + {1 \over {\varphi _{f0}^2}}{{d\varphi _f^2} \over {d\hat P}}\left| {_{_{{{\hat P}_0}}}\Delta \hat P} \right.} \right.} \right] \cr} $

(22)

$ \eqalign{ &{{{N_f}\hat s} \over {{N_{pch}}}}\left( {1 - \alpha } \right){{\hat h}_{f0}}\left[ {{{d{{\hat h}_f}} \over {{{\hat h}_{f0}}d\hat P}}\left| {_{{{\hat P}_0}}\Delta \hat P - {{\Delta \alpha } \over {1 - {\alpha _0}}}} \right.} \right] + {{\Delta \alpha \hat s} \over {{N_{pch}}}} + \cr &{{{N_f}} \over {{N_{pch}}}}{\partial \over {\partial \hat z}}\left\{ {\left( {1 - {x_0}} \right){{\hat h}_{f0}}{{\hat G}_0}\left[ {{{\Delta \hat G} \over {{{\hat G}_0}}} + {{d{{\hat h}_f}} \over {{{\hat h}_{f0}}d\hat P}}} \right]\left| {_{{{\hat P}_0}}\Delta \hat P - {{\Delta x} \over {1 - {x_0}}}} \right.} \right\} \cr & + {{1 - {R_{gl}}} \over {{R_{gl}}{N_{pch}}}}{\partial \over {\partial \hat z}}\left[ {{{\hat G}_0}\Delta x + {x_0}\Delta \hat G} \right] = 0 \cr} $

(23)

同时，将方程(16)代入无量纲的辅助方程，得到线性增量的辅助方程：

$ \eqalign{ &\Delta \alpha = {{{{ - {x_0}} \over {{x_0} + \left( {1 - {x_0}} \right){R_{gl}}}}} \over {{C_0} + \left( {{{{{\hat v}_{gj0}}{R_{gl}}} \over {{x_0} + \left( {1 - {x_0}} \right){R_{gl}}}}} \right){1 \over {{{\hat G}_0}}}}} \cr &\left[ {{{\Delta x} \over {{x_0} + \left( {1 - {x_0}} \right){R_{gl}}}} - {{\Delta x} \over {{x_0}}} - {1 \over {{C_0}}}\left( {{{\Delta \hat G} \over {{{\hat G}_0}}} + {{\Delta x} \over {{x_0} + \left( {1 - {x_0}} \right){R_{gl}}}} - {{\Delta \alpha } \over {{{\hat v}_{gj0}}}}{{d{{\hat v}_{gj}}} \over {d\alpha }}} \right)} \right] \cr} $

(24)

另外，将方程(16)代入无量纲的边界方程中，得到线性增量的边界方程：

加热通道入口位置($ {\hat z} $= 0)：

$ \Delta {\hat P_{c,in}} = - \left( {1 + {\xi _{c,in}}} \right){\hat G_0}\Delta G$

(25)

上升通道入口位置($ {\hat z} $= 1)：

$ \Delta {\hat P_{c,ex}} = \Delta {\hat P_b} - {\xi _b}{\hat G_0}\Delta \hat G$

(26)

上升通道出口位置($ {\hat z} $= 1+$ {\hat L} $_r)

$ \Delta {\hat P_{r,ex}} = - \left( {{\xi _{r,rx}} - 1} \right)\varphi _{f0}^2{\hat G_0}\Delta G$

(27)

2.2 状态空间方程组的建立

根据线性化扰动的假设(即外加扰动$ \Delta \hat G $、$ \Delta \hat P $、Δx等随位置$ {\hat z} $的变化而呈线性变化)，在加热段和上升段(0≤$ {\hat z} $≤1+$ {\hat L} $_r)的范围内对线性增量的控制方程进行积分，其中闪蒸起始点以下区域采用单相的线性增量控制方程，闪蒸点以上区域采用两相的线性增量控制方程。积分时变量的边界条件参照线性增量的边界方程。

同时为了使变量和方程的数量相同，需要将线性增量的辅助方程式(24)代入控制方程，用以消去Δα这个变量，最终得到一组3×3的线性常微分方程组，即描述自然循环闪蒸不稳定动态特性的状态空间表达式：

$ MXs = KX + BU $

(28)

式中： X 为由线性扰动ΔG、ΔP、Δx组成的状态向量，即$，而 M 和 K 为线性常微分方程组的系数矩阵：

$ M = \left[ {\matrix{ {{a_{11}}} &{{a_{12}}} &{{a_{13}}} \cr {{a_{21}}} &{{a_{22}}} &{{a_{23}}} \cr {{a_{31}}} &{{a_{32}}} &{{a_{33}}} \cr } } \right],K = \left[ {\matrix{ {{b_{11}}} &{{b_{12}}} &{{b_{13}}} \cr {{b_{21}}} &{{b_{22}}} &{{b_{23}}} \cr {{b_{31}}} &{{b_{32}}} &{{b_{33}}} \cr } } \right]$

(29)

其中，矩阵 M 和 N 的参数如下所示，为了方便表述，首先定义如下的参数：

$ \eqalign{ &{k_1} = {C_0}\left( {{{{R_{gl}}} \over {{x_0}}} - {R_{gl}} + 1} \right) + {{{{\hat v}_{gj0}}{R_{gl}}} \over {{{\hat G}_0}}}, \cr &{k_2} = {1 \over {{C_0}{{\hat v}_{gj0}}}}{{d{{\hat v}_{gj}}} \over {d\alpha }}\left| {_{{\alpha _0}}} \right. \cr &{k_3} = {{\left( {1 - {x_0}} \right){R_{gl}}{C_0} + {x_0}} \over {\left[ {{x_0} + \left( {1 - {x_0}} \right){R_{gl}}} \right]{C_0}{x_0}}}, \cr &{k_4} = {1 \over {{C_0}{{\hat G}_0}}} \cr &{k_5} = {{\left( {2 + 2{{\hat L}_r} - {{\hat z}_{fla}}} \right){{\hat z}_{fla}}} \over {2 + 2{{\hat L}_r}}}, \cr &{k_6} = {{{{\hat z}_{fla}}{{\left( {1 - {x_0}} \right)}^2}{{\hat G}_0}^2} \over {\left( {1 + {{\hat L}_r}} \right){{\left( {1 - {x_0}} \right)}^2}}} \cr &{k_7} = {{{\xi ^2}{{\left( {1 + {{\hat L}_r} - {{\hat z}_{fla}}} \right)}^2}\hat G} \over {2 + 2{{\hat L}_r}\hat D}}, \cr &{k_8} = {{{x_0}^2{{\hat G}_0}^2} \over {\left( {1 + {{\hat L}_r}} \right){R_{gl}}{\alpha _0}}} \cr &{k_9} = {{\left( {1 + {{\hat L}_r} - {{\hat z}_{fla}}} \right){{\hat z}_{fla}}} \over {\left( {2 + 2{{\hat L}_r}} \right){F_R}}}, \cr &{k_{10}} = {1 \over {\varphi _{f0}^2}}{{d\varphi _f^2} \over {dx}}\left| {_{_{_{_{_{{{\hat p}_0}}}}}}} \right. \cr &{k_{11}} = \left[ {1 + {\xi _{c,in}} - {\xi _b} + \left( {{\xi _{r,ex}} - 1} \right)\varphi _{f0}^2} \right]{{\hat G}_0} \cr &{k_{12}} = {{{{\left( {1 + {{\hat L}_r} - {{\hat z}_{fla}}} \right)}^2}} \over {2{R_{gl}}{N_{pch}}}} \cr &{k_{13}} = {{{N_f}\left( {1 - {\alpha _0}} \right)} \over {{N_{pch}}}} \cdot {{d{{\hat h}_f}} \over {d\hat p}}\left| {_{_{_{_{_{_{{{\hat p}_0}}}}}}}} \right. \cr &{k_{14}} = {{{N_f}{{\left( {1 + {{\hat L}_r} - {{\hat z}_{fla}}} \right)}^2}} \over {\left( {2 + 2{{\hat L}_r}} \right){N_{pch}}}}{{d{{\hat h}_f}} \over {d\hat p}}\left| {_{_{_{_{_{_{{{\hat p}_0}}}}}}}} \right. \cr &\left| {{k_{15}} = {{{N_f}{{\hat h}_{f0}}} \over {{N_{pch}}}}} \right. \cr &{k_{16}} = {{{N_f}\left( {1 + {{\hat L}_r} - {{\hat z}_{fla}}} \right){{\hat G}_0}} \over {{N_{pch}}\left( {2 + 2{{\hat L}_r}} \right)}}{{d{{\hat h}_f}} \over {d\hat p}}\left| {_{_{_{_{_{_{{{\hat p}_0}}}}}}}} \right. \cr &{k_{17}} = {{\xi {{\hat G}_0}^2\varphi _{f0}^2} \over {2D}} \cr &{k_{18}} = {{{N_f}{{\hat h}_{f0}}\left( {1 + {{\hat L}_r} - {{\hat z}_{fla}}} \right){{\hat G}_0}} \over {{N_{pch}}\left( {1 + {{\hat L}_r}} \right)}} \cr &{k_{19}} = {1 \over {\varphi _{f0}^2}}{{d\varphi _f^2} \over {d\hat p}}\left| {_{_{_{_{_{_{{{\hat p}_0}}}}}}}} \right. \cr &{k_{20}} = {{{N_f}{{\hat h}_{f0}}\left( {1 - {x_0}} \right){{\hat G}_0}} \over {{N_{pch}}}}{{{{\hat z}_{fla}}} \over {\left( {1 + {{\hat L}_r}} \right)}} \cr} $

此时矩阵 M 的参数可以表示为

$ \eqalign{ &{a_{11}} = {{{k_3}{k_5}} \over {{k_1} + {k_4}}},{a_{12}} = 0,{a_{13}} = {{{k_2}{k_5}} \over {{k_1} + {k_4}}}, \cr &{a_{21}} = {{1 + {{\hat L}_r}} \over 2},{a_{22}} = 0,{a_{23}} = 0, \cr &{a_{31}} = {k_{12}} + \left( {{{{k_5}} \over {{N_{pch}}}} - {k_5}{k_{15}}} \right){{{k_3}} \over {{k_1} + {k_4}}} \cr &{a_{32}} = {k_{14}} + {k_5}{k_{13}}, \cr &{a_{33}} = \left( {{{{k_5}} \over {{N_{pch}}}} - {k_5}{k_{15}}} \right){{{k_4}} \over {{k_1} + {k_4}}} - {k_{11}}\left( {{k_{14}} + {k_5} \cdot {k_{13}}} \right) \cr} $

而矩阵 N 的参数则可以表示为

$ \eqalign{ &{b_{11}} = 1,{b_{12}} = 0,{b_{13}} = {{{k_2}{k_5}} \over {{k_1} + {k_4}}}, \cr &{b_{21}} = {{{k_3}{k_8}} \over {\left( {{k_1} + {k_4}} \right){\alpha _0}}} - {{2{k_8}} \over {{x_0}}} - {{{k_3}{k_6}} \over {\left( {{k_1} + {k_4}} \right)\left( {1 - {\alpha _0}} \right)}} + \cr &{{{k_3}{k_9}} \over {{k_1} + {k_4}}} + {{2{k_6}} \over {1 - {x_0}}} - {k_5}{k_{10}}{k_{17}}, \cr &{b_{22}} = - 1 - {k_5}{k_{17}}{k_{19}}, \cr &{b_{23}} = {{{k_2}{k_8}} \over {\left( {{k_1} + {k_4}} \right){\alpha _0}}} - {{{k_2}{k_6}} \over {\left( {{k_1} + {k_4}} \right)\left( {1 - {\alpha _0}} \right)}} + \cr &{{{k_2}{k_9}} \over {{k_1} + {k_4}}} + {{2\left( {{k_{17}} - {k_6} - {k_8}} \right)} \over {{{\hat G}_0}}} - {k_7} + {k_{11}} + {k_5}{k_{17}}{k_{19}}{k_{11}}, \cr &{b_{31}} = {{{k_{20}}} \over {1 - {x_0}}} - {{{{\hat G}_0}} \over {{R_{gl}}{N_{pch}}}} \cr &{b_{32}} = - {k_{16}} - {k_{20}}{{d{{\hat h}_f}} \over {{{\hat h}_{f0}}d\hat P}}\left| {_{_{_{_{_{_{{{\hat P}_0}}}}}}}} \right., \cr &{b_{33}} = {k_{11}}{k_{16}} - {k_{18}} - {{{x_0}} \over {{R_{gl}}{N_{pch}}}} - {{{k_{20}}} \over {{{\hat G}_0}}} - {k_{11}}{k_{20}}{{d{{\hat h}_f}} \over {{{\hat h}_{f0}}d\hat P}}\left| {_{_{_{_{_{_{{{\hat P}_0}}}}}}}} \right. \cr} $

2.3 系统空间状态方程稳定域的判断

对于方程(28)所示的状态空间方程组，可以通过判断系数矩阵 N (= M ^T K )的特征根实部的符号来判断系统的稳定性。即求解如下的特征方程：

$ \det \left( {\lambda I - N} \right) = 0$

(30)

其特征根为

$ {\lambda _i} = {\gamma _i} + j{\mu _i} i = 1,2,3$

(31)

当 γ_i 均小于零时，系统可以维持稳定，即表示开式自然循环系统能够形成稳定的闪蒸；而 γ_i 存在大于零的情况时，系统无法维持稳定，即此时的闪蒸流动为闪蒸不稳定工况。

3 模型计算结果及与实验结果比较

由于上述计算过程繁复，较难推导出闪蒸不稳定区域所对应的数学表达形式，因此本文利用MATLAB 7.0程序对方程特征值进行求解和对系统稳定性进行判断，其计算过程如图 2所示。

图 3给出了模型的计算结果，图中实线表示计算得到的闪蒸不稳定边界，数据点表示开式自然循环实验台上得到的闪蒸不稳定区域(实验装置见图 4，实验方法见文献^[11])。通过对比模型计算结果与实验结果发现，两者之间趋势上吻合较好。

其中在高加热功率工况下，计算得到的不稳定上边界与实验结果存在一定偏差。此时实验观察到加热管出口过冷沸腾现象强烈而形成间歇喷泉流动，而本模型中未考虑过冷沸腾现象对闪蒸流动的影响，因此偏差产生的原因极可能是由于过冷沸腾过程的影响。

同样在高加热功率工况下，计算得到的不稳定下边界与实验结果也存在偏差。此时随着加热段入口过冷度降低，使得加热段内沸腾喷发过程增强，并与上升段内闪蒸流动过程共同影响自然循环的流动特性^[11]。而沸腾喷发过程会导致流量、上升段内的温度分布等出现剧烈地振荡，进而导致持续的稳定闪蒸流动无法形成，即可以认为加热段内沸腾喷发过程对于稳定闪蒸过程起到强烈的干扰作用。

因此本实验采用注气诱导法排除沸腾喷发对于稳定闪蒸过程的干扰作用。即在闪蒸不稳定工况时，向上升段内持续注入气体，使得系统形成稳定的空气-水两相流动，此时自然循环流量快速增加，加热段内沸腾喷发过程被抑制。随后逐渐减少注气量直至停止，这一过程中闪蒸流动过程会逐渐主导系统的流动特性，而由于循环流量较高导致该过程中加热段内沸腾喷发不会出现。通过这种方法来排除沸腾喷发对闪蒸过程的影响，再通过观察停止注气后自然循环闪蒸流动的发展演化过程来判断闪蒸流动是否稳定。图 5给出了典型的注气诱导法过程及其结果。

图 5(a)所示的实验工况位于模型计算的闪蒸不稳定区域内(图 3的不稳定区域)，注气前系统处于闪蒸和沸腾共同诱发的流动不稳定状态。当注气逐渐停止后，加热管内沸腾始终被抑制，而上升段内闪蒸流动出现振荡，振幅不断增大直至最终恢复为原来的流动不稳定状态，即说明此工况下闪蒸流动过程为发散过程。而图 5(b)所示实验工况则处于模型计算的稳定闪蒸流动区域(图 3的下部稳定区域)，注气前系统与工况1相似，但注气逐渐停止后，上升段内闪蒸流动首先出现较弱振荡，随后振荡逐渐减弱并形成稳定闪蒸流动，此时稳定闪蒸驱动的流量很大，使得加热段出口为单相流动，即工况2时自然循环闪蒸流动过程为自发收敛过程，也证实了图 3实验中未出现稳定闪蒸流动主要是由于沸腾喷发过程的干扰作用。

图 6给出了通过重复实验得到的注气诱导法形成的稳定闪蒸流动边界，本文认为是准确的闪蒸不稳定的区域下边界，由图可知该边界与计算得到的不稳定下边界符合很好，即证明了本文模型结果具有很好的准确性。另外，可以认为本模型计算得到下稳定区域内的工况，都可以通过注气诱导法将闪蒸不稳定流动转化为持续稳定闪蒸流动。

4 结论

1) 本文运用状态空间法对开式自然循环闪蒸不稳定进行研究，提出了描述该过程的线性均相流模型，并导出了描述系统稳定性的状态空间方程。

2) 通过求解状态空间方程得到了自然循环闪蒸不稳定的边界，将计算结果与实验结果对比发现两者趋势上吻合较好。其中上边界的偏差是由于模型未考虑过冷沸腾的影响而造成的；下边界的偏差是由加热段内沸腾喷发对稳定闪蒸的干扰作用而产生的。

3) 通过注气诱导法排除了沸腾喷发干扰作用，得到了准确的闪蒸不稳定下边界，该边界与计算得到的不稳定下边界吻合很好。

4) 通过实验发现，计算得到的下稳定区域内的工况点，均可通过注气诱导法将闪蒸不稳定工况转化为稳定闪蒸流动工况。