CN106126791A - 一种考虑几何不确定性的高超声速机翼气动力/热分析方法 - Google Patents

Abstract

本发明公开了一种考虑几何不确定性的高超声速机翼气动力/热分析方法。根据伯恩斯坦多项式近似模型构造不确定参数的样本点，在给定飞行工况下，利用气动力/热工程算法计算样本点对应机翼外形的气动力/热响应值，同时应用最小二乘法得到多项式近似模型的拟合系数。在此基础上，利用已建立的多项式近似模型求解气动力/热响应值关于各个不确定参数的极大/小值点，并组合形成最大/小值点的向量，最终得到气动力/热响应值的区间上界和区间下界，实现考虑几何不确定性的高超声速机翼气动力/热分析。本发明方法得到的区间边界与蒙特卡洛方法得到的区间边界吻合较好，并且可以实现区间包络，为高超声速机翼外形的总体设计提供了新思路。

Description

一种考虑几何不确定性的高超声速机翼气动力/热分析方法

技术领域

本发明涉及高超声速飞行器设计技术领域，特别涉及一种考虑几何不确定性的高超声速机翼气动力/热分析方法。

背景技术

气动外形建模方法是高超声速飞行器设计与研制过程中的重要环节。为了提高建模效率和精度，一般需要对设计参数的数目进行控制，并且采用参数化建模方法。以机翼外形的几何建模方法为例，传统的参数化几何建模主要依赖于CAD方法(B样条曲线、NURBS曲线和二次曲线方法)，但该建模方法由于表达参数偏多，不利于概念设计中各学科分析模型的简化。相比而言，基于类函数/形函数转换的CST(Class function and Shape functionTransformation)方法通过修改少量的设计变量，可以生成翼型和机翼几何外形，而且具有快速、准确的拟合效果。

然而，由于加工工艺、制造水平的限制，高超声速飞行器的气动外形存在由加工误差导致的几何不确定性，几何设计参数往往在名义值附近上下波动。对于由加工误差引起的几何不确定性，现有的研究常采用概率方法对不确定因素进行定量化表征，利用标准正态函数描述几何参数的分布特征。但是，当没有足够的数据来验证这些随机变量概率密度的正确性时，概率方法难以可靠地满足精度要求的计算结果。与概率方法相比，非概率区间方法仅需明确不确定参数的分布界限，能够在不确定参数概率密度未知的情况下，利用区间向量对不确定参数进行定量化。非概率区间方法在结构的静、动力特性分析领域已经取得了一定成果，但在高超声速飞行器机翼气动外形建模中还处于起步阶段，相关的研究成果十分有限。同时，当机翼气动外形存在不确定性时，尤其是机翼前缘位置出现的几何波动，会造成机翼气动力/气动热性能存在偏差，这种波动或偏差对飞行安全是不利的。综上所述，亟需发展一种考虑几何不确定性的高超声速机翼气动力/热非概率区间分析方法，以克服传统概率方法对参数大样本容量试验数据的依赖性。

发明内容

本发明要解决的技术问题为：针对高超声速机翼在加工、制造过程中由加工误差所导致的几何不确定性，提出一种考虑几何不确定性的高超声速机翼气动力/热分析方法。该方法在CST参数化建模技术基础上，利用区间向量对机翼几何特征参数进行定量化表征。在给定工况下，利用气动力/热工程算法计算样本点处的气动力/热响应值，通过伯恩斯坦多项式近似模型及最小二乘法建立气动力/热响应值与不确定参数之间的映射关系。在给定区间范围内，利用已建立的多项式近似模型求解气动力/热响应值关于各个不确定参数最大/小值点的向量，计算得到气动力/热响应值的区间上界和区间下界，实现考虑几何不确定性的高超声速机翼气动力/热分析。

本发明解决上述技术问题采用的技术方案为：一种考虑几何不确定性的高超声速机翼气动力/热分析方法，包括以下步骤：

步骤(1)、确定高超声速机翼的平面轮廓参数，包括翼根弦长C_r、翼尖弦长C_t、副翼宽度C_a及机翼半展长L；

步骤(2)、利用类函数/形函数转换方法建立高超声速机翼翼型的参数化表达式，翼型的几何曲线可用下列函数表示：

$\frac{y}{c}\left(\frac{x}{c}\right)=C\left(\frac{x}{c}\right)S\left(\frac{x}{c}\right)+\frac{x}{c}\frac{z_{te}}{c}---\left(1\right)$

