CN113221473B - 发动机燃烧室内气体-液滴两相流动特性的数值模拟方法 - Google Patents

Abstract

本发明公开了一种发动机燃烧室内气体‑液滴两相流动特性的数值模拟方法，所述模拟方法包括如下步骤：构建发动机燃烧室的几何模型；当液滴温度小于液滴的沸点温度时采用光滑离散粒子法对液滴碰撞过程进行模拟，采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟；当液滴温度不小于液滴的沸点温度时，采用光滑离散粒子法对液滴碰撞过程进行模拟，采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟，采用光滑离散粒子法和有限体积耦合法对液滴蒸发与燃烧过程进行模拟。该方法具有计算量小、稳定性高、液滴属性可调整、液滴轨迹可追踪、相间耦合机制灵活等优势。

Description

发动机燃烧室内气体-液滴两相流动特性的数值模拟方法

技术领域

本发明涉及燃油雾化特性研究技术领域，特别涉及发动机燃烧室内气体-液滴两相流动特性的数值模拟方法。

背景技术

航空航天作为世界经济军事强国的重要战略性产业，不仅对一个国家的国防现代化具有重要作用，而且对其经济发展同样具有重要意义。而气体-液滴两相流作为航空发动机及火箭发动机燃烧室中重要的一类现象，对其测量及应用技术的研究，逐渐成为航空航天领域中的关键技术和尖端技术。

在航空发动机或液体火箭发动机中，雾化燃烧特性对于发动机燃烧效率、污染物排放和燃烧稳定性具有重要的影响。航空燃油或液体推进剂经喷嘴喷出后，在自身不稳定波或外部气流场的作用下，破碎成小液滴，小液滴间频繁发生着碰撞聚合、碰撞反弹、碰撞破碎以及液滴的剪切破碎等现象，同时液滴在高温高压环境下蒸发、掺混进而燃烧。经分析发现，液滴与气相之间形成的两相流区别于传统的气-固两相流(气体-固体液滴两相流)和气-液两相流(气体-液体连续界面流)，液滴不仅与气相之间具有连续性界面，液滴在界面张力的作用下变形和破碎，同时液滴又是以分散质点的形式存在，数量众多；液滴间不仅可以像固体液滴碰撞一样发生反弹现象，同时液滴间还可以发生碰撞聚合以及碰撞破碎等现象，液滴粒径在时刻发生改变；液滴与气体相之间不仅具有动量和能量交互，同时液滴也可发生蒸发与燃烧现象，与气体相之间具有质量传递。

对于气体-液滴两相流现有的数值方法有：基于流体体积函数的界面追踪方法、液滴轨道追踪方法、双流体模型方法、流体拟液滴方法。

基于流体体积函数的界面追踪方法的基本思想是，采用体积率函数C(x)描述网格单元中某一流体所占的体积率，当网格单元完全由某一种流体占据，则C＝1，当网格单元中完全不含该流体，则C＝0，对于两种流体的边界有0<C<1。流体体积函数是一种网格固定的表面跟踪技术。该模型用于观察两种及以上互不相融流体间的分界面。流体体积函数模型中，两种流体共用一组动量方程，计算域中各流体的体积分数在每个计算单元上被跟踪。方法中，首先通过求解一组动量方程来模拟两种及两种以上互不相融的流体，在整个计算域跟踪各流体的体积分数；然后对给定的界面流体体积函数进行界面重构；最后再计算下一时刻各网格单元中流体体积函数的值。

液滴轨道追踪方法是基于欧拉-拉格朗日框架所构建，该方法认为流体相是连续的，液滴相是离散的。根据气体-液滴两相系统中液滴运动的特点，该方法将两相流动中液滴的运动过程分解为受冲力支配的瞬时碰撞运动及受流体曳力控制的悬浮运动，从而建立了液滴运动分解模型。该方法中，流体的运动规律用连续介质的N-S方程进行描述，而液滴的行为则通过在拉格朗日坐标系中分析每一个单液滴的运动轨迹而进行描述。其中，在液滴与液滴相互作用的过程中，其运动规律服从碰撞动力学中的动量守恒定律；在流体与液滴相互作用的悬浮过程中，液滴在曳力、重力等力的作用下，其运动规律服从牛顿动力学中的力平衡方程。这样，每个液滴的速度及位移的更新由邻近液滴对它的碰撞作用及流体对它的悬浮作用来确定。与此同时，液滴对流体的瞬时作用则反映在离散相与流体相两相耦合的N-S方程不断进行修正的数值求解过程中。

双流体模型方法是基于欧拉框架进行构建的，是一种相对较为复杂的多相流模型，它将各相看成相互渗透、融合但又具有不同运动特征的连续介质。该模型通过建立一套包含有两组动量方程和连续性方程的方程组来求解每一相，压力项和各界面交换系数是两相共有，通过各相所占据的体积分数比重加权平均得到。不同相之间的动量交换也依赖于混合物的类别。该模型已广泛应用于气泡柱、气泡上浮、液滴悬浮以及流化床等领域。相对于单流体模型，双流体模型考虑了液滴相的湍流输运以及气体-液滴两相间相互滑移引起的阻力，其计算结果在大多数情况下更贴近于实际。该模型可以对各相进行单独的计算，每一相都有独自的守恒方程，具有很大的适应性。但随着相数的增加，方程数相应增加，使得计算逐渐趋于复杂，特别是当液滴相结构尺寸存在较大差别时，各相都要进行独自计算最后进行迭代直至收敛，计算量巨大。此外液滴相的扩散系数、导热系数和粘度系数的确定也仍然存在问题，并且双流体模型在计算中还可能产生伪扩散。

流体拟液滴方法是一种底层的粒子方法，它将流体离散为大量自由运动并同步碰撞的硬球，从而结合了另外两种粒子方法的优势，即分子动力学模拟的可靠性和表现力，加之直接模拟蒙特卡罗的效率和简单性，因此它已能在微机上实现微观液滴流体系统的模拟。后来，为了克服拟液滴模拟对实际系统过于基础从而低效的难题，又以类似SPH中所采用的流体微元间作用代替粒子间的碰撞，称之为宏观拟液滴模拟。

界面追踪方法较难用于大量液滴的界面追踪和蒸发燃烧过程计算：采用界面追踪方法对液滴的碰撞及破碎过程进行直接数值模拟，虽然可以捕捉物理过程的细节，分析物理机理，但不适合于对工程实际问题的研究，因为通常发动机燃烧室中包含有数万甚至数亿个液滴，液滴大小不一，并且有时差别较大，现有计算机硬件远远不能达到要求。同时，界面追踪方法与蒸发燃烧模型的耦合算法复杂，需要在体积分数和质量分数、相界面追踪和物质组分追踪之间做好区别和处理，目前鲜有文献对大量液滴的蒸发燃烧问题采用该方法。

液滴轨道模型依赖于碰撞概率、碰撞结果预测、二次破碎结果预测等模型且计算量大：采用液滴轨道模型计算液滴间的碰撞需要借助直接蒙特卡罗概率判断碰撞与否，再结合液滴碰撞模型直接给出碰撞结果，计算碰撞的概率精度低，同时必须依赖于现有已发展的液滴间碰撞模型和液滴的二次破碎模型，通用性差；该模型计算单液滴的受力基于牛顿运动第二定律，通过计算每一个液滴所受到的所有外部作用力(包括液滴间作用力)跟踪液滴的轨迹，当液滴数量较多时，计算量将迅速增大，虽然可以采用液滴群的概念降低计算量，但对于液滴尺寸分布较广的情形，同样计算量无法有效降低，该模型不适合大规模计算；另外，该方法本质属于微观方法，计算参数需要进行统计才能与宏观可测参数相匹配，不利于实验验证。

双流体模型等对于气体-液滴两相流计算存在一些不足：双流体模型相比较液滴轨道模型，在计算量方面大幅降低，同时可直接获取宏观可测参量，全面考虑液滴的湍流输运作用，但该模型仅能获得网格上面的液滴体积分数信息，无法追踪液滴的运动轨迹，捕捉液滴运动的细节，易产生伪扩散，同时不易加入液滴蒸发、燃烧等物理化学模型。流体拟液滴模型是将气体及液滴的运动均在拉格朗日体系中描述，该模型仅适合于微重力场、低雷诺数下流动过程的模拟，模型的准确程度与“拟液滴”的大小密切相关，“拟液滴”尺寸越小，液滴数越多，计算结果越准确，但工作量也就越大。

发明内容

本发明的目的是一种发动机燃烧室内气体-液滴两相流动特性的数值模拟方法，以在提高模拟结果的准确性的同时减少工作量。

