改进的欧拉-拉格朗日方法预测爆震波后液滴生存距离

《Chinese Journal of Aeronautics》：An improved Eulerian-Lagrangian method to predict the droplet survival distance behind detonation waves

　　

由于简单的几何结构和高体积能量密度，液体燃料爆震发动机引起了广泛关注。然而，其数值模拟仍然具有挑战性，因为喷雾雾化、液滴破碎、湍流混合和爆震波传播等多个过程同时存在并相互作用。特别是在欧拉-拉格朗日（Eulerian-Lagrangian）框架内，现有破碎模型无法准确预测爆震波后的液滴生存距离，这可能导致爆震波传播速度被高估。为此，研究人员首先对液滴破碎过程进行了理论分析。研究发现，由于子液滴直径较小且周围流速较高，液滴破碎时间与松弛时间之比可趋近于1。这表明子液滴的破碎行为可能与母液滴不同。基于理论分析，研究人员提出了一种改进方法，当破碎时间与松弛时间之比超过临界阈值时，抑制破碎过程。随后，研究人员将所提方法纳入OpenFOAM?框架内的Pilch-Erdman和Reitz-Diwakar模型中。爆震波作用下液滴阵列破碎的数值模拟表明，改进方法在不同条件下均具有更好的预测性能。具体而言，对于初始煤油液滴直径为40.7 μm的Reitz-Diwakar模型，液滴生存距离从6.11 mm提升至12.56 mm（实验值：12.73 mm），爆震波传播速度从2153 m/s降至1948 m/s（实验值：1897 m/s）。Pilch-Erdman模型也获得了类似结果。总体而言，所提出的模型有效、稳健，且易于集成到现有框架中。

随着爆震燃烧极高的放热率，爆震被广泛应用于旋转爆震发动机（Rotating Detonation Engines, RDEs）、脉冲爆震发动机（Pulse Detonation Engines, PDEs）、斜爆震发动机（Oblique Detonation Engines, ODEs）以及基于爆震的涡轮发动机（Detonation-Based Turbine Engines, DBTEs）等新兴推进装置中。液体燃料爆震推进系统因其高体积能量密度而受到广泛关注，被认为是下一代航空航天动力系统极具前景的技术。然而，液体喷雾与高速爆震波（传播速度超过1500 m/s）之间的复杂相互作用，使得液体燃料爆震的稳定传播比气体燃料爆震更具挑战性，尤其是爆震波作用下的液滴破碎过程尚未得到充分研究。

大量实验证据一致表明，快速液滴破碎对于液体燃料爆震的成功传播至关重要。Schroeder等人利用激波管研究了RP-2液滴在超声速激波下的变形和剪切剥离破碎机理，指出快速液滴破碎是液体燃料爆震传播的关键物理过程。过去的几十年中，大量实验研究致力于解决两相爆震行为和激波-液滴相互作用问题。经典的喷雾爆震研究表明，液体燃料混合物的爆震速度亏损与液滴尺寸密切相关，较大的液滴导致较长的破碎时间，从而延长反应区。近期利用高速成像、平面激光诱导荧光（Planar Laser Induced Fluorescence, PLIF）和粒子散射技术的研究，绘制了液滴在激波和爆震波作用下的详细变形模式、不稳定性增长和破碎机理。重要的是，Brown、Young、Agee等人对液体燃料爆震的实验揭示，液滴在波后区域持续存在，该区域被视为两相爆震的显著特征。该区域的轴向范围被定义为液滴生存距离，其大小被认为与液滴破碎过程密切相关。这些实验发现强调了液滴破碎在爆震波传播中的重要作用，并提供了评估液滴破碎效果的定量指标，即液滴生存距离。

实验仍是研究爆震发动机的重要方法，但随着计算技术的快速发展，高保真数值模拟使研究人员能够获取爆震发动机内部的详细流场信息，从而更全面地理解相关流动结构和参数影响，加速基于爆震的推进系统研制。然而，当前液体燃料爆震的数值模拟仍面临显著挑战。常见的简化是假设均匀液滴直径而忽略破碎过程。但实验证据表明，即使存在相对较大的液滴，爆震波也能自持传播；相反，在相同液滴尺寸条件下，忽略破碎的模拟往往预测爆震解耦，显示出数值与实验结果之间的根本不一致。为弥补这一差异，传统破碎模型常被引入，但这些模型倾向于预测过快液滴破碎，导致爆震传播速度显著高于实验测量值。更重要的是，近期研究表明，现有欧拉-拉格朗日破碎模型严重低估生存距离。例如，初始煤油液滴直径为40.7 μm时，Pilch-Erdman模型偏差达52.1%，Reitz-Diwakar模型达53.7%，WERT49模型达70.3%，KH-RT模型高达92.7%。这些发现共同表明，现有破碎模型不能满足液体燃料爆震模拟的精度要求。Wang等人提出了一个动态模型以反映两相爆震流中液滴破碎过程的物理图像，但其预测的液滴生存距离未经过实验验证。Musick等人利用精确热力学计算和全局动力学机制模拟了气液烃类爆震，认为破碎过程将决定波后条件下液滴的寿命，但结果无法直接与现有实验比较。

