CN104899387A - 一种蒸汽发生器管束两相横向流弹失稳分析方法 - Google Patents

Abstract

本发明涉及一种蒸汽发生器管束两相横向流弹失稳分析方法。本发明是为了克服核反应堆蒸汽发生器换热管流弹失稳而导致的管束破坏问题，而提出一种蒸汽发生器管束两相横向流弹失稳分析方法。该方法是在Liming提出的单相横向流弹失稳显示分析模型-“流管”模型的基础上，提出针对两相横向流弹失稳的显式解析模型，然后确定求解该解析模型所需要的两相流计算参数，最后求解解析模型得到蒸汽发生器管束临界折合流速和质量阻尼比，从而得到稳定区图。一个确定的流管模型，质量阻尼比是确定，只需要计算出具体折合流速，然后判断两者确定的点，若在稳定区图内，则说明蒸汽发生器管束不会发生流弹失稳，反之，则说明蒸汽发生器管束将发生流弹失稳。

Description

一种蒸汽发生器管束两相横向流弹失稳分析方法

技术领域

本发明涉及一种蒸汽发生器管束两相横向流弹失稳分析方法，提出的两相横向流弹失稳的解析模型，分析核反应堆中蒸汽发生器管束两相横向流弹失稳是否发生的一种方法。

背景技术

流弹失稳机理与航天工业中广泛研究的颤振现象相同，不稳定性的主要参数是系统的阻尼和流体弹性力。横向流作用下的管束，当流速增加到一定值时，系统中流体弹性力所做的功将大于系统阻尼散耗的功，因此管子振动的振幅将急剧增大，这时就称出现了流弹失稳。

从理论上看，流体弹性力应该通过求解一组联立的流体力学和结构动力学方程得到。但由于我们涉及的问题通常十分复杂，如大量的换热管，复杂的边界条件和高雷诺数的流动等，使得从理论上进行求解十分困难，甚至是不可能的。通常的做法是根据大量的实验结果，提出简化的数学物理模型。在这些模型中，流体力一般不再以独立方程的形式出现而是表示成与结构的位移、速度和加速度等有关的函数，这种关系体现了流体作用使结构产生位移和变形，而结构的位移和变形又改变了流体作用这样一个流固耦合振动的事实。

目前，已有分析管阵流弹失稳的模型，主要为拟静力(如Connors模型)、拟定常(如Paidoussis模型)、非定常流模型(如Chen模型)及解析模型(Yetisir流管模型)。分析管阵临界流速时，拟静力、拟定常和非定常流模型的流体弹性力都不同程度地依赖于实验数据，且这些模型中的一些参数是在一定条件下、针对某个管阵形式得到的，有一定的保守性和不确定性，并不能对所有实际问题给出可靠的临界流速。而解析模型则是根据大量的实验结果，对流体弹性系统作一些简化，用比较成熟的数学、力学方法，得到流体弹性力。尽管这类模型与实际情况还有一定差距，有待于进一步完善，但无疑是一类具有普遍意义的数学模型。

现有的解析方法计算过程比较繁琐，甚至需用到数值计算，且都是针对单相流弹失稳的分析，如果用其分析两相流弹失稳问题，还需将气-液两相流等效为单相流，这一方法的有效性也有待考证。

发明内容

本发明的目的是：本发明是为了克服核反应堆蒸汽发生器换热管流弹失稳而导致的管束破坏问题，而提出一种蒸汽发生器管束两相横向流弹失稳分析方法。包括以下步骤：

1两相横向流弹失稳参数确定

热交换器通常有四种典型的管阵几何排列形式，对于所有的几何排列形式，解析模型的基本概念是相同的，排列形式只是影响模型中使用的试验参数。具体以平行三角形管束为例，而对于旋转三角形和旋转正方形管束(P_r>1.7)，流体力F_y具有同样的形式，除了它需要除以2cosα，这是由于所定义的最小间距扰动面积函数a(s_m,t)具有不同的形式。

1.1基本定义

引入“流管”概念，平行三角形管阵流管模型。曲线坐标s用以表示流动路径，流管的平均横截面积在整个流管上假设为常数，且由-s₁位置的最小间隙面积A₀确定。

A₀＝min(Pcosα-D/2,P-D)

瞬时流管面积A(s,t)会随着管子的振动而改变，A(s,t)包含两项：平均项A₀，和波动项a(s,t)。

A(s,t)＝A₀+a(s,t)

同样的，速度和压力定义为：

U(s,t)＝U₀+u(s,t)，P(s,t)＝P₀+p(s,t)

然而，在足够远的上游位置s₁，扰动可以忽略，速度和压力可以分别视为常数U₀和P₀，流管假设在s_a位置开始附着于柔性管，并在s_s位置分离。