为实现上述目的，本发明提供了如下方案：

一种发动机燃烧室内气体-液滴两相流动特性的数值模拟方法，所述数值模拟方法包括如下步骤：

构建发动机燃烧室的几何模型；

判断液滴温度是否小于液滴的沸点温度，获得第一判断结果；

若所述第一判断结果表示是，则基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟；

基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟；

若所述第一判断结果表示否，则基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟；

基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟；

基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴蒸发与燃烧过程进行模拟；

返回步骤“判断液滴温度是否小于液滴的沸点温度，获得第一判断结果”，继续下一个时间步的模拟，直到模拟时间达到预设模拟时间。

可选的，所述基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟，具体包括：

给定初始的液滴分布矩；

根据上一时刻的液滴分布矩，利用公式，

计算每个积分点的当前时刻的权值；

其中，

表示在积分点m处的液滴的粒径的k阶矩，ω_k(t)表示当前时刻t的k阶矩的权值，m_k(t-1)表示前一时刻的液滴尺度分布函数的k阶矩，m＝0,1,2,…,N-1,k＝0,1,2,…,N-1，N表示积分点的数量；

根据当前时刻的权值，利用公式

计算矩量输运方程的聚合与破碎源项

其中，ω_α(t)和ω_γ(t)分别表示当前时刻t的积分点α和积分点γ的权值，L_α和L_γ分别表示在积分点α和积分点γ处的液滴的粒径，β_αγ表示粒径为L_α,L_γ的两液滴碰撞的概率，β_αγ＝β(L_α,L_γ)，β(·)为通过实验总结获得的碰撞概率函数；

表示在积分点α处的液滴破碎的概率，

表示在积分点α处的液滴与其它液滴碰撞的概率；

表示在积分点b处的液滴破碎的概率；

表示在积分点m处的液滴的粒径的k阶矩；

根据矩量输运方程的聚合与破碎源项求解矩量输运方程

获得当前时刻的液滴分布矩；

其中，

为k阶矩的速度；

基于当前时刻的液滴分布矩，根据液滴分布矩与液滴数量、液滴分布的标准差和粒径平均值的对应关系，求解液滴数量、标准差和粒径平均值。

可选的，所述采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟，具体包括：

将以粒子形式存储的网格中心点的信息采用核函数插值的方式插值到SDPH粒子点位置，获得SDPH粒子点位置的网格中心点信息；所述网格中心点信息包括液滴的位置、速度、体积分数、气场压力、湍流耗散量；

根据SDPH粒子点位置的两相体积分数、气相速度和湍流耗散值，求解群体平衡方程

获得不同时刻液滴的空间状态特征；

其中，n(L,t)为t时刻尺寸为L的液滴的数量，

为阶矩的平均速度，B(L,t)和D(L,t)分别表示t时刻尺寸为L的液滴由于聚合和破碎造成的生函数和死函数；B^B(L)表示尺寸为L的液滴由于破碎造成的生函数项

b(L|λ)表示尺寸为L和尺寸为λ的两液滴由于碰撞后破碎形成的子液滴的尺寸分布函数，a(L)表示尺寸为λ的液滴的破碎频率，n(λ)表示尺寸为λ的液滴的数量；B^C(L)表示尺寸为L的液滴由于聚合作用产生的生函数项

β((L³-λ³)^1/3,λ)表示尺寸为(L³-λ³)^1/3和尺寸为λ的两液滴聚合的频率，n((L³-λ³)^1/3)表示尺寸为(L³-λ³)^1/3的液滴的数量；具体的，通过实验数据总结拟合得到尺寸分布函数与We和χ数之间的关系式。液滴碰撞参数韦伯数和偏心度的计算方式为：

韦伯数

偏心度

式中，σ为表面张力系数，D为液滴直径，v_rel为两液滴之间的相对速度，X为两液滴的偏心距离。

根据SDPH粒子点位置的网格中心点信息，计算SDPH粒子点位置所受到的气相曳力、气相压力和热传导量；由网格中心点的信息计算气相曳力就是采用曳力公式R_gp＝β_gp(v_g-v_p)计算。R_gp为气体-颗粒间的曳力，β_gp为气体-颗粒之间的曳力系数，v_g为网格中心点的气相速度矢量，v_p为SDPH粒子点位置的速度矢量。

系数C_D为：

相对雷诺数Re_p定义为：

α_g为网格中心点气相的体积分数；α_p为SDPH粒子的体积分数；d_p为颗粒的直径；μ_g为网格中心点气相的粘度。

将SDPH粒子点位置所受到的气相曳力、气相压力和热传导量作为矩量输运方程的聚合与破碎源项，并根据不同时刻液滴的空间状态特征，求解矩量输运方程获得SDPH粒子的粒径的单位时间变化量，求解能量方程获得SDPH粒子的质量、动量和能量的单位时间变化量；本发明根据聚合与破碎源项的实验拟合公式计算在具体的气相曳力、气相压力和热传导量作用下的源项。

根据SDPH粒子的粒径、质量、动量和能量的单位时间变化量，采用蛙跳法更新SDPH粒子点位置的网格中心点信息；

判断当前仿真时长是否小于FVM时间步长，获得第二判断结果；

若所述第二判断结果表示是，则返回步骤“根据SDPH粒子点位置的两相体积分数、气相速度和湍流耗散值，求解群体平衡方程

获得不同时刻液滴的空间状态特征”；

若所述第二判断结果表示否，则将SDPH粒子点位置的网格中心点信息采用核函数反插值的方式插值到网格中心位置处；

根据网格中心位置处的信息和网格数据，计算气相受到的液滴相曳力和热传导量，作为FVM输运方程的源项；FVM输运方程为：

其中I为单位矩阵，V为流体所占据体积，S为占据体积边界的面积，n为垂直于面S的单位法向量；

采用求解FVM输运方程的方式，求解FVM液滴的质量、动量和能量的单位时间变化量；

FVM输运方程求解形式为

根据FVM液滴的质量、动量和能量的单位时间变化量，采用蛙跳法更新网格中心位置处的信息；

判断FVM是否收敛，获得第三判断结果；

若所述第三判断结果表示否，则返回步骤“根据网格中心位置处的信息和网格数据，计算气相受到的液滴相曳力和热传导量，作为FVM输运方程的源项”；

判断FVM计算总时间是否超过总时间阈值，获得第四判断结果；

若所述第四判断结果表示否，则将SDPH计算时间设置为0，返回步骤“将以粒子形式存储的网格中心点的信息采用核函数插值的方式插值到SDPH粒子点位置，获得SDPH粒子点位置的网格中心点信息”；

若所述第三判断结果表示是，停止仿真。

可选的，所述基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴蒸发与燃烧过程进行模拟，具体包括：

将以粒子形式存储的网格中心点的信息采用核函数插值的方式插值到液滴粒子点位置，获得液滴粒子点位置的网格中心点信息；网格中心点信息包括液滴的位置、速度、体积分数、气场压力、湍流耗散量；

根据液滴粒子点位置的网格中心点信息，利用公式

计算单液滴蒸发量；其中，m_p为液滴质量，N_i为组分i的液滴的数量，A_p为液滴的表面积，M_w,i为组分i的摩尔质量，d_p为液滴直径，α_p为液滴的体积分数；

根据单液滴蒸发量更新液滴的直径、质量、体积分数和光滑离散粒子密度，利用大密度差多相流粒子方法和完全变光滑长度法，求解液滴粒子的密度、速度、光滑长度和位移的单位时间变化量；

根据液滴粒子的密度、速度、光滑长度和位移的单位时间变化量采用蛙跳方式，更新液滴粒子的质量、速度和能量；

将液滴粒子点位置的网格中心点信息和单液滴蒸发量采用核函数反插值的方式插值到网格中心位置处；

根据单液滴蒸发量计算每个网格的蒸发量和体积分数；

利用气相燃烧模型

计算化学反应的平均速率

其中，ω_t为基于k-ε湍流模型计算的湍流化学反应速率，ω_t＝min(ω_t1,ω_t2)；ω_t1表示已燃气体微团破碎速率，ω_t2表示未燃气体微团破碎速率；ω_t1表示已燃气体微团破碎速率，ω_t2表示未燃气体微团破碎速率，

v_i为反应物i的化学计量数；v_R为生成物R的化学计量数；M_i为反应物i的摩尔质量；M_R为生成物R的摩尔质量；A、B为经验常数，A＝4.0，B＝0.5，ε为湍流耗散率，k为湍流强度，ρ为液滴的密度，m_R为生成物的质量分数，m_i为反应物i的质量分数，N为总的反应物的数目。

