一种高精度暂冲型引射式跨声速风洞流场控制结构

摘要

本发明涉及一种高精度暂冲型引射式跨声速风洞流场控制结构及方法。该控制结构将风洞流场控制分解为主引压力控制和总压静压控制两部分，由设定值计算模型、主引控制器、流场控制器三个模块构成。在主引压力控制部分，采用主调压阀控制主引压力。在总压静压控制部分，通过主排气阀和栅指控制总压和静压。设定值计算模型由一维变截面管道非定常欧拉方程组构成，通过计算流体力学软件FLUENT对风洞模型求解主引压力设定值和静压设定值。使用分段增量式PID控制算法构成主引控制器，调节主引射器出口处的主引压力。以模型姿态前馈动态矩阵控制方法为核心构成流场控制器，调节稳定段的总压和试验段的静压。该控制方法实时性好，控制精度高，试验阶段马赫数的变化控制到了0.001的范围。

说明书

技术领域

本发明涉及一种暂冲型引射式跨声速风洞流场的高精度控制结构及方法， 流场马赫数控制精度可达0.001。

背景技术

风洞是一种产生可控均匀气流的管状试验装置，是研究空气流动规律、航 空航天飞行器及其它物体气动特性的重要地面模拟平台。风洞试验作为获取飞 行器气动力数据、评估分析布局设计和预测飞行性能的主要手段，其试验精准 度直接关系到飞行器布局方案的优化选择和对实际飞行性能的准确预测，从而 关系到飞行器的研制质量和研制水平。

作为进行空气动力学试验的主要平台，暂冲型引射式跨声速风洞需要解决 流场高精度控制以及降低能耗等问题。大飞机(我国C919大型客机、ARJ飞机 等)等飞行器的研制对风洞流场马赫数的控制精度提出了严苛的要求，试验过 程中流场马赫数的精度(准确度和精确度)都必须达到0.001，而作为大飞机 这类我国自主在研项目唯一可依赖的跨声速气动力试验平台(2.4m暂冲式跨 声速风洞)的马赫数控制精度仅为0.002～0.003。另外，风洞属于高能耗的设 备。以国内2.4m暂冲式跨声速风洞为例，运行时间一般为45～90秒，试验车 次超过4000次/年。按照一次试验60秒进行计算，能源动力成本高达上万元。

暂冲型引射式风洞主要有以下控制难点：①扰动频繁，主要体现在两个方 面：设定值变化和模型攻角变化。空气动力学试验要求在试验阶段马赫数需要 阶梯或连续变化，从而带来设定值阶跃扰动；攻角阶梯或连续变化也会带来扰 动。所以，如何能够快速准确地消除扰动成为主要解决的控制问题。②时变， 暂冲型引射式风洞以储气罐中的气体作为气源，由于气源容量的限制，气源压 力随时间逐渐下降，导致风洞系统的模型发生变化。如何能够尽量地消除气源 压力对风洞流场的影响也是需要解决的控制问题。③时滞，2.4米风洞的主体 结构是一个全长近200米的环形管道，控制量和被控量在实际物理位置上存在 一定距离，这导致控制量的变化不能直接作用于被控量，例如主排气阀与稳定 段之间的距离有近60米，使得主排气阀对总压的作用存在滞后。④耦合严重， 在风洞系统中，各控制量和被控量之间存在严重的耦合关系。如在使用主排气 阀调节风洞流场时，总压和马赫数都会受到影响；同理，在使用栅指调节风洞 流场时，总压和马赫数也都会发生变化。因此，在设计风洞流场控制器时，必 须同时考虑主排气阀、栅指对总压和马赫数的影响。⑤非线性严重，暂冲型引 射式风洞以总压和马赫数表征风洞流场品质，传统的控制系统都采用总压和马 赫数作为被控参数。但是，由于马赫数不能直接测量，根据亚音速一维等熵流 理论，根据马赫数由总压和静压经公式计算得出，但这种关系是非线性的，不 利于流场精确控制。

