基于collocation方法的Taylor-Culick模型流動(dòng)穩(wěn)定性的特征向量分析①

李 陽(yáng)，劉佩進(jìn)，金秉寧

(西北工業(yè)大學(xué) 燃燒、熱結(jié)構(gòu)與內(nèi)流場(chǎng)重點(diǎn)實(shí)驗(yàn)室，西安 710072)

0 引言

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

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

本文利用譜配置法對(duì)Taylor-Culick流動(dòng)模型進(jìn)行穩(wěn)定性分析時(shí)，利用配置微分矩陣的方式處理了對(duì)稱(chēng)軸處奇性的問(wèn)題，與Griffond[4]利用打靶法相比，不需要選擇初值，減少了人為因素的加入，適應(yīng)性更強(qiáng)，且計(jì)算更為簡(jiǎn)便。國(guó)內(nèi)外對(duì)于固體火箭發(fā)動(dòng)機(jī)內(nèi)部流動(dòng)不穩(wěn)定現(xiàn)象的研究主要關(guān)注于流動(dòng)不穩(wěn)定發(fā)生的可能性及頻率，但對(duì)于流動(dòng)不穩(wěn)定所產(chǎn)生的影響及不穩(wěn)定的局部特征沒(méi)有相應(yīng)的分析。對(duì)于不穩(wěn)定特征研究，可更有效地揭示流動(dòng)不穩(wěn)定現(xiàn)象，并啟發(fā)對(duì)于流動(dòng)不穩(wěn)定的抑制。因此，本文還分析了不同不穩(wěn)定頻率下的特征向量，用來(lái)反映流動(dòng)不穩(wěn)定的局部特征。

1 Taylor-Culick流方程處理

1.1 物理及數(shù)學(xué)模型

本文研究的徑向加質(zhì)Taylor-Culick流幾何構(gòu)型與坐標(biāo)見(jiàn)圖1。

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

·u′=0

(1)

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

(2)

利用分離變量法，假設(shè)擾動(dòng)量為簡(jiǎn)正模態(tài)形式：

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

(3)

圖1 幾何模型

進(jìn)行空間穩(wěn)定性分析，ω為實(shí)數(shù)，代表無(wú)量綱時(shí)間擾動(dòng)頻率，α=αr-iαi為復(fù)數(shù)，實(shí)部αr為軸向擾動(dòng)頻率，負(fù)虛部αi為軸向擾動(dòng)局部增長(zhǎng)率。m為正整數(shù)，代表環(huán)向波數(shù)。將上式代入擾動(dòng)方程中，得到如下關(guān)于α的多項(xiàng)式特征值問(wèn)題：

(4)

矩陣L的元素Lij分別為

Taylor-Culick流的邊界條件為壁面上的速度邊界條件：頭部壁面速度為零、軸向坐標(biāo)軸上的速度對(duì)稱(chēng)條件、側(cè)壁上注入速度為常數(shù)以及側(cè)壁上的無(wú)滑移邊界條件。利用擾動(dòng)量表示物理邊界條件得到擾動(dòng)方程的齊次邊界條件,利用擾動(dòng)方程還可得到壓力邊界條件。

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

(5)

確定時(shí)間擾動(dòng)頻率ω和環(huán)向波數(shù)m后，便可利用譜配置方法求解特征值問(wèn)題，得到任意軸向位置z處的頻率ω對(duì)應(yīng)的增長(zhǎng)率αi和頻率αr。

1.2 微分矩陣的建立

由于配點(diǎn)法所假設(shè)的近似解是在定義域內(nèi)關(guān)于配點(diǎn)對(duì)未知函數(shù)進(jìn)行的高階多項(xiàng)式插值逼近，所以可對(duì)該插值多項(xiàng)式進(jìn)行微分，得到微分矩陣：

(6)

完成對(duì)未知函數(shù)導(dǎo)數(shù)項(xiàng)的逼近，進(jìn)一步建立微分算子矩陣的具體形式(式(6))[10]。

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

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

(7)

