射流阀门开/闭过程无网格法数值模拟

王园丁,郭曼丽,谭俊杰,任登凤

(1.上海空间推进研究所,上海 201112;2.上海空间发动机工程技术研究中心,上海 201112;3.南京理工大学 能源与动力工程学院,江苏 南京 210094)

姿轨控发动机内流场存在着非常复杂的非定常流动现象,尤其当各种射流阀门开/闭时其流场更加复杂。射流阀开/闭过程与射流控制技术密不可分。射流控制理论由美国Harry Diamond实验室于20世纪五十年代提出[1],并已经在工业工程中得到了较为广泛的应用。例如,在石油工业中用于钻采的射流式冲击锤[2]、兵器工业中用于改变导弹姿态的双稳换向阀[3]以及用于化工工业中电解铝温控系统的射流控制系统等[4]。射流阀通过控制端射流与主流的相互作用改变主流的流动特性,使主流流动方向发生偏转,并在某一个方向产生一定的推力,以改变受控体的姿态和轨迹,实现对受控体的姿态控制。

国内外在射流控制方面也做了大量研究工作,大部分工作都是围绕较基础的数值模拟展开的。国外约翰霍普金斯大学应用物理实验室的Roger等[5]、韩国航空航天大学的Heo[6]对射流阀进行了细致研究。北京理工大学是国内最早开展燃气射流阀研制的高校,并成功研制了兵器系统某远程火箭弹修偏用燃气射流阀。文献[7-9]通过数值模拟手段对启动射流阀的舵机压力特性、驱动机构的瞬态响应、喷口遮挡的影响等一系列问题进行了研究。文献[3]对某远程多管火箭弹姿控发动机双稳阀超音速射流流场进行数值模拟,基于A100远程多管火箭弹BK183-5姿控发动机50阀模型以及改进的42.5阀模型,采用隐式算法以及先进的V2F湍流模型对射流阀流场进行求解,通过与已有实验数据的对比,合理解释了实验中发生的现象,验证了利用数值模拟手段求解超音速射流阀的可靠性和准确性,取得了较为满意的结果。

与传统的网格方法相比,无网格方法不需要通过离散节点组成的网格单元进行求解区域的划分,而是通过“点云”结构代替常规的网格单元离散计算域,布点较为灵活,因而国内外开展了大量的研究工作[10-17],主要涉及基础算法研究,如将迎风格式引入无网格算法、求解各项异性点云结构的无网格算法等,但鲜有应用于姿轨控发动机射流阀流场的数值模拟。本文采用耦合Wilcoxk-ω湍流模型的最小二乘无网格,对含拓扑结构改变的姿轨控发动机超音速射流阀流场进行数值模拟,成功模拟出附壁流动、主流的切换等复杂的非定常流动状态。

1 射流阀基本原理

图1为射流换向阀结构示意图[5],主要由主流入口、控制通道、尖劈和排气通道组成,主流通道结构造型为拉瓦尔喷管。图中,b为主流通道宽度,d为喷管厚度,φ为扩张角。

射流换向阀的基本工作过程可以根据图1描述如下:主流从主流入口处进入漏斗形喷管后,经垂直段通道到达喉部,若打开左侧控制端的电磁阀,则控制流从左侧控制通道进入喷管,并在喉部对主流流动产生向右的作用力,主流流动将发生向右偏转,在扩张段形成右侧附壁流动,并经尖劈分流后从右侧排气通道流出。图1中尖劈结构起到限制分流的作用。若要使得主流发生向左偏转,则只需要关闭左侧控制端电磁阀,并打开右侧电磁阀使控制流从右控制端流入喷管,则控制流会在喉部对主流产生向左的作用力而在扩张段形成左侧附壁流,并沿左侧排气通道流出。

图1 射流换向阀结构示意图

2 射流阀开闭过程数值模拟

采用耦合Wilcoxk-ω湍流模型的最小二乘无网格,对含拓扑结构改变的姿轨控发动机超音速射流阀流场进行数值模拟。

2.1 物理模型

图2为不带阀芯装置的计算模型,射流口高度为1.8 mm,喉部主流入口宽度为1.8 mm,喉部出口宽度(扩张段入口)为2.1 mm,扩张半角为12°,z方向为5.0 mm(将z方向中间位置定义为z=0坐标起点)。对此模型进行数值模拟,可以对控制流与主流的相互作用进行初步研究,以观测在控制流作用下的主流切换效果。