发明内容

本发明的目的是针对暂冲型引射式跨声速风洞流场中存在的控制问题，建 立了一种带有解耦、前馈功能的高精度风洞流场控制结构。

一种高精度暂冲型引射式跨声速风洞流场控制结构，由主引压力控制和总 压、静压控制两部分组成，包括设定值计算模型、主引控制器、流场控制器； 所述设定值计算模型的输出端一端输出到主引控制器，主引控制器输出到主调 压阀子系统用于调压，调压后由主引射器输出，且输出的数据传输回主引控制 器；所述设定值计算模型的输出端另一端输出两路到流场控制器，流场控制器， 输出两路，一路输出到主排气阀子系统，一路输出到栅指系统，主排气阀子系 统和栅指系统控制风洞流场分别输出总压和静压，且输出的总压和静压数据传 输到流场控制器。

设定值计算模型由一维变截面管道非定常欧拉方程组构成，该模型用于根 据总压设定值和马赫数设定值计算主引压力设定值和静压设定值。

所述的设定值计算模型的建立步骤如下：

步骤1.建立风洞进气处一维变截面管道非定常欧拉方程组。

利用一维等熵流理论计算主引射器喷嘴处的气体压力，密度和速度，在此 基础上建立主调压阀子系统处的质量方程、动量方程和能量方程。

步骤2.建立风洞排气处一维变截面管道非定常欧拉方程组。

利用一维等熵流理论计算主排气阀处的气体压力，密度和速度，在此基础 上建立主排气阀子系统处的质量方程、动量方程和能量方程。

步骤3.建立风洞气源、集气室气动力方程，建立风洞模型。

所述设定值计算模型的计算过程为：使用计算流体力学软件FLUENT对 风洞模型求解。采用非稳态隐式求解，粘性设置为无粘。压力速度耦合采用 SIMPLE算法，模型输入量为入口边界总压、静压以及出口边界静压，在稳定 段和试验段添加平面作为监测面以获得所需输出量。模型中压力入口边界参数 通过瞬态Profile导入到FLUENT中，软件根据时间步的推进而读取对应数据 作为压力入口边界的参数值。压力出口参数设置为恒定压力(当地大气压)， 动量方程的离散格式选择二阶迎风格式以减少数值耗散的影响，从而计算出主 引压力设定值和静压设定值。

使用分段增量式PID控制算法构成主引控制器，调节主引射器出口处的主 引压力。所述的分段增量式PID控制算法将控制量分成不同的阶段，根据不同 的偏差采用不同的调节参数，在大偏差时加快调节速度，小偏差时减少超调。

以模型姿态前馈动态矩阵控制方法为核心构成流场控制器，调节稳定段的 总压和试验段的静压。所述的模型姿态前馈动态矩阵控制方法的计算步骤如下：

步骤1.根据风洞阶跃响应建立模型向量；

步骤2.构造一步动态矩阵和多步动态矩阵A；

步骤3.计算控制矩阵D；

步骤4.计算k-1时刻的被控量一步预测值

步骤5.采集数据。采集k时刻的总压、静压y₁(k)、y₂(k)；

步骤6.利用控制偏差补偿k时刻的被控量一步预测值

步骤7.利用模型姿态对被控量一步预测值进行前馈补偿修；

按步骤1构造模型姿态变化对总压及模型姿态变化对静压的模型向量，记 为b₁、b₂，其中，b_i＝[b_i(1),b_i(2),…,b_i(N)]^T,i＝1,2；

按步骤2构造一步扰动动态矩阵具体如下：然后计算扰 动前馈补偿量记k时刻的模型姿态变化量为△v(k)，则扰动前馈补偿量 为：

$\Delta {\bar{y}}_d\left(k\right)=\bar{B}\mathrm{\Delta v}\left(k\right)---\left(1\right)$

