CN112069689A - 一种航空发动机燃油雾化特性的仿真方法及系统 - Google Patents

一种航空发动机燃油雾化特性的仿真方法及系统 Download PDF

Info

Publication number
CN112069689A
CN112069689A CN202010944300.0A CN202010944300A CN112069689A CN 112069689 A CN112069689 A CN 112069689A CN 202010944300 A CN202010944300 A CN 202010944300A CN 112069689 A CN112069689 A CN 112069689A
Authority
CN
China
Prior art keywords
particle
ith
density
jth
calculating
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010944300.0A
Other languages
English (en)
Other versions
CN112069689B (zh
Inventor
陈福振
严红
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xi'an Keli Simulation Software Technology Co ltd
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN202010944300.0A priority Critical patent/CN112069689B/zh
Publication of CN112069689A publication Critical patent/CN112069689A/zh
Application granted granted Critical
Publication of CN112069689B publication Critical patent/CN112069689B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种航空发动机燃油雾化特性的仿真方法，首先基于光滑粒子流体动力学方法，建立用于进行航空发动机燃油雾化特性预测的物理模型的离散方程组；然后基于离散方程组进行仿真，根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长、每个粒子的表面张力的单元体积力和壁面处边界力；利用离散方程组，计算每个粒子的密度、速度和位置的单位时间变化量；根据单位时间变化量更新每个粒子的密度、速度和位置。本发明提出了用于航空发动机燃油雾化特性预测的无网格粒子仿真方法，该方法具有计算量小、稳定性高、液体属性可调整、液滴轨迹可追踪等优势，同时具有较好的实用性和可拓展性。

Description

一种航空发动机燃油雾化特性的仿真方法及系统
技术领域
本发明涉及燃油雾化特性研究技术领域，特别涉及一种航空发动机燃油雾化特性的仿真方法及系统。
背景技术
为满足新一代军用飞机对战术机动性、短距起飞、超音速巡航等优异作战性能以及民用飞机对低成本、高清洁、高可靠性的要求，未来航空发动机必然要向高推重比、高压比、高温升、低污染趋势发展，这对燃烧室内的燃油燃烧效率、稳定工作范围、出口温度分布、耗油率、火焰筒冷却以及污染物排放等提出了更为严苛的标准。燃油雾化作为燃烧过程的初始阶段，其性能的好坏对发动机的燃烧效率和燃烧稳定性有重大影响，研究燃油的雾化机理和喷注单元的雾化特性对发动机设计具有重要意义。
早期，航空发动机燃烧室主要以压力雾化喷嘴为主，其存在点火容易、燃烧稳定范围宽等优点，但该雾化模式易引起燃烧室起烟、积碳和火焰筒壁温过高等问题。因此，发展了空气雾化喷嘴，其在发动机工作中排放量低、燃烧效率高，但是燃烧稳定范围较窄、低工况下燃烧效率低、点火困难，满足不了高性能航机要求。这种情况下，组合式空气雾化喷嘴由于兼有压力雾化喷嘴和空气雾化喷嘴的优点，逐渐成为当今航空业发达国家正努力探讨解决的关键技术之一。然而，无论是空气雾化喷嘴还是组合式空气雾化喷嘴，其燃油雾化明显不同于其他动力装置雾化过程，它会受到强湍流扰动、多级旋流和高压强等复杂气动条件的影响，液膜的破碎、流动、蒸发、掺混等过程呈现高度非定常特性与复杂相界面结构，雾化问题逐渐成为航发发动机燃烧技术的重要问题之一。
近年来，研究人员从实验、理论及数值模拟三个方面对航空发动机内部的燃油雾化进行了大量深入系统的研究，其中，实验研究是燃油雾化研究的主要手段，通过改变实验工况、捕捉不同的细节和测量不同的参量，获得雾化实验结果，分析燃油雾化特性；燃油雾化理论目前尚不完善，还主要停留在简单雾化条件下的单一雾化机理揭示，对于复杂气动条件下的航空发动机雾化来说，理论预测和实验结果之间还存在一定的偏差；数值模拟作为实验和理论研究重要的补充，近些年逐渐引起了学者的广泛关注，尤其随着高性能计算平台的发展和新型计算方法的完善，应用数值模拟方法对航空发动机燃油雾化进行了一定的探讨，取得了一些有益的研究成果，但毕竟属于起步阶段，还需进一步快速发展。
根据对燃油相处理方式的不同，目前国内外对燃油雾化过程的数值模拟主要采用两种方法：一种是将燃油设定为连续相，采用基于网格的数值模拟方法，对燃油雾化过程中液膜的形成、波纹的震荡、液滴的剥离以及液丝的破碎等过程进行精确的捕捉，追踪雾化过程中液面的变形破裂过程。另一种研究燃油雾化过程的方法是将燃油直接设定为离散相，在发动机喷嘴入口处，直接给定燃油液滴的初始分布，以颗粒态加入到燃烧室流场中，采用颗粒轨道模型对液滴的运动、变形、碰撞、破碎、蒸发等过程进行计算。
第一种方法中，主要是采用基于流体体积函数的界面追踪方法，它的基本思想是，采用体积率函数C(x)描述网格单元中某一流体所占的体积率，当网格单元完全由某一种流体占据，则C＝1，当网格单元中完全不含该流体，则C＝0，对于两种流体的边界有0<C<1。流体体积函数方法是一种基于固定网格的表面跟踪技术。该模型用于观察两种及以上互不相融流体间的分界面。流体体积函数模型中，两种流体共用一组动量方程，计算域中各流体的体积分数在每个计算单元上被跟踪。方法中，首先通过求解一组动量方程来模拟两种及两种以上互不相融的流体，在整个计算域跟踪各流体的体积分数；然后对给定的界面流体体积函数进行界面重构；最后再计算下一时刻各网格单元中流体体积函数的值。
第二种方法中，颗粒轨道追踪方法是基于欧拉-拉格朗日框架所构建，该方法认为流体相是连续的，液滴相是离散的。根据气体-液滴两相系统中液滴运动的特点，该方法将两相流动中液滴的运动过程分解为受冲力支配的瞬时碰撞运动及受流体曳力控制的悬浮运动，从而建立了液滴运动分解模型。该方法中，流体的运动规律用连续介质的N-S方程进行描述，而液滴的行为则通过在拉格朗日坐标系中分析每一个单液滴的运动轨迹而进行描述。其中，在液滴与液滴相互作用的过程中，其运动规律服从碰撞动力学中的动量守恒定律；在流体与液滴相互作用的悬浮过程中，液滴在曳力、重力等力的作用下，其运动规律服从牛顿动力学中的力平衡方程。这样，每个液滴的速度及位移的更新由邻近液滴对它的碰撞作用及流体对它的悬浮作用来确定。与此同时，液滴对流体的瞬时作用则反映在离散相与流体相两相耦合的N-S方程不断进行修正的数值求解过程中。
第一种界面追踪方法，是基于欧拉网格方法，通过对网格内的流体体积分数值求解，再结合界面定位技术获得实际中的界面位置，虽然可以捕捉界面变形破碎过程的细节，但是随着液膜到液丝再到液滴甚至更小液滴的发展，过大的网格无法捕捉到气液两相的界面，需要采用网格自适应技术不断加密网格，计算量巨大，通常计算单股射流雾化过程便需要五千核以上高性能计算机计算200小时以上，对于实际复杂的航空发动机燃油雾化装置来说，采用该方法进行数值模拟难度巨大。
第二种颗粒轨道追踪方法，忽略了液膜变形破碎细节，直接将燃油雾化后的液滴分布加入到流场中进行计算，虽然适用于大量液滴的追踪模拟，可计算蒸发、燃烧过程，但是无法深入认识燃油雾化的机理，燃油雾化过程中的细节无法有效获取。
实现减少航空发动机燃油雾化特性模拟的计算量，并实现雾化过程的细节的模拟成为一个亟待解决的技术问题。
发明内容
本发明的目的是提供一种航空发动机燃油雾化特性的仿真方法及系统，以减少航空发动机燃油雾化特性模拟的计算量，并实现雾化过程的细节的模拟。
为实现上述目的，本发明提供了如下方案：
一种航空发动机燃油雾化特性的仿真方法，所述仿真方法包括如下步骤：
建立压力旋流雾化喷嘴结构的粒子模型；所述粒子模型包括I个粒子；
建立描述压力旋流雾化喷嘴雾化过程中燃油和气体的气液两相流物理模型；
对所述气液两相流物理模型进行离散化，获得离散方程组；
根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长；
根据上一个仿真周期的每个粒子的密度、速度和位置，计算当前仿真周期的每个粒子的表面张力的单元体积力和壁面处边界力；
根据当前仿真周期的每个粒子的表面张力的单元体积力和壁面处边界力，利用离散方程组，计算每个粒子的密度、速度和位置的单位时间变化量；
根据每个粒子的密度、速度和位置的单位时间变化量和当前仿真周期的时间步长，更新每个粒子的密度、速度和位置，返回步骤“根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长”，等待进行下一个仿真周期的仿真。
可选的，所述建立压力旋流雾化喷嘴结构的粒子模型，具体包括：
采用三维制图软件建立压力旋流雾化喷嘴结构的几何模型；
利用Hypermesh软件对所述几何模型进行网络划分，获得网格划分后的几何模型；
对所述网格划分后的几何模型进行网格到粒子的转化处理，获得压力旋流雾化喷嘴结构的粒子模型。
可选的，所述建立描述压力旋流雾化喷嘴雾化过程中燃油和气体的气液两相流物理模型，具体包括：
建立描述压力旋流雾化喷嘴雾化过程中燃油和气体的气液两相流物理模型为：
Figure BDA0002674700860000051
其中，ρi表示第i个粒子的密度，vi表示第i个粒子的速度，Pi为第i个粒子的压力，μi为第i个粒子的动力粘度，D表示变形张量，g为重力加速度，Fis为第i个粒子的表面张力的单元体积力，fi bp为第i个粒子的壁面处边界力；xi表示第i个粒子的位置。
可选的，所述对所述气液两相流物理模型进行离散化，获得离散方程组，具体包括：
对所述气液两相流物理模型进行离散化，获得离散方程组为：
Figure BDA0002674700860000052
其中，mi表示第i个粒子的质量，Vi和Vj分别表示第i个粒子和第j个粒子的体积；
Figure BDA0002674700860000053
表示第i个粒子与第j个粒子的均值压力，
Figure BDA0002674700860000054
ρi和ρj分别表示第i个粒子和第j个粒子的密度，Pi和Pj分别表示第i个粒子和第j个粒子的压力；
Figure BDA0002674700860000055
表示
Figure BDA0002674700860000056
Wij为第i个粒子和第j个粒子之间的核函数的数值，Wij＝W(xj-xi,hi)，W(·)为核函数，hi为第i个粒子的光滑长度；μi和μj分别表示第i个粒子和第j个粒子的动力粘度，xi和xj分别表示第i个粒子和第j个粒子的位置；rij为第i个粒子和第j个粒子之间的位移；vij为第i个粒子和第j个粒子之间的速度差矢量vij＝vi-vj，vi和vj分别表示第i个粒子和第j个粒子的速度，Fis为第i个粒子的表面张力的单元体积力，fi bp为第i个粒子的壁面处边界力，g为重力加速度。
可选的，所述根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长，具体包括：
利用公式
Figure BDA0002674700860000061
计算第一时间步长Δt1；其中，ci表示第i个粒子的色标；
利用公式
Figure BDA0002674700860000062
计算第二时间步长Δt2；
利用公式
Figure BDA0002674700860000063
计算第三时间步长Δt3；
利用公式
Figure BDA0002674700860000064
计算第四时间步长Δt4；
利用公式Δt＝min(Δt1,Δt2,Δt3,Δt4)，计算当前仿真周期的时间步长Δt；
其中，hi为第i个粒子的光滑长度，fi为第i个粒子所受到的体积力，αΠ和βΠ分别为第一无量纲参量和第二无量纲参量，φij为第i个粒子和第j个粒子之间的粘性，ρi表示第i个粒子的密度，μi为第i个粒子的动力粘度，σ表示表面张力系数，I表示粒子的数量。
可选的，根据上一个仿真周期的每个粒子的密度、速度和位置，计算当前仿真周期的每个粒子的表面张力的单元体积力，具体包括：
根据上一个仿真周期的每个粒子的密度、速度和位置，利用公式
Figure BDA0002674700860000065
计算每个粒子的表面位置；
其中，
Figure BDA0002674700860000066
表示第i个粒子的表面位置，cj表示第j个粒子的色标，mj和ρj分别表示第j个粒子的质量和密度，Wij为第i个粒子和第j个粒子之间的核函数的数值；
根据每个粒子的表面位置
Figure BDA0002674700860000071
，利用公式
Figure BDA0002674700860000072
计算每个粒子的法向量ni＝(nxi,nyi,nzi)；
其中，ni＝(nxi,nyi,nzi)表示第i个粒子的法向量，nαi表示第i个粒子的法向量的α轴分量，α表示坐标方向，α＝x,y,z；nxi、nyi和nzi分别表示第i个粒子的法向量的x轴分量、y轴分量和z轴分量，Wij,x、Wij,y、Wij,z分别表示第i个粒子和第j个粒子之间的光滑核函数值Wij在x、y、z方向的偏导数；
Figure BDA0002674700860000073
和
Figure BDA0002674700860000074
分别表示第i个粒子和第j个粒子的表面位置的α轴分量，
Figure BDA0002674700860000075
和
Figure BDA0002674700860000076
分别表示第i个粒子的位置的x轴分量、y轴分量和z轴分量，
Figure BDA0002674700860000077
和
Figure BDA0002674700860000078
分别表示第j个粒子的位置的x轴分量、y轴分量和z轴分量；mj和ρj分别表示第j个粒子的质量和密度；N表示粒子的数量；
根据每个粒子的法向量，利用公式
Figure BDA0002674700860000079
计算每个粒子的曲率；
其中，ki表示第i个粒子的曲率，
Figure BDA00026747008600000710
表示第i个粒子的法向量的单位向量，
Figure BDA00026747008600000711
和
Figure BDA00026747008600000712
分别表示第i个粒子的法向量的x轴分量、y轴分量和z轴分量的单位向量，
Figure BDA00026747008600000713
和
Figure BDA00026747008600000714
分别表示第i个粒子的法向量的x轴方向分量相对于x轴的单位向量、y轴方向分量相对于y轴的单位向量和z轴方向分量相对于z轴的单位向量，
Figure BDA00026747008600000715
Figure BDA00026747008600000716
为第i个粒子的法向量的α轴方向分量相对于γ轴的单位向量，γ＝x,y,z，
Figure BDA00026747008600000717
和
Figure BDA00026747008600000718
分别表示第i个粒子和第j个粒子的法向量的γ轴分量的单位向量；
根据第i个粒子的法向量和曲率，利用公式Fis＝σkiδsni，计算第i个粒子的表面张力的单元体积力Fis；
其中，σ为表面张力系数，δs为表面狄拉克函数。
可选的，所述根据上一个仿真周期的每个粒子的密度、速度和位置，计算当前仿真周期的每个粒子的壁面处边界力，具体包括：
利用公式
Figure BDA0002674700860000081
计算第i个粒子的壁面处边界力；
其中，fi bp为第i个粒子的壁面处边界力，nj为第j个粒子的法向量，ε为罚参数，B为包含固体边界粒子的空间，vi和vj分别表示第i个粒子和第j个粒子的速度，Aj指第j个粒子与壁面的接触面积，hi是第i个粒子的光滑长度，Wij第i个粒子和第j个粒子之间的光滑核函数值，rij是第i个粒子和第j个粒子之间的位移矢量。
可选的，所述根据每个粒子的密度、速度和位置的单位时间变化量和当前仿真周期的时间步长，更新每个粒子的密度、速度和位置，具体包括：
根据每个粒子的密度、速度和位置的单位时间变化量和当前仿真周期的时间步长，利用公式
Figure BDA0002674700860000082
更新每个粒子的密度、速度和位置；
其中，
Figure BDA0002674700860000083
和
Figure BDA0002674700860000084
分别表示更新后的第i个粒子的密度、速度和位置，
Figure BDA0002674700860000085
和
Figure BDA0002674700860000086
分别表示更新前的第i个粒子的密度、速度和位置，
Figure BDA0002674700860000087
和
Figure BDA0002674700860000088
表示第i个粒子的密度和速度的单位时间变化量，Δt表示当前仿真周期的仿真时长。
一种航空发动机燃油雾化特性的仿真系统，所述仿真系统包括如下步骤：
粒子模型建立模块，用于建立压力旋流雾化喷嘴结构的粒子模型；
气液两相流物理模型建立模块，用于建立描述压力旋流雾化喷嘴雾化过程中燃油和气体的气液两相流物理模型；
离散化模块，用于对所述气液两相流物理模型进行离散化，获得离散方程组；
时间步长计算模块，用于根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长；
单元体积力和壁面处边界力计算模块，用于根据上一个仿真周期的每个粒子的密度、速度和位置，计算当前仿真周期的每个粒子的表面张力的单元体积力和壁面处边界力；
单位时间变化量变化模块，用于根据当前仿真周期的每个粒子的表面张力的单元体积力和壁面处边界力，利用离散方程组，计算每个粒子的密度、速度和位置的单位时间变化量；
仿真更新模块，用于根据每个粒子的密度、速度和位置的单位时间变化量和当前仿真周期的时间步长，更新每个粒子的密度、速度和位置，返回步骤“根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长”，等待进行下一个仿真周期的仿真。
可选的，所述粒子模型建立模块，具体包括：
几何模型建立子模块，用于采用三维制图软件建立压力旋流雾化喷嘴结构的几何模型；
网格划分子模块，用于利用Hypermesh软件对所述几何模型进行网络划分，获得网格划分后的几何模型；
粒子转化子模块，用于对所述网格划分后的几何模型进行网格到粒子的转化处理，获得压力旋流雾化喷嘴结构的粒子模型。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明公开了一种航空发动机燃油雾化特性的仿真方法，首先基于光滑粒子流体动力学方法，建立用于进行航空发动机燃油雾化特性预测的物理模型的离散方程组；然后基于离散方程组进行仿真，根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长、每个粒子的表面张力的单元体积力和壁面处边界力；利用离散方程组，计算每个粒子的密度、速度和位置的单位时间变化量；根据单位时间变化量更新每个粒子的密度、速度和位置。本发明提出了用于航空发动机燃油雾化特性预测的无网格粒子仿真方法，该方法具有计算量小、稳定性高、液体属性可调整、液滴轨迹可追踪等优势，同时具有较好的实用性和可拓展性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明提供的一种航空发动机燃油雾化特性的仿真方法的流程图；
图2为本发明提供的压力旋流雾化喷嘴结构的粒子模型的结构图；
图3为本发明提供的燃油雾化过程中液膜、液丝和液滴的形态变化示意图；其中，图3(a)为燃油雾化过程的侧面视图，图3(b)为燃油雾化过程的斜底部视图；
图4为本发明提供的燃油速度矢量分布变化过程示意图；其中，图4(a)为0.018ms时燃油速度矢量分布图，图4(b)为0.042ms时燃油速度矢量分布图，图4(c)为0.206ms时燃油速度矢量分布图；
图5为本发明提供的液滴粒径分布的概率密度图；
图6为本发明提供的不同压力的SMD曲线图；
图7为本发明提供的不同粘度下液滴的SMD曲线图。
具体实施方式
本发明的目的是提供一种航空发动机燃油雾化特性的仿真方法及系统，以减少航空发动机燃油雾化特性模拟的计算量，并实现雾化过程的细节的模拟。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对发明作进一步详细的说明。
本发明技术方案基于无网格粒子仿真技术(光滑粒子流体动力学方法，SmoothedParticle Hydrodynamics，SPH)模拟航空燃油雾化过程，该技术方法的设计过程如下：(1)航空发动机燃油雾化喷嘴几何模型建立：根据实际燃油雾化喷嘴几何构型，建立喷嘴的三维几何模型，再基于SPH三维高保真技术建立雾化喷嘴的粒子离散模型；(2)燃油-气体两相流物理模型建立：建立气液两相流动物理模型、气液界面张力模型、牛顿流体本构模型；(3)气体-燃油两相材料参数选取：对雾化过程涉及到的气体和燃油的物性参数进行选取和确定；(4)SPH对燃油-气体两相流物理模型的离散：采用SPH数值离散技术对所建立的两相流物理模型进行离散，获得离散方程组；(5)壁面条件选取和施加：选择燃油与喷嘴壁面之间的相互作用力施加方式，建立壁面力施加公式；(6)时间积分格式建立：确定SPH非稳态时间步更新所采用的时间积分格式；(7)数值仿真计算：计算燃油经由管道进入喷嘴内，在喷嘴内部流动，经喷嘴出口喷出形成液膜、液丝再到液滴的整个过程；(8)计算结果后处理：得到燃油雾化过程中液膜、液丝和液滴的形态变化、液滴空间分布、喷嘴内燃油流动规律、燃油速度矢量分布、液滴粒径分布规律等；(9)结果分析和机理揭示：对燃油雾化过程进行分析，揭示气液两相界面演变的物理机理；对不同参数影响下的雾化结果进行分析，获得燃油雾化规律。
如图1所示本发明提供一种航空发动机燃油雾化特性的仿真方法，所述仿真方法包括如下步骤：
步骤101，建立压力旋流雾化喷嘴结构的粒子模型；所述粒子模型包括I个粒子。
步骤101所述建立压力旋流雾化喷嘴结构的粒子模型，具体包括：采用三维制图软件建立压力旋流雾化喷嘴结构的几何模型；利用Hypermesh软件对所述几何模型进行网络划分，获得网格划分后的几何模型；对所述网格划分后的几何模型进行网格到粒子的转化处理，获得压力旋流雾化喷嘴结构的粒子模型。
具体的，燃油雾化喷嘴模型的建立采用商业软件建立三维几何模型，再导入网格划分软件，进行细致均匀的网格划分，最后网格文件导入到程序中按照一个网格转化成一个粒子的转化原则，全部转化为粒子，具体过程为：1)首先建立实际航空发动机雾化装置几何模型：采用SolidWorks软件完成，该软件是达索公司旗下的SolidWorks子公司负责开发的三维软件，无论多复杂的几何模型，都可分解成有限数量的几何构型特征，而每一种构成特征，都可用有限的参数完全约束。同时，基于特征的实体模型化系统，研究人员采用基于智能特性的功能去生产模型，如腔、壳、倒角及圆角，可以任意勾画草图，较易改变模型，这一功能特性给工程设计者提供了设计上的简便和灵活。2)在步骤1)的基础上，对于步骤1)建立的几何模型进行网格划分：对于模型网格生成，采用功能强大的CAE应用软件包—Hypermesh软件完成，该软件具有强大的网格划分处理功能，它使CAE分析工程师可以在高度交互和可视化环境下进行仿真分析工作，与其他的有限元前后处理器比较，Hypermesh的图形用户界面容易理解和学习，特别是它支持直接输入已有的三维CAD几何模型如ProE所建的几何模型，并且导入的效率和生成的模型质量都很高，可以大大减少很多重复性的工作，使得CAE工程师能够投入更多精力和时间到分析计算工作上去。3)将步骤2)划分后形成的网格文件导入到自编程序中进行网格到粒子的转化：按照一个网格对应一个粒子的原则，采用任意六面体体积计算方法计算六面体网格的体积即为SPH粒子的体积，采用任意六面体质心计算方法计算六面体网格的质心即为SPH粒子的质心也就是SPH粒子的初始位置，这样就获得了燃油雾化喷嘴结构的粒子模型。SPH单个粒子的体积决定了SPH单个粒子的质量，SPH粒子的质心，也就是SPH粒子的初始位置，直接决定了物质的初始位置。图2为示例所建的压力旋流雾化喷嘴结构的粒子模型。
步骤102，建立描述压力旋流雾化喷嘴雾化过程中燃油和气体的气液两相流物理模型。
步骤102所述建立描述压力旋流雾化喷嘴雾化过程中燃油和气体的气液两相流物理模型，具体包括：
建立描述压力旋流雾化喷嘴雾化过程中燃油和气体的气液两相流物理模型为：
Figure BDA0002674700860000131
其中，ρi表示第i个粒子的密度，vi表示第i个粒子的速度，Pi为第i个粒子的压力，μi为第i个粒子的动力粘度，D表示变形张量，g为重力加速度，Fis为第i个粒子的表面张力的单元体积力，fi bp为第i个粒子的壁面处边界力；xi表示第i个粒子的位置。
具体的，根据燃油雾化过程中气液两相运动过程，建立非稳态不可压缩Navier-Stokes方程，表面张力作为源项添加在Navier-Stokes方程中，具体控制方程如下：
Figure BDA0002674700860000132
ρ≡ρ(x,t)为燃油液体的密度，u≡(u,v,w)为燃油液体的速度，P为压力，μ≡μ(x,t)为燃油液体的动力粘度，
Figure BDA0002674700860000133
为变形张量，σ为表面张力系数，κ为两相界面的曲率，n为两相界面的法向量，δs表示表面张力仅作用于两相界面处，fbp为壁面对气液两相流体的作用力。
对于表面张力模型，采用基于连续表面力模型(continuum surface force，CSF)的表面张力计算方法，CSF模型将表面张力描述为通过界面的连续的作用力，界面厚度有限，在界面内色函数连续地变化。CSF表面张力模型的思想是从定义色函数出发，通过色函数计算得到表面的法向及曲率，在此基础上求得单元表面力并转换为单元体积力，期间保持转化的动量守恒。
在表面张力系数为常量的条件下，在有限的界面厚度范围内，单元体积力Fs表示为
Fs＝fsδs (4)
δs为表面狄拉克函数，fs为单元表面力，通过下式计算
Figure BDA0002674700860000141
其中σ为表面张力系数，k(x)为界面x处的曲率，
Figure BDA0002674700860000142
为界面的单位法向。法向n可由下式计算
Figure BDA0002674700860000143
则单位法向表示为
Figure BDA0002674700860000144
其中c(x)为针对不同相流体定义的色函数。曲率k为法向的散度，即
Figure BDA0002674700860000145
针对在气液固多物质交界面处，引入Brackbill提出的壁面附着力边界条件处理方法，对在气液固三相交界处的流体粒子及部分固体边界虚粒子的界面法向进行修正，修正前后保持法向模值不变，修正公式如式(9)
Figure BDA0002674700860000146
式中
Figure BDA0002674700860000147
为修正后单位化法向，始终沿液体与壁面接触达到稳定状态时气液交界面的法线方向，
Figure BDA0002674700860000148
为在边界平面内并沿接触线法线方向的单位化法向，
Figure BDA0002674700860000149
为垂直于
Figure BDA00026747008600001410
并指向壁面内的单位化法向，θ为液体与壁面形成的接触角。
对气液两相流物理模型中涉及到的物性参数进行选取，气体密度ρg＝1.228kg/m3、粘度ηg＝1.8×10-5Pa·s，航空燃油液体密度ρl＝780kg/m3、粘度ηl＝3.0×10-3Pa·s，气液界面的表面张力为0.0758N/m。
步骤103，对所述气液两相流物理模型进行离散化，获得离散方程组。
步骤103所述对所述气液两相流物理模型进行离散化，获得离散方程组，具体包括：
对所述气液两相流物理模型进行离散化，获得离散方程组为：
Figure BDA0002674700860000151
其中，mi表示第i个粒子的质量，Vi和Vj分别表示第i个粒子和第j个粒子的体积；
Figure BDA0002674700860000152
表示第i个粒子与第j个粒子的均值压力，
Figure BDA0002674700860000153
ρi和ρj分别表示第i个粒子和第j个粒子的密度，Pi和Pj分别表示第i个粒子和第j个粒子的压力；
Figure BDA0002674700860000154
表示
Figure BDA0002674700860000155
Wij为第i个粒子和第j个粒子之间的核函数的数值，Wij＝W(xj-xi,hi)，W(·)为核函数，hi为第i个粒子的光滑长度；μi和μj分别表示第i个粒子和第j个粒子的动力粘度，xi和xj分别表示第i个粒子和第j个粒子的位置；rij为第i个粒子和第j个粒子之间的位移；vij为第i个粒子和第j个粒子之间的速度差矢量vij＝vi-vj，vi和vj分别表示第i个粒子和第j个粒子的速度，Fis为第i个粒子的表面张力的单元体积力，fi bp为第i个粒子的壁面处边界力，g为重力加速度。
步骤104，根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长。
步骤104所述根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长，具体包括：利用公式
Figure BDA0002674700860000156
计算第一时间步长Δt1；其中，ci表示第i个粒子的色标；
利用公式
Figure BDA0002674700860000157
计算第二时间步长Δt2；
利用公式
Figure BDA0002674700860000161
计算第三时间步长Δt3；
利用公式
Figure BDA0002674700860000162
计算第四时间步长Δt4；
利用公式Δt＝min(Δt1,Δt2,Δt3,Δt4)，计算当前仿真周期的时间步长Δt；
其中，hi为第i个粒子的光滑长度，fi为第i个粒子所受到的体积力，αΠ和βΠ分别为第一无量纲参量和第二无量纲参量，φij为第i个粒子和第j个粒子之间的粘性，ρi表示第i个粒子的密度，μi为第i个粒子的动力粘度，σ表示表面张力系数，I表示粒子的数量。
具体的，对于蛙跳积分，时间步长必须满足稳定性条件，这里应用柯朗－弗里德里奇－列维(Courant-Friedrich-Lewy，简称CFL)条件对时间步长进行估计，具体表达式为：
Figure BDA0002674700860000163
Figure BDA0002674700860000164
Figure BDA0002674700860000165
Figure BDA0002674700860000166
其中，f为作用于单位质量上的外力，μ为流体的动力粘度，σ为界面表面张力系数，最终取式(22)-(25)中最小值作为SPH计算的时间步长。
步骤105，根据上一个仿真周期的每个粒子的密度、速度和位置，计算当前仿真周期的每个粒子的表面张力的单元体积力和壁面处边界力；
步骤105中的根据上一个仿真周期的每个粒子的密度、速度和位置，计算当前仿真周期的每个粒子的表面张力的单元体积力，具体包括：
根据上一个仿真周期的每个粒子的密度、速度和位置，利用公式
Figure BDA0002674700860000167
计算每个粒子的表面位置；
其中，
Figure BDA0002674700860000171
表示第i个粒子的表面位置，cj表示第j个粒子的色标，mj和ρj分别表示第j个粒子的质量和密度，Wij为第i个粒子和第j个粒子之间的核函数的数值；
根据每个粒子的表面位置
Figure BDA0002674700860000172
，利用公式
Figure BDA0002674700860000173
计算每个粒子的法向量ni＝(nxi,nyi,nzi)；
其中，ni＝(nxi,nyi,nzi)表示第i个粒子的法向量，nαi表示第i个粒子的法向量的α轴分量，α表示坐标方向，α＝x,y,z；nxi、nyi和nzi分别表示第i个粒子的法向量的x轴分量、y轴分量和z轴分量，Wij,x、Wij,y、Wij,z分别表示第i个粒子和第j个粒子之间的光滑核函数值Wij在x、y、z方向的偏导数；
Figure BDA0002674700860000174
和
Figure BDA0002674700860000175
分别表示第i个粒子和第j个粒子的表面位置的α轴分量，
Figure BDA0002674700860000176
和
Figure BDA0002674700860000177
分别表示第i个粒子的位置的x轴分量、y轴分量和z轴分量，
Figure BDA0002674700860000178
和
Figure BDA0002674700860000179
分别表示第j个粒子的位置的x轴分量、y轴分量和z轴分量；mj和ρj分别表示第j个粒子的质量和密度；N表示粒子的数量；
根据每个粒子的法向量，利用公式
Figure BDA00026747008600001710
计算每个粒子的曲率；
其中，ki表示第i个粒子的曲率，
Figure BDA00026747008600001711
表示第i个粒子的法向量的单位向量，
Figure BDA00026747008600001712
和
Figure BDA00026747008600001713
分别表示第i个粒子的法向量的x轴分量、y轴分量和z轴分量的单位向量，
Figure BDA00026747008600001714
和
Figure BDA00026747008600001715
分别表示第i个粒子的法向量的x轴方向分量相对于x轴的单位向量、y轴方向分量相对于y轴的单位向量和z轴方向分量相对于z轴的单位向量，
Figure BDA00026747008600001716
Figure BDA00026747008600001717
为第i个粒子的法向量的α轴方向分量相对于γ轴的单位向量，γ＝x,y,z，
Figure BDA00026747008600001718
和
Figure BDA00026747008600001719
分别表示第i个粒子和第j个粒子的法向量的γ轴分量的单位向量；
根据第i个粒子的法向量和曲率，利用公式Fis＝σkiδsni，计算第i个粒子的表面张力的单元体积力Fis；
其中，σ为表面张力系数，δs为表面狄拉克函数。
具体的，首先是表面定位公式
Figure BDA0002674700860000181
cj是粒子j的色标，在定义的流体区域内初始设为1，在流体区域外时初始设为0。
其次是法向n的通用计算式为
Figure BDA0002674700860000182
采用CSPM方法对式(12)进行修正，其核心思想是采用基于Taylor级数展开的校正核估计代替传统方法中的核估计来离散控制方程组，修正后的法向分量计算式(三维)为
Figure BDA0002674700860000183
其中α、β取值为1或2，表示坐标方向，nαi表示粒子i在α方向的法向分量，
Figure BDA0002674700860000184
由公式(11)得出，Wij＝W(xj-xi,h)，
Figure BDA0002674700860000185
此方法在处理尖角等粒子缺失严重的边界问题时将得到相比式(12)精度更高的结果。
然后是进行曲率的计算，曲率即法向的散度，传统计算式
Figure BDA0002674700860000186
采用CSPM方法对曲率进行修正，修正后的曲率分量计算式(三维)为
Figure BDA0002674700860000191
Figure BDA0002674700860000192
为粒子i的法向分量
Figure BDA0002674700860000193
在α方向的偏导数，α、β、γ取值为1或2，表示坐标方向。
Figure BDA0002674700860000194
为粒子i,j在γ方向的法向分量，由公式(13)、(16)、(17)得出
Figure BDA0002674700860000195
Figure BDA0002674700860000196
将公式(15)计算得到的
Figure BDA0002674700860000197
代入曲率计算公式(18)，计算获得曲率值ki
Figure BDA0002674700860000198
有了以上公式(13)和公式(18)计算获得的法向值和曲率值，则根据公式(5)和(4)获得气液两相界面的表面张力数值，再进一步代入Navier-Stokes方程(10)中计算获得流场分布。
步骤105中的根据上一个仿真周期的每个粒子的密度、速度和位置，计算当前仿真周期的每个粒子的壁面处边界力，具体包括：
利用公式
Figure BDA0002674700860000199
计算第i个粒子的壁面处边界力；
其中，fi bp为第i个粒子的壁面处边界力，nj为第j个粒子的法向量，ε为罚参数，B为包含固体边界粒子的空间，vi和vj分别表示第i个粒子和第j个粒子的速度，Aj指第j个粒子与壁面的接触面积，hi是第i个粒子的光滑长度，Wij第i个粒子和第j个粒子之间的光滑核函数值，rij是第i个粒子和第j个粒子之间的位移矢量。
具体的，公式(1)中的壁面边界作用力fbp的计算，需要选取适当的壁面力计算方法，本发明对于壁面边界条件采用罚函数法来对SPH方法施加接触边界条件，计算得到第i个粒子所受到的壁面力fi bp
Figure BDA0002674700860000201
fi bp为壁面处粒子i所受边界力，ε为罚参数，B为包含固体边界粒子的空间，B表示流体和固体交界面，Aj指边界粒子j与壁面的接触面积。
步骤106，根据当前仿真周期的每个粒子的表面张力的单元体积力和壁面处边界力，利用离散方程组，计算每个粒子的密度、速度和位置的单位时间变化量。
有了以上公式(13)和公式(18)计算获得的法向值和曲率值，则根据公式(5)和(4)获得气液两相界面的表面张力数值，再进一步代入Navier-Stokes方程(10)中计算获得流场分布。
步骤107，根据每个粒子的密度、速度和位置的单位时间变化量和当前仿真周期的时间步长，更新每个粒子的密度、速度和位置，返回步骤“根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长”，等待进行下一个仿真周期的仿真。
步骤107所述根据每个粒子的密度、速度和位置的变化量更新每个粒子的密度、速度和位置，具体包括：
根据每个粒子的密度、速度和位置的单位时间变化量和当前仿真周期的时间步长，利用公式
Figure BDA0002674700860000202
更新每个粒子的密度、速度和位置；
其中，
Figure BDA0002674700860000203
和
Figure BDA0002674700860000204
分别表示更新后的第i个粒子的密度、速度和位置，
Figure BDA0002674700860000211
和
Figure BDA0002674700860000212
分别表示更新前的第i个粒子的密度、速度和位置，
Figure BDA0002674700860000213
和
Figure BDA0002674700860000214
表示第i个粒子的密度和速度的单位时间变化量，Δt表示当前仿真周期的仿真时长。
即，利用公式(20)和(21)进行更新。
Figure BDA0002674700860000215
xi(t+δt)＝xi(t)+vi(t+δt/2)δt (21)
式中
Figure BDA0002674700860000216
表示粒子i的密度ρ，速度v；xi为粒子i处的位置坐标。
Figure BDA0002674700860000217
是通过方程组(10)计算获得的dvi(t)和dρi(t)；δt为采用公式(22)-(25)计算获得的时间步长。
在以上步骤1010-106的基础上，自编程序实现以上算法，然后进行编译，采用OpenMP并行计算的方式进行计算，主要计算燃油经由管道进入喷嘴内，在喷嘴内部流动，而后经喷嘴出口喷出形成液膜、液丝再到液滴的整个过程，获得不同时间节点上的流场数据ρ、v、以及位移xi。
计算结果后处理可采用两种方式完成：一种是采用自主研发的数值模拟软件进行数据的读取和结果的显示；第二种是依托商业软件Tecplot，按照程序控制信息提供的数据输出方式，输出所有场变量，生成相关动画。按照程序控制信息提供的粒子/节点编号和变量类型编号，生成相关变量的时间历程曲线。如图3-7所示展示了tecplot软件处理得到的燃油雾化过程中液膜、液丝和液滴的形态变化、液滴空间分布、喷嘴内燃油流动规律、燃油速度矢量分布、液滴粒径分布规律。
本发明还提供一种航空发动机燃油雾化特性的仿真系统，所述仿真系统包括如下步骤：
粒子模型建立模块，用于建立压力旋流雾化喷嘴结构的粒子模型。
所述粒子模型建立模块，具体包括：几何模型建立子模块，用于采用三维制图软件建立压力旋流雾化喷嘴结构的几何模型；网格划分子模块，用于利用Hypermesh软件对所述几何模型进行网络划分，获得网格划分后的几何模型；粒子转化子模块，用于对所述网格划分后的几何模型进行网格到粒子的转化处理，获得压力旋流雾化喷嘴结构的粒子模型。
气液两相流物理模型建立模块，用于建立描述压力旋流雾化喷嘴雾化过程中燃油和气体的气液两相流物理模型。
离散化模块，用于对所述气液两相流物理模型进行离散化，获得离散方程组。
时间步长计算模块，用于根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长。
单元体积力和壁面处边界力计算模块，用于根据上一个仿真周期的每个粒子的密度、速度和位置，计算当前仿真周期的每个粒子的表面张力的单元体积力和壁面处边界力。
单位时间变化量变化模块，用于根据当前仿真周期的每个粒子的表面张力的单元体积力和壁面处边界力，利用离散方程组，计算每个粒子的密度、速度和位置的单位时间变化量。
仿真更新模块，用于根据每个粒子的密度、速度和位置的单位时间变化量和当前仿真周期的时间步长，更新每个粒子的密度、速度和位置，返回步骤“根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长”，等待进行下一个仿真周期的仿真。
(1)本发明基于无网格粒子仿真方法，考虑两相之间的大密度差、气液固三相交界面处理、气液两相表面张力计算、壁面边界条件等因素，计算过程中无需进行自适应，计算量得到了控制，同时该技术属于拉格朗日方法，对于两相界面的追踪比较自然，无需引入特别的界面追踪技术，为解决气液两相界面问题提供了新的途径。
(2)传统对于航空发动机燃烧室内燃油雾化特性预测大多采用实验和理论预测的方式，实验方法费时费力，同时受实验测量技术的制约很多雾化细节无法捕获；理论预测仅能获得最终雾化特性与初始参数之间的关系，无法有效揭示雾化过程机理；而数值模拟可以较好解决实验和理论计算存在的不足，但传统数值模拟大多采用基于网格的数值模拟技术，仅单股射流就需要近万核超三十天周期，因此无法对于航空发动机复杂雾化装置下的雾化过程进行有效模拟，本发明则独辟蹊径引入无网格粒子仿真技术很好的解决该问题。
(3)本发明突破了传统界面追踪技术必须采用网格进行模拟的现状，网格模拟的过程中必须结合网格自适应才能捕获不同尺度的液滴，而无网格粒子方法无需网格，更无需网格自适应，由粒子数量来表征具体液滴大小；同时克服了传统界面追踪技术必须引入界面定位方法才能将界面精确表示的不足，无网格粒子方法属于拉格朗日方法，对于界面的定位属于自然定位，根据两相界面运动的位置直接获取。
(4)本发明突破了传统颗粒轨道追踪技术无法获取液膜、液丝破碎的一次雾化过程的缺陷，传统颗粒轨道追踪技术仅能从雾化后形成液滴之后开始计算，对于一次雾化过程是忽略的，无法采用该方法去认识雾化过程的细节，而本发明的无网格粒子方法是直接从燃油进入雾化装置开始进行计算，包括燃油在喷嘴内部的流动、燃油在喷嘴出口的积聚、液膜的形成与破碎、液丝的形成和断裂、液滴的形成与运动等均可以有效捕捉，克服了颗粒轨道追踪技术的不足。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。
本文中应用了具体个例对发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