式中，x/c为翼型弦向的无量纲坐标值，y/c为翼型法向的无量纲坐标值，C(x/c)和S(x/c)分别为类型函数和形状函数，z_te/c为翼型后缘点的无量纲坐标值，C(x/c)可表示为：

$C\left(\frac{x}{c}\right)={\left(\frac{x}{c}\right)}^{N_1}{(1-\frac{x}{c})}^{N_2},0\≤\frac{x}{c}\≤1---\left(2\right)$

对于Clark-Ys翼型，取指数N₁＝0.5,N₂＝1，S(x/c)可表示为：

$S\left(\frac{x}{c}\right)=\underset{i=0}{\overset{n}{\Σ}}\[{\mathrm{\λ}}_i\·\frac{n!}{i!(n-i)!}\·{\left(\frac{x}{c}\right)}^i\·{(1-\frac{x}{c})}^{n-i}\],0\≤\frac{x}{c}\≤1---\left(3\right)$

S(x/c)与翼型前缘半径R_le/c及后缘倾角β满足以下关系：

$\left\{\begin{array}{c}S\left(0\right)={\mathrm{\λ}}_0=\sqrt{\frac{2R_{le}}{c}},\\ S\left(1\right)={\mathrm{\λ}}_n=tan\mathrm{\β}+\frac{z_{te}}{c}\end{array}\right.---\left(4\right)$

在形状函数S(x/c)中取n＝3，这样共包含等8个几何设计参数，其中，表示翼型前缘半径，β₁和β₂表示上、下翼面的后缘倾角，z_te/c为翼型后缘点的无量纲坐标值，λ₁、λ₂、λ₁'、λ₂'为上、下翼面形状函数多项式的加权系数；

步骤(3)、考虑由加工误差导致翼型几何特征参数存在的波动，将存在不确定性的几何特征参数记为其余参数均视为确定性参数；

步骤(4)、设定翼型几何特征参数α的区间上下界，可表示为：

${\mathrm{\α}}^I=\[\underset{\‾}{\mathrm{\α}},\stackrel{\‾}{\mathrm{\α}}\]---\left(5\right)$

式中，为高超声速机翼翼型几何特征参数的区间下界，为区间上界；

步骤(5)、利用区间[-1,1]内的n阶伯恩斯坦多项式，构造用于计算多项式系数的样本点，将伯恩斯坦多项式模型中子多项式的极值点t记为：

t＝[t₀,t₁,…,t_r,…,t_n-1,t_n]^T (6)

式中：

$t_r=\frac{2r-n}{n},r=0,1,...,n---\left(7\right)$

步骤(6)、将式(6)所表示的极值点映射至区间参数空间α^I内的第s个区间参数α_s，则在第s维用于计算该维近似多项式系数的样本点可以表示为：

${\mathrm{\α}}_{s,r}={\[{\mathrm{\α}}_1^c,...,{\mathrm{\α}}_{s-1}^c,{\mathrm{\α}}_{s,r},{\mathrm{\α}}_{s+1}^c,...,{\mathrm{\α}}_m^c\]}^T,r=0,1,...,n---\left(8\right)$

其中，

$\begin{array}{c}{\mathrm{\α}}_i^c=\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_i+{\underset{\‾}{\mathrm{\α}}}_i}{2}\right),i=1,2,\mathrm{...},s-1,s+1,\mathrm{...},m\\ {\mathrm{\α}}_{s,r}=\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s+{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)+\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s+{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)t_r\end{array}---\left(9\right)$

式中，m为不确定参数的个数；

步骤(7)、根据式(8)中样本点对应翼型的几何特征参数，建立高超声速机翼的气动外形，进而利用高超声速工程算法计算气动力/热响应值，记为Q_s，表示为如下形式：

Q_s＝[Q_s,0,...,Q_s,r,...,Q_s,n]^T (10)

式中，Q_s,r＝Q(α_s,r)表示样本点α_s,r对应的气动力/热响应值；

步骤(8)、利用最小二乘法计算伯恩斯坦多项式近似模型的待定系数ω_s，公式如下：

ω_s＝[ω_s,0,ω_s,1,...,ω_s,n]^T＝(B^TB)^-1.B^T·Q_s (11)

式中，B为伯恩斯坦-范德蒙矩阵，可表示为：

其中，b_n,r(t)为伯恩斯坦多项式的子多项式，记为：

