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

本发明涉及计算流体力学技术领域，尤其涉及一种飞行器绕流数值模拟的多重网格扰动域更新加速方法。

背景技术：

随着计算机技术的飞速发展，计算流体力学数值模拟已成为当代飞行器气动设计的重要技术支撑之一。尽管数值模拟能够精准预测飞行器的静、动态气动特性，但对于飞行器气动设计而言，仍存在数值模拟计算效率难以满足日益增长的实际应用需求的问题。因此，提出适用于飞行器绕流数值模拟的高效方法对提高飞行器设计迭代效率、缩短研制周期具有重要意义和工程应用价值。

多重网格方法是飞行器绕流数值模拟的主流加速技术之一。该方法既可用于飞行器静态气动特性预测中定常绕流流场的迭代求解，也适用于飞行器动态气动特性预测中非定常绕流流场迭代推进中的内迭代求解。多重网格方法是一种在粗、细网格间循环求解流动控制方程的加速收敛技术，其加速原理有两方面：其一，粗网格能以更大的迭代步长推进求解，从而可以带动细网格更快达到收敛态；其二，粗、细网格间的循环求解有助于衰减阻碍细网格收敛的高频误差，从而可以使细网格的整体误差能更快地衰减。

目前，多重网格方法针对飞行器绕流数值模拟时仍存在一个弊端，即每个迭代步均必须对飞行器流场网格预设计算域中的所有网格单元更新求解。这种全局更新的求解方式忽略了求解过程中飞行器绕流流场演化的特点，从而会造成不同程度的无效计算。对于亚声速飞行器，流动的数学性质为椭圆型，因此，为避免数值模拟的有限边界所导致的非物理响应，飞行器流场网格的预设计算域需要取数十倍于飞行器参考长度的范围进行模拟。然而，随着飞行器绕流流场逐渐趋近于定常态，不满足流动控制方程的区域也相应收缩，飞行器流场网格的预设计算域也将趋于冗余。对于超声速飞行器，流动的数学性质为双曲型，即流场中任一点仅受其上游依赖域内流动的影响。因此，在迭代过程中，位于飞行器绕流流场下游的区域不会先于其上游流场收敛，已达到定常态的上游流动也不会受到未收敛下游流动更新的影响。由此可知，在整个迭代求解过程中，始终对飞行器流场网格预设计算域内的所有网格单元全局更新求解将产生大量无效计算，从而显著降低飞行器绕流数值模拟的计算效率。因此，若能消除现有方法中的无效计算，必将显著提升飞行器绕流数值模拟的计算效率，从而有效解决飞行器气动设计中数值模拟计算效率难以满足日益增长的实际应用需求的问题。

技术实现要素：

有鉴于此，本发明提供了一种飞行器绕流数值模拟的多重网格扰动域更新加速方法，用以解决现有多重网格方法在模拟飞行器绕流时存在大量无效计算的问题。

本发明提供的一种飞行器绕流数值模拟的多重网格扰动域更新加速方法，其特征在于，包括如下步骤：

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层单元作为初始单元，细网格取壁面的l层相邻单元作为初始单元，l满足：

l＝2^l(1)

其中，l表示粗网格层数，l为正整数，由步骤s1读入的计算设置给定；较粗层网格的初始单元根据细网格的初始单元生成；对于粘性动态计算域，所有网格层均取紧邻壁面的1层单元作为初始单元；

根据给定流场初始化，对于对流动态计算域，细网格的初始单元为给定流场中流动特性与来流条件不符的已受扰单元，满足：

||w-w∞||/||δw⁽¹⁾||max＞εa,c(2)

其中，w表示守恒量；w∞表示来流条件的守恒量；||δw⁽¹⁾||max表示对流动

态计算域所有单元的守恒量更新量在第1迭代步的最大值；εa,c表示对流新增阈值，取10^-4≤εa,c≤10^-6；对于粘性动态计算域，细网格的初始单元为细网格对流动态计算域中的粘性效应主导单元，满足：

