非定常運動下的激波／邊界層干擾分離特性研究

趙云飛，劉 偉，劉 緒，鄧小剛

（國防科技大學 航天科學與工程學院，湖南 長沙 410073）

0 引 言

在超聲速和高超聲速飛行中，激波與邊界層相互作用、相互影響的現(xiàn)象十分常見。這種干擾不僅產(chǎn)生氣動加熱，而且對機體的控制和發(fā)動機進氣狀態(tài)也有影響。因此，在高速飛行器開發(fā)、設計時，要保持飛行器飛行的高效率、高可靠性就必須深入了解飛行器流場的結構和特性。

最早注意到激波／邊界層干擾現(xiàn)象是在20世紀30年代末，F(xiàn)erri［1］在超聲速風洞試驗中，發(fā)現(xiàn)在機翼后緣附近激波與邊界層干擾后產(chǎn)生了邊界層分離現(xiàn)象。隨后，包括美國、蘇聯(lián)（俄羅斯）、英國、德國等一些最早發(fā)展航空航天的國家相繼開展了對這種干擾現(xiàn)象的研究。幾十年來，人們通過大量的地面風洞試驗、飛行試驗、理論分析及計算研究，對激波與邊界層干擾現(xiàn)象的復雜性及其對飛行穩(wěn)定性及控制的影響已 經(jīng) 有 了 非 常 深 刻 的 認 識 和 理 解［2－9］，Dolling［4］、Knight［5－6］以及 Zheltovodov［7］等的綜述文獻對這些研究工作做了較全面的總結。

要指出的是，目前國內(nèi)外研究主要是集中在飛行器姿態(tài)不變假定下的激波／邊界層干擾特性研究（包括由定姿態(tài)下非定常分離、脈動引起的非定常激波／邊界層干擾）。然而，高速飛行器對于飛行條件的變化是非常敏感的，存在諸如機動動作引起的、大攻角舵效耦合引起的時延非線性、非定常動態(tài)特性問題。另外，諸多隨機干擾因素對飛行器的飛行狀態(tài)有著非常大的影響，對典型的高超聲速飛行器布局而言，長周期模態(tài)是欠阻尼的（或不穩(wěn)定的），短周期模態(tài)是不穩(wěn)定的。這都使得飛行器處于非定常運動狀態(tài)中。目前對于飛行器非定常運動狀態(tài)下的非定常激波／非定常邊界層干擾問題研究，由于受地面實驗技術等因素的限制，缺乏深入了解，這對于確定干擾流場在多大范圍內(nèi)影響空氣動力及氣動加熱環(huán)境，確定干擾域內(nèi)的氣動力載荷及氣動熱載荷，為工程設計提供最嚴重的臨界設計條件還存在不足。因此，在上述研究工作基礎上，進一步開展飛行器非定常運動狀態(tài)下的非定常激波／邊界層干擾特性研究是有重要意義的。

本文以二維壓縮拐角這一經(jīng)典的激波／邊界層干擾問題為研究對象，采用五階精度加權緊致非線性格式（WCNS）和非定?！半p時間步”格式數(shù)值求解非定常雷諾平均N－S方程模擬壓縮拐角在等速抬頭、等速低頭和周期性俯仰振動等不同運動方式下的非定常激波／邊界層干擾現(xiàn)象，分析非定常激波／邊界層干擾引起的分離流區(qū)非定常變化特性，并考察角速度、振幅和頻率等參數(shù)的影響。

1 數(shù)值方法

1.1 流動控制方程和離散方法

流動控制方程為一般曲線坐標系下的非定常雷諾平均N－S方程：

式中Q ＝J－1（ρ，ρu，ρv，ρw，ρe）T是守恒變量，F(xiàn)、G、H和Fv、Gv、Hv分別是無粘通量和粘性通量。湍流模型采用一方程SA模型。

控制方程時間方向采用隱式非定?！半p時間步”格式求解，時間導數(shù)的離散精度為二階。空間項采用五階精度加權緊致非線性格式WCNS［10］離散。該格式采用六階單元中心型差分格式離散無粘導數(shù)項，以ξ方向為例：

