基于collocation方法的Taylor-Culick模型流動穩定性的特征向量分析①

李 陽，劉佩進，金秉寧

(西北工業大學 燃燒、熱結構與內流場重點實驗室，西安 710072)

0 引言

固體火箭發動機(如美國的Space Shuttle SRM、RSRM、Titan系列以及歐洲的Arian 5 P230)工作過程中容易出現軸向、低階的壓強振蕩和推力振蕩[1]。由于發動機幾何尺寸較大，工作時間又很長，在實際預估中，對每一時刻的工作狀態進行數值模擬是不現實的。因此，在流動穩定性機理研究的基礎上，譜方法作為一種快速預估手段，可用來分析流動不穩定特征以及演化規律。

國外自20世紀60～70年代開始，利用穩定性理論分析固體火箭發動機的流動穩定性。Varapaev和Yagodkin最先在平行流的假設下，研究了Taylor平面流的流動穩定性，導出了流函數表示的類似Orr-Sommerfeld方程的擾動方程，得到中性曲線[2〗。Casalis等對Taylor平面流的穩定性做了更加深入的研究，但結果與Varapaev和Yagodkin的結果有一些差異[3]。后來，Griffond發現了上述差異產生的原因是由于所用變量不同所致[4]。由主變量和流函數所表示的擾動方程不同，但所得結果差異不大。Akiki研究可壓縮Taylor平面流的穩定性時出現了幅值的偏差，認為最主要的原因可能是解析方法把有旋量和無旋量分開去求解，而數值方法是整體求解[5]。但Chedevergne發現，解析渦聲解和流動穩定性模態的疊加可準確地重現DNS結果。因此，Akiki的解釋并不令人信服[6]。除了Taylor流模型，Griffond[4]利用Taylor-Culick模型，使用了打靶法，得到了其兩個不穩定頻率，并分析了中性曲線和增長因子。國內對于固體火箭發動機內部流動不穩定性的研究較少，楊尚榮應用局部非平行理論分析了關于主變量和流函數導出的擾動方程的差異，比較了基于空間發展模式的理論結果與大渦模擬計算結果[7]。穩定性預估方法在Taylor平面流的研究中已獲得一定進展，但平面模型顯然與實際發動機差別較大，更合理的模型是Taylor-Culick流。雖然Taylor-Culick流更接近實際發動機，但不可避免的將引入柱坐標，其對稱軸處在數值計算時產生奇性，影響計算結果的準確性。

本文利用譜配置法對Taylor-Culick流動模型進行穩定性分析時，利用配置微分矩陣的方式處理了對稱軸處奇性的問題，與Griffond[4]利用打靶法相比，不需要選擇初值，減少了人為因素的加入，適應性更強，且計算更為簡便。國內外對于固體火箭發動機內部流動不穩定現象的研究主要關注于流動不穩定發生的可能性及頻率，但對于流動不穩定所產生的影響及不穩定的局部特征沒有相應的分析。對于不穩定特征研究，可更有效地揭示流動不穩定現象，并啟發對于流動不穩定的抑制。因此，本文還分析了不同不穩定頻率下的特征向量，用來反映流動不穩定的局部特征。

1 Taylor-Culick流方程處理

1.1 物理及數學模型

本文研究的徑向加質Taylor-Culick流幾何構型與坐標見圖1。

采用不可壓粘性流體N-S方程，利用歸一化參考量：半徑R、徑向加質速度Vinj、密度ρ和運動粘性系數ν進行無量綱化。并將瞬時變量分解為平均量和擾動量的和，得到擾動方程：

·u′=0

(1)

平均量作為基本流，可事先求出解析解[8-9]：

(2)

利用分離變量法，假設擾動量為簡正模態形式：

=(ur,uθ,uz,p)(r)ei(mθ+αz-ωt)

(3)

圖1 幾何模型

進行空間穩定性分析，ω為實數，代表無量綱時間擾動頻率，α=αr-iαi為復數，實部αr為軸向擾動頻率，負虛部αi為軸向擾動局部增長率。m為正整數，代表環向波數。將上式代入擾動方程中，得到如下關于α的多項式特征值問題：

(4)

矩陣L的元素Lij分別為

Taylor-Culick流的邊界條件為壁面上的速度邊界條件：頭部壁面速度為零、軸向坐標軸上的速度對稱條件、側壁上注入速度為常數以及側壁上的無滑移邊界條件。利用擾動量表示物理邊界條件得到擾動方程的齊次邊界條件,利用擾動方程還可得到壓力邊界條件。

(ur,uθ,uz)(±1)=0;Dp(±1)=0

(5)

確定時間擾動頻率ω和環向波數m后，便可利用譜配置方法求解特征值問題，得到任意軸向位置z處的頻率ω對應的增長率αi和頻率αr。

1.2 微分矩陣的建立

由于配點法所假設的近似解是在定義域內關于配點對未知函數進行的高階多項式插值逼近，所以可對該插值多項式進行微分，得到微分矩陣：

(6)

完成對未知函數導數項的逼近，進一步建立微分算子矩陣的具體形式(式(6))[10]。

高階微分矩陣由下面公式計算得到，l階微分矩陣為一階微分矩陣的l次方。

D(l)=(D(1))l,l=1,2,…

(7)

由于計算采用柱坐標，穩定性方程會出現含有1/r和1/r2的項，這些項的值在軸線r=0處無窮大，對稱軸處出現奇性，需要在數值計算時特殊對待。本文考慮到在計算矩陣算子時，邊界r=0在矩陣中始終表現為塊矩陣的一行，而矩陣之間只有加減運算而不做矩陣乘法。也就是說，奇性始終作用在邊界點上，不會進入到內部節點的計算中。

矩陣算子中邊界點上的值最后會被r=0處的邊界條件替換掉。因此，本文直接令r=1參與r=0處的計算，然后用邊界條件替換掉矩陣算子中r=0所在的行向量。這樣既可提高精度，又可解決對稱軸處奇性的問題。圖2陰影所示為L矩陣配置位置，分別為分塊矩陣r的倒數項對應于r=0所在的行。

1.3 計算結果驗證

為驗證算法及程序，計算了長徑比z=10，雷諾數Re=4500，時間波動頻率ω=80，環向波數m=0下的不穩定模態，得到的特征譜如圖3所示(圖中坐標無量綱)，存在兩個不穩定模態。與Griffond的結果[4]比較，結果如表1所示。

圖2 L矩陣配置位置

模態Griffond的結果[4]本文結果誤差/‰1αr6.095 294 565 66.095 295 97αi-1.078 799 814 0-1.078 801 377<12αr3.326 428 536 63.326 428 826αi-0.109 552 558 9-0.109 552 685 1<1

圖3 特征譜

本文所得結果與文獻結果基本吻合，可驗證本文所用方法可行。誤差在1‰以內，滿足精度要求。

2 結果分析

對于3個不同頻率ω，特征值α=0對應的無量綱特征向量見圖4。隨著時間波動頻率的增大，徑向的振蕩周期變小，振動幅度減小。這說明時間波動頻率與徑向波動頻率是同步的，振動頻率越小，振動幅度越大，使得振蕩滲透向軸心(r=0)處的深度也更大。頻率越低的不穩定振動，影響范圍越大。

圖5比較了3種雷諾數下，特征值α=0對應的無量綱特征向量的變化。隨著雷諾數的增加，振動幅度變化不大，但振蕩波向內部的滲透深度變大，振幅的衰減率減小。徑向加質速度越大，雷諾數越大，振動可在幅度不變的情況下影響范圍越大。流動不穩定現象所表現出來的是流體振動的失穩，其特點與一般振動所表現的效果相似。從特征向量可看出，流動失穩首先從壁面開始，逐漸向中心滲透，頻率與振幅呈反比例增長。

(a) ω=30 (b) ω=50 (c) ω=80

(a) Re=900 (b) Re=2000 (c) Re=4000

可從能量的角度進行分析: 當總能量保持不變而頻率變化的情況下，高頻振蕩從壁面注入后很快衰減，低頻振蕩衰減較慢，能夠更深入地滲透到流動區域內；雷諾數的變化同時反映了能量的變化，振蕩從能量注入的位置開始，向內部滲透，能量越大，滲透的深度越深。流動穩定性的特征向量變化表現出阻尼振蕩的特征。整體的能量注入、能量傳遞、能量耗散和能量隨流動流出控制區域的形式和變化，還有待后續工作進行闡述。

3 結論

(1)對于對稱軸處的奇性，可利用配置矩陣的方式進行處理，相比打靶法過分依賴于由經驗所給定的初值，其實用性更強，更適用于計算機計算。

(2)分析了Taylor-Culick不同條件下的特征向量，得到了其流動不穩定局部特征。流動不穩定由能量注入處發生，雷諾數越大且頻率越小，滲透部分越深。這意味著大雷諾數的低頻振蕩輸入對應的模態振型對內部的作用更強烈，影響流場的范圍更大。由于非線性理論可分析有限幅值振蕩，建議在后續工作中考慮非線性因素的影響。