为理解现有破碎模型在爆震波条件下为何不准确，研究人员分析了此类极端环境中控制液滴行为的关键物理因素。首先，传统破碎模型主要针对内燃机和燃气轮机条件开发和验证，其气液相对速度通常在约100 m/s量级；而爆震发动机中液滴承受的载荷远为剧烈，相对速度超过约1500 m/s，显著改变了破碎动力学。其次，爆震波产生的强冲击载荷会产生大量小液滴。这些在高速相对速度下的小液滴（直径<20 μm）会显著加速并紧跟流动，其行为以斯托克斯数（Stokes number, S_t = τ_r/τ_f ? 1）表征，其中τ_f为流动特征时间尺度，τ_r为液滴松弛时间（液滴加速到终端速度所需时间）。较小的斯托克斯数对应更快的液滴加速和更低的破碎敏感性，凸显了松弛时间τ_r在液滴动力学中的关键作用。已有研究注意到液滴松弛时间可能影响液滴的变形与破碎，但主要关注低速场景，如重力下的自由下落。研究人员因此假设，在爆震流特有的大速度差条件下，小液滴的快速加速可能显著影响其破碎速率，但目前尚无研究直接在爆震环境中解决这一问题。

针对现有破碎模型无法准确预测爆震波后液滴生存距离且模拟爆震波传播速度快于实验的问题，研究人员进行了理论分析以研究破碎与小液滴运动之间的关系，进而识别传统破碎模型的局限性，然后提出了适用于爆震波环境的改进液滴破碎方法，并通过实验进行了全面验证。该工作扩展了现有破碎模型的适用性，增进了对爆震波相互作用下液滴破碎行为的理解，对实现未来高性能飞机和航天器具有重要作用。

在计算方法上，研究人员基于OpenFOAM V8?开发了专用求解器LagrDEFoam，采用颗粒-网格双向耦合的欧拉-拉格朗日方法，使用有限体积法离散控制方程，求解质量、动量、能量守恒方程、组分输运方程和理想气体状态方程。液滴运动由相对速度决定的气动力控制，采用Spalding模型和Ranz-Marshall关联式描述液滴蒸发和传热。时间积分采用二阶后向格式，时间步长5 ns；对流通量采用二阶Kurganov型中心差分格式和van Leer限制器；扩散通量采用二阶中心差分格式。研究中修改了Pilch-Erdman和Reitz-Diwakar破碎模型，以Reitz-Diwakar模型为例详细说明了修改流程。计算设置方面，模拟基于实验条件设置计算域长400 mm、宽50 mm，左侧0-300 mm为气态煤油区域以形成稳定爆震波，右侧300-400 mm为核心区域，爆震波进入液态燃料区传播。共设置五个基准工况（C#1-C#5）覆盖小液滴（5-20 μm）和大液滴（40.7 μm）的液体燃料爆震，并针对正辛烷和正十二烷燃料扩展验证工况C#6-C#8。所有工况初始压力1 atm，初始温度300 K。网格无关性验证表明50 μm和100 μm均匀网格结果几乎一致，200 μm网格存在较大偏差，综合考虑计算成本选用100 μm均匀网格。化学机理验证采用10组分/12步机理作为基准，对比6组分/2步简化煤油机理，二者预测的平均爆震胞格尺寸分别为10.55 mm和10.01 mm，偏差较小，故选用6组分/2步机理；正辛烷采用两步机理，正十二烷采用单步机理。

