CN113505551B - 来流变化诱发非定常的模拟方法、系统、存储介质和终端 - Google Patents

来流变化诱发非定常的模拟方法、系统、存储介质和终端 Download PDF

Info

Publication number
CN113505551B
CN113505551B CN202111053024.XA CN202111053024A CN113505551B CN 113505551 B CN113505551 B CN 113505551B CN 202111053024 A CN202111053024 A CN 202111053024A CN 113505551 B CN113505551 B CN 113505551B
Authority
CN
China
Prior art keywords
time
incoming flow
calculating
value
represented
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.)
Active
Application number
CN202111053024.XA
Other languages
English (en)
Other versions
CN113505551A (zh
Inventor
陈琦
谢昱飞
陈坚强
毛枚良
王新光
万钊
华如豪
郭勇颜
钟译漫
孙伟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Original Assignee
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Computational Aerodynamics Institute of China Aerodynamics Research and Development Center filed Critical Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority to CN202111053024.XA priority Critical patent/CN113505551B/zh
Publication of CN113505551A publication Critical patent/CN113505551A/zh
Application granted granted Critical
Publication of CN113505551B publication Critical patent/CN113505551B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computing Systems (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了来流变化诱发非定常的模拟方法、系统、存储介质和终端，方法包括：基于飞行器的外形，生成初始时刻的计算网格；计算初始时刻飞行器的定常流场；计算来流参数变化时的非定常流场；对非定常三维Navier‑stockes方程组进行离散化和线性化处理；计算n+1时刻左端项LHS的值；并在n+1时刻，更新来流参数参考量；根据当前时刻的参考值计算
Figure 592297DEST_PATH_IMAGE002
，根据上述n+1时刻的参考值计算
Figure 978279DEST_PATH_IMAGE004
，同时计算前一时刻的参考值
Figure 100004_DEST_PATH_IMAGE006
；计算n+1时刻的无量纲时间步长
Figure 100004_DEST_PATH_IMAGE008
；计算右端项RHS的值。本发明考虑来流参数变化诱导的气动非定常效应，采用了具有一般普适性的计算右端项RHS的公式，基于理论推导，对程序改动量小，计算量小。

Description

来流变化诱发非定常的模拟方法、系统、存储介质和终端
本发明涉及流体力学领域，尤其涉及来流变化诱发非定常的模拟方法、系统、存储介质和终端。
背景技术
来流参数变化在飞行器运动过程的模拟中广泛存在，在航空航天领域具有广泛的应用需求。目前，对来流参数变化的数值模拟过程，一般只考虑参数自身的变化，没有考虑参数变化产生后附加的非定常效应，如图1所示，导致数值模拟得到的气动力与真实飞行过程存在一定的误差。
因此在本领域中，如何能够提供一种从理论上推导了非定常效应的引入过程，以双时间步方法为基础，详细阐述了如何了对原双时间步方法进行改造，以考虑高度变化导致的气动非定常效应，并且从方法、系统、存储介质和终端的方向进行实现，属于本领域亟待解决的问题。
发明内容
本发明的目的在于克服现有技术的不足，提供来流变化诱发非定常的模拟方法、系统、存储介质和终端。
本发明的目的是通过以下技术方案来实现的：
本发明的第一方面，提供来流变化诱发非定常的模拟方法，包括以下步骤：
基于飞行器的外形，生成初始时刻的计算网格；
计算初始时刻飞行器的定常流场；
基于初始时刻的定常流场，计算来流参数变化时的非定常流场；
对非定常三维Navier-stockes方程组进行离散化和线性化处理后，可以简写成如下形式：
Figure 496166DEST_PATH_IMAGE002
Figure 678886DEST_PATH_IMAGE004
Figure DEST_PATH_IMAGE005
上式中，LHS表示左端时间离散项，RHS表示右端空间离散项，Δτ为伪时间推进步长，Δt为真实时间推进步长，I是单元矩阵，J -1是坐标变换矩阵，n表示时刻，s表示亚迭代指标，Q为守恒变量，Q=[ρ，ρu，ρv，ρw，e]-1，ρ是来流密度，u，v，w分别是三个方向的流动速度，e则是单位质量的总内能；
Figure 927465DEST_PATH_IMAGE006
，
Figure DEST_PATH_IMAGE007
，其中E、F和G为无粘通矢量，E v 、F v 和G v 为粘性通矢量，
Figure 214352DEST_PATH_IMAGE008
为坐标变换系数，
Figure DEST_PATH_IMAGE009
表示基于以无穷远来流参数计算的雷诺数，
Figure 763145DEST_PATH_IMAGE010
τ是剪切应力，q是热流，p为压力；i为张量运算符号，值取1、2、3，分别表示x、y、z三个方向的分量，x、y、z分别表示相应量在x、y、z三个方向的分量；
基于公式（2），计算n+1时刻左端项LHS的值；并在n+1时刻，来流密度、速度、压力或者温度发生变化后，更新来流参数参考量，包括：来流密度参考量
Figure DEST_PATH_IMAGE011
、速度参考量
Figure 382345DEST_PATH_IMAGE012
、压力参考量
Figure DEST_PATH_IMAGE013
、温度参考量
Figure 242854DEST_PATH_IMAGE014
；
根据当前时刻的参考值计算
Figure DEST_PATH_IMAGE015
，根据上述n+1时刻的参考值计算
Figure 176175DEST_PATH_IMAGE016
，同时计算前一时刻的参考值
Figure DEST_PATH_IMAGE017
；
根据
Figure 937064DEST_PATH_IMAGE018
，计算n+1时刻的无量纲时间步长
Figure DEST_PATH_IMAGE019
，L为特征长度；
将
Figure 196007DEST_PATH_IMAGE020
和
Figure 543812DEST_PATH_IMAGE019
代入公式（4），从而替代式（3）计算右端项RHS的值：
Figure 100002_DEST_PATH_IMAGE021
循环上述步骤直到完成整个参数变化过程的计算；
上述公式和参数中，加波浪线的表示有量纲量，无波浪线的为无量纲值，下标∞表示来流参数参考值，上标n表示第n个时间步的值。
进一步地，所述的来流参数变化包括飞行器距离地面高度h的变化，公式为：
Figure 343140DEST_PATH_IMAGE022
式中，t表示时间，T表示模拟的时间尺度，h 0表示初始时刻高度，h 1表示下一时刻高度；对于来流参数中的来流密度参考量
Figure 132105DEST_PATH_IMAGE011
、速度参考量
Figure 30791DEST_PATH_IMAGE012
、压力参考量
Figure 101777DEST_PATH_IMAGE013
、温度参考量
Figure 642480DEST_PATH_IMAGE014
，根据当前的高度h重新计算。
进一步地，所述T=0.09秒，h 0=50。
进一步地，在无量纲时间步长
Figure DEST_PATH_IMAGE023
的计算过程中，物理时间步长
Figure 754792DEST_PATH_IMAGE024
保持不变。
本发明的第二方面，提供来流变化诱发非定常的模拟系统，包括以下模块：
网格初始化模块：用于基于飞行器的外形，生成初始时刻的计算网格；
定常流场初始化模块：用于计算初始时刻飞行器的定常流场；
非定常流场计算模块：基于初始时刻的定常流场，计算来流参数变化时的非定常流场；
方程组转化模块：对非定常三维Navier-stockes方程组进行离散化和线性化处理后，可以简写成如下形式：
Figure 417855DEST_PATH_IMAGE002
Figure 412356DEST_PATH_IMAGE004
Figure 694432DEST_PATH_IMAGE005
上式中，LHS表示左端时间离散项，RHS表示右端空间离散项，Δτ为伪时间推进步长，Δt为真实时间推进步长，I是单元矩阵，J -1是坐标变换矩阵，n表示时刻，s表示亚迭代指标，Q为守恒变量，Q=[ρ，ρu，ρv，ρw，e]-1，ρ是来流密度，u，v，w分别是三个方向的流动速度，e则是单位质量的总内能；
Figure 457989DEST_PATH_IMAGE006
，
Figure 557532DEST_PATH_IMAGE007
，其中E、F和G为无粘通矢量，E v 、F v 和G v 为粘性通矢量，
Figure 977012DEST_PATH_IMAGE008
为坐标变换系数，
Figure 125097DEST_PATH_IMAGE009
表示基于以无穷远来流参数计算的雷诺数，
Figure 71056DEST_PATH_IMAGE010
τ是剪切应力，q是热流，p为压力，i为张量运算符号，值取1、2、3，分别表示x、y、z三个方向的分量，x、y、z分别表示相应量在x、y、z三个方向的分量；
左端项LHS计算模块：用于基于公式（2），计算n+1时刻左端项LHS的值；
来流参数参考量更新模块：用于在n+1时刻，来流密度、速度、压力或者温度发生变化后，更新来流参数参考量，包括：来流密度参考量
Figure 216867DEST_PATH_IMAGE011
、速度参考量
Figure 920380DEST_PATH_IMAGE012
、压力参考量
Figure 190866DEST_PATH_IMAGE013
、温度参考量
Figure 929015DEST_PATH_IMAGE014
；
前后时刻守恒变量计算模块：用于根据当前时刻的参考值计算
Figure 245727DEST_PATH_IMAGE015
，根据上述n+1时刻的参考值计算
Figure 436537DEST_PATH_IMAGE016
，同时计算前一时刻的参考值
Figure 988741DEST_PATH_IMAGE017
；
无量纲时间步长计算模块：用于根据
Figure 315817DEST_PATH_IMAGE018
，计算n+1时刻的无量纲时间步长
Figure 803430DEST_PATH_IMAGE019
，L为特征长度；
右端项RHS计算模块：用于将
Figure 481536DEST_PATH_IMAGE020
和
Figure 837431DEST_PATH_IMAGE019
代入公式（4），从而替代式（3）计算右端项RHS的值：
Figure 487855DEST_PATH_IMAGE021
循环计算步骤：用于循环利用上述模块直到完成整个参数变化过程的计算；
上述公式和参数中，加波浪线的表示有量纲量，无波浪线的为无量纲值，下标∞表示来流参数参考值，上标n表示第n个时间步的值。
进一步地，所述的来流参数变化包括飞行器距离地面高度h的变化，公式为：
Figure 943107DEST_PATH_IMAGE022
式中，t表示时间，T表示模拟的时间尺度，h 0表示初始时刻高度，h 1表示下一时刻高度；对于来流参数中的来流密度参考量
Figure 170826DEST_PATH_IMAGE011
、速度参考量
Figure 2516DEST_PATH_IMAGE012
、压力参考量
Figure 507447DEST_PATH_IMAGE013
、温度参考量
Figure 133600DEST_PATH_IMAGE014
，根据当前的高度h重新计算。
进一步地，所述T=0.09秒，h 0=50。
进一步地，在无量纲时间步长计算模块中的无量纲时间步长
Figure 350080DEST_PATH_IMAGE023
的计算过程中，物理时间步长
Figure 985461DEST_PATH_IMAGE024
保持不变。
本发明的第三方面，提供一种存储介质，其上存储有计算机指令，所述计算机指令运行时执行所述的来流变化诱发非定常的模拟方法的步骤。
本发明的第四方面，提供一种终端，包括存储器和处理器，所述存储器上存储有可在所述处理器上运行的计算机指令，所述处理器运行所述计算机指令时执行所述的来流变化诱发非定常的模拟方法的步骤。
本发明的有益效果是：
在本发明的一示例性实施例中，为考虑来流参数变化诱导的气动非定常效应，采用了具有一般普适性的计算右端项RHS的公式。该方法基于理论推导，对程序改动量小，计算量小，在模拟来流参数剧烈变化时必须使用本方法。
本发明的系统、存储介质和终端也具有相同优点。
附图说明
图1为现有技术计算方法流程示意图；
图2为本发明一示例性实施例中提供的方法流程示意图；
图3为本发明一示例性实施例中生成的计算网格示意图；
图4为本发明一示例性实施例中飞行器飞行高度非定常变化过程中计算得到的轴向力与定常时的气动力比较的示意图。
具体实施方式
下面结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
在本发明的描述中，需要说明的是，属于“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方向或位置关系为基于附图所述的方向或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。
在本发明的描述中，需要说明的是，除非另有明确的规定和限定，属于“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
在本申请使用的术语是仅仅出于描述特定实施例的目的，而非旨在限制本申请。在本申请和所附权利要求书中所使用的单数形式的“一种”、“所述”和“该”也旨在包括多数形式，除非上下文清楚地表示其他含义。还应当理解，本文中使用的术语“和/或”是指并包含一个或多个相关联的列出项目的任何或所有可能组合。
应当理解，尽管在本申请可能采用术语第一、第二、第三等来描述各种信息，但这些信息不应限于这些术语。这些术语仅用来将同一类型的信息彼此区分开。例如，在不脱离本申请范围的情况下，第一信息也可以被称为第二信息，类似地，第二信息也可以被称为第一信息。取决于语境，如在此所使用的词语“如果”可以被解释成为“在……时”或“当……时”或“响应于确定”。此外，属于“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性。
此外，下面所描述的本发明不同实施方式中所涉及的技术特征只要彼此之间未构成冲突就可以相互结合。
参见图2，图2示出了本发明的一示例性实施例提供来流变化诱发非定常的模拟方法（来流参数变化诱发非定常效应的模拟方法），包括以下步骤：
S1：基于飞行器的外形，生成初始时刻的计算网格，在其中一示例性实施例中生成的计算网格如图3所示；
S2：计算初始时刻飞行器的定常流场；
S3：基于初始时刻的定常流场，计算来流参数变化时的非定常流场；
S4：对非定常三维Navier-stockes方程组进行离散化和线性化处理后，可以简写成如下形式：
Figure 79319DEST_PATH_IMAGE002
Figure 141953DEST_PATH_IMAGE004
Figure 78685DEST_PATH_IMAGE005
上式中，LHS表示左端时间离散项，RHS表示右端空间离散项，Δτ为伪时间推进步长，Δt为真实时间推进步长，I是单元矩阵，J -1是坐标变换矩阵，n表示时刻，s表示亚迭代指标，Q为守恒变量，Q=[ρ，ρu，ρv，ρw，e]-1，ρ是来流密度，u，v，w分别是三个方向的流动速度，e则是单位质量的总内能；
Figure 721019DEST_PATH_IMAGE006
，
Figure 997279DEST_PATH_IMAGE007
，其中E、F和G为无粘通矢量，E v 、F v 和G v 为粘性通矢量，
Figure 293131DEST_PATH_IMAGE008
为坐标变换系数，
Figure 654843DEST_PATH_IMAGE009
表示基于以无穷远来流参数计算的雷诺数，
Figure 569709DEST_PATH_IMAGE010
τ是剪切应力，q是热流，p为压力，为张量运算符号，值取1、2、3，分别表示x、y、z三个方向的分量，x、y、z分别表示相应量在x、y、z三个方向的分量；
S5：基于公式（2），计算n+1时刻左端项LHS的值；并在n+1时刻，来流密度、速度、压力或者温度发生变化后，更新来流参数参考量，包括：来流密度参考量
Figure 700476DEST_PATH_IMAGE011
、速度参考量
Figure 167230DEST_PATH_IMAGE012
、压力参考量
Figure 16237DEST_PATH_IMAGE013
、温度参考量
Figure 469215DEST_PATH_IMAGE014
；
其中，需要说明的是：密度、速度（三个方向分量）、压力，这五个变量被称为原始变量，求解方程（1）的最终目的就是要得到这五个变量的空间分布和时间分布情况。所以速度、密度和压力这些物理量在求解过程中一直都需要。
温度可根据空气状态方程
Figure DEST_PATH_IMAGE025
计算得到（其中，ρ是来流密度，R为气体常数，P为压力，T为温度）；相应地，温度变化后对压力和密度也会产生影响；另外温度变化后，声速
Figure 812078DEST_PATH_IMAGE026
也会变化（其中，
Figure DEST_PATH_IMAGE027
表示比热比），相应的速度V也会变化。也就是说温度和压力、密度、速度都相关，如果来流的温度变了，其他的值相应的也要改变。
S6：根据当前时刻的参考值计算
Figure 590678DEST_PATH_IMAGE015
，根据上述n+1时刻的参考值计算
Figure 254878DEST_PATH_IMAGE016
，同时计算前一时刻的参考值
Figure 511547DEST_PATH_IMAGE017
；
S7：根据
Figure 616906DEST_PATH_IMAGE018
，计算n+1时刻的无量纲时间步长
Figure 425462DEST_PATH_IMAGE019
，L为特征长度（一般为飞行器的长度或者直径）；
将
Figure 452324DEST_PATH_IMAGE020
和
Figure 637318DEST_PATH_IMAGE019
代入公式（4），从而替代式（3）计算右端项RHS的值：
Figure 597183DEST_PATH_IMAGE021
循环上述步骤（S5~S7）直到完成整个参数变化过程的计算；
其中，需要说明的是：先计算右端项RHS的值，计算完成后会得到原始变量密度、压力、速度在空间的分布；再计算左端项LHS的值，得到原始变量在不同时刻的值（每一时刻原始变量在空间的分布都不相同）。从t=0开始，重复上述过程，一直到计算结束，得到原始变量的时空分布情况。
更为具体地：先计算右端空间离散项RHS，进行伪时间推进，获得离散空间上的收敛，再计算左端时间离散项LHS进行真实时间推进。
另外，对于RHS的计算值（即公式4）对应原始变量密度、压力、速度在空间的分布：右端空间离散项RHS收敛，获得守恒变量Q的值，Q=[ρ，ρu，ρv，ρw，e]-1，ρ是来流密度，u，v，w分别是三个方向的流动速度，e则是单位质量的总内能，
Figure 717586DEST_PATH_IMAGE028
，从而获得原始变量密度、压力、速度在空间的分布。
另外，计算完成后，得到的是包含了密度、压力、速度的流场数据，这些数据包含了物体表面的数据和物体周围的流动数据。对物体表面的压力等数据积分，就可以得到气动力，气动力包含六个分量：轴向力、法向力、侧向力，俯仰力矩/偏航力矩/滚转力矩（即对应于图4中的轴向力和气动力）。
上述公式和参数中，加波浪线的表示有量纲量，无波浪线的为无量纲值，下标∞表示来流参数参考值，上标n表示第n个时间步的值。
需要说明的是，在上述示例性实施例中，在来流参数无变化时，
Figure DEST_PATH_IMAGE029
，
Figure 857843DEST_PATH_IMAGE030
，则公式（4）与公式（3）等价；在来流参数变化时，两式不相同，为考虑来流参数变化诱导的气动非定常效应，必须采用公式（4）替换公式（3）。公式（4）比公式（3）更具有一般普适性。该方法基于理论推导，对程序改动量小，计算量小，在模拟来流参数剧烈变化时必须使用本方法。
具体地，对于步骤S4~S7的理论推导如下：
描述非定常流动的守恒形式Navier-Stokes方程组为：
Figure DEST_PATH_IMAGE031
其中Q为守恒变量，包含了密度、三方向速度、流体的总能等参数信息。为求解上述非线性方程组，首先需要对各物理量进行无量纲化，无量纲化的参考值一般根据来流的参数进行计算，在常规的CFD（Computational Fluid Dynamics）计算中，来流参数一般保持不变，因此，以方程第一项为例，其一般采用如下方式进行离散：
Figure 721894DEST_PATH_IMAGE032
上式中的各变量均为无量纲的值，将上述无量纲量还原为有量纲的形式，可以得到
Figure DEST_PATH_IMAGE033
上式中加波浪线的表示有量纲量，无波浪线的为无量纲值，下标∞表示来流参考值，上标n表示第n个时间步的值，在来流参数无变化时，不同时刻的参考值相同，即
Figure 864162DEST_PATH_IMAGE034
将（9）式代入（8）式，可以直接回到（7）式，即来流参数不变时，（7）式与（6）式等价；但如果来流参数变化时，（9）式不再成立，（8）式与（7）式不再等价，必须重新求解（8）式，对（8）式整理后得到
Figure DEST_PATH_IMAGE035
与（7）式相比，（10）式项较为复杂，直接求解很不方便，程序改动也较大，不利于推广。考虑到一般在求解非定常问题时，物理时间步长一般不变，也即是
Figure 421045DEST_PATH_IMAGE036
由此，（10）式可以简化为
Figure DEST_PATH_IMAGE037
公式（12）从理论上指出，当来流参数变化时，需要对相关变量的时间离散项进行相应的改造，引入不同时刻的来流参考值。下面以双时间步方法为例，详细介绍具体改造的过程。
为方便描述，（6）式可以写成如下半离散形式
Figure 281554DEST_PATH_IMAGE038
非定常计算一般采用双时间步方法，不考虑高度变化的影响，可写成如下形式
Figure DEST_PATH_IMAGE039
对
Figure 949296DEST_PATH_IMAGE040
项进行线性化处理，
Figure DEST_PATH_IMAGE041
Figure 946070DEST_PATH_IMAGE042
，于是有
Figure 408276DEST_PATH_IMAGE044
根据（12）式，在来流参数变化时，需要对时间离散项进行相应的改造，（11）式不需要修改，（3）式改造为
Figure DEST_PATH_IMAGE045
也就是说，在求解方程（1）的右端项时，使用公式（4）替换了之前普遍采用的公式（3），即用
Figure 254616DEST_PATH_IMAGE046
替换了
Figure DEST_PATH_IMAGE047
。在来流参数变化时，两式不相同，为考虑来流参数变化诱导的气动非定常效应，必须采用公式（4）替换公式（3）。公式（4）比公式（3）更具有一般普适性。该方法基于理论推导，对程序改动量小，计算量小，在模拟来流参数剧烈变化时必须使用本方法。
更优地，在一示例性实施例中，所述的来流参数变化包括飞行器距离地面高度h的变化，公式为：
Figure 460469DEST_PATH_IMAGE048
式中，t表示时间，T表示模拟的时间尺度，h 0表示初始时刻高度，h 1表示下一时刻高度；对于来流参数中的来流密度参考量
Figure 249433DEST_PATH_IMAGE011
、速度参考量
Figure 7174DEST_PATH_IMAGE012
、压力参考量
Figure 248799DEST_PATH_IMAGE013
、温度参考量
Figure 992764DEST_PATH_IMAGE014
，根据当前的高度h重新计算。
更优地，在一示例性实施例中，所述T=0.09秒，h 0=50；同时h 1=70。
更优地，在一示例性实施例中，在无量纲时间步长
Figure 901815DEST_PATH_IMAGE023
的计算过程中，物理时间步长
Figure 564877DEST_PATH_IMAGE024
保持不变。
图4示出了飞行器飞行高度非定常变化过程中计算得到的轴向力与定常时的轴向力（轴向力是气动力的一种，气动力包含六个分量：轴向力、法向力、侧向力，俯仰力矩/偏航力矩/滚转力矩）比较，可以看到，随着飞行高度的上升和攻角的增大，其气动力并非是线性增加的（如图4的实线表示），而采用本发明方法计算的结果可以反应来流参数变化诱导的气动力非定常效应。
与上述示例性实施例具有相同的发明构思，本发明的的又一示例性实施例提供来流变化诱发非定常的模拟系统，包括以下模块：
网格初始化模块：用于基于飞行器的外形，生成初始时刻的计算网格；
定常流场初始化模块：用于计算初始时刻飞行器的定常流场；
非定常流场计算模块：基于初始时刻的定常流场，计算来流参数变化时的非定常流场；
方程组转化模块：对非定常三维Navier-stockes方程组进行离散化和线性化处理后，可以简写成如下形式：
Figure 559378DEST_PATH_IMAGE002
Figure DEST_PATH_IMAGE049
Figure 107034DEST_PATH_IMAGE005
上式中，LHS表示左端时间离散项，RHS表示右端空间离散项，Δτ为伪时间推进步长，Δt为真实时间推进步长，I是单元矩阵，J -1是坐标变换矩阵，n表示时刻，s表示亚迭代指标，Q为守恒变量，Q=[ρ，ρu，ρv，ρw，e]-1，ρ是来流密度，u，v，w分别是三个方向的流动速度，e则是单位质量的总内能；
Figure 932908DEST_PATH_IMAGE006
，
Figure 704554DEST_PATH_IMAGE007
，其中E、F和G为无粘通矢量，E v 、F v 和G v 为粘性通矢量，
Figure 124034DEST_PATH_IMAGE008
为坐标变换系数，
Figure 272119DEST_PATH_IMAGE009
表示基于以无穷远来流参数计算的雷诺数，
Figure 453964DEST_PATH_IMAGE010
τ是剪切应力，q是热流，p为压力，i为张量运算符号，值取1、2、3，分别表示x、y、z三个方向的分量，x、y、z分别表示相应量在x、y、z三个方向的分量；
左端项LHS计算模块：用于基于公式（2），计算n+1时刻左端项LHS的值；
来流参数参考量更新模块：用于在n+1时刻，来流密度、速度、压力或者温度发生变化后，更新来流参数参考量，包括：来流密度参考量
Figure 662091DEST_PATH_IMAGE011
、速度参考量
Figure 303288DEST_PATH_IMAGE012
、压力参考量
Figure 989485DEST_PATH_IMAGE013
、温度参考量
Figure 789950DEST_PATH_IMAGE014
；
前后时刻守恒变量计算模块：用于根据当前时刻的参考值计算
Figure 903400DEST_PATH_IMAGE015
，根据上述n+1时刻的参考值计算
Figure 31893DEST_PATH_IMAGE016
，同时计算前一时刻的参考值
Figure 521780DEST_PATH_IMAGE017
；
无量纲时间步长计算模块：用于根据
Figure 176752DEST_PATH_IMAGE018
，计算n+1时刻的无量纲时间步长
Figure 726682DEST_PATH_IMAGE019
，L为特征长度；
右端项RHS计算模块：用于将
Figure 342471DEST_PATH_IMAGE020
和
Figure 370470DEST_PATH_IMAGE019
代入公式（4），从而替代式（3）计算右端项RHS的值：
Figure 879949DEST_PATH_IMAGE021
循环计算步骤：用于循环利用上述模块直到完成整个参数变化过程的计算；
上述公式和参数中，加波浪线的表示有量纲量，无波浪线的为无量纲值，下标∞表示来流参数参考值，上标n表示第n个时间步的值。
对应地，在一示例性实施例中，所述的来流参数变化包括飞行器距离地面高度h的变化，公式为：
Figure 335201DEST_PATH_IMAGE022
式中，t表示时间，T表示模拟的时间尺度，h 0表示初始时刻高度，h 1表示下一时刻高度；对于来流参数中的来流密度参考量
Figure 500603DEST_PATH_IMAGE011
、速度参考量
Figure 269976DEST_PATH_IMAGE012
、压力参考量
Figure 571644DEST_PATH_IMAGE013
、温度参考量
Figure 289808DEST_PATH_IMAGE014
，根据当前的高度h重新计算。
对应地，在一示例性实施例中，所述T=0.09秒，h 0=50。
对应地，在一示例性实施例中，在无量纲时间步长计算模块中的无量纲时间步长
Figure 880190DEST_PATH_IMAGE023
的计算过程中，物理时间步长
Figure 249991DEST_PATH_IMAGE024
保持不变。
本发明的又一示例性实施例，提供一种存储介质，其上存储有计算机指令，所述计算机指令运行时执行上述任一示例性实施例中所述的来流变化诱发非定常的模拟方法的步骤。
本发明的又一示例性实施例，提供一种终端，包括存储器和处理器，所述存储器上存储有可在所述处理器上运行的计算机指令，所述处理器运行所述计算机指令时执行上述任一示例性实施例中所述的来流变化诱发非定常的模拟方法的步骤。
基于这样的理解，本实施例的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(Read-OnlyMemory，ROM)、随机存取存储器(RandomAccessMemory，RAM)、磁碟或者光盘等各种可以存储程序代码的介质。
显然，上述实施例仅仅是为清楚地说明所作的举例，而并非对实施方式的限定，对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其他不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引申出的显而易见的变化或变动仍处于本发明创造的保护范围之中。

Claims (10)

1.来流变化诱发非定常的模拟方法，其特征在于：包括以下步骤：
基于飞行器的外形，生成初始时刻的计算网格；
计算初始时刻飞行器的定常流场；
基于初始时刻的定常流场，计算来流参数变化时的非定常流场；
对非定常三维Navier-stockes方程组进行离散化和线性化处理后，简写成如下形式：
Figure DEST_PATH_IMAGE002
Figure 227017DEST_PATH_IMAGE003
Figure DEST_PATH_IMAGE004
上式中，LHS表示左端时间离散项，RHS表示右端空间离散项，Δτ为伪时间推进步长，Δt为真实时间推进步长，I是单元矩阵，J -1是坐标变换矩阵，n表示时刻，s表示亚迭代指标，Q为守恒变量，Q=[ρ，ρu，ρv，ρw，e]-1，ρ是来流密度，u，v，w分别是三个方向的流动速度，e则是单位质量的总内能；
Figure 32905DEST_PATH_IMAGE005
，
Figure DEST_PATH_IMAGE006
，其中E、F和G为无粘通矢量，E v 、F v 和G v 为粘性通矢量，
Figure 343801DEST_PATH_IMAGE007
为坐标变换系数，
Figure DEST_PATH_IMAGE008
表示基于以无穷远来流参数计算的雷诺数，
Figure 942272DEST_PATH_IMAGE009
τ是剪切应力，q是热流，p为压力；i为张量运算符号，值取1、2、3，分别表示x、y、z三个方向的分量，x、y、z分别表示相应量在x、y、z三个方向的分量；
基于公式（2），计算n+1时刻左端项LHS的值；并在n+1时刻，来流密度、速度、压力或者温度发生变化后，更新来流参数参考量，包括：来流密度参考量
Figure DEST_PATH_IMAGE010
、速度参考量
Figure 412437DEST_PATH_IMAGE011
、压力参考量
Figure DEST_PATH_IMAGE012
、温度参考量
Figure 500479DEST_PATH_IMAGE013
；
根据当前时刻的参考值计算
Figure DEST_PATH_IMAGE014
，根据上述n+1时刻的参考值计算
Figure 970774DEST_PATH_IMAGE015
，同时计算前一时刻的参考值
Figure DEST_PATH_IMAGE016
；
根据
Figure 326931DEST_PATH_IMAGE017
，计算n+1时刻的无量纲时间步长
Figure DEST_PATH_IMAGE018
，L为特征长度；
将
Figure 261389DEST_PATH_IMAGE019
和
Figure 926857DEST_PATH_IMAGE018
代入公式（4），从而替代式（3）计算右端项RHS的值：
Figure DEST_PATH_IMAGE020
循环上述步骤直到完成整个参数变化过程的计算；
上述公式和参数中，加波浪线的表示有量纲量，无波浪线的为无量纲值，下标∞表示来流参数参考值，上标n表示第n个时间步的值。
2.根据权利要求1所述的来流变化诱发非定常的模拟方法，其特征在于：所述的来流参数变化包括飞行器距离地面高度h的变化，公式为：
Figure DEST_PATH_IMAGE021
式中，t表示时间，T表示模拟的时间尺度，h 0表示初始时刻高度，h 1表示下一时刻高度；对于来流参数中的来流密度参考量
Figure 868137DEST_PATH_IMAGE010
、速度参考量
Figure 339570DEST_PATH_IMAGE011
、压力参考量
Figure 394113DEST_PATH_IMAGE012
、温度参考量
Figure 496062DEST_PATH_IMAGE013
，根据当前的高度h重新计算。
3.根据权利要求2所述的来流变化诱发非定常的模拟方法，其特征在于：所述T=0.09秒，h 0=50。
4.根据权利要求1所述的来流变化诱发非定常的模拟方法，其特征在于：在无量纲时间步长
Figure DEST_PATH_IMAGE022
的计算过程中，物理时间步长
Figure 173905DEST_PATH_IMAGE023
保持不变。
5.来流变化诱发非定常的模拟系统，其特征在于：包括以下模块：
网格初始化模块：用于基于飞行器的外形，生成初始时刻的计算网格；
定常流场初始化模块：用于计算初始时刻飞行器的定常流场；
非定常流场计算模块：基于初始时刻的定常流场，计算来流参数变化时的非定常流场；
方程组转化模块：对非定常三维Navier-stockes方程组进行离散化和线性化处理后，简写成如下形式：
Figure DEST_PATH_IMAGE024
Figure 400094DEST_PATH_IMAGE025
Figure DEST_PATH_IMAGE026
上式中，LHS表示左端时间离散项，RHS表示右端空间离散项，Δτ为伪时间推进步长，Δt为真实时间推进步长，I是单元矩阵，J -1是坐标变换矩阵，n表示时刻，s表示亚迭代指标，Q为守恒变量，Q=[ρ，ρu，ρv，ρw，e]-1，ρ是来流密度，u，v，w分别是三个方向的流动速度，e则是单位质量的总内能；
Figure 918931DEST_PATH_IMAGE005
，
Figure 221474DEST_PATH_IMAGE006
，其中E、F和G为无粘通矢量，E v 、F v 和G v 为粘性通矢量，
Figure 635882DEST_PATH_IMAGE007
为坐标变换系数，
Figure 42592DEST_PATH_IMAGE008
表示基于以无穷远来流参数计算的雷诺数，
Figure 228985DEST_PATH_IMAGE009
τ是剪切应力，q是热流，p为压力，i为张量运算符号，值取1、2、3，分别表示x、y、z三个方向的分量，x、y、z分别表示相应量在x、y、z三个方向的分量；
左端项LHS计算模块：用于基于公式（2），计算n+1时刻左端项LHS的值；
来流参数参考量更新模块：用于在n+1时刻，来流密度、速度、压力或者温度发生变化后，更新来流参数参考量，包括：来流密度参考量
Figure 938315DEST_PATH_IMAGE010
、速度参考量
Figure 810325DEST_PATH_IMAGE011
、压力参考量
Figure 223989DEST_PATH_IMAGE012
、温度参考量
Figure 982998DEST_PATH_IMAGE013
；
前后时刻守恒变量计算模块：用于根据当前时刻的参考值计算
Figure 456704DEST_PATH_IMAGE014
，根据上述n+1时刻的参考值计算
Figure 580125DEST_PATH_IMAGE015
，同时计算前一时刻的参考值
Figure 407266DEST_PATH_IMAGE016
；
无量纲时间步长计算模块：用于根据
Figure 410994DEST_PATH_IMAGE017
，计算n+1时刻的无量纲时间步长
Figure 180236DEST_PATH_IMAGE018
，L为特征长度；
右端项RHS计算模块：用于将
Figure 636625DEST_PATH_IMAGE019
和
Figure 267458DEST_PATH_IMAGE018
代入公式（4），从而替代式（3）计算右端项RHS的值：
Figure 860113DEST_PATH_IMAGE027
循环计算步骤：用于循环利用上述模块直到完成整个参数变化过程的计算；
上述公式和参数中，加波浪线的表示有量纲量，无波浪线的为无量纲值，下标∞表示来流参数参考值，上标n表示第n个时间步的值。
6.根据权利要求5所述的来流变化诱发非定常的模拟系统，其特征在于：所述的来流参数变化包括飞行器距离地面高度h的变化，公式为：
Figure 301721DEST_PATH_IMAGE021
式中，t表示时间，T表示模拟的时间尺度，h 0表示初始时刻高度，h 1表示下一时刻高度；对于来流参数中的来流密度参考量
Figure 510986DEST_PATH_IMAGE010
、速度参考量
Figure 945509DEST_PATH_IMAGE011
、压力参考量
Figure 861513DEST_PATH_IMAGE012
、温度参考量
Figure 582344DEST_PATH_IMAGE013
，根据当前的高度h重新计算。
7.根据权利要求6所述的来流变化诱发非定常的模拟系统，其特征在于：所述T=0.09秒，h 0=50。
8.根据权利要求5所述的来流变化诱发非定常的模拟系统，其特征在于：在无量纲时间步长计算模块中的无量纲时间步长
Figure 403538DEST_PATH_IMAGE022
的计算过程中，物理时间步长
Figure 500807DEST_PATH_IMAGE023
保持不变。
9.一种计算机存储介质，其上存储有计算机指令，其特征在于：所述计算机指令运行时执行权利要求1至4中任一项所述的来流变化诱发非定常的模拟方法的步骤。
10.一种终端，包括存储器和处理器，所述存储器上存储有可在所述处理器上运行的计算机指令，其特征在于，所述处理器运行所述计算机指令时执行权利要求1至4中任一项所述的来流变化诱发非定常的模拟方法的步骤。
CN202111053024.XA 2021-09-09 2021-09-09 来流变化诱发非定常的模拟方法、系统、存储介质和终端 Active CN113505551B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111053024.XA CN113505551B (zh) 2021-09-09 2021-09-09 来流变化诱发非定常的模拟方法、系统、存储介质和终端

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111053024.XA CN113505551B (zh) 2021-09-09 2021-09-09 来流变化诱发非定常的模拟方法、系统、存储介质和终端

Publications (2)

Publication Number Publication Date
CN113505551A CN113505551A (zh) 2021-10-15
CN113505551B true CN113505551B (zh) 2021-12-07

Family

ID=78016676

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111053024.XA Active CN113505551B (zh) 2021-09-09 2021-09-09 来流变化诱发非定常的模拟方法、系统、存储介质和终端

Country Status (1)

Country Link
CN (1) CN113505551B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115840992B (zh) * 2023-02-20 2023-05-26 中国空气动力研究与发展中心计算空气动力研究所 一种弹性飞行器飞行仿真方法、系统、计算机存储介质及终端

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110162822A (zh) * 2019-03-19 2019-08-23 北京机电工程研究所 耦合结构模态的时域快速非定常气动力计算方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103136450B (zh) * 2013-02-07 2016-02-17 南京航空航天大学 超声速下飞行器表面侵蚀量衡量方法
CN105808954A (zh) * 2016-03-11 2016-07-27 中国航天空气动力技术研究院 一种适用于cfd数值模拟的周期非定常流场的预测方法
CN107273625A (zh) * 2017-06-22 2017-10-20 西南科技大学 一种三维不可压缩非定常n‑s方程有限元数值求解方法
CN109184952B (zh) * 2018-08-21 2019-06-18 西安理工大学 一种高超进气道不启动状态分离区自持能力定量分析方法
US11526182B2 (en) * 2019-03-25 2022-12-13 Cbn Nano Technologies Inc. Sensing and operation of devices in viscous flow using derived parameters to reduce data-handling requirements
CN111859530B (zh) * 2020-06-11 2022-04-22 北京航空航天大学 一种飞行器动态气动特性模拟的迭代推进扰动域更新方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110162822A (zh) * 2019-03-19 2019-08-23 北京机电工程研究所 耦合结构模态的时域快速非定常气动力计算方法

Also Published As

Publication number Publication date
CN113505551A (zh) 2021-10-15

Similar Documents

Publication Publication Date Title
CN113111430B (zh) 基于非线性气动力降阶的弹性飞机飞行动力学建模方法
CN108052772A (zh) 一种基于结构降阶模型的几何非线性静气动弹性分析方法
CN104133933A (zh) 一种高超声速飞行器热环境下气动弹性力学特性分析方法
CN114168796B (zh) 一种建立飞行器高空气动力数据库的方法
Zhang et al. High-fidelity aerostructural optimization with integrated geometry parameterization and mesh movement
CN108363843A (zh) 一种基于结构降阶模型的几何非线性静气动弹性全机配平方法
CN109002630B (zh) 一种超弹性材料的快速仿真方法
CN111046615B (zh) 一种黎曼解算器激波不稳定性抑制方法及系统
CN113505551B (zh) 来流变化诱发非定常的模拟方法、系统、存储介质和终端
Castellani et al. Flight loads prediction of high aspect ratio wing aircraft using multibody dynamics
CN115525980A (zh) 一种再入飞行器气动外形的优化方法和优化装置
Im et al. Prediction of a supersonic wing flutter boundary using a high fidelity detached eddy simulation
Oktay et al. Numerical investigation of effects of airspeed and rotational speed on quadrotor UAV propeller thrust coefficient
He et al. HyperFLOW: A structured/unstructured hybrid integrated computational environment for multi-purpose fluid simulation
CN114218672A (zh) 跨超声速大迎角配平翼的抖振响应获取方法及相关装置
Huang et al. Aeroelastic simulation using CFD/CSD coupling based on precise integration method
Wang et al. A novel method for estimating three-domain limit cycles in a 3D wing-aileron model with freeplay in aileron deflection
Verri et al. Static loads evaluation in a flexible aircraft using high-fidelity fluid–structure iteration tool (E2-FSI): extended version
CN103577649B (zh) 运输类飞机货物空投时货舱地板载荷的确定方法
Wang et al. Multi-parametric investigations on aerodynamic force, aeroacoustic, and engine energy utilizations based development of intercity bus associates with various drag reduction techniques through advanced engineering approaches
CN109635370A (zh) 开裂式阻力方向舵静气动弹性特性分析方法
Reasor et al. X-56A Aeroelastic Flight Test Predictions
Xu et al. Aerodynamic optimization based on continuous adjoint method for a flexible wing
Chimakurthi A computational aeroelasticity framework for analyzing flapping wings
Marciniuk et al. Aerodynamic Analysis of Variable Camber-Morphing Airfoils with Substantial Camber Deflections

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant