CN113850008B - 飞行器气动特性预测的自适应网格扰动域更新加速方法 - Google Patents

飞行器气动特性预测的自适应网格扰动域更新加速方法 Download PDF

Info

Publication number
CN113850008B
CN113850008B CN202111455975.XA CN202111455975A CN113850008B CN 113850008 B CN113850008 B CN 113850008B CN 202111455975 A CN202111455975 A CN 202111455975A CN 113850008 B CN113850008 B CN 113850008B
Authority
CN
China
Prior art keywords
grid
domain
unit
dynamic domain
units
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
CN202111455975.XA
Other languages
English (en)
Other versions
CN113850008A (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 CN202111455975.XA priority Critical patent/CN113850008B/zh
Publication of CN113850008A publication Critical patent/CN113850008A/zh
Application granted granted Critical
Publication of CN113850008B publication Critical patent/CN113850008B/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/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F21/00—Security arrangements for protecting computers, components thereof, programs or data against unauthorised activity
    • G06F21/60—Protecting data
    • G06F21/602—Providing cryptographic facilities or services
    • 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
    • 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)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Computer Security & Cryptography (AREA)
  • Software Systems (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Bioethics (AREA)
  • Automation & Control Theory (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种飞行器气动特性预测的自适应网格扰动域更新加速方法，首先，提出自适应网格与动态计算域耦合存储的新型数据格式，从而实现自适应网格上动态计算域信息的高效存储与遍历；其次，建立可标识具有位于脱体强激波区、收敛缓慢等特征单元的自适应探测器；再次，建立既可加速收敛又可最小化计算增量的脱体强激波区各向异性加密网格自适应策略；最后，建立自适应网格加密与动态计算域更新的耦合逻辑及其算法，提出自适应网格扰动域更新加速方法。该方法可进一步提升现有扰动域更新方法的加速效果；相比于传统的全局更新方法，可为飞行器气动特性预测提供同时优化网格与计算域的技术途径。

Description

飞行器气动特性预测的自适应网格扰动域更新加速方法
技术领域
本发明涉及计算流体力学领域，特别涉及飞行器气动特性预测的自适应网格扰动域更新加速方法。
背景技术
计算流体力学数值模拟已成为飞行器气动特性预测的重要手段之一。数值模拟的首要步骤便是指定计算域并将其离散为待求解单元，即生成计算网格。自适应网格技术便是一种根据几何信息与流动特性自动加密或稀疏计算网格的网格生成技术。在飞行器气动特性模拟中引入自适应网格技术，可在大流动梯度区保证数值精度的同时，显著降低计算网格的单元数量。现有基于自适应网格的数值模拟方法均采用静态计算域全局更新求解。由于忽略了数值/物理扰动传播与流动控制方程求解收敛的特征，这种全局更新的求解方式会产生不同程度的无效计算，会因对未受扰区域迭代求解而引入更多数值误差，会因对已收敛区域持续更新而延迟收敛。
专利文献 ZL201810250654.8 提出了一种扰动区域更新方法可通过消除传统方法中的无效计算而获得计算效率的显著提升。这种扰动区域更新方法正好可用于避免现有基于自适应网格数值模拟方法的不足。目前，扰动区域更新方法已在各类可压缩流动的数值模拟中得到了很好地验证，但其也存在计算效率进一步提升的空间。例如，该方法仅能通过优化计算域提升计算效率，却不具备优化网格质量的能力。一方面，网格分布不合理便无法以最小网格量获得最高数值精度。另一方面，若模拟流场中存在脱体强激波区，该区域便可能因网格质量未满足要求而导致动态计算域无法收缩，从而使扰动区域更新方法无法发挥其真正的加速效果。而自适应网格技术正好是在求解中通过改善网格疏密度实现收敛加速或计算精度提升的有效方法。因此，令两种方法优势互补，建立自适应网格技术与扰动区域更新方法的高效耦合方法是进一步提升数值模拟计算效率的有效途径。
发明内容
针对现有基于自适应网格技术的数值模拟方法存在无效计算、现有扰动区域更新方法无法自动优化网格限制加速效果的问题，本发明提出一种飞行器气动特性预测的自适应网格扰动域更新加速方法。通过扰动域区域更新方法解决现有基于自适应网格技术的数值模拟方法中存在无效计算的问题，通过自适应网格技术解决现有扰动区域更新方法中因无法自动优化网格而限制加速效果的问题，并着重提出两种方法耦合时所涉及的新算法与新数据格式。首先，提出自适应网格与动态计算域耦合存储的新型数据格式，从而实现自适应网格上动态计算域信息的高效存储与遍历。其次，建立可标识具有位于脱体强激波区、收敛缓慢等特征单元的自适应探测器。再次，建立既可加速收敛又可最小化计算增量的脱体强激波区各向异性加密网格自适应策略。最后，建立自适应网格加密与动态计算域更新的耦合逻辑及其算法，提出自适应网格扰动域更新加速方法。
一种飞行器气动特性预测的自适应网格扰动域更新加速方法，包括如下步骤：
S1：读入数据；
S2：粗化网格；
S3：初始化流场；
S4：建立动态计算域；
S5：求解控制方程；
S6：判断是否退出当前网格的迭代求解；
若当前网格为步骤S1读入网格的粗化网格，以最大迭代步数构建迭代求解退出判据；判断迭代步数是否达到最大迭代步数，若是，则跳转至步骤S8，若否，则继续步骤S7；
若当前网格即为步骤S1读入的网格，仍以最大迭代步数构建迭代求解退出判据；判断迭代步数是否达到最大迭代步数，若是，则跳转至步骤S9,若否，则继续步骤S7；
若当前网格为步骤S1读入网格的加密网格，以收敛阈值构建迭代求解退出判据；判断所有网格单元守恒量更新量的模值是否均小于收敛阈值，若是，则跳转至步骤S10，若否，则继续步骤S7；
S7：更新动态计算域；
S8：细化当前计算网格；
根据网格粗化准则，以步骤S1读入网格为基础，确定出比当前计算网格少稀疏一倍的计算网格；跳转至步骤S3；
S9：自适应加密网格；
依次进行标识待加密单元、加密标识单元、优化加密区域；跳转至步骤S3；
S10：输出结果，结束。
进一步的，所述步骤S1具体为：读入并存储数值模拟所需的数据，包括计算网格格点坐标、计算网格的边界条件、计算设置参数。
进一步的，所述步骤S2具体为：
将读入的计算网格按2的倍数稀疏，稀疏倍数N由用户指定；令m (0)表示读入的计算网格在某一网格方向上的单元数，若用户指定该方向稀疏N倍，则稀疏i倍后的单元数m (N)表示为
Figure 449012DEST_PATH_IMAGE001
（1）
式中，函数mod表示计算m (0)/2N的余数；
细网格标号ID(0)与粗网格标号ID(N)间的对应关系为
Figure 724266DEST_PATH_IMAGE002
（2）
式中，函数floor表示实数向下取整，函数dble表示将整数转化为实数；
粗网格标号与细网格标号间的对应关系为
Figure 806492DEST_PATH_IMAGE003
（3）
式中，函数ceiling表示向上取整。
进一步的，所述步骤S3具体为：
若由步骤S2进入，将网格单元的守恒量以来流值初始化；若由步骤S8或S9进入，则通过已完成计算的粗网格的守恒量为待求解的密网格初始化。
进一步的，所述步骤S4具体为：
若首次进入步骤S4，则根据网格的壁面边界分别建立对流动态域和粘性动态域；若非首次进入步骤S4，则根据初始化流场分别建立对流动态域与粘性动态域。
进一步的，所述步骤S5具体为：
离散流动控制方程表示为：
Figure 74662DEST_PATH_IMAGE004
（4）
其中，t表示时间，W为守恒量，F c、F v分别表示对流通量与粘性通量，Q c、Q v分别表示源项的无粘部分与粘性部分，
Figure 900446DEST_PATH_IMAGE005
表示网格单元的体积，
Figure 861449DEST_PATH_IMAGE006
表示单元面的面积，N f表示单元面的面数；当计算网格为步骤S1的读入网格及其粗化网格时，N f在二维和三维情况的取值分别为4、6；当计算网格为步骤S1读入网格的加密网格时，N f在二维和三维情况的取值将分别为不小于4和6的整数，具体取值由网格加密情况决定；
求解控制方程包含如下4个子步：
S51：调整存储空间——根据步骤S4或步骤S7所定义的对流动态域范围，调整当地时间步长、守恒量更新量的存储空间；
S52：边界条件处理——遍历步骤S1读入计算网格的所有边界，根据每一边界的条件，为其虚网格的守恒量赋值；
S53：残差估计；
估计式（4）等号右端项；式（4）中，对流通量项与源项无粘部分统称为控制方程残差的无粘项，粘性通量与源项粘性部分统称为控制方程残差的粘性项；遍历对流动态域中的所有单元，计算每一单元的残差无粘项；遍历粘性动态域中的所有单元，计算每一单元的残差粘性项；
在每一单元计算对流通量与粘性通量时，先采用重构格式确定单元界面的左、右状态值，再采用对流离散格式确定单元界面处的通量值；对于结构网格，重构格式按网格标号对模板单元插值；对于非结构网格，重构格式则根据相邻单元的流动梯度与间距进行插值；
S54：时间积分；
通过时间离散格式近似式（2）等号左端项，确定守恒量更新量并更新守恒量W。
进一步的，所述步骤S7包含如下4个子步：
S71：增大对流动态域；
遍历对流动态域的边界单元，对所有单元执行如下步骤：通过守恒量更新量的模值判断该边界单元是否受到扰动；若该边界单元受到扰动，则将其可能受到扰动的紧邻单元加入对流动态域；
S72：缩小对流动态域；
遍历对流动态域的边界单元，对所有单元判断如下条件：（1）周围不存在新增受扰单元，（2）求解已收敛，（3）位于最上游，（4）不再影响对流动态域中其他单元的求解；（5）不再受对流动态域中其他单元的影响；若上述5个条件均满足，则将该单元从对流动态域中移除；
S73：增大粘性动态域；
遍历粘性动态域的边界单元，通过守恒量判断该边界单元是否受到粘性扰动；若受到粘性扰动，则将其所有紧邻单元加入粘性动态域；
S74：缩小粘性动态域；
遍历粘性动态域的边界单元，若其满足以下2个条件：（1）周围不存在新增的粘性效应主导单元；（2）不再是粘性效应主导单元；则将该单元从粘性动态域中移除。
进一步的，所述步骤S8中采用式（1）定义的网格粗化准则。
进一步的，所述步骤S9包含如下3个子步：
S91：标识待加密单元；
遍历对流动态域中所有单元，若其满足以下2个条件：（1）位于最上游，（2）收敛缓慢，则将该单元标识为待加密单元；
S92：加密标识单元；
遍历步骤S91标记的所有待加密单元，以流动梯度为衡量标准，将流动梯度较大的方向作为网格加密方向；将被标记网格沿确定的加密方向二等分，并存储新增网格的格点坐标、格心坐标、面向量和壁面距；
S93：优化加密区域；
遍历步骤S92标记的所有被加密单元，根据被加密单元与其相邻单元的体积比与加密方向长度，对被加密单元的相邻单元进行加密，确保网格平滑过渡。
本发明与现有技术相比所具有的有益效果：
本发明所提飞行器气动特性预测的自适应网格扰动域更新方法可进一步提升现有扰动域更新方法的加速效果；相比于传统的全局更新方法，可为飞行器气动特性预测提供同时优化网格与计算域的技术途径。
本发明提高计算效率的原理主要体现在三个方面。其一，步骤S2将用户的输入网格粗化，使得本发明在根据用户输入的网格求解前，可先在粗化网格上求解，以便采用更大迭代步长使流场接近定常解。其二，步骤S5仅在动态计算域中求解流动控制方程，仅在对流动态域中估计残差的无粘项和时间积分，仅在粘性动态域中估计残差的粘性项，从而通过优化计算域有效消除了传统全局更新方法中的无效计算。其三，步骤S9通过自适应网格加密技术，针对对流动态域上游单元中的难以收敛单元，提高其在流动大梯度方向的网格分辨率。一方面，步骤S9可通过优化网格提高本发明方法的数值精度；另一方面，其又可改善难以收敛区域的收敛性，以避免因网格质量不佳影响动态计算域收缩从而影响扰动域更新方法加速效果的问题。
附图说明
图1为本发明自适应网格扰动域更新加速方法的流程图；
图2(a)为本发明对NACA0012翼型问题生成的自适应加密网格；
图2(b)为本发明对钝头体问题生成的自适应加密网格。
具体实施方式
图1为本发明实施的流程图，下面进行具体介绍。本发明提供了一种飞行器气动特性预测的自适应网格扰动域更新加速方法，具体包括如下步骤：
S1：数据读入；
分配网格格点坐标、网格边界条件和计算设置参数的存储空间，并从用户提供的输入文件中读入网格的节点坐标、边界条件，以及计算设置。
S2：粗化网格；
将读入的计算网格按2的倍数稀疏，稀疏倍数N由用户指定，通常系数2~3倍。令m (0)表示读入的计算网格在某一网格方向上的单元数，若用户指定该方向稀疏N倍，则稀疏N倍后的单元数m (N)可表示为
Figure 837495DEST_PATH_IMAGE007
（5）
式中，函数mod表示计算m (0)/2N的余数。
细网格标号ID(0)与粗网格标号ID(N)间的对应关系为
Figure 643777DEST_PATH_IMAGE008
（6）
式中，函数floor表示实数向下取整，函数dble表示将整数转化为实数；
粗网格标号与细网格标号间的对应关系为
Figure 83986DEST_PATH_IMAGE009
（7）
式中，函数ceiling表示向上取整。
S3：流场初始化；
若由步骤S2进入，流场还未被在任何网格上求解过，将网格单元的守恒量以来流值初始化。若由步骤S8或S9进入，则将已完成计算的粗网格的解插值到待求解的密网格上，为密网格上的流场求解进行初始化。
S4：建立动态计算域；
根据壁面边界建立动态计算域。取所有子网格中壁面边界的10层相邻单元作为对流动态域与非定常动态域的初始单元。
S5：求解控制方程；
离散流动控制方程可表示为：
Figure 153573DEST_PATH_IMAGE010
（8）
其中，t为时间，W为守恒量，F c、F v分别表示对流通量与粘性通量，Q c、Q v分别表示源项的无粘部分与粘性部分，
Figure 164385DEST_PATH_IMAGE011
表示网格单元的体积，
Figure 774358DEST_PATH_IMAGE012
表示单元面的面积，N f表示单元面的面数。当计算网格为结构网格时，即为由步骤S1读入的网格及其粗化网格时，N f在二维、三维情况的取值分别为4和6。当计算网格为步骤S1读入网格的加密网格时，N f在二维、三维情况的取值分别为不小于4和6的整数，具体取值由网格加密情况决定。
步骤S5包含如下4个子步：
S51：调整存储空间——根据步骤S4或步骤S7所定义的对流动态域范围，调整当地时间步长、守恒量更新量的存储空间；
S52：边界条件处理——遍历步骤S1读入计算网格的所有边界，根据计算网格的边界条件，为每一边界的虚网格的守恒量赋值；
S53：残差估计；
采用空间离散格式估计式（8）的等号右端项。首先，遍历对流动态域中的所有单元，计算每一单元在式（8）中的对流通量项与源项的无粘部分；其次，遍历粘性动态域中的所有单元，计算每一单元在式（8）中的粘性通量项与源项的粘性部分。在每一单元计算对流通量与粘性通量时，先采用重构格式确定单元界面的左、右状态值，再采用对流离散格式通过界面的左、右状态值确定出单元界面处的通量值。对于结构网格，重构格式按网格标号对模板单元插值；对于非结构网格，重构格式则根据相邻单元的流动梯度与间距进行插值。
S54：时间积分；
采用时间离散格式离散式（8）的等号左端项，确定守恒量更新量并更新守恒量W。
S6：判断是否退出当前网格的迭代求解；
若当前网格为步骤S1读入网格的粗化网格，以最大迭代步数构建迭代求解退出判据；判断迭代步数是否达到最大迭代步数，若是，则跳转至步骤S8，若否，则继续步骤S7；
若当前网格即为步骤S1读入的网格，仍以最大迭代步数构建迭代求解退出判据；判断迭代步数是否达到最大迭代步数，若是，则跳转至步骤S9,若否，则继续步骤S7；
若当前网格为步骤S1读入网格的加密网格，以收敛阈值构建迭代求解退出判据；判断所有网格单元守恒量更新量的模值是否均小于收敛阈值，若是，则跳转至步骤S10，若否，则继续步骤S7。
S7：更新动态计算域；
步骤S7包含如下4个子步：
S71：增大对流动态域；
遍历所有对流动态域的边界单元，对任一边界单元：
（1）通过守恒量更新量的模值，判断该边界单元是否受到扰动；令||ΔW||表示守恒量更新量的模值，ε a,c表示给定的对流新增阈值，取10-5，则||ΔW||>ε a,c代表单元已受到扰动；
（2）若该单元受到扰动，当处于亚声速流动中，将其所有紧邻单元加入对流动态域和非定常动态域；当处于超声速流动中，令q为该单元格心指向该单元格点的单位向量，u为流动速度矢量，a为声速，则扰动沿q方向传播的条件可表示为
Figure 803494DEST_PATH_IMAGE013
（9）
若该单元格心指向某一格点的单位向量满足式（9）则将与包含该格点的紧邻单元加入对流动态域和非定常动态域。
S72：缩小对流动态域；
遍历所有对流动态域的边界单元，对任一边界单元，若其满足以下5个条件则将该单元从对流动态域中移除。具体判断条件如下：
（1）周围不存在新增无粘受扰单元；
若周围单元均满足||ΔW||<ε a,c，则代表周围不存在新增无粘受扰单元。
（2）求解已收敛；
令ε d表示给定的删除阈值，取10-7，则待删除单元求解已收敛可描述为||ΔW||<ε d。
（3）位于最上游；
令q表示待删除单元格心指向其紧邻单元格心的单位向量，若待删除单元位于最上游，则其所有处于对流动态域中的紧邻单元均应满足
Figure 371879DEST_PATH_IMAGE014
（10）
式中，u表示速度矢量；θ d表示上游单元容差角，超声速流动取10°，亚声速流动取45°。
（4）不再影响对流动态域中其他单元的求解；
以马赫数为判断，马赫数大于0.3的可压缩流动位于最上游，则可不再影响对流动态域中其他单元的求解；而马赫数小于0.3的不可压缩流动，则不可移除。
（5）不再受对流动态域中其他单元的影响；
若考虑紧邻单元对待删除单元守恒量更新量的影响后，其仍满足收敛条件，则可认为该单元不再受其他单元的影响，即满足
Figure 604408DEST_PATH_IMAGE015
（11）
其中
Figure 221334DEST_PATH_IMAGE016
Figure 104976DEST_PATH_IMAGE017
（12）
式中，Δt表示迭代步长，C CFL表示时间推进格式的CFL数，|Ω|表示网格单元的体积；I, J, K表示网格方向；ΔR i 表示对流动态域中的紧邻单元对对流动态域边界单元的残差沿i方向的影响，i=I, J, K；ΔW表示守恒量更新量；
Figure 578683DEST_PATH_IMAGE018
表示对流通量变化量，即当前步与上一步对流通量之差；下标i±1表示对流动态域边界单元沿正、负i方向的紧邻单元；下标i±1/2表示对流动态域边界单元沿正、负i方向的单元面；
Figure 547776DEST_PATH_IMAGE019
表示对流通量Jacobian矩阵沿i方向的谱半径。
S73：增大粘性动态域；
遍历粘性动态域的边界单元，对任一边界单元：
（1）通过守恒量，判断该边界单元是否受到粘性扰动；令ε a,v表示给定的粘性新增单元，粘性效应主导单元应满足：
Figure 515863DEST_PATH_IMAGE020
（13）
式中，Ψ为粘性效应衡量参数，
Figure 519591DEST_PATH_IMAGE021
为第1个迭代步紧邻单元粘性效应衡量参数的最大值；
Figure 164199DEST_PATH_IMAGE022
表示式粘性通量Jacobian矩阵沿i方向的谱半径。
（2）若受到粘性扰动，则将该边界单元的所有紧邻单元加入粘性动态域。
S74：缩小粘性动态域；
遍历粘性动态域的边界单元，若其满足以下2个条件，则将该单元从粘性动态域中移除。具体条件如下：
（1）周围不存在新增的粘性效应主导单元；
判断相邻两层单元是否满足式（13），若均不满足则表示周围不存在新增的粘性效应主导单元；
（2）不再是粘性效应主导单元；
若删除单元不满足式（13），则表示其不再是粘性效应主导单元。
记录流场中速度梯度与压力梯度的平均值，将重叠区域中速度、压力梯度大于平均值的单元标记为大流动梯度单元，并记录其网格节点坐标。
S8：细化当前计算网格；
根据式（5）定义的网格粗化准则，以步骤S1读入网格为基础，确定出比当前计算网格少稀疏一倍的计算网格；
例如，若当前网格为读入网格稀疏N倍后的网格，则下一次计算的应为稀疏N-1倍的网格，即单元数为
Figure 620588DEST_PATH_IMAGE023
（14）
S9：自适应加密网格；
步骤S9包含如下3个子步：
S91：标识待加密单元；
遍历对流动态域中所有单元，若其满足以下2个条件，则将该单元标识为待加密单元：
（1）位于最上游；
是否位于最上游可通过式（10）判断；
（2）收敛缓慢；
计算该单元守恒量更新量随迭代步的下降率，若其小于0.5倍的最大下降率，则可认为其收敛缓慢；
S92：加密标识单元；
遍历步骤S91标记的所有待加密单元，以流动梯度为衡量标准，将流动梯度较大的方向作为网格加密方向；将被标记网格沿确定的加密方向二等分，并存储新增网格的格点坐标、格心坐标、面向量和壁面距；
S93：优化加密区域；
遍历步骤S92标记的所有被加密单元，根据被加密单元与其相邻单元的体积比与加密方向长度，对被加密单元的相邻单元进行加密，确保网格平滑过渡。
跳转至步骤S3。
S10：输出结果，结束。
下面结合附图和实施例进一步描述本发明，应该理解，以下所述实施例旨在便于对本发明的理解，而对其不起任何限定作用。
实施例1：图2(a)为本发明模拟NACA0012翼型问题时生成的自适应加密网格，图2(b)为本发明模拟钝头体问题时生成的自适应加密网格。图中，网格在近壁区依据壁面曲率实现了自适应加密，加密区域可保持光滑过渡。说明本发明方法可实现准确标识需要加密的区域，并合理过渡加密区域，为数值模拟提供分布合理的网格。
对于本领域的普通技术人员来说，在不脱离本发明创造构思的前提下，还可以对本发明的实施例做出若干变型和改进，这些都属于本发明的保护范围。

Claims (7)

1.一种飞行器气动特性预测的自适应网格扰动域更新加速方法，其特征在于，包括如下步骤：
S1：读入数据；
S2：粗化网格；
S3：初始化流场；
S4：建立动态计算域；
S5：求解流动控制方程；
S6：判断是否退出当前网格的迭代求解，包括：
若当前网格为步骤S1读入网格的粗化网格，以最大迭代步数构建迭代求解退出判据；判断迭代步数是否达到最大迭代步数，若是，则跳转至步骤S8，若否，则继续步骤S7；
若当前网格即为步骤S1读入的网格，仍以最大迭代步数构建迭代求解退出判据；判断迭代步数是否达到最大迭代步数，若是，则跳转至步骤S9,若否，则继续步骤S7；
若当前网格为步骤S1读入网格的加密网格，以收敛阈值构建迭代求解退出判据；判断所有网格单元守恒量更新量的模值是否均小于收敛阈值，若是，则跳转至步骤S10，若否，则继续步骤S7；
S7：更新动态计算域；
S8：细化当前计算网格，包括：
根据网格粗化准则，以步骤S1读入网格为基础，确定出比当前计算网格少稀疏一倍的计算网格；跳转至步骤S3；
S9：自适应加密网格；
依次进行标识待加密单元、加密标识单元、优化加密区域；跳转至步骤S3；
S10：输出结果，结束；
所述步骤S1具体为：读入并存储数值模拟所需的数据，包括计算网格格点坐标、计算网格的边界条件、计算设置参数；
所述步骤S2具体为：
将读入的计算网格按2的倍数稀疏，稀疏倍数N由用户指定；令m (0)表示读入的计算网格在某一网格方向上的单元数，若用户指定该方向稀疏N倍，则稀疏i倍后的单元数m (N)表示为
Figure 805126DEST_PATH_IMAGE001
（1）
式中，函数mod表示计算m (0)/2N的余数；
细网格标号ID(0)与粗网格标号ID(N)间的对应关系为
Figure 519004DEST_PATH_IMAGE002
（2）
式中，函数floor表示实数向下取整，函数dble表示将整数转化为实数；
粗网格标号与细网格标号间的对应关系为
Figure 361059DEST_PATH_IMAGE003
（3）
式中，函数ceiling表示向上取整。
2.根据权利要求1所述的飞行器气动特性预测的自适应网格扰动域更新加速方法，其特征在于，所述步骤S3具体为：
若由步骤S2进入S3，将网格单元的守恒量以来流值初始化；若由步骤S8或S9进入S3，则通过已完成计算的粗网格的守恒量为待求解的密网格初始化。
3.根据权利要求1所述的飞行器气动特性预测的自适应网格扰动域更新加速方法，其特征在于，所述步骤S4具体为：
若首次进入步骤S4，则根据网格的壁面边界分别建立对流动态域和粘性动态域；若非首次进入步骤S4，则根据初始化流场分别建立对流动态域与粘性动态域。
4.根据权利要求1所述的飞行器气动特性预测的自适应网格扰动域更新加速方法，其特征在于，所述步骤S5具体为：
离散流动控制方程，则流动控制方程表示为：
Figure 320050DEST_PATH_IMAGE004
（4）
其中，t表示时间，W为守恒量，F c、F v分别表示对流通量与粘性通量，Q c、Q v分别表示源项的无粘部分与粘性部分，
Figure 196739DEST_PATH_IMAGE005
表示网格单元的体积，
Figure 81518DEST_PATH_IMAGE006
表示单元面的面积，N f表示单元面的面数；当计算网格为步骤S1的读入网格及其粗化网格时，N f在二维和三维情况的取值分别为4、6；当计算网格为步骤S1读入网格的加密网格时，N f在二维和三维情况的取值将分别为不小于4和6的整数，具体取值由网格加密情况决定；
求解流动控制方程包含如下4个子步：
S51：调整存储空间——根据步骤S4或步骤S7所定义的对流动态域范围，调整当地时间步长、守恒量更新量的存储空间；
S52：边界条件处理——遍历步骤S1读入计算网格的所有边界，根据每一边界的条件，为其虚网格的守恒量赋值；
S53：残差估计；
式（4）等号右端项中，对流通量项与源项无粘部分统称为流动控制方程残差的无粘项，粘性通量与源项粘性部分统称为流动控制方程残差的粘性项；遍历对流动态域中的所有单元，计算每一单元的残差无粘项；遍历粘性动态域中的所有单元，计算每一单元的残差粘性项；
在每一单元计算对流通量与粘性通量时，先采用重构格式确定单元界面的左、右状态值，再采用对流离散格式确定单元界面处的通量值；对于结构网格，重构格式按网格标号对模板单元插值；对于非结构网格，重构格式则根据相邻单元的流动梯度与间距进行插值；
S54：时间积分；
通过式（4）等号左端项，确定守恒量更新量并更新守恒量W。
5.根据权利要求1所述的飞行器气动特性预测的自适应网格扰动域更新加速方法，其特征在于，所述步骤S7包含如下4个子步：
S71：增大对流动态域；
遍历对流动态域的边界单元，对所有单元执行如下步骤：通过守恒量更新量的模值判断该边界单元是否受到扰动；若该边界单元受到扰动，则将其可能受到扰动的紧邻单元加入对流动态域；
S72：缩小对流动态域；
遍历对流动态域的边界单元，对所有单元判断如下条件：（1）周围不存在新增受扰单元，（2）求解已收敛，（3）位于最上游，（4）不再影响对流动态域中其他单元的求解；（5）不再受对流动态域中其他单元的影响；若上述5个条件均满足，则将该单元从对流动态域中移除；
S73：增大粘性动态域；
遍历粘性动态域的边界单元，通过守恒量判断该边界单元是否受到粘性扰动；若受到粘性扰动，则将其所有紧邻单元加入粘性动态域；
S74：缩小粘性动态域；
遍历粘性动态域的边界单元，若其满足以下2个条件：（1）周围不存在新增的粘性效应主导单元；（2）不再是粘性效应主导单元；则将该单元从粘性动态域中移除。
6.根据权利要求1所述的飞行器气动特性预测的自适应网格扰动域更新加速方法，其特征在于，所述步骤S8中采用式（1）定义的网格粗化准则。
7.根据权利要求1所述的飞行器气动特性预测的自适应网格扰动域更新加速方法，其特征在于，所述步骤S9包含如下3个子步：
S91：标识待加密单元；
遍历对流动态域中所有单元，若其满足以下2个条件：（1）位于最上游，（2）收敛缓慢，则将该单元标识为待加密单元；
S92：加密标识单元；
遍历步骤S91标记的所有待加密单元，以流动梯度为衡量标准，将流动梯度大的方向作为网格加密方向；将被标记网格沿确定的加密方向二等分，并存储新增网格的格点坐标、格心坐标、面向量和壁面距；
S93：优化加密区域；
遍历步骤S92标记的所有被加密单元，根据被加密单元与其相邻单元的体积比与加密方向长度，对被加密单元的相邻单元进行加密，确保网格平滑过渡。
CN202111455975.XA 2021-12-02 2021-12-02 飞行器气动特性预测的自适应网格扰动域更新加速方法 Active CN113850008B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111455975.XA CN113850008B (zh) 2021-12-02 2021-12-02 飞行器气动特性预测的自适应网格扰动域更新加速方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111455975.XA CN113850008B (zh) 2021-12-02 2021-12-02 飞行器气动特性预测的自适应网格扰动域更新加速方法

Publications (2)

Publication Number Publication Date
CN113850008A CN113850008A (zh) 2021-12-28
CN113850008B true CN113850008B (zh) 2022-03-11

Family

ID=78982660

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111455975.XA Active CN113850008B (zh) 2021-12-02 2021-12-02 飞行器气动特性预测的自适应网格扰动域更新加速方法

Country Status (1)

Country Link
CN (1) CN113850008B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114218878B (zh) * 2022-02-17 2022-05-10 北京航空航天大学 一种飞行器机动过程模拟的动网格扰动域更新方法
CN114329799B (zh) * 2022-03-09 2022-06-17 北京航空航天大学 一种飞行器动态特性模拟的时间并行扰动域更新方法
CN114841095B (zh) * 2022-07-05 2022-09-09 北京航空航天大学 一种不可压缩流动的扰动域推进方法
CN114880971B (zh) * 2022-07-13 2022-09-20 中国空气动力研究与发展中心计算空气动力研究所 一种计算流体力学软件采用的隐式方法
CN117252129B (zh) * 2023-11-17 2024-02-20 中国空气动力研究与发展中心高速空气动力研究所 参数化的编队飞行气动干扰快速预测方法
CN117933144B (zh) * 2024-03-22 2024-05-28 中国空气动力研究与发展中心超高速空气动力研究所 一种求解复杂拓扑结构网格流场的多重网格方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102682173A (zh) * 2012-05-13 2012-09-19 北京理工大学 基于自适应径向基函数代理模型的飞行器优化设计方法
CN107273593A (zh) * 2017-06-01 2017-10-20 北京航空航天大学 一种用于高马赫数强激波流场气动热预测的湍流模型及其建立方法
CN107832494A (zh) * 2017-10-13 2018-03-23 南京航空航天大学 高超声速飞行器前缘流‑热‑固一体化计算方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20170116360A1 (en) * 2015-10-27 2017-04-27 Livermore Software Technology Corporation Efficient explicit finite element analysis of a product with a time step size control scheme

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102682173A (zh) * 2012-05-13 2012-09-19 北京理工大学 基于自适应径向基函数代理模型的飞行器优化设计方法
CN107273593A (zh) * 2017-06-01 2017-10-20 北京航空航天大学 一种用于高马赫数强激波流场气动热预测的湍流模型及其建立方法
CN107832494A (zh) * 2017-10-13 2018-03-23 南京航空航天大学 高超声速飞行器前缘流‑热‑固一体化计算方法

Also Published As

Publication number Publication date
CN113850008A (zh) 2021-12-28

Similar Documents

Publication Publication Date Title
CN113850008B (zh) 飞行器气动特性预测的自适应网格扰动域更新加速方法
JP5209298B2 (ja) 流れのシミュレーション計算方法およびシステム
CN111859530B (zh) 一种飞行器动态气动特性模拟的迭代推进扰动域更新方法
Wang A Quadtree-based adaptive Cartesian/Quad grid flow solver for Navier-Stokes equations
Fries et al. The extended/generalized finite element method: an overview of the method and its applications
Praveen et al. Low cost PSO using metamodels and inexact pre-evaluation: Application to aerodynamic shape design
CN111651831B (zh) 飞行器定常粘性可压缩绕流的分区扰动域更新计算方法
CN108563843B (zh) 定常可压缩流动的扰动区域更新方法
CN111859529B (zh) 一种飞行器绕流数值模拟的多重网格扰动域更新加速方法
CN113962030B (zh) 飞行器多体分离模拟的重叠网格扰动域更新方法
Jameson et al. Continuous adjoint method for unstructured grids
CN113609598B (zh) 飞行器气动特性模拟的rans/les扰动域更新方法
WO2012092502A1 (en) Generating inviscid and viscous fluid-flow simulations over an aircraft surface
CN113392472B (zh) 一种飞行器气动特性模拟的OpenMP并行扰动域更新方法
KR101328739B1 (ko) 형상 제어가 가능한 다상유체 시뮬레이션 장치 및 방법
Yamazaki et al. Geometry parameterization and computational mesh deformation by physics-based direct manipulation approaches
CN112613206B (zh) 一种基于各向异性体调和场的附面层网格生成方法
CN111046615B (zh) 一种黎曼解算器激波不稳定性抑制方法及系统
CN113609597B (zh) 超声速飞行器绕流的时间-空间混合推进扰动域更新方法
KR100972624B1 (ko) 다상유체 시뮬레이션 장치 및 방법
WO2014069421A1 (ja) シミュレーション装置、シミュレーション方法およびプログラム
CN113378440A (zh) 一种流固耦合数值模拟计算方法、装置及设备
Ali et al. Optimal multi-block mesh generation for CFD
JP3587827B2 (ja) 翼形性能の推定方法および翼形性能の推定プログラム
CN109726431B (zh) 一种基于平均核函数和迭代密度变化率的自适应sph流体模拟方法

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