2.3.2 基于摄动法的动力学特性求解

2.3.2.1 控制方程组的建立

如图2-7所示，对螺旋槽间隙内流体微元进行流动及受力分析，得间隙内流体的流动控制方程组如下^[88]：

连续性方程：

周向动量方程：

轴向动量方程：

以上控制方程组为摄动法求解各部分流体力及动力学特性系数的基本方程组，对于光滑环形流域、矩形槽流域及螺旋槽流域等，式（2-46）中的切向力可表示为流动速度与摩擦因数的函数：

其中摩擦因数λ由Blasius或Moody摩擦理论模型求得，由于齿顶间隙内流体流动的形态与两平行圆盘之间的流体流动类似^[88]，故在Blasius模型应用时，取经验系数m=0.066，n=-0.25，即Blasius模型下摩擦因数：

Moody模型下摩擦因数：

其中，e_rs为壁面绝对粗糙度，齿顶间隙内流体在z方向的雷诺数为Re_zl=2C_l0u_zl/υ；式（2-46）中齿顶间隙内流体微元所受切向力的轴向分量为

将齿槽射流流域简化为两平行圆盘之间的流动，则与齿顶间隙流域相似，式（2-46）中齿槽内流域流体微元所受切向力的轴向分量为：

其中，摩擦因数λ_f=0.1，Blasius摩擦模型下，Moody摩擦模型下。

将槽内旋流流域简化为管内流动，流域等效水力直径MG如式（2-52）所示。设旋流流域轴向速度为射流流域轴向速度的1/2，且假设流体在齿槽部分的流动能量损失全部来源于槽内射流流域与旋流流域间的摩擦项，则槽内流域流体微元满足式（2-52）^[88]，结合式（2-39）可进一步求得稳态流动状态下，槽内旋流流域稳态周向速度。

其中，Blasius摩擦模型下λ_d=0.0791×（4u_θdMG/ν）^-0.25，Moody摩擦模型下，Re_θd=2C_g0u_θd/υ。

由于螺旋槽间隙内流体微元的轴向速度远大于周向速度，即Re_z＞＞Re_r，则流体微元所受切向力的周向分量可简化为与轴向分量相关的函数^[88]，则式（2-45）中齿顶间隙内流体微元所受切向力周向分量为：

同理，式（2-45）中槽内流体微元所受切向力的周向分量为：

综上所述，由连续性方程、周向动量方程及轴向动量方程组成的原控制方程组可简化为齿顶间隙流域控制方程组及齿槽流域控制方程组，分别如式（2-56）与式（2-57）所示。

齿顶间隙流域控制方程组：

齿槽流域控制方程组：

2.3.2.2 齿顶流域流体微元控制方程的摄动法求解

参考Childs提出的有限长理论分析方法^[74]，分别对齿顶间隙流域、齿槽流域的轴向、周向动量方程及连续性方程组成的控制方程组进行摄动法求解，求取流域内介质压力沿轴向坐标z的分布函数的数值解，并对此数值解进行拟合求取流体力F_x与F_y，结合式（2-2）中小扰动模型下各动力学特性系数的定义，对应求得包括主刚度、交叉刚度、主阻尼、交叉阻尼及主附加质量等5个动力学特性系数。

假设存在某一偏心小量ε，将轴向速度、周向速度、压力分布及环形间隙径向厚度用偏心量ε表示：

将式（2-62）中各参数的零阶摄动表达式代入式（2-56）～式（2-58）组成的控制方程组中，求得齿顶间隙流域控制方程的零阶周向动量方程及轴向动量方程：

忽略零阶轴向速度u_zl0在轴向的变化，即取u_zl0=u_zlm。由螺旋槽齿顶间隙流域的流动状态可知，周向速度零阶摄动值u_θl0在z=0处等于进口周向速度u_θin，且对式（2-63）进行求解，可得齿顶间隙流域周向速度沿轴向的分布函数：

其中，a=-λ_zl/C_l0。

将式（2-65）中各参数的一阶摄动表达式代入式（2-56）至式（2-58）组成的控制方程组中，并忽略摄动量的高阶小量，求得齿顶间隙流域控制方程的一阶连续性方程、周向动量方程及轴向动量方程：

