楊柳青，李 沁，尤延鋮

（廈門大學航空航天學院，廈門 361102）

在超聲速、高超聲速飛行器上的內外流中，廣泛存在因激波/邊界層干擾（Shock-wave/boundary layer interaction，SWBLI）而引起的復雜分離流動現象。SWBLI 所產生的逆壓梯度會引起流動分離和總壓損失，增大流動的不穩(wěn)定性，導致進氣未啟動等情況。這些流動現象將會給飛行器的飛行性能、表面熱防護、結構設計等方面帶來直接影響。

近年來，人們對三維激波/邊界層相互干擾有了相當程度的了解，但目前的研究大多集中于后掠激波與邊界層的相互干擾［1-5］。錐形激波/邊界層的相互干擾（Concial shock-wave/boundary layer interaction，CSWBLI）在高速飛行器上也很常見，例如助推火箭頂部整流罩產生的錐形激波與火箭主體表面的邊界層相互作用，但是這方面的研究相對較少。Migotsky 等［6］首先對CSWBLI 進行了數學分析；Panov［7］開展了錐形激波與平板邊界層干擾的實驗研究；Kussoy 等［8］在來流馬赫數為2.2 的情況下進行了圓錐激波與軸對稱邊界層的實驗研究；Gai 等［9］通過實驗研究了來流馬赫數為2 時，錐形激波與平面湍流邊界層的相互作用；Hale［10］研究了來流馬赫數為2.05 時錐形激波與平面湍流邊界層的相互干擾情況；Zuo 等［11］模擬了Hale［10］的實驗模型和流場條件，進行了高保真的直接數值模擬；之后，Zuo 等［12］又采用雷諾平均NS 方程（Reynolds averaged Navier-Stokes，RANS）的方法，數值模擬了在自由流馬赫數4.0 條件下，由軸對稱壁面發(fā)展而來的內錐形激波/邊界層的相互干擾。這項工作旨在擴展先前Zuo 等［11］對外流CSBLI 的研究。利用直接數值模擬數據庫［11］，Zuo［13］還進一步研究了CSWBLI 中分離泡的尺寸及其與壓升的關系，建立了中等激波強度下壁面壓升與分離泡幾何形狀的標度關系。

在已開展的針對CSWBLI 的研究工作中，干擾流場的實驗結果缺乏足夠的細節(jié)，更多雷諾數變化對干擾結果的影響有待研究；同時，缺乏三維計算和二維計算結果的對比分析以顯示三維效應的影響。為了補充這方面的研究，本文采用兩方程Menter-SST 模型，以Gai 等［9］所采用的實驗模型和流場條件為基礎，對圓錐激波/平板邊界層干擾進行了數值模擬。通過定性與定量分析，進一步研究了入射激波與邊界層交匯后干擾區(qū)的流動情況；同時還分別改變了激波發(fā)生器頂角和來流單位雷諾數以比較不同工況下的干擾結果，并探討其流動分離的變化規(guī)律。此外還進行了二維情形的計算，將其與相應的三維結果進行了對比，分析兩者間的差異，研究三維效應對數值模擬結果的影響。

1 研究對象及數值方法

1.1 數值方法

本文數值模擬使用風雷軟件作為計算平臺，控制方程為RANS 方程，其積分形式為

式中：Q為流動變量，n為控制體表面的外法向向量，FI和FV分別表示對流（無黏）通量和黏性通量。湍流模型選用兩方程Menter-SST 模型，其具體表達式見文獻［14］。時間離散格式采用Jameson［15］提出的雙時間步隱式格式（Lower-upper symmetric-gauss-seidel，LU-SGS）。

1.2 研究對象

1.2.1 幾何外形

本文的研究對象為Gai 等［9］所使用的由圓錐和平板組成的物理模型，其幾何外形如圖1 所示。

圖1 幾何外形示意圖Fig.1 Geometrical sketch map

圓錐為激波發(fā)生器，圓錐半錐角θc可變，圓錐底面直徑為20 mm；平板長279 mm，圓錐中心距平板高度h= 30 mm。為了方便比較，以平板前緣為原點，保證各計算條件下入射激波的入射點始終距離原點xc=169 mm。

1.2.2 網格及邊界條件