在步骤6的基础上，将扰动前馈补偿量叠加到完成对的修正，公式如下：

${\bar{y}}_{\mathrm{cord}}\left(k\right)={\bar{y}}_{\mathrm{cor}}\left(k\right)+\Delta {\bar{y}}_d\left(k\right)---\left(2\right)$

步骤8.计算k时刻的被控量初始预测值

步骤9.提取预测时域P内的被控量初始预测值

步骤10.获取未来P个时刻的总压、静压设定值序列w(k)；

步骤11.计算k时刻的由主排气阀、栅指构成的控制向量u(k)，在将k时 刻的控制量施加到风洞后，转到步骤4准备计算下一时刻的控制量。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

提出一种新的控制结构，通过引入主引压力控制回路，减少气源压力下降 对风洞流场的影响，有效解决了时变问题。

将被控参数从总压、马赫数转换为总压、静压、主引压力，降低系统的非 线性度，减轻了被控参数之间的耦合，有利于提高控制性能。

提出一种适合风洞流场控制的模型姿态前馈动态矩阵控制方法，实现了风 洞流场的抑制模型姿态扰动和优化解耦控制。

风洞流场控制精度得到提高，如马赫数控制精度提高至0.001，总压控制精 度提高至0.1％，满足了大飞机的实验要求。

风洞流场控制精度的提高，使得实验过程中数据测量等待时间减少，从而 加快了实验的进程，缩短实验时间，降低了能耗。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1是暂冲型引射式跨声速风洞控制系统图；

图2是本发明的暂冲型引射式跨声速风洞控制结构图。

具体实施方式

本发明以国内2.4m暂冲型引射式跨声速风洞为例，说明所发明的控制结 构和方法的实施方式。

2.4m暂冲型引射式跨声速风洞流场控制系统结构如图1所示，虚线框中为 广义风洞对象。为消除风洞试验过程中气源压力的降低对风洞流场的影响，引 入主引压力控制回路，将风洞流场控制分解为主引压力控制和总压静压控制两 部分，如图2所示。该控制结构由设定值计算模型、主引控制器、流场控制器 三个模块构成。在主引压力控制部分，采用主调压阀控制主引压力。在总压静 压控制部分，通过主排气阀和栅指控制总压和静压。具体实施步骤如下：

步骤一、建立设定值计算模型

在使用图2所示的控制结构调节风洞之前，首先应建立设定值计算模型， 建立过程为：

步骤1.建立风洞进气处一维变截面管道非定常欧拉方程组

在主调压阀子系统处进气对风洞空气环流的影响相当于给风洞空气环流 注入质量、动量和能量。由于储气罐压力与风洞压力之比总高于临界值，即主 调压阀阀门喉道处(最小面积处)始终处于音速状态，因此相应的源项可表示 为：

注入的质量：$S_1=P_0·A^{*}\sqrt{\frac{K}{{\mathrm{RT}}_0}{\left(\frac{2}{K+1}\right)}^{\frac{K+1}{K-1}}}$

注入的动量：S₂＝S₁U_e+P_eA_e

注入的能量：$S_3=S_1(\frac{K}{K-1}{P}_e+\frac{1}{2}{\rho }_e{U}_e^2)/{\rho }_e$

其中，A_e是主引射器所有喷嘴总面积，P_e，ρ_e，U_e分别是主引射器喷嘴 处的压力，密度和速度，A^*为主调压阀门的喉道截面积。由面积比A_e/A^*，用 下述一维等熵流理论可计算出P_e，ρ_e，U_e。

$\frac{A_e}{A^{*}}=\frac{1}{M}\{\frac{2}{K+1}(1+\frac{K-1}{2}{M}^2){\}}^{\frac{K+1}{2(K-1)}}$

由上式可计算出M数，再由下面三个公式可计算出P_e，ρ_e，U_e。

