CN105808916A - 燃烧室试验台虚拟试验建模方法 - Google Patents

Abstract

本发明公开了一种燃烧室试验台虚拟试验建模方法。其包括以下步骤：步骤1：采用有限体积法将试验模型划分成多个有限控制体，建立所述有限控制体的数学模型方程；步骤2：根据所述数学模型方程建立仿真模型，以进行并行运算，获得所述试验模型的相关参数。采用这种燃烧室试验台虚拟试验建模方法可以较为逼真地模拟发动机工作的真实情况，提高燃烧室试验台虚拟试验的准确度。

Description

燃烧室试验台虚拟试验建模方法

技术领域

本发明涉及发动机虚拟试验技术领域，尤其涉及一种燃烧室试验台虚拟试验建模方法。

背景技术

建立涡轮发动机的虚拟试验可以减少燃烧室实物试验次数，缩短发动机、尤其是燃烧室的研制周期，从而可以极大地降低涡轮发送机的研制成本。

但是，现有技术中所进行的虚拟试验大都是在理想状况下，没有考虑燃烧过程中的实际情况。其实，燃烧室的实际燃烧过程为带化学反应的气粒两相流动过程，由此燃烧室、进气管道等部件的气动热力参数，例如，燃烧室内的总压强、总温度、流量和油气比等的计算较为复杂，并且计算的复杂度也成倍提高，从而对虚拟试验的计算效率提出了更高的要求。因此，需要一种燃烧室试验台虚拟试验建模方法，以解决现有技术中存在的上述技术问题。

发明内容

本发明提供一种燃烧室试验台虚拟试验建模方法。采用这种燃烧室试验台虚拟试验建模方法可以较为逼真地模拟发动机工作的真实情况，提高燃烧室试验台虚拟试验的准确度。

本发明采用的技术方案是：

一种燃烧室试验台虚拟试验建模方法，其包括以下步骤：步骤1：采用有限体积法将试验模型划分成多个有限控制体，建立所述有限控制体的数学模型方程；步骤2：根据所述数学模型方程建立仿真模型，以进行并行运算，获得所述试验模型的相关参数。

优选地，所述试验模型为管道时，所述步骤1包括：对于粘性流动的流体，面积力为法向应力与平行于作用面的切应力的合力：σ_ij＝-pδ_ij+τ_ij，根据高斯公式B表示一个标量场、向量场或张量场，将适用于连续介质的有限控制体的欧拉Euler型积分形式的流动守恒方程组1：

对应地化为方程组2：

式中，单位体积总能量e为单位质量的内能，Q为辐射或化学能释放等因素产生的控制体内单位体积流体热量的增量，为热通量向量，“-”号表示热的流通与外法线方向相反，为作用面外法线方向，为单位质量力，ρf_iV表示体积力，σ_ij为应力张量，其中，i表示作用面外法线方向，j表示面积力方向，n_jσ_jiS表示面积力。

优选地，对于Q＝0的考虑管壁弹性的一维变截面管内流动，只有沿管道轴向的流体运动和管壁弹性变形引起的流体沿管壁外法线方向的膨胀或压缩，从管流中取任一有限控制体积，端面的轴向坐标分别为x₁、x₂，有限控制体的质心处壁面的径向坐标为r_s，当轴向的流动速度和管壁弹性变形引起的流体沿管壁外法线方向的膨胀速度，管壁摩擦取唯象的准稳态描述方式，则方程组2化为如下方程组3：

${\∫\∫\∫}_{V\left(t\right)}{\mathrm{\ρf}}_i\mathrm{dV}={\∫}_{x_1}^{x_2}{\mathrm{\ρf}}_x\mathrm{Adx},f_x=g\cos \mathrm{\θ};$

式中，g为重力加速度，θ为一维管流流动方向至重力加速度方向的旋转角，F_f为管壁摩擦力，f_λ为无量纲摩擦损失系数，D为管道内径，u_s为管壁弹性变形引起的流体沿管壁外法向的膨胀速度，S_s为管壁弹性变形区域的内表面积，λ为流体导热系数，为有限控制体对管壁的单位面积上的热流密度，α为流体与管壁之间的对流换热系数，T_w为管内壁温度，S为有限控制体所在位置管壁的内表面积。

优选地，定义那么f_R相当于单位质量流体受到的管壁摩擦阻力，f_R作用在有限控制体体积上，则：

忽略粘性应力所作的耗散功，即Q＝0的可压缩变截面管流一维瞬变流动守恒方程可积分为方程组4：