1.2 壓縮拐角運動方程

壓縮拐角的運動方式包括等速抬頭、等速低頭和周期性俯仰振動，通過讓模型繞定軸做單自由度旋轉的方式實現(xiàn)。等速抬頭和等速低頭的運動方程為：

式中θ是俯仰角，θm是角振幅，k是無量綱減縮頻率。在本文的坐標系定義中，壓縮拐角低頭時俯仰角θ為正號，抬頭時俯仰角θ為負號，俯仰運動的軸心設在拐角點處。

壓縮拐角運動過程中，網(wǎng)格與物體剛性固聯(lián)，網(wǎng)格運動滿足幾何守恒律（?J－1／?t≡0）。初場采用定常問題的收斂解。壁面速度滿足粘性無滑移條件V＝Vwall，內(nèi)點速度Vi和虛擬邊界速度Vg滿足切向無滑移、法向無穿透的條件Vg＝2Vwall－Vi。

2 壓縮拐角定姿態(tài)風洞實驗模擬

通過模擬Settles等［11－14］的壓縮拐角定姿態(tài)風洞實驗，對所采用的五階WCNS格式進行驗證。表1是實驗條件，包括8°、16°、20°三個壓縮拐角，來流馬赫數(shù)2.85左右，來流邊界層為湍流態(tài)。實驗結果表明：8°壓縮拐角未發(fā)生流動分離，16°壓縮拐角產(chǎn)生了微弱分離，20°壓縮拐角形成了較大范圍的流動分離。表中δ、δ＊和θ是在壓縮拐角上游50.8mm位置處（即表1中的xupstream位置，該處邊界層尚未進入分離區(qū)）測量的邊界層厚度、位移厚度和動量厚度。數(shù)值計算時，為了保證來流邊界層的狀態(tài)與實驗一致，本文采用文獻［15－16］的做法，通過調(diào)整上游平板長度使邊界層發(fā)展到xupstream位置處時動量厚度θ與表1給出的實驗值一致。表2是按照這種方法確定的平板長度，同時給出了各邊界層厚度與實驗［13］和文獻［16］采用SA模型計算結果的對比，整體上符合較好。

表1 壓縮拐角實驗條件和來流邊界層參數(shù)［13］Table1 Freestream and inflow boundary layer data for compression ramp cases［13］

表2 來流邊界層參數(shù)的對比Table2 Inflow matching results

圖1是壓縮拐角計算網(wǎng)格，以20°情況為例，坐標原點位于拐角點處，如前所述拐角前面很長一段平板用來發(fā)展來流邊界層。網(wǎng)格數(shù)量為256×131（流向×法向），物面網(wǎng)格距離為1×10－6m，可保證無量綱壁面距離y＋＜1。

圖1 20°壓縮拐角計算網(wǎng)格Fig.1 Grid of 20°compression corner

圖2是壓縮拐角分離區(qū)陰影圖與實驗［11］的比較。8°壓縮拐角情況下，激波緊貼物面，流動沒有分離；16°壓縮拐角存在微弱分離；20°壓縮拐角形成了較大范圍的流動分離，有明顯的分離激波出現(xiàn)。計算結果與實驗在激波的形狀、位置和流動分離現(xiàn)象等方面與實驗基本吻合。

圖3是壓縮拐角表面壓力（以來流壓力無量綱）與實驗結果［13］的定量對比。8°和16°壓縮拐角情況下計算結果與實驗符合的較好；20°壓縮拐角情況下，分離區(qū)內(nèi)數(shù)值計算的壓力比實驗值偏高（根據(jù)Knight［6］的綜述文獻，這是RANS湍流模型在模擬強干擾、大分離激波／邊界層干擾問題時普遍存在的問題），而在后面恢復段與實驗一致。以上對比驗證情況說明本文采用五階WCNS格式模擬的壓縮拐角激波／邊界層干擾結果整體上是合理的。