其中，ψ表示粘性效应衡量参数，为粘性扰动与无粘扰动的质量流量之比；表示第1迭代步的粘性效应衡量参数；i,j,k分别表示网格方向；表示对流通量jacobian矩阵沿i方向的谱半径，表示粘性通量jacobian矩阵沿i方向的谱半径；εa,v表示粘性新增阈值，取10^-3≤εa,v≤10^-2；所有粗网格层的对流动态计算域依据细网格的对流动态计算域生成，所有粗网格层的粘性动态计算域依据细网格的粘性动态计算域生成。

在一种可能的实现方式中，在本发明提供的上述多重网格扰动域更新加速方法中，步骤s7，将流动控制方程的残差项分为无粘项与粘性项两类；其中，所述无粘项包括对流通量和粗网格的强迫函数，所述无粘项在对流动态计算域中计算；所述粘性项包括粘性通量和湍流模型方程源项，所述粘性项在粘性动态计算域中计算，具体包括：

流动控制方程表示为：

其中，t表示时间，fc表示对流通量，fv表示粘性通量，qf表示粗网格的强迫函数，qt表示湍流模型方程源项，|ω|表示网格单元的体积，nf表示单元面的面数，δs表示单元面的面积；式(5)的右端项统称为流动控制方程的残差项，为无粘项，为粘性项；在对流动态计算域中求解残差项的无粘项，在粘性动态计算域中求解残差项的粘性项。

在一种可能的实现方式中，在本发明提供的上述多重网格扰动域更新加速方法中，步骤s11，对所述对流动态计算域的所有边界单元逐个判断是否已受到无粘扰动，具体包括：

对所述对流动态计算域的所有边界单元逐个判断是否满足：

||δw||/||δw⁽¹⁾||max＞εa,c(6)

其中，||δw||表示对流动态计算域的边界单元在当前迭代步中守恒量更新量的模值；

步骤s12，衡量无粘扰动的传播方向，并将该边界单元的紧邻单元中位于传播方向上的单元纳入所述对流动态计算域，具体包括：

无粘扰动相对于流动以声速传播，沿无粘扰动的传播方向，无粘扰动的传播速度为正；令q表示单位方向向量，则无粘扰动沿q方向传播表示为：

u·q+a＞0(7)

其中，u表示流动速度矢量，a表示声速；q取对流动态计算域边界单元的格心指向边界单元某一格点的单位向量，若q满足式(7)，则将与对流动态计算域边界单元共享该格点的紧邻单元纳入对流动态计算域。

在一种可能的实现方式中，在本发明提供的上述多重网格扰动域更新加速方法中，步骤s13，对所述对流动态计算域的所有边界单元逐个判断是否同时满足已收敛、位于可压缩流动中、位于最上游以及不再受所述对流动态计算域中其他单元的影响四个条件；若是，执行步骤s14；若否，则返回步骤s13，对所述对流动态计算域的下一边界单元进行判断；具体包括如下步骤：

s131：对所述对流动态计算域的所有边界单元逐个判断是否满足：

||δw||/||δw⁽¹⁾||max＜εd(8)

其中，εd表示删除阈值，取εd≤10^-7；若是，则执行步骤s132；若否，则返回步骤s131，对所述对流动态计算域的下一边界单元进行判断；

s132：判断该边界单元的马赫数是否大于0.3；若是，则执行步骤s133；若否，则返回步骤s131，对所述对流动态计算域的下一边界单元进行判断；

s133：判断该边界单元是否满足：

其中，q取对流动态计算域边界单元的格心指向边界单元紧邻单元格心的单位向量，θd表示上游单元容差角，取5°≤θd≤10°；若是，则执行步骤s134；若否，则返回步骤s131，对所述对流动态计算域的下一边界单元进行判断；

s134：对于超声速无粘流动，无需再判断，直接执行步骤s14；

对于亚声速流动和粘性流动，判断该边界单元是否满足：

||δw||+||δ(δw)||＜εd(10)

其中，||δ(δw)||表示边界单元守恒量更新量受边界单元对流动态计算域中紧邻单元影响的模值；δt表示迭代步长，ccfl表示时间推进格式的cfl数，δrⁱ表示对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响；

对于亚声速无粘单元，对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响δrⁱ表示为：