$b_{n,r}\left(t\right)=\frac{n!}{r!(n-r)!}{\left(\frac{t+1}{2}\right)}^r{\left(\frac{-t+1}{2}\right)}^{n-r},r=0,1,\mathrm{...},n---\left(13\right)$

步骤(9)、根据步骤(8)中得到的系数ω_s，建立气动力/热响应值Q关于第s个不确定参数的伯恩斯坦多项式近似模型P_s(t)，表示如下：

$P_s\left(t\right)=\underset{r=0}{\overset{n}{\Σ}}{\mathrm{\ω}}_{s,r}\·b_{n,r}\left(t\right),t\∈\[-1,1\]---\left(14\right)$

步骤(10)、步骤(9)中P_s(t)的最值点可以通过其导函数的零点与自变量的端点产生，其中导函数零点可通过下式计算：

$\frac{{\mathrm{dP}}_s\left(t\right)}{dt}=0---\left(15\right)$

将零点记为t_s＝[t_s,1,t_s,2,...,t_s,n]，根据P_s(t)的自变量范围，可将零点进一步修改为：

其中Im和Re分别表示对应变量的虚部和实部；

步骤(11)、气动力/热响应值Q关于第s个区间参数的最小值点和最大值点可计算为：

$\begin{array}{c}{P}_s\left(t_{s,\min }\right)=\underset{t\∈T_s}{\min }\left\{P_s\left(t\right)\right\}\\ P_s\left(t_{s,\max }\right)=\underset{t\∈T_s}{\min }\left\{P_s\left(t\right)\right\}\end{array}---\left(17\right)$

式中，T_s＝[t_s,-1,1]；

步骤(12)、重复步骤(6)～(11)，以相同的方式计算气动力/热响应值Q关于所有区间参数的最小值点和最大值点，最终可以得到气动力/热响应值在标准区间[-1,1]内的最小值点和最大值点，即：

$\begin{array}{c}{t}_{\min }={\[t_{1,\min },t_{2,\min },\mathrm{...},t_{s,\min },\mathrm{...},t_{m,\min }\]}^T\\ t_{\max }={\[t_{1,\max },t_{2,\max },\mathrm{...},t_{s,\max },\mathrm{...},t_{m,\max }\]}^T\end{array}---\left(18\right)$

步骤(13)、将步骤(12)中标准区间内的最小值点和最大值点映射到实际参数空间有：

$\begin{array}{c}{\mathrm{\α}}_{\min }={\[{\mathrm{\α}}_{1,\min },{\mathrm{\α}}_{2,\min },\mathrm{...},{\mathrm{\α}}_{s,\min },\mathrm{...},{\mathrm{\α}}_{m,\min }\]}^T\\ {\mathrm{\α}}_{\max }={\[{\mathrm{\α}}_{1,\max },{\mathrm{\α}}_{2,\max },\mathrm{...},{\mathrm{\α}}_{s,\max },\mathrm{...},{\mathrm{\α}}_{m,\max }\]}^T\end{array}---\left(19\right)$

式中：

$\begin{array}{c}{\mathrm{\α}}_{s,\min }=\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s+{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)+\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s-{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)t_{s,\min }\\ {\mathrm{\α}}_{s,\max }=\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s+{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)+\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s-{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)t_{s,\max }\end{array}---\left(20\right)$

步骤(14)、根据α_min和α_max建立其对应的机翼参数化外形，利用高超声速工程算法得到气动力/热响应值Q的区间上下界，可表示为：

$\begin{array}{c}\underset{\‾}{Q}=Q\left({\mathrm{\α}}_{\min }\right)\\ \stackrel{\‾}{Q}=Q\left({\mathrm{\α}}_{\max }\right)\end{array}---\left(21\right)$

式中，Q(α_min)和Q(α_max)分别表示几何特征参数α_min和α_max对应的气动力/热响应值；

步骤(15)、将通过上述方法得到的气动力/热响应值Q的区间上下界与蒙特卡洛方法进行对比，验证所建立方法的有效性。

其中，所述步骤(1)中，机翼平面轮廓参数由表1确定。

表1 机翼平面轮廓参数

其中，所述步骤(2)中，所采用的Clark-Ys标准翼型对应的设计参数由表2确定。

表2 Clark-Ys翼型设计参数

其中，所述步骤(4)中，含不确定性的翼型几何特征参数的区间上下界由表3确定。

表3 翼型几何特征参数的区间界限

