基于风雷软件的可压缩气固颗粒两相流计算

雷颖昊南，雷三惠，张生浩，王 萍

(兰州大学 颗粒湍流研究中心，兰州 710000)

0 引 言

可压缩颗粒两相流常见于超燃发动机[1-2]、沙尘或雨中飞行器航行[3-4]和混合物冲击/爆炸[5]等。另一个典型的例子是高超声速飞行器再入大气过程中，热防护材料高温烧蚀“层裂”过程中产生颗粒[6-7]，与飞行器周围气流相互作用。因为单个颗粒可能在稀薄或连续气体中经历亚声速到高超声速的流态，这些类型气流中的颗粒运动建模很复杂。同时，由于两相系统的复杂性和多尺度性，用理论和实验方法进行研究存在诸多困难，数值模拟可直观地展示颗粒内部的物理现象，因此可发挥较大作用。

用于模拟颗粒两相流的方法大体上可以分为欧拉-欧拉、欧拉-拉格朗日两类。在欧拉-欧拉方法中，流体和颗粒均被视为连续介质，采用Navier-Stokes方程组求解[8]，该方法适用于颗粒体积分数较高的情况；而在欧拉-拉格朗日方法中，气体方程在欧拉坐标系下离散，颗粒在拉格朗日坐标系下跟踪，所受到的力可能包括拖曳力、升力、布朗力、热泳力等[9]。假设气相占据整个空间，而颗粒被考虑为具有有限质量和动量的点源。当颗粒体积分数 φp较低时，颗粒相与流体间的动量交换可以忽略，称之为单向耦合；增大φp至颗粒与流体间相互作用不可忽略时，两相通过质量（以化学反应的形式）、动量和能量方程双向耦合；更进一步增大 φp，在密相流的欧拉-拉格朗日方法中需要考虑颗粒间碰撞，即四向耦合[10]。

欧拉-拉格朗日方法由于可以更为准确地模拟单个颗粒行为，最近被广泛用于可压缩多相流研究。如Zhang、Dai、Xiao 等采用湍流直接数值模拟求解器[11-13]，结合点颗粒运动模型研究了颗粒的倾向性聚集和扩散。Dai、Xia 等研究了湍流受颗粒调制的规律[14-15]。Li 等[16]在FLUENT 软件求解流场的基础上，开发颗粒求解器代码，研究超重力作用下亚微米颗粒在超声速层流边界层中的运动。随后，其在此基础上进一步研究了超声速绕楔流中的颗粒沉积[9]。Teh 等[17]基于OpenFOAM 模拟流场研究了颗粒对斜激波与层流边界层相互作用的影响。Palmer 等[4]在美国宇航局艾姆斯研究中心开发的计算流体动力学解算器（DPLR，其中包含有限速率反应动力学、热和化学非平衡、精确的高温传输系数和电离流物理的通用模型）的基础上，根据火星的地面和大气观测，使用一组耦合的常微分方程计算沙尘通过激波层的轨迹，研究沙尘颗粒影响下火星再入飞行器隔热层侵蚀。Davuluri 等[7]采用欧拉-拉格朗日方法数值研究烧蚀层裂现象中表面喷射后剥落颗粒的轨迹，该研究使用高超声速KATS-Kentucky 气动热力学和热响应求解器，用于获得流场。Davuluri 等[18]在此基础上进一步讨论了颗粒控制方程中的拖曳力，利用基于雷诺数、马赫数和克努森数的一系列经验关系，建立了混合阻力系数模型。

本文基于国产自主研发的计算流体力学（Computational Fluid Dynamics，CFD）软件PHengLEI（风雷）[19]，开发可压缩边界层颗粒多相流程序，为国家数值风洞软件开发提供技术储备和集成模块。基于典型算例对所开发的模块进行检验，并初步研究了层流边界层中颗粒的运动行为及其对流动性质的影响。

1 数值模型和方法