為了方便生成網格并減少網格規(guī)模，對物面外形進行了簡化處理：考慮到圓錐激波發(fā)生器上半部分對于激波/邊界層干擾沒有影響，因而生成的網格只包含圓錐下半部分；由于物面外形及其周圍流場沿中心對稱面z= 0 對稱，因此僅生成對稱面一側的網格；激波發(fā)生器在風洞實驗中都是有限長度的，而在簡化網格中，將其底面一直延伸到了計算域的出口處，這樣處理雖然會使激波發(fā)生器在原底面位置處形成的膨脹波系結束角度不同，但是由特征線理論［16］可知，不同轉折角下產生的膨脹波在偏轉同一角度的區(qū)域內的解相同，因此不會給激波/邊界層干擾區(qū)帶來額外影響。

計算域長寬高分別為190、55 和30 mm，起始位置（即來流入口）為平板前緣。根據圓錐高度和激波角，得到三維計算下圓錐頂點至來流入口的距離lc1，見表1。

表1 三維計算下的lc1Table 1 lc1 in three-dimensional cases

以θc=20°的圓錐（即算例4）為例，模型的計算域如圖2 所示。

圖2 算例4 的計算域示意圖Fig.2 Sketch map of computational region for Case 4

計算網格采用結構網格，網格總體圖見圖3（a），網格總數約為540 萬個。對變形較大的位置做了過渡處理且基本確保壁面位置的網格正交性，如圖3（b，c）所示。平板上采用加密網格，平板的物面網格數為198×156（流向×展向），第一層網格厚度為0.001 5 mm（y+為2），網格增長比率為1.064 4。

圖3 算例4 網格示意圖Fig.3 Grid for Case 4

已知受三維效應的影響，入射激波角等效的二維流動與三維流動之間存在差異。為了進一步研究這些差異隨流動參數改變的變化規(guī)律，本文還安排了與三維問題對應的二維算例，其入射激波角度與三維中心對稱面上入射激波的角度相同。二維計算中，圓錐簡化為斜楔。根據激波理論，對于同一激波角，三維計算和二維計算各自需要的頂角大小不同，激波發(fā)生器位置相應發(fā)生變化。此外，在實際計算過程中發(fā)現如果仍采用三維計算時的激波發(fā)生器高度，在二維計算時會出現壅塞現象，因此二維計算時需提高激波發(fā)生器的高度。根據后續(xù)測試，在二維情況下，斜楔高度設置為40 mm。根據斜楔高度和激波角，得到二維計算下相應的斜楔楔角θ2大小及其頂點與來流入口距離lc2，如表2所示。

表2 二維計算下的lc2Table 2 lc2 in two-dimensional cases

計算域入口采用超聲速來流條件；出口為超聲速出流，直接使用外推；平板壁面為無滑移、絕熱邊界條件。

1.2.3 計算工況

本文算例的入口來流狀態(tài)如表3 所示。

表3 入口來流狀態(tài)Table 3 States of incoming flow

2 計算結果與討論

2.1 網格收斂性和數值驗證

2.1.1 網格收斂性

使用算例4-2D，來流狀態(tài)Run 2 來進行網格收斂性分析。參與比較的計算網格（流向網格數×縱向網格數）如表4 所示。

表4 計算網格Table 4 Computational grid

計算結果如圖4 所示，圖4 給出了3 種網格下的壁面壓力分布和壁面摩阻系數分布，可以看到Grid 2 和Grid 3 計算得到的結果基本重合，而Grid 1 的結果與兩者存在一些差距。根據網格收斂性分析，可以認為Grid 2 和Grid 3 的網格計算收斂，因此在后面計算中，中心對稱面上網格采用Grid 2 密度的網格。

圖4 網格收斂性分析Fig.4 Analysis of mesh convergence

2.1.2 數值驗證

本文以文獻［9］提供的實驗數據作為數值驗證的標準，圖5 為流向截面壁面壓力數值模擬結果與實驗結果的比較，使用的是算例1，狀態(tài)Run 2。

圖5 中，實驗結果在分離引起的壓升位置上游存在明顯的波動，推測該情況是實驗測量過程中各種誤差積累的結果。對比結果顯示，使用風雷軟件計算得到的數值模擬結果與實驗測量值吻合得較好，壁面壓力分布的最大誤差控制在2%以內。隨著逐漸遠離中心對稱面（z= 0 mm），各截面壓升頂點逐漸降低，但每次降低的程度并不相等，越遠離對稱面，降低幅度越大；同時，各截面壓升起始位置逐漸向下游移動，與壁面壓升變化情況類似，不同展向位置截面向下游移動的程度并不一致，越遠離對稱面，壓升起始位置移動幅度越大。壓升峰值后的壓降是因為圓錐底面產生的膨脹波傳播到壁面導致的。

圖5 算例1、狀態(tài)Run 2 的流向截面壁面壓力分布Fig.5 Wall pressure distribution in flow direction section for Case 1, Run 2

圖6 為上游影響范圍的數值模擬結果與實驗結果的比較，圖中“ref”指相應的實驗數據，來自文 獻［9］，參 與 比 較 的 算 例 有Case 1、Case 2 和Case 3，初始來流狀態(tài)均為Run 2。上游影響范圍定義為壁面壓升的起始位置，可以表示壁面上的分離激波軌跡。從圖中可以看出，除了θc= 14°時展向位置z= 24 mm 的數值模擬結果比實驗結果有比較明顯的上游偏移外，其余數值模擬結果與實驗結果擬合得較好。由于邊界層的存在，黏性計算下的激波軌跡相較于對應角度的無黏計算下的入射激波軌跡，其位置普遍向上游移動。由于遠離對稱面的激波效應減弱，所有黏性計算下的激波軌跡相較于無黏軌跡，越靠近出口，軌跡的曲率逐漸減??；對于4 個黏性計算下的算例，其激波軌跡向上游移動距離的變化幅度并不相等，14°、16°、18°的 變 化 幅 度 大 致 相 等，而20°相 較 于18°的結果，變化幅度并不大。此外，所有在中心對稱面上的數值結果與實驗結果相比，要略微向下游偏移。

圖6 上游影響范圍與無黏計算下的入射激波軌跡Fig.6 Upstream influence range and incident shock wave trajectory under inviscid calculation

2.2 不同半錐角對錐形激波邊界層干擾的影響

計算不同半錐角時，來流狀態(tài)均為Run 2。

2.2.1 中心對稱面的壁面壓力分布比較

圖7 給出了不同θc角下二維、三維中心對稱面（Z= 0）的壁面壓力分布。

圖7 不同θc 角下二維、三維中心對稱面的壁面壓力分布比較Fig.7 Comparison of wall pressure distribution between 2D cases and 3D cases in center-symmetric plane at different θc

由圖7 中三維數值模擬結果間的對比表明，分離點隨θc逐漸增大而向上游移動，原因是不同θc下產生的分離渦大小不同，導致各自的壓升起始位置不同；壓力比峰值也隨著θc的逐漸提高而增大；此外，隨著θc的增大，壓升過程中逐漸顯現出與分離流有關的壓力平臺，θc越大，壓力平臺越明顯。二維數值模擬結果間的對比顯示出與三維計算類似的規(guī)律。由于在二維計算的計算域內，圓錐底部產生的膨脹波還沒有影響到壁面，因此圖7 的二維計算結果沒有表現出二次壓升后的壓降過程。

考慮到二維計算下沒有橫向流動，理論上二維計算下的壁面壓升比應該會比對應的三維計算結果大，但是對比表明實際計算結果并不全是如此，需要將干擾區(qū)的壓升分為一次壓升和二次壓升來分別討論。對于一次壓升幅度，除了θc=20°的計算結果外，其他3 個較小θc下的二維計算結果都接近對應的三維計算結果。一次壓升是由分離激波引起的［17］，而分離激波的強度與分離渦尺寸密切相關，當θc較小時，三維計算受三維效應的影響并不顯著，使得三維計算下中心對稱面上的分離渦尺寸與對應的二維結果相差并不大（具體見2.2.2 節(jié)），因此兩者的分離激波強度也相似，從而使得一次壓升幅度大致相同；而當θc=20°時，由于受到的三維效應更加顯著，二維計算的分離渦尺寸明顯大于三維結果，因此其產生的分離激波更強，從而使得二維計算的一次壓升幅度大于三維結果。綜合以上分析，可以認為存在一個閾值，當θc大于等于該值時，中心對稱面上會因為受到顯著的三維效應影響而使得二維計算的一次壓升幅度大于三維結果。壓力平臺后的二次壓升則與再附過程有關（受到再附激波等因素的影響）［17］，從圖7 可知，4 組計算在三維計算下的二次壓升增長率都要高于對應的二維結果，而θc=20°時出現三維計算的二次壓升幅度低于對應二維結果的情況，其原因是此時圓錐半錐角增大引起膨脹波向上游移動，使得再附過程更早地受到了膨脹波的影響。

