CN106339534A - 一种用于飞行器定常绕流数值求解的流场初始化方法 - Google Patents
一种用于飞行器定常绕流数值求解的流场初始化方法 Download PDFInfo
- Publication number
- CN106339534A CN106339534A CN201610687769.4A CN201610687769A CN106339534A CN 106339534 A CN106339534 A CN 106339534A CN 201610687769 A CN201610687769 A CN 201610687769A CN 106339534 A CN106339534 A CN 106339534A
- Authority
- CN
- China
- Prior art keywords
- infin
- initial value
- gamma
- pressure
- inlet flow
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Aerodynamic Tests, Hydrodynamic Tests, Wind Tunnels, And Water Tanks (AREA)
Abstract
本发明属于计算流体力学领域,具体涉及一种用于飞行器定常绕流数值求解的流场初始化方法。本发明基于已知飞行状态的绕流定常解,完成飞行器定常绕流数值求解的流场初始化,包括以下步骤:计算待初始化状态来流条件下流场压强初值;计算待初始化状态来流条件下流场马赫数初值并修正流场压强初值;计算待初始化状态来流条件下流场温度初值与速度初值;计算待初始化状态来流条件下流场速度矢量初值;插值得到待初始化状态对应网格上流场压强、温度、速度矢量初值。本发明提出的方法计算量小、通用性强,显著提高了流场初始化质量,大幅提高多状态飞行器定常绕流数值求解效率。
Description
技术领域
本发明属于计算流体力学领域,具体涉及一种用于飞行器定常绕流数值求解的流场初始化方法。
背景技术
飞行器定常绕流数值求解是计算流体力学领域的重要研究方向,能用极少假设给出详尽、准确、广泛的飞行器绕流信息,是获取飞行器绕流特性、气动力/热性能的重要手段,对飞行器气动力/热性能仿真分析与气动设计具有重要意义。流场初始化的核心任务为给定流场数值仿真的初始条件,仿真初始流场的好坏对定常流场数值仿真的数值稳定性、收敛性具有重要影响,好的初始流场可使数值求解在短时间内收敛,而不好的初始流场会导致数值求解收敛速度慢,甚至导致发散。一般而言,初始流场与流场的定常解越接近,求解的数值稳定性、收敛性越好。目前飞行器定常绕流数值求解的流场初始化方法主要有以下几类:(1)根据边界条件给为定值;(2)根据边界条件线性插值;(3)基于拉普拉斯方程数值求解的混合初始化方法。前两类初始化方法计算量小,但难以给出理想的流场初值,第三类方法给出的流场初值优于前两者,但计算量较大。
现有流场初始化方法皆针对特定来流条件单独进行初始化,而飞行器绕流特性、气动力/热性能等的数值求解工作往往需要针对多个不同来流状态(飞行器外形不变)。现有飞行器定常绕流数值求解的流场初始化方法未能利用已知的其它飞行状态的仿真结果,影响了多状态飞行器定常绕流数值求解效率的进一步提高。
发明内容
为进一步提高多状态飞行器定常绕流数值求解效率,本发明提出基于已知状态流场定常解的飞行器定常绕流数值求解的流场初始化方法。流场初始化的主要任务为给出所有计算网格节点的压强、温度、速度矢量。设已知状态流场定常解对应网格节点集合为G1={(x1,y1,z1)},来流条件为“来流条件I”,对应马赫数、压强、温度、攻角、侧滑角依次为M∞1、p∞1、T∞1、α1、β1,流场定常解为“流场I”,其网格节点上的压强、温度、速度矢量依次为p1、T1、待初始化的计算状态对应网格节点集合为G2={(x2,y2,z2)},来流条件为“来流条件II”,对应马赫数、压强、温度、攻角、侧滑角依次为M∞2、p∞2、T∞2、α2、β2。本发明根据G1、M∞1、p∞1、T∞1、α1、β1、p1、T1、给出来流条件II下的飞行器定常绕流数值求解问题的初始条件,即计算网格节点G2上的压强、温度、速度矢量初值p2、T2、
本发明采用的技术方案是:
一种用于飞行器定常绕流数值求解的流场初始化方法,具体包括以下步骤:
第一步:计算来流条件II下流场压强初值p21;
流场中压强系数Cp的计算公式为:
式中,p为压强;p∞、ρ∞、V∞依次为来流压强、密度、速度。
超声速条件下,驻点前存在正激波,正激波后压强满足如下关系:
式中,γ为来流比热比,ps为正激波后压强,M∞为来流马赫数。
综合式(1)、(2)推导得正激波后压强系数Cps:
按式(1)计算得来流条件I下流场的压强系数Cp1:
式中,p1为来流条件I下流场的压强;V∞1为来流条件I对应的来流速度,按下式计算:
式中,γ为来流比热比,R为来流气体常数。
分两种情况计算来流条件II下流场压强系数初值Cp2,
(a)当M∞1>1且M∞2>1时,按下式:
(b)当M∞1≤1或M∞2≤1时,按下式:
Cp2=Cp1;(7)
得来流条件II下压强初值p21计算公式:
式中,V∞2为来流条件II对应的来流速度,按下式计算:
第二步:计算来流条件II下流场马赫数初值M21并修正p21;
设流场总压比为流场总压与来流总压的比值。正激波后总压比rp0为:
流场中总压p0、压强p、马赫数M具有以下关系:
来流条件I下流场总压比rp01按下式计算:
式中,M1为来流条件I下流场的马赫数,按下式计算:
式中,V1为来流条件I下流场的速度,T1为来流条件I下流场的温度。
按下式计算来流条件II下流场总压比初值rp02:
式中,κ≤1,一般可取为0.8~0.99。
根据式(14)、(11)推导得来流条件II下流场马赫数初值M21计算公式:
式中,g为计算M21的中间变量,为来流条件II对应的总压。
当式(15)中g<0,第一步得到的p21不再满足式(11),按式(11)推导重新计算p21:
第三步:计算来流条件II下流场温度初值T21与速度初值V21;
按下式计算来流条件II下流场温度初值T21:
式中,为来流条件II对应的总温。
按下式计算来流条件II下流场温度初值V21:
式中,R为来流气体常数。
第四步:计算来流条件II下流场速度矢量初值
按下式计算来流条件II下流场速度矢量初值
式中:
第五步:插值得到待初始化状态对应网格G2上流场压强、温度、速度矢量初值p2、T2、
前四步得到的压强初值p21、温度初值T21、速度矢量初值对应网格节点为G1,若G1与G2一致,则p21、T21、的集合即为压强、温度、速度矢量初值p2、T2、若网格G1与网格G2不一致,则将网格G1上的压强初值p21、温度初值T21、速度矢量初值插值至网格G2,插值方法取计算量较小的最近邻点插值法,得到压强、温度、速度矢量初值p2、T2、
本发明的有益效果是:
(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下流场压强初值p21:
以标准模型“AGARD HB-2”(参考文献:高超声速风洞气动力试验方法[标准].GJB4399-2002.2002)为例,其外形如图2所示,来流马赫数M∞1=5、压强p∞1=6500Pa、温度T∞1=288.15K、攻角α1=0°、侧滑角β1=0°状态下其绕流定常解在对称面的压强p1、温度T1、速度V1分度依次如图3~图5所示,对应网格G1的对称面网格如图6所示。待初始化的计算状态对应马赫数M∞2、压强p∞2、温度T∞2、攻角α2、侧滑角β2依次为7、3000Pa、288.15K、10°、0°。
按发明内容中第一步描述方法,得来流条件II下流场压强初值p21。
步骤二:计算来流条件II下流场马赫数初值M21并修正p21:
按发明内容中第二步描述方法,得来流条件II下流场马赫数初值M21。根据M21修正后的流场压强初值p21如图7所示。
步骤三:计算来流条件II下流场温度初值T21与速度初值V21:
按发明内容中第三步描述方法,得来流条件II下流场温度初值T21与速度初值V21分布如图8、图9所示。
步骤四:计算来流条件II下流场速度矢量初值
按发明内容中第四步描述方法,得来流条件II下流场速度矢量初值
步骤五:插值得到待初始化状态对应网格G2上流场压强、温度、速度矢量的初值p2、T2、
待初始化状态对应网格G2的对称面网格如图10所示,前四步得到的p21、T21、对应网格节点为G1,由于G2与G1不一致,将G1上的p21、T21、插值至G2,得到待求的压强、温度、速度矢量初值p2、T2、
为分析本发明提出方法的特点,取本发明得到的初始流场与常用的基于来流参数的常值初始流场,分别进行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倍。分析表明,本发明提出的飞行器定常绕流数值求解的流场初始化在计算稳定性、收敛速度等方面皆明显优于常规常值初始化方法。
本发明提出了基于已知飞行状态的绕流定常解的高效、通用的飞行器定常绕流数值求解的流场初始化方法,提高了流场初始化质量,改善了流场数值求解的计算稳定性,加快了收敛速度,大幅提高了多状态飞行器定常绕流数值求解效率,是多状态飞行器定常绕流数值求解的理想初始化方法。
综上所述,虽然本发明已以较佳实施例揭露如上,然其并非用以限定本发明,任何本领域普通技术人员,在不脱离本发明的精神和范围内,当可作各种更动与润饰,因此本发明的保护范围当视权利要求书界定的范围为准。
Claims (1)
1.一种用于飞行器定常绕流数值求解的流场初始化方法,定义已知状态流场定常解对应网格节点集合为G1={(x1,y1,z1)},来流条件为“来流条件I”,对应马赫数、压强、温度、攻角、侧滑角依次为M∞1、p∞1、T∞1、α1、β1,流场定常解为“流场I”,网格节点上的压强、温度、速度矢量依次为p1、T1、待初始化的计算状态对应网格节点集合为G2={(x2,y2,z2)},来流条件为“来流条件II”,对应马赫数、压强、温度、攻角、侧滑角依次为M∞2、p∞2、T∞2、α2、β2,其特征在于:该方法根据G1、M∞1、p∞1、T∞1、α1、β1、p1、T1、给出来流条件II下的飞行器定常绕流数值求解问题的初始条件,即计算网格节点G2上的压强p2、温度T2、速度矢量初值具体步骤如下:
第一步:计算来流条件II下流场压强初值p21:
流场中压强系数Cp的计算公式为:
式中,p为压强;p∞、ρ∞、V∞依次为来流压强、密度、速度;
超声速条件下,驻点前存在正激波,正激波后压强满足如下关系:
式中,γ为来流比热比,ps为正激波后压强,M∞为来流马赫数;
综合式(1)、(2)推导得正激波后压强系数Cps:
按式(1)计算得来流条件I下流场的压强系数Cp1:
式中,p1为来流条件I下流场的压强;V∞1为来流条件I对应的来流速度,按下式计算:
式中,γ为来流比热比,R为来流气体常数;
分两种情况计算来流条件II下流场压强系数初值Cp2:
(a)当M∞1>1且M∞2>1时,按下式:
(b)当M∞1≤1或M∞2≤1时,按下式:
Cp2=Cp1; (7)
得来流条件II下压强初值p21计算公式:
式中,V∞2为来流条件II对应的来流速度,按下式计算:
第二步:计算来流条件II下流场马赫数初值M21并修正p21:
设流场总压比为流场总压与来流总压的比值;正激波后总压比rp0为:
流场中总压p0、压强p、马赫数M具有以下关系:
来流条件I下流场总压比rp01按下式计算:
式中,M1为来流条件I下流场的马赫数,按下式计算:
式中,V1为来流条件I下流场的速度,T1为来流条件I下流场的温度;
按下式计算来流条件II下流场总压比初值rp02:
式中,κ≤1,可取为0.8~0.99;
根据式(14)、(11)推导得来流条件II下流场马赫数初值M21计算公式:
式中,g为计算M21的中间变量,为来流条件II对应的总压;
当式(15)中g<0,第一步得到的p21不再满足式(11),按式(11)推导重新计算p21:
第三步:计算来流条件II下流场温度初值T21与速度初值V21:
按下式计算来流条件II下流场温度初值T21:
式中,为来流条件II对应的总温;
按下式计算来流条件II下流场温度初值V21:
式中,R为来流气体常数;
第四步:计算来流条件II下流场速度矢量初值
按下式计算来流条件II下流场速度矢量初值
式中:
第五步:插值得到待初始化状态对应网格G2上流场压强、温度、速度矢量初值p2、T2、
前四步得到的压强初值p21、温度初值T21、速度矢量初值对应网格节点为G1,若G1与G2一致,则p21、T21、的集合即为压强、温度、速度矢量初值p2、T2、若网格G1与网格G2不一致,则将网格G1上的压强初值p21、温度初值T21、速度矢量初值插值至网格G2,插值方法取计算量较小的最近邻点插值法,得到压强、温度、速度矢量初值p2、T2、
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610687769.4A CN106339534B (zh) | 2016-08-18 | 2016-08-18 | 一种用于飞行器定常绕流数值求解的流场初始化方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610687769.4A CN106339534B (zh) | 2016-08-18 | 2016-08-18 | 一种用于飞行器定常绕流数值求解的流场初始化方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN106339534A true CN106339534A (zh) | 2017-01-18 |
CN106339534B CN106339534B (zh) | 2019-04-12 |
Family
ID=57824846
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201610687769.4A Active CN106339534B (zh) | 2016-08-18 | 2016-08-18 | 一种用于飞行器定常绕流数值求解的流场初始化方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN106339534B (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106991209A (zh) * | 2017-03-01 | 2017-07-28 | 中国航天空气动力技术研究院 | 一种火星大气真实气体环境气动特性预测方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US5070458A (en) * | 1989-03-31 | 1991-12-03 | Honeywell Inc. | Method of analyzing and predicting both airplane and engine performance characteristics |
CN104573296A (zh) * | 2013-10-22 | 2015-04-29 | 北京临近空间飞行器系统工程研究所 | 一种面向相似网格的高超声速流场初始化方法 |
-
2016
- 2016-08-18 CN CN201610687769.4A patent/CN106339534B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US5070458A (en) * | 1989-03-31 | 1991-12-03 | Honeywell Inc. | Method of analyzing and predicting both airplane and engine performance characteristics |
CN104573296A (zh) * | 2013-10-22 | 2015-04-29 | 北京临近空间飞行器系统工程研究所 | 一种面向相似网格的高超声速流场初始化方法 |
Non-Patent Citations (2)
Title |
---|
向敏 等: "超空泡航行器内外流场仿真及性能分析", 《弹箭与制导学报》 * |
张代贤: "大气吸光船推力器工作特性数值计算研究", 《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》 * |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106991209A (zh) * | 2017-03-01 | 2017-07-28 | 中国航天空气动力技术研究院 | 一种火星大气真实气体环境气动特性预测方法 |
CN106991209B (zh) * | 2017-03-01 | 2020-07-14 | 中国航天空气动力技术研究院 | 一种火星大气真实气体环境气动特性预测方法 |
Also Published As
Publication number | Publication date |
---|---|
CN106339534B (zh) | 2019-04-12 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111339681B (zh) | 一种采用空气介质模拟发动机燃气介质喷流气动干扰效应的喷管出口参数匹配方法 | |
CN102650565B (zh) | 风洞模拟实验中涡轮动力模拟器短舱唇口及其设计方法 | |
CN106354903B (zh) | 用于飞行器定常绕流数值求解的计算域外边界确定方法 | |
CN105059530B (zh) | 一种后掠角可控的尖顶点密切锥乘波体 | |
CN109359379B (zh) | 一种自适应的超临界二氧化碳透平设计方法 | |
CN106679925A (zh) | 一种微质量射流流量高精度控制装置和控制方法 | |
CN103366052A (zh) | 一种高超声速飞行器热静气弹分析方法 | |
CN107832530B (zh) | 一种复杂外形的高超声速边界层转捩判定方法 | |
CN112947534A (zh) | 一种高超声速飞行器下压段自适应伪谱法轨迹优化方法 | |
CN104834785A (zh) | 基于单纯形样条函数的航空发动机稳态模型的建模方法 | |
CN104462653A (zh) | 一种发动机仿真设计方法 | |
CN111125972B (zh) | 核电厂破口失水事故水力载荷分析方法 | |
CN106339534A (zh) | 一种用于飞行器定常绕流数值求解的流场初始化方法 | |
Millikan | The boundary layer and skin friction for a figure of revolution | |
CN115079592A (zh) | 一种船舶核动力装置热力系统管网仿真方法 | |
Qu et al. | A parameter-free upwind scheme for all speeds’ simulations | |
Monier et al. | Les investigation of boussinesq constitutive relation validity in a corner separation flow | |
CN115879216A (zh) | 一种内流道强波系干扰控制下的流场重构设计方法 | |
CN102840960B (zh) | 一种用拉法尔管均匀风洞风场的方法 | |
Schwartz et al. | Scaling of distortion caused by crossing shock-wave/boundary-layer interactions | |
INTAKES | Effect of cowl angle in the performance of scramjet air intakes | |
CN204718765U (zh) | 一种可提供多路高压气流的调压系统 | |
Tufts et al. | Stability Equation Based Transition Prediction | |
Paxton et al. | Multi Fidelity Aerodynamic and Aerothermal Boundary Layer and Near Body Characterization of Hypersonic Configurations. | |
CN106202807B (zh) | 判别航天器身部激波/前缘类激波干扰发生条件及类型的方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |