一种用于飞行器定常绕流数值求解的流场初始化方法

摘要

本发明属于计算流体力学领域，具体涉及一种用于飞行器定常绕流数值求解的流场初始化方法。本发明基于已知飞行状态的绕流定常解，完成飞行器定常绕流数值求解的流场初始化，包括以下步骤：计算待初始化状态来流条件下流场压强初值；计算待初始化状态来流条件下流场马赫数初值并修正流场压强初值；计算待初始化状态来流条件下流场温度初值与速度初值；计算待初始化状态来流条件下流场速度矢量初值；插值得到待初始化状态对应网格上流场压强、温度、速度矢量初值。本发明提出的方法计算量小、通用性强，显著提高了流场初始化质量，大幅提高多状态飞行器定常绕流数值求解效率。

说明书

技术领域

本发明属于计算流体力学领域，具体涉及一种用于飞行器定常绕流数值求解的流场初始化方法。

背景技术

飞行器定常绕流数值求解是计算流体力学领域的重要研究方向，能用极少假设给出详尽、准确、广泛的飞行器绕流信息，是获取飞行器绕流特性、气动力/热性能的重要手段，对飞行器气动力/热性能仿真分析与气动设计具有重要意义。流场初始化的核心任务为给定流场数值仿真的初始条件，仿真初始流场的好坏对定常流场数值仿真的数值稳定性、收敛性具有重要影响，好的初始流场可使数值求解在短时间内收敛，而不好的初始流场会导致数值求解收敛速度慢，甚至导致发散。一般而言，初始流场与流场的定常解越接近，求解的数值稳定性、收敛性越好。目前飞行器定常绕流数值求解的流场初始化方法主要有以下几类：(1)根据边界条件给为定值；(2)根据边界条件线性插值；(3)基于拉普拉斯方程数值求解的混合初始化方法。前两类初始化方法计算量小，但难以给出理想的流场初值，第三类方法给出的流场初值优于前两者，但计算量较大。

现有流场初始化方法皆针对特定来流条件单独进行初始化，而飞行器绕流特性、气动力/热性能等的数值求解工作往往需要针对多个不同来流状态(飞行器外形不变)。现有飞行器定常绕流数值求解的流场初始化方法未能利用已知的其它飞行状态的仿真结果，影响了多状态飞行器定常绕流数值求解效率的进一步提高。

发明内容

为进一步提高多状态飞行器定常绕流数值求解效率，本发明提出基于已知状态流场定常解的飞行器定常绕流数值求解的流场初始化方法。流场初始化的主要任务为给出所有计算网格节点的压强、温度、速度矢量。设已知状态流场定常解对应网格节点集合为G₁＝{(x₁，y₁，z₁)}，来流条件为“来流条件I”，对应马赫数、压强、温度、攻角、侧滑角依次为M_∞1、p_∞1、T_∞1、α₁、β₁，流场定常解为“流场I”，其网格节点上的压强、温度、速度矢量依次为p₁、T₁、待初始化的计算状态对应网格节点集合为G₂＝{(x₂，y₂，z₂)}，来流条件为“来流条件II”，对应马赫数、压强、温度、攻角、侧滑角依次为M_∞2、p_∞2、T_∞2、α₂、β₂。本发明根据G₁、M_∞1、p_∞1、T_∞1、α₁、β₁、p₁、T₁、给出来流条件II下的飞行器定常绕流数值求解问题的初始条件，即计算网格节点G₂上的压强、温度、速度矢量初值p₂、T₂、

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

一种用于飞行器定常绕流数值求解的流场初始化方法，具体包括以下步骤：

第一步：计算来流条件II下流场压强初值p₂₁；

流场中压强系数C_p的计算公式为：

$C_p=\frac{p-p_{\infty }}{\frac{1}{2}{\rho }_{\infty }{V}_{\infty }^2}---\left(1\right)$

式中，p为压强；p_∞、ρ_∞、V_∞依次为来流压强、密度、速度。

超声速条件下，驻点前存在正激波，正激波后压强满足如下关系：

