用于挤代式供应系统性能研究的数值分析方法

1.本发明属于能源与动力领域，还涉及用于挤代式供应系统性能研究的数值分析方法。


背景技术：

2.热动力系统是水下航行器的“心脏”，其性能与水下航行器的各项运动指标都密切相关。其中，如图1所示，能源供应系统作为水下航行器热动力系统中的重要组成部分，用以保证将推进剂在规定时间内按一定比例快速稳定的运输至燃烧室中进行燃烧，若不能够在指定时间内在燃料泵前建立指定压强就会导致水下航行器偏离启动工况运行甚至启动失败的情况发生。
3.为保证航行器能够在水下正常运行，航行器制造者必须采用实验测试或计算机数值模拟等方法证明其的可行性。其中，实验测试是最直接的测试手段，完全在自然条件下进行，但这种方法成本较高，且无法模拟达到航行器所遇到的各种极端或突发条件。在此情况下，采用数值模拟软件逐渐成为了计算挤代式供应系统内部流场的支持工具，先进的模拟方法还可以用来低成本模拟各种极端条件下的的流动情况。在此条件下，一类采用fluent、cfx、star-ccm来模拟仿真的流场的商业软件逐渐发展起来。但其对于挤代式供应系统的模拟仿真常常会以较大的网格数量以及高性能的计算机性能为基础，从而导致其很难用于对于大量不同模型的数值仿真测试，以及短时间内快速对于指定截面流体信息的检测。


技术实现要素：