其中，δfc表示对流通量的变化量；下标i+1表示对流动态计算域的边界单元沿正i方向的紧邻单元，下标i-1表示对流动态计算域的边界单元沿负i方向的紧邻单元，δsi+1/2表示对流动态计算域的边界单元沿正i方向单元面的面积；δsi-1/2表示对流动态计算域的边界单元沿负i方向单元面的面积；

对于粘性单元，对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响δrⁱ表示为：

若是，则执行步骤s14；若否，则返回步骤s131，对所述对流动态计算域的下一边界单元进行判断。

在一种可能的实现方式中，在本发明提供的上述多重网格扰动域更新加速方法中，步骤s16，对所述粘性动态计算域的所有边界单元逐个判断是否受粘性效应主导，具体包括：

对所述粘性动态计算域的所有边界单元逐个判断是否满足：

其中，φ表示粘性新增阈值的缩比因子，||δw||max表示对流动态计算域所有单元的守恒量更新量在当前迭代步的最大值，(||δw||max)min表示第1迭代步至当前迭代步中单步守恒量更新量最大值的最小值。

本发明提供的上述飞行器绕流数值模拟的多重网格扰动域更新加速方法，残差估计仅在各层网格的对流动态计算域内计算残差的无粘项，仅在各层网格的粘性动态计算域内计算残差的粘性项，时间积分也仅在各层网格的对流动态计算域内执行，因此，本发明通过在各层网格采用对流、粘性两类动态计算域，可同时降低多重网格方法模拟飞行器绕流问题时各层网格的计算量，显著提升多重网格方法针对飞行器绕流数值模拟的计算效率，从而有效解决数值模拟计算效率难以满足飞行器气动设计日益增长的实际应用需求的问题。将所需存储的数据分为所有单元的固有信息和与求解更新相关的信息两类；其中，对于第一类数据，以静态数据结构存储预设计算域中所有单元的信息；对于第二类数据，以动态数据结构仅存储各网格层对流动态计算域中网格单元的信息，并在每个迭代步及时根据各网格层的对流动态计算域调整第二类数据的存储空间，因此，相比于现有多重网格方法，本发明可有效降低飞行器绕流数值模拟的存储空间需求。针对飞行器流场的细网格，根据扰动传播特性仅将受扰动单元加入对流动态计算域中；通过衡量粘性效应的影响，仅将受粘性效应主导的网格单元加入细网格的粘性动态计算域中；利用求解收敛特征，在不影响收敛速率的前提下，将已收敛单元从对流、粘性两类动态计算域中移除；根据细网格的动态计算域，更新飞行器流场的所有粗网格的对流、粘性两类动态计算域；因此，本发明可始终保持在飞行器流场的所有网格层仅对求解未收敛的受扰单元求解，仅在粘性效应主导单元中考虑粘性，从而可有效避免现有多重网格方法在飞行器绕流数值模拟中的无效计算，解决数值模拟计算效率难以满足飞行器气动设计日益增长的实际应用需求的问题。

附图说明

图1为本发明实施例1提供的多重网格扰动域更新加速方法的流程图；

图2为本发明实施例1中多重网格根据细网格生成粗网格的示意图；

图3为本发明实施例1中多重网格循环示意图；

图4为利用本发明实施例1提供的多重网格扰动域更新加速方法求解跨声速naca0012翼型无粘绕流的流场与对流动态计算域的演化过程图；

图5a为利用本发明实施例1提供的多重网格扰动域更新加速方法求解跨声速naca0012翼型无粘绕流的对流动态计算域网格量变化曲线；

图5b为利用本发明实施例1提供的多重网格扰动域更新加速方法求解跨声速naca0012翼型无粘绕流的收敛曲线；

图6为利用本发明实施例1提供的多重网格扰动域更新加速方法求解跨声速rae2822翼型湍流绕流的流场与对流、粘性两类动态计算域的演化过程图；

图7a为利用本发明实施例1提供的多重网格扰动域更新加速方法求解跨声速rae2822翼型湍流绕流的对流、粘性两类动态计算域网格量变化曲线；

图7b为利用本发明实施例1提供的多重网格扰动域更新加速方法求解跨声速rae2822翼型湍流绕流的收敛曲线。

具体实施方式

下面将结合本发明实施方式中的附图，对本发明实施方式中的技术方案进行清楚、完整的描述，显然，所描述的实施方式仅仅是作为例示，并非用于限制本发明。

