氣動載荷下防熱材料剝離顆粒輸運特性的直接數值模擬研究1)

李婷婷 李 青 涂國華 袁先旭 周 強 ,3)

* (西安交通大學化學工程與技術學院,西安 710049)

? (中國空氣動力研究與發展中心空氣動力學國家重點實驗室,四川綿陽 621000)

引言

高超聲速飛行器以極高速(Ma>10)飛行,將在近壁區產生激波壓縮、黏性摩擦、氣動加熱等作用,需要設計特定的熱防護系統以抵御氣動熱載荷.燒蝕防熱通過飛行器表面材料反應、熱解等將表面熱載荷帶離,是飛行器熱防護的重要形式之一[1-2].熱防護系統常用的纖維編制類、顆粒增強類復合材料(見圖1),在高焓氣動載荷下材料顆粒相或纖維束發生力學剝蝕,對飛行器流態、氣動性能、熱載荷產生影響[3-4].另外,依托電弧加熱類風洞設備開展材料響應評價實驗時,在極高焓電弧作用下,材料將發生微結構變形、蠕化、剝離等現象,形成的燒蝕顆粒進入流場,污染電弧風洞流場,進而干擾材料性能評價[5-6].因此,燒蝕顆粒的剝離動力學問題是高超聲速領域的共性基礎性科學問題,研究氣動載荷下燒蝕剝離顆粒的輸運特性對防熱材料體系評價和熱防護系統性能評估均具有重要而又緊迫的工程價值.

針對飛行器燒蝕及顆粒剝離問題,Li 等[7]研究了碳化復合材料燒蝕面后退的非線性熱解層模型,并預測了燒蝕響應;文獻[8]利用線性穩定理論研究了防熱材料的溫度響應及燒蝕氣體產物對高溫邊界層不穩定性的影響;國義軍等[9]在前人基礎上提出了氧化燒蝕的雙平臺理論和反應控制機理.上述研究重點關注燒蝕后退量工程計算、燒蝕剝離的力學/熱學條件、燒蝕機理等問題,未見有針對燒蝕顆粒在近壁運動的動力學機理以及對邊界層的調制等方面的研究.另一方面,作為一個經典的力學問題,顆粒多相流已在風沙、化工等顆粒流動問題中得到了較為深入的研究,并形成了顆粒解析的相關數值模擬方法和軟件[10-13].文獻[14]綜述了顆粒、液滴和氣泡與湍流的相互作用和調制,及其相關的實驗和模擬技術.Yu 和Shao[15]提出基于虛擬區域方法的顆粒懸浮流直接數值模擬,并證明了其在各種典型顆粒流動均適用并具有較好的精度和魯棒性.文獻[16]通過理論分析得出了球形懸浮顆粒在有界流體中的側向遷移現象.針對氣動載荷下防熱材料剝離顆粒動力學問題,亟需借鑒顆粒解析的模擬方法,進一步開展氣動載荷下防熱材料剝離顆粒輸運特性研究,以期揭示燒蝕顆粒在近壁運動的動力學機理以及對邊界層的調制規律,支撐高超聲速飛行器熱防護系統設計.

本文針對防熱材料燒蝕剝離問題,采用近壁流動量綱分析方法,將氣動環境下材料燒蝕顆粒剝離過程模化為單個圓球慣性燒蝕顆粒在Couette 流動中的動力學問題,采用顆粒解析的直接數值模擬方法(particle resolved-direct numerical simulation,PRDNS),開展燒蝕顆粒剝離和輸運特性原理性分析,為燒蝕質量損失的預測和邊界層調制研究提供理論基礎.

1 顆粒及流動模型

1.1 典型燒蝕顆粒模型