$\frac{P_e}{P_0}={(1+\frac{K-1}{2}{M}^2)}^{-\frac{K}{K-1}}$

$\frac{{\rho }_e}{{\rho }_0}={(1+\frac{K-1}{2}{M}^2)}^{-\frac{1}{K-1}}$

$U_e=M\sqrt{K\frac{P_e}{{\rho }_e}}$

其中，P₀和ρ₀为集气室内的总压与密度。

步骤2.建立风洞排气处一维变截面管道非定常欧拉方程组

在主排气阀子系统排气对风洞空气环流的影响相当于从风洞空气环流中 抽取质量、动量和能量。相应源项表达式为：

$S_1=-P_{0e}·A_e^{*}\sqrt{\frac{K}{{\mathrm{RT}}_{0e}}{\left(\frac{2}{K+1}\right)}^{\frac{K+1}{K-1}}}$

其中，P_0e和T_0e是风洞管道上排气处的总压和总温，是排气系统喉道面积。 上式只有当排气系统喉道面积处达音速时才适用。风洞中压力与排气系统出口 处的压力比未达到临界值时，可按亚音速一维等熵流理论公式计算，具体计算 如下：

$M=\sqrt{\frac{2}{K-1}[{\left(\frac{P_{0e}}{P_b}\right)}^{\frac{K-1}{1}}-1]}$

${\rho }_b={\rho }_e{\left(1\frac{K-1}{2}{M}^2\right)}^{-\frac{1}{K-1}}$

$U_b=M\sqrt{K\frac{P_b}{{\rho }_b}}$

S₁＝ρ_bU_bA_b

其中，P_b，ρ_b，U_b分别是排气系统出口处的压力、密度和速度。ρ_e是风洞 管道上排气处的密度，由于排气方向与风洞气流方向垂直，所以在风洞气流方 向上动量损失为零，即

S₂＝0

$S_3=-S_1(\frac{K}{K-1}{P}_e+\frac{1}{2}{\rho }_e{U}_e^2)/{\rho }_e$

上式中U_e是风洞管道排气处的速度，P_e是风洞管道排气处的静压，并且该 式也只适用于排气阀门喉道面积处为音速时的情况，对亚临界情况，仍然要用 亚音速一维等熵流公式计算出ρ_b、P_b和U_b，计算出S₁，进而利用上式计算出S₃。

步骤3.建立风洞气源、集气室气动力方程，建立风洞模型

假设：

(Ⅰ)气源、集气室容积内气体参数均匀分布。

(Ⅱ)气源、集气室容积内各质点速度U＝0。

对容积V内的气体应用积分形式气体动力学方程组：

应用假设(1)、(2)对气源边界得：

$\frac{{\partial \rho }_s}{\partial t}=-\frac{\stackrel{.}{m}}{V_S}$

U_s＝0

$\frac{\partial }{\partial t}(e_s-P_s)=-\frac{{\stackrel{.}{m}}_A}{V_S}(\frac{K}{K-1}\frac{P_s}{{\rho }_s}+\frac{1}{2}{U}_s^2)$

其中，为流出气源的气体质量流率，值为ρ_sU_sA_s，V_s为气源容积，P_s，ρ_s，U_s， A_s分别为与气源结合处得送气管道截面上的压力、密度、速度和截面积。应用 假设(1)、(2)对集气室边界得：

$\frac{{\partial \rho }_m}{\partial t}=-\frac{{\stackrel{.}{m}}_C-{\stackrel{.}{m}}_B}{V_m}$

U_m＝0

$\frac{\partial }{\partial t}(e_m-P_m)=-\frac{{\stackrel{.}{m}}_C}{V_m}(\frac{K}{K-1}\frac{P_C}{{\rho }_C}+\frac{1}{2}{U}_C^2)-\frac{{\stackrel{.}{m}}_g}{V_m}(\frac{K}{K-1}\frac{P_B}{{\rho }_B}+\frac{1}{2}{U}_B^2)$

