CN113392472A - 一种飞行器气动特性模拟的OpenMP并行扰动域更新方法 - Google Patents

一种飞行器气动特性模拟的OpenMP并行扰动域更新方法 Download PDF

Info

Publication number
CN113392472A
CN113392472A CN202110940149.8A CN202110940149A CN113392472A CN 113392472 A CN113392472 A CN 113392472A CN 202110940149 A CN202110940149 A CN 202110940149A CN 113392472 A CN113392472 A CN 113392472A
Authority
CN
China
Prior art keywords
grid
dynamic
domain
boundary
dynamic domain
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
Application number
CN202110940149.8A
Other languages
English (en)
Other versions
CN113392472B (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 CN202110940149.8A priority Critical patent/CN113392472B/zh
Publication of CN113392472A publication Critical patent/CN113392472A/zh
Application granted granted Critical
Publication of CN113392472B publication Critical patent/CN113392472B/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
    • 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
    • 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]
    • 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)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (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

一种飞行器气动特性模拟的OpenMP并行扰动域更新方法,包括数据读入，静态区域分解，建立动态计算域，动态区域分解，分配存储空间，边界条件处理，残差估计，时间积分，块内增大对流动态域，块内缩小动态计算域，块内增大粘性动态域，块内缩小粘性动态域，块间增大动态计算域，判断求解是否达到收敛条件，结果输出等步骤，可以实现针对现有并行数值模拟方法在模拟飞行器气动特性时存在大量无效计算而现有扰动区域更新方法开展大规模并行时计算效率可能损失的问题。

Description

一种飞行器气动特性模拟的OpenMP并行扰动域更新方法
技术领域
本发明涉及计算流体力学领域，特别涉及一种飞行器气动特性模拟的OpenMP并行扰动域更新方法。
背景技术
计算流体力学数值模拟已改变了当代飞行器的设计流程，成为各类先进飞行器设计中不可或缺的关键技术。通过将现有技术与新技术高效融合，不断提升数值模拟的计算效率，对提高飞行器设计迭代效率、缩短研制周期具有重要的工程应用价值。
并行技术，是开展飞行器气动特性模拟必须使用的加速技术之一。并行技术的加速原理主要是通过区域分解策略将计算任务分解成多个子任务并行处理，从而通过降低串行处理的计算量实现单位时间的求解加速。对于飞行器气动特性的模拟，区域分解策略是将计算网格均匀分配给线程（进程），各线程（进程）同时以相同数值方法求解，并通过网格交界面的数据交换实现整个流场求解的统一。OpenMP是飞行器气动特性模拟中常用的CPU并行编程模式，以并行任务共享内存为特征，具有通信开销低、支撑细粒度循环级并行等优势。
目前，基于传统数值模拟方法发展的OpenMP并行计算主要存在两点影响计算效率的因素。其一，采用全局更新求解的计算策略，产出大量无效计算，从而严重影响计算效率。为避免传统数值模拟方法中的无效计算，专利文献ZL 201810250654.8已建立了名为“扰动区域更新方法”的新型加速技术。不过，在面对需大规模并行计算模拟的问题时，该方式所提技术仍需进一步补充算法并行化、数据交换量最小化等其他关键技术，否则将出现计算效率的明显损失。其二，采用按网格块存储数据、按网格块遍历数据的常规数据结构与循环结构，未利用OpenMP并行共享内存的特征，产生多余的内边界处理操作，从而影响计算效率。
发明内容
本发明的目的在于克服现有技术的不足，提供一种飞行器气动特性模拟的OpenMP并行扰动域更新方法，既可应对现有并行数值模拟方法在模拟飞行器气动特性时存在大量无效计算的问题，又可解决现有扰动区域更新方法开展大规模并行时计算效率可能损失的问题，还可基于新的数据结构与循环结构最少化内边界处理的工作量。
本发明提供了一种飞行器气动特性模拟的OpenMP并行扰动域更新方法，包括如下步骤：
S1：数据读入——分配网格坐标和块信息的存储空间，读入飞行器流场的多块结构化网格、边界条件、计算设置数据；
S2：静态区域分解——将读入的多块结构化网格合并至最小块数；根据该网格的单元数将网格分为网格子块并均匀分配给线程；
S3：流场初始化——按线程并行执行，各线程分配存储流动特性的存储空间，根据来流条件或给定流场对飞行器流场网格的所有单元赋初值；
S4：建立动态计算域——根据壁面边界或指定流场两种方式，建立对流、粘性两类动态计算域；
S5：动态区域分解——根据动态计算域单元数，将计算任务平均分配至各线程；
S6：分配存储空间——按线程并行，各线程根据所包含动态计算域的范围分配守恒量更新量、当地时间步长等与更新相关变量的存储空间；
S7：边界条件处理——按边界并行，为边界虚网格或边界面赋值；物理边界按照其物理定义为虚网格或边界面赋值，内边界需额外存储第一层虚网格的守恒量更新量模值；
S8：残差估计——将流动控制方程的残差分为无粘项和粘性项：按线程并行，在对流动态域中计算残差的无粘项，在粘性动态域中计算残差的粘性项；
S9：时间积分——按线程并行，在对流动态域中，求解守恒量更新量，更新流场变量；
S10：块内增大对流动态域——按网格块并行，将受扰单元加入对流动态域；
S11：块内缩小动态计算域——按网格行并行，判断某一对流动态域边界单元是否可从对流动态域中移除，并相应调整粘性动态域；
S12：块内增大粘性动态域——按网格行并行，将粘性效应主导单元加入粘性动态域；
S13：缩小粘性动态域——按网格行并行，将不受粘性效应主导的粘性动态域边界单元从粘性动态域中移除；
S14：块间增大动态计算域——按网格块并行遍历尚未参与计算的网格块，逐一判断其内边界的第一层虚网格单元；
S15：判断求解是否达到收敛条件；若是，则跳至步骤S16；若否，则跳至步骤S5；
S16：结果输出。
优选的方式中，所述步骤S1、S2、S5和S15为串行执行。
优选的方式中，所述步骤S2具体过程为：
所述将读入的多块结构化网格合并至最小块数时，依据读入的边界条件，遍历所有内边界；对任一内边界的紧邻两网格块，若其在除内边界法向外的另2个网格方向的单元数完全相同，则将该对网格块合并为一块。
优选的方式中，所述步骤S4中根据壁面边界建立时，按边界并行执行；根据指定流场建立时，按线程并行执行。
优选的方式中，步骤S7具体过程如下：
S71：网格块间进行数据交换，获取相邻块的守恒量信息；
S72：计算并存储第一层虚网格未更新守恒量与新传入的相邻块守恒量之差的模值；
S73：更新虚网格内的存储信息。
优选的方式中，所述步骤S8中无粘项与粘性系数无关，粘性项与粘性系数相关。
优选的方式中，步骤S10具体包括：
S101：判断对流动态域边界单元是否已受到无粘扰动；
S102：若该单元已受无粘扰动，衡量无粘扰动的传播方向，将其会受影响的紧邻单元纳入对流动态域。
优选的方式中，所述步骤S11具体为通过以下4个条件判断某一对流动态域边界单元是否可从对流动态域中移除：
条件1：该单元是否已收敛；
条件2：该单元是否位于最上游；
条件3：该单元是否位于可压缩流动中；
条件4：该单元是否不再受对流动态域中其他单元的影响；
若对流动态域某一边界单元同时满足上述4个条件，则该单元将从对流动态域中移除；若该单元同时存在于粘性动态域中，则也将其从粘性动态域中一并移除，以保证粘性动态域始终包含于对流动态域。
优选的方式中，所述步骤S12具体包括：
S121：判断粘性动态域边界单元是否受粘性效应主导；
S122：若该单元受粘性效应主导，将其所有位于对流动态域中的紧邻单元纳入粘性动态域。
优选的方式中，所述步骤S14具体包括：
S141：根据虚网格的守恒量更新量模值、守恒量和坐标，判断紧邻内边界实单元是否将受到无粘扰动；若是，则将该实单元加入对流动态域；
S142：若紧邻内边界实单元被加入对流动态域，则根据虚网格的守恒量，判断紧邻内边界实单元是否受粘性效应主导；若是，则将该实单元加入粘性动态域。
本发明的一种飞行器气动特性模拟的OpenMP并行扰动域更新方法，可解决现有并行数值模拟方法在模拟飞行器气动特性时存在大量无效计算而现有扰动区域更新方法开展大规模并行时计算效率可能损失的问题。一方面，本发明通过采用对流、粘性两类动态计算域并发展动态计算域更新的并行算法，实现仅对未收敛受扰单元求解、仅在局部区域考虑粘性效应的求解思路，从而有效避免传统并行数值模拟方法中的无效计算。另一方面，本发明利用OpenMP并行共享内存的特点，通过在数据结构和循环结构中引入网格子块的概念，使得网格块内边界数据交换量最小并不受并行线程数的影响。
在基于传统数值模拟方法发展的并行计算中，由于全局更新的计算策略存在着大量无效计算。相比之下，本发明仅在动态计算域中执行数值模拟中最为耗时的残差估计和时间积分两步。其中，步骤S8的残差估计仅在对流动态域内计算残差的无粘项，仅在粘性动态域内计算残差的粘性项；步骤S9的时间积分仅在对流动态域内执行。而本发明所采用的动态计算域仅包含必须参与计算的单元：步骤S10、S14仅将受到无粘扰动的单元加入对流动态域，而步骤S12、S14仅将受粘性效应主导的单元加入粘性动态域；步骤S11会将求解已经收敛并不再影响计算的单元从对流、粘性动态计算域中移除，而步骤S13会将误判为受粘性效应主导的单元从粘性动态域中移除。因此，本发明可始终保持仅对求解未收敛的受扰单元求解，仅在粘性效应主导单元中考虑粘性，从而有效避免传统方法中的无效计算，显著降低并行数值模拟的总计算量，获得相比于传统方法更高的计算效率。
其次，除S1、S2、S5、S15这些耗时极少的步骤外，本发明中的其他步骤均通过OpenMP按线程的粗粒度并行或按网格块、网格行、边界的细粒度并行实现了并行执行，从而通过提高算法并行度确保了其可获得相比于传统方法更高的计算效率。
此外，本发明在步骤S2中通过合并网格块最小化了内边界的数量；在步骤S7中采用了特殊的内边界处理方法，使得内边界无需额外增加数据交换量；利用OpenMP并行共享内存的特点，在流场数据存储的数据结构中引入网格子块的概念，使得网格内边界的处理量不随并行线程数的增大而变化，始终保持在最小量。以上特征都进一步确保了本发明可获得相比于传统方法更高的计算效率。
附图说明
图1为OpenMP并行扰动域更新技术的流程图；
图2a)为读入的多块结构化网格的概念示意图；
图2b)为采用2线程并行模拟时使用网格的概念示意图；
图3为采用12线程并行求解跨声速RAE2822翼型湍流绕流的瞬态流场及动态计算域、网格子块的演化过程图；
图4a)为采用12线程并行求解跨声速RAE2822翼型湍流绕流问题中对流动态域、粘性动态域网格量的变化曲线；
图4b)为与传统方法采用12线程并行求解跨声速RAE2822翼型湍流绕流问题的收敛曲线；
图5为采用16线程并行求解超声速三维楔形体绕流的瞬态流场及动态计算域、网格子块的演化过程图；
图6a)为与传统方法采用16线程并行求解超声速三维楔形体绕流问题的收敛曲线；
图6b)为采用16线程并行求解超声速三维楔形体绕流问题中对流动态域网格量的变化曲线。
具体实施方式
下面详细说明本发明的具体实施，有必要在此指出的是，以下实施只是用于本发明的进一步说明，不能理解为对本发明保护范围的限制，该领域技术熟练人员根据上述本发明内容对本发明做出的一些非本质的改进和调整，仍然属于本发明的保护范围。
图1-6b为本发明具体实现的流程图等示意性展示，下面进行具体的介绍。
首先，本发明提供了一种飞行器气动特性模拟的OpenMP并行扰动域更新方法，其基于OpenMP并行扰动域更新方式，结合附图1所示，具体的包括可依次执行的如下步骤：
S1：数据读入；
分配网格坐标和块信息的存储空间，读入飞行器流场的多块结构化网格、边界条件、计算设置等数据；
S2：静态区域分解；
将读入的多块结构化网格合并至最小块数；根据该网格的单元数将网格分为网格子块并均匀分配给线程；此步为串行执行；
S3：流场初始化；
按线程并行执行，各线程分配存储流动特性的存储空间，根据来流条件或给定流场对飞行器流场网格的所有单元赋初值；
S4：建立动态计算域；
根据壁面边界或指定流场两种方式，建立对流、粘性两类动态计算域；根据壁面边界建立时，按边界并行执行；根据指定流场建立时，按线程并行执行；
S5：动态区域分解；
根据动态计算域单元数，将计算任务平均分配至各线程；此步为串行执行；
S6：分配存储空间；
按线程并行，各线程根据所包含动态计算域的范围分配守恒量更新量、当地时间步长等与更新相关变量的存储空间；
S7：边界条件处理；
按边界并行，为边界虚网格或边界面赋值；物理边界按照其物理定义为虚网格或边界面赋值，内边界需额外存储第一层虚网格的守恒量更新量模值；
S8：残差估计；
将流动控制方程的残差分为无粘项和粘性项：无粘项与粘性系数无关，粘性项与粘性系数相关；按线程并行，在对流动态域中计算残差的无粘项，在粘性动态域中计算残差的粘性项；
S9：时间积分；
按线程并行，在对流动态域中，求解守恒量更新量，更新流场变量；
S10：块内增大对流动态域；
按网格块并行，具体包括如下2个子步：
S101：判断对流动态域边界单元是否已受到无粘扰动；
S102：若该单元已受无粘扰动，衡量无粘扰动的传播方向，将其会受影响的紧邻单元纳入对流动态域；
S11：块内缩小动态计算域；
按网格行并行，通过以下4个条件判断某一对流动态域边界单元是否可从对流动态域中移除：
条件1：该单元是否已收敛；
条件2：该单元是否位于最上游；
条件3：该单元是否位于可压缩流动中；
条件4：该单元是否不再受对流动态域中其他单元的影响；
若对流动态域某一边界单元同时满足上述4个条件，则该单元将从对流动态域中移除；若该单元同时存在于粘性动态域中，则也将其从粘性动态域中一并移除，以保证粘性动态域始终包含于对流动态域；
S12：块内增大粘性动态域；
按网格行并行，具体包括如下2个子步：
S121：判断粘性动态域边界单元是否受粘性效应主导；
S122：若该单元受粘性效应主导，将其所有位于对流动态域中的紧邻单元纳入粘性动态域；
S13：块内缩小粘性动态域；
按网格行并行，将不受粘性效应主导的粘性动态域边界单元从粘性动态域中移除；
S14：块间增大动态计算域；
按网格块并行遍历尚未参与计算的网格块，逐一判断其内边界的第一层虚网格单元，具体包括如下2个子步：
S141：根据虚网格的守恒量更新量模值、守恒量和坐标，判断紧邻内边界实单元是否将受到无粘扰动；若是，则将该实单元加入对流动态域；
S142：若紧邻内边界实单元被加入对流动态域，则根据虚网格的守恒量，判断紧邻内边界实单元是否受粘性效应主导；若是，则将该实单元加入粘性动态域；
S15：判断求解是否达到收敛条件；若是，则跳至步骤S16；若否，则跳至步骤S5；
S16：结果输出。
进一步，步骤S2具体过程如下：
为降低内边界数量以提高计算效率，本发明将读入的多块结构化网格合并至最小块数。依据读入的边界条件，遍历所有内边界；对任一内边界的紧邻两网格块，若其在除内边界法向外的另2个网格方向的单元数完全相同，则将该对网格块合并为一块。
飞行器气动特性模拟中，网格块是分配数组、存储数据的基本单元。为避免网格块数随并行线程数增大而增长，从而降低计算效率，本发明在数据结构中引入网格子块的概念。网格子块是网格块的子区域。网格子块的标号范围用于定义其所属线程所需求解的区域，但其仍与该网格块的其他子块共享同一数组存储数据，从而可避免内边界间的数据传递。
进一步，步骤S7具体过程如下：
与传统方法相比，本发明需额外存储内边界第一层虚网格的守恒量更新量模值。为避免增加网格块间的数据交换量，本发明处理内边界需依次执行如下子步：
S71：网格块间进行数据交换，获取相邻块的守恒量信息；
S72：计算并存储第一层虚网格未更新守恒量（即上一步的相邻块守恒量）与新传入的相邻块守恒量之差的模值；
S73：更新虚网格内的存储信息。
本发明所提飞行器气动特性模拟的OpenMP并行扰动域更新技术可显著提高并行数值模拟的计算效率。
基于上述方式，本发明提供具体的实施例进行说明：
实施例1：
如图1所示，本发明的飞行器气动特性模拟的OpenMP并行扰动域更新技术，具体包括如下步骤：
S1：数据读入：
分配网格坐标和块信息的存储空间，读入飞行器流场的多块结构化网格、边界条件、计算设置等数据。
S2：静态区域分解：
将读入的多块结构化网格合并至最小块数；根据该网格的单元数将网格分为网格子块并均匀分配给线程。图2以二维网格为例示意了本发明对网格的处理方法，图中a、b表示沿网格方向的单元数。图2a)示意了读入的多块结构网格，包括3个网格块，2个内边界。对于网格块1，内边界的法向为J方向，判断可知其与网格块2在I方向的单元数相同，则网格块1可与网格块2合并，形成网格块4。继续判断网格块4，可知其与网格块3在I方向的单元不同，则二者无法合并。因此，图2b)所示的本发明模拟中使用的网格中包含了2个网格块，网格块3和4将分别分配数据存储各自网格块内单元的信息。当采用2个线程并行时，网格块4又进一步划分为2个网格子块，具有4ab个单元的右侧子块将独自分配给1个线程。虽然网格块4左、右2个子块由不同线程计算，但由于数据仍存储于同一数组，因此可省去网格块间内边界虚网格信息间的交换。
S3：流场初始化：
按线程并行执行，各线程分配存储流动特性的存储空间，根据来流条件或给定流场对飞行器流场网格的所有单元赋初值。
S4：建立动态计算域：
a.当根据来流条件初始化时，携带着流场变化信息的扰动从物体表面产生，受粘性主导的单元也为紧邻物体表面。因此，对流动态域可取紧邻壁面的10层单元作为初始单元，粘性动态域可取紧邻壁面的1层单元作为初始单元。
b.当根据给定初始流场初始化时，对流、粘性两类动态域根据给定初始流场的流动特性建立。
对于对流动态域，其应包含给定初始流场中的已受扰单元，这些单元的流动特性应与来流条件不符；故对流动态域单元的流动特性应满足
Figure 624310DEST_PATH_IMAGE001
（1）
式中，W表示守恒量；W ∞表示来流守恒量；ε a,c为对流新增阈值，取10-5。
对于粘性动态域，其应包含对流动态域中的粘性效应主导单元。定义粘性效应衡量参数Ψ，令其表示粘性扰动与无粘扰动的质量流量之比。为降低流动特性、数值格式、网格等因素对粘性动态域的影响，以紧邻壁面单元在求解第1步的粘性效应衡量参数
Figure 10292DEST_PATH_IMAGE002
进行标准化。粘性效应主导单元的粘性效应衡量参数应处于较大量级，故粘性动态域中的单元应满足
Figure 524450DEST_PATH_IMAGE003
（2）
式中，I, J, K表示网格方向；
Figure 887036DEST_PATH_IMAGE004
、
Figure 435829DEST_PATH_IMAGE005
分别表示对流通量和粘性通量Jacobian矩阵沿i方向的谱半径，i=I, J, K；ε a,v表粘性新增阈值，取5×10-3。
S5：动态区域分解：
根据动态计算域单元数，将计算任务平均分配至各线程。
S6：分配存储空间：
按线程并行，各线程采用动态数据结构根据所包含动态计算域的范围分配守恒量更新量、当地时间步长等与更新相关变量的存储空间。
S7：边界条件处理：
按边界并行，为边界虚网格或边界面赋值；物理边界按照其物理定义为虚网格或边界面赋值，内边界需额外存储第一层虚网格的守恒量更新量模值。内边界的处理分为三步：首先，网格块间进行数据交换，获取相邻块的守恒量信息；其次，计算并存储第一层虚网格未更新守恒量（即上一步的相邻块守恒量）与新传入的相邻块守恒量之差的模值；最后，更新虚网格内的存储信息。
S8：残差估计：
本发明将控制方程的残差分为无粘性与粘性项；仅在对流动态域中求解残差的无粘项，在粘性动态域中求解残差的粘性项。流动控制方程可表示为
Figure 727133DEST_PATH_IMAGE006
（3）
其中，W表示守恒量；t表示时间；F c、F v分别表示对流通量和粘性通量；Q T表示湍流模型方程源项；|Ω|、N f、ΔS分别表示网格单元的体积、单元面的面数与面积；式（3）中，右端项统称为控制方程的残差，本发明将残差中的对流通量项称为无粘项，而将残差中的粘性通量项、湍流源项统称为粘性项；无粘性、粘性项分别仅在对流动态域和粘性动态域中求解；所采用的计算方法与传统方法保持一致。
S9：时间积分：
在对流动态域中，采用时间推进格式计算式（3）的左端项；所采用的时间积分方法与传统方法完全一致。
S10：块内增大对流动态域：
令单元(I, J, K)表对流动态域的边界单元，则本发明通过以下2个子步实现对流动态域的增大：
S101：判断单元(I, J, K)是否已受到无粘扰动；
受扰单元的守恒量更新量应处于不可忽略的量级，则已受扰单元满足
Figure 728587DEST_PATH_IMAGE007
（4）
式中，||ΔW||表示求解当前迭代步该对流动态域边界单元守恒量更新量的模值；||ΔW (1)||max表示求解第1步所有单元守恒量更新量的最大值；ε a,c表示对流新增阈值，取10-5；
S102：若单元(I, J, K)已受无粘扰动，衡量无粘扰动的传播方向，将其会受影响的紧邻单元纳入对流动态域。
无粘扰动相对于流动以声速传播。沿扰动传播方向，扰动的传播速度应为正；单元(I, J, K)位于该方向的紧邻单元也将受到扰动的影响。令q为单元(I, J, K)格心指向任一方向的单位向量，则扰动沿q方向传播将满足
Figure 396328DEST_PATH_IMAGE008
（5）
式中，u表流动速度矢量，a为声速；单位向量q取单元(I, J, K)格心指向其格点的单位向量；若单元(I, J, K)格心指向其某一格点的单位向量满足式（5），则将与单元(I,J, K)共享该格点的紧邻单元纳入对流动态域。
S11：块内缩小动态计算域：
为在不影响收敛速率的前提下，尽可能地减少单个迭代步的计算量，对于对流动态域的边界单元(I, J, K)，本发明将通过以下4个条件实现对流动态域的缩小：
条件1：单元(I, J, K)是否已收敛；
已收敛单元的守恒量更新量应处于较小量级，即满足
Figure 534049DEST_PATH_IMAGE009
（6）
式中，||ΔW||表示求解当前迭代步该对流动态域边界单元守恒量更新量的模值；||ΔW (1)||max表示求解第1步所有单元守恒量更新量的最大值；ε d表示删除阈值，取10-7。
条件2：单元(I, J, K)是否位于最上游；
令q表单元(I, J, K)格心指向其紧邻单元格心的单位向量，则该紧邻单元为单元(I, J, K)的下游单元应满足
Figure 261833DEST_PATH_IMAGE010
（7）
式中，u表示流动速度矢量；θ d表示上游单元容差角，超声速流动取10°，亚声速流动取45°；若单元(I, J, K)所有位于对流动态域中的紧邻单元均满足式（7）则认为该单元(I, J, K)位于最上游。
条件3：该单元是否位于可压缩流动中；
根据可压缩流动定义，位于可压缩流动中单元的马赫数应大于0.3。
条件4：单元(I, J, K)是否不再受对流动态域中其他单元的影响；
超声速无粘流动中，由于流动控制方程的数学性质呈双曲型，即流场中任一点流动不受其下游流动的影响，因此满足条件3的超声速无粘单元也自然满足条件4。
亚声速流动和粘性流动中，若考虑紧邻单元对单元(I, J, K)守恒量更新量的影响后，其仍满足收敛条件，则认为该单元不再受其他单元的影响，即满足
Figure 750583DEST_PATH_IMAGE011
（8）
其中，Δ(ΔW)表示对流动态域中的紧邻单元对单元(I, J, K)守恒量更新量的影响，可表示为：
Figure 956437DEST_PATH_IMAGE012
（9）
式中，Δt表示迭代步长；C CFL表示时间推进格式的CFL数；|Ω|表示网格单元的体积；I, J, K表示网格方向；ΔR i 表示对流动态域中的紧邻单元对对流动态域边界单元的残差沿i方向的影响，i=I, J, K；
对于亚声速无粘单元，紧邻单元对对流动态域边界单元的残差沿网格方向的影响ΔR i 可表示为
Figure 214243DEST_PATH_IMAGE013
Figure 112929DEST_PATH_IMAGE014
Figure 321931DEST_PATH_IMAGE015
（10）
式中，ΔW表示守恒量更新量；ΔF c表示对流通量变化量，即当前步与上一步对流通量之差；下标i±1表示对流动态域边界单元沿正、负i方向的紧邻单元；下标i±1/2表示对流动态域边界单元沿正、负i方向的单元面；
Figure 331475DEST_PATH_IMAGE016
表示对流通量Jacobian矩阵沿i方向的谱半径。
对于粘性单元
Figure 443788DEST_PATH_IMAGE017
Figure 513375DEST_PATH_IMAGE018
（11）
若单元(I, J, K)同时满足上述4个条件，删除单元操作会将其从对流动态域中移除。若该单元(I, J, K)同时存在于粘性动态域中，为保证粘性动态域始终包含于对流动态域，则删除单元操作也会将其一并从粘性动态域中移除。
S12：块内增大粘性动态域：
按网格行并行，令单元(I, J, K)表粘性动态域的边界单元，新增单元操作通过以下2个子步实现粘性动态域的增大：
S121：判断单元(I, J, K)是否受粘性效应主导；
若单元(I, J, K)的标准化粘性效应衡量参数具有不可忽略的量级，则该单元须考虑粘性效应。单元(I, J, K)受粘性效应主导应满足
Figure 711138DEST_PATH_IMAGE019
（12）
式中，||ΔW||max表示求解当前迭代步所有单元守恒量更新量的最大值；(||ΔW||max)min表示当前时刻中求解第1至当前迭代步中单步守恒量更新量最大值的最小值。
S122：若单元(I, J, K)受粘性效应主导，将其所有位于对流动态域中的紧邻单元纳入粘性动态域。
控制方程粘性项呈椭圆型。因此，在不增大对流动态域的前提下，将已受粘性扰动单元的所有紧邻单元均加入粘性动态域。
S13：缩小粘性动态域：
逐一判断粘性动态域的边界单元，若该单元及其紧邻单元均不满足式（12），则将该边界单元从粘性动态域中移除。
S14：块间增大动态计算域：
若单元(I, J, K)为内边界的虚网格单元，新增单元操作通过以下2个子步实现块间对流、粘性动态域的增大：
S141：通过式（4）判断单元(I, J, K)是否已受到无粘扰动；若式（4）满足，则通过式（5）判断紧邻实单元是否将受无粘扰动，若满足则将该实单元加入对流动态域。
S142：通过式（12）判断单元(I, J, K)是否受粘性效应主导；若满足，则将已加入对流动态域的实单元也加入粘性动态域。
S15：判断求解是否达到收敛条件；若是，则跳至步骤S16；若否，则跳至步骤S5：
S16：结果输出。
在具体的实验过程中，
（1）对马赫数0.729的RAE2822翼型湍流绕流问题进行模拟。图3示意了本发明OpenMP并行扰动域更新技术的动态计算域演化过程；图中η、η v分别表示对流动态域、粘性动态域与预设计算域的网格量之比，N max代表求解收敛所需的总迭代步数。流场以来流条件初始化，扰动从壁面边界产生；第一步，两类动态计算域根据壁面边界建立。在求解初期，对流动态域从壁面逐渐扩展至整个预设计算域。在求解中后期，对流动态域开始随着求解收敛逐渐从远场向壁面、从上游向下游收缩，最后仅剩近壁区及其下游区域。由图4a)展示的动态计算域单元数变化曲线可知，在仅2/3的迭代中对流动态域与预设计算域的单元数都接近或低于0.5，最大粘性动态域仅为预设计算域的0.488。从图4b)展示的本实施例的动态域网格量变化曲线可知，相比于以相同线程数并行的传统方法，本发明可节省51.5%的计算时间。本发明所获得的加速比也远超不考虑并行开销的理想加速情况。
（2）对马赫数6的三维楔形体绕流问题进行模拟，图5给出了本发明OpenMP并行扰动域更新技术的动态计算域演化过程。流场以来流条件初始化，因此对流动态域依据壁面建立。随着扰动传播，对流动态域逐渐向周围流场增大。但受楔形体周围激波的阻挡，对流动态域不会覆盖整个预设计算域，至多仅包含激波附近及波后区域。在对流动态域增大的同时，上游区域的求解逐渐收敛；因此，对流动态域也从上游至下游逐渐收缩。在求解收敛时，仅有位于最下游激波附近的少量单元参与了最后一步的求解。由图6a)所示的动态计算域单元数变化曲线可知，对流动态域至多仅包含了55.5%的预设计算域单元。由图6b)所示的收敛曲线对比可知，相比于以相同线程数并行的传统方法，本发明可节省65.5%的计算时间。对于三维情况，本发明所获得的加速比仍远超不考虑并行开销的理想加速情况。
尽管为了说明的目的，已描述了本发明的示例性实施方式，但是本领域的技术人员将理解，不脱离所附权利要求中公开的发明的范围和精神的情况下，可以在形式和细节上进行各种修改、添加和替换等的改变，而所有这些改变都应属于本发明所附权利要求的保护范围，并且本发明要求保护的产品各个部门和方法中的各个步骤，可以以任意组合的形式组合在一起。因此，对本发明中所公开的实施方式的描述并非为了限制本发明的范围，而是用于描述本发明。相应地，本发明的范围不受以上实施方式的限制，而是由权利要求或其等同物进行限定。