$\left\{\begin{array}{c}\frac{\∂}{\∂t}{\∫}_{x_1}^{x_2}\mathrm{\ρAdx}={\left(\mathrm{\ρuA}\right)}_{x=x_1}-{\left(\mathrm{\ρuA}\right)}_{x=x_2}\\ \frac{\∂}{\∂t}{\∫}_{x_1}^{x_2}\mathrm{\ρuAdx}={({\mathrm{\ρu}}^{2}A+\mathrm{pA})}_{x=x_1}-{({\mathrm{\ρu}}^{2}A+\mathrm{pA})}_{x=x_2}+\\ {\∫}_{x_1}^{x_2}p\frac{\∂A}{\∂x}\mathrm{dx}+{\∫}_{x_1}^{x_2}{\mathrm{\ρf}}_x\mathrm{Adx}-{\∫}_{x_1}^{x_2}{f}_R\mathrm{\ρAdx},f_x=g\cos \mathrm{\θ},f_R=\frac{f_{\mathrm{\λ}}u\left|u\right|}{2D}\\ \frac{\∂}{\∂t}{\∫}_{x_1}^{x_2}\mathrm{EAdx}={(\mathrm{EuA}+\mathrm{puA}-\mathrm{\λA}\frac{\∂T}{\∂x})}_{x=x_1}-{(\mathrm{EuA}+\mathrm{puA}-\mathrm{\λA}\frac{\∂T}{\∂x})}_{x=x_2}-{\left({\mathrm{pu}}_s{S}_s\right)}_{r=r_s}+\\ {\∫}_{x_1}^{x_2}{\mathrm{\ρf}}_x\mathrm{uAdx}-\stackrel{.}{q}S\end{array}\right..$

优选地，对每个有限控制体运用积分形式的守恒方程，方程组4最终可离散为数学模型方程：

$\left\{\begin{array}{c}\frac{d\left({\mathrm{\ρ}}_i{V}_i\right)}{\mathrm{dt}}={\mathrm{\ρ}}_i^{\mathrm{in}}{u}_i{A}_i-{\mathrm{\ρ}}_i^{\mathrm{out}}{u}_{i+1}{A}_{i+1}\\ \frac{d\left(E_i{V}_i\right)}{\mathrm{dt}}=(E_i^{\mathrm{in}}{u}_i+p_i^{\mathrm{in}}{u}_i-{\mathrm{\λ}}_i^{\mathrm{in}}\frac{T_i-T_{i-1}}{{\stackrel{\‾}{x}}_i-{\stackrel{\‾}{x}}_{i-1}})A_i-(E_o^{\mathrm{out}}{u}_{i+1}{p}_i^{\mathrm{out}}{u}_{i+1}-{\mathrm{\λ}}_i^{\mathrm{out}}\frac{T_{i+1}-T_i}{{\stackrel{\‾}{x}}_{i+1}-{\stackrel{\‾}{x}}_i})A_{i+1}-\\ p_i{\left(u_s{S}_s\right)}_i+{\mathrm{\ρ}}_i{f}_{x_i}{V}_i\frac{u_i+u_{i+1}}{2}-{\stackrel{.}{q}}_i{S}_i\\ \frac{{\mathrm{dW}}_j}{\mathrm{dt}}=[{\mathrm{\ρ}}_{j-1}{\left(u_j^{\mathrm{in}}\right)}^2+p_{j-1}]{\stackrel{\‾}{A}}_{j-1}-[{\mathrm{\ρ}}_j{\left(u_j^{\mathrm{out}}\right)}^2+p_j]{\stackrel{\‾}{A}}_j+p_{j-1}(A_j-{\stackrel{\‾}{A}}_{j-1})+p_j({\stackrel{\‾}{A}}_j-A_j)+\\ \frac{{\mathrm{\ρ}}_{j-1}{f}_{x_{j-1}}{V}_{j-1}+{\mathrm{\ρ}}_j{f}_{x_j}{V}_j}{2}-f_{R_j}\frac{{\mathrm{\ρ}}_{j-1}{V}_{j-1}+{\mathrm{\ρ}}_j{V}_j}{2},f_{x_j}=g\cos {\mathrm{\θ}}_j,f_{R_j}=\frac{f_{{\mathrm{\λ}}_j}{u}_j\left|u_j\right|}{{2D}_j}\end{array}\right.$

i＝0,1,…,n-1，u为流速，pⁱⁿ为管道进口的压力，p^out为管道出口的压力，T为管道温，A为管道截面，ρ为流体密度，V为有限控制体体积，对于能量方程中的管壁弹性变形项，其中，可以认为

优选地，所述相关参数包括出口流量、流体压力和管道温度，在步骤2中，顺次求解所述数学模型方程中的三个方程，以获得所述相关参数。

优选地，所述并行计算步骤中采用OpenMP多线程架构。

采用上述技术方案，本发明至少具有下列优点：采用本发明的燃烧室试验台虚拟试验建模方法可以更真实地模拟燃烧过程，提高建模质量和虚拟试验的效率。

附图说明

图1为本发明第1实施例的燃烧室试验台虚拟试验建模方法的一维可压缩管流的有限控制体网格示意图。

具体实施方式

为更进一步阐述本发明为达成预定目的所采取的技术手段及功效，以下结合附图及较佳实施例，对本发明进行详细说明如后。

本发明提供的燃烧室试验台虚拟试验建模方法可以充分考虑实际燃烧过程，提高燃烧室试验台虚拟试验的准确性，并且虚拟试验的运行效率得到提高。

该燃烧室试验台虚拟试验建模方法包括步骤1：采用有限体积法将试验模型划分成多个有限控制体，建立所述有限控制体的数学模型方程和步骤2：根据所述数学模型方程建立仿真模型，以进行并行运算，获得所述试验模型的相关参数。

其中，步骤1：采用有限体积法将试验模型划分成多个有限控制体，建立有限控制体的数学模型方程。

适用于连续介质有限控制体的欧拉Euler型积分形式的流动守恒方程组1为：

对于粘性流动，面积力为法向应力与平行于作用面的切应力的合力：σ_ij＝-pδ_ij+τ_ij，根据高斯公式其中B表示一个标量场、向量场或张量场，则上述公式化为方程组2：

式中，单位体积总能量e为单位质量的内能，Q为辐射或化学能释放等因素产生的有限控制体内单位体积流体热量的增量，为热通量向量，其中，“-”号表示热的流通与外法线方向相反，为作用面外法线方向，为单位质量力，ρf_iV表示体积力，σ_ij为应力张量，i表示作用面外法线方向，j表示面积力方向，n_jσ_jiS表示面积力。

对于Q＝0时，考虑管壁弹性的一维变截面管内流动，只有沿管道轴向(x方向)的流体运动和管壁弹性变形引起的流体沿管壁外法线方向的膨胀或压缩，从管流中取任一有限控制体，端面的轴向坐标分别为x₁、x₂，有限控制体的质心处的壁面的径向坐标为r_s。

由于只有轴向的流动速度和管壁弹性变形引起的流体沿管壁外法线方向的膨胀速度，如果管壁摩擦取唯象的准稳态描述方式，方程组2化为方程组3：

${\∫\∫\∫}_{V\left(t\right)}{\mathrm{\ρf}}_i\mathrm{dV}={\∫}_{x_1}^{x_2}{\mathrm{\ρf}}_x\mathrm{Adx},f_x=g\cos \mathrm{\θ};$

由和

得出

式中，g为重力加速度，θ为一维管流流动方向至重力加速度方向的旋转角，F_f为管壁摩擦力，f_λ为无量纲摩擦损失系数(亦称Darcy摩擦因子)，D为管道内径(通常内径应该用d表示，本文为了区别于全微分符号，采用D表示)，u_s为管壁弹性变形引起的流体沿管壁外法向的膨胀速度，S_s为管壁弹性变形区域的内表面积，λ为流体导热系数，为有限控制体对管壁的单位面积上的热流密度(通常可以按对流换热考虑)，α为流体与管壁之间的对流换热系数，T_w为管内壁温度，S为有限控制体所在位置的管壁的内表面积。

定义那么f_R相当于单位质量流体受到的管壁摩擦阻力，这样一来就好像有一个体积力场f_R作用在有限控制体体积上，则：

考虑到实际流体在管壁处的无滑移条件，又由于准一维流动不考虑剪切变形，因此可忽略粘性应力所作的耗散功，即

方程组3，即Q＝0的可压缩变截面管流一维瞬变流动守恒方程最终可积分为如下形式：

$\left\{\begin{array}{c}\frac{\∂}{\∂t}{\∫}_{x_1}^{x_2}\mathrm{\ρAdx}+{\left(\mathrm{\ρuA}\right)}_{x=x_2}-{\left(\mathrm{\ρuA}\right)}_{x=x_1}=0\\ \frac{\∂}{\∂t}{\∫}_{x_1}^{x_2}\mathrm{\ρuAdx}+{({\mathrm{\ρu}}^2+A)}_{x=x_2}-{\left({\mathrm{\ρu}}^{2}A\right)}_{x=x_1}={\∫}_{x_1}^{x_2}\mathrm{\ρ}{f}_x\mathrm{Adx}+{\left(\mathrm{pA}\right)}_{x=x_1}-{\left(\mathrm{pA}\right)}_{x=x_2}+\\ {\∫}_{x_1}^{x_2}p\frac{\∂A}{\∂x}\mathrm{dx}-{\∫}_{x_1}^{x_2}{\mathrm{\ρf}}_x\mathrm{Adx}\\ \frac{\∂}{\∂t}{\∫}_{x_1}^{x_2}\mathrm{EAdx}+{\left(\mathrm{EuA}\right)}_{x=x_2}-{\left(\mathrm{EuA}\right)}_{x=x_1}={\∫}_{x_1}^{x_2}{\mathrm{\ρf}}_x\mathrm{uAdx}+{\left(\mathrm{puA}\right)}_{x=x_1}-{\left(\mathrm{puA}\right)}_{x=x_2}-{\left({\mathrm{pu}}_s{S}_s\right)}_{r=r_s}\\ -\stackrel{.}{q}S-{(-\mathrm{\λA}\frac{\∂T}{\∂x})}_{x=x_2}+{(-\mathrm{\λA}\frac{\∂T}{\∂x})}_{x=x_1}\end{array}\right.,$

即方程组4：

$\left\{\begin{array}{c}\frac{\∂}{\∂t}{\∫}_{x_1}^{x_2}\mathrm{\ρAdx}={\left(\mathrm{\ρuA}\right)}_{x=x_1}-{\left(\mathrm{\ρuA}\right)}_{x=x_2}\\ \frac{\∂}{\∂t}{\∫}_{x_1}^{x_2}\mathrm{\ρuAdx}={({\mathrm{\ρu}}^{2}A+\mathrm{pA})}_{x=x_1}-{({\mathrm{\ρu}}^{2}A+\mathrm{pA})}_{x=x_2}+\\ {\∫}_{x_1}^{x_2}p\frac{\∂A}{\∂x}\mathrm{dx}+{\∫}_{x_1}^{x_2}{\mathrm{\ρf}}_x\mathrm{Adx}-{\∫}_{x_1}^{x_2}{f}_R\mathrm{\ρAdx},f_x=g\cos \mathrm{\θ},f_R=\frac{f_{\mathrm{\λ}}u\left|u\right|}{2D}\\ \frac{\∂}{\∂t}{\∫}_{x_1}^{x_2}\mathrm{EAdx}={(\mathrm{EuA}+\mathrm{puA}-\mathrm{\λA}\frac{\∂T}{\∂x})}_{x=x_1}-{(\mathrm{EuA}+\mathrm{puA}-\mathrm{\λA}\frac{\∂T}{\∂x})}_{x=x_2}-{\left({\mathrm{pu}}_s{S}_s\right)}_{r=r_s}+\\ {\∫}_{x_1}^{x_2}{\mathrm{\ρf}}_x\mathrm{uAdx}-\stackrel{.}{q}S\end{array}\right.$

设管道周长为C，管壁弹性变形区域的周长为C_s，则方程组4还可化为：

即 $\left\{\begin{array}{c}{\∫}_{x_1}^{x_2}(\frac{\∂\mathrm{\ρA}}{\∂t}+\frac{\∂\mathrm{\ρuA}}{\∂x})\mathrm{dx}=0\\ {\∫}_{x_1}^{x_2}[\frac{\∂\mathrm{\ρuA}}{\∂t}+\frac{\∂({\mathrm{\ρu}}^{2}A+\mathrm{pA})}{\∂x}-p\frac{\∂A}{\∂x}-f_x\mathrm{\ρA}+f_R\mathrm{\ρA}]\mathrm{dx}=0\\ {\∫}_{x_1}^{x_2}[\frac{\∂\mathrm{EA}}{\∂t}+\frac{\∂(\mathrm{EuA}+\mathrm{puA}-\mathrm{\λA}\frac{\∂T}{\∂})}{\∂x}+{\mathrm{pu}}_s{C}_s-{\mathrm{\ρf}}_x\mathrm{uA}+\stackrel{.}{q}C]\mathrm{dx}=0\end{array}\right.:$

由于x₁、x₂是任意选取的，且假定被积函数连续，则可得：

$\left\{\begin{array}{c}\frac{\∂\left(\mathrm{\ρA}\right)}{\∂t}+\frac{\∂\left(\mathrm{\ρuA}\right)}{\∂x}=0\\ \frac{\∂\left(\mathrm{\ρuA}\right)}{\∂t}+\frac{\∂({\mathrm{\ρu}}^{2}A+\mathrm{pA})}{\∂x}=p\frac{\∂A}{\∂x}+{\mathrm{\ρAf}}_x-{\mathrm{\ρAf}}_R\\ \frac{\∂\left(\mathrm{EA}\right)}{\∂t}+\frac{\∂(\mathrm{EuA}+\mathrm{puA}-\mathrm{\λA}\frac{\∂T}{\∂x})}{\∂x}=-{\mathrm{pu}}_s{C}_s+{\mathrm{\ρf}}_x\mathrm{uA}-\stackrel{.}{q}C,f_x=g\cos \mathrm{\θ},f_R=\frac{f_{\mathrm{\λ}}u\left|u\right|}{2D}\end{array}\right.$

上列方程组的矢量形式为：

此即为Q＝0的可压缩变截面管流一维瞬变流动守恒方程的微分形式，其中，

单位体积总能量有限控制体对管壁的单位面积上的热流密度通常可以按对流换热考虑，即

补充流体的状态方程：p＝p(ρ,T)，e＝e(p,ρ)

由上述方程出发可以分析考虑管壁弹性变形、考虑重力场、考虑摩擦力影响、考虑流体轴向热传导、考虑管壁传热的一维可压缩变截面管流的瞬变流动，其中能量方程中的管壁弹性变形项在有限控制体的体积可变时(例如某部分管壁为可动活塞的流体腔)的物理意义为流体向外输出的膨胀功。

以下以试验模型为管道进行详细地说明。

从一维可压缩瞬变流积分形式的守恒方程出发，通过引入空间位置交错的两种有限控制体：状态单元和速度单元，对其分别进行质量、能量和动量守恒方程的空间离散，提出了一维可压缩瞬变流的有限元状态变量模型，解决了以往同类模型的不足之处。与计算流体力学中传统的加权残差解有限元方法不同，该模型只是在几何分割上采用了有限元的思想，最终的常微分方程是通过有限体积离散得到的，实质上属于有限体积法。有限元状态变量模型忽略了重力场和管流轴向热传导，最终得到的常微分方程是针对恒定体积有限控制体推导的，未考虑有限控制体体积可变时的情况。

本节从上一节推导得到的Q＝0的可压缩变截面管流积分形式的一维瞬变流动守恒方程出发，采用有限元状态变量模型的空间交错网格划分方法，推导适用范围更加广泛的有限体积模型。

如图1所示，将流体管道进行网格划分，引入空间位置交错的两种有限控制体单元，在状态单元(图中的实线单元)内，认为流体的压强、密度、温度和单位质量内能等状态参数是瞬时一致的；在速度单元(图中的虚线单元)内，认为流体的速度是瞬时一致的。为了区别这两种有限控制体，分别采用i、j表示状态单元与速度单元的编号，管道两端的半个速度单元贮存管道的边界条件，分别用下标a和b表示。单元长度可以是相等的或不等的，管道可以是等截面或变截面的。采用这样的单元分割方法，单元之间的边界条件易于确定，适应性强，能处理任意的边界条件。

对每个有限控制体积单元运用积分形式的守恒方程，由于在单元内各变量是瞬时一致的，与位置无关，方程组4最终可离散为如下的数学模型方程：

$\left\{\begin{array}{c}\frac{d\left({\mathrm{\ρ}}_i{V}_i\right)}{\mathrm{dt}}={\mathrm{\ρ}}_i^{\mathrm{in}}{u}_i{A}_i-{\mathrm{\ρ}}_i^{\mathrm{out}}{u}_{i+1}{A}_{i+1}\\ \frac{d\left(E_i{V}_i\right)}{\mathrm{dt}}=(E_i^{\mathrm{in}}{u}_i+p_i^{\mathrm{in}}{u}_i-{\mathrm{\λ}}_i^{\mathrm{in}}\frac{T_i-T_{i-1}}{{\stackrel{\‾}{x}}_i-{\stackrel{\‾}{x}}_{i-1}})A_i-(E_o^{\mathrm{out}}{u}_{i+1}{p}_i^{\mathrm{out}}{u}_{i+1}-{\mathrm{\λ}}_i^{\mathrm{out}}\frac{T_{i+1}-T_i}{{\stackrel{\‾}{x}}_{i+1}-{\stackrel{\‾}{x}}_i})A_{i+1}-\\ p_i{\left(u_s{S}_s\right)}_i+{\mathrm{\ρ}}_i{f}_{x_i}{V}_i\frac{u_i+u_{i+1}}{2}-{\stackrel{.}{q}}_i{S}_i\\ \frac{{\mathrm{dW}}_j}{\mathrm{dt}}=[{\mathrm{\ρ}}_{j-1}{\left(u_j^{\mathrm{in}}\right)}^2+p_{j-1}]{\stackrel{\‾}{A}}_{j-1}-[{\mathrm{\ρ}}_j{\left(u_j^{\mathrm{out}}\right)}^2+p_j]{\stackrel{\‾}{A}}_j+p_{j-1}(A_j-{\stackrel{\‾}{A}}_{j-1})+p_j({\stackrel{\‾}{A}}_j-A_j)+\\ \frac{{\mathrm{\ρ}}_{j-1}{f}_{x_{j-1}}{V}_{j-1}+{\mathrm{\ρ}}_j{f}_{x_j}{V}_j}{2}-f_{R_j}\frac{{\mathrm{\ρ}}_{j-1}{V}_{j-1}+{\mathrm{\ρ}}_j{V}_j}{2},f_{x_j}=g\cos {\mathrm{\θ}}_j,f_{R_j}=\frac{f_{{\mathrm{\λ}}_j}{u}_j\left|u_j\right|}{{2D}_j}\end{array}\right.,$ 其中i＝0,1,…,n-1。u为流速，pⁱⁿ为管道进口的压力，p^out为管道出口的压力，T为管道温，A为管道截面，ρ为流体密度，V为有限控制体体积。对于能量方程中的管壁弹性变形项(膨胀功)，通常可以认为

接下来，在步骤2中的相关参数包括出口流量、流体压力和管道温度，在步骤2中，从初始有限控制体开始顺次求解上述数学模型方程中的三个方程，以获得相关参数。进一步地，步骤2中的并行计算采用OpenMP多线程架构。

1)对状态单元：i＝0,1,…,n-1

连续方程： $\frac{d\left({\mathrm{\ρ}}_i{V}_i\right)}{\mathrm{dt}}={\mathrm{\ρ}}_i^{\mathrm{in}}{u}_i{A}_i-{\mathrm{\ρ}}_i^{\mathrm{out}}{u}_{i+1}{A}_{i+1}$

能量方程： $\begin{array}{c}\frac{d\left(E_i{V}_i\right)}{\mathrm{dt}}=(E_i^{\mathrm{in}}{u}_i+p_i^{\mathrm{in}}{u}_i-{\mathrm{\λ}}_i^{\mathrm{in}}\frac{T_i-T_{i-1}}{{\stackrel{\‾}{x}}_i-{\stackrel{\‾}{x}}_{i-1}})A_i-\\ (E_i^{\mathrm{out}}{u}_{i+1}+p_i^{\mathrm{out}}-{\mathrm{\λ}}_i^{\mathrm{out}}\frac{T_{i+1}-T_i}{{\stackrel{\‾}{x}}_{i+1}-{\stackrel{\‾}{x}}_i})A_{i+1}-p_i\frac{{\mathrm{dV}}_i}{\mathrm{dt}}+{\mathrm{\ρ}}_i{f}_{x_i}{V}_i\frac{u_i+u_{i+1}}{2}-{\stackrel{.}{q}}_i{S}_i\end{array};$

式中，单位体积总能量 $E_i={\mathrm{\ρ}}_i{e}_i+\frac{1}{2}{\mathrm{\ρ}}_i{\left(\frac{u_i+u_{i+1}}{2}\right)}^2$

状态单元边界参数采用迎风格式： ${\mathrm{Var}}_i^{\mathrm{in}}=\left\{\begin{array}{ccc}{\mathrm{Var}}_{i-1}& \mathrm{if}& u_i\&GreaterEqual;0\\ {\mathrm{Var}}_i& \mathrm{if}& u_i<0\end{array}\right.,$ ${\mathrm{Var}}_i^{\mathrm{out}}=\left\{\begin{array}{ccc}{\mathrm{Var}}_i& \mathrm{if}& u_{i+1}\&GreaterEqual;0\\ {\mathrm{Var}}_{i+1}& \mathrm{if}& u_{i+1}<0\end{array}\right.,$ Var∈{ρ,p,E,λ,e}，Var_-1＝Var_a，Var_n＝Var_b

对液体和流速不大(Ma<0.3)的气体速度单元，可以忽略动能的影响，即 $E_i={\mathrm{\&rho;}}_i{e}_i+\frac{1}{2}{\mathrm{\&rho;}}_i{\left(\frac{u_i+u_{i+1}}{2}\right)}^2\&ap;{\mathrm{\&rho;}}_i{e}_i,$ 则能量方程简化为方程5：

$\begin{array}{c}\frac{d\left({\mathrm{\&rho;}}_i{e}_i{V}_i\right)}{\mathrm{dt}}=(E_i^{\mathrm{in}}{u}_i+p_i^{\mathrm{in}}{u}_i-{\mathrm{\&lambda;}}_i^{\mathrm{in}}\frac{T_i-T_{i-1}}{{\stackrel{\&OverBar;}{x}}_i-{\stackrel{\&OverBar;}{x}}_{i-1}})A_i-(E_i^{\mathrm{out}}{u}_{i+1}+p_i^{\mathrm{out}}{u}_{i+1}-{\mathrm{\&lambda;}}_i^{\mathrm{out}}\frac{T_{i+1}-T_i}{{\stackrel{\&OverBar;}{x}}_{i+1}-{\stackrel{\&OverBar;}{x}}_i})A_{i+1}-\\ p_i\frac{{\mathrm{dV}}_i}{\mathrm{dt}}+{\mathrm{\&rho;}}_i{f}_{x_i}{V}_i\frac{u_i+u_{i+1}}{2}-{\stackrel{.}{q}}_i{S}_i\end{array}$

a)对于理想气体，引入理想气体的内能状态方程

如果假定气体的比热比γ为常数，方程5化为：

$\begin{array}{c}\frac{1}{\mathrm{\&gamma;}-1}\frac{d\left(p_i{V}_i\right)}{\mathrm{dt}}=[(\frac{p_i^{\mathrm{in}}}{\mathrm{\&gamma;}-1}+\frac{1}{2}{\mathrm{\&rho;}}_i^{\mathrm{in}}{u}_i^2)u_i+p_i^{\mathrm{in}}{u}_i-{\mathrm{\&lambda;}}_i^{\mathrm{in}}\frac{T_i-T_{i-1}}{{\stackrel{\&OverBar;}{x}}_i-{\stackrel{\&OverBar;}{x}}_{i-1}}]A_i-\\ [(\frac{p_i^{\mathrm{out}}}{\mathrm{\&gamma;}-1}+\frac{1}{2}{\mathrm{\&rho;}}_i^{\mathrm{out}}{u}_{i+1}^u)u_{i+1}+p_i^{\mathrm{out}}{u}_{i+1}-{\mathrm{\&lambda;}}_i^{\mathrm{out}}\frac{T_{i+1}-T_i}{{\stackrel{\&OverBar;}{x}}_{i+1}-{\stackrel{\&OverBar;}{x}}_i}]A_{i+1}-p_i\frac{{\mathrm{dV}}_i}{\mathrm{dt}}+{\mathrm{\&rho;}}_i{f}_{x_i}{V}_i\frac{u_i+u_{i+1}}{2}-{\stackrel{\&CenterDot;}{q}}_i{S}_i\end{array}$

即 $\begin{array}{c}\frac{{\mathrm{dp}}_i}{\mathrm{dt}}=\frac{1}{V_i}\{\mathrm{\&gamma;}(p_i^{\mathrm{in}}{u}_i{A}_i-p_i^{\mathrm{out}}{u}_{i+1}{A}_{i+1})+(\mathrm{\&gamma;}-1)[\frac{1}{2}{\mathrm{\&rho;}}_i^{\mathrm{in}}{u}_i^3{A}_i-\frac{1}{2}{\mathrm{\&rho;}}_i^{\mathrm{out}}{u}_{i+1}^3{A}_{i+1}+\\ {\mathrm{\&lambda;}}_i^{\mathrm{in}}\frac{T_{i-1}-T_i}{{\stackrel{\&OverBar;}{x}}_i-{\stackrel{\&OverBar;}{x}}_{i-1}}{A}_i-A_i^{\mathrm{out}}\frac{T_i-T_{i+1}}{{\stackrel{\&OverBar;}{x}}_{i+1}-{\stackrel{\&OverBar;}{x}}_i}{A}_{i+1}+{\mathrm{\&rho;}}_i{f}_{x_i}{V}_i\frac{u_i+u_{i+1}}{2}-{\stackrel{.}{q}}_i{S}_i]-{\mathrm{\&gamma;p}}_i\frac{{\mathrm{dV}}_i}{\mathrm{dt}}\}\end{array};$

b)对不可压缩流体，可采用文阻尼修正后的能量方程：

$\frac{{\mathrm{dp}}_i}{\mathrm{dt}}=\frac{1}{C_i}(Q_i-Q_{i+1})+\mathrm{\&epsiv;}\frac{\mathrm{\&rho;a}}{{\stackrel{\&OverBar;}{A}}_i}(\frac{{\mathrm{dQ}}_i}{\mathrm{dt}}-\frac{d{Q}_{i+1}}{\mathrm{dt}}),C_i=\frac{{\mathrm{\&xi;V}}_i}{{\mathrm{\&rho;}}_i{a}_i^2},Q_i=u_i{A}_i;$

式中，ξ为小分段数时并联导纳的修正系数，ε为阻尼修正系数，音速a可取为常数。

$\mathrm{\&xi;}=\left\{\begin{array}{cc}0.8016& n=1\\ 0.9496& n=2\\ 1& n\&GreaterEqual;3\end{array}\right.,\mathrm{\&epsiv;}=\left\{\begin{array}{cc}0.05& n=\mathrm{1,2}\\ 0.10& 3\&le;n\&le;10\\ 0.5& n>10\end{array}\right.$

以上各式中，Ma为马赫数，e为单位质量的内能，V_i为对应状态单元的体积，γ为气体比热比，Var_a、Var_b为与流体管道相连元件边界单元的状态参数。

2)对速度单元：j＝0,1,…,n

动量方程：

$\left\{\begin{array}{c}\frac{{\mathrm{dW}}_j}{\mathrm{dt}}=[{\mathrm{\&rho;}}_{j-1}{\left(u_j^{\mathrm{in}}\right)}^2+p_{j-1}]{\stackrel{\&OverBar;}{A}}_{j-1}-[{\mathrm{\&rho;}}_j{\left(u_j^{\mathrm{out}}\right)}^2+p_j]{\stackrel{\&OverBar;}{A}}_j+p_{j-1}(A_j-{\stackrel{\&OverBar;}{A}}_{j-1})+p_j({\stackrel{\&OverBar;}{A}}_j-A_j)+\\ \frac{{\mathrm{\&rho;}}_{j-1}{f}_{x_{j-1}}{V}_{j-1}+{\mathrm{\&rho;}}_j{f}_{x_j}{V}_j}{2}-f_{R_j}\frac{{\mathrm{\&rho;}}_{j-1}{V}_{j-1}+{\mathrm{\&rho;}}_j{V}_j}{2},f_{x_j}=g\cos {\mathrm{\&theta;}}_j,f_{R_j}=\frac{f_{{\mathrm{\&lambda;}}_j}{u}_j\left|u_j\right|}{{2D}_j}\end{array}\right.;$

式中，速度单元动量 $W_j=\frac{{\mathrm{\&rho;}}_{j-1}{V}_{j-1}+{\mathrm{\&rho;}}_j{V}_j}{2}{u}_j$

速度单元边界参数采用迎风格式： $u_j^{\mathrm{in}}=\left\{\begin{array}{ccc}{u}_{j-1}& \mathrm{if}& u_j\&GreaterEqual;0\\ u_j& \mathrm{if}& u_j<0\end{array}\right.,$ $u_j^{\mathrm{out}}=\left\{\begin{array}{ccc}{u}_j& \mathrm{if}& u_{j+1}\&GreaterEqual;0\\ u_{j+1}& \mathrm{if}& u_{j+1}<0\end{array}\right.$

也可以采用算术平均： $u_j^{\mathrm{in}}=\frac{1}{2}(u_{j-1}+u_j),$ u_-1＝u₀； $u_j^{\mathrm{out}}=\frac{1}{2}(u_j+u_{j+1}),$ u_n+1＝u_n

V_j为对应状态单元的体积，V_-1＝ξ₀V₀，V_n＝ξ_nV_n-1，ξ₀和ξ_n根据管道的物理边界条件在0～1之间取值，如对管道闭端可取ξ₀＝0。当需要捕捉锐利的激波边缘时，边界速度采用迎风格式比采用算术平均值精度要高一些。

3)状态单元压强计算公式

对理想气体： $p_i=(\mathrm{\&gamma;}-1)[E_i-\frac{1}{2}{\mathrm{\&rho;}}_i{\left(\frac{u_i+u_{i+1}}{2}\right)}^2];$

4)状态单元温度计算公式：

对理想气体：式中，气体常数M_g为分子摩尔质量

5)速度单元速度计算公式：

6)速度单元无量纲摩擦损失系数(Darcy摩擦因子)计算公式：

圆管内的流动通常可分为五个区：第Ⅰ区，层流区，0≤Re≤2300；第Ⅱ区，层流转向紊流的过渡区，2300<Re≤4000；第Ⅲ区，水力光滑区，第Ⅳ区，水力光滑向水力粗糙的过渡区，第Ⅴ区，水力粗糙区(亦称阻力平方区或完全粗糙区)，式中，h_Δ为管壁凸起高度，即绝对粗糙度，Re为雷诺数。

在上述五个区域分别采用布拉休斯(Blasius)公式、布拉休斯公式或普朗特阻力公式、尼古拉兹(Nikuradse)公式或科尔布鲁克(Colebrook)公式、卡门－尼古拉兹公式或适用于整个紊流的通用经验公式。

很多情况下***初始流场速度为零，即Re(0)＝0，认为0≤Re≤10时，f_λ＝6.4。

速度单元雷诺数可按下式计算：μ为流体的动力粘度

7)状态单元对管壁的单位面积上的热流密度计算公式：

按对流换热考虑：

式中，α为管内流体与管壁之间的对流换热系数，T_w为管内壁温度。

8)对流换热系数计算公式：

按雷诺数将圆管内的流动分为五个区域：

a)可以按水平或竖直管内层流混合对流计算，对于贮箱等大容器，由于缺乏专门的有限空间自然对流换热公式，暂时按大空间自然对流换热计算，根据管道布置情况、Pr、Ra值选择合适的计算公式，例如Brown-Gauvin层流混合对流公式、Hyman-Bonila-Ehrlich公式、McAdamas公式等；

b)按层流换热计算，从Hausene公式、Seider-Tate层流换热公式、米海耶夫公式中选择一种；

c)2300<Re≤10⁴，按过渡区换热计算，根据情况从管内过渡区换热公式中选择一种；

d)10⁴<Re≤1.25×10⁵，1.25×10⁵<Re≤1.75×10⁵，按紊流换热计算，根据Re值、流体相态、温差大小、冷却或加热情况选择合适的计算公式，例如Dittus-Boelter公式、Seider-Tate紊流换热公式等；

