一种求解化学非平衡流动化学反应源项的预估修正方法与流程

本发明涉及流体力学数值计算技术领域，尤其涉及一种求解化学非平衡流动化学反应源项的预估修正方法。

背景技术：

在利用计算流体力学技术求解高超声速化学非平衡流动研究中，需要解决因化学反应特征时间尺度与流动不匹配而引起的刚性问题，其中化学反应生成率及各组分密度是所关注的重要参数。时间推进的lusgs隐式算法是可压缩流动计算的常规方法，而将化学反应源项作隐式处理的全隐式算法虽然有较高的稳定性和收敛性，但会造成系数矩阵的非对角化，破坏原有算法的优越性。因此，针对高超声速化学非平衡流动的数值计算，需要发展具有较高收敛性和计算效率的求解化学非平衡流动化学反应源项的方法。

技术实现要素：

本发明所要解决的技术问题是：针对上述存在的问题，提供了一种求解化学非平衡流动化学反应源项的预估修正方法，本方法借鉴流动与化学反应解耦的物理分裂思想，假设化学非平衡过程是由系列局部等压非平衡子过程组成，冻结组分连续性方程的对流和扩散项，分层次对时间步作子过程剖分，并在不同层次对预估值进行修正，进而获得下一时刻的化学反应源项。

本发明采用的技术方案如下：一种求解化学非平衡流动化学反应源项的预估修正方法，包括以下步骤：

步骤1：根据当前时刻的流场，解获得到当前时刻的参数值，设定该计算轮次为初始轮次；

步骤2：根据当前时刻已知的参数值，开展下一轮次的预估，先通过二分法剖分时间离散间隔，然后依次对二分点的参数值进行预估，并对上一轮次各二分点的参数值进行修正；

步骤3：计算第n+1轮和第n轮的预估值之差，若差值大于设定的阈值，则返回步骤2，增加一轮预估修正；否则跳转至下一步骤，其中n为大于等于1的自然数；

步骤4：将第n+1轮计算得到的参数值作为最终结果。

这样，通过综合采用核心控制方程隐式处理和化学反应源项半隐式处理的方式，在保持原有算法优势的情形下，同时解决了化学反应源项的刚性问题，且程序实现简单。

进一步的，所述参数值包括组分密度和化学反应生成率。

进一步的，所述步骤2具体包括以下步骤：

步骤21：根据当前时刻已知的组分密度和化学反应生成率，计算下一时刻的组分密度；

步骤22：进入下一轮次预估修正；

步骤23：对上一轮次2^n-1个子过程作二分，形成2ⁿ个子过程；

步骤24：设定子过程时刻i＝1，开始按子过程时刻进行推进；

步骤25：求解第n+1轮二分点时刻2i-1的组分密度和化学反应生成率，其中，组分密度根据第n+1轮上一子时刻2i-2和第n轮对应下一子时刻i的相应值共同求得；对应的化学反应生成率由化学反应动力学获得；

步骤26：求解第n+1轮整点2i时刻的组分密度和化学反应生成率，其中，组分密度根据第n+1轮上一子时刻2i-1和第n轮子时刻i的相应值共同求得，对应的化学反应生成率由化学反应动力学获得；

步骤27：若当前子过程为该轮次最后子过程，即i＝2^n-1，则进入下一步；否则子过程时刻继续推进，i＝i+1，并返回步骤25；

步骤28：计算完当前轮次中的所有子过程，获得基于本轮次预估修正的下一时刻组分密度和化学反应生成率。

进一步的，所述步骤25中，第n+1轮二分点位置的组分密度计算公式如下：

其中，表示第n+1轮前2i-2时刻的组分密度，表示第(n+1)轮前2i-2时刻的化学反应生成率，表示第n轮对应下一子时刻i的化学反应生成率。

进一步的，所述步骤26中，第n+1轮整点位置的组分密度计算公式如下：

其中，表示第n+1轮上一子时刻2i-1的组分密度，表示第n+1轮上一子时刻2i-1的化学反应生成率，表示第n轮对应下一子时刻i的化学反应生成率。