其中，分别为流入、流出集气室的气体质量流率，值分别为ρ_CU_CA_C和 ρ_BU_BA_B，V_m为集气室容积，X_C为与集气室结合处的送气管道截面上的参数，X_B为与集气室结合处的喷管截面上的参数。

步骤二、计算主引压力设定值和静压设定值

建立好风洞模型后，根据设定值模型计算主引压力设定值和静压设定值， 其计算过程为：使用计算流体力学软件FLUENT对风洞模型求解。采用非稳态 隐式求解，粘性设置为无粘。压力速度耦合采用SIMPLE算法，模型输入量为 入口边界总压、静压以及出口边界静压，在稳定段和试验段添加平面作为监测 面以获得所需输出量。模型中压力入口边界参数通过瞬态Profile导入到 FLUENT中，软件根据时间步的推进而读取对应数据作为压力入口边界的参数值。 压力出口参数设置为恒定压力(当地大气压)，动量方程的离散格式选择二阶 迎风格式以减少数值耗散的影响，从而计算出主引压力设定值和静压设定值。

步骤三、求取控制率

在每个控制周期(假定控制周期为k)内，求取主引射器、主排气阀和栅 指控制率的过程如下：

首先，使用分段增量式PID控制算法求取主调压阀控制率，根据不同的偏 差采用不同的调节参数，在大偏差时加快调节速度，小偏差时减少超调。增量 式PID控制算法如下：

$\left(\begin{array}{c}u\left(k\right)=u(k-1)+K_p[e\left(k\right)-e(k-1)]+\\ K_p\frac{T_s}{T_i}e\left(k\right)+K_p\frac{T_d}{T_S}[e\left(k\right)-2e(k-1)+e(k-2)]\end{array}\right)---\left(3\right)$

其中，e为控制偏差，K_p为比例系数，T_i为积分时间常数，T_d为微分时间常数。

使用模型姿态前馈动态矩阵控制方法求取主排气阀和栅指的控制率，计算 步骤如下：

步骤1.根据风洞阶跃响应建立模型向量。首先对采集到的阶跃响应数据 单位化，然后利用输入输出之间的单位阶跃响应分别构造主排气阀对总压、主 排气阀对静压、栅指对总压和栅指对静压的模型向量，记为a₁₁、a₂₁、a₁₂、a₂₂。 下面以a₁₁为例说明模型向量的构造过程：令a₁₁＝[a₁₁(1)…a₁₁(N)]^T，其中a₁₁(i)为 主排气阀对总压的单位阶跃响应在iT时刻的采样值，i＝1,2,…N₁₁，T为采样周 期，N₁₁为建模时域。为了便于计算，将风洞变量的建模时域统一选取为N。 按照相似的方法，可以构造a₂₁、a₁₂、a₂₂。

步骤2.构造一步动态矩阵和多步动态矩阵A。动态矩阵可以依据模型 向量构造，具体如下：

$\bar{A}=\left(\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right),A=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right)---\left(4\right)$

其中，i＝1,2；j＝1,2；P为预测时域；M为控 制时域。P表示预测未来P个时刻的被控量，M表示所预测的控制量是根据未 来M个时刻的控制增量计算的。

步骤3.计算控制矩阵D。

D＝(A^TQA+R)^-1A^TQ     (5)

其中，Q为误差权矩阵，R为误差权矩阵。

步骤4.计算k-1时刻的被控量一步预测值表示：当 k-1时刻控制有一增量△u(k-1)时，被控量未来N个时刻的预测值。计算公式如 下：

${\bar{y}}_{N1}(k-1)={\bar{y}}_{N0}(k-1)+\bar{A}\mathrm{\Delta u}(k-1)---\left(6\right)$

其中，表示:在k-1时刻，假定控制作用不变时对未来N个时刻的被 控量初始预测值；

步骤5.采集数据。采集k时刻的总压、静压y₁(k)、y₂(k)。