采用多重网格方法模拟飞行器绕流的求解过程具有如下特征：(1)流场的扰动始于控制方程无法满足之处；(2)携带着间断信息的扰动会随时间以有限速度逐渐传向周围流场，尚未受到扰动的区域仍保持初始状态；(3)可压缩流动中，上游流场的收敛不会晚于下游流场；(4)无论是否分离，粘性效应仅在有限区域起主导作用。为充分利用上述4点特征，避免现有多重网格方法中的无效计算，本发明提供了一种飞行器绕流数值模拟的多重网格扰动域更新加速方法，通过在所有网格层仅对求解未收敛的收敛单元求解，仅在粘性效应主导单元中考虑粘性，有效减少单个迭代步的计算量和达到收敛态所需的迭代步数，从而显著提高多重网格方法模拟飞行器绕流的计算效率。

本发明提供的一种飞行器绕流数值模拟的多重网格扰动域更新加速方法，包括如下步骤：

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：输出结果。

下面通过一个具体的实施例对本发明提供的上述多重网格扰动域更新加速方法的具体实施进行详细说明。

实施例1：如图1所示，本发明提供的上述多重网格扰动域更新加速方法具体包括如下步骤：

第一步：数据读入。

读入数据包括飞行器流场的细网格、预设计算域、边界条件和计算设置。

第二步：生成粗网格。

通过稀疏细网格，生成指定层数的粗网格。对于结构网格，常用的粗网格生成方法是将细网格沿所有网格方向稀疏一倍，图2以一维为例示意了该粗网格生成方法。图2中，h代表细网格，2h、4h分别代表第1层、第2层粗网格。即，第1层粗网格的每个单元包含2个细网格单元，第2层粗网格中的每个单元包含2个第1层粗网格单元。

第三步：流场初始化。

流场初始化分为根据来流条件初始化和根据给定流场初始化两种方式。当根据来流条件初始化时，将粗网格和细网格上预设计算域中所有网格单元的守恒量赋为来流值。当根据给定流场初始化时，将细网格上预设计算域中所有网格单元的守恒量赋为给定流场值，粗网格上预设计算域中所有网格单元的守恒量根据细网格上所有网格单元的守恒量插值获得。以一维为例，粗网格单元的守恒量满足：

其中，w2h表示粗网格网格单元的守恒量，wh表示细网格网格单元的守恒量；|ω|2h表示粗网格网格单元的体积，|ω|h表示细网格网格单元的体积；2i-1表示沿网格方向第2i-1个单元，2i表示沿网格方向第2i个单元。

第四步：建立对流动态计算域和粘性动态计算域。

(1)当根据来流条件初始化时，对流动态计算域和粘性动态计算域均根据壁面边界建立。携带着流场变化信息的扰动从物体表面产生，受粘性主导的单元也紧邻物体表面，因此，对流动态计算域取壁面的若干层相邻单元作为初始单元。为避免对流动态计算域的增长速度对多重网格法的收敛速率产生影响，最粗层网格取紧邻壁面的1层单元作为初始单元，则对于结构网格框架下的细网格，其对流动态计算域需取壁面的l层相邻单元作为初始单元，l满足：

l＝2^l(2)

其中，l表示粗网格层数，l为正整数，由步骤s1读入的计算设置给定；较粗层网格的初始单元根据细网格的初始单元生成；对于粘性动态计算域，所有网格层均取紧邻壁面的1层单元作为初始单元，以避免粘性动态计算域增长过快而导致无效计算。

(2)当根据给定流场初始化时，对流动态计算域和粘性动态计算域根据给定流场的流动特性建立。对于对流动态计算域，细网格的初始单元为给定流场中的已受扰单元，这些单元的流动特性与来流条件不符，满足：

||w-w∞||/||δw⁽¹⁾||max＞εa,c(3)

其中，w∞表示来流条件的守恒量；||δw⁽¹⁾||max表示对流动态计算域所有单

元的守恒量更新量在第1迭代步的最大值；εa,c表示对流新增阈值，取10^-4≤εa,c≤10^-6。