管子的振动不能立刻影响到流管的其它位置，因此，采用一个相函数来考虑扰动延迟效应。上游微扰函数a(s,t)可以表示为

f(s)为人工衰减函数，用以表示从柔性管的位置到上游点的微扰衰减。a(s_m,t)是最小间隙位置的面积微扰，是管子几何的函数。

1.2两相横向流弹失稳试验参数确定

唯一依赖于试验输入的参数为开始产生附着和分离点的位置，这可以通过Weaver和AbdRabbo以及Scott的可视化研究进行估计。四种典型的管阵排列的可以归为以下两组：

注：β为流体附着表面的法向矢量角度

1.3两相横向流弹失稳计算参数确定

1.3.1流场计算参数

(1)空泡份额

根据均相流模型(HEM)假设，两相流的两相充分混合，并且密度和温度分布均匀，不存在气相和液相的相对滑移，可将空泡份额定义为：

${\mathrm{\α}}_H=1/(1+\frac{1-X_a}{X_a}\frac{{\mathrm{\ρ}}_g}{{\mathrm{\ρ}}_l})$

上式中，X_a为真实含汽率，ρ_g和ρ_l分别为气相和液相的密度。

(2)流动参数

节距流速U_p：U_p＝U_∞P/(P-D)

均相流密度ρ：ρ＝ρ_l(1-α_H)+ρ_gα_H

自由来流速度U_∞： $U_{\∞}=\frac{{\mathrm{\ρ}}_l{V}_l+{\mathrm{\ρ}}_g{V}_g}{\mathrm{\ρA}}$

自由来流质量流速：

节距质量流速 ${\stackrel{\·}{m}}_p={\stackrel{\·}{m}}_{\∞}\frac{P}{(P-D)}=\mathrm{\ρ}{U}_p$

式中，U_∞是自由来流速度(即假设不存在管子时的速度)，P是管中心距，D是管子直径(对于翅片管，用等效水力直径Dh代替)，A为自由来流的流道截面积，ρ_g为气相密度、ρ_l为气相密度，V_g为气相速度、V_l为气相速度。

1.3.2动水参数

(1)水力直径

对于非翅片管，水力直径就简单的为管子外径D。对于翅片管，水力直径等效为D_h：

D_h＝D_r+R_F(D_o-D_r)

式中，D_r为齿根圆直径，D_o为翅片外径，比值R_F为翅片占据的体积与管子外表面到翅片外径表面的有效体积的比。对于翅片管，D_h要用于振动激励参数的计算，等效节径比P/D也要变为P/D_h。

(2)水动力附加质量

管束中一根管子单位长度的水动力附加质量可以表示为：

$m_h=\left(\frac{\mathrm{\π}}{4}{\mathrm{\ρ}}_{\mathrm{TP}}{D}^2\right)\left[\frac{{(D_e/D)}^2+1}{{(D_e/D)}^2-1}\right]$

式中，D_e是周围管子的等效直径，比值D_e/D是管束对管子约束程度的度量。

对于三角形管阵：

D_e/D＝(0.96+0.5P/D)P/D

对于正方形管阵：

D_e/D＝(1.07+0.56P/D)P/D

单位长度管子总的动力质量m包括水动力附加质量m_h、管子本身的质量m_t和管内流体的质量m_i：

m＝m_h+m_t+m_i

(3)阻尼

两相流中多跨热交换器管的阻尼包括支撑阻尼ζ_S、粘性阻尼ζ_V和两相阻尼ζ_TP：

ζ_T＝ζ_S+ζ_V+ζ_TP

通量很高、空泡份额很大的情况下，支撑就可能是干的，此时只有摩擦阻尼存在，这种情况下，支撑阻尼与气体中的热交换器管阻尼基本相同：

ζ_T＝ζ_S+ζ_V+ζ_TP

管和支撑之间存在液体时，支撑阻尼ζ_S既包括压膜阻尼ζ_SF又包括摩擦阻尼ζ_F：

${\mathrm{\ζ}}_S={\mathrm{\ζ}}_{\mathrm{SF}}+{\mathrm{\ζ}}_F=\left(\frac{N-1}{N}\right)[\frac{1460}{f}\left(\frac{{\mathrm{\ρ}}_l{D}^2}{m}\right)+0.5]{\left(\frac{L}{l_m}\right)}^{1/2}$

粘性阻尼：

${\mathrm{\ζ}}_V=\frac{100\mathrm{\π}}{\sqrt{8}}\left(\frac{{\mathrm{\ρ}}_{\mathrm{TP}}{D}^2}{m}\right){\left(\frac{2{\mathrm{\υ}}_{\mathrm{TP}}}{\mathrm{\π}{\mathrm{fD}}^2}\right)}^{1/2}\left\{\frac{[1+{(D/D_e)}^3]}{[1-{(D/D_e)}^2{]}^2}\right\}$