风雷软件是中国空气动力研究与发展中心研发的面向流体工程的混合CFD 平台。风雷软件的计算范围覆盖低速、亚跨声速和高超声速，是一款具有完全自主知识产权的面向工程应用和学术研究的通用CFD 软件，借鉴了面向对象的大型软件设计理念，采用C++语言编程。风雷软件的整体框架可参考文献[19-20]，这里仅介绍与多相流开发相关的数值模型和方法。

1.1 流体控制方程

本文考虑了颗粒-流体的双向耦合，在求解过程中，通过在流体相控制方程中添加颗粒的动量、能量源项的方式，来体现颗粒对流体相的反作用。实际程序计算中，通过定义欧拉变量，将每一时间步各单元的反馈源项存储到单元中心处。笛卡尔坐标系下，包含颗粒源项的守恒形式流动控制方程为：

其中，ui、 ρ、Tf、E分别表示流体的速度、密度、温度、总能，k为流体传热系数， σij为黏性切应力张量，R为气体常数，计算中取值为 287.974 4 J/(kg·K)。流体计算采用二阶有限体积方法，空间上采用Roe 格式，时间上采用LUSGS 格式。颗粒的动量源项 φui以及能量源项 φei的表达式分别为：

其中，Fp,i是颗粒受到的力，Tp为 颗粒的温度，Qp,i是颗粒释放的热量，Gp,i是 颗粒对流体做的功，dp是颗粒的粒径，Nu是努塞特数，np表示计算网格单元体积内的颗粒数， ΔV表示网格单元的体积，up,i是颗粒速度。

1.2 计算相关参数

本文算例的流体介质为空气，其可压缩计算涉及的黏度由著名的Sutherland 公式给出：

其中，Ts为 萨瑟兰常数，计算中一般取 1 10.4 K，T0、µ0为参考状态的温度与动力学黏度，计算中分别取228.15 K、1.789 4×10−5kg/(m·s)。计算中涉及的无量纲参数：普朗特数Pr=cF,pµk，cF,p为流体的比定压热容；单位雷诺数Re∗∞=ρ∞U∞/µ∞， ρ∞和 µ∞分别为来流密度和黏度；来流马赫数Ma∞=U∞/c∞，c∞为来流声速。

1.3 颗粒控制方程

本文离散颗粒相计算的控制方程分别为颗粒位移、颗粒平动、颗粒转动及颗粒温度控制方程，见式（10～13）：

其中，Xp、vp、 ωp是颗粒的运动学参数，分别是位移、速度、角速度，Tp为 颗粒温度，Ip是颗粒的惯性矩，Fp和Tc分 别是颗粒受到的力及扭矩，cp,p为颗粒相物质的比定压热容。

在颗粒受力方面，Magnaudet 等[21]和Ling 等[22-23]对不可压缩和可压缩流动中颗粒受力及其发展进行了全面概述。然而，Thomas[24]、Tedeschi[25-26]和Hughes[27]等的研究发现，对于粒径在0.1～100.0 μm范围内的颗粒，非定常力并不显著。因此，由于本文研究的颗粒-流体密度比例较大且颗粒尺度较低，非定常力对结果影响较小。这里考虑的颗粒受力模型参考了Li[9]的工作，包括拖曳力FD、 升力FSaff、布朗力FBi、热泳力Ft等。颗粒受力Fp表达式见式（14）：

本文研究的问题中，流动范围涉及亚声速到高超声速、连续流到过渡流区，颗粒从低雷诺数到高雷诺数，因此，所采用的拖曳力FD模型需要适用较大的颗粒雷诺数Rep（以颗粒直径为参考长度的雷诺数）、克努森数Kn和 马赫数Map（ 定义为 |vr|/af,af为颗粒位置处流体的声速，vr为流体与颗粒之间的相对速度）范围。而颗粒常见的Stokes[28]均匀不可压缩蠕流条件下单颗颗粒的阻力公式（FD=−6πµavr，a=dp/2为颗粒半径， µ是流体相介质的动力黏度）已不适用。因此，参考Tedeschi等[29]的工作，本文计算所用的阻力公式（15、16）考虑了可压缩性影响和稀薄效应：

