邊界層表面熱流的三維數(shù)值計(jì)算

邊界層表面熱流的三維數(shù)值計(jì)算

0熱流方程的數(shù)值模擬正確預(yù)測(cè)高超速飛機(jī)的表面溫度，并選擇適當(dāng)?shù)目篃岽胧?。這是確保飛機(jī)安全飛行并研究高超速飛機(jī)維修的一項(xiàng)重要技術(shù)之一。預(yù)測(cè)氣動(dòng)熱方法主要有地面與飛行試驗(yàn)、工程計(jì)算、數(shù)值求解Navier-Stokes方程及其近似形式等。工程算法為獲得表面熱流,在計(jì)算激波形狀和壓力分布時(shí)引入了不同程度的近似。基于跟蹤流線的軸對(duì)稱(chēng)比擬法是目前國(guó)外應(yīng)用最廣的一種計(jì)算有攻角三維物體繞流的氣動(dòng)加熱表面熱流的工程算法。其基本原理是:物體表面邊界層內(nèi)流體流動(dòng)的方向基本與繞物體的三維無(wú)黏流表面流線方向一致,在與流線垂直并平行于物面方向的流動(dòng)速度與主流速度相比很小,即小橫向流。作為近似若略去小橫向流,采用Manger變換,在物面以無(wú)黏表面流線作為正交坐標(biāo)系的一個(gè)坐標(biāo)軸,三維邊界層方程就可簡(jiǎn)化為軸對(duì)稱(chēng)比擬邊界層方程,故可用軸對(duì)稱(chēng)零攻角物體邊界層內(nèi)熱流近似計(jì)算公式計(jì)算有攻角三維物體沿流線的熱流分布。該方法由DEJARNETTE,HAMILTON在20世紀(jì)70年代提出,THOMPSON,RILEY,DEJARNETTE等不斷加以發(fā)展和完善,至20世紀(jì)90年代該工程算法更實(shí)用,計(jì)算結(jié)果更有效。但研究表明,該純工程算法在預(yù)測(cè)復(fù)雜外形飛行器表面熱流、保證局部區(qū)域預(yù)測(cè)精度等方面存在局限。隨著計(jì)算流體動(dòng)力學(xué)(CFD)研究的進(jìn)展,以及計(jì)算機(jī)技術(shù)發(fā)展實(shí)現(xiàn)的不斷增長(zhǎng)的速度和內(nèi)存,已經(jīng)發(fā)展了Navier-Stokes方程及其近似形式的計(jì)算程序。單個(gè)計(jì)算狀態(tài)下,該數(shù)值算法較工程算法計(jì)算時(shí)間長(zhǎng),占用內(nèi)存更多。在飛行器設(shè)計(jì)的方案階段,常需計(jì)算飛行器在整個(gè)軌道上不同時(shí)間各種飛行高度、馬赫數(shù)、攻角的多個(gè)狀態(tài)氣動(dòng)熱流分布,因此完全的數(shù)值模擬方法在工程上的實(shí)際應(yīng)用還存在一定限制。本文基于跟蹤流線的軸對(duì)稱(chēng)比擬法,采用數(shù)值求解Euler方程與邊界層內(nèi)工程算法相結(jié)合的方法,即由于高超聲速下邊界層很薄,忽略黏性項(xiàng)的影響,在邊界層外緣無(wú)黏流場(chǎng)數(shù)值求解Euler方程,得到更準(zhǔn)確、適用性廣的物面壓力分布,將其作為外緣條件用于邊界層內(nèi)跟蹤流線的軸對(duì)稱(chēng)比擬工程算法,計(jì)算高超聲速有攻角再入鈍頭體的表面熱流。1理論和公式1.1尺度因子h的控制方程為將軸對(duì)稱(chēng)比擬用于三維邊界層,先確定物面無(wú)黏流線,以及相應(yīng)沿每條流線的尺度因子。在物體幾何坐標(biāo)系(一般用柱坐標(biāo))中物面無(wú)黏流線方程可表示為dθdS=-(psρs(v∞)2)(ρsρ(v∞)2v2)×[-sinθcosΓ??x(pps)+(cosθcosδφ+sinθsinδφsinΓ)f??φ(pps)]-sinΓ[cosθcosΓ?σ?x+(sinθcosδφ-cosθsinδφsinΓ)f?σ?φ].(1)dθdS=?(psρs(v∞)2)(ρsρ(v∞)2v2)×[?sinθcosΓ??x(pps)+(cosθcosδφ+sinθsinδφsinΓ)f??φ(pps)]?sinΓ[cosθcosΓ?σ?x+(sinθcosδφ?cosθsinδφsinΓ)f?σ?φ].(1)式中:θ,φ為流線幾何坐標(biāo);S為弧長(zhǎng);ps,ρs分別為駐點(diǎn)壓力和密度;v∞為自由來(lái)流速度;v,p,ρ分別為來(lái)流速度、壓力和密度;Γ,δφ,σ為物體角;f為物體半徑。流線上的幾何坐標(biāo)x隨流線弧長(zhǎng)S的變化規(guī)律為dxdS=cosθcosΓ;(2)dφdS=sinθcosδφ-cosθsinδφsinΓf.(3)dxdS=cosθcosΓ;(2)dφdS=sinθcosδφ?cosθsinδφsinΓf.(3)尺度因子h的控制方程為1hd2hdS2=-[psρs(v2∞)2ρsρ(v∞)2v21h??β(pps)]2(3-(Μa∞)2)+psρs(v∞)2ρsρ(v∞)2v21h×??β[1h??β(pps)]+cos2Γcosδφf(shuō)×[?Γ?x?σ?φ-?σ?x?Γ?φ].(4)1hd2hdS2=?[psρs(v2∞)2ρsρ(v∞)2v21h??β(pps)]2(3?(Ma∞)2)+psρs(v∞)2ρsρ(v∞)2v21h×??β[1h??β(pps)]+cos2Γcosδφf(shuō)×[?Γ?x?σ?φ??σ?x?Γ?φ].(4)式中:Ma∞為自由來(lái)流馬赫數(shù);β為物面上流線的法向坐標(biāo)。由式(1)～(4)可知:在物體形狀(物體角Γ,δφ,σ;物體半徑f)、自由來(lái)流條件(ρ∞,v∞,Ma∞)和物面壓力分布(p,ps)已知后,只需指定流線步長(zhǎng)dS,沿每條流線數(shù)值積分,即可唯一確定流線離開(kāi)駐點(diǎn)后的空間走向(柱坐標(biāo)系中流線幾何坐標(biāo)θ,x,φ)。所得h相當(dāng)于等效旋成體在(x,φ)處對(duì)應(yīng)的半徑。由等效旋成體半徑沿軸向的變化規(guī)律,即可利用軸對(duì)稱(chēng)邊界層結(jié)果進(jìn)行計(jì)算。1.2基于修正牛頓理論的求解無(wú)黏表面流線和熱流的計(jì)算均依賴(lài)于物面壓力分布。提供表面壓力分布的方法有修正牛頓理論、無(wú)黏Euler方程求解等多種。修正牛頓理論易于使用,對(duì)某些簡(jiǎn)單幾何形狀的壓力分布預(yù)測(cè)相當(dāng)準(zhǔn)確,但不能用于有平直段的物形,對(duì)更復(fù)雜的物體幾何外形,其應(yīng)用更受限制。本文采用數(shù)值求解三維可壓縮定常Euler方程的方法計(jì)算物面壓力分布,該法適于復(fù)雜的幾何外形,所得壓力分布數(shù)據(jù)更為詳盡準(zhǔn)確。1.3超聲速條件下邊界層以打造聲速地層、聲速ae、再生聲速力;為計(jì)算熱流,需獲知邊界層外緣無(wú)黏參數(shù)(壓力pe、密度ρe、焓he、速度Ue、聲速ae、黏性系數(shù)μe)。一階邊界層近似允許邊界層外緣壓力ps等于相應(yīng)的物面壓力pw,在高超聲速條件下,邊界層很薄,可假設(shè)邊界層外緣熵等于正激波后的熵。由表面壓力分布和駐點(diǎn)滯止參數(shù)計(jì)算其他參數(shù)。1.4壓力分布和駐點(diǎn)滯止參數(shù)邊界層外緣參數(shù)的計(jì)算僅依賴(lài)于表面壓力分布和駐點(diǎn)滯止參數(shù)。利用來(lái)流條件p∞,ρ∞,v∞,溫度T∞和正激波關(guān)系可求駐點(diǎn)滯止參數(shù)ps,ρs。1.5壓力公式本文取LEES根據(jù)冷壁軸對(duì)稱(chēng)體加熱率關(guān)系獲得的駐點(diǎn)熱流qw,s=1.02×10-3√B+12×[1+0.527ˉβ0.686s1.116+0.411ˉβ0.686s](Ρr)-0.6√[?vΤ?SΤ]s×(ρwμw)as(ρeμe)bs(hs-hw).(5)qw,s=1.02×10?3B+12????√×[1+0.527βˉ0.686s1.116+0.411βˉ0.686s](Pr)?0.6[?vT?ST]s??????√×(ρwμw)as(ρeμe)bs(hs?hw).(5)式中:ˉβs為駐點(diǎn)區(qū)壓力梯度參數(shù),對(duì)球形頭部ˉβs=0.5;B為駐點(diǎn)處物面當(dāng)?shù)貦M、縱向曲率半徑之比;Pr為普朗特?cái)?shù);μw,μe分別為邊界層外緣和物面黏性系數(shù);hw為物面焓;a=0.1-0.08(ˉβs-0.5);(6)b=0.5-a;(7)[?vΤ?SΤ]s=1RΤ√2(ps-p∞)ρs.(8)此處:RT為駐點(diǎn)處物面當(dāng)?shù)貦M向曲率半徑;下標(biāo)s,T分別表示駐點(diǎn)和橫向。1.6壓力梯度參數(shù)對(duì)非駐點(diǎn)區(qū),因重要的峰值加熱區(qū),邊界層為層流,而湍流邊界層的理論和半經(jīng)驗(yàn)理論的發(fā)展常要借鑒層流區(qū)的處理方法。本文研究給定的流動(dòng)狀態(tài)為層流,主要討論層流邊界層的加熱。對(duì)等溫壁,有ζ′wζ′w,s=[1.116+0.411(ˉβs)0.6861+0.527(ˉβs)0.686]×[1+0.527ˉβ0.6861.116+0.411ˉβ0.686]×(1.1-0.1625te+0.0625(te)2)×0.85+0.15te-ζw(1-ζw,s);(9)qwqw,s=ppsUev∞hζ′wζ′w,s[2(B+1)v∞(?vΤ?SΤ)s∫s0ppsUev∞h2ds]1/2.(10)式中:ˉβ為壓力梯度參數(shù);te=hehs;ζw=ζw,s=hwhs。2遠(yuǎn)場(chǎng)邊界條件及基本控制方程以文獻(xiàn)的鈍雙錐為算例,幾何形狀如圖1所示。鈍雙錐的網(wǎng)格劃分如圖2所示。采用結(jié)構(gòu)化網(wǎng)格,可更易使用高精度格式,網(wǎng)格節(jié)點(diǎn)間關(guān)系更直接。物面附近網(wǎng)格加密,調(diào)整網(wǎng)格使其盡量正交、光滑,整個(gè)解域網(wǎng)格疏密變化均勻。設(shè)參數(shù)為:狀態(tài)1,Ma∞=9.86,p∞=64.08Pa,溫度T∞=48.96K,v∞=1383m/s,攻角(自由來(lái)流與后錐軸線夾角)α=16°,p∞=0.004560kg/m3,壁面溫度Tw=300K,每米雷諾數(shù)Re/m=2.472×105;狀態(tài)2,Ma∞=9.86,p∞=60.58Pa,T∞=49.23K,v∞=1387m/s,α=4°,ρ∞=0.004287kg/m3,Tw=300K,Re/m=2.315×105。邊界區(qū)域標(biāo)號(hào)如圖2所示。對(duì)稱(chēng)面邊界(區(qū)域1)上,平行于對(duì)稱(chēng)面的速度分量和標(biāo)量是關(guān)于對(duì)稱(chēng)面對(duì)稱(chēng)。物面邊界(區(qū)域2)條件滿足無(wú)穿透條件,即物面的法向速度分量為零。遠(yuǎn)場(chǎng)邊界條件采用應(yīng)用廣泛的Riemann不變量關(guān)系處理。對(duì)高超聲速流動(dòng),進(jìn)口邊界條件(區(qū)域3)指定為自由來(lái)流狀態(tài);出口邊界條件(區(qū)域4)由內(nèi)場(chǎng)外插值獲得?？刂品匠虨槿S可壓縮定常Euler方程。采用有限體積格心控制法,將計(jì)算空間離散為有限體積的小單元,并將算得的變量值存儲(chǔ)在單元中心?？臻g離散選用高階Roe’sFDS格式(FDS是典型的線性化Riemann解,對(duì)線性波有高分辨率),因而對(duì)黏性分辨率較高。采用Osher-C(L)限制器提高格式的精度,定熵參數(shù)對(duì)線性與非線性波都設(shè)為0.2。時(shí)間離散采用PointJacobi全隱格式及Backward-Euler格式。球頭鈍錐壓力等值線計(jì)算結(jié)果如圖3所示。數(shù)值求解Euler方程方法所得鈍雙錐表面迎風(fēng)子午線(φ=180°)和背風(fēng)子午線(φ=0°)上的無(wú)量綱壓力沿軸向分布分別如圖4、5所示,描述了球頭與錐體交接、前錐與后錐交接處的氣流的過(guò)膨脹和隨后再壓區(qū)后的錐體上的壓力分布。該方法采用數(shù)值方法求解高超聲速流動(dòng)繞流復(fù)雜幾何外形物體的流場(chǎng),可以得到很好的壓力分布結(jié)果。以此壓力分布作為跟蹤流線的軸對(duì)稱(chēng)比擬方法計(jì)算中的壓力分布輸入條件,計(jì)算鈍雙錐表面的熱流分布,所得鈍雙錐表面沿迎風(fēng)子午線(φ=180°)、背風(fēng)子午線(φ=0°)上的無(wú)量綱熱流沿軸向分布分別如圖6、7所示。由圖可知:準(zhǔn)確刻畫(huà)了交接區(qū)域膨脹波引起的熱流變化,給出的熱流計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合很好,表明方法正確、有效。3實(shí)驗(yàn)結(jié)果分析本文采用數(shù)值方法求解壓力分