4.本发明的目的是提供用于挤代式供应系统性能研究的数值分析方法，能有效提高计算效率，更高效的检测流动过程中挤代式供应系统中各个位置流体参数的变化。
5.本发明所采用的技术方案是，用于挤代式供应系统性能研究的数值分析方法，具体按照以下步骤实施：
6.步骤1、根据给定的挤代式供应系统模型尺寸初始化各部分网格，得到网格单元，并给对应的网格单元内初始化赋予储存信息，并设定分析时间；
7.步骤2、通过ausmdv方法计算得到网格单元相邻界面上的状态通量；
8.步骤3、由连续性方程、动量方程、能量方程计算得到网格单元内部的流体信息；
9.步骤4、根据某时刻网格单元内部的流体信息采用欧拉预测矫正公式在给定的时间范围内计算下一时刻的流体信息；
10.步骤5、返回步骤2，直到计算总时间步达到设定分析时间，得到各个位置流体信息随时间的变化。
11.本发明的特点还在于：
12.步骤1中根据给定的挤代式供应系统模型尺寸初始化各部分网格具体过程为：采用有限体积法将挤代式供应系统模型离散为多个网格单元。
13.给对应的网格单元内初始化赋值过程为：根据每个网格单元的位置信息赋予储存
信息，储存信息包括边界条件、时间步长、中心位置坐标、交界面位置坐标、流体信息的平均值以及网格单元质量动量和能量的时间导数。
14.步骤2具体过程为：假设网格单元内流体信息为线性变化，结合网格单元内初始化赋予的流体信息利用网格单元内的状态平均值插值到网格单元内表面上，得到每个网格单元内表面左端和右端两侧的质量、动量和能量流动状态，即网格单元相邻界面上的状态通量，插值是针对网格单元相邻界面上的流体信息单独执行的，得到的流体信息经过整合后变为相邻界面上的状态通量。
15.流体信息包括网格单元密度、速度、温度、压力、内能、音速。
16.步骤3中连续性方程表达式为：
[0017][0018]
式中，u
x
为轴线方向的速度分量；t为时间；ρ为密度，表示网格单元内随时间变化的质量部分，表示轴线方向上网格单元相邻界面流入流出的质量。
[0019]
步骤3中动量方程表达式为：
[0020][0021]
式中，p为网格单元内部的压强，τ是因分子黏性作用而产生的作用在微元体表面上的黏性应力τ的分量；f为单位质量力。
[0022]
步骤3中能量方程表达式为：
[0023][0024]
式中，e为网格单元内部的总能，包含内能、动能和势能之和，h为焓；其中t
ref
＝298.15k；k
eff
为有效热传导系数，k
eff
＝k k
t
，k
t
为湍流热传导系数。
[0025]
步骤4中欧拉预测矫正公式表示为：
[0026][0027][0028]
其中，x是网格单元的空间步长，包含控制体j单元内的质量、动量、能量；包含控制体前后j
±
1/2单元壁面上的质量通量、动量通量、能量通量；是由上一次迭代计算出的控制体内部的状态参数得到的控制体内储存的质量、动量、能量。
[0029]
本发明有益效果是：
[0030]
本发明用于挤代式供应系统性能研究的数值分析方法，能够快速准确的计算出系
统各个位置的流体信息参数，并且可以通过更改系统尺寸，达到快速计算的目的，相比于传统的流体仿真软件极大的节约了时间成本。
附图说明
[0031]
图1是本发明所提供的用于挤代式供应系统的原理框图；
[0032]
图2是本发明所提供的模拟方法的实现流程图；
[0033]
图3是本发明所提供的用于挤代式供应系统的二维结构图；
[0034]
图4是本发明所提供的基于matlab构建的用户交换界面；
[0035]
图5是本发明仿真得到的高压气舱末端压力变化曲线；
[0036]
图6是本发明仿真得到的推进剂舱末端压力变化曲线。
具体实施方式
[0037]
下面结合附图和具体实施方式对本发明进行详细说明。
[0038]
本发明用于挤代式供应系统性能研究的数值分析方法，具体按照以下步骤实施：
[0039]
步骤1、根据给定的挤代式供应系统模型尺寸，采用有限体积法将挤代式供应系统模型离散为多个网格单元，得到网格单元，并根据每个网格单元的位置信息赋予储存信息，储存信息包括边界条件、时间步长、中心位置坐标、交界面位置坐标、流体信息的平均值以及网格单元质量、动量和能量的时间导数，并设定分析时间，其中流体信息包括网格单元密度、速度、温度、压力、内能、音速；
[0040]
步骤2、通过ausmdv方法计算得到网格单元相邻界面上的状态通量；具体过程为：假设网格单元内流体信息为线性变化，结合网格单元内初始化赋予的流体信息利用网格单元内的状态平均值插值到网格单元内表面上，得到每个网格单元内表面左端和右端两侧的质量、动量和能量流动状态，即网格单元相邻界面上的状态通量，插值是针对网格单元相邻界面质上的流体信息单独执行的，得到的流体信息经过整合后变为相邻界面上的状态通量。例如，使用下面的表达式通过非线性插值得到内表面左右任意一侧的密度：
[0041][0042][0043][0044][0045]
其中，(δ-)j和(δ )j是对网格(控制体)j密度斜率的两种估计方式，x
j-1
、xj、x
j 1
是网格(控制体)j-1、j、j 1相对应的中心点位置。如果这两个密度斜率具有相同的符号，则minmod限制器函数将选择最小幅度的密度函数，否则返回零。
[0046]
步骤3、由连续性方程、动量方程、能量方程计算得到网格单元内部的密度、速度、温度、压力、内能、音速；
[0047]
连续性方程表达式为：
[0048][0049][0050]
式中，u
x
为轴线方向的速度分量；t为时间；ρ为密度，表示网格单元内随时间变化的质量部分，表示轴线方向上网格单元相邻界面流入流出的质量。
[0051]
动量方程表达式为：
[0052][0053]
式中，p为网格单元内部的压强，τ是因分子黏性作用而产生的作用在微元体表面上的黏性应力τ的分量；f为单位质量力。
[0054]
能量方程表达式为：
[0055][0056]
式中，e为网格单元内部的总能，包含内能、动能和势能之和，h为焓；其中t
ref
＝298.15k；.k
eff
为有效热传导系数，.k
eff
＝k k
t
，k
t
为湍流热传导系数。
[0057]
根据动量方程、能量方程、连续性方程计算出的网格单元的内的质量、动量和能量，能够用来更新网格单元内的流体信息：
[0058]
ρj＝u
mass,j
[0059][0060][0061]
其中u
mass,j
、u
momentum,j
、u
energy,j
分别为j网格单元上的质量，动量、能量；vj为单元网格上流体的速度，ej为网格单元上流体的比内能。挤代系统中氮气可被认为是热完全气体，则其比内能与温度成正比，表示为：
[0062][0063]
其中，cv是气体的定容比热容，继而根据已求解出的密度与温度通过理想气体状态方程可求出压力：
[0064]
pj＝ρjrtj[0065]
其中，r是气体常数。
[0066]
声速为：
[0067]
[0068]
其中，aj为网格单元上的音速，γ为比热比。
[0069]
步骤4、根据某时刻网格单元内部的流体信息(流体信息包括网格单元密度、速度、温度、压力、内能、音速)采用欧拉预测矫正公式在给定的时间范围内计算下一时刻的网格单元密度、速度、温度、压力、内能、音速；
[0070]
欧拉预测矫正公式表示为：
[0071][0072][0073]
其中，x是网格单元的空间步长，包含控制体j单元内的质量、动量、能量；包含控制体前后j
±
1/2单元壁面上的质量通量、动量通量、能量通量；是由上一次迭代计算出的控制体内部的状态参数得到的控制体内储存的质量、动量、能量。
[0074]
步骤5、返回步骤2，直到计算总时间步达到设定分析时间，得到各个位置密度、速度、温度、压力、内能、音速随时间的变化。
[0075]
实施例
[0076]
如图1所示，能源供应系统作为水下航行器热动力系统中的重要组成部分，用以保证将推进剂在规定时间内按一定比例快速稳定的运输至燃烧室中进行燃烧，若不能够在指定时间内在燃料泵前建立指定压强就会导致水下航行器偏离启动工况运行甚至启动失败的情况发生。
[0077]
用于挤代式供应系统的快速数值模拟方法，这种数值模拟技术可以用matlab编程实现，并通过计算机运行来模拟挤代式供应系统内部的流动状态。
[0078]
模拟的条件是：
[0079]
高压气舱直径：35e-3m高压气舱长度：1.6m
[0080]
活门直径：6e-3m活门长度：6e-3m
[0081]
管路1直径：15e-3m管路1长度：50e-3m
[0082]
管路2直径：20e-3m管路2长度：500e-3m
[0083]
推进剂舱容积：0.024l
[0084]
初始压力：0.1mpa初始温度：300k
[0085]
高压气舱压力：20mpa高压气舱温度：300k
[0086]
如图2所示在确定模型尺寸后，根据模型结构特点进行对方程进行简化，将如图3所示模型离散为相邻的网格单元，并将初始化的流体信息赋值在其网格中心点，通过ausmdv方法计算出网格交界面的状态通量后，再通过连续性方程、动量方程、能量方程计算出各个网格单元内的流体信息，最后通过欧拉预测矫正公式计算得到下一时刻的网格单元内的流体信息参数。
[0087]
根据此程序开发了图4所示简易的用户界面，用于修改得到不同尺寸下系统内部流体流动的过程，通过计算后最终得到如图5和图6所示的系统模型高压气舱末端和推挤机舱中部的压力曲线。
[0088]
图5和图6是根据程序计算得到的不同位置压力随时间的变化曲线，相比于传统商
业软件计算迅速，数值准确。
[0089]
通过上述方式，本发明用于挤代式供应系统性能研究的数值分析方法，能够快速准确的计算出系统各个位置的流体信息参数，并且可以通过更改系统尺寸，达到快速计算的目的，相比于传统的流体仿真软件极大的节约了时间成本。