邊界層表面熱流的三維數值計算

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