其中，所述步骤(7)中，工程算法的计算流程为：首先根据机翼的外形参数建立几何模型，对几何模型进行非结构表面网格划分，然后分别利用切劈法和达黑姆巴克法计算机翼迎风面和背风面的表面压力系数。在此基础上，利用费雷德经验公式计算机翼驻点加热，从而得到机翼表面的最大热流密度、升力系数、阻力系数等三个气动力/热响应值。

其中，所述步骤(15)中，蒙特卡洛方法的计算流程为：在由不确定参数组成的区间范围内，随机选择1000个样本点，逐次计算所有样本点对应的气动力/热响应值，并将响应值中的最大值和最小值作为区间上界和区间下界。

本发明的有益效果是：

本发明利用CST参数化建模方法实现机翼几何外形的快速表征，同时考虑机翼加工过程中由加工误差导致的几何不确定性，建立含区间参数的高超声速机翼气动外形。利用伯恩斯坦多项式近似模型和最小二乘法，建立气动力/热响应值与不确定参数之间的函数关系。通过求解气动力/热响应值关于各个不确定参数的最大/小值点的向量，最终得到气动力/热响应值的区间上界和区间下界，实现考虑几何不确定性的高超声速机翼气动力/热分析。本发明方法得到的区间边界与蒙特卡洛方法得到的区间边界吻合较好，并且可以实现区间包络，为高超声速机翼外形的总体设计提供了新思路。

附图说明

图1为高超声速机翼平面轮廓参数示意图；

图2为Clark-Ys标准翼型示意图；

图3为翼型几何特征参数示意图；

图4为高超声速机翼气动外形示意图；

图5为高超声速机翼表面非结构网格示意图；

图6为机翼表面最大热流密度随马赫数的变化曲线图；

图7为机翼升力系数随马赫数的变化曲线图；

图8为机翼阻力系数随马赫数的变化曲线图；

图9为本发明的方法实现流程图。

具体实施方式

以下将参照附图，对本发明的设计实例进行详细描述。应当理解，所选实例仅为了说明本发明，而不是限制本发明的保护范围。

(1)以类X-37B高超声速飞行器机翼作为研究对象，首先确定高超声速机翼的平面轮廓参数，包括翼根弦长C_r、翼尖弦长C_t、副翼宽度C_a及机翼半展长L。机翼平面轮廓如图1所示，相关参数取值见表4。

表4 机翼平面轮廓参数

(2)以Clark-Ys翼型作为标准翼型(见图2)，利用类型函数/形状函数转换CST方法设定参数化翼型曲线的表达式，其中部分参数的几何含义如图3所示，几何参数由表5给出。

表5 Clark-Ys翼型设计参数

(3)考虑由加工误差导致翼型几何特征参数存在的波动，设定几何特征参数的区间上下界，如表6所示。

表6 翼型几何特征参数的区间界限

(4)利用区间[-1,1]内的n阶伯恩斯坦多项式，构造用于计算多项式系数的样本点，将伯恩斯坦多项式模型中子多项式的极值点t记为：

t＝[t₀,t₁,…,t_r,…,t_n-1,t_n]^T (22)

式中：

$t_r=\frac{2r-n}{n},r=0,1,...,n---\left(23\right)$

本实例中，采用10阶伯恩斯坦多项式，故这里取n＝10。

(5)将式(22)所表示的极值点映射至区间参数空间α^I内的第s个区间参数α_s，则在第s维用于计算该维近似多项式系数的样本点可以表示为：

${\mathrm{\α}}_{s,r}={\[{\mathrm{\α}}_1^c,...,{\mathrm{\α}}_{s-1}^c,{\mathrm{\α}}_{s,r},{\mathrm{\α}}_{s+1}^c,...,{\mathrm{\α}}_m^c\]}^T,r=0,1,...,n---\left(24\right)$

其中，

$\begin{array}{c}{\mathrm{\α}}_i^c=\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_i+{\underset{\‾}{\mathrm{\α}}}_i}{2}\right),i=1,2,\mathrm{...},s-1,s+1,\mathrm{...},m\\ {\mathrm{\α}}_{s,t}=\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s+{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)+\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s+{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)t_r\end{array}---\left(25\right)$

本实例中，共包含3个不确定参数，故这里取m＝3。

(6)根据式(24)中各个样本点α_s,r对应的翼型几何特征参数，利用商业软件CATIA的参数化建模功能建立高超声速机翼的气动外形，如图4所示。将气动外形的几何文件导入商业软件ICEM中进行机翼表面非结构网格的自由划分，如图5所示。得到表面网格节点信息后，在给定工况条件下(如表7所示)，分别利用切劈法和达黑姆巴克法计算机翼迎风面和背风面的表面压力系数。在此基础上，利用费雷德经验公式计算机翼驻点加热情况，从而得到机翼表面的最大热流密度q_max、升力系数C_l、阻力系数C_d等三个气动力/热响应值。同时，将气动力/热响应值以向量Q_s表示为如下形式：