在结果与讨论部分，研究人员首先进行了理论分析。液滴破碎受相对速度U_rel影响，强气动扰动如爆震波使液滴倾向于剥离破碎，较大的U_rel导致更快的破碎过程。然而，气动力使液滴沿气流方向加速，加速所需时间称为松弛时间τ_r。液滴加速在破碎时间τ_b内增大其速度U_d，从而减小U_rel。通过数值模拟获取波后气相速度、密度和液滴速度等合理参数进行分析，采用代表性液滴初始速度200 m/s，液滴直径0-5 μm，煤油密度738.62 kg/m³，波后气相密度3 kg/m³。基于Reitz-Diwakar模型的计算表明，小液滴松弛时间极短，加速强烈，速度快速上升使相对速度迅速降低；直径小于0.5 μm的液滴加速后相对速度低于临界破碎速度，破碎无法持续而提前终止。通过定义f = τ_b/τ_r，参数f代表液滴在破碎时间内的累积增速能力，理论分析得出当液滴直径小于0.5 μm时f小于0.6，因而f取值约为0.6附近，实际动态过程应视为一个范围。

基于理论分析，研究人员提出了改进破碎方法，明确考虑破碎期间液滴加速的影响。核心思想是根据松弛时间τ_r和破碎时间τ_b的相对大小判断液滴运动是否会在破碎过程中显著降低气液相对速度并提前抑制破碎。当τ_b?τ_r时，液滴加速太慢无法显著改变U_rel，破碎可维持；当τ_b与τ_r相当或更大时，破碎期间加速会快速减小U_rel，使破碎标准在完成前不再满足，破碎可能提前终止。因此，改进模型仅在τ_b/τ_r < f（f取0.4）时允许破碎发生，否则禁用破碎模型。研究人员将Pilch-Erdman和Reitz-Diwakar模型分别修改为YPilch-Erdman和YReitz-Diwakar模型。

在液滴生存距离定义方面，将计算域沿x方向划分均匀1 mm间隔，统计每个间隔内的总液滴质量，将原始液体质量蒸发99.5%的特定距离定义为液滴生存距离，同时计算每个间隔的几何平均直径D10。在C#5工况下，原始Pilch-Erdman和Reitz-Diwakar模型结果差异很小，改进模型使液滴生存距离增至12.56 mm，与实验值12.73 mm高度吻合。液滴平均直径沿爆震波传播反方向减小，且靠近爆震锋面处下降更陡。根据破碎速率将整个过程分为快速破碎阶段（312-318 mm）和较慢破碎阶段（305-312 mm）。在快速破碎阶段，大液滴受爆震波载荷产生显著变形和破碎，改进模型对大液滴正常破碎几乎没有限制，四种模型结果无明显差异，证实了改进的合理性；在慢速破碎阶段，液滴直径和质量同时减小，改进模型限制小液滴破碎。尽管YReitz-Diwakar给出的波后液滴直径大于YPilch-Erdman，但两种改进模型获得的生存距离基本相同。

研究人员还考察了f值对生存距离的影响。理论分析基于恒定特征速度得出f约为0.6，但真实爆震条件下液滴速度和气体速度在空间和时间上变化（0-1200 m/s），因此基于实验与模拟对比进一步调整f值。通过比较C#5条件下不同f值的模拟结果，当f设为0.4时与实验观测最一致。

在爆震波与液滴相互作用及模型验证方面，无破碎模型模拟显示激波与化学反应区解耦，而加入改进破碎模型后激波稳定传播并与反应区强耦合，证明破碎模型对于液体燃料爆震模拟的必要性。针对小液滴工况（C#1，5 μm），原始破碎模型严重低估生存距离，改进模型预测值与实验吻合良好，且对不同燃料均表现稳健。针对大初始液滴直径工况，煤油40.7 μm和正十二烷41 μm，原始Reitz-Diwakar模型预测生存距离仅6.11 mm和6.93 mm，与实验值12.73 mm和11.20 mm偏差显著；改进YReitz-Diwakar模型预测为12.56 mm和11.69 mm，相对偏差均降至5%以内。图12显示液滴直径、速度和压力波相对位置的关系：波后液滴直径逐渐减小，随着液滴变小，载气与液滴之间相对速度降低，破碎速率减慢，验证了理论分析的合理性。进一步在C#2-C#4工况（初始液滴直径8 μm、12 μm、20 μm）下验证，实验生存距离分别为（3.19 ± 0.38）mm、（8.52 ± 0.30）mm和（15.72 ± 0.30）mm，对应爆震速度为（1748 ± 46）m/s、（1666 ± 54）m/s和（1646 ± 75）m/s；YReitz-Diwakar预测的生存距离和爆震速度均与实验高度一致，而改进前模拟值偏差较大。图15显示原始模型预测液滴即时破碎汽化，波后几乎无液相，与实验观察的液滴云区不符；改进模型显著延长生存距离，且存在明确的尺寸依赖趋势：初始直径越大，液滴生存距离越长。

