CN115618171B - 一种基于同伦算法的推进剂燃烧平衡产物求解方法 - Google Patents

Abstract

本发明公开了一种基于同伦算法的推进剂燃烧平衡产物求解方法，包括以下步骤：根据推进剂的燃烧的实际情况，选择合适的非线性燃烧模型，获取推进剂燃烧组分模型；根据燃烧过程需求，构造满足条件的同伦方程，获得具有边界条件的非线性常微分求解方程；选取满足要求的设置求解燃烧产物的初始值；设置迭代过程的推进步长、迭代精度与最大步数；选取合适的迭代求解代数方程，通过循环算法，求得满足要求的最优近似解。本发明基于同伦算法的思想，使用基于拓扑学同伦原理的数值解法，在不改变方程组形式的情况下直接求解，并且同伦解法求解过程对初值不敏感，从定义域的任意初始值开始，通过迭代最终都可求得收敛解，是一种稳定快速的求解方法。

Description

一种基于同伦算法的推进剂燃烧平衡产物求解方法

技术领域

本发明涉及固体推进剂技术领域，具体涉及一种基于同伦算法的推进剂燃烧平衡产物求解方法。

背景技术

随着现代战争对武器系统在远程化、机动性等方面的要求越来越高，以固体火箭发动机为动力源的火箭、导弹武器发挥着重要的作用。推进剂作为控制固体火箭发动机的重要动力来源，燃烧过程的分析对发动机工作具有重要的意义。

传统的推进剂燃烧平衡产物求解方法是将对数项进行泰勒级数展开，进行线性化近似处理后通过牛顿迭代法或最速梯度下降法进行求解。使用这种方法可以使求解方程得到简化，但是通过泰勒级数展开取前两项近似，将要丢失一定的精度。并且这两种数值方法对求解初始值都具有一定的敏感性，若初值设置不当则有可能造成求解发散。

发明内容

有鉴于此，本发明提供了一种基于同伦算法的推进剂燃烧平衡产物求解方法，是一种稳定快速的求解方法。

为达到上述目的，本发明的技术方案包括如下步骤：

S1根据推进剂的燃烧的实际情况，确定推进剂组员构成以及燃烧产物组分。

S2根据吉布斯最小自由能法，建立燃烧平衡状态的守恒控制方程。

S3根据燃烧过程需求，构造满足条件的同伦方程，转化为常微分线性控制方程组。

S4选取满足要求的设置求解燃烧产物的初始值x(0)；设置迭代过程的推进步长h、迭代精度e与最大步数N；

S5选取合适的迭代求解代数方程，通过循环算法，求得满足要求的最优近似解。

进一步地，S1根据推进剂的燃烧的实际情况，确定推进剂组员构成以及燃烧产物组分，具体为：

所述推进剂为lkg固体推进剂，其化学式为摩尔质量为1000g/mol，燃烧反应平衡时组分数量为下标N，第i种组分的摩尔数为n_i，在此种组分中第j种元素的摩尔原子数为A_ij。

则根据质量守恒定律得式中，M为推进剂中含有的元素种类总数。

设燃烧产物中第i种组分的摩尔吉布斯自由能为G_m,i，则在反应达到平衡时N种组分的吉布斯自由能表达式为

进一步地，根据吉布斯最小自由能法，建立燃烧平衡状态的守恒控制方程，具体为：选取合适的燃烧状态平衡方程，构建吉布斯自由能平衡方程，设置求解燃烧平衡产物，得到燃烧平衡状态的守恒控制方程。

进一步地，根据燃烧过程需求，构造满足条件的同伦方程，转化为常微分线性控制方程组，依据所述吉布斯自由能平衡方程，设置x为包含所有燃烧产物组分物质的量n_i，以及拉格朗日乘子λ_i的向量，FF(x)为包含所有组分平衡方程的向量函数，表述为