Claims (10)

1.一种飞行器气动特性模拟的OpenMP并行扰动域更新方法，其特征在于，包括如下步骤：
S1：数据读入——分配网格坐标和块信息的存储空间，读入飞行器流场的多块结构化网格、边界条件、计算设置数据；
S2：静态区域分解——将读入的多块结构化网格合并至最小块数；根据该网格的单元数将网格分为网格子块并均匀分配给线程；
S3：流场初始化——按线程并行执行，各线程分配存储流动特性的存储空间，根据来流条件或给定流场对飞行器流场网格的所有单元赋初值；
S4：建立动态计算域——根据壁面边界或指定流场两种方式，建立对流、粘性两类动态计算域；
S5：动态区域分解——根据动态计算域单元数，将计算任务平均分配至各线程；
S6：分配存储空间——按线程并行，各线程根据所包含动态计算域的范围分配守恒量更新量、当地时间步长等与更新相关变量的存储空间；
S7：边界条件处理——按边界并行，为边界虚网格或边界面赋值；物理边界按照其物理定义为虚网格或边界面赋值，内边界需额外存储第一层虚网格的守恒量更新量模值；
S8：残差估计——将流动控制方程的残差分为无粘项和粘性项：按线程并行，在对流动态域中计算残差的无粘项，在粘性动态域中计算残差的粘性项；
S9：时间积分——按线程并行，在对流动态域中，求解守恒量更新量，更新流场变量；
S10：块内增大对流动态域——按网格块并行，将受扰单元加入对流动态域；
S11：块内缩小动态计算域——按网格行并行，判断某一对流动态域边界单元是否可从对流动态域中移除，并相应调整粘性动态域；
S12：块内增大粘性动态域——按网格行并行，将粘性效应主导单元加入粘性动态域；
S13：缩小粘性动态域——按网格行并行，将不受粘性效应主导的粘性动态域边界单元从粘性动态域中移除；
S14：块间增大动态计算域——按网格块并行遍历尚未参与计算的网格块，逐一判断其内边界的第一层虚网格单元；
S15：判断求解是否达到收敛条件；若是，则跳至步骤S16；若否，则跳至步骤S5；
S16：结果输出。
2.如权利要求1所述的方法，其特征在于：所述步骤S1、S2、S5和S15为串行执行。
3.如权利要求1所述的方法，其特征在于：所述步骤S2具体过程为：所述将读入的多块结构化网格合并至最小块数时，依据读入的边界条件，遍历所有内边界；对任一内边界的紧邻两网格块，若其在除内边界法向外的另2个网格方向的单元数完全相同，则将该对网格块合并为一块。
4.如权利要求3所述的方法，其特征在于：所述步骤S2、S4中根据壁面边界建立时，按边界并行执行；根据指定流场建立时，按线程并行执行。
5.如权利要求4所述的方法，其特征在于：步骤S7具体过程如下：
S71：网格块间进行数据交换，获取相邻块的守恒量信息；
S72：计算并存储第一层虚网格未更新守恒量与新传入的相邻块守恒量之差的模值；
S73：更新虚网格内的存储信息。
6.如权利要求1所述的方法，其特征在于：所述步骤S8中无粘项与粘性系数无关，粘性项与粘性系数相关。
7.如权利要求1所述的方法，其特征在于：步骤S10具体包括：
S101：判断对流动态域边界单元是否已受到无粘扰动；
S102：若该单元已受无粘扰动，衡量无粘扰动的传播方向，将其会受影响的紧邻单元纳入对流动态域。
8.如权利要求1所述的方法，其特征在于：所述步骤S11具体为通过以下4个条件判断某一对流动态域边界单元是否可从对流动态域中移除：
条件1：该单元是否已收敛；
条件2：该单元是否位于最上游；
条件3：该单元是否位于可压缩流动中；
条件4：该单元是否不再受对流动态域中其他单元的影响；
若对流动态域某一边界单元同时满足上述4个条件，则该单元将从对流动态域中移除；若该单元同时存在于粘性动态域中，则也将其从粘性动态域中一并移除，以保证粘性动态域始终包含于对流动态域。
9.如权利要求1所述的方法，其特征在于：所述步骤S12具体包括：
S121：判断粘性动态域边界单元是否受粘性效应主导；
S122：若该单元受粘性效应主导，将其所有位于对流动态域中的紧邻单元纳入粘性动态域。
10.如权利要求1所述的方法，其特征在于：所述步骤S14具体包括：
S141：根据虚网格的守恒量更新量模值、守恒量和坐标，判断紧邻内边界实单元是否将受到无粘扰动；若是，则将该实单元加入对流动态域；
S142：若紧邻内边界实单元被加入对流动态域，则根据虚网格的守恒量，判断紧邻内边界实单元是否受粘性效应主导；若是，则将该实单元加入粘性动态域。
CN202110940149.8A 2021-08-17 2021-08-17 一种飞行器气动特性模拟的OpenMP并行扰动域更新方法 Active CN113392472B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110940149.8A CN113392472B (zh) 2021-08-17 2021-08-17 一种飞行器气动特性模拟的OpenMP并行扰动域更新方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110940149.8A CN113392472B (zh) 2021-08-17 2021-08-17 一种飞行器气动特性模拟的OpenMP并行扰动域更新方法