两相阻尼比：

${\mathrm{\ζ}}_{\mathrm{TP}}=4.0\left(\frac{{\mathrm{\ρ}}_l{D}^2}{m}\right)\left[f\left({\mathrm{\α}}_H\right)\right]\left\{\frac{[1+{(D/D_e)}^3]}{{[1-{(D/D_e)}^2]}^2}\right\}$

式中，ρ_l是管和支撑间液体的密度，D_e是周围管子的等效直径，D是管子的外径，m是单位长度管子的总质量，包括用两相流均相密度计算的水动力附加质量，l_m是管子的特征长度，L是支承板宽度，f是管子的频率，ρ_TP是气液两相等效密度，υ_TP是等效运动粘度，υ_TP＝υ_l/[1+α_H(υ_l/υ_g-1)]，υ_g和υ_l分别为液相和气相的运动黏度，而f(α_H)的确定

f(α_H)＝α_H/40         α_H<40％

f(α_H)＝1              40％≤α_H≤70％_。

f(α_H)＝1-(α_H-70)/30  α_H>70％

2在Liming提出的单相横向流弹失稳显式分析模型-“流管”模型的基础上，提出针对两相横向流弹失稳的解析模型。热交换器通常有四种典型的管阵几何排列形式，对于所有的几何排列形式，解析模型的基本概念是相同的，排列形式只是影响模型中使用的试验参数。对于平行三角形管束模型：

2.1不稳定流速

式中Y，很多研究已经证明流体弹性不稳定通常首先在升力方向(y方向)产生。

假设柔性管在两自由度方向以相同的固有频率做简谐运动

$m\stackrel{\&CenterDot;\&CenterDot;}{x}+c\stackrel{\&CenterDot;}{x}+\mathrm{kx}=F_x{a}_1(s_m,t)+F_x{a}_2(s_m,t)$

$m\stackrel{\&CenterDot;\&CenterDot;}{y}+c\stackrel{\&CenterDot;}{y}+\mathrm{ky}=F_y{a}_1(s_m,t)-F_y{a}_2(s_m,t)$

其中，a₁(s_m,t)和a₂(s_m,t)分别是流管1和流管2在最小间隙处的面积微扰(面积变小为正，变大为负。)，它们由管子几何参数定义。

$a_1(s_m,t)=-\frac{1}{2}y\left(t\right),a_2(s_m,t)=\frac{1}{2}y\left(t\right)$

方程的解为y(t)＝Ye^iωt

式中，ω是振动复频率。U_r＝U₀/ω_nl₀是缩减速度，l₀＝2s_l，频率比ω_r＝ω/ω_n，ω_n是管子在静水中的固有频率，s_l是从振动管到压力扰动可以忽略的位置的距离。

由于扰动衰减和时间衰减仅发生在“流管”路径上，以下给出两组针对平行三角形管束的面积扰动方程

$f\left(s\right)=\frac{s-s_1}{s_a-s_1},s_1\&le;s\&le;s_a$

f(s)＝1,s_a≤s≤s_s

而管子运动和流体力之间的时间延迟通过相位函数来表示

为了便于使用，在接下来的公式中将写为

$\stackrel{\&OverBar;}{s}=\frac{s-s_a}{s_1-s_a},$ 则衰减函数 $f\left(s\right)=1-\stackrel{\&OverBar;}{s}$

$U_r=\frac{U_p}{\mathrm{fD}}\frac{1}{l_0{A}_0}\frac{\cos \mathrm{\&alpha;}}{2\mathrm{\&pi;}}(\mathrm{Pr}-1)$

其中：

$A_4=-\frac{{\mathrm{\&omega;}}_{r}i}{l_0{A}_0{U}_r}(s_l-s_a)$

2.2非稳态压力的计算

各分项计算：

$\begin{array}{c}\underset{s_a}{\overset{s}{\&Integral;}}{u}_2(s,U_r)\mathrm{ds}=(s_1-s_a)\underset{0}{\overset{\stackrel{\&OverBar;}{s}}{\&Integral;}}(A_4\stackrel{\&OverBar;}{s}+A_5)d\stackrel{\&OverBar;}{s}\\ =(s_1-s_a)(\frac{A_4}{2}{\stackrel{\&OverBar;}{s}}^2+A_5\stackrel{\&OverBar;}{s})\\ =(B_2{\stackrel{\&OverBar;}{s}}^2+B_3\stackrel{\&OverBar;}{s})\end{array}$

B₂＝(s₁-s_a)A₄/2

B₃＝(s₁-s_a)A₅