根据每个网格的蒸发量、体积分数和化学反应的平均速率，采用求解组分输运方程和N-S方程的方式求解网格的质量、速度和能量的单位时间变化量；具体的，利用组分运输方程获得每个网格上的组分的分布数据，将这些数据作为N-S方程中连续性方程的源项，计算N-S方程，获得每个网格质量、速度和能量的时间变化量，进而采用蛙跳时间更新获得新的时刻每个网格的质量、速度和能量。

根据质量、速度和能量的单位时间变化量，采用蛙跳时间更新获得新的时刻每个网格的质量、速度和能量。

返回步骤“将以粒子形式存储的网格中心点的信息采用核函数插值的方式插值到液滴粒子点位置，获得液滴粒子点位置的网格中心点信息”，当总的计算时间大于设定的终止时间时，计算终止。

可选的，所述利用大密度差多相流粒子方法和完全变光滑长度法，求解液滴粒子的密度、速度、光滑长度和位移的单位时间变化量，具体包括：

大密度差多相流粒子方法的公式为：

完全变光滑长度方法的公式为：

其中，ρ_i表示组分i和组分j的密度，m_i和m_j分别表示组分i和组分j的密度，V_i和V_j分别表示组分i组分j的体积，v_i表示组分i的速度，v_ij表示组分i与组分j的相对速度，x_i表示组分i的位置，W_ij表示粒子i和粒子j之间的核函数数值，P_i和P_j分别表示粒子i和粒子j的压力值，μ_i和μ_j分别表示粒子i和粒子j的动力粘度值，r_ij表示粒子i和粒子j之间的距离，

表示粒子i和粒子j之间的平均压力，

表示粒子i和粒子j之间的人工应力值，Π_ij表示粒子i和粒子j之间的人工粘性值，h_i和h_j分别表示组分i和组分j的能量，f_ij表示粒子i与粒子j之间的人工应力系数，用公式表示为f_ij＝W_ij/W(Δp)，Δp为粒子之间的初始间距，W(Δp)为粒子之间的距离为初始间距时的粒子之间的核函数数值。

可选的，所述组分输运方程为：

其中，

为组分i的化学反应的平均速率，S_i为离散相及自定义的源项产生的额外物质质量生成速率，J_i为组分i的扩散通量；Y_i表示组分i的主流蒸汽浓度，ρ和v分别表示液滴的密度和速度；

ρ为液滴密度；D_i,m为蒸汽扩散系数；Sc为施密特数；μ为液滴的动力粘度系数；Y_i表示组分i的主流蒸汽浓度。

本发明还提供一种发动机燃烧室内气体-液滴两相流动特性的数值模拟系统，所述数值模拟系统包括：

几何模型判断模块，用于构建发动机燃烧室的几何模型；

第一判断模块，用于判断液滴温度是否小于液滴的沸点温度，获得第一判断结果；

第一液滴碰撞过程模拟模块，用于若所述第一判断结果表示是，则基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟；

第一液滴剪切破碎过程模拟模块，用于基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟；

第二液滴碰撞过程模拟模块，用于若所述第一判断结果表示否，则基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟；

第二液滴剪切破碎过程模拟模块，用于基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟；

液滴蒸发与燃烧过程模拟模块，用于基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴蒸发与燃烧过程进行模拟；

返回模块，用于返回步骤“判断液滴温度是否小于液滴的沸点温度，获得第一判断结果”，继续下一个时间步的模拟，直到模拟时间达到预设模拟时间。

可选的，所述第一液滴碰撞过程模拟模块，具体包括：

液滴分布矩初始化子模块，用于给定初始的液滴分布矩；

权值计算子模块，用于根据上一时刻的液滴分布矩，利用公式

计算每个积分点的当前时刻的权值；

其中，

表示在积分点m处的液滴的粒径的k阶矩，ω_k(t)表示当前时刻t的k阶矩的权值，m_k(t-1)表示前一时刻的液滴尺度分布函数的k阶矩，m＝0,1,2,…,N-1,k＝0,1,2,…,N-1，N表示积分点的数量；

矩量输运方程的聚合与破碎源项计算子模块，用于根据当前时刻的权值，利用公式

计算矩量输运方程的聚合与破碎源项

其中，ω_α(t)和ω_γ(t)分别表示当前时刻t的积分点α和积分点γ的权值，L_α和L_γ分别表示在积分点α和积分点γ处的液滴的粒径，β_αγ表示粒径为L_α,L_γ的两液滴碰撞的概率；

表示在积分点m处的液滴的粒径的k阶矩；

表示在积分点α处的液滴破碎的概率，

表示在积分点α处的液滴与其它液滴碰撞的概率；

表示在积分点b处的液滴破碎的概率；

液滴分布矩求解子模块，用于根据矩量输运方程的聚合与破碎源项求解矩量输运方程

获得当前时刻的液滴分布矩；

其中，

为k阶矩的速度；

分布数据求解子模块，用于基于当前时刻的液滴分布矩，根据液滴分布矩与液滴数量、液滴分布的标准差和粒径平均值的对应关系，求解液滴数量、标准差和粒径平均值。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明公开了一种发动机燃烧室内气体-液滴两相流动特性的数值模拟方法，所述模拟方法包括如下步骤：构建发动机燃烧室的几何模型；当液滴温度小于液滴的沸点温度时采用光滑离散粒子法对液滴碰撞过程进行模拟；采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟；当液滴温度不小于液滴的沸点温度时，采用光滑离散粒子法对液滴碰撞过程进行模拟；采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟，采用光滑离散粒子法和有限体积耦合法对液滴蒸发与燃烧过程进行模拟。本发明采用无网格粒子法与网格法耦合的数值模拟方法，进行发动机燃烧室内气体-液滴两相流动特性的数值模拟，该方法具有计算量小、稳定性高、液滴属性可调整、液滴轨迹可追踪、相间耦合机制灵活等优势，同时具有较好的实用性和可拓展性，对于解决气体-液滴两相流问题具有一定的优势。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明提供的一种发动机燃烧室内气体-液滴两相流动特性的数值模拟方法的流程图；

图2为本发明提供的三维几何模型的结构图；图2(a)为喷嘴剖面图，

图2(b)为喷嘴外部整体图，图2(c)为燃烧室火焰筒结构图；

图3为本发明提供的液滴碰撞过程实验结果图；

图4为本发明提供的采用光滑粒子方法对液滴碰撞过程的数值模拟结果图；

图5为本发明提供的液滴碰撞结果与液滴碰撞参数之间的关系图；

图6为本发明提供的采用光滑离散粒子法对液滴碰撞过程进行模拟的流程图。

具体实施方式

本发明的目的是一种发动机燃烧室内气体-液滴两相流动特性的数值模拟方法，以在提高模拟结果的准确性的同时减少工作量。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对发明作进一步详细的说明。

本发明方案的核心技术包括：建立了描述液滴间碰撞聚合、碰撞破碎和单液滴剪切破碎行为的宏观方程式，解决了液滴粒径改变光滑粒子方法模拟的关键科学问题，突破光滑粒子离散求解部分微积分方程关键技术，实现光滑粒子无网格法和有限体积法对相间传质耦合作用的描述，形成一套可完整求解液滴与气相复杂化学反应及含液滴碰撞和破碎过程的新型光滑粒子法-有限体积耦合方法，并在典型发动机燃烧室内两相雾化燃烧问题上实现应用示范。

如图1所示，本发明提供一种发动机燃烧室内气体-液滴两相流动特性的数值模拟方法，所述数值模拟方法包括如下步骤：

步骤101，构建发动机燃烧室的几何模型。

根据实际发动机燃烧室几何构型，建立燃烧室壳体、燃油雾化喷嘴、旋流器、机匣等结构的三维几何模型如图2所示。

步骤102，判断液滴温度是否小于液滴的沸点温度，获得第一判断结果；

步骤103，若所述第一判断结果表示是，则基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟。

步骤103所述基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟，具体包括：给定初始的液滴分布矩；根据上一时刻的液滴分布矩，利用公式，

计算每个积分点的当前时刻的权值；其中，

表示在积分点m处的液滴的粒径的k阶矩，ω_k(t)表示当前时刻t的k阶矩的权值，m_k(t-1)表示前一时刻的液滴尺度分布函数的k阶矩，m＝0,1,2,…,N-1,k＝0,1,2,…,N-1，N表示积分点的数量；根据当前时刻的权值，利用公式

计算矩量输运方程的聚合与破碎源项