当颗粒位于流场强剪切区域时，会受到Saffman（萨夫曼）升力。Saffman[30]基于无界线性剪切流情况，对单个球体受力问题进行研究，通过理论推导发现了这种垂向流体诱导的剪切升力，这一升力因此被命名为Saffman 升力：

其中，升力系数J=1， ω为颗粒所在位置处流体涡量。

Li 等[9]指出，热泳力Ft在弓形激波附近和温度梯度较高的热边界层中可能很重要, 对于亚微米颗粒，布朗运动及布朗力FBi也可能不可忽视。本文参考Li 等的工作，在颗粒运动模型中也增加了热泳力和布朗力，其中，热泳力采用了Yamamoto 等[31]的无量纲形式的公式：

式中：Aw、A0、Hw及H0均 为克努森数的函数，kˆ 是颗粒和流体的导热系数之比。布朗力可以被考虑为一个高斯白噪声随机过程，其谱强度Snij的表达式为：

式中：Tf是流体温度， δij是 克罗内克符号，kB是玻尔兹曼常数，Cc是斯托克斯-坎宁安修正项。最终，布朗力表示为：

其中， ξi是零均值的单位方差无关的高斯随机数[17]。注意到，在以往的欧拉-拉格朗日模拟中，大部分研究认为拖曳力是占主导地位的，因此，往往在模拟颗粒轨迹时也仅考虑了拖曳力[20]。Li 等[9]在平板边界层层流颗粒轨迹的模拟中发现，热泳力在其模拟的参数范围内（自由来流马赫数小于3.01）是不重要的。同时，他们的模拟结果表明小粒径（dp=0.05µm）颗粒的布朗运动比大粒径（dp=1.0µm）颗粒的布朗运动严重得多，且随着主流马赫数的增大，布朗运动将变得更加剧烈。但无论如何，升力的影响非常重要，当颗粒进入边界层时，考虑萨夫曼升力的颗粒轨迹显著低于未考虑萨夫曼升力的颗粒轨迹。

1.4 流场速度插值方法

Gimenez 等[32]提出了颗粒位置处的插值方法，将流场离散化，并将单元中的流场变量 ϕ存储在单元中心xc处 ，即 ϕc=ϕ(xc)。 考 虑点xc周 围 任 意 位 置 处 ϕ的Taylor 级数展开，并略去高阶项，得到任意位置xp处ϕ(xp)的值：

1.5 颗粒边界条件

所开发代码中，颗粒边界条件包括入口、出口边界条件，周期边界条件和固壁完全反弹边界条件三类。颗粒壁面碰撞和反弹示意图见图1。颗粒的入口条件包括颗粒的位置坐标和速度信息。当颗粒中心超出出口边界或者颗粒中心距出口边界距离小于半径时，从存储颗粒的数组中删除该颗粒信息，如果将颗粒出口速度信息保留并重置其位置坐标，可实现周期边界条件。注意后者仅适用于具有周期边界条件的流场。

图1 计算网格点内颗粒壁面碰撞和反弹示意图Fig.1 Schematic of particle-wall collision and rebound in the calculation cell

颗粒中心与壁面距离小于颗粒半径时，在碰撞点所在壁面网格单元平面上将速度矢量进行反射变换，更新颗粒的位置信息与速度信息，并继承其他物理参数，反射过程无动量损失。

2 颗粒求解可靠性检验

本节将用3 个算例对代码进行验证测试，测试算例分别为：1）泊肃叶流动中颗粒运动算例，用于检验颗粒插值方法；2）超声速绝热平板边界层中单颗颗粒运动算例，用于检验不同受力模型对颗粒运行轨迹的影响；3）等温平板边界层颗粒群对热流的影响算例，验证代码双向耦合部分的准确性。