2.2.2 中心對稱面的旋渦結構和渦量分布比較

圖8 給出了不同θc角下二維、三維中心對稱面上的截面流線和展向渦量分布比較。以下從分離渦的尺寸、拓撲結構和渦量分布方面對圖8 分別進行分析討論。

圖8 不同θc角下二維、三維中心對稱面上的截面流線和展向渦量分布比較Fig.8 Comparison of streamline and vorticity between 2D cases and 3D cases in center-symmetric plane at different θc

（1）分離渦尺寸

三維和二維計算結果比較表明，由于保持不同θc下激波入射點位置不變，分離渦的渦核位置基本一致（在x= 166～167 mm），而分離點隨θc的減小而逐漸向下游移動。結合表5 分析，分離渦的尺寸隨θc增大而逐漸增大，其中寬度變化比較平均（二維計算中按2～3 倍增長；三維計算中每次增長1.5～2 mm），而分離渦的高度在θc=16°、18°、20°時是按照2～4 倍逐漸增大的，但是θc=16°下的分離渦高度相較于14°，增長了10 倍。文獻［9］中，實驗結果表明θc=14°時沒有觀察到分離現象，但數值模擬結果表明，分離還是存在的，只是尺度非常小。此外，圖7 中的壓力平臺范圍顯然與分離渦的尺寸有直接關系，分離渦尺寸越大，壓力平臺就越明顯。

表5 不同半錐角（楔角）的二維/三維分離渦尺寸對比Table 5 Comparison of separated vortex between 2D cases and 3D cases at different θc

考慮到氣流的橫向逸出，二維計算得到的分離渦應該要比對應的三維計算結果大，但是在目前的4 組算例中，只有θc=20°的結果符合預期。其他二維計算下的分離渦頂點大多位于相應三維計算下的分離流線和再附流線之間（或者接近再附流線頂點），而且分離點位置變化不大。

（2）分離渦內流線的拓撲結構

二維計算中，流線在分離渦內圍繞渦核形成系列封閉流線，其邊界是分離流線。而在三維流線中，分離渦不再是類似的封閉形狀，而是開放形式；流線也不再是封閉的曲線，而是圍繞渦核的螺旋形曲線；起于分離點的分離流線和滯止于再附點的再附流線是不同的流線，在圖8 中可以看到兩者間存在明顯的間距。渦核處，流動向側向逸出，這是一種典型的三維效應。

（3）渦量分布

二維結果和三維結果反映出同樣的渦量大小分布規(guī)律，即分離渦內的渦量并不大，渦量較大的區(qū)域集中在邊界層之中，邊界層內渦量較大是因為邊界層內受黏性作用形成剪切流。此外，由于邊界層會受到分離渦的影響而發(fā)生抬升，因此對于不同θc的計算結果，分離渦附近的渦量分布與分離渦的尺寸有直接關系。而比較二維結果與相應的三維結果可以發(fā)現，兩者間的差異與各自的分離渦大小有直接關系，因此可以認為渦量分布也會受到三維效應的影響。

2.2.3x=169 mm 位置的展向截面密度梯度分布比較

圖9 給出了不同θc角下x=169 mm 的展向截面密度梯度云圖。從圖中可以分辨出較為清晰的來流邊界層（黑框）、入射激波（藍框）、反射激波（紅框）等流動結構。截面上，密度梯度的大小隨θc增大而增大；入射激波、反射激波等流動結構的尺寸也與θc密切相關，當θc增大時，截面上入射激波、反射激波、膨脹波系的尺寸均增大。

圖9 不同θc 角下x = 169 mm 的展向截面密度梯度Fig.9 Density gradient of spanwise section (x = 169 mm)at different θc

2.3 不同雷諾數對錐形激波邊界層干擾的影響

計算不同雷諾數時，使用的激波發(fā)生器外形均為Case 4。

2.3.1 中心對稱面的壁面壓力分布比較

