可压缩流大涡模拟耦合概率密度函数方法中的滤波压力模型.pdf

第 卷 第 期国防科技大学学报 年 月 ：可压缩流大涡模拟耦合概率密度函数方法中的滤波压力模型张林，梁剑寒，王翼，孙明波（国防科技大学 空天科学学院，湖南 长沙 ）摘要：亚格子组分 温度关联项显著影响反应流大涡模拟精度。利用概率密度函数方法中概率等效的特点，发展一种新的滤波压力模型，可以良好封闭亚格子组分 温度关联项。介绍概率密度函数及其耦合求解方法，在已有模型基础上推导建立新的滤波压力模型，并在三维超声速氢气 空气时间发展反应混合层中对不同的滤波压力模型进行数值测试。结果表明，与传统的滤波压力模型相比，新的滤波压力模型可以明显改善反应混合层的模拟准确度。特别地，基于新的滤波压力模型，大涡模拟耦合概率密度函数方法可以较好地模拟链式听语音聊科研与作者互动反应中间微小组分如超氧化氢基等，有望更有效地再现自点火等复杂现象。关键词：可压缩流；大涡模拟；概率密度函数方法；滤波压力；亚格子组分 温度关联项中图分类号：文献标志码：开放科学（资源服务）标识码（）：文章编号：（），（，）：，：；以超燃冲压发动机为动力的吸气式高超声速飞行器同时在军用和民用方面展现巨大潜力，多年来受到广泛关注 。当前，超燃冲压发动机正在向宽域、高性能、强鲁棒性的方向发展，其中的宽域湍流燃烧过程极为复杂，存在来流条件变化大、燃烧模式复杂、稳焰方式多样等特点，湍流 化学反应相互作用强烈，需持续开展探索。随着计算流体力学和计算燃烧学的发展，数值模拟逐渐成为理解认识宽域湍流燃烧过程不可或缺的 重 要 手 段。当 前，大 涡 模 拟（，）已广泛用于高速湍流燃烧模拟 。虽取得了巨大成功，但仍然有一些问题尚待明确，特别是部分亚格子模型的封闭。其中，最为重要的模型之一就是湍流燃烧模型，通常用于封闭滤波的化学反应源项。概率密度函数（，）方法的基础为随机假设 ，即假设湍流流场中的物理量为随机变量。从随机方法出发，无论过程多么复杂，所有的单点单时刻过程均严格封闭。因此，湍流 化学反应相互作用（即滤波的化学反应源项）可以由 方法精确封闭。为了提高 方法的鲁棒性，通常将其与 耦合求解，称之为 方法，该方法在湍流反应收稿日期：基金项目：国家自然科学基金资助项目（）作者简介：张林（），男，陕西富平人，副研究员，博士，：国 防 科 技 大 学 学 报第 卷流中得到广泛应用 。方法发展于低速流，为了将其拓展至高速流，众多学者做出了大量杰出工作。等 基于焓的控制方程，在 能量方程中引入了平均压力的物质导数项，使得 可以考虑高速流中的可压缩效应。然而，该方法并未考虑黏性耗散。为了处理超声速流场中广泛存在的间断，等 在 能量方程的计算中引入了 关系式。在 等 的工作中，可压缩效应和黏性耗散效应基于 场计算，然后通过插值的形式反馈给 粒子。后来，学者们提出了能量 压力 速度 标量联合 以及速度 标量 压力 方法。尽管理论上更加完备，但却需要封闭额外的条件滤波项，这些条件滤波项封闭极其困难，而已有的模型大多缺乏充分的验证。在传统 中，滤波压力通常基于理想气体状态方程，由滤波的组分、温度、密度直接得到，并未考虑亚格子组分 温度关联项（本文也称为亚格子压力）。由于该项建模极为困难，已有大涡模拟中大都忽略了它。研究 表明，忽略亚格子组分 温度关联项在燃烧释热较强时将带来显著误差。而对于 方法，由于 粒子信息更加全面，有望更好封闭滤波压力。提出了一种新方法，平均压力通过名义粒子上的瞬时参数值直接计算。考虑到 系综平均场包含统计噪声，直接反馈给 可能引入数值误差，甚至引起数值不稳定，因此采用移动时间平均方法以减小统计噪声，但该处理仅适用于统计定常流 。在 方法用于非定常高速流动燃烧问题的模拟时，滤波压力模型对于该方法的数值解质量、湍流 化学反应相互作用过程的模拟精度、耦合方法鲁棒性等具有显著影响，因此迫切需要发展高分辨率的滤波压力模型。为了适用于非定常强释热条件下的可压缩流模拟，基于 方法，提出了一种新的滤波压力模型，充分利用 和 各自优势，并在考虑亚格子组分 温度关联项的同时降低统计噪声的影响。在三维超声速氢气 空气时间发展反应混合层中进行了数值测试，并与直接数值模拟（，）结果进行对比，验证新模型的优势。方法 标量 输运方程为了改善间断流场中的 能量精度及 耦合求解精度，采用文献 的一致守恒型 方法。对于气相多组分可压缩流，混合物总能的能量方程可以表述为 （）式（）满足爱因斯坦求和约定（单个式子中成对出现的下标为哑标，表示求和，与具体的哑标字母无关，如 ）。其中，为密度，为速度分量，为黏性应力张量分量，为压力，为时间，、均为坐标分量，为混合物的总能，为能量通量项，由式（）确定。（）其中，为导热系数，为温度，为混合物包含的组分数模，、和分别为组分 的焓、分子扩散系数和质量分数。经过简单变化，混合物总能的能量方程可以改写为 ()()()（）式中：为分子黏性系数；为普朗特数，即分子黏性系数与热扩散系数之比；为气体常数。上述能量方程和各组分的质量方程可以写成统一的形式，称为标量（含组分质量分数和总能）守恒方程，即 （）其中，表示标量，和 分别表示标量通量和标量源项，分别由式（）（）计算。（）（）（），（）其中：（）为组分施密特数，即黏性系数与分子扩散系数之比，下标中包含括号与不含括号的意义相同，旨在区分该下标不进行爱因斯坦求和；表示组分质量生成率；表示燃烧释热，表示可压缩效应和黏性耗散效应源项，其计算公式为 ()()（）基于标量守恒方程（），结合 的基本性第 期张林，等：可压缩流大涡模拟耦合概率密度函数方法中的滤波压力模型质 ，可以推导得到精确的标量 输运方程为 珘 （）()（珘）（）其中，表示标量 ，表示标量，表示对应的状态空间变量，表示组成的状态空间矢量，“”表示密度加权空间滤波，表示条件滤波。参考常用的封闭模型 ，模化的标量 输运方程为 珘 珋（珚）（珋）（）（，珋，珘）（）其中，变量上方的“”表示空间滤波，为平均分析扩散系数，为湍流扩散系数，为混合频率，能量方程中可压缩效应和黏性耗散效应的影响隐含在标量源项 中。耦合求解 方程的模化过程和粒子求解格式等详细信息，可以参考文献 ，。简单起见，这里不再赘述。由于粒子方法为随机方法，包含统计噪声，在求解过程中将随着时间的累积带来较大误差，甚至引起计算发散。为了提高数值解质量，通常将拉格朗日粒子 方法与欧拉网格 方法相结合，采用耦合方法求解。一般认为，耦合方法及其求解过程对数值结果有重要影响，其核心在于充分利用粒子方法对单点非线性过程高精度模拟的优势，同时尽可能规避或减小随机方法统计噪声带来的误差。在不可压缩流中，为了解决 直接反馈粒子质量密度带来的数值不稳定问题，一般通过额外求解变量（如比体积、焓、等效温度等）的输运方程，此时粒子的影响以该输运方程源项的形式来体现。本文在此基础上，发展了适用于可压缩流的 耦合方法：求解滤波的质量、动量、能量以及组分输运方程，而 求解标量（即组分和能量）的输运、混合和反应。两种求解器冗余求解组分和能量，该冗余量一方面可以评估 与 的一致性，也可用于改善数值解质量。将化学反应源项信息传递给 ，则将滤波速度、密度、扩散系数、可压缩效应和黏性耗散效应等信息传递给 用于粒子迭代求解。由于 反馈给 的唯一纽带为组分源项，该量为单点过程，不包含时空导数信息，同时 有限差分求解器对于反应源项具备较好的适应性，因此该求解方法可以改善 耦合计算的鲁棒性。详细的迭代耦合求解过程如图 所示。前期研究结果显示，该耦合方法具备良好的数值稳定性。图 耦合求解示意图 国 防 科 技 大 学 学 报第 卷 滤波压力模型 滤波压力模型多组分理想气体的压力为 （）式中，表示通用气体常数，表示组分的分子质量。对式（）进行空间滤波处理，可以得到珋 珋 槇()（）式中，括号中的第二项称之为亚格子组分 温度关联项。一般认为，当燃烧释热较弱时，该项的影响较小，可以忽略；当燃烧释热较强时，其影响将变得显著，忽略该项将带来明显误差 。但由于封闭困难，通常 中均对其进行了忽略，即模化为珋 珋 （）其中，称之为 可解压力。为了方便对比，本文将该模型称之为压力模型 。滤波压力模型对于 拉 格 朗 日 粒 子 方 法，近 期 提出采用如下模型计算滤波压力珋 珋 （）（）（）()（）珋 （）式中，珋 表示 滤波压力，上标（）表示系综域内第 个粒子的参数，表示粒子质量。类似地，本文将该模型称之为压力模型 。由于实际计算中粒子数量有限，滤波压力通过拉格朗日随机粒子求解时包含统计噪声，引入了随机误差。若该结果直接反馈给 进行迭代计算，可能会降低数值解质量，甚至可能导致计算发散。滤波压力新模型考虑 可解压力，定义为 珋 （）其中，上标 表示 进行系综统计得到的结果，即 （）（）（）（）将式（）代入方程（），则 平均压力可以改写为珋 珋 槇（）考虑到 可以有效地完整表征单点物理过程，即 系综得到的 滤波压力和 可解压力的相互关系与流体物理本身是相符的。因此，假设 部分和 部分的滤波压力和可解压力之间保持尺度相似，即珋 珋 （）因此，方法中的滤波压力可以采用 滤波压力来模化，即珋 珋 珋 （）式（）即所提出的 滤波压力新模型，为方便对比称之为压力模型 。可以看出，该滤波压力模型包含了传统的 可解压力。借助 可解压力的数值稳定性，可以有效改善数值计算的鲁棒性。此外，还考虑了 滤波压力中包含的亚格子组分 温度关联项信息，当燃烧释热较强时仍具有良好分辨率。算例验证 算例及计算条件设置为了验证本文提出的滤波压力模型，采用超声速氢气 空气反应时间发展混合层算例进行测试。该算例构型及边界条件相对简单，计算量小，又包含诸如激波、大梯度、湍流 化学反应相互作用等超声速燃烧流场典型特征，被广泛用于 方法和模型的验证。时间混合层包含上下两组平行流相向运动，如图 所示。初始时刻的速度满足双曲正切剖面。给定初始条件，在周期边界条件的作用下，混合层会逐渐失稳转捩，逐步发展为充分发展的湍流。为加速转捩过程的进行，通常需要在平均速度剖面上添加一个小的初始扰动。本文采用给定湍流能谱的带宽扰动方式，采用各向同性湍流能谱生成初始扰动速度信息。该初始能谱具有如下形式（）()()（）其中，为波数，为能谱密度最大时对应的波数。式（）为一维能谱，计算中采用的三维能谱密度为第 期张林，等：可压缩流大涡模拟耦合概率密度函数方法中的滤波压力模型（，）（）（）其中，槡。三维能谱向初始速度场的转换可以参考文献 。图 时间混合层及边界条件示意图 初始的湍流强度取自由来流速度即参考速度的 。这种弱可压缩扰动可以有效抑制可压缩转捩，降低计算耗费。的取值满足流向计算域的网格，具有约 个峰值波长，于是 ()（）其中，表示流向计算网格点数，（）表示取整函数。为了将速度脉动限定在混合层中心附近的有限区域内，在基于式（）获得的扰动速度再乘以一个沿着 方向的指数衰减函数获得最终的扰动速度。该衰减函数表达式为 （）（）式中，为扰动的限制宽度。为了模拟超声速氢气 空气化学反应，采用 等 的组分 步简化反应机理。时间混合层计算参数的设置如表 所示。为保持文章简洁并突出重点，其他计算参数和计算条件的详细设置（如 数值格式、亚格子模型、相关模型及求解格式等），可以参考文献 。为了进行对比验证，采用 、基于均匀搅拌反应器（，）模型的大涡模拟（）以及具有不同滤波压力模型的 分别对该算例进行计算，并对结果进行对比分析。结果与讨论为验证本文的数值格式，采用 模拟了 等 的氮气 氧气超声速惰性时间混合层，对流马赫数为 。需要强调的是，等 的计算设置并不详尽，部分详细参数尚未可表 超声速反应混合层初始化参数 参数燃料流（上层流体）氧化剂流（下层流体）压力 温度 密度（）速度（）声速（）马赫数 组分质量分数 图 无反应时间混合层动量厚度变化对比 知。尽管如此，本文计算所得的混合层动量厚度随时间的变化与 等的结果很吻合，如图 所示。这表明当前方法及数值格式对超声速流场捕捉良好，结果具备可信度，可以作为基准结果评估 的滤波压力模型。需要强调的是，以下所有结果均为无量纲结果。本文无量纲处理中，长度、密度、速度、温度的单位为基准单位，其他单位由以上基准单位组合得到。除长度单位采用 进行无量纲处理外，其他所国 防 科 技 大 学 学 报第 卷有参数均采用氧化剂流（即下层流体）的参数及其组合结果进行无量纲化处理。需要强调的是，在 方法用于超声速反应流时间混合层计算时，对于滤波压力模型 ，由于直接反馈的 滤波压力包含统计噪声，同时本文实际计算中选取的粒子数目有限（每个网格 个粒子），统计噪声在 方程的求解中不断累积，最终造成计算发散。要提高该模型的鲁棒性，可以采用更多的粒子，但粒子数目的增加会显著提高该方法的计算量；同时，根据统计理论 ，有限样本的统计值相较于真实值的误差与样本数量的二分之一次方成反比，因此提高样本数量即粒子数目的方法成效并不显著。兼顾以上两点，本文不再讨论更多的粒子数目下该模型的鲁棒性。此外，参考 的讨论，此时 滤波压力由粒子瞬时参数经系综平均获得，需要配合采用移动时间平均方法等以减小统计噪声，但目前这些处理仅适用于统计定常流，而发展适用于非定常流的降噪方法并非本文的主要关切。鉴于以上两方面的考虑，在后面结果的展示中，不再展示压力模型 的结果。不同模型得到的混合层动量厚度随时间动态增长情况如图 所示。其中，表示基于 数据进行滤波处理得到的结果，表示采用基于压力模型 的 方法计算结果，表示采用基于压力模型 的 方法计算结果。总的来说，采用压力模型 后，对于湍流燃烧释热过程的捕捉相比于 模型以及传统的压力模型 更加精确，混合层动量厚度增长与 结果吻合更好。在混合层发展初始阶段（无量纲时间 时间段），由于燃料和氧化剂尚未充分混合，同时有限的混合气受高温作用时间较短，整体反应十分微弱，燃烧释热基本可以忽略，此时不同 方法得到的动量厚度基本保持一致。值得注意的是，得到的动量厚度与 存在一定误差，原因在于初始化时给定初始能谱的波数不同，使得初始湍流分布存在一定差别，再加上网格分辨率大不相同，从而使得 和 的混合层初始发展过程显著不同。随着混合层的继续发展，燃料与高温氧化剂得到初步的混合，链式反应开始变得重要，燃烧释热的影响开始凸显。此后（无量纲时间 时刻后）模型和 模型得到的混合层厚度发展便不再一致。大约在无量纲时间 时，混合层开始进入自相似的充分发展阶段，此时整体化学反应更加剧烈，燃烧释热引起的条件亚格子项的作用开始凸显，从而导致不同滤波压力模型的 方法计算所得的混合层厚度增长开始出现明显不同。图 不同方法计算得到的动量厚度 图 给出了无量纲时间 时刻的涡量等值面，并分别采用以可解压力进行无量纲的亚格子压力以及无量纲燃烧释热率进行着色。可以看出，燃烧放热反应主要在氧化剂一侧进行，同时燃烧释热的增强会使得亚格子组分 温度关联项（）亚格子压力（）的影响增加。尽管相比于可解压力小 个量级，但由于每步迭代都会产生作用，通过时间累积效应对整体反应流场带来较大影响，如混合层动量厚度所示。为了进行定量对比，本文选取充分发展阶段无量纲时间 时刻的结果进行统计，不同方法所得的平均结果如图 所示。平均处理沿着周期性边界方向（和 方向），具体来说，就是针对第 期张林，等：可压缩流大涡模拟耦合概率密度函数方法中的滤波压力模型（）无量纲释热率（）图 无量纲涡量 等值面 每个无量纲 值，计算该 值平行于 平面的参数平均值。可以看出，在充分发展湍流阶段，湍流脉动十分强烈，湍流 化学反应相互作用不能忽（）温度分布（）（）氢氧基分布（）（）氢自由基分布（）（）超氧化氢基分布（）（）过氧化氢分布（）图 充分发展阶段（时刻）统计结果对比 （）国 防 科 技 大 学 学 报第 卷略，模型由于忽略了湍流 化学反应相互作用带来了较大误差，特别是在高温的氧化剂一侧。混合层靠氧化剂一侧温度较高，此时自点火效应变得显著。数值计算中一般采用超氧化氢基（）作为自点火过程的标识。可以看出，在氧化剂一侧，混合气在高温作用下，生成 的相关关链式反应以较快速度进行，使得 快速累积并达到峰值。即使是在燃料流一侧，模型得到的 也存在明显偏差。总体来说，方法由于忽略了亚格子尺度湍流 化学反应的相互作用，不能准确地描述包含自点火特征的湍流燃烧过程。如前文所述，有限样本的统计值相较于真实值的误差与样本数量的二分之一次方成反比，因此提高粒子数目，预期可以同时改善匹配压力模型 和压力模型 的 方法计算结果。但由于压力模型 采用的压力实际为 可解压力，本身不直接受粒子数目的影响，而是通过流场演化间接作用，因此粒子数目对于模型的改进主要体现为粒子化学反应计算本身及流场发展。而对于压力模型 ，粒子数目的增加还同时直接改善压力本身的计算结果，实现双向正反馈耦合。因此可以认为，增加粒子数目，压力模型 的改善可能会比压力模型 更加明显。由于粒子数目的影响并非关注的重点，同时针对本算例，本文采用的粒子数目可以获得较好的结果及不同模型的对比效果，也满足当前 方法中的通用粒子数目设置，额外考虑计算量，不再对不同粒子数目开展进一步计算分析。图 时刻统计结果 为了探索不同滤波压力模型的时间累积误差，在充分发展阶段的较早时刻（），也对时均流场进行了统计，对比如图 所示。可以看出，压力模型 在中心峰值之外出现了另一小的峰值，这主要是由于压力模型 中仍包含一定程度的统计噪声。同时，相比于压力模型 ，压力模型 的改善效果并不十分明显，这主要是由于流动进入充分发展阶段的时间还比较短，亚格子组分 温度关联项相对可解压力较小，需要经过长时间的积累才能展现优势。结论面向超声速反应流 ，发展了一种新的滤波压力模型。该模型综合了有限差分 和拉格朗日粒子 的优势，可以合理充分考虑亚格子组分 温度关联项，同时避免了直接使用 滤波压力带来的严重统计噪声问题，具备较好的数值稳定性。从 与流体保持概率等效的基本特征出发，和 中亚格子压力与可解压力的比值理论上具有一致性，基于此对 中的亚格子压力项进行建模。在三维超声速氢气 空气时间发展反应混合层中进行了数值测试，分别采用不同压力模型的 和 以及 对该算例进行模拟并对结果进行对比分析。研究显示，相比于传统的压力封闭模型，新发展的滤波压力模型可以明显改善 方法的计算精度。相比于传统的 方法，基于新滤波压力模型的 方法可以较好地捕捉时间发展反应混合层中亚格子尺度强烈的湍流 化学反应相互作用和亚格子尺度的分子 湍流混合、扩散。特别地，方法可以较好地模拟链式反应过程中的中间微小组分如过氧化氢的反应变化过程，可以更有效地模拟如自点火等复杂的高速湍流燃烧过程中的典型现象。致谢相关数值模拟在国家超级计算长沙中心开展，谨致谢意！参考文献（），：，：，：，第 期张林，等：可压缩流大涡模拟耦合概率密度函数方法中的滤波压力模型 ，（）：岳连捷，张旭，张启帆，等高马赫数超燃冲压发动机技术研究进展 力学学报，（）：，（）：（），（）：（），（）：，（）：，（）：，（）：，（）：，（）：，：，（）：，（），（）：王方，窦力，魏观溢，等基于 模型的凹腔支板火焰稳定器模拟 工程热物理学报，（）：，（）：（），（）：?，（）：，（）：，汪洪波超声速气流中凹腔稳定的射流燃烧模式及振荡机制研究 长沙：国防科技大学，：，（），：，（）：，（）：，（），（）：张林高速湍流燃烧 方法及其应用研究 长沙：国防科技大学，：，（）?，（）：，：