Claims (10)

1.一种航空发动机燃油雾化特性的仿真方法，其特征在于，所述仿真方法包括如下步骤：
建立压力旋流雾化喷嘴结构的粒子模型；所述粒子模型包括I个粒子；
建立描述压力旋流雾化喷嘴雾化过程中燃油和气体的气液两相流物理模型；
对所述气液两相流物理模型进行离散化，获得离散方程组；
根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长；
根据上一个仿真周期的每个粒子的密度、速度和位置，计算当前仿真周期的每个粒子的表面张力的单元体积力和壁面处边界力；
根据当前仿真周期的每个粒子的表面张力的单元体积力和壁面处边界力，利用离散方程组，计算每个粒子的密度、速度和位置的单位时间变化量；
根据每个粒子的密度、速度和位置的单位时间变化量和当前仿真周期的时间步长，更新每个粒子的密度、速度和位置，返回步骤“根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长”，等待进行下一个仿真周期的仿真。
2.根据权利要求1所述的航空发动机燃油雾化特性的仿真方法，其特征在于，所述建立压力旋流雾化喷嘴结构的粒子模型，具体包括：
采用三维制图软件建立压力旋流雾化喷嘴结构的几何模型；
利用Hypermesh软件对所述几何模型进行网络划分，获得网格划分后的几何模型；
对所述网格划分后的几何模型进行网格到粒子的转化处理，获得压力旋流雾化喷嘴结构的粒子模型。
3.根据权利要求1所述的航空发动机燃油雾化特性的仿真方法，其特征在于，所述建立描述压力旋流雾化喷嘴雾化过程中燃油和气体的气液两相流物理模型，具体包括：
建立描述压力旋流雾化喷嘴雾化过程中燃油和气体的气液两相流物理模型为：
Figure FDA0002674700850000021
其中，ρi表示第i个粒子的密度，vi表示第i个粒子的速度，Pi为第i个粒子的压力，μi为第i个粒子的动力粘度，D表示变形张量，g为重力加速度，Fis为第i个粒子的表面张力的单元体积力，
Figure FDA0002674700850000026
为第i个粒子的壁面处边界力；xi表示第i个粒子的位置。
4.根据权利要求1所述的航空发动机燃油雾化特性的仿真方法，其特征在于，所述对所述气液两相流物理模型进行离散化，获得离散方程组，具体包括：
对所述气液两相流物理模型进行离散化，获得离散方程组为：
Figure FDA0002674700850000022
其中，mi表示第i个粒子的质量，Vi和Vj分别表示第i个粒子和第j个粒子的体积；
Figure FDA0002674700850000023
表示第i个粒子与第j个粒子的均值压力，
Figure FDA0002674700850000024
ρi和ρj分别表示第i个粒子和第j个粒子的密度，Pi和Pj分别表示第i个粒子和第j个粒子的压力；
Figure FDA0002674700850000027
表示
Figure FDA0002674700850000025
Wij为第i个粒子和第j个粒子之间的核函数的数值，Wij＝W(xj-xi,hi)，W(·)为核函数，hi为第i个粒子的光滑长度；μi和μj分别表示第i个粒子和第j个粒子的动力粘度，xi和xj分别表示第i个粒子和第j个粒子的位置；rij为第i个粒子和第j个粒子之间的位移；vij为第i个粒子和第j个粒子之间的速度差矢量vij＝vi-vj，vi和vj分别表示第i个粒子和第j个粒子的速度，Fis为第i个粒子的表面张力的单元体积力，
Figure FDA0002674700850000037
为第i个粒子的壁面处边界力，g为重力加速度。
5.根据权利要求1所述的航空发动机燃油雾化特性的仿真方法，其特征在于，所述根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长，具体包括：
利用公式
Figure FDA0002674700850000031
计算第一时间步长Δt1；其中，ci表示第i个粒子的色标；
利用公式
Figure FDA0002674700850000032
计算第二时间步长Δt2；
利用公式
Figure FDA0002674700850000033
计算第三时间步长Δt3；
利用公式
Figure FDA0002674700850000034
计算第四时间步长Δt4；
利用公式Δt＝min(Δt1,Δt2,Δt3,Δt4)，计算当前仿真周期的时间步长Δt；
其中，hi为第i个粒子的光滑长度，fi为第i个粒子所受到的体积力，αΠ和βΠ分别为第一无量纲参量和第二无量纲参量，φij为第i个粒子和第j个粒子之间的粘性，ρi表示第i个粒子的密度，μi为第i个粒子的动力粘度，σ表示表面张力系数，I表示粒子的数量。
6.根据权利要求1所述的航空发动机燃油雾化特性的仿真方法，其特征在于，根据上一个仿真周期的每个粒子的密度、速度和位置，计算当前仿真周期的每个粒子的表面张力的单元体积力，具体包括：
根据上一个仿真周期的每个粒子的密度、速度和位置，利用公式
Figure FDA0002674700850000035
计算每个粒子的表面位置；
其中，
Figure FDA0002674700850000036
表示第i个粒子的表面位置，cj表示第j个粒子的色标，mj和ρj分别表示第j个粒子的质量和密度，Wij为第i个粒子和第j个粒子之间的核函数的数值；
根据每个粒子的表面位置
Figure FDA0002674700850000041
利用公式
Figure FDA0002674700850000042
计算每个粒子的法向量ni＝(nxi,nyi,nzi)；
其中，ni＝(nxi,nyi,nzi)表示第i个粒子的法向量，nαi表示第i个粒子的法向量的α轴分量，α表示坐标方向，α＝x,y,z；nxi、nyi和nzi分别表示第i个粒子的法向量的x轴分量、y轴分量和z轴分量，Wij,x、Wij,y、Wij,z分别表示第i个粒子和第j个粒子之间的光滑核函数值Wij在x、y、z方向的偏导数；
Figure FDA0002674700850000043
和
Figure FDA0002674700850000044
分别表示第i个粒子和第j个粒子的表面位置的α轴分量，
Figure FDA0002674700850000045
和
Figure FDA0002674700850000046
分别表示第i个粒子的位置的x轴分量、y轴分量和z轴分量，
Figure FDA0002674700850000047
和
Figure FDA0002674700850000048
分别表示第j个粒子的位置的x轴分量、y轴分量和z轴分量；mj和ρj分别表示第j个粒子的质量和密度；N表示粒子的数量；
根据每个粒子的法向量，利用公式
Figure FDA0002674700850000049
计算每个粒子的曲率；
其中，ki表示第i个粒子的曲率，
Figure FDA00026747008500000410
表示第i个粒子的法向量的单位向量，
Figure FDA00026747008500000411
和
Figure FDA00026747008500000412
分别表示第i个粒子的法向量的x轴分量、y轴分量和z轴分量的单位向量，
Figure FDA00026747008500000413
和
Figure FDA00026747008500000414
分别表示第i个粒子的法向量的x轴方向分量相对于x轴的单位向量、y轴方向分量相对于y轴的单位向量和z轴方向分量相对于z轴的单位向量，
Figure FDA00026747008500000415
Figure FDA00026747008500000416
为第i个粒子的法向量的α轴方向分量相对于γ轴的单位向量，γ＝x,y,z，
Figure FDA00026747008500000417
和
Figure FDA00026747008500000418
分别表示第i个粒子和第j个粒子的法向量的γ轴分量的单位向量；
根据第i个粒子的法向量和曲率，利用公式Fis＝σkiδsni，计算第i个粒子的表面张力的单元体积力Fis；
其中，σ为表面张力系数，δs为表面狄拉克函数。
7.根据权利要求1所述的航空发动机燃油雾化特性的仿真方法，其特征在于，所述根据上一个仿真周期的每个粒子的密度、速度和位置，计算当前仿真周期的每个粒子的壁面处边界力，具体包括：
利用公式
Figure FDA0002674700850000051
计算第i个粒子的壁面处边界力；
其中，
Figure FDA0002674700850000057
为第i个粒子的壁面处边界力，nj为第j个粒子的法向量，ε为罚参数，B为包含固体边界粒子的空间，vi和vj分别表示第i个粒子和第j个粒子的速度，Aj指第j个粒子与壁面的接触面积，hi是第i个粒子的光滑长度，Wij第i个粒子和第j个粒子之间的光滑核函数值，rij是第i个粒子和第j个粒子之间的位移矢量。
8.根据权利要求7所述的航空发动机燃油雾化特性的仿真方法，其特征在于，所述根据每个粒子的密度、速度和位置的单位时间变化量和当前仿真周期的时间步长，更新每个粒子的密度、速度和位置，具体包括：
根据每个粒子的密度、速度和位置的单位时间变化量和当前仿真周期的时间步长，利用公式
Figure FDA0002674700850000052
更新每个粒子的密度、速度和位置；
其中，
Figure FDA0002674700850000053
和
Figure FDA0002674700850000054
分别表示更新后的第i个粒子的密度、速度和位置，
Figure FDA0002674700850000055
和
Figure FDA0002674700850000056
分别表示更新前的第i个粒子的密度、速度和位置，
Figure FDA0002674700850000061
和
Figure FDA0002674700850000062
表示第i个粒子的密度和速度的单位时间变化量，Δt表示当前仿真周期的仿真时长。
9.一种航空发动机燃油雾化特性的仿真系统，其特征在于，所述仿真系统包括如下步骤：
粒子模型建立模块，用于建立压力旋流雾化喷嘴结构的粒子模型；所述粒子模型包括I个粒子；
气液两相流物理模型建立模块，用于建立描述压力旋流雾化喷嘴雾化过程中燃油和气体的气液两相流物理模型；
离散化模块，用于对所述气液两相流物理模型进行离散化，获得离散方程组；
时间步长计算模块，用于根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长；
单元体积力和壁面处边界力计算模块，用于根据上一个仿真周期的每个粒子的密度、速度和位置，计算当前仿真周期的每个粒子的表面张力的单元体积力和壁面处边界力；
单位时间变化量变化模块，用于根据当前仿真周期的每个粒子的表面张力的单元体积力和壁面处边界力，利用离散方程组，计算每个粒子的密度、速度和位置的单位时间变化量；
仿真更新模块，用于根据每个粒子的密度、速度和位置的单位时间变化量和当前仿真周期的时间步长，更新每个粒子的密度、速度和位置，返回步骤“根据上一个仿真周期的每个粒子的密度和速度，计算当前仿真周期的时间步长”，等待进行下一个仿真周期的仿真。
10.根据权利要求9所述的航空发动机燃油雾化特性的仿真系统，其特征在于，所述粒子模型建立模块，具体包括：
几何模型建立子模块，用于采用三维制图软件建立压力旋流雾化喷嘴结构的几何模型；
网格划分子模块，用于利用Hypermesh软件对所述几何模型进行网络划分，获得网格划分后的几何模型；
粒子转化子模块，用于对所述网格划分后的几何模型进行网格到粒子的转化处理，获得压力旋流雾化喷嘴结构的粒子模型。
CN202010944300.0A 2020-09-10 2020-09-10 一种航空发动机燃油雾化特性的仿真方法及系统 Active CN112069689B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010944300.0A CN112069689B (zh) 2020-09-10 2020-09-10 一种航空发动机燃油雾化特性的仿真方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010944300.0A CN112069689B (zh) 2020-09-10 2020-09-10 一种航空发动机燃油雾化特性的仿真方法及系统