圖2 壓縮拐角分離區(qū)陰影圖的比較（上：實驗結果［11］，下：本文計算結果）Fig.2 Comparison of shadowgraghs of compression corner with experiment（up：experiment［11］；down：present）

圖3 壓縮拐角表面壓力分布與實驗結果的對比Fig.3 Surface pressure distribution of the 8°，16°and 20°compression corner

3 動姿態(tài)壓縮拐角計算結果

3.1 計算網(wǎng)格

本節(jié)以20°壓縮拐角為模型研究其在等速抬頭、等速低頭和周期性俯仰振動等運動方式下的非定常激波／邊界層干擾問題，自由來流條件與表1相同，來流馬赫數(shù)2.85。與圖1定姿態(tài)情況的計算網(wǎng)格相比，姿態(tài)運動情況下對壓縮拐角計算網(wǎng)格做了以下調(diào)整（如圖4）：

（1）為了更好模擬壓縮拐角俯仰運動時平板前緣的非定常繞流情況，采用10°楔面補充了平板下方的計算區(qū)域，并在楔前緣采用采用很小的球頭鈍化（球頭直徑與平板長度的比值約為1／10000，測試情況表明這種鈍化處理不影響計算結果的準確性）。

（2）由于壓縮拐角運動時分離區(qū)大小會產(chǎn)生非定常變化，因此在分離區(qū)附近流向方向上進行了更寬范圍的加密。網(wǎng)格總數(shù)為流向503×法向151。

3.2 等速抬頭、等速低頭運動

讓20°壓縮拐角按照方程（3）做等速抬頭和等速低頭運動，無量綱角速度取ω＝±0.001745，正號代表低頭即逆時針旋轉，負號代表抬頭即順時針旋轉，旋轉的軸心位于拐角處。初始時刻模型的俯仰角為零，以定常收斂解作為零時刻的初場。

圖4 姿態(tài)運動情況下的壓縮拐角計算網(wǎng)格Fig.4 Grid of compression corner for a movement condition

圖5和圖6分別是壓縮拐角等速抬頭和等速低頭運動時的流場情況，從上到下包括5個不同的角位置，抬頭：θ＝（0°，－4°，－8°，－12°，－16°），低頭：θ＝（0°，4°，8°，12°，16°）。可以看出壓縮拐角俯仰運動時，分離點、再附點和分離激波等都處于非定常運動中，且分離區(qū)的大小變化具有明顯規(guī)律：等速抬頭運動時分離區(qū)減小，等速低頭運動時分離區(qū)增大。

模型運動引起分離區(qū)大小變化的原因可從兩方面考慮：一是模型運動時姿態(tài)角發(fā)生了變化（引起了流動條件的變化），二是壁面運動速度等非定常效應的影響。為此將壓縮拐角分別固定在θ＝±（0°，4°，8°，12°，16°）這些位置上不動，在同樣的來流條件下做定常計算，進行對比分析。本文分離區(qū)大小的定義方法如圖7，其中Δs代表分離點到拐角的距離，Δr代表再附點到拐角的距離（分離點和再附點根據(jù)表面摩阻系數(shù)Cf＝0確定），以Δs與Δr的和作為分離區(qū)的大小Δ，即Δ＝Δs＋Δr。分離區(qū)大小Δ隨俯仰角θ的變化曲線如圖8和圖9。

圖5 壓縮拐角等速抬頭過程中5個不同角位置上的流場情況，θ＝（0°，－4°，－8°，－12°，－16°）Fig.5 Flow of compression corner at five pitching angles during the pitching up course，θ＝（0°，－4°，－8°，－12°，－16°）

圖6 壓縮拐角等速低頭過程中5個不同角位置上的流場情況，θ＝（0°，4°，8°，12°，16°）Fig.6 Flow of compression corner at five pitching angles during the pitching down course，θ＝（0°，4°，8°，12°，16°）

圖7 壓縮拐角分離區(qū)大小的定義方式Fig.7 The definition of separation size

圖8 靜止不動與抬頭運動的壓縮拐角分離區(qū)大小對比Fig.8 Comparison of separation size between static and pitching up compression corner