将以上代入非稳态压力的积分计算表达式得：

其中：

C₁＝-ω_rl₀U_riB₂

$C_2=-{\mathrm{\&omega;}}_r{l}_0{U}_{r}i{B}_3-2l_0^2{U}_r^2{A}_4-\frac{{\mathrm{\&omega;}}_r{l}_0{U}_{r}i}{A_0}(s_1-s_a)$

$\begin{array}{c}{C}_3=-{\mathrm{\&omega;}}_r{l}_0{U}_{r}i{B}_1-2l_0^2{U}_r^2{A}_5-\frac{l_0^2{U}_r^2}{A_0}\\ -\frac{{\mathrm{\&omega;}}_r{l}_0{U}_r}{A_0}{B}_4\end{array}$

2.3作用在管子上的流体力

$\begin{array}{c}{F}_y=\frac{Y}{D}{e}^{\mathrm{i\&omega;t}}\cos \mathrm{\&alpha;}\underset{s_a}{\overset{s_s}{\&Integral;}}{p}_2(s,U_r)\cos \left(\frac{2s}{\mathrm{Pr}}\right)\mathrm{ds}\\ =\frac{Y}{D}{e}^{\mathrm{i\&omega;t}}\cos \mathrm{\&alpha;}\underset{s_a}{\overset{s_s}{\&Integral;}}(C_1{\stackrel{\&OverBar;}{s}}^2+C_2\stackrel{\&OverBar;}{s}+C_3)\cos \left(\frac{2s}{\mathrm{Pr}}\right)\mathrm{ds}\end{array}$

计算各分项：

为便于计算令

$\begin{array}{c}\underset{s_a}{\overset{s_s}{\&Integral;}}\stackrel{\&OverBar;}{s}\cos \left(\frac{2s}{\mathrm{Pr}}\right)\mathrm{ds}=\underset{s_a}{\overset{s_s}{\&Integral;}}(\frac{s}{s_1-s_a}-\frac{s_a}{s_1-s_a})\cos \left(\mathrm{bs}\right)\mathrm{ds}\\ =\frac{1}{s_1-s_a}(\frac{1}{b^2}\cos \mathrm{bs}+\frac{s}{b}\sin \mathrm{bs}){|}_{s_a}^{s_s}-\frac{s_a}{s_1-s_a}\frac{1}{b}\sin \mathrm{bs}{|}_{s_a}^{s_s}=D_2\end{array}$

$\underset{s_a}{\overset{s_s}{\&Integral;}}\cos \left(\frac{2s}{\mathrm{Pr}}\right)\mathrm{ds}=\underset{s_a}{\overset{s_s}{\&Integral;}}\cos \left(\mathrm{bs}\right)\mathrm{ds}=\frac{1}{b}\sin \mathrm{bs}{|}_{s_a}^{s_s}=D_3$

将以上各项代入流体力的计算公式可得：

$\begin{array}{c}{F}_y=\frac{Y}{D}{e}^{\mathrm{i\&omega;t}}\cos \mathrm{\&alpha;}\underset{s_a}{\overset{s_s}{\&Integral;}}{p}_2(s,U_r)\cos \left(\frac{2s}{\mathrm{Pr}}\right)\mathrm{ds}\\ =\frac{Y}{D}{e}^{\mathrm{i\&omega;t}}\cos \mathrm{\&alpha;}\underset{s_a}{\overset{s_s}{\&Integral;}}(C_1{\stackrel{\&OverBar;}{s}}^2+C_2\stackrel{\&OverBar;}{s}+C_3)\cos \left(\frac{2s}{\mathrm{Pr}}\right)\mathrm{ds}\\ =\frac{Y}{D}{e}^{\mathrm{i\&omega;t}}\cos \mathrm{\&alpha;}\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}{C}_i{D}_i\end{array}$

F_y＝F_R+F_Ii是复数形式表示的流体力，将上式代入公式

$\frac{m}{\mathrm{\&rho;}{D}^2}=\frac{F_y(U_r,\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_n})}{1-{\left(\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_n}\right)}^2+i\frac{\mathrm{\&delta;}}{\mathrm{\&pi;}}\left(\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_n}\right)}$

$F_y={2l}_0{U}_r\cos \mathrm{\&alpha;}{F}_y\left(s\right)=F_R(U_r,\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_n})+i{F}_l(U_r,\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_n})$

中就可得到折合流速和质量阻尼比之间的关系，从而得到稳定区图。对于一个确定的流管模型，质量阻尼比是确定，只需要计算出具体折合流速，然后判断两者确定的点，若在稳定区图内，则说明蒸汽发生器管束不会发生流弹失稳，若两者确定的点在稳定区图外，则说明蒸汽发生器管束将发生流弹失稳。值得注意的是以上推导过程是根据平行三角形管阵推导而来，对于旋转三角形管阵和旋转方阵，流体力具有相同的形式，只需将上式流体力除以2cos(α)即可。