Publications (2)

Publication Number Publication Date
CN113392472A true CN113392472A (zh) 2021-09-14
CN113392472B CN113392472B (zh) 2021-11-09

Family

ID=77622691

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110940149.8A Active CN113392472B (zh) 2021-08-17 2021-08-17 一种飞行器气动特性模拟的OpenMP并行扰动域更新方法

Country Status (1)

Country Link
CN (1) CN113392472B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113962030A (zh) * 2021-12-20 2022-01-21 北京航空航天大学 飞行器多体分离模拟的重叠网格扰动域更新方法
CN115017604A (zh) * 2022-04-28 2022-09-06 北京航空航天大学 一种可压/不可压混合流动的高效扰动域推进方法
CN116225722A (zh) * 2023-05-08 2023-06-06 中国空气动力研究与发展中心计算空气动力研究所 流场变量的通信方法、装置、终端设备及存储介质
WO2023168772A1 (zh) * 2022-03-09 2023-09-14 北京航空航天大学 一种飞行器动态特性模拟的时间并行扰动域更新方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130198426A1 (en) * 2011-12-22 2013-08-01 Airbus Operations S.L. Heterogeneous parallel systems for accelerating simulations based on discrete grid numerical methods
CN104881510A (zh) * 2015-02-13 2015-09-02 南京航空航天大学 一种直升机旋翼/尾桨气动干扰数值仿真方法
CN108595277A (zh) * 2018-04-08 2018-09-28 西安交通大学 一种基于OpenMP/MPI混合编程的CFD仿真程序的通信优化方法
CN110096838A (zh) * 2019-05-16 2019-08-06 杭州电子科技大学 一种基于n-s方程的直升机流场数值并行隐式求解方法
CN111859530A (zh) * 2020-06-11 2020-10-30 北京航空航天大学 一种飞行器动态气动特性模拟的迭代推进扰动域更新方法
CN112632874A (zh) * 2020-12-31 2021-04-09 杭州电子科技大学 一种直升机流场数值模拟的优化方法及系统

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130198426A1 (en) * 2011-12-22 2013-08-01 Airbus Operations S.L. Heterogeneous parallel systems for accelerating simulations based on discrete grid numerical methods
CN104881510A (zh) * 2015-02-13 2015-09-02 南京航空航天大学 一种直升机旋翼/尾桨气动干扰数值仿真方法
CN108595277A (zh) * 2018-04-08 2018-09-28 西安交通大学 一种基于OpenMP/MPI混合编程的CFD仿真程序的通信优化方法
CN110096838A (zh) * 2019-05-16 2019-08-06 杭州电子科技大学 一种基于n-s方程的直升机流场数值并行隐式求解方法
CN111859530A (zh) * 2020-06-11 2020-10-30 北京航空航天大学 一种飞行器动态气动特性模拟的迭代推进扰动域更新方法
CN112632874A (zh) * 2020-12-31 2021-04-09 杭州电子科技大学 一种直升机流场数值模拟的优化方法及系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
孙学功 等: "高超声速飞行器并行仿真方法研究", 《系统仿真学报》 *
罗磊 等: "CFD 技术发展及其在航空领域中的应用进展", 《航空制造技术》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113962030A (zh) * 2021-12-20 2022-01-21 北京航空航天大学 飞行器多体分离模拟的重叠网格扰动域更新方法
WO2023168772A1 (zh) * 2022-03-09 2023-09-14 北京航空航天大学 一种飞行器动态特性模拟的时间并行扰动域更新方法
CN115017604A (zh) * 2022-04-28 2022-09-06 北京航空航天大学 一种可压/不可压混合流动的高效扰动域推进方法
CN115017604B (zh) * 2022-04-28 2024-06-07 北京航空航天大学 一种可压/不可压混合流动的高效扰动域推进方法
CN116225722A (zh) * 2023-05-08 2023-06-06 中国空气动力研究与发展中心计算空气动力研究所 流场变量的通信方法、装置、终端设备及存储介质