2.1 泊肃叶流动中颗粒运动

计算二维稳态泊肃叶流动中单个颗粒在流场中的运动速度和轨迹的模拟示意图见图2。图中，T∞为来流总温和等温壁的壁面温度，U∞为来流速度，pin和pout分别为入口和出口总压强。

图2 计算域与边界条件示意图Fig.2 Schematic diagram of computational domain and boundary conditions

根据解析解，该流动中速度均具有抛物线形式解：

其中Ly为半槽道高度。当惯性颗粒（运动方程中仅考虑Stokes 拖曳力）以流向初速度up0=0和垂向初速度vp0≠0在 初始位置处 (xp0,yp0)处释放进入流动中时，其速度 (up(t),vp(t))和 位置坐标 (xp(t),yp(t))的理论 解如下：

计算域总长度设置为 0 mm ≤Lx≤20 mm，高度−0.2 mm ≤Ly≤0.2 mm。在流向上采用均匀网格，垂向网格采用沿壁面指数加密并关于x轴对称。流向200 个网格节点，垂向21 个网格节点。来流温度T∞=334.5 K ， 来流速度U∞=11.7 m/s，来流马赫数Ma∞=0.031 9 ， 入口压强pin=101.625 kPa，普朗特数Pr=0.72 ， 比 热比 γ=1.4 。 颗 粒 密度 ρp=2 000 kg/m3，粒径dp=1×10−5m。 颗粒释放的初始速度为up0=0，vp0=1.0 m/s， 初始位置为xp0=11 mm，yp0=0.01 mm，该位置处流动不受入口影响，充分发展，平均压力梯度为 dp/dx=−10.149 kPa。图3 给出了颗粒轨迹和颗粒速度随时间变化的数值解和解析解的对比，其中图3（a）的纵坐标表示无量纲的颗粒速度，即up/U∞和vp/U∞；图3（b）的纵坐标表示为无量纲的颗粒位置坐标，即xp/L和yp/L，从图中可以看出颗粒解析解和数值解验证良好。

图3 数值模拟的颗粒速度和轨迹与理论结果对比Fig.3 Comparison of particle velocities and trajectories from numerical simulations with theoretical results

2.2 超声速绝热平板边界层中颗粒运动

Li 等[9]对超声速平板层流边界层的亚微米颗粒在不同受力模型下的运动进行了分析，并详细讨论了颗粒初始位置、流动马赫数、颗粒尺寸和密度对于轨迹的影响。本文在与文献[9]相同的条件下模拟颗粒运动，具体流动参数见表1。其中，以单位长度表示的有量纲流体松弛时间为 τ∗f,∞=1/U∞，以单相流边界层流动出口处边界层厚度 δout表示的有量纲流体松弛时间为 τf,out=δout/U∞=δoutτ∗f,∞， 以 δout表示的流体雷诺数为Reout=ρ∞U∞δout/µ∞=δoutRe∗∞。

表1 平板边界层流动参数Table 1 Parameters of the flat-plate boundary layer flow

二维超声速平板层流边界层流动计算域大小为50 mm×50 mm，在x、y方向均采用330 个网格节点，并在近壁处以及平板前缘激波位置处进行加密。与文献[9]不同的是，模拟采用理想气体状态方程，而非Redlich–Kwong 方程。图4 给出了模拟流场速度廓线与Crocco[33]解析解的对比，其中横坐标表示为无量纲高度，纵坐标为无量纲的流体速度，可以看到，数值模拟结果与理论解吻合得很好。

图4 流场速度廓线计算值与解析解对比Fig.4 Comparison between calculated values and analytical solutions of velocity profile