本发明的有益效果是：具有显式表达、输入试验参数少、不依赖数值积分的优点。

附图说明

图1为四种典型的管阵排列形式；

图2为流管模型；

图3为根据以上横向流弹失稳的原理编制计算程序和聂清德等人的试验参数绘制出正方形管阵(p/d＝1.25)的稳定区图；

图4为根据以上横向流弹失稳的原理编制计算程序和聂清德等人的试验参数绘制出旋转三角形管阵(p/d＝1.25)的稳定区图；

图5为根据以上横向流弹失稳的原理编制计算程序和聂清德等人的试验参数绘制出旋转方正方形管阵(p/d＝1.25)的稳定区图；

图6为根据以上横向流弹失稳的原理编制计算程序和聂清德等人的试验参数绘制出三角形管阵(p/d＝1.25)的稳定区图。

具体实施方式

为了使本发明的目的，技术方案及优点更加清晰，下面以聂清德等人风洞实验为实施例结合附图对本发明作进一步的描述。

本发明提供的一种蒸汽发生器管束两相横向流弹失稳分析方法，其中主要包括以下步骤：1确定管子的基本参数，见表1、表2。

表1 各种管束的试验数据

注：m＝1.24kg/m管径d＝0.022m

表2 实验管束的几何参数

2确定两相横向流弹失稳参数确定

2.1两相横向流弹失稳试验参数确定

唯一依赖于试验输入的参数为开始产生附着和分离点的位置，这可以通过Weaver和AbdRabbo以及Scott的可视化研究进行估计。四种典型的管阵排列的可以归为以下两组：

注：β为流体附着表面的法向矢量角度

2.2两相横向流弹失稳计算参数确定

2.2.1流场计算参数

(1)空泡份额

根据均相流模型(HEM)假设，两相流的两相充分混合，并且密度和温度分布均匀，不存在气相和液相的相对滑移，可将空泡份额定义为：

${\mathrm{\&alpha;}}_H=1/(1+\frac{1-X_a}{X_a}\frac{{\mathrm{\&rho;}}_g}{{\mathrm{\&rho;}}_l})$

式中，X_a为真实含汽率，ρ_g和ρ_l分别为气相和液相的密度。

(2)流动参数

节距流速Up：U_p＝U_∞P/(P-D)

均相流密度ρ：ρ＝ρ_l(1-α_H)+ρ_gα_H

自由来流速度U_∞： $U_{\&infin;}=\frac{{\mathrm{\&rho;}}_l{V}_l+{\mathrm{\&rho;}}_g{V}_g}{\mathrm{\&rho;A}}$

自由来流质量流速：

节距质量流速： ${\stackrel{\&CenterDot;}{m}}_p={\stackrel{\&CenterDot;}{m}}_{\&infin;}\frac{P}{(P-D)}=\mathrm{\&rho;}{U}_p$

式中，U_∞是自由来流速度(即假设不存在管子时的速度)，P是管中心距，D是管子直径(对于翅片管，用等效水力直径D_h代替)。A为自由来流的流道截面积；ρ_g为气相密度、ρ_l为气相密度；Vg为气相速度、V_l为气相速度。

2.2.2动水参数

(1)水力直径

对于非翅片管，水力直径就简单的为管子外径D。对于翅片管，水力直径等效为D_h：

D_h＝D_r+R_F(D_o-D_r)

式中，D_r为齿根圆直径，D_o为翅片外径，比值R_F为翅片占据的体积与管子外表面到翅片外径表面的有效体积的比。对于翅片管，D_h要用于振动激励参数的计算，等效节径比P/D也要变为P/D_h。

(2)水动力附加质量

管束中一根管子单位长度的水动力附加质量可以表示为：

$m_h=\left(\frac{\mathrm{\&pi;}}{4}{\mathrm{\&rho;}}_{\mathrm{TP}}{D}^2\right)\left[\frac{{(D_e/D)}^2+1}{{(D_e/D)}^2-1}\right]$

式中，D_e是周围管子的等效直径，比值D_e/D是管束对管子约束程度的度量。

对于三角形管阵：

D_e/D＝(0.96+0.5P/D)P/D

对于正方形管阵：

D_e/D＝(1.07+0.56P/D)P/D

单位长度管子总的动力质量m包括水动力附加质量m_h、管子本身的质量m_t和管内流体的质量m_i：

m＝m_h+m_t+m_i

(3)阻尼

两相流中多跨热交换器管的阻尼包括支撑阻尼ζ_S、粘性阻尼ζ_V和两相阻尼ζ_TP：