常用防熱燒蝕材料為纖維編織/顆粒填充增強復合材料.燒蝕過程分為熱化學燒蝕和機械剝蝕兩部分,前者指燒蝕材料與高溫氣體組分之間的化學反應以及燒蝕材料在高溫情況下的升華過程[17-18],后者指基體和基質密度不同造成材料燒蝕速率差異而暴露出的纖維束、黏結劑優先燒蝕而暴露出的顆粒填充物等,在氣流壓力以及剪切力作用下被剝離吹走[19-20],如圖1 所示.

圖1 纖維增強防熱材料在高溫氣流作用下的機械剝落示意圖Fig.1 Schematic diagram of mechanical spalling of fiber reinforced ablative composite materials under high temperature air flow

單根纖維的直徑一般1～ 10 μm,纖維斷裂一般是多根一起發生,產生的斷裂產物一般幾微米到幾十微米;顆粒填充物反應產物的典型粒徑為0.1～10 μm,因此,實際燒蝕顆粒的尺寸在微納米量級0.1～ 10 μm[21-22].當顆粒粒徑大于1 μm 時,只需考慮顆粒與流體相互作用力;當顆粒直徑小于1 μm時,還需要考慮顆粒之間的分子間作用力、顆粒的團聚作用等.實際燒蝕顆粒不僅粒徑尺寸多變,而且具有不同的長徑比,且由于燒蝕掉的材料不同,顆粒密度也不同.本文以真實數據為參考,研究單個圓球形慣性顆粒,直徑量級(1～10) μm (本文上標*代表有量綱物理量,其他為無量綱物理量,除了角度量θ,α,β).

1.2 典型流動模型

表1 激波前后來流參數范圍Table 1 Flow parameters beside the shock

飛行器頭部位置流動可模化為經典的楔形物體繞流問題,其控制方程為Falkner-Skan 方程,是一個典型的層流邊界層相似解問題.Falkner-Skan 方程[26]表述為

邊界條件為

式中 β 為壓力梯度,且 β 與等效半錐角 α 的關系為

根據 β 取不同的值,飛行器頭部流動可模化為三種不同的典型標模流動,如圖2 所示.

圖2 飛行器頭部流動示意簡圖Fig.2 Flow around the aircraft nose

(1) β=1,Hiemenz 流動

Hiemenz 流動區域有一個向下的擠壓力,正壓力梯度,加速流動.

(2) 0 <β<1,Falkner-Skan 流動

Falkner-Skan 流動區域無向下擠壓力,正壓力梯度,加速流動.

(3) β=0,Blasius 流動

Blasius 流動區域無向下擠壓力,無壓力梯度,平行剪切流,本文重點研究該區域.

邊界層外邊界處的勢流速度為

其中,x*為沿壁面位移,c為常數.

不同壓力梯度的Falkner-Skan 流動速度分布見圖3,本文研究的壁面燒蝕顆粒在微米量級,因此不妨先關注距離壁面10～100 倍顆粒直徑,即y*≤0.1 m 范圍內顆粒的啟動和輸運情況.

圖3 不同壓力梯度Falkner-Skan 流動速度分布圖Fig.3 Velocity profiles of Falkner-Skan flow at different half cone angles

在β=0 的平行剪切流近壁區y*≤0.1 m,相似變量

此時,剪切率B*接近常數,如圖3 紅色虛線所示,因此,作為一種原理性研究,可用剪切率為常數的Couette 流動來研究微米尺度顆粒在近壁面的輸運特性.根據定義,剪切率B*和速度u*的關系式為

自相似解求解過程中,有

其中,φ*為流函數,η 為相似變量,將式(8)代入式(7),并經過變換得剪切率

流體推動顆粒運動,顆粒運動的能量由流體慣性提供,表征顆粒周圍的流體慣性的特征參數為顆粒雷諾數

表征顆粒慣性的特征參數為斯托克斯數

即可保證數值實驗模擬真實物理過程的合理性.

真實情況下,顆粒直徑與特征長度的比為

顆粒密度與氣體密度比為

2 數值方法及驗證算例

