CN112562793B - 一种针对燃料爆震燃烧的两步反应模型计算方法 - Google Patents

一种针对燃料爆震燃烧的两步反应模型计算方法 Download PDF

Info

Publication number
CN112562793B
CN112562793B CN202011439758.7A CN202011439758A CN112562793B CN 112562793 B CN112562793 B CN 112562793B CN 202011439758 A CN202011439758 A CN 202011439758A CN 112562793 B CN112562793 B CN 112562793B
Authority
CN
China
Prior art keywords
reaction model
detonation wave
detonation
reaction
step reaction
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.)
Active
Application number
CN202011439758.7A
Other languages
English (en)
Other versions
CN112562793A (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.)
Beijing Institute of Technology BIT
Beijing Power Machinery Institute
Original Assignee
Beijing Institute of Technology BIT
Beijing Power Machinery Institute
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 Beijing Institute of Technology BIT, Beijing Power Machinery Institute filed Critical Beijing Institute of Technology BIT
Priority to CN202011439758.7A priority Critical patent/CN112562793B/zh
Publication of CN112562793A publication Critical patent/CN112562793A/zh
Application granted granted Critical
Publication of CN112562793B publication Critical patent/CN112562793B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/10—Analysis or design of chemical reactions, syntheses or processes
    • 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
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Chemical & Material Sciences (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Analytical Chemistry (AREA)
  • Chemical Kinetics & Catalysis (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Computing Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提供一种针对燃料爆震燃烧的两步反应模型计算方法，包括根据已有基元反应模型计算燃料爆震波关键参数，选取典型关键参数构建燃料爆震燃烧两步诱导‑放热反应模型，采用实验及数值模拟验证所述两步反应模型；通过建模使两步反应模型对应燃料稳定性参数与基元反应模型一致，保证两步反应模型在爆震波稳定性描述方面与复杂基元反应模型一致，可根据实际问题，灵活选择重点关注的爆震燃烧特征参数，使两步反应模型计算得到的特征参数与基元反应模型一致，解决两步反应模型对实际燃料爆震燃烧物理过程模拟的问题，实现两步反应模型在较高爆震燃烧数值模拟效率基础上兼顾模拟精度，可推动两步反应模型在爆震燃烧领域及工程设计领域的广泛应用。

Description

一种针对燃料爆震燃烧的两步反应模型计算方法
技术领域
本发明涉及航空航天技术领域，特别是一种针对燃料爆震燃烧的两步反应模型计算方法。
背景技术
燃烧是燃料化学能转化为热能的主要形式，一般可分为缓燃燃烧和爆震燃烧两种模式。对于缓燃燃烧，燃烧波传播速度为米每秒量级，燃烧过程通常可近似为等压燃烧过程。而对于爆震燃烧，强激波与化学反应紧密耦合在一起，其传播速度可达千米每秒量级，燃烧过程通常可近似为等容燃烧过程。相比于缓燃燃烧，爆震燃烧放热快，产生的熵增较小，具有更高的热效率。因此，与基于缓燃燃烧的传统发动机相比，基于爆震燃烧的爆震发动机具有自增压、熵增低、热循环效率高的优点。
然而针对爆震发动机的试验系统非常复杂，并且由于试验观测手段有限，往往无法获得爆震发动机内部爆震燃烧过程以及爆震波等流场参数分布的详细信息。另外由于试验所需周期长且试验成本和费用昂贵等诸多因素的限制，目前通常采用数值模拟的手段研究爆震发动机内爆震燃烧详细过程，进而弥补上述试验存在的局限性。
目前对于爆震燃烧进行数值模拟所采用的化学动力学模型主要可以分为基元反应模型和分步反应模型。
基元反应模型是化学反应流体力学中最基本、最重要的化学反应模型之一。基元反应模型主要通过若干基本粒子组元之间的化学反应实现，可以描述宏观的燃烧、离解等过程。由于涉及的化学反应较多，相比分步反应更能反映燃料的燃烧特性。通常的基元反应模型是非常复杂的，例如高链烃燃烧的全反应涉及几十种基本粒子和上千个基元反应，即使通过详细化学反应模型简化的反应机制也达上百个反应之多，而且其中很多反应的化学动力学过程目前还不清楚。常见的氢气/氧气简化基元反应模型相对简单，通常涉及到大约十种基本粒子和二十个左右的基元反应，因此这方面数值模拟工作比较多。基元反应模型能够完整系统的描述整个燃烧系统的化学反应过程，但是其较低的模拟效率制约了其应用。在采用基元反应模型的氢氧爆震数值模拟中，通常求解化学反应源项的时间要占整个数值计算时间的80％以上。如果模拟碳氢燃料的爆震燃烧，随着模型的复杂程度增加，越来越多的时间被用来求解源项。然而，这些求解得到的关于流场各部分组元分布的信息并不是我们所需要的，只有极少量的数据是有用的。由于源项的求解耗费了大量的时间，网格密度受到限制从而可能导致影响数值结果的精度。因此，有必要对化学反应模型进行简化，从而提高数值模拟的效率。
分步反应模型根据其对爆震波的不同简化程度，可以分为单步反应模型和两步反应模型等。
单步反应模型通常简单的考虑化学反应过程中热量释放的影响，将爆震燃烧对应反应放热过程通过一个参数(即化学反应进度变量)和能量守恒以及流动方程联系起来。同时化学反应进度变量是通过对化学反应速率的积分得到的，因此建立了流动和化学反应之间的联系。单步反应模型虽然简单，可通过该模型定性的模拟爆震燃烧的一些宏观现象，但是进一步关于燃料爆震燃烧的详细特性和一些定量参数就无法获得。比如爆震波一维稳定性研究表明，数值模拟采用单步反应模型不能得到临界起爆能量、爆震极限等反应爆震波起爆与传播的关键特征参数。
两步反应模型与单步反应模型思想相同，但其在单步反应模型基础上进一步增加了对于爆震波结构的描述，即用两个无量纲化学反应进度变量ξ(1→0)和λ(0→1)分别描述爆震燃烧过程中包含的诱导反应过程和放热反应过程。两个反应进度变量同时也可以看作反应物和燃烧产物的质量分数。当诱导反应没有发生时ξ＝1，随着诱导反应进行，ξ逐渐减小，当诱导反应结束时，即反应气体已全部转化为活性基团(ξ＝0)。同样，在放热反应阶段，用无量纲变量λ来表示放热反应进行的程度。放热反应没有发生时λ＝0，此时没有燃烧产物的产生；当诱导反应进度变量ξ＝0时，放热反应开始进行，并且随着放热反应的进行，λ逐渐增大，当放热反应结束时λ＝1，此时反应气体全部转化为燃烧产物。图1为两步化学反应模型示意图，爆震波从右向左传播。未反应气体经前导激波压缩后，温度和压力都相对较高，未反应气体吸收一部分能量使得气体分子开始活化，诱导反应开始进行，当气体分子活化完成后，放热反应开始进行，大量的热量释放出来，驱动爆震波自持传播。
分步反应模型通过简单的反应模拟爆震燃烧放热过程，计算效率高，通用性强，因此在过去几十年得到了广泛的应用。其中，在早期支链反应模型基础上改进的两步诱导-放热反应模型在爆震波一维稳定性、爆震波胞格理论、斜爆震波起爆区结构及波面稳定性等爆震基础机理研究方面已经取得了比较好的效果，具有广泛的应用前景。然而上述应用两步反应模型对爆震燃烧的研究中所采用的燃料大都不是实际燃料，而是满足研究中给定的两步反应模型参数的假想燃料。这是因为目前尚缺乏针对实际燃料的两步反应模型参数确定方法，即针对实际燃料的两步反应模型建模技术，由此大大限制了两步反应模型在爆震燃烧领域的广泛应用。因此，急需一种针对实际燃料的两步反应模型建模计算新方法。
发明内容
本发明的目的在于提供一种针对燃料爆震燃烧的两步反应模型计算方法，可以实现针对实际燃料的两步反应模型参数确定，及在爆震燃烧领域的广泛应用。
本发明提供的一种针对燃料爆震燃烧的两步反应模型计算方法是这样实现的：
S1、根据基元反应模型确定爆震波关键参数；
S2、根据所述基元反应模型爆震波关键参数，构建两步反应模型；
S3、采用实验及数值模拟验证所述两步反应模型。
进一步的，S1包括以下步骤：
S11、首先采用基元反应模型计算给定初始状态燃料/氧化剂均匀混合气体中C-J(Chapman–Jouguet)爆震波速
Figure GDA0003817505060000031
C-J爆震波马赫数
Figure GDA0003817505060000032
C-J状态参数及von Neumann状态参数。其中C-J状态参数包括：温度
Figure GDA0003817505060000033
压力
Figure GDA0003817505060000034
及密度
Figure GDA0003817505060000035
von Neumann状态参数包括：温度
Figure GDA0003817505060000036
压力
Figure GDA0003817505060000037
密度
Figure GDA0003817505060000038
比热比
Figure GDA0003817505060000039
S12、利用基元反应模型计算混气C-J爆震波对应的ZND(Zel’dovich-von
Figure GDA00038175050600000310
)结构，确定诱导区实际物理长度ΔI，同时计算爆震波一维稳定性参数
Figure GDA00038175050600000311
进一步的，S2两步反应模型为诱导-放热两步反应模型；其使用两个反应速率控制方程分别模拟爆震燃烧的两个过程，第一步诱导反应代表的是诱导区内点火过程，第二步放热反应代表的是放热反应区内热量释放过程，其具体形式如下：
Figure GDA00038175050600000312
Figure GDA00038175050600000313
其中ξ是诱导反应进度变量，λ是放热反应进度变量，H(1-ξ)是一个阶跃函数：
Figure GDA00038175050600000314
kI是诱导区化学反应速率常数，通常取kI＝-uvn；其中uvn是爆震波坐标系下，C-J爆震波中前导激波后气流的速度；kR是放热反应区化学反应速率常数；EI和ER分别是诱导区和放热反应区的活化能，Ts是前导激波后气流的温度；上述物理量均为无量纲量，长度采用的是诱导区长度进行无量纲化，时间采用诱导区长度与参考速度的比值进行无量纲化；
用C-J爆震波前导激波后气流温度对诱导反应和放热反应活化能进行归一化处理，即：
Figure GDA0003817505060000041
进一步的，S2包括以下步骤：
S21、取基元反应模型计算得到的混气C-J爆震波速
Figure GDA0003817505060000042
和von Neumann状态下比热比
Figure GDA0003817505060000043
分别作为两步反应模型的C-J爆震波速VCJ和比热比γ；
S22、从基元反应模型计算得到的C-J状态参数(温度
Figure GDA0003817505060000044
压力
Figure GDA0003817505060000045
密度
Figure GDA0003817505060000046
)、vonNeumann状态参数(温度
Figure GDA0003817505060000047
压力
Figure GDA0003817505060000048
密度
Figure GDA0003817505060000049
)及C-J爆震波马赫数
Figure GDA00038175050600000410
共七个爆震波特征参数中任选其一，作为两步反应模型对应的特征参数，再结合比热比γ和vCJ，根据固定比热比爆震理论，唯一确定两步反应模型的无量纲放热量Q和气体常数R。其中，若选取基元反应模型的C-J爆震波马赫数
Figure GDA00038175050600000411
作为两步反应模型的C-J爆震波马赫数MCJ，则两步反应模型无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure GDA00038175050600000412
Figure GDA00038175050600000413
其中
Figure GDA00038175050600000414
为爆震波波前混气静温。
若选取基元反应模型爆震波C-J状态温度
Figure GDA00038175050600000415
作为两步反应模型爆震波C-J状态温度TCJ，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure GDA00038175050600000416
Figure GDA00038175050600000417
Figure GDA00038175050600000418
其中
Figure GDA00038175050600000419
为爆震波波前混气静温。
若选取基元反应模型爆震波C-J状态压力
Figure GDA00038175050600000420
作为两步反应模型爆震波C-J状态压力PCJ，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure GDA0003817505060000051
Figure GDA0003817505060000052
Figure GDA0003817505060000053
其中
Figure GDA0003817505060000054
为爆震波波前混气静温，
Figure GDA0003817505060000055
为爆震波波前混气静压。
若选取基元反应模型爆震波C-J状态密度
Figure GDA0003817505060000056
作为两步反应模型爆震波C-J状态密度ρCJ，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure GDA0003817505060000057
Figure GDA0003817505060000058
Figure GDA0003817505060000059
其中
Figure GDA00038175050600000510
为爆震波波前混气静温，
Figure GDA00038175050600000511
为爆震波波前混气密度。
若选取基元反应模型爆震波von Neumann状态温度
Figure GDA00038175050600000512
作为两步反应模型爆震波von Neumann状态温度Tvn，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure GDA00038175050600000513
Figure GDA00038175050600000514
Figure GDA00038175050600000515
其中
Figure GDA00038175050600000516
为爆震波波前混气静温。
若选取基元反应模型爆震波von Neumann状态压力
Figure GDA00038175050600000517
作为两步反应模型爆震波von Neumann状态压力Pvn，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure GDA0003817505060000061
Figure GDA0003817505060000062
Figure GDA0003817505060000063
其中
Figure GDA0003817505060000064
为爆震波波前混气静温，
Figure GDA0003817505060000065
为爆震波波前混气静压。
若选取基元反应模型爆震波von Neumann状态密度
Figure GDA0003817505060000066
作为两步反应模型爆震波von Neumann状态密度ρvn，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure GDA0003817505060000067
Figure GDA0003817505060000068
Figure GDA0003817505060000069
其中
Figure GDA00038175050600000610
为爆震波波前混气静温，
Figure GDA00038175050600000611
为爆震波波前混气密度。
S23、根据确定的两步反应模型的无量纲放热量Q和比热比γ，根据固定比热比爆震理论计算两步反应模型对应的C-J状态参数(温度tCJ、压力pCJ、密度ρCJ)和von Neumann状态参数(温度tvn、压力pvn、密度ρvn、速度uvn)；
S24、假设定容爆炸的点火延迟时间τi具有Arrhenius形式：
Figure GDA00038175050600000612
其中T为温度，R为气体常数，EI为活化能，A为指前因子。
给定不同扰动爆震波速度±1％
Figure GDA00038175050600000613
采用基元反应计算经过爆震波前导激波压缩之后气体温度T1和T2，以及定容爆炸对应点火延迟时间τ1和τ2；根据定容爆炸理论计算无量纲诱导反应活化能εI：
Figure GDA00038175050600000614
S25、取诱导区化学反应速率常数kI＝-uvn，同时取定εR＝1.0；通过调整放热反应区化学反应速率常数kR，使得两步反应模型计算得到的爆震波一维稳定性参数χ与基元反应模型计算得到的爆震波一维稳定性参数
Figure GDA0003817505060000071
相等，从而确定放热反应区化学反应速率常数kR。
进一步的，计算爆震波一维稳定性参数χ方法如下：
Figure GDA0003817505060000072
其中εI为诱导反应活化能；ΔI为诱导区长度，定义为前导激波与反应热释放速率最大的位置之间的长度；ΔR放热区长度，定义为：
Figure GDA0003817505060000073
其中uCJ为爆震波坐标系下，C-J状态下当地声速；
Figure GDA0003817505060000074
为化学反应热释放速率，可通过下式计算：
Figure GDA0003817505060000075
其中c为当地声速。
本发明的有益效果：
1、爆震波稳定性参数是描述爆震波稳定性的重要特征参数，本方法通过建模使得两步反应模型对应燃料稳定性参数χ与基元反应模型
Figure GDA0003817505060000076
保持一致，保证了简单的两步反应模型在爆震波的稳定性描述方面与复杂的基元反应模型一致。
2、本建模方法可根据实际研究问题需要，灵活选择七个爆震波特征参数，研究重点关注的爆震燃烧特征参数，使得两步反应模型计算得到的特征参数与基元反应模型一致。
3、本计算方法可以保证两步反应模型对燃料爆震燃烧典型特征参数的准确模拟，解决两步反应模型对实际燃料爆震燃烧物理过程模拟的问题，实现两步反应模型在较高爆震燃烧数值模拟效率基础上兼顾模拟精度，可推动两步反应模型在爆震燃烧领域及工程设计领域的广泛应用。
附图说明
图1是本发明提供的实施例一P0＝1atm,T0＝300K,化学恰当比乙烯/空气两步反应模型与基元反应模型ZND结构对比；
图2是本发明提供的实施例二P0＝1atm,T0＝500K，化学恰当比JP-10/空气两步反应模型与基元反应模型ZND结构对比图；
图3是采用两步反应模型模拟得到的化学恰当比氢气/空气预混气旋转爆震燃烧温度和压力流场；
图4是两步反应模型数值模拟不同时刻爆震波传播速度与基元反应模型计算的C-J速度比值。
具体实施方式
为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本申请保护的范围。
本发明提供了一种针对燃料爆震燃烧的两步反应模型计算方法：
S1、根据基元反应模型确定爆震波关键参数；
S2、根据所述基元反应模型爆震波关键参数，构建两步反应模型；
S3、采用实验及数值模拟验证所述两步反应模型。
进一步的，S1包括以下步骤：
S11、首先采用基元反应模型计算给定初始状态燃料/氧化剂均匀混合气体中C-J爆震波速
Figure GDA0003817505060000081
C-J爆震波马赫数
Figure GDA0003817505060000082
C-J状态参数及von Neumann状态参数。其中C-J状态参数包括：温度
Figure GDA0003817505060000083
压力
Figure GDA0003817505060000084
及密度
Figure GDA0003817505060000085
von Neumann状态参数包括：温度
Figure GDA0003817505060000086
压力
Figure GDA0003817505060000087
密度
Figure GDA0003817505060000088
比热比
Figure GDA0003817505060000089
S12、利用基元反应模型计算混气C-J爆震波对应的ZND(Zel’dovich-von
Figure GDA00038175050600000810
)结构，确定诱导区实际物理长度ΔI，同时计算爆震波一维稳定性参数
Figure GDA00038175050600000811
进一步的，S2两步反应模型为诱导-放热两步反应模型；其使用两个反应速率控制方程分别模拟爆震燃烧的两个过程，第一步诱导反应代表的是诱导区内点火过程，第二步放热反应代表的是放热反应区内热量释放过程，其具体形式如下：
Figure GDA00038175050600000812
Figure GDA00038175050600000813
其中ξ是诱导反应进度变量，λ是放热反应进度变量，H(1-ξ)是一个阶跃函数：
Figure GDA0003817505060000091
kI是诱导区化学反应速率常数，通常取kI＝-uvn；其中uvn是爆震波坐标系下，C-J(Chapman–Jouguet)爆震波中前导激波后气流的速度；kR是放热反应区化学反应速率常数；EI和ER分别是诱导区和放热反应区的活化能，Ts是前导激波后气流的温度；上述物理量均为无量纲量，长度采用的是诱导区长度进行无量纲化，时间采用诱导区长度与参考速度的比值进行无量纲化；
用C-J爆震波前导激波后气流温度对诱导反应和放热反应活化能进行归一化处理，即：
Figure GDA0003817505060000092
进一步的，S2包括以下步骤：
S21、取基元反应模型计算得到的混气C-J爆震波速VCJ和von Neumann状态下比热比
Figure GDA0003817505060000093
分别作为两步反应模型的C-J爆震波速VCJ和比热比γ；
S22、从基元反应模型计算得到的C-J状态参数(温度
Figure GDA0003817505060000094
压力
Figure GDA0003817505060000095
密度
Figure GDA0003817505060000096
)、vonNeumann状态参数(温度
Figure GDA0003817505060000097
压力
Figure GDA0003817505060000098
密度
Figure GDA0003817505060000099
)及C-J爆震波马赫数
Figure GDA00038175050600000910
共七个爆震波特征参数中任选其一，作为两步反应模型对应的特征参数，再结合比热比γ和vCJ，根据固定比热比爆震理论，唯一确定两步反应模型的无量纲放热量Q和气体常数R；
其中，若选取基元反应模型的C-J爆震波马赫数
Figure GDA00038175050600000911
作为两步反应模型的C-J爆震波马赫数MCJ，则两步反应模型无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure GDA00038175050600000912
Figure GDA00038175050600000913
其中
Figure GDA00038175050600000914
为爆震波波前混气静温。
若选取基元反应模型爆震波C-J状态温度
Figure GDA00038175050600000915
作为两步反应模型爆震波C-J状态温度TCJ，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure GDA00038175050600000916
Figure GDA0003817505060000101
Figure GDA0003817505060000102
其中
Figure GDA0003817505060000103
为爆震波波前混气静温。
若选取基元反应模型爆震波C-J状态压力
Figure GDA0003817505060000104
作为两步反应模型爆震波C-J状态压力PCJ，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure GDA0003817505060000105
Figure GDA0003817505060000106
Figure GDA0003817505060000107
其中
Figure GDA0003817505060000108
为爆震波波前混气静温，
Figure GDA0003817505060000109
为爆震波波前混气静压。
若选取基元反应模型爆震波C-J状态密度
Figure GDA00038175050600001010
作为两步反应模型爆震波C-J状态密度ρCJ，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure GDA00038175050600001011
Figure GDA00038175050600001012
Figure GDA00038175050600001013
其中
Figure GDA00038175050600001014
为爆震波波前混气静温，
Figure GDA00038175050600001015
为爆震波波前混气密度。
若选取基元反应模型爆震波von Neumann状态温度
Figure GDA00038175050600001016
作为两步反应模型爆震波von Neumann状态温度Tvn，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure GDA00038175050600001017
Figure GDA0003817505060000111
Figure GDA0003817505060000112
其中
Figure GDA0003817505060000113
为爆震波波前混气静温。
若选取基元反应模型爆震波von Neumann状态压力
Figure GDA0003817505060000114
作为两步反应模型爆震波von Neumann状态压力Pvn，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure GDA0003817505060000115
Figure GDA0003817505060000116
Figure GDA0003817505060000117
其中
Figure GDA0003817505060000118
为爆震波波前混气静温，
Figure GDA0003817505060000119
为爆震波波前混气静压。
若选取基元反应模型爆震波von Neumann状态密度
Figure GDA00038175050600001110
作为两步反应模型爆震波von Neumann状态密度ρvn，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure GDA00038175050600001111
Figure GDA00038175050600001112
Figure GDA00038175050600001113
其中
Figure GDA00038175050600001114
为爆震波波前混气静温，
Figure GDA00038175050600001115
为爆震波波前混气密度。
S23、根据确定的两步反应模型的无量纲放热量Q和比热比γ，根据固定比热比爆震理论计算两步反应模型对应的C-J状态参数(温度tCJ、压力pCJ、密度ρCJ)和von Neumann状态参数(温度tvn、压力pvn、密度ρvn、速度uvn)。
S24、假设定容爆炸的点火延迟时间τi具有Arrhenius形式：
Figure GDA00038175050600001116
其中T为温度，R为气体常数，EI为活化能，A为指前因子。
给定不同扰动爆震波速度±1％
Figure GDA0003817505060000121
采用基元反应计算经过爆震波前导激波压缩之后气体温度T1和T2，以及定容爆炸对应点火延迟时间τ1和τ2；根据定容爆炸理论计算无量纲诱导反应活化能εI：
Figure GDA0003817505060000122
S25、取诱导区化学反应速率常数kI＝-uvn，同时取定εR＝1.0；通过调整放热反应区化学反应速率常数kR，使得两步反应模型计算得到的爆震波一维稳定性参数χ与基元反应模型计算得到的爆震波一维稳定性参数
Figure GDA0003817505060000127
相等，从而确定放热反应区化学反应速率常数kR。
进一步的，计算爆震波一维稳定性参数χ方法如下：
Figure GDA0003817505060000123
其中εI为诱导反应活化能；ΔI为诱导区长度，定义为前导激波与反应热释放速率最大的位置之间的长度；ΔR放热区长度，定义为：
Figure GDA0003817505060000124
其中uCJ为爆震波坐标系下，CJ状态下当地声速；
Figure GDA0003817505060000125
为化学反应热释放速率，可通过下式计算：
Figure GDA0003817505060000126
其中c为当地声速。
具体实施例一：
下表给出了按照本发明建模方法计算的P0＝1atm,T0＝300K,化学恰当比乙烯/空气均匀预混气两步反应参数及两步反应模型与基元反应模型爆轰燃烧特征参数对比。同时图1给出了两步反应模型混气C-J爆轰波ZND结构温度与基元反应模型计算结果对比。其中，建模中保证两步反应模型与基元反应模型计算得到的爆轰波马赫数保持一致。可以看出，本发明所述建模方法得到的两步反应模型，能够很好描述爆轰波的ZND结构，两步反应温度曲线变化趋势与基元反应结果定性十分接近，说明了本发明的可行性和有效性。
Q γ ε<sub>I</sub> ε<sub>R</sub> k<sub>R</sub> R(J/(kg·K)) χ Δ<sub>I</sub>(m)
27.8 1.26 13.53 1.00 6.95 317.0 12.90 6.07E-4
P0＝1atm,T0＝300K,化学恰当比乙烯/空气两步反应模型参数
Figure GDA0003817505060000131
P0＝1atm,T0＝300K,化学恰当比乙烯/空气两步反应模型与基元反应模型爆轰燃烧特征对比
具体实施例二：
进一步，下表及图2给出了采用本发明建模方法对大分子碳氢燃料JP-10两步反应建模结果与基元反应模型对比。同样在建模中保证两步反应模型与基元反应模型计算得到的爆轰波马赫数保持一致。建模对应的状态及燃料组分参数为P0＝1atm,T0＝500K，化学恰当比JP-10/空气。大分子碳氢燃料燃烧过程十分复杂，典型的基元反应模型包含几十上百种组分，对应的反应在几百到几千量级，直接利用其进行数值模拟十分困难。而采用本文建模方法得到的两步反应模型，能够十分简单并且相对准确的描述其发生爆轰燃烧时对应的激波诱导及燃烧放热过程，充分说明了本建模方法的特点和优势。
Q γ ε<sub>I</sub> ε<sub>R</sub> k<sub>R</sub> R(J/(kg·K)) χ Δ<sub>I</sub>(m)
15.8 1.27 6.87 1.00 5.27 286.9 5.02 8.62E-4
P0＝1atm,T0＝500K,化学恰当比JP-10/空气两步反应模型参数
Figure GDA0003817505060000132
P0＝1atm,T0＝500K,化学恰当比JP-10/空气两步反应模型与基元反应模型爆轰燃烧特征对比
具体实施例三：
为进一步验证本发明两步反应模型计算方法对实际燃料爆震燃烧模拟情况，采用两步反应模型对对旋转爆震发动机燃烧室流场进行数值模拟，并将模拟得到的爆震波特征参数与基元反应模型结果进行对比。同样在建模中保证两步反应模型与基元反应模型计算得到的爆震波速度和马赫数保持一致。建模对应的状态及燃料组分参数为P0＝5500Pa,T0＝216K，化学恰当比氢气/空气预混气，其中上述状态参数对应旋转爆震发动机典型工作高度H＝20km来流静压和静温条件。下表给出了采用本发明计算得到的两步反应模型参数。
Q γ ε<sub>I</sub> ε<sub>R</sub> k<sub>R</sub> R(J/(kg·K)) χ Δ<sub>I</sub>(m)
24.9 1.32 5.77 1.00 2.63 422 1.59 2.49E-3
P0＝5500Pa,T0＝216K,化学恰当比氢气/空气预混气两步反应模型参数
旋转爆震发动机的燃烧室构型主要为等直的圆环，当燃烧室的径向厚度与内径相比很小时，可忽略其径向的影响。把燃烧室沿其母线展开得到二维矩形计算区域，对上下边界采用周期性边界条件以实现爆震波的周期性传播，左侧边界为预混气超声速入口边界条件，右侧出口采用外推边界条件。二维数值模拟算例计算条件设置为：周向长度Y＝240ΔI，轴向长度X＝160ΔI，反应物为化学恰当比的氢气/空气预混气，预混气填充总压Pt＝29P0，预混气总温Tt＝3.6T0。数值模拟收敛之后，旋转爆震流场温度和压力分布云图如图3所示。从图中可以看出典型的旋转爆震流场中爆震波-斜激波-滑移线典型结构特征，并且压力云图中可以明显看到爆震波面存在沿波面法向左右运动的三波点和横波结构。上述结果表明两步反应模型能够较为精细的模拟典型爆震波波系结构特征。
图4所示为监测到的不同时刻爆震波面传播速度的变化。为方便对比，采用基元反应模型计算的C-J爆震波速度对监测到的两步反应模型数值模拟爆震波传播速度进行归一化处理。可以看出，横波结构的存在导致爆震波面传播速度出现周期性的变化，爆震波传播速度整体在0.8到1.0左右上下波动，平均速度大约C-J速度的90％。爆震波面运动横波结构的存在导致爆震波传播速度发生波动，而爆震波后高温高压燃烧产物向燃烧室下游的侧向膨胀导致实际爆震波传播速度一定程度小于C-J速度，上述数值模拟速度与C-J速度之间的差异是正常的物理现象。上述结果说明本发明计算方法针对实际燃料建模结果，能够实现对燃料爆震燃烧特征参数的准确模拟，说明了本发明的可行性和有效性。
以上所述，仅为本发明的较佳的具体实施方式，但本发明的保护范围并不局限于此，所有熟悉本技术领域的技术人员在本发明公开的技术范围内，根据本发明的技术方案及其本发明的构思加以等同替换或改变均应涵盖在本发明的保护范围之内。

Claims (3)

1.一种针对燃料爆震燃烧的两步反应模型计算方法，其特征在于，包括：
S1、根据基元反应模型确定爆震波关键参数；
S2、根据所述基元反应模型爆震波关键参数，构建两步反应模型；
所述S2两步反应模型为诱导-放热两步反应模型；其使用两个反应速率控制方程分别模拟爆震燃烧的两个过程，第一步诱导反应代表的是诱导区内点火过程，第二步放热反应代表的是放热反应区内热量释放过程，其具体形式如下：
Figure FDA0003888976890000011
Figure FDA0003888976890000012
其中ξ是诱导反应进度变量，λ是放热反应进度变量，H(1-ξ)是一个阶跃函数：
Figure FDA0003888976890000013
kI是诱导区化学反应速率常数，取kI＝-uvn；其中uvn是爆震波坐标系下，C-J(Chapman–Jouguet)爆震波中前导激波后气流的速度；kR是放热反应区化学反应速率常数；EI和ER分别是诱导区和放热反应区的活化能，Ts是前导激波后气流的温度；上述物理量均为无量纲量，长度采用的是诱导区长度进行无量纲化，时间采用诱导区长度与参考速度的比值进行无量纲化；
用C-J爆震波前导激波后气流温度对诱导反应和放热反应活化能进行归一化处理，即：
Figure FDA0003888976890000014
所述S2包括以下步骤：
S21、取基元反应模型计算得到的混气C-J爆震波速
Figure FDA0003888976890000015
和von Neumann状态下比热比
Figure FDA0003888976890000021
分别作为两步反应模型的C-J爆震波速VCJ和比热比γ；
S22、从基元反应模型计算得到的C-J状态参数：温度
Figure FDA0003888976890000022
压力
Figure FDA0003888976890000023
密度
Figure FDA0003888976890000024
vonNeumann状态参数：温度
Figure FDA0003888976890000025
压力
Figure FDA0003888976890000026
密度
Figure FDA0003888976890000027
及C-J爆震波马赫数
Figure FDA0003888976890000028
共七个爆震波特征参数中任选其一，作为两步反应模型对应的特征参数，再结合比热比γ和vCJ，根据固定比热比爆震理论，唯一确定两步反应模型的无量纲放热量Q和气体常数R；
其中，若选取基元反应模型的C-J爆震波马赫数
Figure FDA0003888976890000029
作为两步反应模型的C-J爆震波马赫数MCJ，则两步反应模型无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure FDA00038889768900000210
Figure FDA00038889768900000211
其中
Figure FDA00038889768900000212
为爆震波波前混气静温；
若选取基元反应模型爆震波C-J状态温度
Figure FDA00038889768900000213
作为两步反应模型爆震波C-J状态温度TCJ，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure FDA00038889768900000214
Figure FDA00038889768900000215
Figure FDA00038889768900000216
其中
Figure FDA0003888976890000031
为爆震波波前混气静温；
若选取基元反应模型爆震波C-J状态压力
Figure FDA0003888976890000032
作为两步反应模型爆震波C-J状态压力PCJ，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure FDA0003888976890000033
Figure FDA0003888976890000034
Figure FDA0003888976890000035
其中
Figure FDA0003888976890000036
为爆震波波前混气静温，
Figure FDA0003888976890000037
为爆震波波前混气静压；
若选取基元反应模型爆震波C-J状态密度
Figure FDA0003888976890000038
作为两步反应模型爆震波C-J状态密度ρCJ，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure FDA0003888976890000039
Figure FDA00038889768900000310
Figure FDA00038889768900000311
其中
Figure FDA00038889768900000312
为爆震波波前混气静温，
Figure FDA00038889768900000313
为爆震波波前混气密度；
若选取基元反应模型爆震波von Neumann状态温度
Figure FDA00038889768900000314
作为两步反应模型爆震波vonNeumann状态温度Tvn，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure FDA0003888976890000041
Figure FDA0003888976890000042
Figure FDA0003888976890000043
其中
Figure FDA0003888976890000044
为爆震波波前混气静温；
若选取基元反应模型爆震波von Neumann状态压力
Figure FDA0003888976890000045
作为两步反应模型爆震波vonNeumann状态压力Pvn，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure FDA0003888976890000046
Figure FDA0003888976890000047
Figure FDA0003888976890000048
其中
Figure FDA0003888976890000049
为爆震波波前混气静温，
Figure FDA00038889768900000410
为爆震波波前混气静压；
若选取基元反应模型爆震波von Neumann状态密度
Figure FDA00038889768900000411
作为两步反应模型爆震波vonNeumann状态密度ρvn，则两步反应无量纲放热量Q和气体常数R可通过以下方程组求解得到：
Figure FDA00038889768900000412
Figure FDA00038889768900000413
Figure FDA0003888976890000051
其中
Figure FDA0003888976890000052
为爆震波波前混气静温，
Figure FDA0003888976890000053
为爆震波波前混气密度；
S23、根据确定的两步反应模型的无量纲放热量Q和比热比γ，根据固定比热比爆震理论计算两步反应模型对应的C-J状态参数：温度tCJ、压力pCJ、密度ρCJ；和von Neumann状态参数：温度tvn、压力pvn、密度ρvn、速度uvn；
S24、假设定容爆炸的点火延迟时间τi具有Arrhenius形式：
Figure FDA0003888976890000054
其中T为温度，R为气体常数，EI为活化能，A为指前因子；
给定不同扰动爆震波速度
Figure FDA0003888976890000055
采用基元反应计算经过爆震波前导激波压缩之后气体温度T1和T2，以及定容爆炸对应点火延迟时间τ1和τ2；根据定容爆炸理论计算无量纲诱导反应活化能εI：
Figure FDA0003888976890000056
S25、取诱导区化学反应速率常数kI＝-uvn，同时取定εR＝1.0；通过调整放热反应区化学反应速率常数kR，使得两步反应模型计算得到的爆震波一维稳定性参数χ与基元反应模型计算得到的爆震波一维稳定性参数
Figure FDA0003888976890000057
相等，从而确定放热反应区化学反应速率常数kR；
S3、采用实验及数值模拟验证所述两步反应模型。
2.根据权利要求1所述的一种针对燃料爆震燃烧的两步反应模型计算方法，其特征在于，所述S1包括以下步骤：
S11、首先采用基元反应模型计算给定初始状态燃料/氧化剂均匀混合气体中C-J爆震波速
Figure FDA0003888976890000061
C-J爆震波马赫数
Figure FDA0003888976890000062
C-J状态参数及von Neumann状态参数；其中C-J状态参数包括：温度
Figure FDA0003888976890000063
压力
Figure FDA0003888976890000064
及密度
Figure FDA0003888976890000065
von Neumann状态参数包括：温度
Figure FDA0003888976890000066
压力
Figure FDA0003888976890000067
密度
Figure FDA0003888976890000068
比热比
Figure FDA0003888976890000069
S12、利用基元反应模型计算混气C-J爆震波对应的
Figure FDA00038889768900000610
Figure FDA00038889768900000611
结构，确定诱导区实际物理长度ΔI，同时计算爆震波一维稳定性参数
Figure FDA00038889768900000612
3.根据权利要求1所述的一种针对燃料爆震燃烧的两步反应模型计算方法，其特征在于，所述计算爆震波一维稳定性参数χ方法如下：
Figure FDA00038889768900000613
其中εI为诱导反应活化能；ΔI为诱导区长度，定义为前导激波与反应热释放速率最大的位置之间的长度；ΔR放热区长度，定义为：
Figure FDA00038889768900000614
其中uCJ为爆震波坐标系下，CJ状态下当地声速；
Figure FDA00038889768900000615
为化学反应热释放速率，可通过下式计算：
Figure FDA00038889768900000616
其中c为当地声速。
CN202011439758.7A 2020-12-10 2020-12-10 一种针对燃料爆震燃烧的两步反应模型计算方法 Active CN112562793B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011439758.7A CN112562793B (zh) 2020-12-10 2020-12-10 一种针对燃料爆震燃烧的两步反应模型计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011439758.7A CN112562793B (zh) 2020-12-10 2020-12-10 一种针对燃料爆震燃烧的两步反应模型计算方法

Publications (2)

Publication Number Publication Date
CN112562793A CN112562793A (zh) 2021-03-26
CN112562793B true CN112562793B (zh) 2023-01-03

Family

ID=75060811

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011439758.7A Active CN112562793B (zh) 2020-12-10 2020-12-10 一种针对燃料爆震燃烧的两步反应模型计算方法

Country Status (1)

Country Link
CN (1) CN112562793B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113294264B (zh) * 2021-04-16 2022-08-19 中国人民解放军战略支援部队航天工程大学 基于针栓喷注器的双组元变推力旋转爆震火箭发动机
WO2023003685A1 (en) 2021-07-20 2023-01-26 Combustion Research And Flow Technology, Inc. Method for designing a combustion system with reduced environmentally-harmful emissions

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104870780A (zh) * 2012-10-12 2015-08-26 阿卜杜拉国王科技大学 爆震发动机
WO2018076403A1 (zh) * 2016-10-25 2018-05-03 浙江邦业科技股份有限公司 一种预测水泥回转窑熟料质量的一维仿真方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP4036138B2 (ja) * 2003-05-02 2008-01-23 日産自動車株式会社 火花点火式内燃機関の燃焼制御装置
CN110516310B (zh) * 2019-07-31 2022-10-04 中国空气动力研究与发展中心 旋转爆震反压的非定常数值模拟方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104870780A (zh) * 2012-10-12 2015-08-26 阿卜杜拉国王科技大学 爆震发动机
WO2018076403A1 (zh) * 2016-10-25 2018-05-03 浙江邦业科技股份有限公司 一种预测水泥回转窑熟料质量的一维仿真方法

Also Published As

Publication number Publication date
CN112562793A (zh) 2021-03-26

Similar Documents

Publication Publication Date Title
CN112562793B (zh) 一种针对燃料爆震燃烧的两步反应模型计算方法
Grasso et al. Three-dimensional computations of flows in a stratified-charge rotary engine
CN112417596B (zh) 一种航空发动机燃烧室通流模型并行网格仿真方法
CN110765698B (zh) 一种燃气轮机燃烧室变工况排放性能预测方法
Kyprianidis et al. EVA: A tool for environmental assessment of novel propulsion cycles
CN109408915B (zh) 固体火箭超燃冲压发动机燃烧流场仿真方法
Lietz et al. Parametric investigation of rotating detonation rocket engines using large eddy simulations
Lietz et al. Numerical investigation of rotating detonation rocket engines
Ristori et al. Numerical simulation of ducted rocket motor
Iannetti et al. The effect of spray initial conditions on heat release and emissions in LDI CFD calculations
CN110057536B (zh) 发动机燃烧条件下的吸气式飞行器内外流耦合模拟方法
Roy et al. A physics based reactor network model of a rotating detonation engine combustor
CN110442934A (zh) 一种考虑固体发动机尾喷流辐射的高精度气动热计算方法
Lopes et al. Are the available data from laboratory spray burners suitable for CFD modelling validations? A review
CN113987694A (zh) 一种基于空间推进算法的旋转爆震发动机流场参数分布预测方法
Zhu et al. Conceptual design and optimization of scramjet engines using the exergy method
Dubovas et al. Research on characteristics of turbo jet engine combustion chamber
CN116738894B (zh) 一种火箭发动机燃气流动数值仿真的方法
Kalina et al. Report on the implementation of the POIG project „turbine engine with a detonation combustion chamber”
Nyong et al. Prediction of Ignition Delay Behavior of Aviation Jet Fuel Model in a Constant Volume Adiabatic Reactor Relevant to Gas Turbine Conditions
Koshelev et al. Application of reactingCentralFOAM for modeling processes in combustion test chamber
Shunn et al. GPU-Accelerated High-Fidelity Rotating Detonation Engine Simulations Using an Extended Flamelet Progress Variable Approach
Liu et al. Generic Modeling Method of Quasi-One-Dimensional Flow for Aeropropulsion System Test Facility. Symmetry 2022, 14, 1161
Andrew et al. Numerical Investigation of Geometry Effect Channel on Fuel Source Flow Distribution in Micro Gas Turbine Annular Combustor
Nguyen et al. Aerodynamics and combustion of a realistic annular gas turbine combustor: A simulation study

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