圖9 靜止不動與低頭運動的壓縮拐角分離區(qū)大小對比Fig.9 Comparison of separation size between static and pitching down compression corner

對于等速抬頭運動的情況，如圖8，static是將壓縮拐角固定在各角位置上不動時的定常計算結果，這時不存在壁面運動效應的影響，分離區(qū)大小的變化完全是由于姿態(tài)角的不同造成的，俯仰角（絕對值）越大分離區(qū)越小，這與圖5壓縮拐角運動時分離區(qū)大小隨俯仰角的變化趨勢相同，說明壓縮拐角運動時姿態(tài)角的變化是影響分離區(qū)大小變化的主要因素之一。從圖8中也可以看出，與靜止不動的情況相比，壓縮拐角等 速 抬 頭 運 動 時 （ω ＝ －0.000872、－0.001745、－0.003491）分離區(qū)減小的相對更慢，說明非定常遲滯作用的影響減緩了分離區(qū)的變化。隨著角速度減小曲線向靜止不動時的曲線收斂。

對圖9等速低頭的情況進行類似的分析后得到的結論是：壓縮拐角低頭運動時俯仰角的增大使分離區(qū)增大；受非定常遲滯作用的影響，低頭運動時分離區(qū)的增大速度也比靜止不動時慢。

下面考察不同姿態(tài)角對分離區(qū)大小的影響。注意到在圖5（a）和圖6（a）的流場整體馬赫數(shù)云圖中，壓縮拐角抬頭或低頭運動時來流在楔的上表面發(fā)生了膨脹或壓縮，這樣對于拐角區(qū)的流動來說其來流邊界層的馬赫數(shù)相應發(fā)生了變化。因此考慮上述姿態(tài)角不同引起壓縮拐角分離區(qū)大小變化也可能是由于姿態(tài)角變化后引起來流馬赫數(shù)的變化而造成的。為此，在其他來流條件不變的情況下，將模型的俯仰角固定在θ＝0°，只改變來流馬赫數(shù)，單方面考察來流馬赫數(shù)對壓縮拐角分離區(qū)大小的影響。表3是5個來流馬赫數(shù)條件，分別按照壓縮拐角在俯仰角為8°、4°、0°、－4°和－8°時的波后馬赫數(shù)取值，以便參考。圖10是分離區(qū)大小隨馬赫數(shù)的變化關系，當馬赫數(shù)增大時壓縮拐角的分離區(qū)減小。因此可以將壓縮拐角運動時引起來流馬赫數(shù)變化對分離區(qū)大小的影響機制歸納為：

抬頭角度增大→來流在楔的上表面發(fā)生膨脹→馬赫數(shù)增大→分離區(qū)減??；

低頭角度增大→來流在楔的上表面發(fā)生壓縮→馬赫數(shù)減小→分離區(qū)增大。

表3 來流馬赫數(shù)的取值情況Table3 Value of freestream Mach number

圖10 壓縮拐角分離區(qū)大小隨馬赫數(shù)的變化曲線Fig.10 Variation of separation size with Mach number

3.3 周期性俯仰振動

令壓縮拐角按照方程（4）做周期性俯仰振動，振幅取θm＝4°，無量綱減縮頻率k＝0.0314，對應的振動周期為T＝200。

圖11是前3.5個振動周期內(nèi)俯仰角、分離點Δs和再附點Δr隨時間的非定常變化過程。壓縮拐角做周期性俯仰振動時，分離點和再附點與俯仰角同樣形成了周期性變化，而且頻率相同。由于非定常遲滯作用的影響，分離點和再附點的相位比俯仰角落后約3π／10。圖12是一個振動周期內(nèi)8個不同相位角上（φ＝0，π／4，π／2，3π／4，π，5π／4，3π／2，7π／4）分離區(qū)的非定常流動情況。整體上模型下俯的過程中分離區(qū)變大，模型上仰的過程中分離區(qū)變小，這與前面等速抬頭和等速低頭運動得到的結論相符。

