一种飞行器绕流数值模拟的多重网格扰动域更新加速方法技术

本发明专利技术公开了一种飞行器绕流数值模拟的多重网格扰动域更新加速方法，针对现有多重网格方法在模拟飞行器绕流时存在大量无效计算的问题，通过在细网格和粗网格上采用对流、粘性两类动态计算域，使得流动控制方程在每个迭代步中均在仅包含未收敛受扰单元的对流动态计算域中求解，而粘性效应也仅在只包含粘性效应主导单元的粘性动态计算域内考虑，从而实现了仅对未收敛受扰单元求解、仅在局部区域考虑粘性效应的求解思路。通过有效减少单个迭代步的计算量和达到收敛态所需的迭代步数，本发明专利技术可显著提高多重网格方法模拟飞行器绕流的计算效率。

【技术实现步骤摘要】
一种飞行器绕流数值模拟的多重网格扰动域更新加速方法
本专利技术涉及计算流体力学
，尤其涉及一种飞行器绕流数值模拟的多重网格扰动域更新加速方法。
技术介绍
随着计算机技术的飞速发展，计算流体力学数值模拟已成为当代飞行器气动设计的重要技术支撑之一。尽管数值模拟能够精准预测飞行器的静、动态气动特性，但对于飞行器气动设计而言，仍存在数值模拟计算效率难以满足日益增长的实际应用需求的问题。因此，提出适用于飞行器绕流数值模拟的高效方法对提高飞行器设计迭代效率、缩短研制周期具有重要意义和工程应用价值。多重网格方法是飞行器绕流数值模拟的主流加速技术之一。该方法既可用于飞行器静态气动特性预测中定常绕流流场的迭代求解，也适用于飞行器动态气动特性预测中非定常绕流流场迭代推进中的内迭代求解。多重网格方法是一种在粗、细网格间循环求解流动控制方程的加速收敛技术，其加速原理有两方面：其一，粗网格能以更大的迭代步长推进求解，从而可以带动细网格更快达到收敛态；其二，粗、细网格间的循环求解有助于衰减阻碍细网格收敛的高频误差，从而可以使细网格的整体误差能更快地衰减。目前，多重网格方法针对飞行器绕流数值模拟时仍存在一个弊端，即每个迭代步均必须对飞行器流场网格预设计算域中的所有网格单元更新求解。这种全局更新的求解方式忽略了求解过程中飞行器绕流流场演化的特点，从而会造成不同程度的无效计算。对于亚声速飞行器，流动的数学性质为椭圆型，因此，为避免数值模拟的有限边界所导致的非物理响应，飞行器流场网格的预设计算域需要取数十倍于飞行器参考长度的范围进行模拟。然而，随着飞行器绕流流场逐渐趋近于定常态，不满足流动控制方程的区域也相应收缩，飞行器流场网格的预设计算域也将趋于冗余。对于超声速飞行器，流动的数学性质为双曲型，即流场中任一点仅受其上游依赖域内流动的影响。因此，在迭代过程中，位于飞行器绕流流场下游的区域不会先于其上游流场收敛，已达到定常态的上游流动也不会受到未收敛下游流动更新的影响。由此可知，在整个迭代求解过程中，始终对飞行器流场网格预设计算域内的所有网格单元全局更新求解将产生大量无效计算，从而显著降低飞行器绕流数值模拟的计算效率。因此，若能消除现有方法中的无效计算，必将显著提升飞行器绕流数值模拟的计算效率，从而有效解决飞行器气动设计中数值模拟计算效率难以满足日益增长的实际应用需求的问题。
技术实现思路
有鉴于此，本专利技术提供了一种飞行器绕流数值模拟的多重网格扰动域更新加速方法，用以解决现有多重网格方法在模拟飞行器绕流时存在大量无效计算的问题。本专利技术提供的一种飞行器绕流数值模拟的多重网格扰动域更新加速方法，其特征在于，包括如下步骤：S1：读入数据，包括飞行器流场的细网格、预设计算域、边界条件和计算设置；S2：通过稀疏所述细网格，生成指定层数的粗网格；S3：根据来流条件或根据给定流场，对所述粗网格和所述细网格上预设计算域中的流场进行初始化；S4：根据流场初始化方式，建立对流动态计算域和粘性动态计算域；S5：将所有层网格求解所需存储的数据分为两类；第一类数据是单元的固有信息，采用静态数据结构存储所述预设计算域中所有网格单元的信息，包括网格坐标和流场变量；第二类数据是与求解更新相关的信息，采用动态数据结构仅存储所述对流动态计算域内网格单元的信息，包括守恒量更新量和当地迭代步长；S6：根据边界条件类型，为当前层网格的边界虚网格的守恒量赋值；S7：将流动控制方程的残差项分为无粘项与粘性项两类；其中，所述无粘项包括对流通量和粗网格的强迫函数，所述无粘项在对流动态计算域中计算；所述粘性项包括粘性通量和湍流模型方程源项，所述粘性项在粘性动态计算域中计算；S8：在所述对流动态计算域中，求解守恒量更新量，更新流场变量；S9：判断当前网格层是否为细网格；若是，则执行步骤S10；若否，则执行步骤S21；S10：判断细网格是否达到收敛条件；若否，则执行步骤S11；若是，则执行步骤S23；S11：对所述对流动态计算域的所有边界单元逐个判断是否已受到无粘扰动；若是，则执行步骤S12；若否，则返回步骤S11，对所述对流动态计算域的下一边界单元进行判断；遍历所述对流动态计算域的所有边界单元后，执行步骤S13；S12：衡量无粘扰动的传播方向，并将该边界单元的紧邻单元中位于传播方向上的单元纳入所述对流动态计算域，返回步骤S11，对所述对流动态计算域的下一边界单元进行判断；S13：对所述对流动态计算域的所有边界单元逐个判断是否同时满足已收敛、位于可压缩流动中、位于最上游以及不再受所述对流动态计算域中其他单元的影响四个条件；若是，则执行步骤S14；若否，则返回步骤S13，对所述对流动态计算域的下一边界单元进行判断；遍历所述对流动态计算域的所有边界单元后，执行步骤S16；S14：将该边界单元从所述对流动态计算域中移除，并判断该边界单元是否存在于所述粘性动态计算域中；若是，则执行步骤S15；若否，则返回步骤S13，对所述对流动态计算域的下一边界单元进行判断；S15：将该边界单元从所述粘性动态计算域中移除，返回步骤S13，对所述对流动态计算域的下一边界单元进行判断；S16：对所述粘性动态计算域的所有边界单元逐个判断是否受粘性效应主导；若是，则执行步骤S17；若否，则返回步骤S16，对所述粘性动态计算域的下一边界单元进行判断；遍历所述粘性动态计算域的所有边界单元后，执行步骤S18；S17：将该边界单元的紧邻单元中位于所述对流动态计算域中的单元纳入所述粘性动态计算域，返回步骤S16，对所述粘性动态计算域的下一边界单元进行判断；S18：根据细网格的对流动态计算域和粘性动态计算域，更新所有层粗网格的对流动态计算域和粘性动态计算域；S19：根据所述对流动态计算域，重新分配所有粗网格和细网格中第二类数据的存储空间；S20：将守恒量插值至更粗网格，返回步骤S6，进入更粗网格的计算；S21：判断是否需要计算更粗网格；若是，则返回步骤S20；若否，则执行步骤S22；S22：将守恒量修正量插值回更细网格，并判断更细网格是否为细网格；若是，则返回步骤S6，进入下一迭代步的计算；若否，则返回步骤S6，进入更细网格的计算；S23：输出结果。在一种可能的实现方式中，在本专利技术提供的上述多重网格扰动域更新加速方法中，步骤S3，根据来流条件或根据给定流场，对所述粗网格和所述细网格上预设计算域中的流场进行初始化，具体包括：根据来流条件初始化，将所述粗网格和所述细网格上预设计算域中所有网格单元的守恒量赋为来流值；根据给定流场初始化，将所述细网格上预设计算域中所有网格单元的守恒量赋为给定流场值，所述粗网格上预设计算域中所有网格单元的守恒量根据细网格上所有网格单元的守恒量插值获得。在一种可能的实现方式中，在本专利技术提供的上述多重网格扰动域更新加速方法中，步骤S4，根据流场初始化方式，建立对流动态计算域和粘性动态计算域，具体包括：本文档来自技高网...