傳統的PP-DNS (point particle-DNS)是對顆粒動力學過程的一種近似,其忽略了有限尺寸帶來的尾流效應,且需要的流體信息是基于未擾流動信息,而真實PP-DNS 求解器中的流場都是受擾動流場[27-28].

而PR-DNS 是基于第一性原理的數值模擬方法:由于網格尺寸比顆粒尺寸小至少一個數量級,使得流體求解器可以獲得顆粒周圍的小尺度流體動力學信息,然后根據柯西定理,對有限尺寸顆粒表面的流體動力學信息進行積分,從而獲得基于拉格朗日坐標下的顆粒瞬時受力,即

其中,Σ為有限尺寸顆粒控制體的表面積.

本研究采用浙江大學開發的基于并行計算的PR-DNS 代碼來進行數值實驗[29-30],研究單個慣性顆粒在流動近壁區的動力學問題,從而揭示其剝離及輸運特性.

2.1 中性顆粒在平面Couette 流中的驗證

驗證算例選取平面Couette 流中的中性浮力顆粒.首先在一個長方體計算域中構造一個數值Couette 流動,如圖4 所示.計算域為OX×OY×OZ=16rp× 16rp× 8rp,球形顆粒半徑rp=1 (代碼中所有量均為無量綱量),則上下壁面距離H=16,上壁面速度Uw=8,κ=rp/H.在計算域中構造流動,參數(詳見表2)如下:特征流體剪切率B=0.5,特征流體密度=1,特征運動黏度ν=1,特征長度為,特征速度,特征時間t=δ/U=2,來與理論解對比,驗證PR-DNS 的準確性和合理性.

圖4 長方體計算域及其相關參數示意圖Fig.4 Cuboid computing domain and parameters

表2 數值計算參數設置Table 2 Numerical calculation parameters settings

初始流場全場施加一個平面Couette流U=By,V=0,W=0,之后每個時刻,O Y施加Dirichlet 邊界條件,OX和OZ方向邊界施加周期邊界條件,均勻網格尺寸為 ?=rp/16 .已驗證結果與網格尺寸 ?、時間步長dt、計算域大小OX×OY×OZ無關.

計算結果圖5 表明,PR-DNS 數值方法預測的中性顆粒在Couette 流中的法向速度與文獻[31-32]理論解、Wang[33]數值模擬結果吻合良好.

圖5 Couette 流中單個中性顆粒法向速度隨顆粒到壁面位移的變化Fig.5 The normal velocity of a single neutral particle in Couette flow varies with the particle displacement to the wall

2.2 顆粒解析的直接數值模擬方法

構造一個數值Couette 流動,如表2 所示,其特征流體剪切率B=1,特征流體密度 ρf=1,特征運動黏度 ν=1,則特征長度為δ=,特征速度,特征時間t=δ/U=1 .為保證和S t與真實物理過程一致,數值實驗中的,dp/δ ～O(1),因此,ρr取值分別為10000,20000,30000,dp/δ 取值分別為0.5,0.75,1,1.5,2,4,6,8,10,12,14.

本文考慮數值計算域是一個長方體,如上述圖4 所示,計算域OX×OY×OZ=10rp× 10rp× 5rp,OX×OY×OZ=20rp×20rp×10rp,其中rp為無量綱顆粒半徑.網格尺寸為 ?=rp/6.4,?=rp/12.8,時間步長為dt=0.0025t,0.005t.已驗證PR-DNS 的數值模擬結果與計算域、網格尺寸、時間步長無關,其中一組網格無關驗證如圖6 所示,本文使用均勻網格.初始流場全場施加一個平面Couette 流U=By,V=0,W=0,之后每個時刻,OY施加Dirichlet 邊界條件,OX和OZ方向邊界施加周期邊界條件.

圖6 網格無關性驗證算例(ρr=10000,rp=0.25,OX×OY×OZ=5rp×5rp×2.5rp)Fig.6 Grid independent verification(ρr=10000,rp=0.25,OX×OY×OZ=5rp×5rp×2.5rp)