Q_s＝[Q_s,0,...,Q_s,r,...,Q_s,n]^T (26)

式中，Q_s,r＝Q(α_s,r)表示样本点α_s,r对应的气动力/热响应值。

表7 飞行工况

注：表7中，H表示飞行高度，T_∞表示来流温度，ρ_∞表示来流密度，p_∞表示压强，alpha表示攻角，Ma表示飞行马赫数。

(7)根据式(22)中极值点分布，计算得到伯恩斯坦-范德蒙矩阵B，如下所示：

其中，b_n,r(t)为伯恩斯坦多项式的子多项式，记为：

$b_{n,r}\left(t\right)=\frac{n!}{r!(n-r)!}{\left(\frac{t+1}{2}\right)}^r{\left(\frac{-t+1}{2}\right)}^{n-r},r=0,1,...,n---\left(28\right)$

(8)利用最小二乘法计算伯恩斯坦多项式近似模型的待定系数ω_s，公式如下：

ω_s＝[ω_s,0,ω_s,1,...,ω_s,n]^T＝(B^TB)^-1.B^T·Q_s (29)

由此得到气动力/热响应值Q关于第s个不确定参数的伯恩斯坦多项式近似模型P_s(t)，表示如下：

$P_s\left(t\right)=\underset{r=0}{\overset{n}{\Σ}}{\mathrm{\ω}}_{s,r}\·b_{n,r}\left(t\right),t\∈\[-1,1\]---\left(30\right)$

(9)P_s(t)的最值点可以通过其导函数的零点与自变量的端点产生，其中导函数零点可通过下式计算：

$\frac{{\mathrm{dP}}_s\left(t\right)}{dt}=0---\left(31\right)$

将零点记为t_s＝[t_s,1,t_s,2,...,t_s,n]，根据P_s(t)的自变量范围，可将零点进一步修改为：

其中Im和Re分别表示对应变量的虚部和实部。

(10)气动力/热响应值Q关于第s个区间参数的最小值点和最大值点可计算为：

$\begin{array}{c}{P}_s\left(t_{s,\min }\right)=\underset{t\∈T_s}{\min }\left\{P_s\left(t\right)\right\}\\ P_s\left(t_{s,\max }\right)=\underset{t\∈T_s}{\min }\left\{P_s\left(t\right)\right\}\end{array}---\left(33\right)$

式中，T_s＝[t_s,-1,1]。

(11)重复步骤(5)～(10)，以相同的方式计算气动力/热响应值Q关于所有区间参数的最小值点和最大值点，最终可以得到气动力/热响应值在标准区间[-1,1]内的最小值点和最大值点，即：

$\begin{array}{c}{t}_{\min }={\[t_{1,\min },t_{2,\min },\mathrm{...},t_{s,\min },\mathrm{...},t_{m,\min }\]}^T\\ t_{\max }={\[t_{1,\max },t_{2,\max },\mathrm{...},t_{s,\max },\mathrm{...},t_{m,\max }\]}^T\end{array}---\left(34\right)$

(12)将步骤(11)中的标准区间内的最小值点和最大值点映射到实际参数空间有：

$\begin{array}{c}{\mathrm{\α}}_{\min }={\[{\mathrm{\α}}_{1,\min },{\mathrm{\α}}_{2,\min },\mathrm{...},{\mathrm{\α}}_{s,\min },\mathrm{...},{\mathrm{\α}}_{m,\min }\]}^T\\ {\mathrm{\α}}_{\max }={\[{\mathrm{\α}}_{1,\max },{\mathrm{\α}}_{2,\max },\mathrm{...},{\mathrm{\α}}_{s,\max },\mathrm{...},{\mathrm{\α}}_{m,\max }\]}^T\end{array}---\left(35\right)$

式中：

$\begin{array}{c}{\mathrm{\α}}_{s,\min }=\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s+{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)+\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s-{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)t_{s,\min }\\ {\mathrm{\α}}_{s,\max }=\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s+{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)+\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s-{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)t_{s,\max }\end{array}---\left(36\right)$

(13)根据α_min和α_max建立其对应的机翼参数化外形，利用高超声速工程算法得到气动力/热响应值Q的区间上下界，表示为：