e)Re>1.75×10⁵，属于超高雷诺数紊流换热，暂时采用Dittus-Boelter公式。

式中，普朗特数格拉晓夫数瑞利数运动粘性系数u＝m/r，μ、c_P、β分别为流体的动力粘度、定压比热、容积膨胀系数，g为重力加速度，特征尺度l对竖圆柱取高度、对横圆柱取直径。

9)动态仿真全场时间步长的确定

根据一维流场CFL条件确定的状态单元局部时间步长为：

${\mathrm{\&Delta;t}}_i=\frac{\mathrm{CFL}}{\frac{|(u_i+u_{i+1})/2|}{x_{i+1}-x_i}+\frac{a_i}{x_{i+1}-x_i}},$ 其中0<CFL≤1，可取CFL＝0.8

动态仿真时必须保证全场(包括温度场)在时间推进上的同步性，因此全场时间步长取为：Δt＝min{min{Δt_i},min{Δτ_i,k}}；

式中，Δτ_i,k为下一节将会讨论到的根据管壁温度场网格数值计算稳定性准则确定的时间步长，a为状态单元处的音速。

10)状态单元音速的计算公式：

对理想气体： $a_i=\sqrt{{\mathrm{\&gamma;RT}}_i}=\sqrt{\frac{{\mathrm{\&gamma;p}}_i}{{\mathrm{\&rho;}}_i}}$