ζ_T＝ζ_S+ζ_V+ζ_TP

通量很高、空泡份额很大的情况下，支撑就可能是干的，此时只有摩擦阻尼存在，这种情况下，支撑阻尼与气体中的热交换器管阻尼基本相同：

ζ_T＝ζ_S+ζ_V+ζ_TP

管和支撑之间存在液体时，支撑阻尼ζ_S既包括压膜阻尼ζ_SF又包括摩擦阻尼ζ_F：

${\mathrm{\&zeta;}}_S={\mathrm{\&zeta;}}_{\mathrm{SF}}+{\mathrm{\&zeta;}}_F=\left(\frac{N-1}{N}\right)[\frac{1460}{f}\left(\frac{{\mathrm{\&rho;}}_l{D}^2}{m}\right)+0.5]{\left(\frac{L}{l_m}\right)}^{1/2}$

粘性阻尼：

${\mathrm{\&zeta;}}_V=\frac{100\mathrm{\&pi;}}{\sqrt{8}}\left(\frac{{\mathrm{\&rho;}}_{\mathrm{TP}}{D}^2}{m}\right){\left(\frac{2{\mathrm{\&upsi;}}_{\mathrm{TP}}}{\mathrm{\&pi;}{\mathrm{fD}}^2}\right)}^{1/2}\left\{\frac{[1+{(D/D_e)}^3]}{[1-{(D/D_e)}^2{]}^2}\right\}$

两相阻尼比：

${\mathrm{\&zeta;}}_{\mathrm{TP}}=4.0\left(\frac{{\mathrm{\&rho;}}_l{D}^2}{m}\right)\left[f\left({\mathrm{\&alpha;}}_H\right)\right]\left\{\frac{[1+{(D/D_e)}^3]}{{[1-{(D/D_e)}^2]}^2}\right\}$

式中，ρ_l是管和支撑间液体的密度，D_e是周围管子的等效直径，D是管子的外径，m是单位长度管子的总质量，N是折流板的数目，包括用两相流均相密度计算的水动力附加质量，l_m是管子的特征长度，L是支承板宽度，f是管子的频率，ρ_TP是气液两相等效密度，υ_TP是等效运动粘度，υ_TP＝υ_l/[1+α_H(υ_l/υ_g-1)]，υ_g和υ_l分别为液相和气相的运动黏度，而f(α_H)的确定

f(α_H)＝α_H/40         α_H<40％

f(α_H)＝1              40％≤α_H≤70％。

f(α_H)＝1-(α_H-70)/30  α_H>70％

3建立并求解两相横向流弹失稳的解析模型

聂清德等人风洞实验：风洞宽0.25m，高0.35m。试验段中的自由流速可达40m/s，其上游的紊流度接近1％。放在试验段中的管束分别按正方形、转角三角形、转角正方形、三角形四种典型的方式布置，节径比P/d为1.25与1.41。管子材料为T2紫铜，管长0.39m，直径0.022m。上端以螺纹连接的方法与上管板固定后呈自由悬挂状，使管子的一阶固有频率为90Hz。在距离管子的固定端20mm处，沿圆周且相应于流向以及与流动垂直的方向各粘贴两片箔式应变片，采用半桥双臂互补法联接到动态应变仪的电桥盒上。管子振动时，应变片电阻变化反映了管子自由端位移的大小，由此得到流速与管子振幅的关系曲线。相应于振幅陡然急剧增大时的流速，即为发生流体弹性不稳定时的临界流速V_C，实验数据参见表1、表2。

在Liming提出的单相横向流弹失稳显式分析模型-“流管”模型的基础上，提出针对两相横向流弹失稳的解析模型。热交换器通常有四种典型的管阵几何排列形式，对于所有的几何排列形式，解析模型的基本概念是相同的，排列形式只是影响模型中使用的试验参数。对于平行三角形管束模型：

3.1不稳定流速

式中Y，很多研究已经证明流体弹性不稳定通常首先在升力方向(y方向)产生。

假设柔性管在两自由度方向以相同的固有频率做简谐运动

$m\stackrel{\&CenterDot;\&CenterDot;}{x}+c\stackrel{\&CenterDot;}{x}+\mathrm{kx}=F_x{a}_1(s_m,t)+F_x{a}_2(s_m,t)$

$m\stackrel{\&CenterDot;\&CenterDot;}{y}+c\stackrel{\&CenterDot;}{y}+\mathrm{ky}=F_y{a}_1(s_m,t)-F_y{a}_2(s_m,t)$

其中，a₁(s_m,t)和a₂(s_m,t)分别是流管1和流管2在最小间隙处的面积微扰(面积变小为正，变大为负。)，它们由管子几何参数定义。