【技术保护点】
1.一种飞行器绕流数值模拟的多重网格扰动域更新加速方法，其特征在于，包括如下步骤：/nS1：读入数据，包括飞行器流场的细网格、预设计算域、边界条件和计算设置；/nS2：通过稀疏所述细网格，生成指定层数的粗网格；/nS3：根据来流条件或根据给定流场，对所述粗网格和所述细网格上预设计算域中的流场进行初始化；/nS4：根据流场初始化方式，建立对流动态计算域和粘性动态计算域；/nS5：将所有层网格求解所需存储的数据分为两类；第一类数据是单元的固有信息，采用静态数据结构存储所述预设计算域中所有网格单元的信息，包括网格坐标和流场变量；第二类数据是与求解更新相关的信息，采用动态数据结构仅存储所述对流动态计算域内网格单元的信息，包括守恒量更新量和当地迭代步长；/nS6：根据边界条件类型，为当前层网格的边界虚网格的守恒量赋值；/nS7：将流动控制方程的残差项分为无粘项与粘性项两类；其中，所述无粘项包括对流通量和粗网格的强迫函数，所述无粘项在对流动态计算域中计算；所述粘性项包括粘性通量和湍流模型方程源项，所述粘性项在粘性动态计算域中计算；/nS8：在所述对流动态计算域中，求解守恒量更新量，更新流场变量；/nS9：判断当前网格层是否为细网格；若是，则执行步骤S10；若否，则执行步骤S21；/nS10：判断细网格是否达到收敛条件；若否，则执行步骤S11；若是，则执行步骤S23；/nS11：对所述对流动态计算域的所有边界单元逐个判断是否已受到无粘扰动；若是，则执行步骤S12；若否，则返回步骤S11，对所述对流动态计算域的下一边界单元进行判断；遍历所述对流动态计算域的所有边界单元后，执行步骤S13；/nS12：衡量无粘扰动的传播方向，并将该边界单元的紧邻单元中位于传播方向上的单元纳入所述对流动态计算域，返回步骤S11，对所述对流动态计算域的下一边界单元进行判断；/nS13：对所述对流动态计算域的所有边界单元逐个判断是否同时满足已收敛、位于可压缩流动中、位于最上游以及不再受所述对流动态计算域中其他单元的影响四个条件；若是，则执行步骤S14；若否，则返回步骤S13，对所述对流动态计算域的下一边界单元进行判断；遍历所述对流动态计算域的所有边界单元后，执行步骤S16；/nS14：将该边界单元从所述对流动态计算域中移除，并判断该边界单元是否存在于所述粘性动态计算域中；若是，则执行步骤S15；若否，则返回步骤S13，对所述对流动态计算域的下一边界单元进行判断；/nS15：将该边界单元从所述粘性动态计算域中移除，返回步骤S13，对所述对流动态计算域的下一边界单元进行判断；/nS16：对所述粘性动态计算域的所有边界单元逐个判断是否受粘性效应主导；若是，则执行步骤S17；若否，则返回步骤S16，对所述粘性动态计算域的下一边界单元进行判断；遍历所述粘性动态计算域的所有边界单元后，执行步骤S18；/nS17：将该边界单元的紧邻单元中位于所述对流动态计算域中的单元纳入所述粘性动态计算域，返回步骤S16，对所述粘性动态计算域的下一边界单元进行判断；/nS18：根据细网格的对流动态计算域和粘性动态计算域，更新所有层粗网格的对流动态计算域和粘性动态计算域；/nS19：根据所述对流动态计算域，重新分配所有粗网格和细网格中第二类数据的存储空间；/nS20：将守恒量插值至更粗网格，返回步骤S6，进入更粗网格的计算；/nS21：判断是否需要计算更粗网格；若是，则返回步骤S20；若否，则执行步骤S22；/nS22：将守恒量修正量插值回更细网格，并判断更细网格是否为细网格；若是，则返回步骤S6，进入下一迭代步的计算；若否，则返回步骤S6，进入更细网格的计算；/nS23：输出结果。/n...