对不可压缩流体：为液体体积弹性模量，取为常数。

OpenMP是一种面向共享存储体系结构的多线程并行编程语言，是可被用于显示指导多线程、共享内存并行的API，支持多种平台，包括类UNIX***和WindowsNT***(Windows2000、WindowsXP和WindowsVista等)。OpenMP由编译制导(CompilerDirective)、运行时库函数(RuntmieLibrary)和环境变量(EnvironmentVariables)构成，同时支持C、C++语言或Fortran语言，可选择任意一种语言以及支持OpenMP的编译器编写程序。OpenMP的优点是实现线程级并行，线程间通信通过读/写共享变量实现，避免了消息传递的开销；将串行程序转换为并行程序时无须对代码作大的改动，编程简单；允许运行时调度，提供了细粒度和粗粒度并行机制。

OpenMP使用派生/连接Fork-Join并行执行模型。一个OpenMP程序从一个单个线程开始执行，在程序某点需要并行时程序派生(Fork)出一些额外(可能为零个)线程组成线程组，被派生出来的线程称为组的从属线程，并行区域中的代码在不同的线程中并行执行，程序执行到并行区域末尾，线程将会等待直到整个线程组到达，然后将它们连接Join在一起。在该点处线程组中的从属线程终止而初始主线程继续执行直到下一个并行区域到来。一个程序中可以定义任意数目的并行块，因此，在一个程序的执行中可派生/连接Fork-Join若干次。

