气动载荷下防热材料剥离颗粒输运特性的直接数值模拟研究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 代码支持.