一种燃烧模拟降维提速方法及装置、稳态计算方法与流程

1.本发明属于航空发动机燃烧室流动燃烧仿真模拟技术领域，具体涉及一种基于火焰面生成流形(fgm)模型的燃烧模拟降维提速方法。


背景技术：

2.燃烧室作为航空发动机的核心部件之一，其设计直接影响到发动机的整机性能。在航空发动机的研发中，燃烧室是一个包含多组分的燃烧化学反应系统，所发生的湍流燃烧过程极其复杂，可能包含上千步中间反应，直接求解详细的化学反应机理以及完整的组分输运方程所需要的计算量通常难以承受，对于湍流燃烧模型的研究带来挑战。
3.目前湍流燃烧模型工程常用的主要有：有限速率类模型和火焰面类模型。有限速率类模型对多组分输运方程进行求解，计算量随组分数量增多而增加；火焰面类模型引入少量控制变量，建立这些控制变量与化学热力学状态之间的对应关系，从而获得火焰面数据库，通过对这些少量控制变量的输运方程求解，再对火焰面数据库查表，即可得到全部化学热力学状态参数，这种方法由于求解输运方程的数量减小，从而减小了燃烧模拟的计算量。火焰面生成流形(fgm)模型是火焰面类模型的一种。在fgm模型中，有三方程模型、四方程模型等不同模式。基于平均混合分数混合分数方差平均进度变量三变量输运方程组的模式，只考虑了的均值，没有考虑进度变量的波动影响，精度相对较低。基于四变量输运方程组的模式，同时考虑了的均值和波动，计算结果更准确，但由于多求解一个输运方程，计算速度明显降低。


技术实现要素：