对于粘性动态计算域，细网格的初始单元为细网格对流动态计算域中的粘性效应主导单元。定义粘性效应衡量参数ψ，令其表示粘性扰动与无粘扰动的质量流量之比。粘性效应主导单元的粘性效应衡量参数应处于较大量级，故粘性动态计算域中的单元应满足：

其中，ψ表示粘性效应衡量参数，为粘性扰动与无粘扰动的质量流量之比；表示第1迭代步的粘性效应衡量参数；i,j,k分别表示网格方向；表示对流通量jacobian矩阵沿i方向的谱半径，表示粘性通量jacobian矩阵沿i方向的谱半径；εa,v表示粘性新增阈值，取10^-3≤εa,v≤10^-2。

所有粗网格层的对流动态计算域依据细网格的对流动态计算域生成，所有粗网格层的粘性动态计算域依据细网格的粘性动态计算域生成。

第五步：分配存储空间。

将所有层网格求解所需存储的数据分为两类；第一类数据是单元的固有信息，采用静态数据结构存储预设计算域中所有网格单元的信息，包括网格坐标和流场变量；第二类数据是与求解更新相关的信息，采用动态数据结构仅存储对流动态计算域内网格单元的信息，包括守恒量更新量和当地迭代步长。

对于第一类数据，采用静态数据结构存储预设计算域中所有网格单元的信息。例如，令i,j,k表示网格方向，imax,jmax,kmax分别表示各网格方向的单元数，对于一块imax×jmax×kmax的网格，若每个单元有10个量需要记录，则可采用一个实数型四维数组a1(10,imax,jmax,kmax)进行存储。

对于第二类数据，采用动态数据结构仅存储对流动态计算域内网格单元的信息。例如，令i,j,k表示网格方向，具有相同j,k标号的单元统称为一行，定义一个特殊数据类型用于存储对流、粘性两类动态计算域中一行的信息，该特殊数据类型包括：a.一条单项链表——在两类动态计算域更新中记录当前迭代步该行单元的i标号范围；b.一个整型二维数组——在两类动态计算域更新时存储上一迭代步该行单元的i标号范围，在两类动态计算域更新后存储当前迭代步该行单元的i标号范围；c.一个实数型二维数组——存储该行位于对流动态计算域中单元的第二类数据，与整型二维数组的i标号范围相对应。对于一块imax×jmax×kmax的网格，若每个单元有10个量需要记录，则采用一个特殊数据类型二维数组a2(jmax,kmax)，a2中实数型二维数组为b(10,i1:i2)，i1、i2分别表示行(j,k)中单元的i标号的最小值与最大值，1＜i1,i2＜imax。

现有多重网格方法中，所有网格层中所有网格单元的所有数据均存储在静态数据结构中。相比之下，本发明提供的多重网格扰动域更新加速方法，将所需存储的数据分为所有单元的固有信息和与求解更新相关的信息两类；其中，对于第一类数据，以静态数据结构存储预设计算域中所有网格单元的信息；对于第二类数据，以动态数据结构仅存储各网格层对流动态计算域中网格单元的信息，并且，第十五步会在每个迭代步及时根据各网格层的对流动态计算域调整第二类数据的存储空间。因此，相比于现有多重网格方法，本发明可有效降低存储空间需求。

第六步：边界条件处理。

根据边界条件类型，为当前层网格的边界虚网格的守恒量赋值。

第七步：残差估计。

流动控制方程表示为：

其中，t表示时间，fc表示对流通量，fv表示粘性通量，qf表示粗网格的强迫函数，qt表示湍流模型方程源项，nf表示单元面的面数，δs表示单元面的面积；式(6)的右端项统称为流动控制方程的残差项，为无粘项，为粘性项；在对流动态计算域中求解残差项的无粘项，在粘性动态计算域中求解残差项的粘性项。

第八步：时间积分。

在对流动态计算域中，采用时间推进格式计算式(6)的左端项；所采用的时间积分方法与现有方法一致。

现有多重网格方法中，最为耗时的残差估计、时间积分两步需对静态的预设计算域内各层网格的所有网格单元执行，约占99％的总求解时间。相比之下，本发明提供的多重网格扰动域更新加速方法，第七步的残差估计仅在各层网格的对流动态计算域内计算残差的无粘项，仅在各层网格的粘性动态计算域内计算残差的粘性项；第八步的时间积分也仅在各层网格的对流动态计算域内执行。因此，本发明可同时降低各层网格的计算量，显著提升多重网格方法的计算效率。