$a_1(s_m,t)=-\frac{1}{2}y\left(t\right),a_2(s_m,t)=\frac{1}{2}y\left(t\right)$

方程的解为y(t)＝Ye^iωt

其中ω是振动复频率。U_r＝U₀/ω_nl₀是缩减速度，l₀＝2s_l，频率比ω_r＝ω/ω_n，ω_n是管子在静水中的固有频率，s_l是从振动管到压力扰动可以忽略的位置的距离。

由于扰动衰减和时间衰减仅发生在“流管”路径上，以下给出两组针对平行三角形管束的面积扰动方程

$f\left(s\right)=\frac{s-s_1}{s_a-s_1},s_1\&le;s\&le;s_a$

f(s)＝1,s_a≤s≤s_s

而管子运动和流体力之间的时间延迟通过相位函数来表示

为了便于使用，在接下来的公式中将写为

$\stackrel{\&OverBar;}{s}=\frac{s-s_a}{s_1-s_a},$ 则衰减函数 $f\left(s\right)=1-\stackrel{\&OverBar;}{s}$

$U_r=\frac{U_p}{\mathrm{fD}}\frac{1}{l_0{A}_0}\frac{\cos \mathrm{\&alpha;}}{2\mathrm{\&pi;}}(\mathrm{Pr}-1)$

其中：

$A_4=-\frac{{\mathrm{\&omega;}}_{r}i}{l_0{A}_0{U}_r}(s_l-s_a)$

3.2非稳态压力的计算

各分项计算：

B₂＝(s₁-s_a)A₄/2

B₃＝(s₁-s_a)A₅

将以上代入非稳态压力的积分计算表达式得：

其中：

C₁＝-ω_rl₀U_riB₂

$C_2=-{\mathrm{\&omega;}}_r{l}_0{U}_{r}i{B}_3-2l_0^2{U}_r^2{A}_4-\frac{{\mathrm{\&omega;}}_r{l}_0{U}_{r}i}{A_0}(s_1-s_a)$

$\begin{array}{c}{C}_3=-{\mathrm{\&omega;}}_r{l}_0{U}_{r}i{B}_1-2l_0^2{U}_r^2{A}_5-\frac{l_0^2{U}_r^2}{A_0}\\ -\frac{{\mathrm{\&omega;}}_r{l}_0{U}_r}{A_0}{B}_4\end{array}$

3.3作用在管子上的流体力

$\begin{array}{c}{F}_y=\frac{Y}{D}{e}^{\mathrm{i\&omega;t}}\cos \mathrm{\&alpha;}\underset{s_a}{\overset{s_s}{\&Integral;}}{p}_2(s,U_r)\cos \left(\frac{2s}{\mathrm{Pr}}\right)\mathrm{ds}\\ =\frac{Y}{D}{e}^{\mathrm{i\&omega;t}}\cos \mathrm{\&alpha;}\underset{s_a}{\overset{s_s}{\&Integral;}}(C_1{\stackrel{\&OverBar;}{s}}^2+C_2\stackrel{\&OverBar;}{s}+C_3)\cos \left(\frac{2s}{\mathrm{Pr}}\right)\mathrm{ds}\end{array}$

计算各分项：

为便于计算令

将以上各项代入流体力的计算公式可得：

$\begin{array}{c}{F}_y=\frac{Y}{D}{e}^{\mathrm{i\&omega;t}}\cos \mathrm{\&alpha;}\underset{s_a}{\overset{s_s}{\&Integral;}}{p}_2(s,U_r)\cos \left(\frac{2s}{\mathrm{Pr}}\right)\mathrm{ds}\\ =\frac{Y}{D}{e}^{\mathrm{i\&omega;t}}\cos \mathrm{\&alpha;}\underset{s_a}{\overset{s_s}{\&Integral;}}(C_1{\stackrel{\&OverBar;}{s}}^2+C_2\stackrel{\&OverBar;}{s}+C_3)\cos \left(\frac{2s}{\mathrm{Pr}}\right)\mathrm{ds}\\ =\frac{Y}{D}{e}^{\mathrm{i\&omega;t}}\cos \mathrm{\&alpha;}\underset{i=1}{\overset{3}{\mathrm{\&Sigma;}}}{C}_i{D}_i\end{array}$

F_y＝F_R+F_Ii是复数形式表示的流体力，将上式代入公式

$\frac{m}{\mathrm{\&rho;}{D}^2}=\frac{F_y(U_r,\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_n})}{1-{\left(\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_n}\right)}^2+i\frac{\mathrm{\&delta;}}{\mathrm{\&pi;}}\left(\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_n}\right)}$