$\begin{array}{c}\underset{\‾}{Q}=Q\left({\mathrm{\α}}_{\min }\right)\\ \stackrel{\‾}{Q}=Q\left({\mathrm{\α}}_{\max }\right)\end{array}---\left(37\right)$

式中，Q(α_min)和Q(α_max)分别表示几何参数α_min和α_max对应的气动力/热响应值。

(14)将通过上述方法得到的气动力/热响应值Q的区间上下界与蒙特卡洛方法进行对比，验证所建立方法的有效性。蒙特卡洛方法的计算流程为：在由不确定参数组成的区间范围内，随机选择1000个样本点，逐次计算所有样本点对应的气动力/热响应值，并将响应值中的最大值和最小值作为区间上界和区间下界。

(15)将本发明方法发生得到的机翼表面最大热流密度q_max、升力系数C_l及阻力系数C_d的区间边界与通过蒙特卡洛方法得到的区间边界进行对比，计算结果如表8-10及图6-8所示。从结果对比中可以看出，本发明方法得到的区间边界与蒙特卡洛方法得到的区间边界吻合较好，并且可以包络蒙特卡洛计算结果，从而验证了本发明方法的有效性及可行性。

表8 机翼表面最大热流密度q_max(Kw/m²)的对比

表9 机翼升力系数C_l的对比

表10 机翼升力系数C_d的对比

综上所述，本发明提出了一种考虑几何不确定性的高超声速机翼气动力/热分析方法。该方法充分考虑外形设计中存在的不确定因素，在不确定参数概率密度未知的情况下，利用区间向量实现不确定参数的定量化表征。在给定飞行工况下，根据伯恩斯坦多项式近似模型构造不确定参数的样本点，利用气动力/热工程算法计算样本点处的气动力/热响应值，同时基于最小二乘策略得到多项式近似模型的待定系数。在给定区间范围内，利用已建立的多项式近似模型求解气动力/热响应值关于各个不确定参数的极大/小值点，从而组合形成最大/小值点的向量，并最终得到气动力/热响应值的区间上界和区间下界，实现考虑几何不确定性的高超声速机翼气动力/热分析。与蒙特卡洛方法相比，本发明方法得到的气动力/热响应区间边界与蒙特卡洛方法得到的区间边界误差较小，并且可以包络蒙特卡洛计算结果，从而验证了本方法的有效性及合理性，为高超声速机翼外形的总体设计提供了新思路。

以上仅是本发明的具体步骤，对本发明的保护范围不构成任何限制，其可扩展应用于高超声速飞行器外形设计领域，凡采用等同转换或者等效替换而形成的技术方案，均落在本发明权利保护范围之内。

Claims (6)

1.一种考虑几何不确定性的高超声速机翼气动力/热分析方法，其特征在于实现步骤如下：

步骤(1)、确定高超声速机翼的平面轮廓参数，包括翼根弦长C_r、翼尖弦长C_t、副翼宽度C_a及机翼半展长L；

步骤(2)、利用类函数/形函数转换方法建立高超声速机翼翼型的参数化表达式，翼型的几何曲线可用下列函数表示：

$\frac{y}{c}\left(\frac{x}{c}\right)=C\left(\frac{x}{c}\right)S\left(\frac{x}{c}\right)+\frac{x}{c}\frac{z_{te}}{c}---\left(1\right)$

式中，x/c为翼型弦向的无量纲坐标值，y/c为翼型法向的无量纲坐标值，C(x/c)和S(x/c)分别为类型函数和形状函数，z_te/c为翼型后缘点的无量纲坐标值，C(x/c)可表示为：

$C\left(\frac{x}{c}\right)={\left(\frac{x}{c}\right)}^{N_1}{(1-\frac{x}{c})}^{N_2},0\≤\frac{x}{c}\≤1---\left(2\right)$

对于Clark-Ys翼型，取指数N₁＝0.5,N₂＝1，S(x/c)可表示为：

$S\left(\frac{x}{c}\right)=\underset{i=0}{\overset{n}{\Σ}}\[{\mathrm{\λ}}_i\·\frac{n!}{i!(n-i)!}\·{\left(\frac{x}{c}\right)}^i\·{(1-\frac{x}{c})}^{n-i}\],0\≤\frac{x}{c}\≤1---\left(3\right)$

S(x/c)与翼型前缘半径R_le/c及后缘倾角β满足以下关系：