$\frac{p_s-p_{\infty }}{p_{\infty }}=\frac{\gamma }{\gamma +1}(M_{\infty }^2-1)---\left(2\right)$

式中，γ为来流比热比，p_s为正激波后压强，M_∞为来流马赫数。

综合式(1)、(2)推导得正激波后压强系数C_ps：

$C_{ps}=C_{ps}\left(M_{\infty }\right)=\frac{4(M_{\infty }^2-1)}{M_{\infty }^2(\gamma +1)}---\left(3\right)$

按式(1)计算得来流条件I下流场的压强系数C_p1：

$C_{p1}=\frac{p_1-p_{\infty 1}}{\frac{1}{2}{\rho }_{\infty 1}{V}_{\infty 1}^2}---\left(4\right)$

式中，p₁为来流条件I下流场的压强；V_∞1为来流条件I对应的来流速度，按下式计算：

$V_{\infty 1}=M_{\infty 1}\sqrt{{\mathrm{\gamma RT}}_{\infty 1}}---\left(5\right)$

式中，γ为来流比热比，R为来流气体常数。

分两种情况计算来流条件II下流场压强系数初值C_p2，

(a)当M_∞1＞1且M_∞2＞1时，按下式：

$C_{p2}=\left(\begin{array}{cc}{C}_{p1}& C_{p1}=0\\ C_{p1}\frac{C_{ps}\left(M_{\infty 2}\right)}{C_{ps}\left(M_{\infty 1}\right)}& C_{p1}\ne 0\end{array}\right);---\left(6\right)$

(b)当M_∞1≤1或M_∞2≤1时，按下式：

C_p2＝C_p1；(7)

得来流条件II下压强初值p₂₁计算公式：

$p_{21}=\frac{1}{2}{C}_{p2}{\rho }_{\infty 2}{V}_{\infty 2}^2+p_{\infty 2};---\left(8\right)$

式中，V_∞2为来流条件II对应的来流速度，按下式计算：

$V_{\infty 2}=M_{\infty 2}\sqrt{{\mathrm{\gamma RT}}_{\infty 2}};---\left(9\right)$

第二步：计算来流条件II下流场马赫数初值M₂₁并修正p₂₁；

设流场总压比为流场总压与来流总压的比值。正激波后总压比r_p0为：

$r_{p0}=\frac{{\left(\frac{\frac{\gamma +1}{2}{M}_{\infty }^2}{1+\frac{\gamma -1}{2}{M}_{\infty }^2}\right)}^{\frac{\gamma }{\gamma -1}}}{{(\frac{2\gamma }{\gamma +1}{M}_{\infty }^2-\frac{\gamma -1}{\gamma +1})}^{\frac{1}{\gamma -1}}}---\left(10\right)$

流场中总压p₀、压强p、马赫数M具有以下关系：

$p_0=p_0(p,M)=p{(1+\frac{\gamma -1}{2}{M}^2)}^{\frac{\gamma }{\gamma -1}}---\left(11\right)$

来流条件I下流场总压比r_p01按下式计算：

$r_{p01}=\frac{p_0(p_1,M_1)}{p_0(p_{\infty },M_{\infty 1})}---\left(12\right)$

式中，M₁为来流条件I下流场的马赫数，按下式计算：

$M_1=V_1/\sqrt{{\mathrm{\gamma RT}}_1}---\left(13\right)$

式中，V₁为来流条件I下流场的速度，T₁为来流条件I下流场的温度。

按下式计算来流条件II下流场总压比初值r_p02：

$r_{p02}=\left(\begin{array}{cc}{r}_{p01}& r_{p01}\ge \kappa \\ r_{p01}\frac{r_{p0}\left(M_{\infty 2}\right)}{r_{p0}\left(M_{\infty 1}\right)}& r_{p01}<\kappa \end{array}\right)---\left(14\right)$

式中，κ≤1，一般可取为0.8～0.99。

根据式(14)、(11)推导得来流条件II下流场马赫数初值M₂₁计算公式：