图2 不带阀芯装置的计算模型(z=0截面模型尺寸)

图3为带阀芯装置的计算模型。主流流入射流阀后,控制流在喉部对主流作用,使得主流发生偏转,并对阀芯产生作用力驱动阀芯运动到另一侧,从而起到密封的作用。

图3 带阀芯装置的计算模型

2.2 控制方程

可压缩雷诺平均N-S控制方程可表示为如下守恒形式:

(1)

式中:Q为守恒变矢量;E,F,G分别为x,y,z坐标方向的无粘通量;Ev,Fv,Gv分别为x,y,z坐标方向的粘性通量。

Q=(ρuρvρwρe)T
E=(ρuρu2+pρuvρuw(ρe+p)u)T
F=(ρvρuvρv2+pρvw(ρe+p)v)T
G=(ρwρuwρvwρw2+p(ρe+p)w)T
Ev=(0τxxτxyτxzuτxx+vτxy+wτxz-qx)T
Fv=(0τxyτyyτyzuτxy+vτyy+wτyz-qy)T
Gv=(0τxzτyzτzuτxx+vτyz+wτzz-qz)T

ρ,p,e分别为流体密度、压强、单位质量的气体总内能;u,v,w分别为x,y,z方向的速度分量;γ为比热比,对于空气γ=1.4;τ为应力项;q为热通量项。具体表达式见文献[12]。

2.3 数值计算方法

2.3.1 空间导数求解

基于无网格方法进行流场计算首要工作便是确定空间导数的拟合方法,这也是无网格算法能否成功的基础。本文所采用的空间导数求解方法是基于Taylor级数展开而进行的。可以简单描述为:首先对中心点处的基本变量结合卫星点处的具体变量值进行泰勒展开,利用二次平方极小曲面逼近求解空间导数,进而利用最小二乘求导方法计算数值通量。

假设流场中任意离散点i处的流场值可以表示为函数f(x,y,z),且设函数f(x,y,z)在点云Ci内n阶连续可导,对点云Ci内任意一个卫星点s处的流场值f(xs,ys,zs),可以通过Taylor级数展开得到,即f(xs,ys,zs)可表示为

(2)

式中:i表示中心点,s表示以点i为中心点的点云结构中的卫星点,s∈Ci。由点云结构的构造方法可知,卫星点个数大于3,而式(2)中只有3个未知数,所以式(2)是一个矛盾方程组,此时2个未知的偏导数可以通过最小二乘法拟合得到。最小二乘法需要使得拟合曲线尽可能地逼近物理解,即需要满足监测数据与拟合曲线之间的偏差(亦可理解为式(2)中的高阶小量)的平方和最小:

f(O(h2))=

(3)

(4)

显式求解式(4),可将空间导数表示成如下形式:

(5)

方程组系数αis,βis和γis的具体表达形式为

(6)

式(6)中的系数具体表达式为

(7)

(8)

2.3.2 数值通量求解

本文借鉴文献[10]的思想,在点云结构中心点与卫星点连线的中点处构造一个虚拟界面J,将MUSCL(monotone upstream-centred schemes for conservation laws,MUSCL)格式引入无网格方法中,用以重构该界面左、右两侧的原始变量值,进而虚拟界面处的无粘通量可通过迎风格式进行求解。本文采用的迎风格式为AUSM+-up格式。

对于粘性通量中的速度和温度的梯度可以直接应用式(5)求解。为了避免因出现奇偶失联导致的计算不准或者错误的情况,对虚拟中点处的梯度值进行修正:

(9)

通过上述方法对控制方程进行离散之后,控制方程将变成常微分方程,本文采用三阶SSP型Runge-Kutta法,并使用当地时间步长进行时间推进求解离散后的控制方程。湍流模型采用Wilcoxk-ω两方程模型,针对各项异性点云结构,采用经过验证的点云重构技术[15]进行病态点云结构的修正,以提高模拟粘性流动的准确性。

2.4 计算条件

