边界层表面热流的三维数值计算

边界层表面热流的三维数值计算

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实验结果分析本文采用数值方法求解压力分