3 顆粒剝離和輸運特性模擬結果分析

3.1 密度比對顆粒輸運的影響

相同直徑不同密度比顆粒的輸運軌跡如圖7 所示.從圖7 可以看出,顆粒水平位移約是法向位移的20 倍,水平位移遠大于法向位移,即大密度比球形顆粒在Couette 流動近壁區以水平輸運為主.除此之外,密度比越大的顆粒,法向位移越大,即越重的顆粒被吹的越高.為了揭示個中原因,進一步研究了顆粒輸運速度的變化規律.

圖7 不同密度比顆粒輸運軌跡 (dp/δ =1)Fig.7 Particle transport trajectories with different density ratios(dp/δ =1)

相同直徑不同密度比顆粒的水平輸運速度沿流向變化規律如圖8 所示.從圖8 可以看出,隨著密度比 ρr的增大,顆粒水平輸運速度Up呈現減小趨勢.而相同直徑不同密度比顆粒水平滑移速度沿流向變化規律如圖9 所示.由圖9 可知,隨著密度比 ρr的增大,顆粒/流體之間的水平滑移速度呈增加趨勢.圖8 和圖9 共同佐證了越重的顆粒越難以被吹遠,這是符合認知規律的,根據St～ ρr(式(12)),在相同體積下,密度比 ρr越大,顆粒慣性St也越大,沿流向越難以被吹遠.

圖8 不同密度比顆粒水平速度沿流向變化規律 (dp/δ=1)Fig.8 The horizontal velocity of particles with different density ratios varies along the flow direction (dp/δ=1)

圖9 不同密度比顆粒水平滑移速度沿流向變化規律 (dp/δ=1)Fig.9 The horizontal slip velocity of particles with different density ratios varies along the flow direction (dp/δ=1)

進一步考察法向運動情況.相同直徑不同密度比顆粒的法向速度隨流向變化規律如圖10 所示.由圖可知,顆粒密度ρr越大,顆粒法向速度越小.相同直徑不同密度比顆粒法向位移隨時間變化曲線如圖11 所示.由圖可知,顆粒密度ρr越大,顆粒法向位移越小.上述現象說明相同體積下,越重的顆粒,其慣性St越大,沿法向也越難以被抬高.

圖10 不同密度比顆粒法向速度沿流向變化規律(dp/δ=1)Fig.10 Normal velocity of particles with different density ratios varies along the flow direction (dp/δ=1)

圖11 不同密度比顆粒法向位移隨時間變化規律(dp/δ=1)Fig.11 The normal displacement of particles with different density ratios varies with time (dp/δ=1)

然而,有趣的是從上述圖7 可觀察到大密度比的顆粒輸運軌跡反而在較高的位置,這可能是由于密度比對顆粒水平輸運速度的影響遠大于對法向速度的影響,即密度比的增加對顆粒水平輸運速度的抑制作用遠大于對法向速度的抑制作用.

3.2 直徑對顆粒輸運的影響

由于顆粒直徑dp越大,則顆粒慣性St(～dp)越大,而大顆粒的重心離開壁面的距離也越高,其重心處流體速度越大(Uf=Byp),顆粒受到流體的推力也越大.因而,直徑對顆粒輸運的影響必然是顆粒慣性和流體慣性競爭的結果.相同密度比、不同直徑的顆粒水平輸運速度如圖12 所示.從圖12 可以看出,直徑dp越大的顆粒,水平輸運速度Up越小,體現了顆粒慣性作用優勢.

圖12 不同直徑顆粒水平速度沿流向變化規律(ρr=10000)Fig.12 The horizontal velocity of particles with different diameters varies along the flow direction (ρr=10000)

然而,相同密度比、不同直徑的顆粒水平滑移速度沿流向變化規律,如圖13 所示,越大的顆粒,其周圍流體速度Uf=Byp也越大,這導致顆粒/流體之間的滑移速度就越大,體現出流體慣性作用優勢.