Also Published As

Publication number Publication date
CN113392472B (zh) 2021-11-09

Similar Documents

Publication Publication Date Title
CN113392472B (zh) 一种飞行器气动特性模拟的OpenMP并行扰动域更新方法
CN111859530B (zh) 一种飞行器动态气动特性模拟的迭代推进扰动域更新方法
CN113505443B (zh) 一种任意外形的三维绕流问题自适应笛卡尔网格生成方法
CN114329799B (zh) 一种飞行器动态特性模拟的时间并行扰动域更新方法
CN113850008B (zh) 飞行器气动特性预测的自适应网格扰动域更新加速方法
CN111651831B (zh) 飞行器定常粘性可压缩绕流的分区扰动域更新计算方法
CN111859529B (zh) 一种飞行器绕流数值模拟的多重网格扰动域更新加速方法
CN110516316B (zh) 一种间断伽辽金法求解欧拉方程的gpu加速方法
Fietz et al. Optimized hybrid parallel lattice Boltzmann fluid flow simulations on complex geometries
CN113609598B (zh) 飞行器气动特性模拟的rans/les扰动域更新方法
CN113609597B (zh) 超声速飞行器绕流的时间-空间混合推进扰动域更新方法
Shershnev et al. HyCFS, a high-resolution shock capturing code for numerical simulation on hybrid computational clusters
Williamschen et al. Parallel anisotropic block-based adaptive mesh refinement algorithm for three-dimensional flows
Menshov et al. Efficient parallel shock-capturing method for aerodynamics simulations on body-unfitted cartesian grids
Deng et al. A hybrid aerodynamic optimization algorithm based on differential evolution and RBF response surface
Jung et al. An implementation of self-organizing maps for airfoil design exploration via multi-objective optimization technique
Antepara et al. Parallel adaptive mesh refinement simulation of the flow around a square cylinder at Re= 22000
CN114218878B (zh) 一种飞行器机动过程模拟的动网格扰动域更新方法
Liu et al. Massively parallel CFD simulation software: CCFD development and optimization based on Sunway TaihuLight
Shchur et al. Evolution of time horizons in parallel and grid simulations
Biazewicz et al. Problems Related to Parallelization of CFD Algorithms on GPU, Multi‐GPU and Hybrid Architectures
Liu et al. Implementation of the HLL-GRP solver for multidimensional ideal MHD simulations based on finite volume method
Schwing et al. Parallelization of unsteady adaptive mesh refinement for unstructured Navier-Stokes solvers
US20220222399A1 (en) Data structure, numerical calculation method, and numerical calculation program
Liu et al. Feature-aligned poly-square mapping of large-scale 2D geometries for semi-structured quad mesh generation

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