【技术特征摘要】
1.一种飞行器绕流数值模拟的多重网格扰动域更新加速方法，其特征在于，包括如下步骤：
S1：读入数据，包括飞行器流场的细网格、预设计算域、边界条件和计算设置；
S2：通过稀疏所述细网格，生成指定层数的粗网格；
S3：根据来流条件或根据给定流场，对所述粗网格和所述细网格上预设计算域中的流场进行初始化；
S4：根据流场初始化方式，建立对流动态计算域和粘性动态计算域；
S5：将所有层网格求解所需存储的数据分为两类；第一类数据是单元的固有信息，采用静态数据结构存储所述预设计算域中所有网格单元的信息，包括网格坐标和流场变量；第二类数据是与求解更新相关的信息，采用动态数据结构仅存储所述对流动态计算域内网格单元的信息，包括守恒量更新量和当地迭代步长；
S6：根据边界条件类型，为当前层网格的边界虚网格的守恒量赋值；
S7：将流动控制方程的残差项分为无粘项与粘性项两类；其中，所述无粘项包括对流通量和粗网格的强迫函数，所述无粘项在对流动态计算域中计算；所述粘性项包括粘性通量和湍流模型方程源项，所述粘性项在粘性动态计算域中计算；
S8：在所述对流动态计算域中，求解守恒量更新量，更新流场变量；
S9：判断当前网格层是否为细网格；若是，则执行步骤S10；若否，则执行步骤S21；
S10：判断细网格是否达到收敛条件；若否，则执行步骤S11；若是，则执行步骤S23；
S11：对所述对流动态计算域的所有边界单元逐个判断是否已受到无粘扰动；若是，则执行步骤S12；若否，则返回步骤S11，对所述对流动态计算域的下一边界单元进行判断；遍历所述对流动态计算域的所有边界单元后，执行步骤S13；
S12：衡量无粘扰动的传播方向，并将该边界单元的紧邻单元中位于传播方向上的单元纳入所述对流动态计算域，返回步骤S11，对所述对流动态计算域的下一边界单元进行判断；
S13：对所述对流动态计算域的所有边界单元逐个判断是否同时满足已收敛、位于可压缩流动中、位于最上游以及不再受所述对流动态计算域中其他单元的影响四个条件；若是，则执行步骤S14；若否，则返回步骤S13，对所述对流动态计算域的下一边界单元进行判断；遍历所述对流动态计算域的所有边界单元后，执行步骤S16；
S14：将该边界单元从所述对流动态计算域中移除，并判断该边界单元是否存在于所述粘性动态计算域中；若是，则执行步骤S15；若否，则返回步骤S13，对所述对流动态计算域的下一边界单元进行判断；
S15：将该边界单元从所述粘性动态计算域中移除，返回步骤S13，对所述对流动态计算域的下一边界单元进行判断；
S16：对所述粘性动态计算域的所有边界单元逐个判断是否受粘性效应主导；若是，则执行步骤S17；若否，则返回步骤S16，对所述粘性动态计算域的下一边界单元进行判断；遍历所述粘性动态计算域的所有边界单元后，执行步骤S18；
S17：将该边界单元的紧邻单元中位于所述对流动态计算域中的单元纳入所述粘性动态计算域，返回步骤S16，对所述粘性动态计算域的下一边界单元进行判断；
S18：根据细网格的对流动态计算域和粘性动态计算域，更新所有层粗网格的对流动态计算域和粘性动态计算域；
S19：根据所述对流动态计算域，重新分配所有粗网格和细网格中第二类数据的存储空间；
S20：将守恒量插值至更粗网格，返回步骤S6，进入更粗网格的计算；
S21：判断是否需要计算更粗网格；若是，则返回步骤S20；若否，则执行步骤S22；
S22：将守恒量修正量插值回更细网格，并判断更细网格是否为细网格；若是，则返回步骤S6，进入下一迭代步的计算；若否，则返回步骤S6，进入更细网格的计算；
S23：输出结果。