圖11 壓縮拐角周期性俯仰振動過程中俯仰角、分離點和再附點隨時間的變化過程Fig.11 Variation of pitching angle，separation point and reattachment point during the forced oscillation

圖13是分離區(qū)大小Δ隨俯仰角θ的變化曲線，可以看出周期性振動情況下曲線具有明顯的遲滯環(huán)特征。與圖中模型固定在各角位置上不動時的定常計算結果（static）相比，周期性俯仰振動時（periodical oscillation）分離區(qū)的大小盡管整體趨勢上仍然是模型下俯（down）的過程中分離區(qū)增大、上仰（up）的過程中分離區(qū)減小，但Δ的最大值和最小值并不是剛好出現(xiàn)在最大俯仰角（θ＝4°）和最小俯仰角（θ＝－4°）位置上，而是滯后了一些。根據(jù)圖11，這個相位差大約是3π／10。

圖12 壓縮拐角一個俯仰振動周期內(nèi)8個不同相位角上的分離區(qū)流動情況Fig.12 Flow of separation region in 8phases during aperiod of pitching oscillation

圖13 周期性振動情況下分離區(qū)大小隨俯仰角的變化曲線Fig.13 Variation of separation size with pitching angle

考察周期性俯仰振動的振幅θm對分離區(qū)大小和相位的影響，令方程（4）中減縮頻率k＝0.0314不變，角振幅分別取θm＝2°、4°、8°和12°。圖14是不同角振幅下分離區(qū)大小Δ隨俯仰角θ的變化曲線。可以看出振幅的改變對曲線相位基本沒有影響，但是對分離區(qū)的大小變化影響很大。振幅增大時，Δ的最大值變得更大，同時最小值也變得更小，也就是說分離區(qū)整體的變化范圍變得更寬。

另外考察振動頻率k的影響，令運動方程（4）中角振幅θm＝4°不變，無量綱減縮頻率分別取k＝0.0157、0.0314、0.0628和0.1256。圖15是不同頻率下分離區(qū)大小Δ隨俯仰角θ的變化曲線。不同頻率下分離區(qū)大小的最大值和最小值差異不大，但不同頻率的曲線間的相位差異十分明顯。在低頻率下（k＝0.0157），模型的運動速度較小，非定常遲滯效應相對較弱，與模型靜止不動時（static）的情況類似，最大和最小分離區(qū)分別出現(xiàn)在最大和最小俯仰角上；當振動頻率增大一些后（k＝0.0314），模型的運動速度加快，遲滯效應增強，流場的變化開始落后于模型的運動，所以最大分離區(qū)和最小分離區(qū)沒有出現(xiàn)在最大俯仰角和最小俯仰角上，而是滯后了一些；當振動頻率進一步增大后（k＝0.0628、0.1256），遲滯效應增強，曲線間的相位差越來越大，最大、最小分離區(qū)的角位置甚至與模型不動時完全顛倒?？梢娬駝宇l率的改變對分離區(qū)大小的相位變化影響很大。

圖14 振幅θm對分離區(qū)大小的影響Fig.14 Effect of amplitude on separation zone

圖15 振動頻率k對分離區(qū)大小的影響Fig.15 Effect of frequency on separation zone

4 結 論

本文采用五階精度WCNS格式和非定?！半p時間步”格式數(shù)值求解非定常雷諾平均N－S方程研究了二維壓縮拐角模型在姿態(tài)運動情況下的非定常激波／邊界層干擾問題，主要結論如下：

（1）壓縮拐角等速抬頭運動時分離區(qū)減小，等速低頭運動時分離區(qū)增大；周期性俯仰振動時分離點和再附點的位置以及分離區(qū)的大小形成與模型同頻率的周期性變化。

（2）通過與相同來流條件下定姿態(tài)計算結果進行對比分析表明，壓縮拐角運動時模型姿態(tài)角的變化（引起了流動條件的變化）是導致分離區(qū)大小變化的主要因素，而模型運動的非定常遲滯作用會使這種變化減慢。