Publications (2)

Publication Number Publication Date
CN112069689A true CN112069689A (zh) 2020-12-11
CN112069689B CN112069689B (zh) 2022-02-08

Family

ID=73663309

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010944300.0A Active CN112069689B (zh) 2020-09-10 2020-09-10 一种航空发动机燃油雾化特性的仿真方法及系统

Country Status (1)

Country Link
CN (1) CN112069689B (zh)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112800700A (zh) * 2021-04-13 2021-05-14 中国空气动力研究与发展中心计算空气动力研究所 低温表面干模态结霜模拟方法、装置、电子设备和介质
CN113343450A (zh) * 2021-05-26 2021-09-03 重庆长安汽车股份有限公司 一种内燃机气缸内单燃油液滴蒸发速率的确定方法
CN114218674A (zh) * 2021-12-16 2022-03-22 西北工业大学太仓长三角研究院 一种用于航空发动机燃油雾化全过程性能预测方法及系统
CN114492147A (zh) * 2022-02-15 2022-05-13 西南石油大学 一种颗粒运动全过程的跟踪方法、系统及可存储介质
CN115186570A (zh) * 2022-07-11 2022-10-14 中国人民解放军国防科技大学 一种低成本超声速液体射流喷注雾化数值仿真方法
CN117195776A (zh) * 2023-10-08 2023-12-08 江苏大学 一种基于ibm-voset-dem的三相流动传热耦合仿真方法及系统
CN117252128A (zh) * 2023-11-17 2023-12-19 中国空气动力研究与发展中心计算空气动力研究所 一种旋流喷嘴雾化过程模拟方法、装置、设备及存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104143027A (zh) * 2014-08-01 2014-11-12 北京理工大学 一种基于sph算法的流体热运动仿真系统
CN105183965A (zh) * 2015-08-27 2015-12-23 中国人民解放军国防科学技术大学 用于预测雾化过程的大涡模拟方法
CN106650064A (zh) * 2016-12-09 2017-05-10 华东师范大学 一种基于粒子模型的凝结现象仿真方法
CN108897902A (zh) * 2018-04-04 2018-11-27 上海大学 喷雾干燥塔中物料蒸发的数值模拟方法
US20190024899A1 (en) * 2017-07-21 2019-01-24 General Electric Company Fuel nozzle for a gas turbine engine

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104143027A (zh) * 2014-08-01 2014-11-12 北京理工大学 一种基于sph算法的流体热运动仿真系统
CN105183965A (zh) * 2015-08-27 2015-12-23 中国人民解放军国防科学技术大学 用于预测雾化过程的大涡模拟方法
CN106650064A (zh) * 2016-12-09 2017-05-10 华东师范大学 一种基于粒子模型的凝结现象仿真方法
US20190024899A1 (en) * 2017-07-21 2019-01-24 General Electric Company Fuel nozzle for a gas turbine engine
CN108897902A (zh) * 2018-04-04 2018-11-27 上海大学 喷雾干燥塔中物料蒸发的数值模拟方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
PEREIRA GERALD G等: "Prediction of fluid flow through and jet formation from a high pressure nozzle using Smoothed Particle Hydrodynamics", 《CHEMICAL ENGINEERING SCIENCE》 *
侯燕: "多喷嘴喷雾冷却实验研究与数值模拟", 《中国博士学位论文全文数据库 工程科技Ⅱ辑》 *
强洪夫等: "基于大密度差多相流SPH方法的二维液滴碰撞数值模拟", 《物理学报》 *

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112800700A (zh) * 2021-04-13 2021-05-14 中国空气动力研究与发展中心计算空气动力研究所 低温表面干模态结霜模拟方法、装置、电子设备和介质
CN112800700B (zh) * 2021-04-13 2021-06-25 中国空气动力研究与发展中心计算空气动力研究所 低温表面干模态结霜模拟方法、装置、电子设备和介质
CN113343450A (zh) * 2021-05-26 2021-09-03 重庆长安汽车股份有限公司 一种内燃机气缸内单燃油液滴蒸发速率的确定方法
CN113343450B (zh) * 2021-05-26 2022-08-09 重庆长安汽车股份有限公司 一种内燃机气缸内单燃油液滴蒸发速率的确定方法
CN114218674A (zh) * 2021-12-16 2022-03-22 西北工业大学太仓长三角研究院 一种用于航空发动机燃油雾化全过程性能预测方法及系统
CN114218674B (zh) * 2021-12-16 2023-11-10 西北工业大学太仓长三角研究院 一种用于航空发动机燃油雾化全过程性能预测方法及系统
CN114492147A (zh) * 2022-02-15 2022-05-13 西南石油大学 一种颗粒运动全过程的跟踪方法、系统及可存储介质
CN115186570A (zh) * 2022-07-11 2022-10-14 中国人民解放军国防科技大学 一种低成本超声速液体射流喷注雾化数值仿真方法
CN117195776A (zh) * 2023-10-08 2023-12-08 江苏大学 一种基于ibm-voset-dem的三相流动传热耦合仿真方法及系统
CN117252128A (zh) * 2023-11-17 2023-12-19 中国空气动力研究与发展中心计算空气动力研究所 一种旋流喷嘴雾化过程模拟方法、装置、设备及存储介质
CN117252128B (zh) * 2023-11-17 2024-01-26 中国空气动力研究与发展中心计算空气动力研究所 一种旋流喷嘴雾化过程模拟方法、装置、设备及存储介质