其中，ω_α(t)和ω_γ(t)分别表示当前时刻t的积分点α和积分点γ的权值，L_α和L_γ分别表示在积分点α和积分点γ处的液滴的粒径，β_αγ表示粒径为L_α,L_γ的两液滴碰撞的概率，β_αγ＝β(L_α,L_γ)，β(·)为通过实验总结获得的碰撞概率函数；

表示在积分点α处的液滴破碎的概率，

表示在积分点α处的液滴与其它液滴碰撞的概率；

表示在积分点b处的液滴破碎的概率；

表示在积分点m处的液滴的粒径的k阶矩；根据矩量输运方程的聚合与破碎源项求解矩量输运方程

获得当前时刻的液滴分布矩；其中，

为k阶矩的速度；基于当前时刻的液滴分布矩，根据液滴分布矩与液滴数量、液滴分布的标准差和粒径平均值的对应关系，求解液滴数量、标准差和粒径平均值。

具体的，该方案中包含三个大的步骤，首先建立描述液滴在碰撞过程中发生粒径变化和体积变化的输运方程式(群体平衡模型)，然后对方程式中的液滴碰撞聚合和碰撞破碎进行理论研究，建立宏观描述的物理模型，最后采用光滑离散粒子数值方法结合直接矩积分方法进行求解。具体每个大的步骤中包含的细节描述如下：

(1)液滴碰撞过程中粒径变化和体积变化的有效描述

通过描述液滴在碰撞过程中发生粒径变化和体积变化的方程式来实现这一步骤，首先引入液滴数目密度的连续形式(公式(1))，表述为在点x＝(x₁；x₂；x₃)附近体积dx＝dx₁dx₂dx₃液滴尺寸(L,L+dL)范围内液滴单元的数量

n(L；x,t)dxdL (1)

在公式(1)的基础上，引入液滴碰撞的概率模型函数(下面步骤(2)详细描述)，包括：描述液滴碰撞聚合的概率模型函数β(L,λ)，表征尺寸为L和λ的两液滴聚合的频率；描述液滴碰撞后破碎的概率模型函数b(L|λ)，表征尺寸为和λ的两液滴由于碰撞后破碎形成的子液滴的尺寸分布函数；描述单液滴气动破碎的概率模型函数a(L)L，表征尺寸为L的液滴的破碎频率，建立尺寸为L的液滴由于聚合和破碎造成的生函数和死函数B(L；x,t)和D(L；x,t)。其中，由于破碎造成的生函数项为

由于聚合作用产生的生函数项为

由于破碎和聚合作用造成的液滴死函数为

D^B(L)＝a(L)n(L) (4)

在上述获得的生死函数的基础上，根据公式(1)的液滴数目密度表达式，结合物质输运表达式，建立描述液滴间碰撞微观行为及粒径变化的宏观方程式(6)，也就是群体平衡方程

(2)液滴碰撞聚合和碰撞破碎的宏观模型描述

上述步骤(1)公式(6)中右端项中的生函数和死函数，如公式(2)、(3)、(4)、(5)，均表示成了液滴碰撞聚合、碰撞破碎概率密度函数的形式，那么具体的概率密度函数是什么，生函数和死函数与实际碰撞过程中的韦伯数、雷诺数、偏心度等参量之间的关系是什么，需要本步骤进行推导获得，具体细节描述如下：

图3为国外学者针对液滴碰撞问题开展的实验获得的结果，图4为发明人采用光滑粒子方法针对液滴碰撞问题开展的数值模拟工作。在这些结果的基础上，总结获得不同条件下的液滴碰撞结果分布图，即液滴碰撞结果与液滴碰撞参数之间的关系图，如图5所示，有了该结果便可在已知液滴碰撞参数的基础上(碰撞的韦伯数We和偏心度χ)，明晰液滴碰撞后的结果，也就是获得了液滴由于聚合和破碎造成的生函数B(L；x,t)和死函数D(L；x,t)与液滴碰撞参数之间的具体关系式B(L；x,t)＝f_B(We,χ)，D(L；x,t)＝f_D(We,χ)。

(3)液滴碰撞物理模型的数值离散求解

通过步骤(1)和(2)建立了描述液滴间碰撞问题的物理模型(公式6)，本步骤是采用一定的方法对该物理模型进行数值离散，从而能够进行数值求解，获得不同时刻液滴的空间状态特征。

首先采用矩方法对液滴数目密度表达式(1)进行矩量积分，如式(7)，得到液滴尺度分布函数的k阶矩

当k取不同值时，液滴分布矩m_k的含义不同，并且低阶矩与表征液滴分布的特定物理性质相关。当k＝0时，m₀表征系统中液滴的总数量；当k＝2时，m₂与面积形状因子k_A的乘积表示液滴的总表面积，即A_f＝k_Am₂(t)；当k＝3时，m₃与体积形状因子k_V的乘积表示液滴的总体积，即V_t＝k_Vm₃(t)。因此，根据液滴的平均索太尔直径d₃₂的定义：液滴总体积与总表面积的比值，可获得液滴群的平均索太尔直径为

根据光滑液滴离散动力学(SDPH)方法的定义可知，SDPH单粒子代表具有一定分布形态的液滴群，液滴群的粒径均值(平均索太尔直径)、方差及液滴数量由SDPH粒子计算表征，基于该定义本发明就建立SDPH粒子所表征的液滴的粒径参量与粒数密度的低阶矩量(公式(7))之间的对应关系，比如假设液滴群的初始分布为对数正态分布，则液滴数密度可写为

f₁为液滴数密度与均值粒径、方差及液滴数量的函数关系。

接下来选择不同的数值密度函数，得到不同的矩方法，为降低计算量，本发明对数值密度函数采用如下连续Dirac Delta函数加权的形式，同时结合方程(9)可得

其中，N为积分点的个数，δ(L)为Dirac Delta函数，ω_i为Dirac Delta函数权值，随时间、均值粒径、方差及液滴数量变化。L_i为液滴的特征尺寸。

将方程(10)代入(7)得

利用方程(11)计算得到权值ω_i：取k＝0,…,N-1，方程(11)写成矩阵形式为

方程(12)为Vandermonde方程，首先取[0,∞]上N个拉塞尔高斯点，将值赋予L_i，i＝0,…,N-1，然后采用针对Vandermonde线性系统的求解算法进行求解，实现任意数目矩的计算，得到权值ω_i。

接下来将矩量方程(11)代入式(6)得矩量输运方程

其中，

为k阶矩的速度，它对于液滴尺度分布的各阶矩理论上不同。本发明中仅计算一种液滴速度分布，

对于各阶矩均相同。将方程(13)右端聚合与破碎源项写为矩量形式

进一步写为直接矩积分方式

其中，

当涉及子液滴的分布函数时，使用以下公式

m和n表征两个破碎液滴间的质量比的关系。如果m＝1，n＝1，表示两个破碎液滴具有相同的体积，因此考虑为对称破碎形式；如果m≠n，破碎非对称，特定情况如m远大于n(或n远大于m)则为通常所熟知的侵蚀破碎。在本发明中，仅考虑对称破碎过程，即m＝n＝1。

由上述阐述的由公式(12)计算得到的权值ω_i，特征长度L_i及聚合破碎源项(19)求解矩量输运方程(10)，得到下一时刻各个矩量值，再转到针对Vandermonde线性系统的求解算法求解方程(12)，循环计算，直至模拟结束。

最后再由SDPH假定液滴群服从的对数正态分布的性质，根据粒径均值

粒径方差σ和总液滴数N为k阶矩m_k(t)的函数关系

进而由各阶矩量输运方程求解出

σ及N的数值，这样便实现了采用SDPH数值方法计算出液滴在不同时刻具体空间分布的目标，采用后处理对数据结果进行分析，从而预测发动机燃烧室内气体-液滴两相流在不同时刻的流动特性，为发动机部件的设计提供数据支撑。

步骤104，基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟。

步骤104所述采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟，具体包括：

将以粒子形式存储的网格中心点的信息采用核函数插值的方式插值到SDPH粒子点位置，获得SDPH粒子点位置的网格中心点信息；所述网格中心点信息包括液滴的位置、速度、体积分数、气场压力、湍流耗散量；

根据SDPH粒子点位置的两相体积分数、气相速度和湍流耗散值，求解群体平衡方程

获得不同时刻液滴的空间状态特征；

其中，n(L,t)为t时刻尺寸为L的液滴的数量，

为阶矩的平均速度，B(L,t)和D(L,t)分别表示t时刻尺寸为L的液滴由于聚合和破碎造成的生函数和死函数；B^B(L)表示尺寸为L的液滴由于破碎造成的生函数项