$\left\{\begin{array}{c}S\left(0\right)={\mathrm{\λ}}_0=\sqrt{\frac{2R_{le}}{c}},\\ S\left(1\right)={\mathrm{\λ}}_n=tan\mathrm{\β}+\frac{z_{te}}{c}\end{array}\right.---\left(4\right)$

在形状函数S(x/c)中取n＝3，这样共包含等8个几何设计参数，其中，表示翼型前缘半径，β₁和β₂表示上、下翼面的后缘倾角，z_te/c为翼型后缘点的无量纲坐标值，λ₁、λ₂、λ₁'、λ₂'为上、下翼面形状函数多项式的加权系数；

步骤(3)、考虑由加工误差导致翼型几何特征参数存在的波动，将存在不确定性的几何特征参数记为其余参数均视为确定性参数；

步骤(4)、设定翼型几何特征参数α的区间上下界，可表示为：

${\mathrm{\α}}^I=\[\underset{\‾}{\mathrm{\α}},\stackrel{\‾}{\mathrm{\α}}\]---\left(5\right)$

式中，为高超声速机翼翼型几何特征参数的区间下界，为区间上界；

步骤(5)、利用区间[-1,1]内的n阶伯恩斯坦多项式，构造用于计算多项式系数的样本点，将伯恩斯坦多项式模型中子多项式的极值点t记为：

t＝[t₀,t₁,…,t_r,…,t_n-1,t_n]^T (6)

式中：

$t_r=\frac{2r-n}{n},r=0,1,...,n---\left(7\right)$

步骤(6)、将式(6)所表示的极值点映射至区间参数空间α^I内的第s个区间参数α_s，则在第s维用于计算该维近似多项式系数的样本点可以表示为：

${\mathrm{\α}}_{s,r}={\[{\mathrm{\α}}_1^c,...,{\mathrm{\α}}_{s-1}^c,{\mathrm{\α}}_{s,r},{\mathrm{\α}}_{s+1}^c,...,{\mathrm{\α}}_m^c\]}^T,r=0,1,...,n---\left(8\right)$

其中，

$\begin{array}{c}{\mathrm{\α}}_i^c=\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_i+{\underset{\‾}{\mathrm{\α}}}_i}{2}\right),i=1,2,...,s-1,s+1,...,m\\ {\mathrm{\α}}_{s,r}=\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s+{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)+\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s-{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)t_r\end{array}---\left(9\right)$

式中，m为不确定参数的个数；

步骤(7)、根据式(8)中样本点对应翼型的几何特征参数，建立高超声速机翼的气动外形，进而利用高超声速工程算法计算气动力/热响应值，记为Q_s，表示为如下形式：

Q_s＝[Q_s,0,...,Q_s,r,...,Q_s,n]^T (10)

式中，Q_s,r＝Q(α_s,r)表示样本点α_s,r对应的气动力/热响应值；

步骤(8)、利用最小二乘法计算伯恩斯坦多项式近似模型的待定系数ω_s，公式如下：

ω_s＝[ω_s,0,ω_s,1,...,ω_s,n]^T＝(B^TB)^-1·B^T·Q_s (11)

式中，B为伯恩斯坦-范德蒙矩阵，可表示为：

其中，b_n,r(t)为伯恩斯坦多项式的子多项式，记为：

$b_{n,r}\left(t\right)=\frac{n!}{r!(n-r)!}{\left(\frac{t+1}{2}\right)}^r{\left(\frac{-t+1}{2}\right)}^{n-r},r=0,1,...,n---\left(13\right)$

步骤(9)、根据步骤(8)中得到的系数ω_s，建立气动力/热响应值Q关于第s个不确定参数的伯恩斯坦多项式近似模型P_s(t)，表示如下：

$P_s\left(t\right)=\underset{r=0}{\overset{n}{\Σ}}{\mathrm{\ω}}_{s,r}\·b_{n,r}\left(t\right),t\∈\[-1,1\]---\left(14\right)$

步骤(10)、步骤(9)中P_s(t)的最值点可以通过其导函数的零点与自变量的端点产生，其中导函数零点可通过下式计算：

$\frac{{\mathrm{dP}}_s\left(t\right)}{dt}=0---\left(15\right)$

将零点记为t_s＝[t_s,1,t_s,2,...,t_s,n]，根据P_s(t)的自变量范围，可将零点进一步修改为：

其中Im和Re分别表示对应变量的虚部和实部；

步骤(11)、气动力/热响应值Q关于第s个区间参数的最小值点和最大值点可计算为：