第九步：判断当前网格层是否为细网格；若是，则执行第十步；若否，则执行第十七步。

第十步：判断细网格是否达到收敛条件；若否，则执行第十一步；若是，则执行第十九步。

第十一步：增大对流动态计算域。

对对流动态计算域的所有边界单元逐个判断是否已受到无粘扰动，即对对流动态计算域的所有边界单元逐个判断是否满足：

||δw||/||δw⁽¹⁾||max＞εa,c(7)

其中，||δw||表示对流动态计算域的边界单元在当前迭代步中守恒量更新量的模值；若对流动态计算域的边界单元满足式(7)，则衡量无粘扰动的传播方向，并将该边界单元的紧邻单元中位于传播方向上的单元纳入对流动态计算域，返回第十一步，对对流动态计算域的下一边界单元进行判断；若对流动态计算域的边界单元不满足式(7)，则返回第十一步，对对流动态计算域的下一边界单元进行判断；遍历对流动态计算域的所有边界单元后，执行第十二步。

无粘扰动相对于流动以声速传播，沿无粘扰动的传播方向，无粘扰动的传播速度为正；令q表示单位方向向量，则无粘扰动沿q方向传播表示为：

u·q+a＞0(8)

其中，u表示流动速度矢量，a表示声速；q取对流动态计算域边界单元的格心指向边界单元某一格点的单位向量，若q满足式(8)，则将与对流动态计算域边界单元共享该格点的紧邻单元纳入对流动态计算域。

第十二步：缩小对流动态计算域和粘性动态计算域。

为了在不影响收敛速率的前提下，尽可能地减少单个迭代步的计算量，可以缩小对流动态计算域的边界单元。对对流动态计算域的所有边界单元逐个判断以下条件实现对流动态计算域的缩小：

(1)对对流动态计算域的所有边界单元逐个判断是否已收敛，即对对流动态计算域的所有边界单元逐个判断是否满足：

||δw||/||δw⁽¹⁾||max＜εd(9)

其中，εd表示删除阈值，取εd≤10^-7；若该边界单元满足式(9)，则继续判断(2)；若否，则对对流动态计算域的下一边界单元进行判断；

(2)判断该边界单元是否位于可压缩流动中，即判断该边界单元的马赫数是否大于0.3；若是，则继续判断(3)；若否，则返回判断(1)，对对流动态计算域的下一边界单元进行判断；

(3)判断该边界单元是否位于最上游，即判断该边界单元是否满足：

其中，q取对流动态计算域边界单元的格心指向边界单元紧邻单元格心的单位向量，θd表示上游单元容差角，取5°≤θd≤10°；若该边界单元满足式(10)，则继续判断(4)；若否，则返回判断(1)，对对流动态计算域的下一边界单元进行判断；

(4)判断该边界单元是否不再受对流动态计算域中其他单元的影响，对于超声速无粘流动，则无需判断、自动满足；对于亚声速流动和粘性流动，判断该边界单元是否满足：

||δw||+||δ(δw)||＜εd(11)

其中，||δ(δw)||表示边界单元守恒量更新量受边界单元对流动态计算域中紧邻单元影响的模值；δt表示迭代步长，ccfl表示时间推进格式的cfl数，δrⁱ表示对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响；对于亚声速无粘单元，对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响δrⁱ表示为：

其中，δfc表示对流通量的变化量；下标i+1表示对流动态计算域的边界单元沿正i方向的紧邻单元，下标i-1表示对流动态计算域的边界单元沿负i方向的紧邻单元，δsi+1/2表示对流动态计算域的边界单元沿正i方向单元面的面积；δsi-1/2表示对流动态计算域的边界单元沿负i方向单元面的面积；对于粘性单元，对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响δrⁱ表示为：