（3）模型姿態(tài)角的改變引起來流馬赫數(shù)發(fā)生變化，從而對分離區(qū)大小產(chǎn)生了較大影響，其規(guī)律是：壓縮拐角抬頭（低頭）→氣流發(fā)生膨脹（壓縮）→來流馬赫數(shù)增大（減小）→分離區(qū)減?。ㄔ龃螅?。

（4）周期性俯仰振動的振幅增大后，分離區(qū)大小的變化范圍明顯擴大，對相位影響不大；振動頻率的改變對分離區(qū)大小的變化范圍影響不大，但對分離區(qū)大小變化的相位的影響十分明顯。

［1］FERRI A.Experimental results with airfoils tested in the high speed tunnel at Guidonia［R］.NACA TM－946，1940.

［2］LI S X.The complex flow dominated by shock wave and boundary layer［M］.Beijing：Science Press，2007.（in Chinese）李素循.激波與邊界層主導的復雜流動［M］.北京：科學出版社，2007.

［3］ZHANG H J，MA H Y，TONG B G.Turbulence modeling validation in hypersonic complex flows［J］.ACTA Aerodynamica Sinica，2001，19（2）：210－216.（in Chinese）張紅杰，馬暉揚，童秉綱.高超聲速復雜流動中湍流模式應用的評估［J］.空氣動力學學報，2001，19（2）：210－216.

［4］DOLLING D S.Fifty years of shock－wave／boundary－layer interaction research：What next？［J］.AIAA Journal，2001，39（8）：1517－1531.

［5］KNIGHT D，JOSE L，DIMITRIS D，et al.Assessment of CFD capability for prediction of hypersonic shock interactions［J］.Progress in Aerospace Sciences，2012，48－49：8－26.

［6］KNIGHT D，HONG Y，ARGYRIS G P，et al.Advances in CFD prediction of shockwave turbulent boundary layer interactions［J］.Progress in Aerospace Sciences，2003，39（2－3）：121－184.

［7］ZHELTOVODOV A.Some advances in research of shock wave turbulent boundary layer interactions［R］.AIAA 2006－496.

［8］JACK R E.Numerical simulations of shock／boundary layer interactions using time dependent modeling techniques：a survey of recent results［J］.Progress in Aerospace Sciences，2008，44（6）：447－465.

［9］WU Y，YI S H，CHEN Z，et al.Experimental investigations on structures of supersonic laminar／turbulent flow over a compression ramp［J］.Acta Phys.Sin.，2013，62（18）：1－12.（in Chinese）武宇，易仕和，陳植，等.超聲速層流／湍流壓縮拐角流動結構的實驗研究［J］.物理學報，2013，62（18）：1－12.

［10］DENG X G，ZHANG H X.Developing high－order weighted compact nonlinear schemes［J］.Journal of Computational Physics，2000，165：22－44.

［11］SETTLES G S，F(xiàn)ITZPATRICK T J，BOGDONOFF S M.Detailed study of attached and separated compression corner flowfields in high reynolds number supersonic flow［J］.AIAA Journal，1979，17（6）：579－585.

［12］SETTLES G S，DODSON L J.Hypersonic shock／boundarylayer interaction database［R］.NASA CR－177577，1991.

［13］SETTLES G S，DODSON L J.Hypersonic shock／boundarylayer interaction database：new and corrected data［R］.NASA CR－177638，1994.

［14］SMITS A J，MUCK K C.Experimental study of three shock wave／turbulent boundary layer interactions［J］.Journal of Fluid Mechanics，1987，130：291－314.

［15］ZHAO R，YAN C.Evaluation of engineering turbulence models for complex supersonic flows［J］.Journal of Beijing University of Aeronautics and Astronautics，2011，37（2）：202－205.（in Chinese）趙瑞，閻超.超聲速復雜流動中湍流模型的評估［J］.北京航空航天大學學報，2011，37（2）：202－205.

［16］OLIVER A B，LILLARD R P，BLAISDELL G A，et al.Validation of high－speed turbulent boundary layer and shock－boundary layer interaction computations with the overflow code［R］.AIAA Paper 2006－894.