$\begin{array}{c}{P}_s\left(t_{s,\min }\right)=\underset{t\∈T_s}{min}\left\{P_s\left(t\right)\right\}\\ P_s\left(t_{s,\max }\right)=\underset{t\∈T_s}{max}\left\{P_s\left(t\right)\right\}\end{array}---\left(17\right)$

式中，T_s＝[t_s,-1,1]；

步骤(12)、重复步骤(6)～(11)，以相同的方式计算气动力/热响应值Q关于所有区间参数的最小值点和最大值点，最终可以得到气动力/热响应值在标准区间[-1,1]内的最小值点和最大值点，即：

$\begin{array}{c}{t}_{\min }={\[t_{1,\min },t_{2,\min },...,t_{s,\min },...,t_{m,\min }\]}^T\\ t_{\max }={\[t_{1,\max },t_{2,\max },...,t_{s,\max },...,t_{m,\max }\]}^T\end{array}---\left(18\right)$

步骤(13)、将步骤(12)中标准区间内的最小值点和最大值点映射到实际参数空间有：

$\begin{array}{c}{\mathrm{\α}}_{\min }={\[{\mathrm{\α}}_{1,\min },{\mathrm{\α}}_{2,\min },...,{\mathrm{\α}}_{s,\min },...,{\mathrm{\α}}_{m,\min }\]}^T\\ {\mathrm{\α}}_{\max }={\[{\mathrm{\α}}_{1,\max },{\mathrm{\α}}_{2,\max },...,{\mathrm{\α}}_{s,\max },...,{\mathrm{\α}}_{m,\max }\]}^T\end{array}---\left(19\right)$

式中：

$\begin{array}{c}{\mathrm{\α}}_{s,\min }=\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s+{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)+\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s-{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)t_{s,\min }\\ {\mathrm{\α}}_{s,\max }=\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s+{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)+\left(\frac{{\stackrel{\‾}{\mathrm{\α}}}_s-{\underset{\‾}{\mathrm{\α}}}_s}{2}\right)t_{s,\max }\end{array}---\left(20\right)$

步骤(14)、根据α_min和α_max建立其对应的机翼参数化外形，利用高超声速工程算法得到气动力/热响应值Q的区间上下界，可表示为：

$\begin{array}{c}\underset{\‾}{Q}=Q\left({\mathrm{\α}}_{min}\right)\\ \stackrel{\‾}{Q}=Q\left({\mathrm{\α}}_{max}\right)\end{array}---\left(21\right)$

式中，Q(α_min)和Q(α_max)分别表示几何特征参数α_min和α_max对应的气动力/热响应值；

步骤(15)、将通过上述方法得到的气动力/热响应值Q的区间上下界与蒙特卡洛方法进行对比，验证所建立方法的有效性。

2.根据权利要求1所述的一种考虑几何不确定性的高超声速机翼气动力/热分析方法，其特征在于：所述步骤(1)中，

表1 机翼平面轮廓参数

机翼平面轮廓参数由表1确定。

3.根据权利要求1所述的一种考虑几何不确定性的高超声速机翼气动力/热分析方法，其特征在于：所述步骤(2)中，

表2 Clark-Ys翼型设计参数

所采用的Clark-Ys标准翼型对应的设计参数由表2确定。

4.根据权利要求1所述的一种考虑几何不确定性的高超声速机翼气动力/热分析方法，其特征在于：所述步骤(4)中，

表3 翼型几何特征参数的区间界限

含不确定性的翼型几何特征参数的区间上下界由表3确定。

5.根据权利要求1所述的一种考虑几何不确定性的高超声速机翼气动力/热分析方法，其特征在于：所述步骤(7)中，工程算法的计算流程为：首先根据机翼的几何参数建立气动外形，对几何外形进行非结构网格自由划分，然后分别利用切劈法和达黑姆巴克法计算机翼迎风面和背风面的表面压力系数，在此基础上，利用费雷德经验公式计算机翼驻点加热，从而得到机翼表面的最大热流密度、升力系数、阻力系数等三个气动力/热响应值。

6.根据权利要求1所述的一种考虑几何不确定性的高超声速机翼气动力/热分析方法，其特征在于：所述步骤(15)中，蒙特卡洛方法的计算原理为：在由不确定参数组成的区间范围内，随机选择1000个样本点，逐次计算所有样本点对应的气动力/热响应值，并将响应值中的最大值和最小值作为区间上界和区间下界。