改进模型还揭示液滴破碎行为与局部波结构强烈相关（C#2-C#4）。在马赫杆（Mach Stem, MS）后方，激波强度和波后压力最高，液滴因强烈动量交换发生均匀而快速的空气动力破碎；在三波点（Triple Points, TP）附近，破碎最剧烈，横波与剪切层复杂相互作用产生高相对速度，导致致密细雾云的灾难性破碎；在入射波（Incident Wave, IW）下游，激波强度相对较弱，破碎过程也较弱，较大液滴在完全破碎前持续更久。韦伯数（Weber number, We）和欧内索格数（Ohnesorge number, Oh）分布显示波后液滴Oh数极小（Oh < 0.1），内部粘性力影响可忽略，破碎主要由惯性力与表面张力竞争控制，即由韦伯数主导。Reitz-Diwakar模型在We≈6-200内保持高精度，随We接近1000捕捉灾难性破碎能力逐渐下降但仍可接受；本研究中波后大部分液滴We在0-9000范围，通过成功复现不同燃料和工况下的实验生存距离，证明了该范围内数值模拟具有高可信度。

与气相爆震不同，液体燃料爆震具有更分散的放热区和更长的生存距离。研究人员通过计算放热率（Heat Release Rate, HRR）超过阈值2×10¹¹ J/(m³·s)的总面积A_s定量评估反应区尺寸。液相爆震的反应区面积达113.3 mm²，约为气相值36.9 mm²的3倍，这种显著扩张归因于两相流的异质性，离散液滴燃烧形成分布式反应区而非气相爆震典型薄而紧凑的放热层。关于破碎模型的影响，改进YReitz-Diwakar模型预测的高HRR面积始终大于原始模型，例如初始直径20 μm时，改进模型面积为113.3 mm²，较原始模型预测的69.0 mm²增加64.2%。原始模型低估液滴生存距离导致快速汽化，而改进模型捕捉到延迟破碎过程，使液滴更深地穿透波后区域，从而增加强反应区的有效体积。

最后，研究人员详细讨论了爆震波载荷下的液滴破碎过程。小液滴受爆震波作用后，气体速度U_g急剧增大而液滴速度U_d近似为零，产生超过临界阈值U_rel,crit的大相对速度并启动变形。但满足破碎标准并不意味着瞬时完全破碎，破碎需要有限发展时间，变形必须持续到破碎产物从母液滴分离。由于小液滴对气动载荷响应迅速，快速加速使U_rel下降，一旦U_rel降至U_rel,crit以下，气动载荷不足以维持破碎，液滴仅发生变形而未能完成破碎；随U_rel进一步降低，表面张力占优，可能出现部分形状恢复。已有实验报道高马赫数条件下小液滴特别容易发生此类不完全破碎。大液滴虽然也受波后气流加速，但较大惯性导致响应较慢，U_rel维持在U_rel,crit以上的时间更长，从而维持大尺度变形并推进破碎，最终从母液滴剥离小液滴。

在结论部分，研究人员指出，液滴破碎是爆震燃烧数值模拟中的重要过程。传统破碎模型无法准确预测液滴生存距离，且模拟的爆震波传播速度快于实验结果。通过引入破碎时间和松弛时间，研究人员建立了改进破碎模型，有效解决了上述问题。主要结论如下：（1）爆震波与不同直径液滴相互作用产生不同破碎行为。大液滴松弛时间长、加速缓慢，破碎时间远短于松弛时间，破碎期间始终满足破碎标准，运动加速对破碎影响很小；小液滴松弛时间极短、加速快，当松弛时间等于或小于破碎时间时，液滴可能没有足够时间完成破碎，液滴加速降低液滴与载气间的相对速度从而抑制破碎，某些条件下液滴甚至可能变形后恢复原状而不解体。（2）在欧拉-拉格朗日框架内发展了一种适用于爆震诱导破碎的改进液滴破碎方法，修改了两种现有破碎模型以更好地捕捉爆震波载荷下的液滴行为。基于改进模型的模拟与实验结果高度一致，爆震波传播速度偏差低于8%，液滴生存距离偏差约8.4%。该方法不仅限于爆震或激波下的液滴破碎，还可适用于涉及小液滴尺寸和大相对速度的广泛场景。重要的是，该修改保留了原始破碎模型的适用性和通用性。此外，该方法可类似应用于其他破碎模型，如TAB模型和KH-RT模型，不仅适用于欧拉-拉格朗日框架，也可用于界面模拟，相关研究将在下一步开展。研究发现小液滴破碎与液滴运动强耦合，有望启发更精细破碎模型的发展。