式中，向量x属于空间X，此空间是N+M维实数空间的子集，在空间X中应满足条件/>M为推进剂中含有的元素种类总数；λ_j表示第j种元素的拉格朗日乘子，n_i表示第i种组分的摩尔数，向量FF(x)属于实数空间/>x₁,x₂,...,x_N+M为第1～第N+M个燃烧产物组分；p为压强，单位为一个标准大气压(atm)；n_g为所有气体产物的摩尔数总和；R_u为摩尔气体常数。

进一步地，选取合适的迭代求解代数方程，通过循环算法，求得满足要求的最优近似解，具体为：

根据S4设置的设置求解燃烧产物的初始值x(0)；设置迭代过程λ的推进步长h、迭代精度e与最大步数N；选用四阶Runge-Kutta方法进行迭代推进求解，包括如下步骤：

SS1：进入Runge-Kutta方法计算步骤。

SS2：计算当前步的Jacobian矩阵J(x(λ_i))。

SS3：使用Gauss-Siedel迭代求解代数方程J(x(λ_i))k_j＝FF(x(0))；k_j为当前迭代步的反应速率，FF(x(0))为包含所有组分平衡方程的向量函数初始值。

SS4：判断是否满足迭代精度e，若是则执行SS5。

否则判断是否达到最大迭代步数N，若是则执行SS5，否则返回SS3。

SS4：根据计算结果k_j更新x。

SS5：j自增1；判断j是否大于4，若是则执行SS6，否则返回SS2。

SS6：判断λ_i是否为1；若是则输出x(1)，结束本流程；否则λ_i+1＝λ_i+h，返回SS1。

有益效果：

1、本发明提出了一种基于同伦算法的推进剂燃烧平衡产物求解方法，使用基于拓扑学同伦原理的数值解法，在不改变方程组形式的情况下可以直接求解，并且同伦解法求解过程对初值不敏感，从定义域的任意初始值开始，通过迭代最终都可以求得收敛解，是一种稳定快速的求解方法。

2、本发明基于同伦算法，提出了一种新的推进剂燃烧平衡产物求解方法，本方法在求解非线性方程过程中无需降阶，选取初值对结果影响较小，结果精度更高，求解过程更稳定。

3、本方面选用四阶Rungc-Kutta方法对上进行迭代推进求解，由于向量x中的元素较多在每一步中使用Gauss-Sicdel迭代求解代数方程，实现了快速稳定的求解。

附图说明

图1为本发明实施例提供的一种基于同伦算法的推进剂燃烧平衡产物求解方法流程图。

图2为本发明实施例中的基于同伦算法的火箭发动机燃烧平衡产物求解方法流程图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本发明提供了一种基于同伦算法的推进剂燃烧平衡产物求解方法，具体流程如图1所示，包括以下步骤：

步骤1：根据推进剂的燃烧的实际情况，选择合适的非线性燃烧模型，获取推进剂燃烧组分模型；

本实施例中，设推进剂为lkg固体推进剂，其假定化学式为摩尔质量为1000g/mol，燃烧反应平衡时组分数量为下标N，第i种组分的摩尔数为n_i，在此种组分中第j种元素的摩尔原子数为A_ij，则根据质量守恒定律可得/>式中，M为推进剂中含有的元素种类总数。设燃烧产物中第i种组分的摩尔吉布斯自由能为G_m,i，则在反应达到平衡时N种组分的吉布斯自由能表达式为/>将式

带入上式，可得式

式中对于1mol气体是标准自由生成焓，即在标准状态下由化学单质合成该物质时自由烩的变化。G_m是在标准压强下该气体的自由焓，是与温度相关的函数。式中，p_i为第i种组分的分压其量纲为atm，根据理想气体状态方程有/>式中，n_g为所有气体产物的摩尔数总和，有n_g＝∑n_i，代入式

可得式

引入函数Φ并令当反应达到平衡状态时，产物组分的吉布斯自由能之和达到最小值，所以求平衡时各组分的摩尔数就等价于求解目标函数中最小值的最优化问题，约束条件为式

使用拉格朗日乘子法构建新的目标函数F，将有条件极值问题转化为无条件极值问题则可得式