这就是标准的并行模式fork/join式并行模式，标准并行模式执行代码的基本思想是，程序开始时只有一个主线程，程序中的串行部分都由主线程执行，并行的部分是通过派生其他线程来执行，但是如果并行部分没有结束时是不会执行串行部分的。

通过具体实施方式的说明，应当可对本发明为达成预定目的所采取的技术手段及功效得以更加深入且具体的了解，然而所附图示仅是提供参考与说明之用，并非用来对本发明加以限制。

Claims (7)

1.一种燃烧室试验台虚拟试验建模方法，其特征在于，包括以下步骤：

步骤1：采用有限体积法将试验模型划分成多个有限控制体，建立所述有限控制体的数学模型方程；

步骤2：根据所述数学模型方程建立仿真模型，以进行并行运算，获得所述试验模型的相关参数。

2.根据权利要求1所述的燃烧室试验台虚拟试验建模方法，其特征在于，所述试验模型为管道时，所述步骤1包括：

对于粘性流动的流体，面积力为法向应力与平行于作用面的切应力的合力：σ_ij＝-pδ_ij+τ_ij，根据高斯公式B表示一个标量场、向量场或张量场，将适用于连续介质的有限控制体的欧拉Euler型积分形式的流动守恒方程组1：

对应地化为方程组2：