若该边界单元满足式(11)和式(12)，则该边界单元已同时满足已收敛、位于可压缩流动中、位于最上游以及不再受所述对流动态计算域中其他单元的影响四个条件，则将该边界单元从对流动态计算域中移除，若该边界单元同时存在于粘性动态计算域中，则将该边界单元从粘性动态计算域中移除，返回第十二步，对对流动态计算域的下一边界单元进行判断；若该边界单元不满足式(11)和式(12)，则返回第十二步，对对流动态计算域的下一边界单元进行判断。遍历对流动态计算域的所有边界单元后，执行第十三步。

第十三步：增大粘性动态计算域。

对粘性动态计算域的所有边界单元逐个判断是否受粘性效应主导，即对粘性动态计算域的所有边界单元逐个判断是否满足：

其中，φ表示粘性新增阈值的缩比因子，||δw||max表示对流动态计算域所有单元的守恒量更新量在当前迭代步的最大值，(||δw||max)min表示第1迭代步至当前迭代步中单步守恒量更新量最大值的最小值；若粘性动态计算域的边界单元满足式(15)和式(16)，则将该边界单元的紧邻单元中位于对流动态计算域中的单元纳入粘性动态计算域，返回第十三步，对粘性动态计算域的下一边界单元进行判断；若粘性动态计算域的边界单元不满足式(15)和式(16)，则返回第十三步，对粘性动态计算域的下一边界单元进行判断；遍历粘性动态计算域的所有边界单元后，执行第十四步。流动控制方程的粘性项呈椭圆型，因此，可以在不增大对流动态计算域的前提下，将已受粘性扰动单元的所有紧邻单元均加入粘性动态计算域。

第十四步：更新粗网格动态计算域。

根据细网格的对流动态计算域更新所有层粗网格的对流动态计算域，根据细网格的粘性动态计算域更新所有层粗网格的粘性动态计算域。对于结构网格，粗网格与细网格间动态计算域的对应关系如图2所示，若某一细网格单元包含于细网格的动态计算域中，则对应于该单元的粗网格单元也会被加入粗网格的动态计算域中。

本发明在第十一步中，针对细网格，根据扰动传播特性仅将受扰动单元加入对流动态计算域中；在第十三步中，通过衡量粘性效应的影响，仅将受粘性效应主导的网格单元加入细网格的粘性动态计算域中；在第十二步中，利用求解收敛特征，在不影响收敛速率的前提下，及时将已收敛单元从对流、粘性两类动态计算域中移除；在第十四步中，根据细网格的动态计算域，更新所有粗网格的对流、粘性两类动态计算域。因此，本发明可始终保持在所有网格层仅对求解未收敛的受扰单元求解，仅在粘性效应主导单元中考虑粘性，从而可有效避免现有多重网格方法求解中的无效计算。

第十五步：重新分配第二类数据的存储空间。

根据对流动态计算域，重新分配所有粗网格和细网格中第二类数据的存储空间。

第十六步：将守恒量插值至更粗网格，返回第六步，进入更粗网格的计算。

守恒量从细网格至粗网格的插值方式与第三步中根据给定流场初始化的插值方式相同。

第十七步：判断是否需要计算更粗网格；若需要计算更粗网格，则返回第十六步；若不需要计算更粗网格，则执行第十八步；

多重网格常用的循环方式如图3所示，左侧子图代表“v”型循环，右侧子图代表“w”型循环，图3中h代表细网格，2h代表第1层粗网格，4h代表第2层粗网格。以“v”型循环为例，若当前计算层为细网格，则第十七步便选择需要计算更粗网格；若当前计算层为最粗网格，则第十七步便选择无需计算更粗网格。

第十八步：将守恒量修正量插值回更细网格，并判断更细网格是否为细网格；若更细网格是细网格，则返回第六步，进入下一个迭代步的计算；若更细网格不是细网格，则返回第六步，进入更细网格的计算。

第十九步：输出结果。

下面分别对马赫数0.755的naca0012翼型无粘绕流问题和马赫数0.729的跨声速rae2822翼型湍流绕流问题进行模拟。

(1)对马赫数0.755的naca0012翼型无粘绕流问题进行模拟。