b(L|λ)表示尺寸为L和尺寸为λ的两液滴由于碰撞后破碎形成的子液滴的尺寸分布函数，a(L)表示尺寸为λ的液滴的破碎频率，n(λ)表示尺寸为λ的液滴的数量；B^C(L)表示尺寸为L的液滴由于聚合作用产生的生函数项

β((L³-λ³)^1/3,λ)表示尺寸为(L³-λ³)^1/3和尺寸为λ的两液滴聚合的频率，n((L³-λ³)^1/3)表示尺寸为(L³-λ³)^1/3的液滴的数量；具体的，通过实验数据总结拟合得到尺寸分布函数与We和χ数之间的关系式。液滴碰撞参数韦伯数和偏心度的计算方式为：

韦伯数

偏心度

式中，η为粘度系数，σ为表面张力系数，D为液滴直径，v_rel为两液滴之间的相对速度，X为两液滴的偏心距离。

根据SDPH粒子点位置的网格中心点信息，计算SDPH粒子点位置所受到的气相曳力、气相压力和热传导量；由网格中心点的信息计算气相曳力就是采用曳力公式R_gp＝β_gp(v_g-v_p)计算。R_gp为气体-颗粒间的曳力，β_gp为气体-颗粒之间的曳力系数，v_g为网格中心点的气相速度矢量，v_p为SDPH粒子点位置的速度矢量。

系数C_D为：

相对雷诺数Re_p定义为：

α_g为网格中心点气相的体积分数；α_p为SDPH粒子的体积分数；d_p为颗粒的直径；μ_g为网格中心点气相的粘度。

将SDPH粒子点位置所受到的气相曳力、气相压力和热传导量作为矩量输运方程的聚合与破碎源项，并根据不同时刻液滴的空间状态特征，求解矩量输运方程获得SDPH粒子的粒径的单位时间变化量，求解能量方程获得SDPH粒子的质量、动量和能量的单位时间变化量；本发明根据聚合与破碎源项的实验拟合公式计算在具体的气相曳力、气相压力和热传导量作用下的源项。

根据SDPH粒子的粒径、质量、动量和能量的单位时间变化量，采用蛙跳法更新SDPH粒子点位置的网格中心点信息；

判断当前仿真时长是否小于FVM时间步长，获得第二判断结果；

若所述第二判断结果表示是，则返回步骤“根据SDPH粒子点位置的两相体积分数、气相速度和湍流耗散值，求解群体平衡方程

获得不同时刻液滴的空间状态特征”；

若所述第二判断结果表示否，则将SDPH粒子点位置的网格中心点信息采用核函数反插值的方式插值到网格中心位置处；

根据网格中心位置处的信息和网格数据，计算气相受到的液滴相曳力和热传导量，作为FVM输运方程的源项；FVM输运方程为：

其中I为单位矩阵，V为流体所占据体积，S为占据体积边界的面积，n为垂直于面S的单位法向量；

采用求解FVM输运方程的方式，求解FVM液滴的质量、动量和能量的单位时间变化量；

FVM输运方程求解形式为

根据FVM液滴的质量、动量和能量的单位时间变化量，采用蛙跳法更新网格中心位置处的信息；

判断FVM是否收敛，获得第三判断结果；

若所述第三判断结果表示否，则返回步骤“根据网格中心位置处的信息和网格数据，计算气相受到的液滴相曳力和热传导量，作为FVM输运方程的源项”；

判断FVM计算总时间是否超过总时间阈值，获得第四判断结果；

若所述第四判断结果表示否，则将SDPH计算时间设置为0，返回步骤“将以粒子形式存储的网格中心点的信息采用核函数插值的方式插值到SDPH粒子点位置，获得SDPH粒子点位置的网格中心点信息”；

若所述第三判断结果表示是，停止仿真。

具体的，湍流气流场中液滴剪切破碎过程模拟与液滴的碰撞过程模拟存在着一定的共性问题，如均需要建立液滴生成与破碎的生函数和死函数(如上述1中的步骤(2))，再建立考虑液滴粒径改变的光滑离散粒子方法(如上述1中的步骤(3))，同时研究思路与描述液滴碰撞的方法相同，在现有实验数据和液滴二次破碎问题上数值模拟结果的基础上开展模型研究，这里不再重复介绍。

除此之外，在建立考虑液滴粒径改变的光滑离散粒子方法的基础上，需要考虑液滴和气体之间的耦合效应，这是与上述1存在差别的地方，上述1中不需要液滴和气体之间的数据交换。因此，这里将离散液滴所采用的光滑粒子方法和离散气流场所采用的有限体积方法之间的耦合过程详述如下：基于光滑离散粒子方法和有限体积方法之间的基础数据交换，根据剪切破碎后的粒径改变，计算液滴与气相场之间的曳力，更新两相各自的体积分数，实现光滑离散粒子法与有限体积方法之间的耦合。具体耦合过程如下：

(1)对问题建模和初始化，首先进行网格离散和初始化设置，而后SDPH模块进行粒子离散和初始化设置；

(2)根据系统的初始分布，利用方程(7)求解出初始分布矩，并根据方程(12)确定初始分布的权值ω_i，利用方程(22)计算液滴群初始平均粒径沿空间分布。

(3)将网格中心点的信息以粒子的形式进行存储，包括位置、速度、体积分数、气场压力、湍流耗散等，然后开始SDPH程序计算，进行临近粒子搜索和核函数计算，FVM等待；

(4)将以粒子形式存储的网格中心处的信息以核函数插值的方式，插值到SDPH粒子点位置处，计算SDPH粒子所受到气相曳力、压力、热传导量等，作为源项添加到SDPH输运方程中，进而计算SDPH质量、动量和能量守恒方程；

(5)利用步骤(4)插值到SDPH粒子处的两相体积分数、气相速度以及湍流耗散值，求解群体平衡方程；

(6)由CEL条件计算SDPH下一步的时间步长，采用蛙跳更新SDPH粒子的密度、速度、温度、矩量以及位置等信息，根据确定的矩量，利用方程(22)确定该时刻平均粒径沿空间分布；

(7)判定SDPH计算的时间t是否超出FVM时间步长ΔT_FVM，假如未超出，则转到步骤(3)循环计算，直至计算总时间大于ΔT_FVM；

(8)重置SDPH计算的时间t＝0，更新两相的体积分数值，将SDPH粒子信息同样采用核函数反插值到网格中心位置处，开始FVM程序计算，SDPH等待；

(9)FVM利用插值到网格位置处的信息和自身网格数据，计算气相受到液滴相的曳力、热传导量，作为源项添加到FVM输运方程中，同时结合SDPH计算得到的液滴相粒径分布及各相体积分数值，计算FVM质量、动量、能量和湍流输运方程；

(10)更新气相速度、温度、压力、湍流耗散等参量，判定FVM计算是否收敛，收敛则更新FVM计算总时间，未收敛则转到步骤(8)循环计算，直至收敛；

(11)由FVM计算总时间判定是否完成计算，未完成计算则转到步骤(3)更新以粒子形式存储的网格信息，开始SDPH计算，若完成计算则停止所有程序计算。

(12)计算过程中输出不同时刻液滴群和气流场的数据，通过后处理软件显示获得不同时刻液滴群的空间分布、液滴大小、液滴速度、气流场的速度等信息，预测发动机燃烧室内气体-液滴两相流动特性，为发动机部件的优化设计以及复杂问题的解决提供数据支撑。

步骤105，若所述第一判断结果表示否，则基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟；

步骤106，基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟；

步骤107，基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴蒸发与燃烧过程进行模拟。

步骤107所述基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴蒸发与燃烧过程进行模拟，具体包括：

将以粒子形式存储的网格中心点的信息采用核函数插值的方式插值到液滴粒子点位置，获得液滴粒子点位置的网格中心点信息；网格中心点信息包括液滴的位置、速度、体积分数、气场压力、湍流耗散量；

根据液滴粒子点位置的网格中心点信息，利用公式

计算单液滴蒸发量；其中，m_p为液滴质量，N_i为组分i的液滴的数量，A_p为液滴的表面积，M_w,i为组分i的摩尔质量，d_p为液滴直径，α_p为液滴的体积分数；

根据单液滴蒸发量更新液滴的直径、质量、体积分数和光滑离散粒子密度，利用大密度差多相流粒子方法和完全变光滑长度法，求解液滴粒子的密度、速度、光滑长度和位移的单位时间变化量；

根据液滴粒子的密度、速度、光滑长度和位移的单位时间变化量采用蛙跳方式，更新液滴粒子的质量、速度和能量；