直接数值模拟平板边界层达到稳定后，在入口边界层发展前缘处释放单颗颗粒（见图4 中插图）。释放颗粒的具体参数见表2。其中，颗粒的Stokes 数表示为S t=τp/τf， τp表示颗粒松弛时间。这里以流体来流黏度 µ∞和 δout表示的有量纲颗粒松弛时间尺度为τp,out=ρpdp2/(18µ∞)，对 应 的 颗 粒Stokes 数 表 示 为S tout=τp,out/τf,out，颗粒在入口边界初始释放高度为yp0=0.15 mm，颗粒释放的初始速度为当地流场速度。

表2 模拟的颗粒参数Table 2 Parameters of a single particle released on plate

颗粒受拖曳力、萨夫曼力、热泳力和布朗力的共同作用。模拟得到的颗粒运动轨迹见图5，可以看到，考虑萨夫曼升力时，颗粒运动轨迹呈向下的趋势，这与Li 等[9]的结果一致。注意到，由于气体状态方程的不同，颗粒轨迹会和参考文献存在一定区别，但从定性上分析，颗粒在受到不同力的情况下，总体运动趋势是一致的。

图5 平板边界层中运动颗粒y 坐标随时间变化Fig.5 Temporal variation of coordinate values of moving particles along wall-normal direction in boundary layer flow

同时，这一流动情况下其他作用力对轨迹的影响见图6。可以从图中看出，在来流马赫数为2.05 时，在拖曳力和升力的基础上，进一步考虑热泳力以及布朗力，对粒径为1 um 的颗粒的运动轨迹并无显著影响，因此，热泳力以及布朗力可以忽略。这也与文献[9]中的结论相吻合。

图6 不同受力情况下颗粒轨迹对比Fig.6 Comparison of particle trajectories under different force conditions

2.3 等温平板边界层中颗粒群对热流的影响

Ching 等[35]在采用欧拉-拉格朗日方法进行双向耦合模拟时，与Wang 等[34]进行了对比。在本研究中，参考Ching 等[35]的方法，模拟Wang 等[34]文献中的大滑移区域，并与理论结果进行比较，模拟参数见表3。其中，Ts为萨瑟兰常数，颗粒与流体的质量比β 定 义 为ρ¯d,∞=m¯dnd是来流颗粒相的体密度，nd是 单位体积的颗粒数，m¯d是颗粒的平均质量。为了与理论值进行对比，颗粒仅受斯托克斯阻力作用。二维计算域总流向长度设置为 [−0.1, 1.0]mm，高度[0, 0.2] mm 。 在流向上采用靠近x=0的指数加密网格，在垂向上采用靠近y=0的指数加密网格，网格节点数Nx×Ny=91×56。本算例颗粒初始释放参考Ching 等[35]方式，见图7。初始颗粒从x=0、y=0～δout处均匀释放，假设初始颗粒速度与无限远来流速度相同并按照来流的质量比 β控制入口处释放的颗粒数。每隔 ΔT=2.71×10−9s重复释放；当颗粒超出出口时将颗粒移除。

表3 双向耦合平板边界参数Table 3 Parameters of flow and particles in two-way coupling

图7 颗粒释放示意图Fig.7 Schematic of particle release

图8 单相与两相流中无量纲壁面热通量Fig.8 Non-dimensionl wall heat flux in particle-free and particle-laden flows

3 颗粒与层流边界层之间的相互影响

3.1 颗粒对流体的调制研究

基于以上的测试结果，进一步分析了层流边界层中的颗粒-流场相互作用。由2.2 节可以得知萨夫曼升力在超声速平板边界层的颗粒运动中是不可忽略的（萨夫曼升力对于颗粒的运动轨迹具有一定影响），因此，在本算例中，我们将对2.3 节的双向耦合进行验证，具体流动参数见表3。在验证双向耦合的基础上将流向计算域延长至 2.5 mm，给颗粒附加萨夫曼升力，考虑在更真实颗粒受力情况下的双向耦合调制规律及其沿程信息。