式中，单位体积总能量e为单位质量的内能，Q为辐射或化学能释放等因素产生的控制体内单位体积流体热量的增量，为热通量向量，“-”号表示热的流通与外法线方向相反，为作用面外法线方向，为单位质量力，ρf_iV表示体积力，σ_ij为应力张量，其中，i表示作用面外法线方向，j表示面积力方向，n_jσ_jiS表示面积力。

3.根据权利要求2所述的燃烧室试验台虚拟试验建模方法，其特征在于，对于Q＝0的考虑管壁弹性的一维变截面管内流动，只有沿管道轴向的流体运动和管壁弹性变形引起的流体沿管壁外法线方向的膨胀或压缩，从管流中取任一有限控制体积，端面的轴向坐标分别为x₁、x₂，有限控制体的质心处壁面的径向坐标为r_s，当轴向的流动速度和管壁弹性变形引起的流体沿管壁外法线方向的膨胀速度，管壁摩擦取唯象的准稳态描述方式，则方程组2化为如下方程组3：

${\&Integral;\&Integral;\&Integral;}_{V\left(t\right)}\mathrm{\&rho;}{f}_i\mathrm{dV}={\&Integral;}_{x_1}^{x_2}\mathrm{\&rho;}{f}_x\mathrm{Adx},f_x=g\cos \mathrm{\&theta;};$