步骤6.利用控制偏差补偿k时刻的被控量一步预测值

${\bar{y}}_{\mathrm{cor}}\left(k\right)={\bar{y}}_{N1}(k-1)+\mathrm{He}\left(k\right)---\left(7\right)$

其中，e(k)为k时刻的预测误差，e(k)为k时刻被控量真实值与k时刻被控量预 测值的差值，H为误差权矩阵。

步骤7.利用模型姿态对被控量一步预测值进行前馈补偿修。

按步骤1构造模型姿态变化对总压及模型姿态变化对静压的模型向量，记 为b₁、b₂。其中，b_i＝[b_i(1),b_i(2),…,b_i(N)]^T,i＝1,2。按步骤2构造一步扰动动 态矩阵具体如下：然后计算扰动前馈补偿量记k时刻 的模型姿态变化量为△v(k)，则扰动前馈补偿量为：

$\Delta {\bar{y}}_d\left(k\right)=\bar{B}\mathrm{\Delta v}\left(k\right)---\left(8\right)$

在步骤6的基础上，将扰动前馈补偿量△y_d(k)叠加到完成对的修正，公式如下：

${\bar{y}}_{\mathrm{cord}}\left(k\right)={\bar{y}}_{\mathrm{cor}}\left(k\right)+\Delta {\bar{y}}_d\left(k\right)---\left(9\right)$

步骤8.计算k时刻的被控量初始预测值计算公式如下：

${\bar{y}}_{N0}\left(k\right)=S_0{\bar{y}}_{\mathrm{cord}}\left(k\right)---\left(10\right)$

其中，S₀为移位矩阵，$S_0=\left(\begin{array}{cc}S& 0\\ 0& S\end{array}\right),$

步骤9.提取预测时域P内的被控量初始预测值计算公式如下：

${\bar{y}}_{P0}\left(k\right)=I_0{\bar{y}}_{N0}\left(k\right)---\left(11\right)$

其中，I₀为移位矩阵，$I_0={\left(\begin{array}{cc}I& 0\\ 0& I\end{array}\right)}_{2P\times 2N},$

步骤10.获取未来P个时刻的总压、静压设定值序列w(k)。

w(k)＝[sP₀(k),sP_S(k)]^T     (12)

其中，

sP₀(k)＝[sP₀(k+1),…,sP₀(k+P)]^T为总压设定值，

sP_S(k)＝[sP_S(k+1),…,sP_S(k+P)]^T为静压设定值。

步骤11.由主排气阀、栅指构成的控制向量u(k)。

$\mathrm{\Delta u}\left(k\right)=D[w\left(k\right)-{\bar{y}}_{P0}\left(k\right)]---\left(13\right)$

u(k)＝u(k-1)+△u(k)      (14)

在将k时刻的控制量施加到风洞后，转到步骤4准备计算下一时刻的控制 量。

试验效果：在国内2.4m暂冲型引射式跨声速风洞上进行控制效果的测试。 考察在吹风实验过程中的马赫数和总压的控制精度。风洞试验的不同工况对应 不同的马赫数、总压设定值，设定值计算模型求出的静压和主引压力的设定值 如表1。

表1不同工况对应的被控变量设定值

分别选择两种阻塞度不同的飞行器模型进行控制效果对比试验，模型A的 阻塞度为1％，模型B的阻塞度为0.8％。实验过程中，飞行器模型从-4°到10° 阶梯变攻角，攻角变化速度2°/s，阶梯4°。选择四种吹风工况进行对比， 采用控制过程中被控量的平均值和标准差两个指标衡量风洞控制系统的准度 和精度。控制效果如表2所示。

表2实验阶段应用验证控制性能对比表

表2显示，基于所发明的控制结构及方法，风洞试验过程中马赫数控制精 度提高至0.001，总压控制精度提高至0.1％。

本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中 披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何 新的组合。