Also Published As

Publication number Publication date
CN112069689B (zh) 2022-02-08

Similar Documents

Publication Publication Date Title
CN112069689B (zh) 一种航空发动机燃油雾化特性的仿真方法及系统
CN113221473B (zh) 发动机燃烧室内气体-液滴两相流动特性的数值模拟方法
Minier et al. PDF model based on Langevin equation for polydispersed two-phase flows applied to a bluff-body gas-solid flow
Wei The development and application of CFD technology in mechanical engineering
CN114218674B (zh) 一种用于航空发动机燃油雾化全过程性能预测方法及系统
CN108009383A (zh) 一种自然层流短舱外形的确定方法及系统
CN114168796B (zh) 一种建立飞行器高空气动力数据库的方法
CN113792432A (zh) 基于改进型fvm-lbfs方法的流场计算方法
CN112632709A (zh) 一种基于fluent仿真的连续激光推力器工质分析方法
CN109408836A (zh) 利用Boltzmann方程进行流体仿真的方法
Suryaprakash et al. Secondary breakup of drops
Rodriguez et al. Formulation and implementation of inflow/outflow boundary conditions to simulate propulsive effects
Cummings et al. Supersonic, turbulent flow computation and drag optimization for axisymmetric afterbodies
Hao et al. Novel design method for inward-turning inlets with non-uniform inflow
CN104516999A (zh) 基于双混合分数的jp5000超音速火焰喷涂的仿真方法
CN110414141B (zh) 可压流体跨音速流动过程中的液滴雾化三维数值模拟方法
Qiao et al. Numerical simulation of distributed propulsion systems using CFD
Waters et al. Modeling turbulent reactive flow in internal combustion engines with an LES in a semi-implicit/explicit finite element projection method
Zhang et al. Simulation of powder transport in plasma jet via hybrid Lattice Boltzmann method and probabilistic algorithm
Stokes et al. Multi-fidelity Computational Investigations of Hypersonic Shock Wave-Boundary Layer Interactions in a Multi-compression Scramjet Inlet
Macleod et al. Modelling Aircraft Fuel Jettison Using Smoothed Particle Hydrodynamics (SPH)
Cai et al. Parallel Numerical Simulation of Complex Unsteady Multi-Component Three-Dimensional Flow Field of Nonequilibrium Chemical Reaction
Pydakula Narayanan High-fidelity simulations of a rotary bell atomizer with electrohydrodynamic effects
Tong et al. Numerical Simulation of Heat Flow Field on the Outer Surface of the Airframe
Bai Study and analysis of lift to drag ratio performance of supercritical wing based on computational fluid dynamics method

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right

Effective date of registration: 20230417

Address after: Room 10505-2, Building 6, Qidi Science and Technology Park, No. 65 Keji Second Road, High tech Zone, Xi'an City, Shaanxi Province, 710075

Patentee after: Xi'an Automaison Information Technology Co.,Ltd.

Address before: 710000 No. 127 Youyi West Road, Shaanxi, Xi'an

Patentee before: Northwestern Polytechnical University

TR01 Transfer of patent right
TR01 Transfer of patent right

Effective date of registration: 20230727

Address after: Room 031, F1903, 19th Floor, Building 4-A, Xixian Financial Port, Fengdong New City Energy Jinmao District, Xi'an City, Shaanxi Province, 710086

Patentee after: Xi'an Keli Simulation Software Technology Co.,Ltd.

Address before: Room 10505-2, Building 6, Qidi Science and Technology Park, No. 65 Keji Second Road, High tech Zone, Xi'an City, Shaanxi Province, 710075

Patentee before: Xi'an Automaison Information Technology Co.,Ltd.

TR01 Transfer of patent right