$F_y={2l}_0{U}_r\cos \mathrm{\&alpha;}{F}_y\left(s\right)=F_R(U_r,\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_n})+i{F}_l(U_r,\frac{\mathrm{\&omega;}}{{\mathrm{\&omega;}}_n})$

中就可得到折合流速和质量阻尼比之间的关系，及可得到的折合流速与质量阻尼比的关系曲线得到管束的稳定区图。对于一个确定的流管模型，质量阻尼比是确定，只需要计算出具体折合流速，然后判断两者确定的点，若在稳定区图内，则说明蒸汽发生器管束不会发生流弹失稳，若两者确定的点在稳定区图外，则说明蒸汽发生器管束将发生流弹失稳。值得注意的是以上推导过程是根据平行三角形管阵推导而来，对于旋转三角形管阵和旋转方阵，流体力具有相同的形式，只需将上式流体力除以2cos(α)即可。

根据以上横向流弹失稳的原理编制计算程序，并根据聂清德等人的试验参数绘制出各种管束的稳定区图。由图可见，对常见的四种管束，显式解析模型方法的结果与试验结果基本一致。当mδ/ρD²>5，基于当前方法的结果相比试验结果更保守。当mδ/ρD²<5，计算结果相对试验结果偏大，但是两者之间的误差并不是很显著。结果表明，该显式解析模型继承了Y-W“流管”模型的优点，即具有较好的计算精度及需要较少的试验数据。

Claims (5)

1.一种蒸汽发生器管束两相横向流弹失稳分析方法，其特征在于，蒸汽发生器管束两相横向流弹失稳分析方法包括：建立两相横向流弹失稳的显式解析模型；确定求解该解析模型所需要的两相流计算参数；最后求解解析模型得到管束的稳定区图。

2.如权利要求1所述的一种蒸汽发生器管束两相横向流弹失稳分析方法，其特征在于，建立两相横向流弹失稳的显式解析模型中，为了考虑从中心管向外的面积微扰，采用以下线性衰减函数

s₁≤s≤s_a

f(s)＝1,s_a≤s≤s_s

式中，曲线坐标s用以表示流动路径；足够远的上游位置为s_l；流管在s_a位置开始附着于柔性管，并在s_s位置分离。

3.如权利要求1所述的一种蒸汽发生器管束两相横向流弹失稳分析方法，其特征在于，建立两相横向流弹失稳的显式解析模型中，扰动衰减和时间衰减仅发生在“流管”路径上，以下给出两组针对平行三角形管束的面积扰动方程

  s₁≤s≤s_a

a₂＝cosα,f(s)＝1.0,  s_a≤s≤s_s

式中，f(s)为线性衰减函数，管子运动和流体力之间的时间延迟通过相位函数，为扰动衰减和时间衰减的相位函数，为扰动衰减的相位函数。

4.如权利要求1所述的一种蒸汽发生器管束两相横向流弹失稳分析方法，其特征在于，确定求解该解析模型所需要的两相流计算参数，根据均相流模型(HEM)假设，两相流的两相充分混合，并且密度和温度分布均匀，不存在气相和液相的相对滑移，可将空泡份额定义为：

式中，X_a为真实含气率，ρ_g和ρ_l分别为气相和液相的密度。

5.如权利要求1所述的一种蒸汽发生器管束两相横向流弹失稳分析方法，其特征在于， 确定求解该解析模型所需要的两相流计算参数，两相流中多跨热交换器管的阻尼包括支撑阻尼ζ_S、粘性阻尼ζ_V和两相阻尼ζ_TP：

ζ_T＝ζ_S+ζ_V+ζ_TP

通量很高、空泡份额很大的情况下，支撑就可能是干的，此时只有摩擦阻尼存在，这种情况下，支撑阻尼与气体中的热交换器管阻尼基本相同：

ζ_T＝ζ_S+ζ_V+ζ_TP

管和支撑之间存在液体时，支撑阻尼ζ_S既包括压膜阻尼ζ_SF又包括摩擦阻尼ζ_F：

粘性阻尼：

两相阻尼比：

式中，ρ_l是管和支撑间液体的密度，D_e是周围管子的等效直径，D是管子的外径，m是单位长度管子的总质量，N是折流板的数目，包括用两相流均相密度计算的水动力附加质量，l_m是管子的特征长度，L是支承板宽度，f是管子的频率，ρ_TP是气液两相等效密度，υ_TP是等效运动粘度，υ_TP＝υ_l/[1+α_H(υ_l/υ_g-1)]，υ_g和υ_l分别为液相和气相的运动黏度，而f(α_H)的确定

f(α_H)＝α_H/40  α_H<40％

f(α_H)＝1  40％≤α_H≤70％。

f(α_H)＝1-(α_H-70)/30  α_H>70％。