$\left(\begin{array}{c}g=\frac{2}{\gamma -1}[{\left(r_{p02}\frac{p_{02}}{p_{21}}\right)}^{\frac{\gamma -1}{\gamma }}-1]\\ M_{21}=\left(\begin{array}{cc}0& g<0\\ \sqrt{g}& g\ge 0\end{array}\right)\end{array}\right)---\left(15\right)$

式中，g为计算M₂₁的中间变量，为来流条件II对应的总压。

当式(15)中g＜0，第一步得到的p₂₁不再满足式(11)，按式(11)推导重新计算p₂₁：

$p_{21}=p_{02}/{(1+\frac{\gamma -1}{2}{M}_{21}^2)}^{\frac{\gamma }{\gamma -1}}=p_{02};---\left(16\right)$

第三步：计算来流条件II下流场温度初值T₂₁与速度初值V₂₁；

按下式计算来流条件II下流场温度初值T₂₁：

$T_{21}=T_{02}/(1+\frac{\gamma -1}{2}{M}_{21}^2)---\left(17\right)$

式中，为来流条件II对应的总温。

按下式计算来流条件II下流场温度初值V₂₁：

$V_{21}=M_{21}\sqrt{{\mathrm{\gamma RT}}_{21}}---\left(18\right)$

式中，R为来流气体常数。

第四步：计算来流条件II下流场速度矢量初值

按下式计算来流条件II下流场速度矢量初值

${\overrightarrow{V}}_{21}=\frac{V_{21}}{V_1}{M}_3^{\alpha }{M}_2^{\beta }{\overrightarrow{V}}_1---\left(19\right)$

式中：

$M_3^{\alpha }=\left(\begin{array}{ccc}cos({\alpha }_1-{\alpha }_2)& sin({\alpha }_1-{\alpha }_2)& 0\\ -sin({\alpha }_1-{\alpha }_2)& cos({\alpha }_1-{\alpha }_2)& 0\\ 0& 0& 1\end{array}\right),M_2^{\beta }=\left(\begin{array}{ccc}1& 0& 0\\ 0& cos({\beta }_2-{\beta }_1)& \sin ({\beta }_2-{\beta }_1)\\ 0& -\sin ({\beta }_2-{\beta }_1)& cos({\beta }_2-{\beta }_1)\end{array}\right)---\left(20\right)$

第五步：插值得到待初始化状态对应网格G₂上流场压强、温度、速度矢量初值p₂、T₂、

前四步得到的压强初值p₂₁、温度初值T₂₁、速度矢量初值对应网格节点为G₁，若G₁与G₂一致，则p₂₁、T₂₁、的集合即为压强、温度、速度矢量初值p₂、T₂、若网格G₁与网格G₂不一致，则将网格G₁上的压强初值p₂₁、温度初值T₂₁、速度矢量初值插值至网格G₂，插值方法取计算量较小的最近邻点插值法，得到压强、温度、速度矢量初值p₂、T₂、

本发明的有益效果是：

(1)基于已知飞行状态的绕流定常解，完成飞行器定常绕流数值求解的流场初始化，提高了流场初始化质量，大幅提高多状态飞行器定常绕流数值求解效率；

(2)本发明提出的流场初始化方法仅需代数运算，计算量小，耗时短，效率高；

(3)本发明提出的流场初始化方法通用性强，在来流为亚、跨、超声速时皆可用，且适用于已知流场状态网格与待初始化状态网格不一致的情形。

附图说明

图1为本发明提出的用于飞行器定常绕流数值求解的流场初始化方法的流程图；

图2为标准模型“AGARD HB-2”外形；

图3为已知飞行状态定常流场对称面的压强分布(单位Pa)；

图4为已知飞行状态定常流场对称面的温度分布(单位K)；

图5为已知飞行状态定常流场对称面的速度分布(单位m/s)；

图6为已知飞行状态定常流场对应网格(对称面)；

图7为本发明方法得到的初始流场对称面的压强分布(单位Pa)；

图8为本发明方法得到的初始流场对称面的温度分布(单位K)；

图9为本发明方法得到的初始流场对称面的速度分布(单位m/s)；

图10为待初始化状态对应的网格(对称面)；

图11为不同求解过程的连续性残差在计算迭代过程的变化情况；