将液滴粒子点位置的网格中心点信息和单液滴蒸发量采用核函数反插值的方式插值到网格中心位置处；

根据单液滴蒸发量计算每个网格的蒸发量和体积分数；

利用气相燃烧模型

计算化学反应的平均速率

其中，ω_t为基于k-ε湍流模型计算的湍流化学反应速率，ω_t＝min(ω_t1,ω_t2)；ω_t1表示已燃气体微团破碎速率，ω_t2表示未燃气体微团破碎速率；ω_t1表示已燃气体微团破碎速率，ω_t2表示未燃气体微团破碎速率，

v_i为反应物i的化学计量数；v_R为生成物R的化学计量数；M_i为反应物i的摩尔质量；M_R为生成物R的摩尔质量；A、B为经验常数，A＝4.0，B＝0.5，ε为湍流耗散率，k为湍流强度，ρ为液滴的密度，m_R为生成物的质量分数，m_i为反应物i的质量分数，N为总的反应物的数目。

根据每个网格的蒸发量、体积分数和化学反应的平均速率，采用求解组分输运方程和N-S方程的方式求解网格的质量、速度和能量的单位时间变化量；具体的，利用组分运输方程获得每个网格上的组分的分布数据，将这些数据作为N-S方程中连续性方程的源项，计算N-S方程，获得每个网格质量、速度和能量的时间变化量，进而采用蛙跳时间更新获得新的时刻每个网格的质量、速度和能量。

根据质量、速度和能量的单位时间变化量，采用蛙跳时间更新获得新的时刻每个网格的质量、速度和能量。

返回步骤“将以粒子形式存储的网格中心点的信息采用核函数插值的方式插值到液滴粒子点位置，获得液滴粒子点位置的网格中心点信息”，当总的计算时间大于设定的终止时间时，计算终止。

可选的，所述利用大密度差多相流粒子方法和完全变光滑长度法，求解液滴粒子的密度、速度、光滑长度和位移的单位时间变化量，具体包括：

大密度差多相流粒子方法的公式为：

完全变光滑长度方法的公式为：

其中，ρ_i表示组分i和组分j的密度，m_i和m_j分别表示组分i和组分j的密度，V_i和V_j分别表示组分i组分j的体积，v_i表示组分i的速度，v_ij表示组分i与组分j的相对速度，x_i表示组分i的位置，W_ij表示粒子i和粒子j之间的核函数数值，P_i和P_j分别表示粒子i和粒子j的压力值，μ_i和μ_j分别表示粒子i和粒子j的动力粘度值，r_ij表示粒子i和粒子j之间的距离，

表示粒子i和粒子j之间的平均压力，

表示粒子i和粒子j之间的人工应力值，Π_ij表示粒子i和粒子j之间的人工粘性值，h_i和h_j分别表示组分i和组分j的能量。

可选的，所述组分输运方程为：

其中，

为组分i的化学反应的平均速率，S_i为离散相及自定义的源项产生的额外物质质量生成速率，J_i为组分i的扩散通量；Y_i表示组分i的主流蒸汽浓度，ρ和v分别表示液滴的密度和速度；

ρ为液滴密度；

D_i,m为蒸汽扩散系数；

Sc为施密特数；

μ为液滴的动力粘度系数；

Y_i表示组分i的主流蒸汽浓度。

发动机燃烧室内液滴的燃烧主要是以扩散燃烧为主，即液滴从外界吸收热量，自身发生相变，液体燃料转变为燃料蒸气向周围扩散，与外界氧化剂接触，达到燃烧条件后，发生化学反应。因此，将这里的液滴的燃烧考虑为扩散燃烧，包括液滴蒸发和气相燃烧两个主要过程，要在光滑离散粒子法-有限体积耦合方法中实现这两个过程，具体流程如图6所示：

1)当液滴温度达到或超过相变临界值T_vap后，光滑离散粒子所表征的液滴进行蒸发反应，所采用的蒸发模型和组分输运模型如下：

采用液滴蒸发传质定律计算液滴蒸发过程，并在达到沸腾温度之前并且可蒸发部分完全析出之前一直保持蒸发状态，即

T_vap≤T_p≤T_bp且m_p＞(1-f_v,0)m_p,0 (23)

本发明中假设液滴可完全蒸发，即f_v,0＝1.0。液滴的沸腾温度相对较高，液滴无法达到该温度值。液滴的蒸发量由梯度扩散值确定，即，液滴向气相中的扩散率和液滴与气流主流之间的蒸汽浓度梯度相关联

N_i＝k_i(C_i,p-C_i,g) (24)

式中，N_i(kg·mol^-1/m²s^-1)为蒸汽的摩尔流率，k_c(m/s)为传质系数，C_i,p(kg·mol/m³)为液滴表面的蒸汽浓度，C_i,g(kg·mol/m³)为气相主流的蒸汽浓度。如果N_i为正值，表征液滴处于蒸发状态，如果N_i为负值，则表征液滴处于凝结状态，此时，将该液滴视为惯性液滴，将N_i调整为0。液滴表面的蒸汽分压假定为液滴温度T_p所对应的饱和压力p_sat，而此时的蒸汽浓度为对应此分压的浓度：

对于第i个组分，主流蒸汽浓度由组分输运方程求解得到，式(26)

其中R_i为化学反应的净生成速率，S_i为离散相及自定义的源项产生的额外物质质量生成速率。J_i为组分i的扩散通量，取决于物质的浓度梯度。对于气体混合物进行组分输运方程的计算，气相的计算采用有限体积方法完成。

对于非预混及部分预混燃烧，蒸汽浓度通过查找概率密度函数表格得到：

X_i为第i个组分的当地体积摩尔分数，p_op为操作压强，T_g为当地气体体积平均温度。

公式(25)中的传质系数由努塞尔关联式求得：

其中D_i,m(m²/s)为蒸汽扩散系数，Sc为施密特数，Sc＝μ/ρD_i,m，d_p为液滴直径，由此得液滴的质量消耗为：

M_w,i(kg/kgmol)为i组分的摩尔质量，m_p(kg)为液滴质量，A_p(m²)为液滴的表面积。

由于液滴蒸发会造成相应的动量和能量的变化，由此得到的变化值作为源项施加到相应的动量方程和能量方程中，动量的变化量为

能量的变化量为

h_fg(J/kg)为液滴的汽化潜热。

2)液滴蒸发后，SDPH粒子的质量减小，质量减小值为公式(29)，动量和能量同样发生相应改变(公式(30)和(31))。而这些变化量反作用于有限体积网格上。在计算相间传输量的同时更新SDPH粒子所表征的液滴的粒径分布，根据质量减小量公式(29)和液滴密度值获得新的液滴粒径分布；更新液滴相以及气体相的体积分数值，液滴相的体积分数值为SDPH粒子的密度与液滴相密度的比值，气体相的体积分数与液滴相体积分数值和为1，由此获得液滴相体积分数值。

3)SDPH离散粒子属性发生改变后会直接造成计算的不稳定，尤其是液滴蒸发造成光滑离散粒子密度减小，流场密度不均匀，传统的方法计算易发散，引入大密度差多相光滑离散粒子方法解决该问题，方程组如式(32)

针对密度改变后，光滑长度改变的问题，采用完全变光滑长度解决，方程组如式(33)

3)根据单液滴的蒸发量(公式(29))，以及FVM网格中所包含的液滴情况，计算出每个网格所获取的液滴蒸发物质含量，进而更新计算气体相物质组分输运方程(公式(26))，得到各物质的空间分布。

4)网格上的燃料蒸汽在高温火源的作用下，开始燃烧化学反应，计算EBU-Arrhenius气相湍流燃烧模型，化学反应的平均速率仅取决于低温的反应物和高温的燃烧产物之间的湍流混合作用，而与化学反应动力学无关。平均反应速率由ω_t和ω_A两者中较小的量决定，如下式

其中，ω_t为基于k-ε湍流模型计算的湍流化学反应速率，由已燃和未燃气体微团破碎速率中的较小值决定，两者分别为

ω_A为基于Arrhenius公式计算的平均化学反应速率，如下式

燃烧消耗燃料和氧气，产成新的生成物，直至反应物燃烧耗尽。

步骤108，返回步骤“判断液滴温度是否小于液滴的沸点温度，获得第一判断结果”，继续下一个时间步的模拟，直到模拟时间达到预设模拟时间。

本发明还提供一种发动机燃烧室内气体-液滴两相流动特性的数值模拟系统，所述数值模拟系统包括：

几何模型判断模块，用于构建发动机燃烧室的几何模型。