圖10 給出了不同Re數下二維、三維中心對稱面的壁面壓力分布。圖10 顯示三維計算下分離點隨著Re增大而逐漸向下游移動，并且一次壓升后的壓力平臺也變得越來越不明顯，直至Re=3×109時，壓力平臺完全消失（說明該雷諾數下沒有發(fā)生分離）。這表明入射激波與邊界層相互作用產生的分離渦隨Re增大而逐漸減小，一方面是因為邊界層厚度減小，另一方面是因為高雷諾數下，邊界層內速度虧損小，剖面含有更高的動能，對逆壓梯度阻滯的抵抗能力更強。但是壓升峰值隨著Re增大而逐漸提高，且提高過程不完全是一個線性過程，Re= 3×106、3×107、3×108下的增長幅度大致相同，而Re= 3×109的壓升峰值相較于Re= 3×108的結果，增長幅度很大。二維計算也顯示出了與三維計算類似的規(guī)律，中心對稱面的壁面壓力變化情況總體上與分離渦的尺寸相匹配。

圖10 不同Re 數下二維、三維中心對稱面的壁面壓力分布比較Fig.10 Comparison of wall pressure distribution between 2D cases and 3D cases in center-symmetric plane at different Re

比較三維計算結果和對應的二維計算結果發(fā)現，除了Re= 3×109時的計算外，其余3 個雷諾數下的二維計算結果相較于對應的三維計算結果，分離點都顯著地向上游移動，且壓力比基本上都大于三維結果。Re= 3×109時，不僅沒有發(fā)生流動分離，而且二維計算下的壁面壓升峰值也不及三維計算。Re= 3×106、3×107、3×108時三維計算與二維計算之間的差異可以用三維效應來解釋，而Re= 3×109的計算結果則表明高雷諾數可能會減弱三維效應帶來的影響。

2.3.2 中心對稱面的旋渦結構和渦量分布比較

圖11 給出了不同Re數下中心對稱面的截面流線和展向渦量分布。以下從分離渦的尺寸、拓撲結構和渦量分布方面對圖11 分別進行分析討論。

（1）分離渦尺寸和拓撲結構

圖11 的結果表明，隨著Re增大，分離渦尺寸顯著減小，分離點向下游移動，同時渦核位置基本保持不變（因為控制了激波的入射位置）；直至雷諾數增長至Re= 3×109時，入射激波與邊界層的相互作用未發(fā)生分離。該規(guī)律對于二維計算和三維計算都成立，其原因與2.3.1 節(jié)中的分析相似。

圖11 不同Re 數下中心對稱面的截面流線和展向渦量分布比較Fig.11 Comparison of streamline and vorticity between 2D cases and 3D cases in center-symmetric plane at different Re

對比二維和三維的數值模擬結果，同樣可以發(fā)現二維計算中分離渦內的流線是封閉的，而三維情況則是螺旋形曲線；同時結合表6，二維計算下的分離渦尺寸要大于相應的三維計算結果，這些都體現了三維效應的影響。但是進一步比較二維、三維分離渦尺寸的變化，注意到隨著Re的增大，二維分離渦與三維分離渦之間的差距在逐漸縮小，直至Re= 3×109時都沒有發(fā)生流動分離，流線幾乎一致，這表明高雷諾數下中心對稱面上二維計算與三維計算的流動結構差異減小。

表6 不同Re 數的二維/三維分離渦尺寸對比Table 6 Comparison of separated vortex between 2D cases and 3D cases at different Re

另外，Re= 3×109時，二維、三維計算都沒有產生分離渦，考慮到高雷諾數對三維效應的抑制作用，是否存在一個雷諾數的閾值，可以使得三維計算時無分離渦，而二維計算時有分離渦，還需要進一步探索。

（2）渦量分布

與2.2.2 節(jié)的計算結果類似，二維、三維計算下渦量較大的區(qū)域均集中于邊界層之中，而分離渦內的渦量并不大，且渦量分布與分離渦的尺寸（影響邊界層狀態(tài)）有直接關系，因此渦量分布也受到雷諾數變化和三維效應的影響。

2.3.3x=169 mm 位置的展向截面密度梯度分布比較