图12为不同求解过程的轴向力系数在计算迭代过程的变化情况。

具体实施方式

本发提出的用于飞行器定常绕流数值求解的流场初始化方法的流程如图1所示。下面结合附图，对本发明的具体实施方式作进一步的说明。

步骤一：计算来流条件II下流场压强初值p₂₁：

以标准模型“AGARD HB-2”(参考文献：高超声速风洞气动力试验方法[标准].GJB4399-2002.2002)为例，其外形如图2所示，来流马赫数M_∞1＝5、压强p_∞1＝6500Pa、温度T_∞1＝288.15K、攻角α₁＝0°、侧滑角β₁＝0°状态下其绕流定常解在对称面的压强p₁、温度T₁、速度V₁分度依次如图3～图5所示，对应网格G₁的对称面网格如图6所示。待初始化的计算状态对应马赫数M_∞2、压强p_∞2、温度T_∞2、攻角α₂、侧滑角β₂依次为7、3000Pa、288.15K、10°、0°。

按发明内容中第一步描述方法，得来流条件II下流场压强初值p₂₁。

步骤二：计算来流条件II下流场马赫数初值M₂₁并修正p₂₁：

按发明内容中第二步描述方法，得来流条件II下流场马赫数初值M₂₁。根据M₂₁修正后的流场压强初值p₂₁如图7所示。

步骤三：计算来流条件II下流场温度初值T₂₁与速度初值V₂₁：

按发明内容中第三步描述方法，得来流条件II下流场温度初值T₂₁与速度初值V₂₁分布如图8、图9所示。

步骤四：计算来流条件II下流场速度矢量初值

按发明内容中第四步描述方法，得来流条件II下流场速度矢量初值

步骤五：插值得到待初始化状态对应网格G₂上流场压强、温度、速度矢量的初值p₂、T₂、

待初始化状态对应网格G₂的对称面网格如图10所示，前四步得到的p₂₁、T₂₁、对应网格节点为G₁，由于G₂与G₁不一致，将G₁上的p₂₁、T₂₁、插值至G₂，得到待求的压强、温度、速度矢量初值p₂、T₂、

为分析本发明提出方法的特点，取本发明得到的初始流场与常用的基于来流参数的常值初始流场，分别进行AGARD HB-2定常绕流仿真，仿真过程库朗数取0.75与2.5两个值，比较二者的计算稳定性与收敛速度。实施的四个仿真过程依次为：

仿真过程I：基于本发明初始化方法，库朗数取2.5；

仿真过程II：基于本发明初始化方法，库朗数取0.75；

仿真过程III：基于常值初始化方法，库朗数取2.5；

仿真过程IV：基于常值初始化方法，库朗数取0.75。

不同仿真过程的连续性残差、轴向力系数在计算迭代过程的变化情况分别如图11～图12所示。由图11中连续性残差曲线得：四个仿真过程中，仿真过程III在迭代至第72步时，残差激增，计算发散，其余三个仿真过程皆收敛；仿真过程I与仿真过程II的残差数值小于且残差下降速度快于仿真过程III与仿真过程IV。由图12中轴向力系数曲线得：仿真过程I在迭代至约240步时轴向力系数收敛；仿真过程II在迭代至约450步时轴向力系数收敛；仿真过程IV在迭代至约835步时轴向力系数收敛；仿真过程I的收敛速率约为仿真过程IV的3.5倍。分析表明，本发明提出的飞行器定常绕流数值求解的流场初始化在计算稳定性、收敛速度等方面皆明显优于常规常值初始化方法。

本发明提出了基于已知飞行状态的绕流定常解的高效、通用的飞行器定常绕流数值求解的流场初始化方法，提高了流场初始化质量，改善了流场数值求解的计算稳定性，加快了收敛速度，大幅提高了多状态飞行器定常绕流数值求解效率，是多状态飞行器定常绕流数值求解的理想初始化方法。

综上所述，虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明，任何本领域普通技术人员，在不脱离本发明的精神和范围内，当可作各种更动与润饰，因此本发明的保护范围当视权利要求书界定的范围为准。