圖13 不同直徑顆粒水平滑移速度沿流向變化規律(ρr=10000)Fig.13 The horizontal slip velocity of particles with different diameters varies along the flow direction (ρr=10000)

圖14 為不同直徑顆粒法向速度沿流向變化規律,表明直徑越大的顆粒,法向速度越大,即越大的顆粒,越容易被吹高.圖15 為不同直徑顆粒的輸運軌跡,從中看出,顆粒直徑越大,法向位移越大.這是由于剪切流中的顆粒會受到向上的Saffman 升力FL*[34-35]

圖14 不同直徑顆粒法向速度沿流向變化規律(ρr=10000)Fig.14 The normal velocity of particles with different diameters varies along the flow direction (ρr=10000)

由圖15 還可以看出,相同密度比不同直徑的顆粒,水平位移是法向位移的10～ 20 倍,顆粒也是以水平輸運為主.圖16 進一步給出了受顆粒影響的流場壓力云圖,由圖16 可知,水平方向壓差明顯大于法向壓差,顆粒主要在較大的水平壓差作用下運動,以水平輸運為主.

圖15 不同直徑顆粒的運動軌跡 (ρr=10000)Fig.15 Trajectories of particles of different diameters (ρr=10000)

圖16 含顆粒的Couette 流場壓力云圖Fig.16 Pressure contour of particle influenced flow field

3.3 歸一化的顆粒輸運規律

由圖17 可知,無量綱的顆粒啟動速度隨著顆粒雷諾數Rep增加而增加.通過擬合,進一步得到顆粒啟動速度的歸一化表達式

圖17 不同密度比顆粒的非定常歸一化啟動速度隨顆粒雷諾數變化規律Fig.17 Dimensionless unsteady start-up velocity of particles along the particle Reynolds number for different density ratios

其中,無量綱剪切率B=B?/B∞=1 .式(19)清晰地表明了顆粒和流體慣性作用對顆粒輸運的影響.右邊第一項為流體微團或中性顆粒作為示蹤粒子在流體中的啟動速度[37-40];右邊第二項為慣性修正項,體現了流體慣性和顆粒慣性的耦合修正效應,Up0=,其中顆粒慣性St對顆粒啟動速度Up0起到了阻滯作用.

4 結論

本文通過單個球形慣性顆粒在Couette 流動的PR-DNS 模擬,獲得了不同顆粒/流體密度比、不同粒徑條件下燒蝕顆粒在近壁流動中的輸運規律.

(1)獲得了燒蝕顆粒關鍵特征參量對顆粒輸運的動力學規律:輸運速度與顆粒/流體密度比ρr和顆粒直徑dp有關.燒蝕顆粒密度比ρr增大,輸運速度因顆粒慣性St增大而減小;顆粒粒徑dp增大,水平輸運速度因顆粒慣性St增大而減小,法向速度和位移均因受到的Saffman 升力的增大而增大.燒蝕顆粒法向位移遠小于水平位移,顆粒以水平輸運為主.

(2) PR-DNS 得到燒蝕顆粒啟動速度的歸一化表達式,與理論分析一致,歸一化的顆粒啟動速度是顆粒和流體慣性的函數,這與前人的經典理論和數值模擬結果定性吻合[33,35,37].具體地,顆粒水平輸運速度等于流體微團或中性顆粒的速度減去顆粒和流體的慣性修正項.

本文研究對燒蝕顆粒動力學問題作了單個顆粒、球形顆粒、Couette 流動等假設,未考慮顆粒間相互作用、復雜顆粒形貌、近壁真實流動條件下的顆粒動力學機理.后續研究將進一步重點關注可壓縮性、顆粒群作用、橢球形等復雜顆粒形狀、真實流動等更符合實際物理過程的燒蝕顆粒動力學問題.

致謝

感謝浙江大學航空航天學院余釗圣教授提供Fortran 并行 PR-DNS 代碼支持.