4.针对现有技术中存在不足，本发明提供了一种基于火焰面生成流形(fgm)模型的燃烧模拟降维提速方法，在精度误差允许范围内减少输运方程数量，实现燃烧计算过程的提速，具有较大的工程应用价值。
5.本发明是通过以下技术手段实现上述技术目的的。
6.一种燃烧模拟降维提速方法，用于火焰面数据库类燃烧模型，所述火焰面数据库类燃烧模型的变量至少包括平均混合分数混合分数方差平均进度变量和进度变量方差的四个变量；所述方法包括：对所述火焰面数据库类燃烧模型在燃烧室流体域的每个网格单元上按照时间步依次求解计算；其中，对于任一时间步的求解计算包括：
7.分别求解每个网格单元的各变量的值；
8.基于火焰面数据库，获取每个变量的列向量；
9.对于任一变量，根据每个网格单元的所述变量的值和所述变量的列向量，确定每个网格单元中所述变量在火焰面数据库中相应变量维度方向上的位置信息和权重信息；
10.根据每个变量对应的位置信息和权重信息，基于火焰面数据库确定当前时间步对应的化学热力学相关状态参数，并更新当前时间步的流体密度。
11.进一步地，当所述变量为平均混合分数混合分数方差或平均进度变量时，求解每个网格单元的任一变量的值，包括：基于所述火焰面数据库，求解所述变量的输运方程，获得每个网格单元的所述变量的值。
12.更进一步地，当所述变量为混合分数方差时，每个网格单元的的列向量为
13.其中，值为当地网格单元的混合分数的值；是从火焰面数据库中读取的一个列向量，所述列向量中元素按照从小到大排列，且每个元素的值均位于[0，1]范围内。
[0014]
更进一步地，当所述变量为进度变量方差时，基于火焰面数据库，获取进度变量方差的列向量，包括：
[0015]
根据网格单元的值和列向量fc(：,n)，获得位置信息(i
temp
，i
temp
1)和权重信息(pi，p
i 1
)，结合二维矩阵构造的列向量中第n个元素为的列向量中第n个元素为循环遍历维度上的每一个位置，获取的列向量；
[0016]
其中：fc(：,n)为二维矩阵fc(m,n)在维度上取第n个数、在维度上分别取从第1到第m个数时形成的列向量。
[0017]
更进一步地，求解每个网格单元的进度变量方差的值，包括：
[0018]
构造求解进度变量方差的代数式：其中：θ为预估值的比例系数，θ∈[0，1]；为当地网格中的列向量最大值与最小值和的二分之一；为进度变量方差预估值，且其中：t表示第t时间步，t 1表示第t 1时间步，μ
t
为湍流黏度，σ
t
为施密特数，为流体密度，k为湍动能，ε为湍动能耗散率，rf为系数，δt为时间步长，为梯度运算符。
[0019]
更进一步地，当所述变量为进度变量方差时，基于火焰面数据库，获取平均进度变量的列向量，包括：
[0020]
根据进度变量方差的值及其列向量，获取的位置信息(n，n 1)和权重信息(fn，f
n 1
)，结合二维矩阵fc(m,n)，构造的列向量中第m个元素为fc(m,n)*fn fc(m,n 1)*f
n 1
，循环遍历维度上每一个位置，获取的列向量；
[0021]
其中：fc(m,n)表示在维度上取第m个数、在维度上取第n个数时，对应值的大小。
[0022]
更进一步地，对于任一变量，根据每个网格单元的所述变量的值和所述变量的列向量，确定每个网格单元中所述变量在火焰面数据库中相应变量维度方向上的位置信息和权重信息，包括：
[0023]
根据每个网格单元的所述变量的值，利用二分法查询并获得所述变量在火焰面数据库中相应变量维度方向上的位置信息和权重信息；
[0024]
其中，所述位置信息是使所述变量的值x位于该变量对应的列向量中第i个元素和
第i 1个元素之间；所述权重信息是存在一对实数(α，β)，使所述变量的值x＝α*xi β*x
i 1
，且α β＝1，α、β均位于[0，1]范围内。
[0025]
进一步地，根据每个变量对应的位置信息和权重信息，基于火焰面数据库确定当前时间步对应的化学热力学相关状态参数，包括：温度、组分质量分数、比热容、导热率、密度、混合气体的分子粘度。
[0026]
一种燃烧模拟降维提速装置，包括：
[0027]
获取进度变量方差模块，用于构造求解进度变量方差的代数式，并得到进度变量方差；
[0028]
位置权重信息获取模块，用于得到平均混合分数混合分数方差平均进度变量以及进度变量方差的位置信息和权重信息；
[0029]
更新模块，基于获取的位置信息和权重信息，更新化学热力学物性场，并更新流体域密度场。
[0030]
一种稳态计算方法，所述方法包括：
[0031]
采用上述燃烧模拟降维提速方法，进行基于伪时间方法的稳态计算，计算结果稳定收敛后,采用四方程模型方法进行基于伪时间的稳态计算。
[0032]
本发明的有益效果为：
[0033]
(1)本发明通过构造进度变量方差的代数式，对采取代数求解，将四方程fgm模型中关于平均混合分数混合分数方差平均进度变量进度变量方差四个变量的输运方程组降为关于三个变量的输运方程组，在精度误差允许范围内实现湍流燃烧过程的更快仿真模拟计算，同时保证了计算结果的精确度和效率，具有较高的工程适用性，为发动机燃烧室内的湍流燃烧过程的快速模拟提供新的方法。
[0034]
(2)本发明根据进度变量方差列向量，求取列向量的最大值和最小值列向量的最大值和最小值对每个网格的进度变量方差预估值进行截断，确保其在合理的预估范围内。
附图说明
[0035]
图1为本发明所述燃烧模拟降维提速方法流程图；
[0036]
图2为本发明所述s1的具体实施流程图；
[0037]
图3为本发明所述s2的具体实施流程图；
[0038]
图4为本发明所述火焰面数据库查表流程图；
[0039]
图5为本发明所述s3的具体实施流程图；
[0040]
图6为本发明所述s4的具体实施流程图；
[0041]
图7(a)为气相燃烧四输运方程温度分布图；
[0042]
图7(b)为气相燃烧本发明方法的温度分布图；
[0043]
图7(c)为气相燃烧四输运方程y_co2分布图；
[0044]
图7(d)为气相燃烧本发明方法的y_co2分布图；
[0045]
图7(e)为气相燃烧四输运方程y_h2o分布图；
[0046]
图7(f)为气相燃烧本发明方法的y_h2o分布图；
[0047]
图7(g)为气相燃烧四输运方程分布图；
[0048]
图7(h)为气相燃烧本发明方法的分布图；
[0049]
图8(a)为气液两相燃烧四输运方程温度分布图；
[0050]
图8(b)为气液两相燃烧本发明方法的温度分布图；
[0051]
图8(c)为气液两相燃烧四输运方程y_co2分布图；
[0052]
图8(d)为气液两相燃烧本发明方法的y_co2分布图；
[0053]
图8(e)为气液两相燃烧四输运方程y_h2o分布图；
[0054]
图8(f)为气液两相燃烧本发明方法的y_h2o分布图；
[0055]
图8(g)为气液两相燃烧四输运方程分布图；
[0056]
图8(h)为气液两相燃烧本发明方法的分布图。
具体实施方式
[0057]
下面结合附图以及具体实施例对本发明作进一步的说明，但本发明的保护范围并不限于此。
[0058]
如图1所示，本发明一种基于火焰面生成流形(fgm)模型的燃烧模拟降维提速方法，具体包括如下步骤：
[0059]
s1，四变量火焰面数据库构建，的输运方程变量求解，以及值的预估计。
[0060]
在进行流场计算之前，需要先完成火焰面数据库的构建和输运方程变量的求解。在此，为了降维提速，只求解三个变量的输运方程，对于进度变量方差通过构造基于脉动生成和耗散效应的代数表达式进行预估计。具体参见图2：
[0061]
s11，构建四变量的湍流火焰面数据库。
[0062]
本发明中具体的构建火焰面数据库的方法为现有技术，获取四变量在湍流火焰面数据库中的位置信息和权重信息，从而确定对应的状态参数；特别指出的，在湍流火焰面数据库状态参数表的最后两列增加列出了进度变量方差和平均进度变量的值，以备后续构建这两个量的列向量时使用。其中状态参数表中包含温度、组分质量分数、比热容、导热率、密度、混合气体的分子粘度。
[0063]
s12，求解输运方程，对火焰面模型在燃烧室流体域中每个网格单元上进行求解计算，获得每个网格单元的数值。
[0064]
对三个变量的输运方程进行数值求解，方程形式为现有技术。特别指出的，在计算平均进度变量输运方程的过程中，将计算出的的梯度场数据在内存中进行保存，以备后续构造的代数表达式时使用。
[0065]
s13，遍历每个流体域网格单元，构造基于脉动生成和耗散效应的进度变量方差预估值的代数表达式。
[0066]
考虑到进度变量方差的产生主要是由于湍流脉动造成的标量方差的生成和耗散
效应，在此首先构造进度变量方差预估值后续再对预估值进行加权修正，得到最终的进度变量方差的代数表达式。
[0067]
进度变量方差的生成率为耗散率为(生成率和耗散率的表达式为现有技术)，其中μ
t
为湍流黏度，σ
t
为施密特数，通常取0.7，ρ为流体密度，k为湍动能，ε为湍动能耗散率，系数rf通常取为0.5，为梯度运算符。
[0068]
忽略对流和扩散效应，构造表达式其中δt为时间步长。
[0069]
忽略气体密度在极小时间步长内变化的影响，得到进度变量方差预估值忽略气体密度在极小时间步长内变化的影响，得到进度变量方差预估值其中t表示第t时间步，t 1表示第t 1时间步。在计算时已完成t 1步的速度、密度、压力的计算，以及湍流模型计算，因此，等式右边均为已知量，对于初始场其进度变量方差设为0。通过预估值表达式可以看出，当值过大时，值会由于耗散率的增加而减小，不会在值过大的情况下继续增大；当值过小时，导致耗散率接近0，值会由于生成率的原因而增大，不会在值过小的情况下继续减小，由此可以反映出表达式构造的合理性。
[0070]
s2，获得和的列向量，查询和的位置信息与权重信息。
[0071]
针对每个网格单元的和值，在获得和的列向量的基础上，采用二分法确定和值在其列向量中所处的位置和权重信息。具体参见图3：
[0072]
s21，获得列向量并查询的位置信息和权重信息。
[0073]
为方便描述，这里把由k个控制变量构建的火焰面数据库，称为k维火焰面数据库。控制变量的列向量用于查询该控制变量在火焰面数据库中相应维度方向上的位置信息和权重信息。列向量中每个元素按照从小到大的顺序排列。遍历流体域的每个网格单元之前，已通过输运方程求解获得网格单元的值，通过读取火焰面数据库，得到的列向量(每个网格上的的列向量均相同)。在遍历流体域的每个网格单元时，根据每个网格计算的值，利用二分法查询其在列向量中的位置信息和权重信息。其中，值的位置信息指一对整数值(i，i 1)，使得的数值大小在列向量中位于第i个元素和第i 1个元素之间，使得的数值大小在列向量中位于第i个元素和第i 1个元素之间，使得值的权重信息指一对的实数值(ai，a
i 1
)，使得且ai a
i 1
＝1，0《＝ai，a
i 1
《＝1，如图4所示。
[0074]
s22，构造列向量并查询的位置信息和权重信息。
[0075]
首先，在遍历每个网格单元之前，从火焰面数据库中读取一个列向量该向量每个元素数值都在0到1之间，且各不相同，从小到大排列，第一元素值为0，最后一个元素值为1。在遍历流体域每个网格单元时，构造该网格单元的的列向量为其中值
为该网格单元的混合分数的数值；而后通过二分法查找值在其列向量中的位置信息(j，j 1)和权重信息(bj，b
j 1
)，具体方法与值查找方法相同。
[0076]
s3，构造的列向量、值的代数表达式。
[0077]
构造的列向量，先后查询和值的位置信息和权重信息。
[0078]
其中，值的代数表达式的构造为其中为当地网格中的列向量最大值与最小值和的二分之一，为基于脉动生成与耗散效应构造的预估值，θ为预估值的比例系数。通过调整θ值大小，可以用于标定不同算例的参数设置，θ值越大，考虑生成和耗散效应的预估值权重越大，θ值越小越，接近于值可能产生波动大小的中间水平。具体参见图5：
[0079]
s31，在遍历流体域网格时，根据和值查询信息，通过火焰面数据库f(i,j,p,q)，分别构建和的二维矩阵fc(m,n)和其中，i表示二维矩阵中维度上的位置序号，j表示二维矩阵中维度上的位置序号，m表示二维矩阵中维度上的位置序号，n表示二维矩阵中维度上的位置序号；火焰面数据库函数f(i,j,m,n)表示在维度上取第i个数、在维度上取第j个数、在维度上取第m个数、在维度上取第n个数时，对应某个化学热力学状态量(比如温度f
t
、组分质量分数fy、比热容f
cp
、导热率f
λ
或在s11中指出最后两列列出的值fc、值)的大小；fc(m,n)表示在维度上取第m个数、在维度上取第n个数时，对应值的大小；表示在维度上取第m个数、在维度上取第n个数时，对应值的大小。利用计算式和计算式可获得fc(m,n)和两个二维矩阵，其中i、i 1、ai、a
i 1
、j、j 1、bj、b
j 1
分别是s21和s22计算获得的值位置信息权重信息和值位置信息、权重信息。
[0080]
s32，根据值、fc(m,n)和构造的列向量。
[0081]
已知输运方程求解得到的值，而fc(m,n)和二维矩阵均为m行n列的矩阵。循环遍历维度上每一个位置，即n从1到n。根据网格的值和列向量fc(：,n)，二分法查找获得位置信息(i
temp
，i
temp
1)和权重信息(pi，p
i 1
)，其中，列向量fc(：,n)为二维矩阵fc(m,n)在维度上取第n个数、在维度上分别取从第1到第m个数时形成的列向量。构造的列向量中第n个元素大小为的列向量中第n个元素大小为当结束n从1到n的遍历循环时，完成列向量的构造。
[0082]
s33，对进度变量方差预估值进行截断。
[0083]
根据s32得到的列向量，进一步求取列向量的最大值和最小值对s13获得的每个网格的进度变量方差预估值进行截断，确保其在合理的预估范围
内。
[0084]
s34，构造的代数表达式。
[0085]
进度变量方差取时，可用于表征值可能产生波动大小的中间水平，其中，结合s13和s33截断后得到的进度变量方差预估值将值代数表达式构造为其中θ为预估值的比例系数，且θ∈[0，1]，根据不同的算例构型和计算工况，在[0，1]任取一个数值，改变θ值的大小，标定不同算例的参数设置；考虑生成和耗散效应的预估值权重越大时，θ值越接近于1；需越接近于值可能产生波动大小的中间水平时，θ值越接近于0，。
[0086]
s35，查询值在列向量中的位置信息和权重信息。
[0087]
通过s34构造的值和s32获得的的列向量，利用二分法查询出值位置信息(n，n 1)和权重信息(fn，f
n 1
)。
[0088]
s36，根据值、s35获得的位置信息(n，n 1)、权重信息(fn，f
n 1
)和s31获得的二维矩阵fc(m,n)，构造的列向量。
[0089]
循环遍历维度上每一个位置，即m从1到m。构造的列向量中第m个元素大小为fc(m,n)*n c(m,n 1)*f
n 1
。当结束m从1到m的遍历循环时，完成列向量的构造。
[0090]
s37，查询值在列向量中的位置信息和权重信息。
[0091]
通过输运方程已求解的值和s36获得的的列向量，利用二分法查询出值位置信息(m，m 1)和权重信息(em，e
m 1
)。
[0092]
s4，查表更新化学热力学物性场，更新流体域密度场，完成湍流燃烧模型计算流程。
[0093]
通过获得的四个变量的位置信息和权重信息，查询火焰面数据库得到化学热力学相关状态参数，并更新流体密度，实现燃烧过程的物性更新，从而完成湍流燃烧模型的计算流程。具体参照图6：
[0094]
s41，根四个量的位置信息和权重信息，查表更新化学热力学物性场。基于以上得到的四个变量在其所对应列向量中所处的位置信息和权重信息，查询湍流火焰面数据库f(k1,k2,k3,k4)，利用计算式得到流场全部化学热力学相关状态参数值，其中表示火焰面数据库中任一化学热力学状态参数，如温度、组分质量分数、比热容、导热率、密度、混合气体的分子粘度。
[0095]
s42，更新流体域密度场。
[0096]
根据s41更新的温度和算例设置的参考压力值，利用气体状态方程，计算得到密度场数据ρ
计算
；通过与前一步密度场数据的松弛计算(ρ
n 1
＝α*ρ
计算
(1-α)ρn，第一步计算时，ρ1＝ρ
计算
；其中α为松弛因子，n为第n时间步，也就是“前一步”)，更新流体域密度场，完成湍流燃烧模型计算流程。
[0097]
本发明通过构建进度变量方差的函数来替换对输运方程的求解，实现了对火焰面生成流形(fgm)模型应用于工程实际航空发动机燃烧室湍流燃烧仿真模拟计算过程的降维提速，同时保证了计算结果的精确度和效率，具有较高的工程适用性。所取得的具体效果如下：
[0098]
(1)通过对的代数计算式构造的降维提速方法，与四输运方程求解的结果比较上看，温度场、组分场等场分布趋势一致，场均值误差较小，可以达到四输运方程求解的计算精度水平。例如，在分别采用基于四变量输运方程求解方法和基于三变量输运方程结合代数表达式计算的求解方法，针对网格数量为5190的台阶流气相燃烧算例的仿真计算结果显示，与四变量输运方程的求解结果相比，降维提速方法计算结果的流场温度、组分、进度变量方差等变量的分布趋势均一致，且平均值的误差较小。如图7(a)、(b)所示，通过对比四输运方程和本发明方法的温度分布图，得到温度t的流体域结果平均值的误差为0.14％；如图7(c)、(d)所示，通过对比四输运方程和本发明方法的y_co2分布图，得到组分y_co2的流体域结果平均值的误差为0.32％；如图7(e)、(f)所示，通过对比四输运方程和本发明方法的y_h2o分布图，得到组分y_h2o的流体域结果平均值的误差为0.34％；如图7(g)、(h)所示，通过对比四输运方程和本发明方法的分布图，得到进度变量方差的流体域结果平均值的误差为0.29％。针对网格数量为461530的模型燃烧室气液两相燃烧算例计算结果中，流场温度、组分等变量的分布趋势也一致，且平均值的误差较小。如图8(a)、(b)所示，通过对比四输运方程和本发明方法的温度分布图，得到温度t的流体域结果平均值的误差为0.71％；如图7(c)、(d)所示，通过对比四输运方程和本发明方法的y_co2分布图，得到组分y_co2的流体域结果平均值的误差为0.61％；如图7(e)、(f)所示，通过对比四输运方程和本发明方法的y_h2o分布图，得到组分y_h2o的流体域结果平均值的误差为1.5％；如图7(g)、(h)所示，通过对比四输运方程和本发明方法的分布图，得到进度变量方差的流体域结果平均值的误差为2.8％。
[0099]
(2)相比于四变量输运求解方法，本发明的方法可明显缩短燃烧模型的计算时间。在发动机燃烧室湍流燃烧模拟中，通常燃烧模型的求解在整个湍流燃烧模拟过程中占很大的比重，仅次于n-s方程求解，因此缩短燃烧模型求解时间，是提高燃烧模拟计算效率的重要途径之一。通过模拟测试结果对比，可以发现，与fgm模型中基于要途径之一。通过模拟测试结果对比，可以发现，与fgm模型中基于四变量输运方程求解结果相比，台阶流气相燃烧算例中对于燃烧模型求解的计算时间缩短24.7％，而模型燃烧室气液两相燃烧算例中对于燃烧模型求解的计算时间缩短22.6％。
[0100]
(3)本发明的方法不仅可单独使用，也可用于配合四方程模型切换使用：在相关技术中基于伪时间方法的稳态计算通常采用四方程模型方法进行求解计算，在本发明中，进行基于伪时间方法的稳态计算时，计算开始使用本发明方法，计算结果稳定收敛后，切换为四方程模型进行稳态计算，这样既可以获得稳态时的四方程计算结果，又可加速四方程求解的迭代收敛速度，在工程实际中具有较大的应用价值。
[0101]
(4)本发明的方法也可用于其它包含但不仅限于四个控制变量的火焰面类燃烧模型(如变量的火焰面进度变量fpv模型，
总焓h五变量的非绝热fgm模型)的降维提速。通过去掉输运方程的求解，改为代数表达式计算的方法，可以获得与fgm燃烧模型改进前后相似的计算精度对比效果与提速对比效果。
[0102]
一种燃烧模拟降维提速装置，包括：
[0103]
获取进度变量方差模块，用于构造求解进度变量方差的代数式，并得到进度变量方差；
[0104]
位置权重信息获取模块，用于得到平均混合分数混合分数方差平均进度变量和进度变量方差的位置信息和权重信息；
[0105]
更新模块，基于获取的位置信息和权重信息，更新化学热力学物性场，并更新流体域密度场。
[0106]
基于与一种燃烧模拟降维提速方法相同的发明构思，本技术还提供了一种电子设备，该电子设备包括一个或多个处理器和一个或多个存储器，存储器中存储了计算机可读代码，其中，计算机可读代码当由一个或多个处理器执行时，进行本发明一种燃烧模拟降维提速方法的实施。其中，存储器可以包括非易失性存储介质和内存储器；非易失性存储介质可存储操作系统和计算机可读代码。该计算机可读代码包括程序指令，该程序指令被执行时，可使得处理器执行任意一种燃烧模拟降维提速方法。处理器用于提供计算和控制能力，支撑整个电子设备的运行。存储器为非易失性存储介质中的计算机可读代码的运行提供环境，该计算机可读代码被处理器执行时，可使得处理器执行任意一种燃烧模拟降维提速方法。
[0107]
应当理解的是，处理器可以是中央处理单元(central processing unit，cpu)，该处理器还可以是其他通用处理器、数字信号处理器(digital signal processor，dsp)、专用集成电路(application specific integrated circuit，asic)、现场可编程门阵列(field-programmable gate array，fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。其中，通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。
[0108]
其中，所述计算机可读存储介质可以是前述实施例所述电子设备的内部存储单元，例如所述计算机设备的硬盘或内存。所述计算机可读存储介质也可以是所述电子设备的外部存储设备，例如所述电子设备上配备的插接式硬盘、智能存储卡(smartmedia card，smc)、安全数字(secure digital，sd)卡、闪存卡(flash card)等。
[0109]
所述实施例为本发明的优选的实施方式，但本发明并不限于上述实施方式，在不背离本发明的实质内容的情况下，本领域技术人员能够做出的任何显而易见的改进、替换或变型均属于本发明的保护范围。