由于計(jì)算采用柱坐標(biāo)，穩(wěn)定性方程會(huì)出現(xiàn)含有1/r和1/r2的項(xiàng)，這些項(xiàng)的值在軸線r=0處無(wú)窮大，對(duì)稱(chēng)軸處出現(xiàn)奇性，需要在數(shù)值計(jì)算時(shí)特殊對(duì)待。本文考慮到在計(jì)算矩陣算子時(shí)，邊界r=0在矩陣中始終表現(xiàn)為塊矩陣的一行，而矩陣之間只有加減運(yùn)算而不做矩陣乘法。也就是說(shuō)，奇性始終作用在邊界點(diǎn)上，不會(huì)進(jìn)入到內(nèi)部節(jié)點(diǎn)的計(jì)算中。

矩陣算子中邊界點(diǎn)上的值最后會(huì)被r=0處的邊界條件替換掉。因此，本文直接令r=1參與r=0處的計(jì)算，然后用邊界條件替換掉矩陣算子中r=0所在的行向量。這樣既可提高精度，又可解決對(duì)稱(chēng)軸處奇性的問(wèn)題。圖2陰影所示為L(zhǎng)矩陣配置位置，分別為分塊矩陣r的倒數(shù)項(xiàng)對(duì)應(yīng)于r=0所在的行。

1.3 計(jì)算結(jié)果驗(yàn)證

為驗(yàn)證算法及程序，計(jì)算了長(zhǎng)徑比z=10，雷諾數(shù)Re=4500，時(shí)間波動(dòng)頻率ω=80，環(huán)向波數(shù)m=0下的不穩(wěn)定模態(tài)，得到的特征譜如圖3所示(圖中坐標(biāo)無(wú)量綱)，存在兩個(gè)不穩(wěn)定模態(tài)。與Griffond的結(jié)果[4]比較，結(jié)果如表1所示。

圖2 L矩陣配置位置

模態(tài)Griffond的結(jié)果[4]本文結(jié)果誤差/‰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 特征譜

本文所得結(jié)果與文獻(xiàn)結(jié)果基本吻合，可驗(yàn)證本文所用方法可行。誤差在1‰以?xún)?nèi)，滿(mǎn)足精度要求。

2 結(jié)果分析

對(duì)于3個(gè)不同頻率ω，特征值α=0對(duì)應(yīng)的無(wú)量綱特征向量見(jiàn)圖4。隨著時(shí)間波動(dòng)頻率的增大，徑向的振蕩周期變小，振動(dòng)幅度減小。這說(shuō)明時(shí)間波動(dòng)頻率與徑向波動(dòng)頻率是同步的，振動(dòng)頻率越小，振動(dòng)幅度越大，使得振蕩滲透向軸心(r=0)處的深度也更大。頻率越低的不穩(wěn)定振動(dòng)，影響范圍越大。

圖5比較了3種雷諾數(shù)下，特征值α=0對(duì)應(yīng)的無(wú)量綱特征向量的變化。隨著雷諾數(shù)的增加，振動(dòng)幅度變化不大，但振蕩波向內(nèi)部的滲透深度變大，振幅的衰減率減小。徑向加質(zhì)速度越大，雷諾數(shù)越大，振動(dòng)可在幅度不變的情況下影響范圍越大。流動(dòng)不穩(wěn)定現(xiàn)象所表現(xiàn)出來(lái)的是流體振動(dòng)的失穩(wěn)，其特點(diǎn)與一般振動(dòng)所表現(xiàn)的效果相似。從特征向量可看出，流動(dòng)失穩(wěn)首先從壁面開(kāi)始，逐漸向中心滲透，頻率與振幅呈反比例增長(zhǎng)。

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

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

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

3 結(jié)論

(1)對(duì)于對(duì)稱(chēng)軸處的奇性，可利用配置矩陣的方式進(jìn)行處理，相比打靶法過(guò)分依賴(lài)于由經(jīng)驗(yàn)所給定的初值，其實(shí)用性更強(qiáng)，更適用于計(jì)算機(jī)計(jì)算。

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