当目标函数F达到最小值时，应满足以下条件:

整理之后可得:

上式即为燃烧产物的组分平衡方程，通过求解该方程可得在不同压强和温度下的推进剂表面燃烧产物的组分构成。

选取推进剂燃烧产物的组分平衡方程为

式中，G_m,i表示燃烧产物中第i种组分的摩尔吉布斯自由能，R_u表示普适气体常量，T表示温度，n_i表示第i种组分的摩尔数，p表示发动机内部压强，n_g表示所有气体产物的摩尔数总和，λ_j表示第j种元素的拉格朗日乘子，A_ji表示第j种元素的摩尔原子数，N_j表示第j种元素燃烧反应平衡时的组分数量，M为推进剂中含有的元素种类总数。

步骤2：根据燃烧过程需求，构造满足条件的同伦方程，获得具有边界条件的非线性常微分求解方程；

在本实施例中，设x为包含所有燃烧产物组分物质的量n_i；以及拉格朗日乘子λ_i的向量，FF(x)为包含所有组分平衡方程的向量函数，则式

可以表述为式

向量x属于空间X，此空间是N+M维实数空间的子集，在空间X中应满足条件

向量FF(x)属于实数空间

步骤3：选取满足要求的设置求解燃烧产物的初始值x(0)。

引入参数λ建立从空间X映射到空间中的参数函数G(λ,x)，则有式

式中，x(0)表示初始值，当λ＝0时有G(0,x)＝FF(x)-FF(x(0))，当λ＝1时有G(1,x)＝FF(x)，通过参数λ的连续变化函数G就可以从函数FF(x)-FF(x(0))变换为函数FF(x)，因此函数G是函数FF(x)-FF(x(0))与函数FF(x)之间的同伦。令G(λ,x)＝0则x是λ取[0,1]区间任意值时G＝0的解，对于待求解的组分平衡方程，因为当吉布斯自由能达到最小值时系统的平衡状态有且仅有一种情况，所以对于任意λ值时G＝0的解都是唯一的。因此可以看出解x与λ的取值有关，可以将x表示为λ的函数x(λ)，当λ＝0时x＝x(0)是初始值，当λ＝1时x＝x(1)，x(1)是函数FF(x)＝0的解即是待求解值。构成函数FF(x)的子函数f_i(x)在空间X内都是连续且可微的，因此函数FF(x)以及由此构成的函数G(λ,x)也是连续且可微的。对函数G(λ,x)＝0取λ的导数，则可得式

式中，x'(λ)是函数x(λ)关于λ的导数，将x(λ)看做在空间[0,1]×X内的曲线，因为函数FF(x)-FF(x(0))到函数FF(x)的变换是连续的，所以曲线x(λ)是连续且光滑的对参数入可微。对式

进行变化可得式

根据上式对于G(λ,x)的定义，以及向量微分对Jacobian的定义，可得式

函数G(λ,x)对λ的偏微分为

整理后得式

x′(λ)＝-[J(x(λ))]^-1FF(x(0))

步骤4：设置迭代过程的推进步长h、迭代精度e与最大步数N；

本实施例中，设置推进步长为h＝0.01，迭代精度e＝1×10^-6，最大迭代步数N＝1×10⁶。

步骤5：选取合适的迭代求解代数方程，通过循环算法，求得满足要求的最优近似解。

本文选用四阶Rungc-Kutta方法对x′(λ)＝-[J(x(λ))]^-1FF(x(0))进行迭代推进求解，由于向量x中的元素较多在每一步中使用Gauss-Sicdel迭代求解代数方程，整体解算流程如图1所示。

综上，本发明公开了一种基于同伦算法的推进剂燃烧平衡产物求解方法，该方法包括以下步骤：

步骤1：根据推进剂的燃烧的实际情况，选择合适的非线性燃烧模型，获取推进剂燃烧组分模型；

步骤2：根据燃烧过程需求，构造满足条件的同伦方程，获得具有边界条件的非线性常微分求解方程；

步骤3：选取满足要求的设置求解燃烧产物的初始值x(0)。