2.如权利要求1所述的多重网格扰动域更新加速方法，其特征在于，步骤S3，根据来流条件或根据给定流场，对所述粗网格和所述细网格上预设计算域中的流场进行初始化，具体包括：
根据来流条件初始化，将所述粗网格和所述细网格上预设计算域中所有网格单元的守恒量赋为来流值；
根据给定流场初始化，将所述细网格上预设计算域中所有网格单元的守恒量赋为给定流场值，所述粗网格上预设计算域中所有网格单元的守恒量根据细网格上所有网格单元的守恒量插值获得。


3.如权利要求1所述的多重网格扰动域更新加速方法，其特征在于，步骤S4，根据流场初始化方式，建立对流动态计算域和粘性动态计算域，具体包括：
根据来流条件初始化，对于对流动态计算域，最粗层网格取紧邻壁面的1层单元作为初始单元，细网格取壁面的L层相邻单元作为初始单元，L满足：
L＝2l(1)
其中，l表示粗网格层数，l为正整数，由步骤S1读入的计算设置给定；较粗层网格的初始单元根据细网格的初始单元生成；对于粘性动态计算域，所有网格层均取紧邻壁面的1层单元作为初始单元；
根据给定流场初始化，对于对流动态计算域，细网格的初始单元为给定流场中流动特性与来流条件不符的已受扰单元，满足：
||W-W∞||/||ΔW(1)||max＞εa,c(2)
其中，W表示守恒量；W∞表示来流条件的守恒量；||ΔW(1)||max表示对流动态计算域所有单元的守恒量更新量在第1迭代步的最大值；εa,c表示对流新增阈值，取10-4≤εa,c≤10-6；对于粘性动态计算域，细网格的初始单元为细网格对流动态计算域中的粘性效应主导单元，满足：

...

【专利技术属性】
技术研发人员：蒋崇文，胡姝瑶，高振勋，李椿萱，
申请(专利权)人：北京航空航天大学，
类型：发明
国别省市：北京;11