CN114357363A - 一种计算零维氢气-空气燃烧反应的自适应步长方法 - Google Patents

一种计算零维氢气-空气燃烧反应的自适应步长方法 Download PDF

Info

Publication number
CN114357363A
CN114357363A CN202111335099.7A CN202111335099A CN114357363A CN 114357363 A CN114357363 A CN 114357363A CN 202111335099 A CN202111335099 A CN 202111335099A CN 114357363 A CN114357363 A CN 114357363A
Authority
CN
China
Prior art keywords
chemical reaction
ordinary differential
equation set
differential equation
time
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.)
Pending
Application number
CN202111335099.7A
Other languages
English (en)
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 Aerospace Testing Technology
Original Assignee
Beijing Institute of Aerospace Testing Technology
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 Aerospace Testing Technology filed Critical Beijing Institute of Aerospace Testing Technology
Priority to CN202111335099.7A priority Critical patent/CN114357363A/zh
Publication of CN114357363A publication Critical patent/CN114357363A/zh
Pending legal-status Critical Current

Links

Images

Landscapes

  • Feedback Control In General (AREA)

Abstract

本发明公开了一种计算零维氢气‑空气燃烧反应的自适应步长方法，能够随着反应的进行不断调整时间步长，节省计算时间。在经典的自适应步长策略的基础上，根据氢气‑氧气化学反应模型进行调整，随着反应的进行不断调整时间步长，节省计算时间。相比于龙格库塔法，本发明方法可以节省大量的计算量，易于实现自适应步长的计算。本发明方法基于Gear方法，即一种求解常微分方程组的隐式多步法，根据氢氧燃烧的特征进行了一些改进，使得氢氧燃烧的计算更加稳定、高效。

Description