步骤4：设置迭代过程的推进步长h、迭代精度e与最大步数N。

步骤5：选取合适的迭代求解代数方程，通过循环算法，求得满足要求的最优近似解。

本实施例中，对非线性方程组

的求解通过拓扑同伦理论变换为了对式

x′(λ)＝-[J(x(λ))]^-1FF(x(0))

常微分方程初值问题的求解。本文选用四阶Rungc-Kutta方法对上进行迭代推进求解，由于向量x中的元素较多在每一步中使用Gauss-Sicdel迭代求解代数方程，整体解算流程如图2所示。

与传统的火箭发动机燃烧平衡产物求解方法相比，本发明的优势在于：本发明中基于同伦算法的思想，使用基于拓扑学同伦原理的数值解法，在不改变方程组形式的情况下可以直接求解，并且同伦解法求解过程对初值不敏感，从定义域的任意初始值开始，通过迭代最终都可以求得收敛解，是一种稳定快速的求解方法。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示意性实施例”、“示例”、“具体示例”或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改﹑替换和变型，本发明的范围由权利要求及其等同物限定。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于同伦算法的推进剂燃烧平衡产物求解方法，其特征在于，包括如下步骤：

S1根据推进剂的燃烧的实际情况，确定推进剂组员构成以及燃烧产物组分；具体为：

所述推进剂为lkg固体推进剂，其化学式为摩尔质量为1000g/mol，燃烧反应平衡时组分数量为下标N，第i种组分的摩尔数为n_i，在此种组分中第j种元素的摩尔原子数为A_ij；

则根据质量守恒定律得式中，M为推进剂中含有的元素种类总数；

设燃烧产物中第i种组分的摩尔吉布斯自由能为G_m,i，则在反应达到平衡时N种组分的吉布斯自由能表达式为

S2根据吉布斯最小自由能法，建立燃烧平衡状态的守恒控制方程；具体为：

选取合适的燃烧状态平衡方程，构建吉布斯自由能平衡方程，设置求解燃烧平衡产物，得到燃烧平衡状态的守恒控制方程；

S3根据燃烧过程需求，构造满足条件的同伦方程，转化为常微分线性控制方程组，依据所述吉布斯自由能平衡方程，设置x为包含所有燃烧产物组分物质的量n_i，以及拉格朗日乘子λ_i的向量，FF(x)为包含所有组分平衡方程的向量函数，表述为

式中，向量x属于空间X，此空间是N+M维实数空间的子集，在空间X中应满足条件M为推进剂中含有的元素种类总数；λ_j表示第j种元素的拉格朗日乘子，n_i表示第i种组分的摩尔数，向量FF(x)属于实数空间/>x₁,x₂,...,x_N+M为第1～第N+M个燃烧产物组分；p为压强，单位为一个标准大气压(atm)；n_g为所有气体产物的摩尔数总和；R_u为摩尔气体常数；

S4选取满足要求的设置求解燃烧产物的初始值x(0)；设置迭代过程的推进步长h、迭代精度e与最大步数N；

S5选取合适的迭代求解代数方程，通过循环算法，求得满足要求的最优近似解；

选用四阶Runge-Kutta方法进行迭代推进求解，包括如下步骤：

SS1：进入Runge-Kutta方法计算步骤；

SS2：计算当前步的Jacobian矩阵J(x(λ_i))；

SS3：使用Gauss-Siedel迭代求解代数方程J(x(λ_i))k_j＝FF(x(0))；k_j为当前迭代步的反应速率，FF(x(0))为包含所有组分平衡方程的向量函数初始值；

SS4：判断是否满足迭代精度e，若是则执行SS5；

否则判断是否达到最大迭代步数N，若是则执行SS5，否则返回SS3；

SS4：根据计算结果k_j更新x；

SS5：j自增1；判断j是否大于4，若是则执行SS6，否则返回SS2；

SS6：判断λ_i是否为1；若是则输出x(1)，结束本流程；否则λ_i+1＝λ_i+h，返回SS1。