图9 给出了不同情况下模拟得到的无量纲平板壁面热通量沿程变化。可以看出，两相流中颗粒导致壁面热通量显著增加。注意到，图中仅考虑拖曳力的颗粒两相流相比纯气体的无量纲热流已经有显著增大。而相比不考虑萨夫曼升力的理论结果，萨夫曼升力对壁面热流有进一步的增强作用，并且增强效应沿流向不断增长。在x∗=0.35处，不含萨夫曼升力的算例相对于纯流体壁面，壁面热通量增强达到了52.8%，附加萨夫曼升力后，相对于不加萨夫曼升力的情况又增加了 9.3%。

图9 不同情况中无量纲壁面热通量沿平板流向变化Fig.9 Variation of non-dimensional wall heat flux along flow direction under different conditions

接下来通过边界层厚度的变化来分析热流和摩阻变化的原因。图11 给出了不同模拟中的边界层位移厚度，其中，位移厚度 δ′的 无量纲形式为δ∗=δ′/(λ∞Re1∞/2)。从图11 可以看出，由于颗粒以来流速度进入边界层，边界层内的流体通过和颗粒的相互作用从颗粒相获得动量，进而加速，导致边界层厚度减小，从而进一步导致边界层内的流体速度梯度和温度梯度增加，这是导致图9 和图10 含颗粒流相较于纯流体摩阻与热流剧烈增加的主要原因。

图10 不同情况中壁面摩阻系数沿平板流向变化Fig.10 Variation of skin friction coefficient along flow direction under different conditions

图11 不同情况下边界层位移厚度对比Fig.11 Variation of displacement thickness along flow direction under different conditions

3.2 流体对颗粒输运的影响研究

首先，讨论萨夫曼升力对颗粒运动的影响，由式（22）可知，萨夫曼升力的方向主要由颗粒-流体间的滑移速度以及涡量确定。因此，图12 给出了典型单个颗粒在进入边界层后，颗粒位置处流向上的滑移速度 ΔU∗及 流体涡量 ω∗x随 时间的变化（其中 ΔU∗和 ω∗x均用来流速度无量纲化）。可以看到颗粒以初始来流速度进入边界层后，会受到一个向下的萨夫曼升力，这是由于颗粒与流体流向滑移速度为负导致的，如图12（a）。

图12 单个颗粒流向滑移速度和涡量随时间步变化Fig.12 Variation of individual particle slip velocity and vorticity with time steps

其次，为分析流场对颗粒运动的影响，进一步统计了流场中颗粒流向、垂向速度分布u¯p、v¯p和颗粒浓度 ϕp。由于颗粒对流体影响的沿程变化，且在大滑移区含颗粒与不含颗粒情况下的流动性质差别不显著，为简化分析，首先沿流向将计算域划分为四个部分，分 别 为x*= [0, 0.087 5],x*= [0.087 5, 0.175 0],x*=[0.175 0, 0.262 5],x*= [0.262 5, 0.350 0]，对这 四个区域的颗粒分别进行统计平均，u¯p、v¯p和 ϕp的统计结果分别见图13、图14 和图15。图中，纵坐标y均以单相流边界层流动出口处边界层最大厚度 δout进行无量纲化，而速度u¯p和v¯p分 别采用U∞无 量纲化，表示为u¯∗p、v¯∗p。为了更直观地观察颗粒与流体的动量交换情况，将相同流向段流体的速度u¯∗f也同时绘制在图中，并用虚线表示。

图13 不同流向段颗粒和流体平均流向速度分布廓线Fig.13 Profiles of time-averaged streamwise veloicity of particle and fluid in different streamwise segments

图14 不同流向段颗粒和流体平均垂向速度分布廓线Fig.14 Profiles of time-averaged vertical veloicity of particle and fluid in different streamwise segments

图15 不同流向段颗粒浓度沿高度分布Fig.15 Distribution of particle concentration along wall-normal direction in different streamwise segments