圖12 給出了不同Re數下x= 169 mm 的展向截面密度梯度云圖。從圖12 可以看到較為清晰的來流邊界層、入射激波、反射激波等流動結構。密度梯度大小隨Re增大而增大；來流邊界層、反射激波、膨脹波系的尺寸變化情況則與之相反，當Re增大時，這些流動結構的尺寸均明顯減小。

圖12 不同Re 數下x = 169 mm 的展向截面密度梯度Fig.12 Density gradient of spanwise section (x = 169 mm)at different Re

2.3.4 壁面極限流線比較

圖13給出了不同Re數的壁面極限流線，極限流線圖取第一層網格高度位置的水平截面進行繪制。

圖13 不同Re 數的壁面極限流線Fig.13 Limiting streamlines of wall at different Re

除了Re= 3×109時的計算外，其他雷諾數下的計算都產生了反射激波誘導分離現象，相應的極限流線圖上都存在相似的鞍點（紅圈）、結點（藍圈）、分離流線（紅色虛線框）和再附流線（藍色虛線框）等典型特征結構。Re= 3×109的計算結果則呈現典型的無分離弱干擾流動特征。對于Re=3×106、3×107、3×108的 結 果，隨 著 雷 諾 數 增大，鞍點向下游移動，結點向上游移動，整體上分離流線和再附流線之間的間距在逐漸縮小，這一現象也反映出不同雷諾數下入射激波與邊界層相互作用而產生的分離渦尺寸的變化。

2.3.5 壁面上游影響比較

圖14 給出了不同Re數的上游影響范圍和無黏計算下的入射激波軌跡。

圖14 不同Re 數的上游影響范圍與無黏計算下的入射激波軌跡Fig.14 Upstream influence and inviscid shock trace at different Re

與圖6 顯示的結果類似，黏性條件計算下的分離激波軌跡相較于無黏激波軌跡，軌跡位置均向上游移動，且越靠近出口位置，軌跡的曲率越小。隨著Re增大，黏性計算下的激波軌跡向上游移動的距離逐漸減小，這與分離渦尺寸的變化情況相一致。

3 結 論

本文針對圓錐激波與平板邊界層的相互干擾，使用兩方程Menter-SST 模型進行了定常數值模擬。計算模擬了不同半錐角和雷諾數下的干擾流場，總結了不同參數下計算結果的變化規(guī)律，擴展了先前的實驗和數值研究內容；此外，將三維計算結果與對應的二維計算結果進行了對比，分析了三維效應對計算結果的影響。

（1）隨著半錐角的增大，入射激波與邊界層干擾產生的分離渦尺寸增大，上游影響范圍朝上游移動；壁面壓升峰值提高，一次壓升后的壓力平臺越發(fā)明顯；干擾區(qū)的密度梯度大小增大，顯示出的反射激波和膨脹波系等流動結構的尺寸也在擴大；中心對稱面上的渦量分布則與分離渦尺寸直接相關。受三維效應的影響，中心對稱面上的分離渦內流線的拓撲結構在二維計算中是封閉的，而在三維計算中則是開放的螺旋形曲線。此外，在分離渦尺寸和壁面一次壓升方面，除了θc=20°的結果，其他3 個較小θc下的二維計算結果均接近對應的三維結果，因此可能存在一個閾值，當θc大于該值時，三維效應的影響才會明顯顯現出來。

（2）隨著雷諾數的增大，入射激波與邊界層干擾產生的分離渦尺寸減小，直至不發(fā)生流動分離，上游影響范圍朝下游移動；壁面壓升峰值提高，但是一次壓升后的壓力平臺愈發(fā)不明顯；干擾區(qū)的密度梯度大小增大，但是顯示出的反射激波和膨脹波系等流動結構的尺寸減?。恢行膶ΨQ面上的渦量分布同樣與分離渦尺寸直接相關。除了拓撲結構上的區(qū)別，受三維效應的影響，二維計算下的分離渦尺寸要大于相應的三維計算結果，但同時二維分離渦與三維分離渦之間的尺寸差距會隨著Re的增大而逐漸縮小，這表明高雷諾數下中心對稱面上兩者流動結構的差異會減小。而在中心對稱面的二次壓升峰值方面，較小的3 個雷諾數下的二維結果大于三維結果，但是兩者間的差距會隨著Re的增大而逐漸縮??；當Re= 3×109時，三維情況下的二次壓升峰值大于二維結果。