CN114218878B - 一种飞行器机动过程模拟的动网格扰动域更新方法 - Google Patents

一种飞行器机动过程模拟的动网格扰动域更新方法 Download PDF

Info

Publication number
CN114218878B
CN114218878B CN202210144150.4A CN202210144150A CN114218878B CN 114218878 B CN114218878 B CN 114218878B CN 202210144150 A CN202210144150 A CN 202210144150A CN 114218878 B CN114218878 B CN 114218878B
Authority
CN
China
Prior art keywords
grid
dynamic
domain
unit
boundary
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
CN202210144150.4A
Other languages
English (en)
Other versions
CN114218878A (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.)
Beihang University
Original Assignee
Beihang University
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 Beihang University filed Critical Beihang University
Priority to CN202210144150.4A priority Critical patent/CN114218878B/zh
Publication of CN114218878A publication Critical patent/CN114218878A/zh
Application granted granted Critical
Publication of CN114218878B publication Critical patent/CN114218878B/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
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • 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
    • 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)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Evolutionary Computation (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computational Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Physics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明涉及计算流体力学技术领域，提供了一种飞行器机动过程模拟的动网格扰动域更新方法。该方法，首先在流动控制方程求解中考虑网格单元速度和动网格几何守恒律的影响，在动网格技术的框架下建立仅对非定常未收敛受扰单元求解、仅在局部区域考虑粘性效应的求解方法，从而有效避免传统基于动网格单元的数值模拟方法中离散方程的无效计算。其次，在动网格单元的坐标、面向量、体积、网格单元速度等变量的更新中引入在非定常效应主导区域内局部更新的思路，从而有效避免传统基于动网格单元的数值模拟方法中更新网格单元参数的无效计算。

Description

一种飞行器机动过程模拟的动网格扰动域更新方法
技术领域
本发明涉及计算流体力学领域，尤其涉及一种飞行器机动过程模拟的动网格扰动域更新方法。
背景技术
飞行器在飞行任务执行过程中，需通过俯仰、滚转、偏航等运动不断调整飞行姿态，从而实现既定轨道飞行。当飞行器调整姿态时，气流分离、旋涡运动等流动现象将会出现由物体运动导致的非定常效应，产生远超定常状态预测的空气动力学效应，甚至产生危及飞行安全的颠覆性气动特性。因此，飞行器机动过程中的气动特性预测对飞行器关键子系统的设计至关重要。
计算流体力学数值模拟是飞行器机动过程气动特性预测的重要技术之一。在保证数值精度的同时尽可能提高数值模拟的计算效率一直是计算流体力学的研究热点；对于飞行器机动过程模拟这类非定常问题，提高数值模拟的计算效率更是工程应用中的难题。传统数值模拟方法普遍采用固定的预设计算域求解动态问题，却忽略了“实际所需计算域会随流场演化与求解收敛而不断变化”的特征，因而可能导致大量无效计算。为避免传统方法中的无效计算，专利文献 202010528415.1 已针对由流动自身导致非定常问题的数值模拟建立了一种名为“迭代推进扰动域更新方法”的新型加速技术。不过，在模拟飞行器机动过程时，由于飞行姿态或气动部件位置不断变化，数值模拟须进一步结合动网格技术方可正确预测气动特性。一方面，在与动网格技术高效结合时，该专利所提技术仍需考虑网格单元速度、几何守恒量等问题而进一步补充算法。另一方面，局部更新的计算策略也可在动网格单元更新中加以利用，从而挖掘出数值算法更高的计算效率。
发明内容
针对上述现有技术中的不足，本发明提出一种飞行器机动过程模拟的动网格扰动域更新方法。该方法，首先在流动控制方程求解中考虑网格单元速度和动网格几何守恒律的影响，在动网格技术的框架下建立仅对非定常未收敛受扰单元求解、仅在局部区域考虑粘性效应的求解方法，从而有效避免传统基于动网格单元的数值模拟方法中离散方程的无效计算。其次，在动网格单元的坐标、面向量、体积、网格单元速度等变量的更新中引入在非定常效应主导区域内局部更新的思路，从而有效避免传统基于动网格单元的数值模拟方法中更新网格单元参数的无效计算。
本发明的技术方案具体如下：
S1：数据读入：
分配存储网格坐标、守恒量的静态数组，并读入飞行器流场的计算网格坐标、边界条件、计算设置；
S2：计算初始化
计算并存储计算域内所有网格单元的网格单元参数，并将所有网格单元速度归零；
根据来流条件或给定流场为计算域内所有网格单元的守恒量赋值；
S3：建立动态计算域
建立对流、粘性、非定常三类动态计算域；
S4：在动态计算域内求解流动控制方程
S4-1：边界条件处理；
S4-2：残差估计；
S4-3：时间积分；
S5：动态计算域更新
S5-1：增大对流、非定常动态域；
S5-2：缩小对流动态域；
S5-3：增大粘性动态域；
S5-4：缩小粘性动态域；
S6：判断内迭代是否收敛；
判断动态计算域内所有网格单元的守恒量更新量模值是否均小于给定的收敛阈值，若满足，则继续步骤S7；若不满足，则跳转至步骤S4，进入当前时刻内迭代的下一个迭代步；
S7：判断计算是否完成
若当前时刻即为末态时刻，则执行步骤S11；若否，则继续步骤S8；
S8：更新动网格单元
在非定常动态域内，采用动网格技术移动网格单元节点，并更新网格单元的格心坐标、体积、面向量、网格单元速度和距飞行器壁面最小距离；
S9：时间推进中更新动态计算域
S9-1：缩小非定常动态域；
S9-2：重置对流、粘性动态域；
S9-3：分配存储空间；
S10：返回步骤S4，进入下一时刻；
S11：输出结果。
在一些实施例中，优选的，所述步骤S2中的网格单元参数包括格心坐标、体积、面向量和距飞行器壁面最小距离。
在一些实施例中，优选的，所述网格单元参数的计算方式如下：
令x i,j,k 表示网格单元(I, J, K)标号最小节点的坐标，则网格单元(I, J, K)的格心坐标表示为
Figure 662307DEST_PATH_IMAGE001
网格单元(I, J, K)的面向量表示为
Figure 807112DEST_PATH_IMAGE002
网格单元(I, J, K)的体积|Ω|表示为
Figure 500261DEST_PATH_IMAGE003
距飞行器壁面最小距离通过壁面距有效单元计算方法获得。
在一些实施例中，优选的，所述步骤S2中，根据来流条件是将计算域内所有网格单元的守恒量赋为来流值，根据给定流场是按给定流场插值获得计算域内所有网格单元的守恒量。
在一些实施例中，优选的，所述步骤S3中，当步骤S2中根据来流条件初始化时，取紧邻飞行器壁面的10层网格单元为对流和非定常动态域的初始网格单元、紧邻飞行器壁面的1层网格单元为粘性动态域的初始网格单元；当步骤S2中根据给定流场初始化时，筛选守恒量与给定流场条件不同的网格单元为对流和非定常动态域的初始网格单元，对流动态域中的粘性效应主导网格单元为粘性动态域的初始网格单元。
在一些实施例中，优选的，所述步骤S4-1中的边界条件处理方法为：以虚网格方法处理计算域的边界条件。
在一些实施例中，优选的，所述步骤S4-2中的残差估计方法如下：
动网格单元上，基于双时间步法离散的非定常流动控制方程为
Figure 30600DEST_PATH_IMAGE004
其中，R(W*)表示残差，定义为
Figure 783792DEST_PATH_IMAGE005
Figure 914428DEST_PATH_IMAGE006
Figure 411268DEST_PATH_IMAGE007
动网格单元上的对流通量
Figure 327272DEST_PATH_IMAGE008
定义为
Figure 985786DEST_PATH_IMAGE009
Figure 354451DEST_PATH_IMAGE010
其中，W表示网格单元的守恒量；τ表示虚拟时间，t表示物理时间；|Ω|、N f、ΔS分别表示网格单元的体积、网格单元面的面数与面积；F v表示粘性通量；Q T表示湍流模型方程源项；ρ、p、H分别表示密度、压力、总焓，ui表示速度分量，ni表示网格单元面外法向分量；V t表示网格单元速度在网格单元面外法向上的分量，V r表示网格单元面外法向上的逆变速度，上标“*”代表当前时刻内迭代的求解值，“n+1”、“n”、“n-1”分别代表当前时刻、上一时刻和上上一时刻，
Figure 402785DEST_PATH_IMAGE011
为无粘项，仅在对流动态域中计算；
Figure 173295DEST_PATH_IMAGE012
为粘性项，仅在粘性动态域中计算。
在一些实施例中，优选的，所述步骤S4-3中，在对流动态域中，采用时间格式离散虚拟时间项，并更新W*。
在一些实施例中，优选的，所述步骤S5-1具体包括：遍历对流动态域的边界单元，判断每一个边界单元是否受到扰动，若受到扰动，则将该边界单元周围所有紧邻单元加入对流和非定常动态域。
在一些实施例中，优选的，通过边界单元守恒量更新量的模值，判断该边界单元是否受到扰动。
在一些实施例中，优选的，所述步骤S5-2具体包括：遍历对流动态域的边界单元，将满足以下5个条件的边界单元从对流动态域中移除，若该边界单元也属于粘性动态域，则也将其从粘性动态域中移除；（1）周围不存在新增无粘受扰单元；（2）求解已收敛；（3）位于最上游；（4）不再影响对流动态域中其他单元的求解；（5）不再受对流动态域中其他单元的影响。
在一些实施例中，优选的，所述步骤S5-3具体包括：遍历粘性动态域的边界单元，判断每一个边界单元是否受到扰动，若受到扰动，则将该边界单元的所有紧邻单元加入粘性动态域。
在一些实施例中，优选的，通过边界单元的守恒量判断该边界单元是否受到粘性扰动。
在一些实施例中，优选的，所述步骤S5-4具体包括：遍历粘性动态域的边界单元，将满足以下2个条件的边界单元从粘性动态域中移除；（1）周围不存在新增的粘性效应主导单元；（2）不再是粘性效应主导单元。
在一些实施例中，优选的，所述步骤S9-1具体包括：遍历非定常动态域的边界单元，将满足以下2个条件的边界单元从非定常动态域中移除；（1）不再是非定常效应主导单元；（2）位于最上游。
在一些实施例中，优选的，所述步骤S9-2具体包括：根据非定常动态域范围，重新建立对流和粘性动态域；其中，对流动态域范围与非定常动态域相同，粘性动态域包含非定常动态域内所有粘性效应主导单元。
在一些实施例中，优选的，所述步骤S9-3具体包括：按照非定常动态域的范围，为需要更新的变量分配存储空间，并更新W (n)和W (n-1)。
相比于现有技术，本发明具有以下有益效果：
在模拟飞行器机动过程的传统数值方法中，为便于程序实现，往往不考虑非定常流动实际所需求解的区域，而是采用全局更新所有计算网格单元的策略，从而可能导致大量无效计算。所产生的无效计算主要源自两方面：其一，离散控制方程求解中的无效计算；其二，网格单元参数更新中的无效计算。本发明提出的飞行器机动过程的动网格扰动域更新方法可显著降低飞行器机动过程数值模拟的计算量，提升飞行器机动过程数值模拟的计算效率，主要体现在以下几个方面：
（1）针对第一类无效计算，本发明建立了对流、粘性、非定常等三类动态计算域。其中，仅在对流动态域中计算流动控制方程残差的无粘项，仅在粘性动态域内计算残差的粘性项，仅在对流动态域中更新当前时刻的守恒量更新量，从而有效的减少了离散控制方程求解中的无效计算。
（2）针对第二类无效计算，本发明仅在非定常动态域内采用动网格技术更新网格节点，并只重新计算非定常动态域内单元的格心坐标、体积、面向量、网格单元速度和距飞行器壁面最小距离等网格单元参数。且在更新动态计算域时，确保三类动态计算域在数值模拟每个物理时刻、每个内迭代步中的范围都保持在保证正确计算的最小求解范围，从而有效的减少了网格单元参数更新中的无效计算。
附图说明
通过参考附图可更好地理解本发明。图中的构件不应视作按比例绘制，重点应放在示出本发明的原理上。
图1为本发明动网格扰动域更新方法的流程图；
图2为本发明动网格扰动域更新方法模拟NACA0012翼型俯仰运动所生成的动网格；
图3为本发明动网格扰动域更新方法与传统方法模拟跨声速NACA0012翼型俯仰运动问题的升力系数对比；
图4为本发明动网格扰动域更新方法与传统方法模拟超声速导弹俯仰运动问题的气动力系数对比；
图5为本发明动网格扰动域更新方法模拟超声速导弹俯仰运动问题的动态计算域网格量变化曲线。
具体实施方式
以下结合说明书附图和具体实施例来进一步说明本发明，以使本领域的技术人员可以更好地理解本发明并能予以实施，但本发明并不局限于此，凡在本发明的精神和原则之内所做的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
如图1所示，本发明的具体实施例中提供了一种飞行器机动过程的动网格扰动域更新方法的流程，具体包括如下步骤：
S1：数据读入
分配存储网格坐标、守恒量的静态数组，并读入飞行器流场的计算网格坐标、边界条件、计算设置；
S2：计算初始化
计算并存储计算域内所有网格单元的格心坐标、体积、面向量和距飞行器壁面最小距离，所有单元网格单元速度归零。
令x i,j,k 表示单元(I, J, K)标号最小节点的坐标，则网格单元(I, J, K)的格心坐标表示为
Figure 2710DEST_PATH_IMAGE013
（1）
网格单元(I, J, K)的面向量表示为
Figure 858671DEST_PATH_IMAGE014
（2）
单元(I, J, K)的体积V表示为
Figure 946582DEST_PATH_IMAGE015
（3）
距飞行器壁面最小距离通过壁面距有效单元计算方法获得。
根据来流条件或给定流场为计算域内所有网格单元的守恒量赋值。根据来流条件初始化是将计算域内所有单元的守恒量赋为来流值。根据给定流场初始化是按给定流场插值获得计算域内所有单元的守恒量。
S3：建立动态计算域
根据初始化方式，分别建立对流、粘性、非定常三类动态计算域。当根据来流条件初始化时，取紧邻壁面的10层单元为对流和非定常动态域的初始单元、紧邻壁面的1层单元为粘性动态域的初始单元。当根据给定流场初始化时，筛选守恒量与来流条件不同的单元为对流和非定常动态域的初始单元，对流动动态域中的粘性效应主导单元为粘性动态域的初始单元。
S4：在动态计算域内求解流动控制方程
S4-1：边界条件处理
以虚网格方法处理计算域的边界条件；其中，对称边界计算壁面法向的逆变速度时需考虑网格单元速度，壁面边界处的流动速度应为网格单元速度。
S4-2：残差估计
动网格单元上，基于飞行器机动过程常用方法双时间步法离散的非定常流动控制方程，表示为
Figure 571598DEST_PATH_IMAGE016
（4）
其中，R(W *)表示残差，定义为
Figure 837494DEST_PATH_IMAGE017
Figure 915172DEST_PATH_IMAGE018
Figure 823085DEST_PATH_IMAGE019
（5）
动网格单元上的对流通量
Figure 53340DEST_PATH_IMAGE020
定义为
Figure 490138DEST_PATH_IMAGE021
（6）
Figure 320690DEST_PATH_IMAGE022
（7）
其中，W表示守恒量；τ表示虚拟时间，t表示物理时间；|Ω|、N f、ΔS分别表示网格单元的体积、网格单元面的面数与面积；F v表示粘性通量；Q T表示湍流模型方程源项；ρ、p、H分别表示密度、压力、总焓，u i 表示速度分量，n i 表示网格单元面外法向分量；V t表示网格单元速度在网格单元面外法向上的分量，V r表示网格单元面外法向上的逆变速度，V r=n i u i-V t；上标“*”代表当前时刻内迭代的求解值，“n+1”、“n”、“n-1”分别代表当前时刻、上一时刻和上上一时刻。本发明将残差中的对流通量项、二阶后向差分近似物理时间导数产生的源项，将其中的粘性通量项、湍流模型方程源项统称为粘性项，即
Figure 501136DEST_PATH_IMAGE023
。无粘项
Figure 100745DEST_PATH_IMAGE024
仅在对流动态域中计算，粘性项仅在粘性动态域中计算。
S4-3：时间积分
在对流动态域中，采用时间格式离散虚拟时间项，并更新W *。
S5：动态计算域更新
S5-1：增大对流、非定常动态域
遍历对流动态域的边界单元，对任一边界单元：
（1）通过守恒量更新量的模值，判断该边界单元是否受到扰动；令||ΔW||表示守恒量更新量的模值，ε a,c表示给定的对流新增阈值，则该边界单元受到扰动可描述为
Figure 957711DEST_PATH_IMAGE025
（8）
（2）若该边界单元满足式（8），则将其可能受到扰动的紧邻单元加入对流和非定常动态域。对于处于亚声速流动中的单元，其可能受到扰动的紧邻单元是四周所有单元。对于处于超声速流动中的单元，令q为该单元格心指向该单元格点的单位向量，则扰动沿q方向传播将满足
Figure 275560DEST_PATH_IMAGE026
（9）
式中，u代表流动速度矢量，a为声速。若该单元格心指向某一格点的单位向量满足上式，则将与包含该格点的紧邻单元加入对流和非定常动态域。
S5-2：缩小对流动态域
遍历对流动态域的边界单元，对任一边界单元，若其满足以下5个条件则将该单元从对流动态域中移除；若该单元也属于粘性动态域，则也将其从粘性动态域中移除。具体判断条件如下：
（1）周围不存在新增无粘受扰单元
新增无粘受扰单元的紧邻单元满足式（8），因此若周围单元均不满足式（8）即表示待删除单元周围不存在新增无粘受扰单元。
（2）求解已收敛
网格单元的收敛条件可通过守恒量更新量的模值描述。令ε d表示给定的删除阈值，则待删除单元求解已收敛可描述为
Figure 994117DEST_PATH_IMAGE027
（10）
（3）位于最上游
通过流动方向及单元的几何关系即可判断是否位于最上游。令q表示待删除单元格心指向其紧邻单元格心的单位向量，若待删除单元位于最上游，则其所有处于对流动态域中的紧邻单元均应满足
Figure 713812DEST_PATH_IMAGE028
（11）
式中，θ d表示上游单元容差角，超声速流动取10°，亚声速流动取45°。
（4）不再影响对流动态域中其他单元的求解
数值实验表明，马赫数小于0.3的不可压缩流动内，上、下游流动的收敛没有先后顺序。因此，马赫数大于0.3的可压缩流动位于最上游，则可不再影响对流动态域中其他单元的求解；而马赫数小于0.3的不可压缩流动，则不可移除。
（5）不再受对流动态域中其他单元的影响
对于超声速无粘流动中，由于控制方程的数学性质呈双曲型，即流场中任一点流动不受其下游流动的影响，因此位于最上游的超声速无粘单元也自然不再受对流动态域中其他单元的影响。
对于亚声速和粘性流动，若考虑紧邻单元对待删除单元守恒量更新量的影响后，其仍满足收敛条件，则可认为该单元不再受其他单元的影响，即满足
Figure 492412DEST_PATH_IMAGE029
（12）
式中，Δt表示迭代步长，C CFL表示时间推进格式的CFL数，|Ω|表示网格单元的体积；I, J, K表示网格单元方向；ΔR i 表示对流动态域中的紧邻单元对对流动态域边界单元的残差沿i方向的影响，i=I, J, K。
对于亚声速无粘单元，紧邻单元对对流动态域边界单元的残差沿网格单元方向的影响ΔR i 可表示为
Figure 297557DEST_PATH_IMAGE030
Figure 591045DEST_PATH_IMAGE031
（13）
式中，ΔW表示守恒量更新量的模值；DF c表示对流通量变化量，即当前步与上一步对流通量之差；下标i±1表示对流动态域边界单元沿正、负i方向的紧邻单元；下标i±1/2表示对流动态域边界单元沿正、负i方向的单元面；
Figure 899667DEST_PATH_IMAGE032
表示对流通量Jacobian矩阵沿i方向的谱半径。
对于粘性单元
Figure 114747DEST_PATH_IMAGE033
Figure 876030DEST_PATH_IMAGE034
（14）
S5-3：增大粘性动态域
遍历粘性动态域的边界单元，对任一边界单元：
（1）通过守恒量，判断该边界单元是否受到粘性扰动；令ε a,v表示给定的粘性新增单元，粘性效应主导单元应满足：
Figure 467548DEST_PATH_IMAGE035
（15）
式中，Ψ为粘性效应衡量参数，
Figure 879944DEST_PATH_IMAGE036
为第1个迭代步紧邻单元粘性效应衡量参数的最大值。
（2）若受到粘性扰动，则将该边界单元的所有紧邻单元加入粘性动态域。
S5-4：缩小粘性动态域
遍历粘性动态域的边界单元，若其满足以下2个条件，则将该单元从粘性动态域中移除。具体条件如下：
（1）周围不存在新增的粘性效应主导单元；
判断相邻两层单元是否满足式（19），若均不满足则表示周围不存在新增的粘性效应主导单元；
（2）不再是粘性效应主导单元；
若待删除单元不满足式（19），则表示其不再是粘性效应主导单元。
S6：判断内迭代是否收敛
通过所有网格单元守恒量更新量模值的最大值判断当前时刻内迭代是否收敛；若该值小于给定的收敛阈值，则当前时刻的求解已收敛，继续步骤S7；若该值大于给定的收敛阈值，则跳转至步骤S4，进入当前时刻内迭代的下一个迭代步。
S7：判断计算是否完成
若当前时刻即为末态时刻，则继续步骤S11；若否，则继续步骤S8。
S8：更新动网格单元；
在非定常动态域内，采用动网格技术移动网格单元节点。并更新移动单元的格心坐标、体积、面向量、网格单元速度和距飞行器壁面最小距离等网格单元参数。更新移动单元的格心坐标、体积可分别采用式（1）、（3）计算。当采用动网格单元数值模拟时，须考虑几何守恒率以避免网格单元变形而产生的误差。几何守恒律要求：网格单元变形过程中控制体的体积增量应与面运动形成的体积相等。由流动控制方程中的质量守恒方程可知，
Figure 265926DEST_PATH_IMAGE037
（16）
式中，u表示流动速度，g表示网格单元速度。
对于密度均匀、流速均匀的封闭控制体，质量守恒方程可化简为
Figure 780084DEST_PATH_IMAGE038
（17）
若以二阶后向差分离散上式的时间项，则单元面m在n+1时刻的网格单元速度可近似为
Figure 909714DEST_PATH_IMAGE039
Figure 927349DEST_PATH_IMAGE040
（18）
式中，下标“m”表单元面；
Figure 500543DEST_PATH_IMAGE041
表单元面m上四个顶点在时刻n和n+1所在位置所组成六面体的体积；
Figure 236418DEST_PATH_IMAGE042
为修正的单元面面积，定义为
Figure 169739DEST_PATH_IMAGE043
（19）
S9：时间推进中更新动态计算域
S9-1：缩小非定常动态域
遍历非定常动态域的边界单元，若其满足以下2个条件，则将该单元从非定常动态域中移除。具体条件如下：
（1）不再是非定常效应主导单元；
若单元已不受非定常效应主导，则其物理时间导数应趋近于零。因此，单元不再受非定常效应主导的条件可表示为
Figure 41880DEST_PATH_IMAGE044
Figure 35244DEST_PATH_IMAGE045
（20）
（2）位于最上游，可由式（11）判断。
S9-2：重置对流、粘性动态域
根据非定常动态域范围，重新建立对流、粘性动态域；其中，对流动态域范围与非定常动态域相同，粘性动态域应包含非定常动态域内所有粘性效应主导单元。
S9-3：分配存储空间
按照非定常动态域的范围为存储守恒量更新量、当地时间步长等与更新相关的变量分配空间，并更新W (n)、W (n-1)。
S10：返回步骤S4，进入下一时刻；
S11：输出结果，结束。
实施例1：
采用本发明的方法对马赫数0.755的NACA0012翼型俯仰运动问题进行模拟。图2展示了本发明模拟该实施例在攻角最大和最小状态时的网格，其中α表示迎角。图中，由于翼型做俯仰运动，数值模拟中的计算网格坐标也随翼型一起运动。图3将本发明所得升力曲线与传统方法的数值结果、风洞实验数据进行了对比。对比结果表明，本发明所得非定常气动力结果与传统方法、风洞实验吻合良好，相对气动力偏差均低于10-3。而相比于传统方法，本发明由于有效避免了传统方法中的无效计算，可缩短27%的计算时间，具有显著的加速效果。
实施例2：
采用本发明的方法对马赫数1.58的有翼导弹俯仰运动进行模拟。图4对比了本发明与传统方法所获得的有翼导弹气动力系数。图中，两种方法所得结果完全吻合，所得阻力系数相对偏差至多仅10-5量级。图5展示了本发明在求解该实施例中的动态域网格单元量变化曲线；其中，η u表示非定常动态域与预设计算域的网格单元量之比，η c,min表示内迭代达到收敛时的对流动态域与预设计算域的网格单元量之比。图中，非定常动态域的范围会随导弹的俯仰运动呈周期性变化；而在每个时刻的内迭代步中，对流动态域都会从非定常动态域的范围缩小至壁面附近。得益于动态计算域的变化，本发明相比于传统方法可节省39%的计算时间，加速效果明显。
应理解，前述仅说明了一些实施方式，可进行改变、修改、增加和/或变化而不偏离所公开的实施方式的范围和实质，该实施方式是示意性的而不是限制性的。此外，所说明的实施方式涉及当前考虑为最实用和最优选的实施方式，其应理解为实施方式不应限于所公开的实施方式，相反地，旨在覆盖包括在该实施方式的实质和范围内的不同的修改和等同设置。此外，上述说明的多种实施方式可与其它实施方式共同应用，如，一个实施方式的方面可与另一个实施方式的方面结合而实现再另一个实施方式。另外，任何给定组件的各独立特征或构件可构成另外的实施方式。
以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围，其均应涵盖在本发明的权利要求和说明书的范围当中。

Claims (6)

1.一种飞行器机动过程模拟的动网格扰动域更新方法，包括以下步骤：
S1：数据读入
分配存储网格坐标、守恒量的静态数组，并读入飞行器流场的计算网格坐标、边界条件、计算设置；
S2：计算初始化
计算并存储计算域内所有网格单元的网格单元参数，并将所有网格单元速度归零；
根据来流条件或给定流场为计算域内所有网格单元的守恒量赋值；
S3：建立动态计算域
建立对流、粘性、非定常三类动态计算域；
S4：在动态计算域内求解流动控制方程
S4-1：边界条件处理；
S4-2：残差估计；
S4-3：时间积分；
S5：动态计算域更新
S5-1：增大对流、非定常动态域；
S5-2：缩小对流动态域；
S5-3：增大粘性动态域；
S5-4：缩小粘性动态域；
S6：判断内迭代是否收敛；
判断动态计算域内所有网格单元的守恒量更新量模值是否均小于给定的收敛阈值，若满足，则继续步骤S7；若不满足，则跳转至步骤S4，进入当前时刻内迭代的下一个迭代步；
S7：判断计算是否完成
若当前时刻即为末态时刻，则执行步骤S11；若否，则继续步骤S8；
S8：更新动网格单元
在非定常动态域内，采用动网格技术移动网格单元节点，并更新网格单元的格心坐标、体积、面向量、网格单元速度和距飞行器壁面最小距离；
S9：时间推进中更新动态计算域
S9-1：缩小非定常动态域；
S9-2：重置对流、粘性动态域；
S9-3：分配存储空间；
S10：返回步骤S4，进入下一时刻；
S11：输出结果；
所述步骤S2中的网格单元参数包括格心坐标、体积、面向量和距飞行器壁面最小距离；
所述网格单元参数的计算方式如下：
令x i,j,k 表示网格单元(I, J, K)标号最小节点的坐标，则网格单元(I, J, K)的格心坐标表示为
Figure 585243DEST_PATH_IMAGE001
网格单元(I, J, K)的面向量表示为
Figure 161718DEST_PATH_IMAGE002
单元(I, J, K)的体积V表示为
Figure 88085DEST_PATH_IMAGE003
距飞行器壁面最小距离通过壁面距有效单元计算方法获得；
所述步骤S3中，当步骤S2中根据来流条件初始化时，取紧邻飞行器壁面的10层网格单元为对流和非定常动态域的初始网格单元、紧邻飞行器壁面的1层网格单元为粘性动态域的初始网格单元；当根据给定流场初始化时，筛选守恒量与给定流场条件不同的网格单元为对流和非定常动态域的初始网格单元，对流动态域中的粘性效应主导网格单元为粘性动态域的初始网格单元。
2.根据权利要求1所述的飞行器机动过程模拟的动网格扰动域更新方法，其特征在于，所述步骤S2中，根据来流条件是将计算域内所有网格单元的守恒量赋为来流值，根据给定流场是按给定流场插值获得计算域内所有网格单元的守恒量。
3.根据权利要求1所述的飞行器机动过程模拟的动网格扰动域更新方法，其特征在于，
所述步骤S4-1中的边界条件处理方法为：以虚网格方法处理计算域的边界条件；
所述步骤S4-2中的残差估计方法如下：
动网格单元上，基于双时间步法离散的非定常流动控制方程为
Figure 685813DEST_PATH_IMAGE004
其中，R(W*)表示残差，定义为
Figure 39434DEST_PATH_IMAGE005
Figure 572047DEST_PATH_IMAGE006
Figure 302105DEST_PATH_IMAGE007
动网格单元上的对流通量
Figure 252875DEST_PATH_IMAGE008
定义为
Figure 511818DEST_PATH_IMAGE009
Figure 797306DEST_PATH_IMAGE010
其中，W表示网格单元的守恒量；τ表示虚拟时间，t表示物理时间；|Ω|、N f、ΔS分别表示网格单元的体积、网格单元面的面数与面积；F v表示粘性通量；Q T表示湍流模型方程源项；ρ、p、H分别表示密度、压力、总焓，ui表示速度分量，ni表示网格单元面外法向分量；V t表示网格单元速度在网格单元面外法向上的分量，V r表示网格单元面外法向上的逆变速度，上标“*”代表当前时刻内迭代的求解值，“n+1”、“n”、“n-1”分别代表当前时刻、上一时刻和上上一时刻，
Figure 65476DEST_PATH_IMAGE011
为无粘项，仅在对流动态域中计算；
Figure 120020DEST_PATH_IMAGE012
为粘性项，仅在粘性动态域中计算；
所述步骤S4-3中，在对流动态域中，采用时间格式离散虚拟时间项，并更新W*。
4.根据权利要求1所述的飞行器机动过程模拟的动网格扰动域更新方法，其特征在于，
所述步骤S5-1具体包括：遍历对流动态域的边界单元，判断每一个边界单元是否受到扰动，若受到扰动，则将该边界单元周围所有紧邻单元加入对流和非定常动态域；
所述步骤S5-2具体包括：遍历对流动态域的边界单元，将满足以下5个条件的边界单元从对流动态域中移除，若该边界单元也属于粘性动态域，则也将其从粘性动态域中移除；（1）周围不存在新增无粘受扰单元；（2）求解已收敛；（3）位于最上游；（4）不再影响对流动态域中其他单元的求解；（5）不再受对流动态域中其他单元的影响；
所述步骤S5-3具体包括：遍历粘性动态域的边界单元，判断每一个边界单元是否受到扰动，若受到扰动，则将该边界单元的所有紧邻单元加入粘性动态域；
所述步骤S5-4具体包括：遍历粘性动态域的边界单元，将满足以下2个条件的边界单元从粘性动态域中移除；（1）周围不存在新增的粘性效应主导单元；（2）不再是粘性效应主导单元。
5.根据权利要求4所述的飞行器机动过程模拟的动网格扰动域更新方法，其特征在于，所述步骤S5-1中通过边界单元守恒量更新量的模值，判断该边界单元是否受到扰动；所述步骤S5-3中通过边界单元的守恒量判断该边界单元是否受到粘性扰动。
6.根据权利要求1所述的飞行器机动过程模拟的动网格扰动域更新方法，其特征在于，
所述步骤S9-1具体包括：遍历非定常动态域的边界单元，将满足以下2个条件的边界单元从非定常动态域中移除；（1）不再是非定常效应主导单元；（2）位于最上游；
所述步骤S9-2具体包括：根据非定常动态域范围，重新建立对流和粘性动态域；其中，对流动态域范围与非定常动态域相同，粘性动态域包含非定常动态域内所有粘性效应主导单元；
所述步骤S9-3具体包括：按照非定常动态域的范围，为需要更新的变量分配存储空间，并更新W (n)和W (n-1)。
CN202210144150.4A 2022-02-17 2022-02-17 一种飞行器机动过程模拟的动网格扰动域更新方法 Active CN114218878B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210144150.4A CN114218878B (zh) 2022-02-17 2022-02-17 一种飞行器机动过程模拟的动网格扰动域更新方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210144150.4A CN114218878B (zh) 2022-02-17 2022-02-17 一种飞行器机动过程模拟的动网格扰动域更新方法

Publications (2)

Publication Number Publication Date
CN114218878A CN114218878A (zh) 2022-03-22
CN114218878B true CN114218878B (zh) 2022-05-10

Family

ID=80709286

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210144150.4A Active CN114218878B (zh) 2022-02-17 2022-02-17 一种飞行器机动过程模拟的动网格扰动域更新方法

Country Status (1)

Country Link
CN (1) CN114218878B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114841095B (zh) * 2022-07-05 2022-09-09 北京航空航天大学 一种不可压缩流动的扰动域推进方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111859530A (zh) * 2020-06-11 2020-10-30 北京航空航天大学 一种飞行器动态气动特性模拟的迭代推进扰动域更新方法
CN113850008A (zh) * 2021-12-02 2021-12-28 北京航空航天大学 飞行器气动特性预测的自适应网格扰动域更新加速方法
CN113962030A (zh) * 2021-12-20 2022-01-21 北京航空航天大学 飞行器多体分离模拟的重叠网格扰动域更新方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2014132138A2 (en) * 2013-02-07 2014-09-04 King Abdullah University Of Science And Technology Method and sensors for estimating and predicting airflow around air vehicles
US11551346B2 (en) * 2018-11-15 2023-01-10 Airbus (S.A.S.) Systems and methods of ultrasonic data evaluation of composite aircraft components
US20200410147A1 (en) * 2019-06-28 2020-12-31 Viettel Group Aerodynamic derivatives calculation method for flight vehicle
US11022524B2 (en) * 2019-11-05 2021-06-01 Dalian University Of Technology Stochastic configuration network based turbofan engine health parameter estimation method
CN111651831B (zh) * 2020-04-13 2022-04-08 北京航空航天大学 飞行器定常粘性可压缩绕流的分区扰动域更新计算方法
CN111859529B (zh) * 2020-06-11 2022-04-08 北京航空航天大学 一种飞行器绕流数值模拟的多重网格扰动域更新加速方法
CN113609597B (zh) * 2021-10-09 2022-01-21 北京航空航天大学 超声速飞行器绕流的时间-空间混合推进扰动域更新方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111859530A (zh) * 2020-06-11 2020-10-30 北京航空航天大学 一种飞行器动态气动特性模拟的迭代推进扰动域更新方法
CN113850008A (zh) * 2021-12-02 2021-12-28 北京航空航天大学 飞行器气动特性预测的自适应网格扰动域更新加速方法
CN113962030A (zh) * 2021-12-20 2022-01-21 北京航空航天大学 飞行器多体分离模拟的重叠网格扰动域更新方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Disturbance region update method for steady compressible flows;Shuyao Hu等;《Computer Physics Communications》;20180831;第229卷;全文 *
Zonal disturbance region update method for steady compressible viscous flows;Shuyao Hu等;《Computer Physics Communications》;20191130;第244卷;全文 *
可压缩湍流边界层燃烧减阻研究综述;刘宏鹏等;《空气动力学学报》;20200615;第38卷(第3期);全文 *
混流式转轮与导叶相互干涉的非定常数值模拟;杨昌明等;《农业机械学报》;20080225(第02期);全文 *

Also Published As

Publication number Publication date
CN114218878A (zh) 2022-03-22

Similar Documents

Publication Publication Date Title
CN111859530B (zh) 一种飞行器动态气动特性模拟的迭代推进扰动域更新方法
CN108563843B (zh) 定常可压缩流动的扰动区域更新方法
CN111859529B (zh) 一种飞行器绕流数值模拟的多重网格扰动域更新加速方法
CN108052772A (zh) 一种基于结构降阶模型的几何非线性静气动弹性分析方法
CN111651831B (zh) 飞行器定常粘性可压缩绕流的分区扰动域更新计算方法
CN113609597B (zh) 超声速飞行器绕流的时间-空间混合推进扰动域更新方法
CN113158338B (zh) 一种基于粗网格快速湍流壁面函数气动力预测方法
CN114444214B (zh) 一种基于舵面效率的飞行器控制方法
CN108363843A (zh) 一种基于结构降阶模型的几何非线性静气动弹性全机配平方法
CN113392472B (zh) 一种飞行器气动特性模拟的OpenMP并行扰动域更新方法
CN113609598B (zh) 飞行器气动特性模拟的rans/les扰动域更新方法
Marques et al. Multifidelity method for locating aeroelastic flutter boundaries
CN114218878B (zh) 一种飞行器机动过程模拟的动网格扰动域更新方法
WO2023168772A1 (zh) 一种飞行器动态特性模拟的时间并行扰动域更新方法
Bourgault-Côté et al. Multi-layer icing methodologies for conservative ice growth
Asada et al. FFVHC-ACE: fully automated Cartesian-grid-based solver for compressible large-eddy simulation
Goizueta et al. Parametric Krylov-based order reduction of aircraft aeroelastic models
CN115329458A (zh) 基于Fluent的海空多栖航行器空中飞行智能控制仿真方法
Murman et al. A vortex wake capturing method for potential flow calculations
Xiong et al. Study of Mach 0.8 Transonic Truss-Braced Wing Aircraft Wing-Strut Interference Effects
Görtz et al. Variable-fidelity and reduced-order models for aero data for loads predictions
Tantaroudas et al. An adaptive aeroelastic control approach using non linear reduced order models
Xiong et al. Simulations of Mach 0.8 Transonic Truss-Braced Wing Aircraft Aerodynamics at High Angles of Attack
Tamaki et al. Wall-Modeled Large-Eddy Simulation of Transonic Buffet over NASA-CRM Using FFVHC-ACE
Sang et al. An unstructured/structured multi‐layer hybrid grid method and its application

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