一种计算零维氢气-空气燃烧反应的自适应步长方法
技术领域
本发明涉及计算化学技术领域，具体涉及一种计算零维氢气-空气燃烧反应的自适应步长方法。
背景技术
氢气-空气燃烧反应可分为三个阶段：诱导阶段、放热阶段与平衡阶段。反应初期称之为诱导阶段，此时各组分变化极小，但相互之间的反应速率可能相差很大。之后化学反应进入放热阶段，组分反应速率大幅加快；最后的阶段是平衡阶段，各组分缓慢趋于最终的极限值。
在各个反应阶段，不同化学组分的反应速率相差很大，因此可取得的最大时间步长也相差很大。反应初期的诱导阶段的数值稳定性是最差的，因此能取得的最大时间步长最小。在计算氢气-空气燃烧反应的化学反应常微分方程组时，随着反应的进行，最大时间步长总是变化，需要不断调整，延长了反应时间。
目前常用的化学反应常微分方程组求解方法是显式的龙格库塔法。龙格库塔法是一种单步法，具有计算简单，稳定性好等特点。但是多级的龙格库塔法误差估计较为复杂，不太适合进行自适应步长的计算。
为了使得求解更加高效，目前亟需一种自适应时间步长的方法，能够随着反应的进行不断调整时间步长，节省计算时间。
发明内容
有鉴于此，本发明提供了一种计算零维氢气-空气燃烧反应的自适应步长方法，能够自主调节时间步长，节省计算时间。
为实现上述发明目的，本发明的技术方案为：
一种计算零维氢气-空气燃烧反应的自适应步长方法，具体步骤包括：
步骤一、设定时间步长h初始值为h0、精度阶k的初始值为2、期望精度阶X、组分质量分数yn初始为y0和历史最大截断误差O初始值为-1；其中，n为步数。
步骤二、判断精度阶k是否等于期望精度阶X；是，则进入步骤三；否，则校正氢气-空气燃烧反应的化学反应常微分方程组Zn后，使精度阶k自增1；时间推进，更新Zn并使步数n自增1，重新执行步骤二。
步骤三、判断历史最大截断误差O是否小于0；若是，则进入步骤四；否则更新最大截断误差o，使历史最大截断误差O取值为o，重新执行步骤三。
步骤四、更新化学反应常微分方程组Zn，在当前的时间步数n的值上加1；更新o，判断o是否大于20O或小于0.5O，是则进入步骤五，否则回到步骤三。
步骤五、判断时间步长h是否大于预设的最大时间步长；若是，结束本流程并输出当前计算得到的组分质量分数yn；否则变更时间步长h，并更新Zn的第2到k+1行后，输出当前计算得到的组分质量分数yn，回到步骤四。
进一步的，化学反应常微分方程组Zn为：
Figure BDA0003350282390000021
其中，y'n为组分质量分数yn的一阶导数，是化学反应常微分方程组Zn的第1行，以此类推，yn (k)为yn的k阶导数，yn (k)为化学反应常微分方程组Zn的第k+1行，Δt为时间间隔，每一行包含ns列；其中ns为组分个数。
进一步的，校正迭代化学反应常微分方程组Zn，使精度阶k自增1，具体方法为：
将化学反应常微分方程组Zn进行泰勒展开后，将Zn中的元素更新，得到泰勒展开估计的方程组化学反应常微分方程组
Figure BDA0003350282390000022
化学反应常微分方程组Zn中的每个元素泰勒展开，其中泰勒展开后第k+1行元素第j列元素更新为：
Figure BDA0003350282390000031
其中，
Figure BDA0003350282390000032
为在第n步的化学反应常微分方程组Zn的第k+1行，第j列元素；
Figure BDA0003350282390000033
为在第n-1步的化学反应常微分方程组Zn的第k+1行，第j列元素；
Figure BDA0003350282390000034
为在第n步下的化学反应常微分方程组Zn的第k-1行，第j列元素；hk-1为时间步长h的k-1次方；j为化学反应常微分方程组Zn的列数且1≤j≤ns。
进一步的，更新最大截断误差o，具体采用如下方法为：
Figure BDA0003350282390000035
进一步的，更新化学反应常微分方程组Zn的具体方法为：
Figure BDA0003350282390000036
其中，Zn为n+k的时间步长下的化学反应常微分方程组，
Figure BDA0003350282390000037
为采用泰勒展开估计的化学反应常微分方程组，L为化学反应常微分方程组Zn的第k+1行的向量；w为化学反应常微分方程组Zn的第ns列的向量。
进一步的，变更时间步长h的求解公式为：
Figure BDA0003350282390000038
其中，hn+1为时间步长h的n+1次方，hn为时间步长h的n次方。
有益效果：
1、本发明提供了一种计算零维氢气-空气燃烧反应的自适应步长方法，运用Gear方法(一种求解常微分方程组的隐式多步法)，在经典的自适应步长策略的基础上，根据氢气-氧气化学反应模型进行调整，随着反应的进行不断调整时间步长，节省计算时间。
2、相比于龙格库塔法，本发明方法可以节省大量的计算量，易于实现自适应步长的计算。本发明方法基于Gear方法，根据氢氧燃烧的特征进行了一些改进，使得氢氧燃烧的计算更加稳定、高效。
附图说明
图1为本发明的方法流程图。
图2为时间步长增长曲线图。
图3为四种方法计算得到的氢气摩尔分数随时间变化的曲线与实验值的对比图。
具体实施方式
下面结合附图并举实施例，对本发明进行详细描述。
如图1所示，本发明提供了一种计算零维氢气-空气燃烧反应的自适应步长方法，运用Gear方法(一种求解常微分方程组的隐式多步法)，在经典的自适应步长策略的基础上，根据氢气-氧气化学反应模型进行调整。
具体步骤包括：
步骤一、设定时间步长h初始值为h0、精度阶k的初始值为2、期望精度阶X、组分质量分数yn初始为y0和历史最大截断误差O初始值为-1；其中，n为步数。
步骤二、判断精度阶k是否等于期望精度阶X；是，则进入步骤三；否，则校正氢气-空气燃烧反应的化学反应常微分方程组Zn后，使精度阶k自增1；时间推进，更新Zn并使步数n自增1，重新执行步骤二。
化学反应常微分方程组Zn为：
Figure BDA0003350282390000041
其中，y'n为组分质量分数yn的一阶导数，是化学反应常微分方程组Zn的第1行，以此类推，yn (k)为yn的k阶导数，yn (k)为化学反应常微分方程组Zn的第k+1行，Δt为时间间隔，每一行包含ns列；其中ns为组分个数。
校正迭代常微分方程组z的具体方法为：将常微分方程组zn进行泰勒展开后，将zn中的元素更新：
将化学反应常微分方程组Zn进行泰勒展开后，将Zn中的元素更新，得到泰勒展开估计的方程组化学反应常微分方程组
Figure BDA0003350282390000051
化学反应常微分方程组Zn中的每个元素泰勒展开，其中泰勒展开后第k+1行元素第j列元素更新为：
Figure BDA0003350282390000052
其中，
Figure BDA0003350282390000053
为在第n步的化学反应常微分方程组Zn的第k+1行，第j列元素；
Figure BDA0003350282390000054
为在第n-1步的化学反应常微分方程组Zn的第k+1行，第j列元素；
Figure BDA0003350282390000055
为在第n步下的化学反应常微分方程组Zn的第k-1行，第j列元素；hk-1为时间步长h的k-1次方；j为化学反应常微分方程组Zn的列数且1≤j≤ns。
步骤三、判断历史最大截断误差O是否小于0；若是，则进入步骤四；否则更新最大截断误差o，使历史最大截断误差O取值为o，重新执行步骤三。
更新最大截断误差o的具体方法为：
Figure BDA0003350282390000056
步骤四、更新化学反应常微分方程组Zn，在当前的时间步数n的值上加1；更新o，判断o是否大于20O或小于0.5O，是则进入步骤五，否则回到步骤三。
更新化学反应常微分方程组Zn的具体方法为：
Figure BDA0003350282390000057
其中，Zn为n+k的时间步长下的化学反应常微分方程组，
Figure BDA0003350282390000061
为采用泰勒展开估计的化学反应常微分方程组，L为化学反应常微分方程组Zn的第k+1行的向量；w为化学反应常微分方程组Zn的第ns列的向量。
步骤五、判断时间步长h是否大于预设的最大时间步长；若是，结束本流程并输出当前计算得到的组分质量分数yn；否则变更时间步长h，并更新Zn的第2到k+1行后，输出当前计算得到的组分质量分数yn，回到步骤四。
变更时间步长h的求解公式为：
Figure BDA0003350282390000062
其中，hn+1为时间步长h的n+1次方，hn为时间步长h的n次方。
本发明实施例中，使用的氢气-氧气化学反应模型的原理如下。对于零维的氢气-氧气燃烧反应，随时间变化的组成成分(以下简称“组分”)可用一个常微分方程组描述：
Figure BDA0003350282390000063
Figure BDA0003350282390000064
其中，ys为每个组分的质量分数，
Figure BDA0003350282390000065
为每个组分的反应速率，E为总能量，ρ为反应时所有组分的混合密度。
本发明实施例中，反应速率
Figure BDA0003350282390000066
的计算方法：
设整个化学反应系统由ns种组分组成，其化学动力学过程由nr个基元反应构成，则该化学反应系统的化学反应方程式可表示为：
Figure BDA0003350282390000067
其中，r＝1,2,...,nr。
其中，As为参加化学反应的组分s，v'rs为第r个化学反应方程式中反应物的化学计量系数，v″rs为第r个化学反应方程式中生成物的化学计量系数。
Figure BDA0003350282390000071
为第r个化学反应方程式的正向反应速率常数，
Figure BDA0003350282390000072
为第r个化学反应方程式的逆向反应速率常数。化学反应速率常数k可由分子碰撞理论导出的Arrhenius公式确定：
Figure BDA0003350282390000073
其中，A、b和C为实验拟合系数。需要特别注意实验拟合系数的量纲及单位，对于二体反应及三体反应，其量纲不同。
化学反应方程式以及相应的A、b和C参数设置如表1所示：
表1实验拟合系数表
Figure BDA0003350282390000074
化学反应源项
Figure BDA0003350282390000075
表示组分s在单位时间单位体积内的质量生成率(即反应物的消耗率和生成物的生成率)。根据质量作用定律，对于上述nr个基元反应，组分s的化学反应源项可表达为：
Figure BDA0003350282390000081
其中，Rr和R-r分别表第r个基元反应的正反应速率和逆反应速率，其具体形式为：
Figure BDA0003350282390000082
Figure BDA0003350282390000083
其中，
Figure BDA0003350282390000084
为三体碰撞项，表征三体碰撞对基元反应的影响，Crk为第k种组分在第r个反应中作为三体碰撞的碰撞系数。Lr为是否发生三体碰撞反应的控制参数。当出现三体碰撞时Lr设为1，三体碰撞项起作用；否则Lr为0，此时三体碰撞项的值为1。
基于上述模型，本发明采用的常微分方程组为：
Figure BDA0003350282390000085
其中，yn为待求向量在第n个时间步长的值，y'n为yn的一阶导数，以此类推，yn (k)为yn的k阶导数，n为时间步长数，Δt为时间间隔。
因为Gear方法最多达到6阶，Passcal矩阵A为7×7的矩阵。
由泰勒展开表达式有：
Zn+1＝AZn
泰勒展开表达式未用到原微分方程，因此在数值分析上它是不稳定的，还需要增加一个校正迭代，添加后的泰勒展开表达式如下：
Figure BDA0003350282390000086
L和F的值或表达式由Gear方法以及隐式多步法中对第n+1步初值的预估方法所确定。这里直接给出L和F的表达式：
L为k+1行1列的向量，其值见表2。
表2L的取值表
Figure BDA0003350282390000091
F为1行ns列的向量，其表达式为：
Figure BDA0003350282390000092
综上，算法可以写成：
Figure BDA0003350282390000093
Figure BDA0003350282390000094
这个迭代格式可以采用简单迭代法，也可以采用牛顿迭代法。本文中采取的是牛顿迭代法。
如果上式收敛，那么
Figure BDA0003350282390000095
将收敛到：
Figure BDA0003350282390000096
其中w为1行ns列的向量，并且满足：
Figure BDA0003350282390000097
上式是一个关于w非线性方程组，可以运用牛顿迭代法进行求解：
Figure BDA0003350282390000101
Figure BDA0003350282390000102
为避免矩阵求逆，将迭代公式化为求解线性方程组：
Figure BDA0003350282390000103
基于上述模型和常微分方程组，本发明使用的自适应步长策略为：
化学反应起始阶段，即诱导阶段初期，方程组的稳定性限制是最大的，可以取得的步长很小，因此在初期很难实现步长的增长。当反应完成了诱导阶段初期，稳定区域会大幅增长，此时需要增大时间步长以提升计算效率。当经过了诱导阶段初期之后，实时监测截断误差估计值，记历史最大截断误差值为o。当某一步的最大截断误差O远远大于o时，则变更步长(此时会使得步长变小)；当最大截断误差O远小于o时，变更步长(此时会使得步长稍微变大)；最大截断误差O与历史最大截断误差值o相差不大时，不改变时间步长。当修正后的时间步长大于可接受的最大时间步长时，则取消修正。
判断完成诱导阶段初期的标准：
||yn-y0||∞>10-5
变更时间步长判断标准：
Figure BDA0003350282390000104
或
Figure BDA0003350282390000105
其中，on为第n时间步长的最大截断误差，O为历史最大截断误差值，变更步长的方法为：
Figure BDA0003350282390000106
本发明实施例中，如图2给出的时间步长增长曲线。可以看出，经过诱导阶段后，时间步长有大幅度增长。
如图3所示，图3给出二级龙格库塔法、定步长四阶Gear方法、变步长四阶Gear方法计算得到的氢气摩尔分数随时间变化的曲线与实验值的对比，可以看到，变步长方法在提升计算效率的同时，维持了原有的计算精度。
如表3所示，表3给出了二级龙格库塔法、定步长四阶Gear方法与变步长四阶Gear方法的CPU耗时，可以看到自适应步长算法可以大幅缩减计算耗时。综上，可以认为本专利的自适应算法是有效的。
表3四种方法的CPU耗时时间表
数值方法 初始步长 CPU耗时
RK2 1d-5 16.81
BDF4(定步长) 1d-5 27.78
BDF4(变步长) 1d-5 7.62
综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种计算零维氢气-空气燃烧反应的自适应步长方法，其特征在于，具体步骤包括：
步骤一、设定时间步长h初始值为h0、精度阶k的初始值为2、期望精度阶X、组分质量分数yn初始为y0和历史最大截断误差O初始值为-1；其中，n为步数；
步骤二、判断所述精度阶k是否等于期望精度阶X；是，则进入步骤三；否，则校正氢气-空气燃烧反应的化学反应常微分方程组Zn后，使精度阶k自增1；时间推进，更新Zn并使步数n自增1，重新执行步骤二；
步骤三、判断历史最大截断误差O是否小于0；若是，则进入步骤四；否则更新最大截断误差o，使历史最大截断误差O取值为o，重新执行步骤三；
步骤四、更新化学反应常微分方程组Zn，在当前的时间步数n的值上加1；更新o，判断o是否大于20O或小于0.5O，是则进入步骤五，否则回到步骤三；
步骤五、判断时间步长h是否大于预设的最大时间步长；若是，结束本流程并输出当前计算得到的组分质量分数yn；否则变更时间步长h，并更新Zn的第2到k+1行后，输出当前计算得到的组分质量分数yn，回到步骤四。
2.如权利要求1所述的方法，其特征在于，所述化学反应常微分方程组Zn为：
Figure FDA0003350282380000011
其中，y′n为组分质量分数yn的一阶导数，是化学反应常微分方程组Zn的第1行，以此类推，yn (k)为yn的k阶导数，yn (k)为化学反应常微分方程组Zn的第k+1行，Δt为时间间隔，每一行包含ns列；其中ns为组分个数。
3.如权利要求2所述的方法，其特征在于，所述校正迭代化学反应常微分方程组Zn，使精度阶k自增1，具体方法为：
将化学反应常微分方程组Zn进行泰勒展开后，将Zn中的元素更新，得到泰勒展开估计的方程组化学反应常微分方程组
Figure FDA0003350282380000021
化学反应常微分方程组Zn中的每个元素泰勒展开，其中泰勒展开后第k+1行元素第j列元素更新为：
Figure FDA0003350282380000022
其中，
Figure FDA0003350282380000023
为在第n步的化学反应常微分方程组Zn的第k+1行，第j列元素；
Figure FDA0003350282380000024
为在第n-1步的化学反应常微分方程组Zn的第k+1行，第j列元素；
Figure FDA0003350282380000025
为在第n步下的化学反应常微分方程组Zn的第k-1行，第j列元素；hk-1为时间步长h的k-1次方；j为化学反应常微分方程组Zn的列数且1≤j≤ns。
4.如权利要求3所述的方法，其特征在于，所述更新最大截断误差o，具体采用如下方法为：
Figure FDA0003350282380000026
5.如权利要求3所述的方法，其特征在于，所述更新化学反应常微分方程组Zn的具体方法为：
Figure FDA0003350282380000027
其中，所述Zn为n+k的时间步长下的化学反应常微分方程组，
Figure FDA0003350282380000028
为采用泰勒展开估计的化学反应常微分方程组，L为化学反应常微分方程组Zn的第k+1行的向量；w为化学反应常微分方程组Zn的第ns列的向量。
6.如权利要求1所述的方法，其特征在于，所述变更时间步长h的求解公式为：
Figure FDA0003350282380000029
其中，hn+1为时间步长h的n+1次方，hn为时间步长h的n次方。
CN202111335099.7A 2021-11-11 2021-11-11 一种计算零维氢气-空气燃烧反应的自适应步长方法 Pending CN114357363A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111335099.7A CN114357363A (zh) 2021-11-11 2021-11-11 一种计算零维氢气-空气燃烧反应的自适应步长方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111335099.7A CN114357363A (zh) 2021-11-11 2021-11-11 一种计算零维氢气-空气燃烧反应的自适应步长方法