第一判断模块，用于判断液滴温度是否小于液滴的沸点温度，获得第一判断结果。

第一液滴碰撞过程模拟模块，用于若所述第一判断结果表示是，则基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟。

第一液滴剪切破碎过程模拟模块，用于基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟。

第二液滴碰撞过程模拟模块，用于若所述第一判断结果表示否，则基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟。

第二液滴剪切破碎过程模拟模块，用于基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟。

液滴蒸发与燃烧过程模拟模块，用于基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴蒸发与燃烧过程进行模拟。

返回模块，用于返回步骤“判断液滴温度是否小于液滴的沸点温度，获得第一判断结果”，继续下一个时间步的模拟，直到模拟时间达到预设模拟时间。

所述第一液滴碰撞过程模拟模块，具体包括：液滴分布矩初始化子模块，用于给定初始的液滴分布矩；权值计算子模块，用于根据上一时刻的液滴分布矩，利用公式

计算每个积分点的当前时刻的权值；其中，

表示在积分点m处的液滴的粒径的k阶矩，ω_k(t)表示当前时刻t的k阶矩的权值，m_k(t-1)表示前一时刻的液滴尺度分布函数的k阶矩，m＝0,1,2,…,N-1,k＝0,1,2,…,N-1，N表示积分点的数量；矩量输运方程的聚合与破碎源项计算子模块，用于根据当前时刻的权值，利用公式

计算矩量输运方程的聚合与破碎源项

其中，ω_α(t)和ω_γ(t)分别表示当前时刻t的积分点α和积分点γ的权值，L_α和L_γ分别表示在积分点α和积分点γ处的液滴的粒径，β_αγ表示粒径为L_α,L_γ的两液滴碰撞的概率；

表示在积分点m处的液滴的粒径的k阶矩；

表示在积分点α处的液滴破碎的概率，

表示在积分点α处的液滴与其它液滴碰撞的概率；

表示在积分点b处的液滴破碎的概率；液滴分布矩求解子模块，用于根据矩量输运方程的聚合与破碎源项求解矩量输运方程

获得当前时刻的液滴分布矩；其中，

为k阶矩的速度；分布数据求解子模块，用于基于当前时刻的液滴分布矩，根据液滴分布矩与液滴数量、液滴分布的标准差和粒径平均值的对应关系，求解液滴数量、标准差和粒径平均值。

本发明获得了不同时刻发动机燃烧室内燃料组分分布信息、燃烧后热量的空间分布信息、燃烧产物的空间分布信息，通过数据后处理直接获得整个发动机燃烧室的燃烧性能，为发动机燃烧室的部件设计和性能改善提供支撑。

(1)本发明提供了一种新的用于预测发动机燃烧室内两相流动与燃烧特性的技术，采用该技术可直接获得实际发动机燃烧室内的流动状态和化学反应状态(实际实验很难测量得到)，从而预估出发动机燃烧室的燃烧性能(点、熄火性能、温度场分布、燃烧不稳定等)，为解决实际发动机燃烧室设计研发过程中存在的难点问题提供强有力的支撑。

(2)本发明提出了计算液滴碰撞和液滴剪切破碎的新思路，并将新的思路在无网格粒子方法上进行了实现。不仅可对大量液滴表现出的宏观特性进行描述，同时又对液滴间相互作用的微观行为进行了刻画；不仅可以进行轨道追踪，同时计算量又得到了控制，实现了一条由微观机理分析到宏观拟流体模型建立再到微观粒子仿真实现的新路径，为液滴和气泡这类特殊液滴问题的计算提供了一种新方法。

(3)本发明解决了包括液滴粒径改变和相间传质耦合计算等科学问题，克服了现有方法在燃烧室两相流计算方面遇到的难题，不仅计算量大幅降低，而且计算精度也将有明显提高，为有效认识发动机燃烧室内的气体-液滴两相流动规律，揭示不同物质和过程之间的相互作用机理，解决发动机燃烧室问题提供了有效的计算工具。

(4)本发明一方面突破了适用于大量液滴求解的光滑离散粒子方法，突破了传统光滑粒子法仅用于离散连续性物质的现状：基于液滴动力学理论，通过增加光滑粒子所表征的液滴相的材料属性，推导光滑粒子属性与液滴属性间的关系式，将光滑粒子改造成适于离散相求解的光滑离散粒子法，建立光滑粒子与真实液滴间的一一对应关系。该研究突破了传统光滑粒子方法材料属性固定，粒子仅作为几何质点的局限，拓展至具有物理质点特性的范畴，为后续进一步拓展光滑粒子方法的应用提供了方向。

(5)本发明建立了光滑离散粒子法-有限体积耦合新方法，为气体-液滴两相流动问题的求解提供了一条新途径：基于双流体模型提出了光滑离散粒子法-有限体积耦合新方法，新方法思想不仅适用于光滑离散粒子法与有限体积之间的耦合，对于其他相类似的拉格朗日粒子法与欧拉网格法之间的耦合同样适用，为气体-液滴两相流问题的研究提供了一类更加有效的新的数值方法。同时，液滴动力学模型与光滑离散粒子方法相结合求解液滴相的思想，实现了离散相由微观因素综合分析到宏观拟流体模型建立再到微观粒子方法求解的模拟过程，为工业乃至自然界中所有具有离散流动特性的问题提供了一条新的有效求解途径。

(6)本发明建立了考虑液滴与气相复杂化学反应及含液滴碰撞聚合、破碎过程的新型光滑离散粒子法-有限体积耦合方法：通过引入液滴的蒸发模型、气相燃烧模型以及描述液滴间碰撞聚合与破碎等微观行为的群体平衡模型，建立了考虑液滴与气相复杂化学反应及含液滴碰撞聚合、破碎过程的新型光滑离散粒子法-有限体积耦合方法。新方法解决了传统液滴轨道模型和双流体模型在求解该问题时遇到的难题。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。

本文中应用了具体个例对发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

Claims (6)

1.一种发动机燃烧室内气体-液滴两相流动特性的数值模拟方法，其特征在于，所述数值模拟方法包括如下步骤：

构建发动机燃烧室的几何模型；

判断液滴温度是否小于液滴的沸点温度，获得第一判断结果；

若所述第一判断结果表示是，则基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟；

基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟；

若所述第一判断结果表示否，则基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟；

基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟；

基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴蒸发与燃烧过程进行模拟；

返回步骤“判断液滴温度是否小于液滴的沸点温度，获得第一判断结果”，继续下一个时间步的模拟，直到模拟时间达到预设模拟时间；

所述基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟，具体包括：

给定初始的液滴分布矩；

根据上一时刻的液滴分布矩，利用公式

计算每个积分点的当前时刻的权值；

其中，

表示在积分点m处的液滴的粒径的k阶矩，ω_k(t)表示当前时刻t的k阶矩的权值，m_k(t-1)表示前一时刻的液滴尺度分布函数的k阶矩，m＝0,1,2,…,N-1,k＝0,1,2,…,N-1，N表示积分点的数量；

根据当前时刻的权值，利用公式

计算矩量输运方程的聚合与破碎源项

其中，ω_α(t)和ω_γ(t)分别表示当前时刻t的积分点α和积分点γ的权值，L_α和L_γ分别表示在积分点α和积分点γ处的液滴的粒径，β_αγ表示粒径为L_α,L_γ的两液滴碰撞的概率；

表示在积分点m处的液滴的粒径的k阶矩；

表示在积分点α处的液滴破碎的概率，

表示在积分点α处的液滴与其它液滴碰撞的概率；

表示在积分点b处的液滴破碎的概率；

根据矩量输运方程的聚合与破碎源项求解矩量输运方程

获得当前时刻的液滴分布矩；

其中，

为k阶矩的速度；

基于当前时刻的液滴分布矩，根据液滴分布矩与液滴数量、液滴分布的标准差和粒径平均值的对应关系，求解液滴数量、标准差和粒径平均值。

2.根据权利要求1所述的发动机燃烧室内气体-液滴两相流动特性的数值模拟方法，其特征在于，所述采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟，具体包括：

将以粒子形式存储的网格中心点的信息采用核函数插值的方式插值到SDPH粒子点位置，获得SDPH粒子点位置的网格中心点信息；所述网格中心点信息包括液滴的位置、速度、体积分数、气场压力、湍流耗散量；

根据SDPH粒子点位置的两相体积分数、气相速度和湍流耗散值，求解群体平衡方程

获得不同时刻液滴的空间状态特征；

其中，n(L,t)为t时刻尺寸为L的液滴的数量，

