56 工程科学学报，第42卷，增刊1 当地径向速度在x或y方向的偏导数，s;x和y分 式中：Ce、Cc、oK及oe为计算常数，其数值分别 别表示轴向和径向方向：符号“±”代表特征线左侧 为1.44、1.92、1.0及1.341：Gx及G分别代表因 或右侧方向；y为氧气比热容，数值为1.4. 平均速度梯度或耗散率产生的湍流动能源项， 2数值模拟模型建立 kgms3:μ为层流黏度系数，Pas;S为自定义 附加湍流动能源项，kgm1s3:浮力剪切层内主流 2.1控制方程 方向与重力方向一致时，C3e常数为1；浮力剪切层 本文选取雷诺时均方法(Reynolds Averaged 内主流方向垂直于重力方向时，C3c常数为0.YM为 Navier--Stokes,简称RANS)求解流场控制方程，该 可压缩湍流的脉动值，kgms,其计算公式如式 方法将湍流分为平均运动和脉动运动两部分，并用 (12)所示 YM=2ps K 时间平均值和脉动值共同计算流场内变量.雷诺时 yRT (12) 均方法在直角坐标系下，以张量形式表示的可压缩 式中：R为气态常数，JmoK,数值为8.314；T为 流时均Navier--Stokes控制微分方程如下所示l): 气相温度，K 连续性方程 2.2网格划分及过程算法 dp d(pvi) dt+ =0 (6) 图2为超音速氧枪的二维网格结构示意图 0xi 动量方程 考虑到拉瓦尔喷管具有几何对称与流动对称特 点，为降低求解过程计算量，选取喷管半截面及轴 apvi）,pviv）_ap, 0(tij-pv'iv'j) (7) 对称处理进行模拟.计算域均采用结构化网格，为 Ot +0x对j 0xi oxi 提高模拟结果准确性，对喷管近壁面、喉口附近及 能量方程 射流速度梯度较大部分的区域进行网格加密处 a(pE)a(pEvi+pvi) (qa +Cppiv;T) 理.初始状态下，计算域充满静止不动的空气，表1 Ot 0xi Oxi 及表2分别为模型边界条件及相关物化参数. O(tijvj-Pavv'j) (8) (a) 8xi 式中：1为射流时间；及号分别为在x或x方向 的平均速度分量，ms;和分别为：或x的平 0.02mm 均速度分量波动值，ms;E为总能量，J:t为黏性 应力(Nm:p为压强，Pa;Cp为比定压热容 kgK:T为脉动温度，K;qa为物质a的热流， Js.pvy代表雷诺应力，由Boussinesq涡黏性假 设建立雷诺应力与时均运动速度的关系为： 4+m-2o 20mm -p啊=（成+a)pk+4 δij (9) xk (b) 式中：分别为在x:方向的平均速度分量，ms: 山为湍流黏度系数，Pas;K为湍流动能，m2s2; 6为克罗内克符号，当j,6l;当，=0 500mm 选取可实现湍流模型确定模型内流场，其中 图2模型网格结构图.()拉瓦尔喷管网格示意图：(b)计算域网格 湍流动能K及湍流耗散率ε(m2s)计算公式如下 示意图 所示： Fig.2 Schematic of the mesh distribution of the simulation model: (a)Laval nozzle;(b)computational domain ax(Kv)= μ+ GK+Gs-ps-YM 模型中空气与氧气的定压比热容C,(小kgK), (10) 分别按式(13)和(14)计算 pe2-+e） Cp.air=1161.48-2.36×T+1.49×10-2×r2- 5.03×10-5×T3+9.93×10-8×T4-1.11×10-10× pC2eK+ (11) T5+6.54×10-l4×T6-1.57×10-17×T7(13)当地径向速度在 x 或 y 方向的偏导数，s −1 ；x 和 y 分 别表示轴向和径向方向；符号“±”代表特征线左侧 或右侧方向；γ 为氧气比热容，数值为 1.4. 2    数值模拟模型建立 2.1    控制方程 本文选取雷诺时均方法（Reynolds Averaged Navier–Stokes，简称 RANS）求解流场控制方程，该 方法将湍流分为平均运动和脉动运动两部分，并用 时间平均值和脉动值共同计算流场内变量. 雷诺时 均方法在直角坐标系下，以张量形式表示的可压缩 流时均 Navier–Stokes 控制微分方程如下所示[15] ： 连续性方程 ∂ρ ∂t + ∂(ρvi) ∂xi = 0 （6） 动量方程 ∂(ρvi) ∂t + ∂(ρvivj) ∂xj =− ∂p ∂xi + ∂(τi j −ρv ′ iv ′ j) ∂xj （7） 能量方程 ∂(ρE) ∂t + ∂(ρEvi + pvi) ∂xi = − ∂(qa +Cpρiv ′ i T ′ ) ∂xi + ∂(τi jvj −ρav ′ i v ′ j vj) ∂xi （8） v ′ i v ′ j T ′ ρv ′ i v ′ j 式中：t 为射流时间；vi 及 vj 分别为在 xi 或 xj 方向 的平均速度分量，m·s−1 ； 和 分别为 xi 或 xj 的平 均速度分量波动值，m·s−1 ；E 为总能量，J；τij 为黏性 应 力 (N·m−2)； p 为 压 强 ， Pa； Cp 为 比 定 压 热 容 (J·kg−1·K−1)； 为脉动温度，K；qa 为物质 a 的热流， J·s−1 . 代表雷诺应力，由 Boussinesq 涡黏性假 设建立雷诺应力与时均运动速度的关系为： −ρv ′ i v ′ j = µt ( ∂vi ∂xj + ∂vj ∂xi ) − 2 3 ( ρK +µt ∂vk ∂xk ) δi j （9） 式中：vk 分别为在 xk 方向的平均速度分量，m·s−1 ； μt 为湍流黏度系数，Pa·s−1 ；K 为湍流动能，m 2 ·s−2 ； δij 为克罗内克符号，当 i=j，δij=1；当 i≠j，δij=0. 选取可实现湍流模型确定模型内流场，其中 湍流动能 K 及湍流耗散率 ε (m2 ·s−3) 计算公式如下 所示[15] ： ∂ ∂xi (ρKvi) = ∂ ∂x j [(µ+ µt σK ) ∂K ∂xi ] −GK +Gε −ρε−YM （10） ∂(ρεvi) ∂xi = ∂ ∂xj ( µ+ µt σε ∂ε ∂xi ) +ρC1εS ε− ρC2ε ε 2 K + √ vε +C1ε ε K C3εGb （11） 式中：C1ε、C2ε、σK 及 σε 为计算常数，其数值分别 为 1.44、1.92、1.0 及 1.3[14–15] ；GK 及 Gb 分别代表因 平均速度梯度或耗散率产生的湍流动能源项 ， kg·m−1·s−3 ；μ 为层流黏度系数，Pa·s−1 ；Sε 为自定义 附加湍流动能源项，kg·m−1·s−3；浮力剪切层内主流 方向与重力方向一致时，C3ε 常数为 1；浮力剪切层 内主流方向垂直于重力方向时，C3ε 常数为 0. YM 为 可压缩湍流的脉动值，kg·m−1·s−3，其计算公式如式 （12）所示. YM = 2ρε K γRT （12） 式中：R 为气态常数，J·mol−1·K−1，数值为 8.314；T 为 气相温度，K. 2.2    网格划分及过程算法 图 2 为超音速氧枪的二维网格结构示意图. 考虑到拉瓦尔喷管具有几何对称与流动对称特 点，为降低求解过程计算量，选取喷管半截面及轴 对称处理进行模拟. 计算域均采用结构化网格，为 提高模拟结果准确性，对喷管近壁面、喉口附近及 射流速度梯度较大部分的区域进行网格加密处 理. 初始状态下，计算域充满静止不动的空气，表 1 及表 2 分别为模型边界条件及相关物化参数. 模型中空气与氧气的定压比热容 Cp (J·kg−1·K−1)， 分别按式（13）和（14）计算. Cp,air = 1161.48−2.36×T +1.49×10−2 ×T 2− 5.03×10−5 ×T 3 +9.93×10−8 ×T 4 −1.11×10−10× T 5 +6.54×10−14 ×T 6 −1.57×10−17 ×T 7 （13） 0 (a) (b) 0.02 0 5 0 20 0 500 mm mm mm mm 图 2    模型网格结构图. （a）拉瓦尔喷管网格示意图；（b）计算域网格 示意图 Fig.2     Schematic  of  the  mesh  distribution  of  the  simulation  model: (a) Laval nozzle; (b) computational domain · 56 · 工程科学学报，第 42 卷，增刊 1