本发明公开了一种飞行器气动特性预测的自适应网格扰动域更新加速方法，首先，提出自适应网格与动态计算域耦合存储的新型数据格式，从而实现自适应网格上动态计算域信息的高效存储与遍历；其次，建立可标识具有位于脱体强激波区、收敛缓慢等特征单元的自适应探测器；再次，建立既可加速收敛又可最小化计算增量的脱体强激波区各向异性加密网格自适应策略；最后，建立自适应网格加密与动态计算域更新的耦合逻辑及其算法，提出自适应网格扰动域更新加速方法。该方法可进一步提升现有扰动域更新方法的加速效果；相比于传统的全局更新方法，可为飞行器气动特性预测提供同时优化网格与计算域的技术途径。
1.一种飞行器气动特性预测的自适应网格扰动域更新加速方法，其特征在于，包括如下步骤：
S6：判断是否退出当前网格的迭代求解，包括：若当前网格为步骤S1读入网格的粗化网格，以最大迭代步数构建迭代求解退出判据；
判断迭代步数是否达到最大迭代步数，若是，则跳转至步骤S8，若否，则继续步骤S7；
若当前网格即为步骤S1读入的网格，仍以最大迭代步数构建迭代求解退出判据；判断迭代步数是否达到最大迭代步数，若是，则跳转至步骤S9,若否，则继续步骤S7；
若当前网格为步骤S1读入网格的加密网格，以收敛阈值构建迭代求解退出判据；判断所有网格单元守恒量更新量的模值是否均小于收敛阈值，若是，则跳转至步骤S10，若否，则继续步骤S7；
S8：细化当前计算网格，包括：根据网格粗化准则，以步骤S1读入网格为基础，确定出比当前计算网格少稀疏一倍的计算网格；跳转至步骤S3；
依次进行标识待加密单元、加密标识单元、优化加密区域；跳转至步骤S3；
所述步骤S1具体为：读入并存储数值模拟所需的数据，包括计算网格格点坐标、计算网格的边界条件、计算设置参数；
将读入的计算网格按2的倍数稀疏，稀疏倍数N由用户指定；令m 表示读入的计算网格(N)
在某一网格方向上的单元数，若用户指定该方向稀疏N倍，则稀疏i倍后的单元数m 表示为（1）
细网格标号ID 与粗网格标号ID 间的对应关系为（2）
式中，函数floor表示实数向下取整，函数dble表示将整数转化为实数；
2.根据权利要求1所述的飞行器气动特性预测的自适应网格扰动域更新加速方法，其特征在于，所述步骤S3具体为：若由步骤S2进入S3，将网格单元的守恒量以来流值初始化；若由步骤S8或S9进入S3，则通过已完成计算的粗网格的守恒量为待求解的密网格初始化。
3.根据权利要求1所述的飞行器气动特性预测的自适应网格扰动域更新加速方法，其特征在于，所述步骤S4具体为：若首次进入步骤S4，则根据网格的壁面边界分别建立对流动态域和粘性动态域；若非首次进入步骤S4，则根据初始化流场分别建立对流动态域与粘性动态域。
4.根据权利要求1所述的飞行器气动特性预测的自适应网格扰动域更新加速方法，其特征在于，所述步骤S5具体为：离散流动控制方程，则流动控制方程表示为：（4）
其中，t表示时间，W为守恒量，Fc、Fv分别表示对流通量与粘性通量，Qc、Qv分别表示源项的无粘部分与粘性部分， 表示网格单元的体积， 表示单元面的面积，Nf表示单元面的面数；当计算网格为步骤S1的读入网格及其粗化网格时，Nf在二维和三维情况的取值分别为4、6；当计算网格为步骤S1读入网格的加密网格时，Nf在二维和三维情况的取值将分别为不小于4和6的整数，具体取值由网格加密情况决定；
求解流动控制方程包含如下4个子步：S51：调整存储空间——根据步骤S4或步骤S7所定义的对流动态域范围，调整当地时间步长、守恒量更新量的存储空间；
S52：边界条件处理——遍历步骤S1读入计算网格的所有边界，根据每一边界的条件，为其虚网格的守恒量赋值；
式（4）等号右端项中，对流通量项与源项无粘部分统称为流动控制方程残差的无粘项，粘性通量与源项粘性部分统称为流动控制方程残差的粘性项；遍历对流动态域中的所有单元，计算每一单元的残差无粘项；遍历粘性动态域中的所有单元，计算每一单元的残差粘性项；
在每一单元计算对流通量与粘性通量时，先采用重构格式确定单元界面的左、右状态值，再采用对流离散格式确定单元界面处的通量值；对于结构网格，重构格式按网格标号对模板单元插值；对于非结构网格，重构格式则根据相邻单元的流动梯度与间距进行插值；
通过式（4）等号左端项，确定守恒量更新量并更新守恒量W。
5.根据权利要求1所述的飞行器气动特性预测的自适应网格扰动域更新加速方法，其特征在于，所述步骤S7包含如下4个子步：S71：增大对流动态域；
遍历对流动态域的边界单元，对所有单元执行如下步骤：通过守恒量更新量的模值判断该边界单元是否受到扰动；若该边界单元受到扰动，则将其可能受到扰动的紧邻单元加入对流动态域；
遍历对流动态域的边界单元，对所有单元判断如下条件：（1）周围不存在新增受扰单元，（2）求解已收敛，（3）位于最上游，（4）不再影响对流动态域中其他单元的求解；（5）不再受对流动态域中其他单元的影响；若上述5个条件均满足，则将该单元从对流动态域中移除；
遍历粘性动态域的边界单元，通过守恒量判断该边界单元是否受到粘性扰动；若受到粘性扰动，则将其所有紧邻单元加入粘性动态域；
遍历粘性动态域的边界单元，若其满足以下2个条件：（1）周围不存在新增的粘性效应主导单元；（2）不再是粘性效应主导单元；则将该单元从粘性动态域中移除。
6.根据权利要求1所述的飞行器气动特性预测的自适应网格扰动域更新加速方法，其特征在于，所述步骤S8中采用式（1）定义的网格粗化准则。
7.根据权利要求1所述的飞行器气动特性预测的自适应网格扰动域更新加速方法，其特征在于，所述步骤S9包含如下3个子步：S91：标识待加密单元；
遍历对流动态域中所有单元，若其满足以下2个条件：（1）位于最上游，（2）收敛缓慢，则将该单元标识为待加密单元；
遍历步骤S91标记的所有待加密单元，以流动梯度为衡量标准，将流动梯度大的方向作为网格加密方向；将被标记网格沿确定的加密方向二等分，并存储新增网格的格点坐标、格心坐标、面向量和壁面距；
遍历步骤S92标记的所有被加密单元，根据被加密单元与其相邻单元的体积比与加密方向长度，对被加密单元的相邻单元进行加密，确保网格平滑过渡。
飞行器气动特性预测的自适应网格扰动域更新加速方法
技术领域
[0001] 本发明涉及计算流体力学领域，特别涉及飞行器气动特性预测的自适应网格扰动域更新加速方法。
背景技术
[0002] 计算流体力学数值模拟已成为飞行器气动特性预测的重要手段之一。数值模拟的首要步骤便是指定计算域并将其离散为待求解单元，即生成计算网格。自适应网格技术便
是一种根据几何信息与流动特性自动加密或稀疏计算网格的网格生成技术。在飞行器气动
特性模拟中引入自适应网格技术，可在大流动梯度区保证数值精度的同时，显著降低计算
网格的单元数量。现有基于自适应网格的数值模拟方法均采用静态计算域全局更新求解。
由于忽略了数值/物理扰动传播与流动控制方程求解收敛的特征，这种全局更新的求解方
式会产生不同程度的无效计算，会因对未受扰区域迭代求解而引入更多数值误差，会因对
已收敛区域持续更新而延迟收敛。
[0003] 专利文献 ZL201810250654.8 提出了一种扰动区域更新方法可通过消除传统方法中的无效计算而获得计算效率的显著提升。这种扰动区域更新方法正好可用于避免现有
基于自适应网格数值模拟方法的不足。目前，扰动区域更新方法已在各类可压缩流动的数
值模拟中得到了很好地验证，但其也存在计算效率进一步提升的空间。例如，该方法仅能通
过优化计算域提升计算效率，却不具备优化网格质量的能力。一方面，网格分布不合理便无
法以最小网格量获得最高数值精度。另一方面，若模拟流场中存在脱体强激波区，该区域便
可能因网格质量未满足要求而导致动态计算域无法收缩，从而使扰动区域更新方法无法发
挥其真正的加速效果。而自适应网格技术正好是在求解中通过改善网格疏密度实现收敛加
速或计算精度提升的有效方法。因此，令两种方法优势互补，建立自适应网格技术与扰动区
域更新方法的高效耦合方法是进一步提升数值模拟计算效率的有效途径。
发明内容
[0004] 针对现有基于自适应网格技术的数值模拟方法存在无效计算、现有扰动区域更新方法无法自动优化网格限制加速效果的问题，本发明提出一种飞行器气动特性预测的自适
应网格扰动域更新加速方法。通过扰动域区域更新方法解决现有基于自适应网格技术的数
值模拟方法中存在无效计算的问题，通过自适应网格技术解决现有扰动区域更新方法中因
无法自动优化网格而限制加速效果的问题，并着重提出两种方法耦合时所涉及的新算法与
新数据格式。首先，提出自适应网格与动态计算域耦合存储的新型数据格式，从而实现自适
应网格上动态计算域信息的高效存储与遍历。其次，建立可标识具有位于脱体强激波区、收
敛缓慢等特征单元的自适应探测器。再次，建立既可加速收敛又可最小化计算增量的脱体
强激波区各向异性加密网格自适应策略。最后，建立自适应网格加密与动态计算域更新的
耦合逻辑及其算法，提出自适应网格扰动域更新加速方法。
[0005] 一种飞行器气动特性预测的自适应网格扰动域更新加速方法，包括如下步骤：
[0006] S1：读入数据；
[0007] S2：粗化网格；
[0008] S3：初始化流场；
[0009] S4：建立动态计算域；
[0010] S5：求解控制方程；
[0011] S6：判断是否退出当前网格的迭代求解；
[0012] 若当前网格为步骤S1读入网格的粗化网格，以最大迭代步数构建迭代求解退出判据；判断迭代步数是否达到最大迭代步数，若是，则跳转至步骤S8，若否，则继续步骤S7；
[0013] 若当前网格即为步骤S1读入的网格，仍以最大迭代步数构建迭代求解退出判据；判断迭代步数是否达到最大迭代步数，若是，则跳转至步骤S9,若否，则继续步骤S7；
[0014] 若当前网格为步骤S1读入网格的加密网格，以收敛阈值构建迭代求解退出判据；判断所有网格单元守恒量更新量的模值是否均小于收敛阈值，若是，则跳转至步骤S10，若
否，则继续步骤S7；
[0015] S7：更新动态计算域；
[0016] S8：细化当前计算网格；
[0017] 根据网格粗化准则，以步骤S1读入网格为基础，确定出比当前计算网格少稀疏一倍的计算网格；跳转至步骤S3；
[0018] S9：自适应加密网格；
[0019] 依次进行标识待加密单元、加密标识单元、优化加密区域；跳转至步骤S3；
[0020] S10：输出结果，结束。
[0021] 进一步的，所述步骤S1具体为：读入并存储数值模拟所需的数据，包括计算网格格点坐标、计算网格的边界条件、计算设置参数。
[0022] 进一步的，所述步骤S2具体为：
[0023] 将读入的计算网格按2的倍数稀疏，稀疏倍数N由用户指定；令m(0)表示读入的计算(N)
网格在某一网格方向上的单元数，若用户指定该方向稀疏N倍，则稀疏i倍后的单元数m 表
示为
[0024] （1）
[0025] 式中，函数mod表示计算m(0)/2N的余数；
[0026] 细网格标号ID(0)与粗网格标号ID(N)间的对应关系为
[0027] （2）
[0028] 式中，函数floor表示实数向下取整，函数dble表示将整数转化为实数；
[0029] 粗网格标号与细网格标号间的对应关系为
[0030] （3）
[0031] 式中，函数ceiling表示向上取整。
[0032] 进一步的，所述步骤S3具体为：
[0033] 若由步骤S2进入，将网格单元的守恒量以来流值初始化；若由步骤S8或S9进入，则通过已完成计算的粗网格的守恒量为待求解的密网格初始化。
[0034] 进一步的，所述步骤S4具体为：
[0035] 若首次进入步骤S4，则根据网格的壁面边界分别建立对流动态域和粘性动态域；若非首次进入步骤S4，则根据初始化流场分别建立对流动态域与粘性动态域。
[0036] 进一步的，所述步骤S5具体为：
[0037] 离散流动控制方程表示为：
[0038] （4）
[0039] 其中，t表示时间，W为守恒量，Fc、Fv分别表示对流通量与粘性通量，Qc、Qv分别表示源项的无粘部分与粘性部分， 表示网格单元的体积， 表示单元面的面积，Nf表示单
元面的面数；当计算网格为步骤S1的读入网格及其粗化网格时，Nf在二维和三维情况的取
值分别为4、6；当计算网格为步骤S1读入网格的加密网格时，Nf在二维和三维情况的取值将
分别为不小于4和6的整数，具体取值由网格加密情况决定；
[0040] 求解控制方程包含如下4个子步：
[0041] S51：调整存储空间——根据步骤S4或步骤S7所定义的对流动态域范围，调整当地时间步长、守恒量更新量的存储空间；
[0042] S52：边界条件处理——遍历步骤S1读入计算网格的所有边界，根据每一边界的条件，为其虚网格的守恒量赋值；
[0043] S53：残差估计；
[0044] 估计式（4）等号右端项；式（4）中，对流通量项与源项无粘部分统称为控制方程残差的无粘项，粘性通量与源项粘性部分统称为控制方程残差的粘性项；遍历对流动态域中
的所有单元，计算每一单元的残差无粘项；遍历粘性动态域中的所有单元，计算每一单元的
残差粘性项；
[0045] 在每一单元计算对流通量与粘性通量时，先采用重构格式确定单元界面的左、右状态值，再采用对流离散格式确定单元界面处的通量值；对于结构网格，重构格式按网格标
号对模板单元插值；对于非结构网格，重构格式则根据相邻单元的流动梯度与间距进行插
值；
[0046] S54：时间积分；
[0047] 通过时间离散格式近似式（2）等号左端项，确定守恒量更新量并更新守恒量W。
[0048] 进一步的，所述步骤S7包含如下4个子步：
[0049] S71：增大对流动态域；
[0050] 遍历对流动态域的边界单元，对所有单元执行如下步骤：通过守恒量更新量的模值判断该边界单元是否受到扰动；若该边界单元受到扰动，则将其可能受到扰动的紧邻单
元加入对流动态域；
[0051] S72：缩小对流动态域；
[0052] 遍历对流动态域的边界单元，对所有单元判断如下条件：（1）周围不存在新增受扰单元，（2）求解已收敛，（3）位于最上游，（4）不再影响对流动态域中其他单元的求解；（5）不
再受对流动态域中其他单元的影响；若上述5个条件均满足，则将该单元从对流动态域中移
除；
[0053] S73：增大粘性动态域；
[0054] 遍历粘性动态域的边界单元，通过守恒量判断该边界单元是否受到粘性扰动；若受到粘性扰动，则将其所有紧邻单元加入粘性动态域；
[0055] S74：缩小粘性动态域；
[0056] 遍历粘性动态域的边界单元，若其满足以下2个条件：（1）周围不存在新增的粘性效应主导单元；（2）不再是粘性效应主导单元；则将该单元从粘性动态域中移除。
[0057] 进一步的，所述步骤S8中采用式（1）定义的网格粗化准则。
[0058] 进一步的，所述步骤S9包含如下3个子步：
[0059] S91：标识待加密单元；
[0060] 遍历对流动态域中所有单元，若其满足以下2个条件：（1）位于最上游，（2）收敛缓慢，则将该单元标识为待加密单元；
[0061] S92：加密标识单元；
[0062] 遍历步骤S91标记的所有待加密单元，以流动梯度为衡量标准，将流动梯度较大的方向作为网格加密方向；将被标记网格沿确定的加密方向二等分，并存储新增网格的格点
坐标、格心坐标、面向量和壁面距；
[0063] S93：优化加密区域；
[0064] 遍历步骤S92标记的所有被加密单元，根据被加密单元与其相邻单元的体积比与加密方向长度，对被加密单元的相邻单元进行加密，确保网格平滑过渡。
[0065] 本发明与现有技术相比所具有的有益效果：
[0066] 本发明所提飞行器气动特性预测的自适应网格扰动域更新方法可进一步提升现有扰动域更新方法的加速效果；相比于传统的全局更新方法，可为飞行器气动特性预测提
供同时优化网格与计算域的技术途径。
[0067] 本发明提高计算效率的原理主要体现在三个方面。其一，步骤S2将用户的输入网格粗化，使得本发明在根据用户输入的网格求解前，可先在粗化网格上求解，以便采用更大
迭代步长使流场接近定常解。其二，步骤S5仅在动态计算域中求解流动控制方程，仅在对流
动态域中估计残差的无粘项和时间积分，仅在粘性动态域中估计残差的粘性项，从而通过
优化计算域有效消除了传统全局更新方法中的无效计算。其三，步骤S9通过自适应网格加
密技术，针对对流动态域上游单元中的难以收敛单元，提高其在流动大梯度方向的网格分
辨率。一方面，步骤S9可通过优化网格提高本发明方法的数值精度；另一方面，其又可改善
难以收敛区域的收敛性，以避免因网格质量不佳影响动态计算域收缩从而影响扰动域更新
方法加速效果的问题。
附图说明
[0068] 图1为本发明自适应网格扰动域更新加速方法的流程图；
[0069] 图2(a)为本发明对NACA0012翼型问题生成的自适应加密网格；
[0070] 图2(b)为本发明对钝头体问题生成的自适应加密网格。
具体实施方式
[0071] 图1为本发明实施的流程图，下面进行具体介绍。本发明提供了一种飞行器气动特性预测的自适应网格扰动域更新加速方法，具体包括如下步骤：
[0072] S1：数据读入；
[0073] 分配网格格点坐标、网格边界条件和计算设置参数的存储空间，并从用户提供的输入文件中读入网格的节点坐标、边界条件，以及计算设置。
[0074] S2：粗化网格；
[0075] 将读入的计算网格按2的倍数稀疏，稀疏倍数N由用户指定，通常系数2 3倍。令m(0)~
表示读入的计算网格在某一网格方向上的单元数，若用户指定该方向稀疏N倍，则稀疏N倍
(N)
后的单元数m 可表示为
[0076] （5）
[0077] 式中，函数mod表示计算m(0)/2N的余数。
[0078] 细网格标号ID(0)与粗网格标号ID(N)间的对应关系为
[0079] （6）
[0080] 式中，函数floor表示实数向下取整，函数dble表示将整数转化为实数；
[0081] 粗网格标号与细网格标号间的对应关系为
[0082] （7）
[0083] 式中，函数ceiling表示向上取整。
[0084] S3：流场初始化；
[0085] 若由步骤S2进入，流场还未被在任何网格上求解过，将网格单元的守恒量以来流值初始化。若由步骤S8或S9进入，则将已完成计算的粗网格的解插值到待求解的密网格上，
为密网格上的流场求解进行初始化。
[0086] S4：建立动态计算域；
[0087] 根据壁面边界建立动态计算域。取所有子网格中壁面边界的10层相邻单元作为对流动态域与非定常动态域的初始单元。
[0088] S5：求解控制方程；
[0089] 离散流动控制方程可表示为：
[0090] （8）
[0091] 其中，t为时间，W为守恒量，Fc、Fv分别表示对流通量与粘性通量，Qc、Qv分别表示源项的无粘部分与粘性部分， 表示网格单元的体积， 表示单元面的面积，Nf表示单元
面的面数。当计算网格为结构网格时，即为由步骤S1读入的网格及其粗化网格时，Nf在二
维、三维情况的取值分别为4和6。当计算网格为步骤S1读入网格的加密网格时，Nf在二维、
三维情况的取值分别为不小于4和6的整数，具体取值由网格加密情况决定。
[0092] 步骤S5包含如下4个子步：
[0093] S51：调整存储空间——根据步骤S4或步骤S7所定义的对流动态域范围，调整当地时间步长、守恒量更新量的存储空间；
[0094] S52：边界条件处理——遍历步骤S1读入计算网格的所有边界，根据计算网格的边界条件，为每一边界的虚网格的守恒量赋值；
[0095] S53：残差估计；
[0096] 采用空间离散格式估计式（8）的等号右端项。首先，遍历对流动态域中的所有单元，计算每一单元在式（8）中的对流通量项与源项的无粘部分；其次，遍历粘性动态域中的
所有单元，计算每一单元在式（8）中的粘性通量项与源项的粘性部分。在每一单元计算对流
通量与粘性通量时，先采用重构格式确定单元界面的左、右状态值，再采用对流离散格式通
过界面的左、右状态值确定出单元界面处的通量值。对于结构网格，重构格式按网格标号对
模板单元插值；对于非结构网格，重构格式则根据相邻单元的流动梯度与间距进行插值。
[0097] S54：时间积分；
[0098] 采用时间离散格式离散式（8）的等号左端项，确定守恒量更新量并更新守恒量W。
[0099] S6：判断是否退出当前网格的迭代求解；
[0100] 若当前网格为步骤S1读入网格的粗化网格，以最大迭代步数构建迭代求解退出判据；判断迭代步数是否达到最大迭代步数，若是，则跳转至步骤S8，若否，则继续步骤S7；
[0101] 若当前网格即为步骤S1读入的网格，仍以最大迭代步数构建迭代求解退出判据；判断迭代步数是否达到最大迭代步数，若是，则跳转至步骤S9,若否，则继续步骤S7；
[0102] 若当前网格为步骤S1读入网格的加密网格，以收敛阈值构建迭代求解退出判据；判断所有网格单元守恒量更新量的模值是否均小于收敛阈值，若是，则跳转至步骤S10，若
否，则继续步骤S7。
[0103] S7：更新动态计算域；
[0104] 步骤S7包含如下4个子步：
[0105] S71：增大对流动态域；
[0106] 遍历所有对流动态域的边界单元，对任一边界单元：
[0107] （1）通过守恒量更新量的模值，判断该边界单元是否受到扰动；令||ΔW||表示守‑5
恒量更新量的模值，εa,c表示给定的对流新增阈值，取10 ，则||ΔW||>εa,c代表单元已受到
扰动；
[0108] （2）若该单元受到扰动，当处于亚声速流动中，将其所有紧邻单元加入对流动态域和非定常动态域；当处于超声速流动中，令q为该单元格心指向该单元格点的单位向量，u为
流动速度矢量，a为声速，则扰动沿q方向传播的条件可表示为
[0109] （9）
[0110] 若该单元格心指向某一格点的单位向量满足式（9）则将与包含该格点的紧邻单元加入对流动态域和非定常动态域。
[0111] S72：缩小对流动态域；
[0112] 遍历所有对流动态域的边界单元，对任一边界单元，若其满足以下5个条件则将该单元从对流动态域中移除。具体判断条件如下：
[0113] （1）周围不存在新增无粘受扰单元；
[0114] 若周围单元均满足||ΔW||<εa,c，则代表周围不存在新增无粘受扰单元。
[0115] （2）求解已收敛；
[0116] 令εd表示给定的删除阈值，取10‑7，则待删除单元求解已收敛可描述为||ΔW||<εd。
[0117] （3）位于最上游；
[0118] 令q表示待删除单元格心指向其紧邻单元格心的单位向量，若待删除单元位于最上游，则其所有处于对流动态域中的紧邻单元均应满足
[0119] （10）
[0120] 式中，u表示速度矢量；θd表示上游单元容差角，超声速流动取10°，亚声速流动取45°。
[0121] （4）不再影响对流动态域中其他单元的求解；
[0122] 以马赫数为判断，马赫数大于0.3的可压缩流动位于最上游，则可不再影响对流动态域中其他单元的求解；而马赫数小于0.3的不可压缩流动，则不可移除。
[0123] （5）不再受对流动态域中其他单元的影响；
[0124] 若考虑紧邻单元对待删除单元守恒量更新量的影响后，其仍满足收敛条件，则可认为该单元不再受其他单元的影响，即满足
[0125] （11）
[0126] 其中
[0127]
[0128] （12）
[0129] 式中，Δt表示迭代步长，CCFL表示时间推进格式的CFL数，|Ω|表示网格单元的体i
积；I, J, K表示网格方向；ΔR 表示对流动态域中的紧邻单元对对流动态域边界单元的残
差沿i方向的影响，i=I, J, K；ΔW表示守恒量更新量； 表示对流通量变化量，即当前
步与上一步对流通量之差；下标i±1表示对流动态域边界单元沿正、负i方向的紧邻单元；
下标i±1/2表示对流动态域边界单元沿正、负i方向的单元面； 表示对流通量Jacobian
矩阵沿i方向的谱半径。
[0130] S73：增大粘性动态域；
[0131] 遍历粘性动态域的边界单元，对任一边界单元：
[0132] （1）通过守恒量，判断该边界单元是否受到粘性扰动；令εa,v表示给定的粘性新增单元，粘性效应主导单元应满足：
[0133] （13）
[0134] 式中，Ψ为粘性效应衡量参数， 为第1个迭代步紧邻单元粘性效应衡量参数的最大值； 表示式粘性通量Jacobian矩阵沿i方向的谱半径。
[0135] （2）若受到粘性扰动，则将该边界单元的所有紧邻单元加入粘性动态域。
[0136] S74：缩小粘性动态域；
[0137] 遍历粘性动态域的边界单元，若其满足以下2个条件，则将该单元从粘性动态域中移除。具体条件如下：
[0138] （1）周围不存在新增的粘性效应主导单元；
[0139] 判断相邻两层单元是否满足式（13），若均不满足则表示周围不存在新增的粘性效应主导单元；
[0140] （2）不再是粘性效应主导单元；
[0141] 若删除单元不满足式（13），则表示其不再是粘性效应主导单元。
[0142] 记录流场中速度梯度与压力梯度的平均值，将重叠区域中速度、压力梯度大于平均值的单元标记为大流动梯度单元，并记录其网格节点坐标。
[0143] S8：细化当前计算网格；
[0144] 根据式（5）定义的网格粗化准则，以步骤S1读入网格为基础，确定出比当前计算网格少稀疏一倍的计算网格；
[0145] 例如，若当前网格为读入网格稀疏N倍后的网格，则下一次计算的应为稀疏N‑1倍的网格，即单元数为
[0146] （14）
[0147] S9：自适应加密网格；
[0148] 步骤S9包含如下3个子步：
[0149] S91：标识待加密单元；
[0150] 遍历对流动态域中所有单元，若其满足以下2个条件，则将该单元标识为待加密单元：
[0151] （1）位于最上游；
[0152] 是否位于最上游可通过式（10）判断；
[0153] （2）收敛缓慢；
[0154] 计算该单元守恒量更新量随迭代步的下降率，若其小于0.5倍的最大下降率，则可认为其收敛缓慢；
[0155] S92：加密标识单元；
[0156] 遍历步骤S91标记的所有待加密单元，以流动梯度为衡量标准，将流动梯度较大的方向作为网格加密方向；将被标记网格沿确定的加密方向二等分，并存储新增网格的格点
坐标、格心坐标、面向量和壁面距；
[0157] S93：优化加密区域；
[0158] 遍历步骤S92标记的所有被加密单元，根据被加密单元与其相邻单元的体积比与加密方向长度，对被加密单元的相邻单元进行加密，确保网格平滑过渡。
[0159] 跳转至步骤S3。
[0160] S10：输出结果，结束。
[0161] 下面结合附图和实施例进一步描述本发明，应该理解，以下所述实施例旨在便于对本发明的理解，而对其不起任何限定作用。
[0162] 实施例1：图2(a)为本发明模拟NACA0012翼型问题时生成的自适应加密网格，图2(b)为本发明模拟钝头体问题时生成的自适应加密网格。图中，网格在近壁区依据壁面曲率
实现了自适应加密，加密区域可保持光滑过渡。说明本发明方法可实现准确标识需要加密
的区域，并合理过渡加密区域，为数值模拟提供分布合理的网格。
[0163] 对于本领域的普通技术人员来说，在不脱离本发明创造构思的前提下，还可以对本发明的实施例做出若干变型和改进，这些都属于本发明的保护范围。