图13 给出了不同流向段颗粒流向速度分布廓线。从图中可以观察到，颗粒速度沿着流向不断衰减，并且越靠近壁面，衰减程度越大，这一速度变化和流场的速度廓线相类似。在计算域的四个流向段中，颗粒速度均大于流体速度，因此当颗粒进入边界层时，由于速度高于流体速度，会产生一个方向向下的萨夫曼升力，可能会使颗粒向壁面运动，这一点也与图12 观察到的结果一致。

同时，流体的垂向速度也可能对颗粒的垂向输运造成影响。图14 给出了不同流向段颗粒相和流体相的平均垂向速度分布廓线。可以看到，在靠近边界层前缘段（第一段），由于存在激波，流场垂向速度沿高度增加；在第二段，由于远离了激波区域，流场的垂向速度骤减并且随着沿程逐渐降低。根据图14 流场的垂向沿程速度分布，可以观察到流场的垂向速度对颗粒垂向速度起着加速作用。同时，由于在近壁处颗粒会受到更强的升力作用，颗粒在近壁区域的垂向速度会存在负值。

对图15 中颗粒沿程浓度进行分析可以发现，在近壁区域，颗粒会存在壁面的累积效应，这是由于近壁处边界层内流体速度梯度使得颗粒产生靠近壁面的垂向速度导致的（见图14），并且这种累积效应会沿流向逐渐加剧（见图15）。在远离壁面区域（ y ＞0.5δout）可以观察到，颗粒垂向速度场相比于近壁区域（y＜0.5δout）沿高度并没有明显的变化，因此在垂向统计域内不会出现类似于近壁区域的明显聚集现象。其次，分析颗粒沿流向的分布，可以观察到，颗粒平均浓度沿流向逐渐增加，这是由于在颗粒进入边界层后，边界层内的低速流体在流向上不断将颗粒减速，导致沿程颗粒平均流向速度降低（见图13），从而造成沿程颗粒沿流向累积。

4 结 论

基于PHengLEI 软件开源框架，采用紧耦合的模式开发了点颗粒模拟求解器。文中介绍了点颗粒模型的具体求解方法，包括颗粒控制方程和受力模型、颗粒追踪算法和边界条件、颗粒位置处流动信息插值等。分别采用泊肃叶流动中颗粒运动轨迹和速度、超声速绝热平板边界层中颗粒的运动和颗粒对等温平板边界层热流的影响三个算例，对所设计的求解器进行验证。计算表明，颗粒求解器对可压缩流动中颗粒的运动轨迹、双向耦合两相流中壁面热流的预测与文献中的计算和理论结果吻合，证明了求解器的可靠性以及对于可压缩两相流的模拟能力。

在同时考虑拖曳力和升力的情况下，双向耦合模拟了马赫数为1.5 的等温无滑移壁面上的层流边界层两相流。结果表明，颗粒和边界层流动存在显著的相互作用，具体表现为当颗粒以来流速度和温度进入边界层时，通过和流体的能量交换，颗粒的存在会增强边界层的壁面热通量，通过和边界层的动量交换，极大提高边界层的壁面摩阻。颗粒的存在还会减少边界层的位移厚度。同时，超声速层流边界层流场沿流向变化也会影响颗粒沿程分布。边界层内颗粒速度大于流体速度而产生向下的剪切升力，使颗粒垂向速度为负，向壁面堆积，这种堆积效应沿流向不断增强，导致近壁颗粒浓度增加。

后续，将进一步持续开展颗粒求解器功能的扩充、开发和优化。同时，基于求解器研究可压缩边界层两相流的气动特性和颗粒相运动、分布规律。讨论在不同的流动参数（马赫数、雷诺数、普朗特数、壁温等）和颗粒参数（体积分数、混合粒径、密度比、释放方式、重力等）情况下颗粒-流体的相互作用，为高超声速飞行器烧蚀颗粒边界层特性等问题提供相应参考。