式中，g为重力加速度，θ为一维管流流动方向至重力加速度方向的旋转角，F_f为管壁摩擦力，f_λ为无量纲摩擦损失系数，D为管道内径，u_s为管壁弹性变形引起的流体沿管壁外法向的膨胀速度，S_s为管壁弹性变形区域的内表面积，λ为流体导热系数，为有限控制体对管壁的单位面积上的热流密度，α为流体与管壁之间的对流换热系数，T_w为管内壁温度，S为有限控制体所在位置管壁的内表面积。

4.根据权利要求3所述的燃烧室试验台虚拟试验建模方法，其特征在于，定义那么f_R相当于单位质量流体受到的管壁摩擦阻力，f_R作用在有限控制体体积上，则：

忽略粘性应力所作的耗散功，即则将方程组3，即Q＝0的可压缩变截面管流一维瞬变流动守恒方程可积分为方程组4：

$\left\{\begin{array}{c}\frac{\&PartialD;}{\&PartialD;t}{\&Integral;}_{x_1}^{x_2}\mathrm{\&rho;Adx}={\left(\mathrm{\&rho;uA}\right)}_{x=x_1}-{\left(\mathrm{\&rho;uA}\right)}_{x=x_2}\\ \frac{\&PartialD;}{\&PartialD;t}{\&Integral;}_{x_1}^{x_2}\mathrm{\&rho;uAdx}={({\mathrm{\&rho;u}}^{2}A+\mathrm{pA})}_{x=x_1}-{({\mathrm{\&rho;u}}^{2}A+\mathrm{pA})}_{x=x_2}+\\ {\&Integral;}_{x_1}^{x_2}p\frac{\&PartialD;A}{\&PartialD;x}\mathrm{dx}+{\&Integral;}_{x_1}^{x_2}\mathrm{\&rho;}{f}_x\mathrm{Adx}-{\&Integral;}_{x_1}^{x_2}{f}_R\mathrm{\&rho;Adx},f_x=g\cos \mathrm{\&theta;},f_R=\frac{f_{\mathrm{\&lambda;}}u\left|u\right|}{2D}\\ \frac{\&PartialD;}{\&PartialD;t}{\&Integral;}_{x_1}^{x_2}\mathrm{EAdx}={(\mathrm{EuA}+\mathrm{puA}-\mathrm{\&lambda;A}\frac{\&PartialD;T}{\&PartialD;x})}_{x=x_1}-{(\mathrm{EuA}+\mathrm{puA}-\mathrm{\&lambda;A}\frac{\&PartialD;T}{\&PartialD;x})}_{x=x_2}-{\left({\mathrm{pu}}_s{S}_s\right)}_{r=r_s}+\\ {\&Integral;}_{x_1}^{x_2}\mathrm{\&rho;}{f}_x\mathrm{uAdx}-\stackrel{.}{q}S\end{array}\right..$

5.根据权利要求4所述的燃烧室试验台虚拟试验建模方法，其特征在于，对每个有限控制体运用积分形式的守恒方程，方程组4最终可离散为数学模型方程：

$\left\{\begin{array}{c}\frac{d\left({\mathrm{\&rho;}}_i{V}_i\right)}{\mathrm{dt}}={\mathrm{\&rho;}}_i^{\mathrm{in}}{u}_i{A}_i-{\mathrm{\&rho;}}_i^{\mathrm{out}}{u}_{i+1}{A}_{i+1}\\ \frac{d\left(E_i{V}_i\right)}{\mathrm{dt}}=(E_i^{\mathrm{in}}{u}_i+p_i^{\mathrm{in}}{u}_i-{\mathrm{\&lambda;}}_i^{\mathrm{in}}\frac{T_i-T_{i-1}}{{\stackrel{\&OverBar;}{x}}_i-{\stackrel{\&OverBar;}{x}}_{i-1}})A_i-(E_i^{\mathrm{out}}{u}_{i+1}+p_i^{\mathrm{out}}{u}_{i+1}-{\mathrm{\&lambda;}}_i^{\mathrm{out}}\frac{T_{i+1}-T_i}{{\stackrel{\&OverBar;}{x}}_{i+1}-{\stackrel{\&OverBar;}{x}}_i})A_{i+1}-\\ p_i\left(u_s{S}_s\right)+{\mathrm{\&rho;}}_i{f}_{x_i}{V}_i\frac{u_i+u_{i+1}}{2}{\stackrel{.}{q}}_i{S}_i\\ \frac{{\mathrm{dW}}_j}{\mathrm{dt}}=[{\mathrm{\&rho;}}_{j-1}{\left(u_j^{\mathrm{in}}\right)}^2+p_{j-1}]{\stackrel{\&OverBar;}{A}}_{j-1}-[{\mathrm{\&rho;}}_j{\left(u_j^{\mathrm{out}}\right)}^2+p_j]{\stackrel{\&OverBar;}{A}}_j+p_{j-1}(A_j-{\stackrel{\&OverBar;}{A}}_{j-1})+p_j({\stackrel{\&OverBar;}{A}}_j-A_j)+\\ \frac{{\mathrm{\&rho;}}_{j-1}{f}_{x_{j-1}}{V}_{j-1}+{\mathrm{\&rho;}}_j{f}_{x_j}{V}_j}{2}-f_{R_j}\frac{{\mathrm{\&rho;}}_{j-1}{V}_{j-1}{+\mathrm{\&rho;}}_j{V}_j}{2},f_{x_j}=g\cos {\mathrm{\&theta;}}_j,f_{R_j}=\frac{f_{{\mathrm{\&lambda;}}_j}{u}_j\left|u_j\right|}{2D_j}\end{array}\right.$

i＝0,1,…,n-1，u为流速，pⁱⁿ为管道进口的压力，p^out为管道出口的压力，T为管道温，A为管道截面，ρ为流体密度，V为有限控制体体积，对于能量方程中的管壁弹性变形项，其中，可以认为

6.根据权利要求5所述的燃烧室试验台虚拟试验建模方法，其特征在于，所述相关参数包括出口流量、流体压力和管道温度，在步骤2中，顺次求解所述数学模型方程中的三个方程，以获得所述相关参数。

7.根据权利要求1至6中任一项所述的燃烧室试验台虚拟试验建模方法，其特征在于，所述并行计算步骤中采用OpenMP多线程架构。