针对上述2种计算物理模型,设置主流入口压力为10 MPa,温度300 K;左右控制端控制流入口压力为6.0 MPa,温度300 K;不带阀芯装置的计算模型,出口压力为3.5 MPa,温度300 K,带阀芯装置下部左右两侧出口压力为1个大气压,温度300 K。不带阀芯装置的计算模型,整个计算域离散点分布为90万;带阀芯装置的整个计算域离散点分布为286万,其中阀芯周围离散点为56万。先打开主流,将流场计算稳定,再打开一侧控制端使控制流进入流场,进行流场计算。

2.5 计算结果分析

2.5.1 不带阀芯装置的数值模拟

1)稳态计算结果。

图4给出了不带阀芯装置模型稳态计算过程中马赫数的变化情况。此时左右通道控制端做固壁处理。

图4 稳态计算结果(马赫数云图)

从图4中可以看出,随着主流的流入,开始会形成左右对称的流场结构,如图4(a)所示,随着时间的推进,根据射流控制理论及主流自身紊流特性,主流会向右侧发生偏转,如图4(b)所示,并最终形成稳定的右侧附壁流动。此时左右控制端的压力为4.65 MPa左右。在稳态计算结果的基础之上,如图4(c)所示,打开右端控制阀进行非稳态计算。

2)右侧控制端打开的计算结果。

图5为在稳态基础之上打开右端控制阀后,计算得到的马赫数随时间变化的云图,此时左通道控制端做固壁处理。

图5 右侧控制阀打开后马赫数变化历程

从图5中可以看出,在t=0.4 ms时,主流便发生了明显偏转;在t=0.8 ms时,主流已经在下部经过分流尖劈从右出口逐渐转向左出口;在t=1.2 ms时,大部分的主流已经经由下部左侧出口流出;t=1.6 ms时,主流已经实现从右侧到左侧的切换,只是流场还没有最终稳定下来;在2.5 ms时气流实现了从右侧附壁流到左侧附壁流的完全切换。完全切换后的左通道控制端附近的压力为5.6 MPa,右通道控制端附近的压力为6.0 MPa。在此基础之上,打开左通道控制端,关闭右通道控制端(做固壁处理),进行二次切换的非稳态计算。

3)左侧控制端打开的计算结果。

在气流完全切换至稳定的左侧附壁流动后,关闭右侧控制端,打开左侧控制端,对流场进行计算,所得的马赫数随时间变化的云图如图6所示。

图6 左侧控制阀打开后马赫数变化历程

从图6中可以看出,在6.0 MPa的控制压力下,主流也实现了偏转,但因为此时左右通道压差较小,所以需要更长的时间实现主流的完全切换。主流实现完全切换的时间为6.5 ms。

对不带阀芯装置的模型进行计算,在右通道控制端6.0 MPa压力的作用下,主流从右侧附壁流在2.5 ms内成功切换至左侧附壁流。在此基础之上,关闭右通道控制端,打开左通道控制端让6.0 MPa控制流对主流进行作用,成功计算出了二次切换的流场结构图,所需要的时间为6.5 ms。说明在现有计算模型和计算条件下可以实现射流阀的双向切换。

2.5.2 带阀芯装置的数值模拟

针对带阀芯装置的模型,设置主流入口压力为10.0 MPa,将主流控制端阀门打开,关闭左右控制端(做固壁处理),计算流场至稳定状态,速度云图如图7所示。由于阀芯和阀体之间间隙小,在最小的喉道位置处气体流动速度将为音速。由于阀芯与阀体完全接触,左侧排气通道被完全关闭,气流无法通过左侧排气通道流出,只能通过右侧排气通道流出,此时左右通道端部压力为4.35 MPa,因此要实现左右控制端气流顺利进入通道驱动主流发生偏转,控制端气流需大于4.35 MPa。

图7 带阀芯流场图

3 结束语

采用耦合Wilcoxk-ω湍流模型的最小二乘无网格,对含拓扑结构改变的姿轨控发动机超音速射流阀流场进行数值模拟。成功模拟出附壁流动、主流的切换等复杂的非定常流动状态,捕捉到了较为清晰的流场结构与流动特征,能够客观地反应射流阀的基本切换过程。数值模拟结果为姿轨控发动机射流阀门的设计提供了参考,也扩大了耦合两方程湍流模型的最小二乘无网格方法数值模拟的范围。