图4示意了本发明提供的上述多重网格扰动域更新加速方法的对流动态计算域演化过程。图4中，ηc表示对流动态计算域与预设计算域的网格量之比，nmax代表达到收敛态所需的总迭代步数。流场以来流条件初始化，扰动从壁面边界产生。在求解初期，对流动态计算域从壁面逐渐扩展至整个预设计算域。在求解中后期，对流动态计算域开始随着求解收敛逐渐从远场向壁面、从上游向下游收缩，最后仅剩近壁区及其下游区域。从图5a所示的多重网格扰动域更新加速方法的对流动态计算域网格量变化曲线可知，本发明提供的上述多重网格扰动域更新加速方法适用于不同网格层数的多重网格，不同网格层数情况下的对流动态计算域演化过程也基本相同。

图5b对比了本发明提供的上述多重网格扰动域更新加速方法与现有全局更新的多重网格方法的收敛曲线。对比结果表明，对于不同层数的多重网格方法，相比于现有方法，本发明提供的上述多重网格扰动域更新加速方法通过减少单个迭代步的计算量和总迭代步数，均可显著提升数值模拟效率。其中，当本发明采用1层粗网格时，相比于现有多重网格方法可节省52.1％的计算时间；当本发明采用2层粗网格时，相比于现有多重网格方法可节省45.6％的计算时间。

(2)对马赫数0.729的跨声速rae2822翼型湍流绕流问题进行模拟。

图6给出了本发明提供的上述多重网格扰动域更新加速方法的动态计算域演化过程。流场以来流条件初始化，因此，对流动态计算域和粘性动态计算域仍依据壁面建立。对流动态计算域的增大速率远大于粘性动态计算域的增大速率。由图6中1％nmax的情况可知，当对流动态计算域已覆盖整个预设计算域时，粘性动态计算域仍保持在壁面附近，仅包含近壁区及其下游区域。由图7a给出的多重网格扰动域更新加速方法的对流、粘性两类动态计算域网格量变化曲线可知，本发明提供的上述多重网格扰动域更新加速方法适用于粘性流动，对流动态计算域和粘性动态计算域均遵循先增大后收缩的演化规律。

图7b对比了本发明提供的上述多重网格扰动域更新加速方法与现有全局更新的多重网格方法的收敛曲线。对比结果表明，本发明提供的上述多重网格扰动域更新加速方法与现有的多重网格方法一样，均可获得收敛解。得益于对流、粘性两类动态计算域的缩小和仅在局部考虑粘性，本发明提供的上述多重网格扰动域更新加速方法获得了显著的计算效率提升，可节省40.7％的计算时间。

本发明提供的上述飞行器绕流数值模拟的多重网格扰动域更新加速方法，残差估计仅在各层网格的对流动态计算域内计算残差的无粘项，仅在各层网格的粘性动态计算域内计算残差的粘性项，时间积分也仅在各层网格的对流动态计算域内执行，因此，本发明通过在各层网格采用对流、粘性两类动态计算域，可同时降低多重网格方法模拟飞行器绕流问题时各层网格的计算量，显著提升多重网格方法针对飞行器绕流数值模拟的计算效率，从而有效解决数值模拟计算效率难以满足飞行器气动设计日益增长的实际应用需求的问题。将所需存储的数据分为所有单元的固有信息和与求解更新相关的信息两类；其中，对于第一类数据，以静态数据结构存储预设计算域中所有单元的信息；对于第二类数据，以动态数据结构仅存储各网格层对流动态计算域中网格单元的信息，并在每个迭代步及时根据各网格层的对流动态计算域调整第二类数据的存储空间，因此，相比于现有多重网格方法，本发明可有效降低飞行器绕流数值模拟的存储空间需求。针对飞行器流场的细网格，根据扰动传播特性仅将受扰动单元加入对流动态计算域中；通过衡量粘性效应的影响，仅将受粘性效应主导的网格单元加入细网格的粘性动态计算域中；利用求解收敛特征，在不影响收敛速率的前提下，将已收敛单元从对流、粘性两类动态计算域中移除；根据细网格的动态计算域，更新飞行器流场的所有粗网格的对流、粘性两类动态计算域；因此，本发明可始终保持在飞行器流场的所有网格层仅对求解未收敛的受扰单元求解，仅在粘性效应主导单元中考虑粘性，从而可有效避免现有多重网格方法在飞行器绕流数值模拟中的无效计算，解决数值模拟计算效率难以满足飞行器气动设计日益增长的实际应用需求的问题。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。