其中，L′=L_l/cosα。以上方程定义了齿顶间隙流域内流体微元的一阶压力摄动量p_l1、周向速度摄动量u_θl1及轴向速度摄动量u_zl1均为轴向坐标z、周向坐标θ及时间t的函数，即以上三个参数变量可表示为p_l1（z，t，θ）、u_θl1（z，t，θ）及u_zl1（z，t，θ）。值得注意的是，式（2-66）、式（2-67）及式（2-68）组成的控制方程组均满足周向连续性条件：

结合以上周向连续性方程的要求，一阶压力摄动量p_l1、周向速度摄动量u_θl1及轴向速度摄动量u_zl1可表示为以下形式：

引入式（2-71）所示复数变量表示以上三组未知数及密封转子的位移一阶摄动量，对控制方程组一阶摄动方程进行换算，即可得由轴向、周向速度及压力一阶摄动量及其对轴向坐标z的偏导组成的运算方程，如式（2-71）、式（2-72）及式（2-73）。

轴向动量方程一阶摄动形式：

周向动量方程一阶摄动形式：

连续性方程一阶摄动形式：

将式（2-75）代入式（2-72）～式（2-74）中，对各式中的周向速度、轴向速度、压力及位移等变量函数进行无量纲化处理，则式（2-72）、式（2-73）及式（2-74）可整合为齿顶间隙流域流体微元的一阶微分方程组，如式（2-76）所示。

式（2-76）可化为轴向速度、周向速度及压力一阶摄动量u_zl1，u_θl1及p_l1关于轴向坐标z的一阶微分方程组：

其中，E₁₁=E₁₃=E₃₃=0，，，，，，。

以上方程组需满足以下边界条件：

1）齿顶间隙流域流体出口端面一阶压力摄动量为0，即=0；

2）齿顶间隙流域流体进口端面一阶周向速度摄动量为0，即=0；

3）由进口压力损失的定义，齿顶间隙流域流体进口端面一阶压力摄动量与一阶轴向速度摄动量间满足关系：。

结合边界条件采用打靶法对式（2-77）所示微分方程组进行求解，介于进口与出口位置一阶压力摄动量沿轴向的分布均与轴向速度有关，故在求解中假设轴向速度为基础变量γ_k，将进口与出口位置压力值用基础变量表示，即：

且设，则。

式（2-77）中对γ_k求偏导数，可得：

结合原控制方程组的边界条件及式（2-78）可知，式（2-79）的数值求解边界条件为M₁（0）=1，M₂（0）=0，M₃（0）=k。设存在一函数F，且F（γ_k）=。求解中，首先给定γ_k一固定初值，即可求得u_zl1，u_θl1及p_l1随z分布函数的数值解，判断此时p_l1（1）的大小，若p_l1（1）小于设定的残差，如10^-6，则表示函数收敛，停止计算。若p_l1（1）大于设定的残差，则利用函数F，采用牛顿法对γ_k的初值进行修正以加速收敛，修正方法：

综上，原控制方程组式（2-66）的求解可化为对初值γ_k的不断改进过程，并验证原边界条件是否满足的迭代求解过程。最终残差满足计算精度要求，迭代结束，将得到环形间隙内流体压力沿轴向坐标z分布函数的数值解，将此数值解沿z坐标进行拟合，可得压力沿轴向坐标z的分布函数的解析解，其形式如下：

对于该研究对象，某一具有I_s头，长度为L的螺旋槽转子迷宫密封，齿顶流域流体作用于转子上面的反作用力可表示为

基于式（2-2）中对光滑环形间隙内流体等效动特性系数的定义，当t=0时，密封转子外壁面及定子内壁面所受径向与周向流体力可表示为5个动特性系数与涡动转速的多项式，如下：