Publications (1)

Publication Number Publication Date
CN114357363A true CN114357363A (zh) 2022-04-15

Family

ID=81095373

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111335099.7A Pending CN114357363A (zh) 2021-11-11 2021-11-11 一种计算零维氢气-空气燃烧反应的自适应步长方法

Country Status (1)

Country Link
CN (1) CN114357363A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115618171A (zh) * 2022-06-06 2023-01-17 北京理工大学 一种基于同伦算法的推进剂燃烧平衡产物求解方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115618171A (zh) * 2022-06-06 2023-01-17 北京理工大学 一种基于同伦算法的推进剂燃烧平衡产物求解方法
CN115618171B (zh) * 2022-06-06 2023-10-24 北京理工大学 一种基于同伦算法的推进剂燃烧平衡产物求解方法

Similar Documents

Publication Publication Date Title
Gupta et al. On a statistic which arises in selection and ranking problems
CN114357363A (zh) 一种计算零维氢气-空气燃烧反应的自适应步长方法
Li et al. Robust H∞ control of TS fuzzy systems with input time-varying delays: a delay partitioning method
CN104036331A (zh) 一种基于改进粒子群算法的电力系统动态经济调度方法
Kaczorek Positive stable realizations with system Metzler matrices
CN110898454A (zh) 精馏塔控制方法、存储介质及电子设备
Chen et al. Stability analysis of systems with timevarying delay via a novel Lyapunov functional
Soloklo et al. Chebyshev rational functions approximation for model order reduction using harmony search
Haus et al. Long-time tail effects on particle diffusion in a disordered system
CN117608327A (zh) 一种用于反应釜温度控制系统的控制器优化方法及装置
Lin et al. Event-triggered stabilisation of sampled-data singular systems: a hybrid control approach
CN112507627A (zh) 基于混沌黏菌算法的多元非线性函数求解方法
CN112596557A (zh) 半导体温控装置输出量的控制方法及装置
CN111399376A (zh) 一种t-s模糊系统的二维重复控制器设计优化方法
CN108333947B (zh) 基于智能寻优的单调整系数预测函数控制参数整定方法
Zhou et al. A perturbed risk model with dependence between premium rates and claim sizes
CN113885325A (zh) 基于阶跃响应的一阶惯性加纯延迟环节闭环辨识方法
CN103401650A (zh) 一种(n,1,m)有误码卷积码的盲识别方法
Rabah et al. State feedback with fractional PIλDμ control structure for genesio-tesi chaos stabilization
CN115659792A (zh) 多时段多场景scuc解耦方法、系统、设备及存储介质
CN106950834A (zh) 一种高阶非线性系统快速有限时间稳定的控制方法
CN107968444B (zh) 一种新能源集群协调优化控制方法
Foscoliano et al. Improving the wastewater treatment plant performance through model predictive control strategies
CN108021999B (zh) 一种快速逼近最大负荷功率点的方法及装置
Shamgah et al. Decomposition via QP

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