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%。由此現象分析，計算中，對于周向速度的處理方法將直接影響動特性系數的計算結果及其隨轉速的變化趨勢，本書提出的分析方法在原流體微元控制方程組中加入周向動量方程，對流場的描述更加準確，一定程度上提高了計算精度，但方法中涉及的周向速度的簡化方式仍有優化的空間。