与现有技术相比，采用上述技术方案的有益效果为：本发明的求解化学反应源项的预估修正方法，对核心控制方程隐式处理、化学反应源项半隐式处理，保持原有lusgs算法的优势，同时高效准确地求解化学反应源项，程序实现简单。

附图说明

图1是本发明的求解化学反应源项的预估修正实施示意图。

图2是本发明的求解化学非平衡流动化学反应源项的预估修正方法的流程图。

具体实施方式

本发明提供一种求解化学非平衡流动化学反应源项的预估修正方法。本发明的求解化学非平衡流动化学反应源项的预估修正方法借鉴流动与化学反应解耦的物理分裂思想，假设化学非平衡过程是由系列局部等压非平衡子过程组成，冻结组分连续性方程的对流和扩散项，分层次对时间步作子过程剖分，并在不同层次对预估值进行修正，进而获得下一时刻的化学反应源项。

下面结合附图进一步对本发明的求解化学非平衡流动化学反应源项的预估修正方法做详细说明。

如图1所示，在已知当前时刻t0的流场，借鉴流动与化学反应解耦的物理分裂思想，假设化学非平衡过程是由系列局部等压非平衡子过程组成，并按照二分法对时间间隔δt做均匀剖分，第1轮生成1个子过程(2^1-1)，第2轮生成2个子过程(2^2-1)，……，第n轮生成2^n-1个子过程。按照预估修正规则对各轮次的子过程作计算。

对照如图2所示的流程图给出具体实施步骤。

步骤01：基于当前时刻流场准备求解化学反应源项所需的变量。

步骤02：获取当前时刻t0的组分密度ρ⁰和化学反应生成率ω⁰。

步骤03：设定初始轮次n＝1。

步骤04：计算下一时刻t0+δt的组分密度ρ1⁽¹⁾，即

对应的化学反应生成率ω1⁽¹⁾由化学反应动力学获得。

步骤05：在已知第n轮各子过程的解情况下，可对第n+1轮各子过程进行预估修正。所以，进入下一轮预估修正，n＝n+1。

步骤06：对上一轮次2^n-1个δt子过程作二分，形成2ⁿ个δt子过程。

步骤07：按子过程时刻顺序开始推进，设定子过程初始时刻i＝1。

步骤08：已获得第(n+1)轮前(2i-2)时刻的组分密度和化学反应生成率求解第n+1轮二分点时刻(2i-1)的组分密度和化学反应生成率，其中组分密度根据第(n+1)轮上一子时刻(2i-2)和第n轮对应下一子时刻(i)的相应值共同求得，如公式(1)

对应的化学反应生成率由化学反应动力学获得。

步骤09：对于第n+1轮整点2i时刻，其组分密度根据第n+1轮上一子时刻(2i-1)和第n轮该子时刻(i)的相应值共同求得，如公式(2)

对应的化学反应生成率由化学反应动力学获得。

步骤10：进行判断

若当前子过程为该轮次最后子过程，即i＝2^n-1，则进入步骤11；否则，i＝i+1，并返回步骤08计算下一子过程。

步骤11：计算完成当前轮次中所有子过程后，获得基于本轮次预估修正的下一时刻组分密度和化学反应生成率。

步骤12：进行判断

若第n+1轮和第n轮的化学反应生成率之差足够小，则进入步骤13；否则，n＝n+1，增加一轮预估修正。

步骤13：最终获得下一时刻的组分密度和化学反应生成率，将其作为最终计算结果，结束化学反应源项的预估修正计算。

本发明的求解化学非平衡流动化学反应源项的预估修正方法综合采用核心控制方程隐式处理、化学反应源项半隐式处理，保持原有算法的优势，同时解决化学反应源项的刚性问题，程序实现简单。本发明的求解化学非平衡流动化学反应源项的预估修正方法可快速高效求解化学非平衡流动问题，具有较好的计算精度和效率。

本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。如果本领域技术人员，在不脱离本发明的精神所做的非实质性改变或改进，都应该属于本发明权利要求保护的范围。