令式（2-82）、式（2-83）与式（2-84）、式（2-85）分别相等，则可求得动力学特性系数组成的两个多项式，值得注意的是，以上周向与轴向力的求解仅针对某一固定涡动频率Ω，为求得五个动力学特性系数，可针对某一固定工作转速n，取涡动频率为0、0.5、1.0、1.5、2.0倍的工作转速，组成5组五元一次方程组，每三组方程求解出一组动特性系数，5组方程共求解10组动特性系数，将求得每5个动特性系数取平均值，即为所求螺旋槽转子迷宫密封齿顶流域的等效动力学特性系数。

2.3.2.3 槽内流域流体微元的控制方程组摄动法求解

槽内流域流体微元控制方程组摄动法求解中，忽略槽内流体周向速度随z坐标的变化，即取周向速度的零阶摄动量u_θg0=u_θgm。则，将式（2-70）中各参数的一阶摄动形式代入式（2-56）、式（2-57）及式（2-58）中，槽内流域流体微元一阶连续性方程、周向动量方程及轴向动量方程的一阶摄动形式可简化为

其中，L″=L_g/cosα。结合式（2-69）、式（2-70）及式（2-71）的参数处理方法，将以上三式化为关于槽内流域压力、周向速度及轴向速度一阶摄动量关于轴向坐标z的微分方程组。如下：

其中：

E₁₁=E₁₃=E₃₃=0

参考齿顶间隙流域流体微元一阶摄动方程的求解方法，对式（2-89）进行求解，即可求得齿槽内流体压力沿轴向坐标z分布函数的数值解，将此数值解沿z坐标进行拟合，可得压力沿轴向坐标z的分布函数的解析解，其形式如下：

对于具有I_s头，长度为L的螺旋槽转子迷宫密封，槽内流域作用于转子上的反作用力可表示为

以上两式结合齿顶间隙流域等效动力学计算方法，可进一步求得齿槽流域等效动力学特性系数，包括主刚度系数K_g、交叉刚度系数k_g、主阻尼系数C_g、交叉阻尼系数c_g及主附加质量系数M_g。

则，整个螺旋槽转子迷宫密封的动力学特性系数如下：

以文献^[157]中的试件几何尺寸、实验工况为计算模型，将泄漏量与动力学参数计算结果与Iwatsubo的分析方法^[88]及实验结果^[157]做对比。计算模型的具体几何参数及操作工况见表2-1。

图2-8及图2-9给出了Iwatsubo实验模型为计算模型的泄漏量计算实验与计算对比情况。由图中可以看出，在对螺旋形转子迷宫密封稳态流场及泄漏量的求解中，本书所述分析方法及Iwatsubo所述分析方法计算结果均略有偏大，但与实验结果基本一致，泄漏量均随转速的增加而减小。本书及Iwatsubo的求解结果均随转速的增加呈线性减小，实验结果随转速的变化呈抛物线函数形态，减小趋势随转速的增加急剧减小。本书修正的稳态流场分析方法的计算误差小于Iwatsubo的方法，在转速低于3500r/min的实验工况下，误差小于6%，在转速4500r/min的实验工况下，计算误差为8.9%，远小于Iwatsubo的计算误差。图2-9分别给出了采用两种分析方法计算出的螺旋槽转子迷宫密封的主刚度系数、交叉刚度系数、主阻尼系数及交叉阻尼系数四个动力学特性系数及实验结果。由图中可以看出，本书提出的分析方法在主刚度系数、交叉刚度系数及交叉阻尼系数的计算中，误差远小于Iwatsubo的求解方法。特别是针对主刚度系数的计算，在6组实验工况下计算误差均小于10%。两种方法对主刚度系数的计算均偏大，交叉刚度系数及主阻尼系数的计算均偏小，且4个动特性系数的实验结果对转速变化的敏感程度均高于两种计算方法的计算结果。对环形间隙内非定常流体激励力径向分量的计算精度的提高较大，以0.588MPa下的500r/min及600r/min为例，计算精度提高约30%。由此现象分析，计算中，对于周向速度的处理方法将直接影响动特性系数的计算结果及其随转速的变化趋势，本书提出的分析方法在原流体微元控制方程组中加入周向动量方程，对流场的描述更加准确，一定程度上提高了计算精度，但方法中涉及的周向速度的简化方式仍有优化的空间。