为阶矩的平均速度，B(L,t)和D(L,t)分别表示t时刻尺寸为L的液滴由于聚合和破碎造成的生函数和死函数；

根据SDPH粒子点位置的网格中心点信息，计算SDPH粒子点位置所受到的气相曳力、气相压力和热传导量；

将SDPH粒子点位置所受到的气相曳力、气相压力和热传导量作为矩量输运方程的聚合与破碎源项，并根据不同时刻液滴的空间状态特征，求解矩量输运方程获得SDPH粒子的粒径的单位时间变化量，求解能量方程获得SDPH粒子的质量、动量和能量的单位时间变化量；

根据SDPH粒子的粒径、质量、动量和能量的单位时间变化量，采用蛙跳法更新SDPH粒子点位置的网格中心点信息；

判断当前仿真时长是否小于FVM时间步长，获得第二判断结果；

若所述第二判断结果表示是，则返回步骤“根据SDPH粒子点位置的两相体积分数、气相速度和湍流耗散值，求解群体平衡方程

获得不同时刻液滴的空间状态特征”；

若所述第二判断结果表示否，则将SDPH粒子点位置的网格中心点信息采用核函数反插值的方式插值到网格中心位置处；

根据网格中心位置处的信息和网格数据，计算气相受到的液滴相曳力和热传导量，作为FVM输运方程的源项；

采用求解FVM输运方程的方式，求解FVM液滴的质量、动量和能量的单位时间变化量；

根据FVM液滴的质量、动量和能量的单位时间变化量，采用蛙跳法更新网格中心位置处的信息；

判断FVM是否收敛，获得第三判断结果；

若所述第三判断结果表示否，则返回步骤“根据网格中心位置处的信息和网格数据，计算气相受到的液滴相曳力和热传导量，作为FVM输运方程的源项”；

判断FVM计算总时间是否超过总时间阈值，获得第四判断结果；

若所述第四判断结果表示否，则将SDPH计算时间设置为0，返回步骤“将以粒子形式存储的网格中心点的信息采用核函数插值的方式插值到SDPH粒子点位置，获得SDPH粒子点位置的网格中心点信息”；

若所述第三判断结果表示是，则停止仿真。

3.根据权利要求1所述的发动机燃烧室内气体-液滴两相流动特性的数值模拟方法，其特征在于，所述基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴蒸发与燃烧过程进行模拟，具体包括：

将以粒子形式存储的网格中心点的信息采用核函数插值的方式插值到液滴粒子点位置，获得液滴粒子点位置的网格中心点信息；网格中心点信息包括液滴的位置、速度、体积分数、气场压力、湍流耗散量；

根据液滴粒子点位置的网格中心点信息，利用公式

计算单液滴蒸发量；其中，m_p为液滴质量，N_i为组分i的液滴的数量，A_p为液滴的表面积，M_w,i为组分i的摩尔质量，d_p为液滴直径，α_p为液滴的体积分数；

根据单液滴蒸发量更新液滴的直径、质量、体积分数和光滑离散粒子密度，利用大密度差多相流粒子方法和完全变光滑长度法，求解液滴粒子的密度、速度、光滑长度和位移的单位时间变化量；

根据液滴粒子的密度、速度、光滑长度和位移的单位时间变化量采用蛙跳方式，更新液滴粒子的质量、速度和能量；

将液滴粒子点位置的网格中心点信息和单液滴蒸发量采用核函数反插值的方式插值到网格中心位置处；

根据单液滴蒸发量计算每个网格的蒸发量和体积分数；

利用气相燃烧模型

计算化学反应的平均速率

其中，ω_t为基于k-ε湍流模型计算的湍流化学反应速率，ω_t＝min(ω_t1,ω_t2)；ω_t1表示已燃气体微团破碎速率，ω_t2表示未燃气体微团破碎速率，ω_A表示基于阿仑尼乌斯模型计算的化学反应速率；

根据每个网格的蒸发量、体积分数和化学反应的平均速率，采用求解组分输运方程和N-S方程的方式求解网格的质量、速度和能量的单位时间变化量；

根据质量、速度和能量的单位时间变化量，采用蛙跳时间更新获得新的时刻每个网格的质量、速度和能量；

返回步骤“将以粒子形式存储的网格中心点的信息采用核函数插值的方式插值到液滴粒子点位置，获得液滴粒子点位置的网格中心点信息”。

4.根据权利要求3所述的发动机燃烧室内气体-液滴两相流动特性的数值模拟方法，其特征在于，利用大密度差多相流粒子方法和完全变光滑长度法，求解液滴粒子的密度、速度、光滑长度和位移的单位时间变化量，具体包括：

大密度差多相流粒子方法的公式为：

完全变光滑长度方法的公式为：

其中，ρ_i表示组分i和组分j的密度，m_i和m_j分别表示组分i和组分j的密度，V_i和V_j分别表示组分i组分j的体积，v_i表示组分i的速度，v_ij表示组分i与组分j的相对速度，x_i表示组分i的位置，W_ij表示粒子i和粒子j之间的核函数数值，

表示W_ij的梯度，P_i和P_j分别表示粒子i和粒子j的压力值，μ_i和μ_j分别表示粒子i和粒子j的动力粘度值，r_ij表示粒子i和粒子j之间的距离，

表示粒子i和粒子j之间的平均压力，

表示粒子i和粒子j之间的人工应力值，Π_ij表示粒子i和粒子j之间的人工粘性值，h_i和h_j分别表示组分i和组分j的能量，f_ij表示粒子i与粒子j之间的人工应力系数。

5.根据权利要求3所述的发动机燃烧室内气体-液滴两相流动特性的数值模拟方法，其特征在于，所述组分输运方程为：

其中，

为组分i的化学反应的平均速率，S_i为离散相及自定义的源项产生的额外物质质量生成速率，J_i为组分i的扩散通量；Y_i表示组分i的主流蒸汽浓度，ρ和v分别表示液滴的密度和速度。

6.一种发动机燃烧室内气体-液滴两相流动特性的数值模拟系统，其特征在于，所述数值模拟系统包括：

几何模型判断模块，用于构建发动机燃烧室的几何模型；

第一判断模块，用于判断液滴温度是否小于液滴的沸点温度，获得第一判断结果；

第一液滴碰撞过程模拟模块，用于若所述第一判断结果表示是，则基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟；

第一液滴剪切破碎过程模拟模块，用于基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟；

第二液滴碰撞过程模拟模块，用于若所述第一判断结果表示否，则基于所述几何模型，采用光滑离散粒子法对液滴碰撞过程进行模拟；

第二液滴剪切破碎过程模拟模块，用于基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴剪切破碎过程进行模拟；

液滴蒸发与燃烧过程模拟模块，用于基于所述几何模型，采用光滑离散粒子法和有限体积耦合法对液滴蒸发与燃烧过程进行模拟；

返回模块，用于返回步骤“判断液滴温度是否小于液滴的沸点温度，获得第一判断结果”，继续下一个时间步的模拟，直到模拟时间达到预设模拟时间；

所述第一液滴碰撞过程模拟模块，具体包括：

液滴分布矩初始化子模块，用于给定初始的液滴分布矩；

权值计算子模块，用于根据上一时刻的液滴分布矩，利用公式

计算每个积分点的当前时刻的权值；

其中，

表示在积分点m处的液滴的粒径的k阶矩，ω_k(t)表示当前时刻t的k阶矩的权值，m_k(t-1)表示前一时刻的液滴尺度分布函数的k阶矩，m＝0,1,2,…,N-1,k＝0,1,2,…,N-1，N表示积分点的数量；

矩量输运方程的聚合与破碎源项计算子模块，用于根据当前时刻的权值，利用公式

计算矩量输运方程的聚合与破碎源项

其中，ω_α(t)和ω_γ(t)分别表示当前时刻t的积分点α和积分点γ的权值，L_α和L_γ分别表示在积分点α和积分点γ处的液滴的粒径，β_αγ表示粒径为L_α,L_γ的两液滴碰撞的概率；

表示在积分点m处的液滴的粒径的k阶矩；

表示在积分点α处的液滴破碎的概率，

表示在积分点α处的液滴与其它液滴碰撞的概率；

表示在积分点b处的液滴破碎的概率；

液滴分布矩求解子模块，用于根据矩量输运方程的聚合与破碎源项求解矩量输运方程

获得当前时刻的液滴分布矩；

其中，

为k阶矩的速度；

分布数据求解子模块，用于基于当前时刻的液滴分布矩，根据液滴分布矩与液滴数量、液滴分布的标准差和粒径平均值的对应关系，求解液滴数量、标准差和粒径平均值。

Images