一种飞行器动态气动特性模拟的迭代推进扰动域更新方法与流程
本发明涉及计算流体力学技术领域，尤其涉及一种飞行器动态气动特性模拟的迭代推进扰动域更新方法。
背景技术：
动态气动特性是飞行器设计的关键性能参数。得益于计算机硬件性能的飞速发展，非定常流动数值模拟已广泛应用于各类飞行器的动态气动特性预测中。然而，由于流动问题复杂导致计算量大、数值模拟方法求解效率低等问题，飞行器动态气动特性的数值模拟仍需耗费大量计算资源和超长求解周期，存在数值模拟计算效率难以满足飞行器气动设计日益增长的实际应用需求的问题。因此，提出适用于飞行器动态气动特性预测的高效数值模拟技术一直是计算流体力学的研究热点。
当前，迭代推进方法是飞行器动态气动特性数值模拟的主流方法之一，例如基于二阶后向差分的双时间步法。该方法通过在飞行器动态气动特性数值模拟的时间推进求解中引入内迭代，既可保证飞行器动态气动特性的数值模拟精度，又可将成熟的飞行器静态气动特性数值模拟方法推广至动态气动特性的数值模拟中。
目前，适用于迭代推进类飞行器动态气动特性数值模拟的加速技术主要通过迭代格式隐式化、空间并行化和网格生成合理化三种途径实现。迭代格式隐式化可打破迭代步长的稳定性条件限制，通过大幅增大迭代步长减少内迭代所需的步数。空间并行化是通过将飞行器流场网格划分为多个网格块，即将计算任务分解成多个子任务并行处理，从而可降低串行处理的计算量。网格生成合理化则是通过合理控制飞行器流场网格的疏密度，从而在保证空间精度的同时可通过减小网格量来降低计算量。
上述技术从离散格式、求解算法、网格等方面实现了飞行器动态气动特性数值模拟的加速，但在每一个物理时刻、每个迭代步中却在飞行器流场网格的静态计算域中对所有网格单元进行全局更新求解。这种更新方式存在两点弊端：一方面，飞行器流场网格的静态计算域对动态气动特性数值模拟的精度与效率至关重要，但其选取却完全依赖于经验，所选计算域不足会导致求解失败，所选计算域过大又浪费计算资源；另一方面，在飞行器流场网格的静态计算域中数值模拟无法利用求解中所需求解范围随流动演化与求解收敛不断变化的特点，会造成大量无效计算，甚至引入更多数值误差。因此，若能消除现有迭代推进类飞行器动态气动特性数值模拟方法中因采用飞行器流场网格的静态计算域模拟所造成的无效计算，必将显著降低飞行器动态气动特性数值模拟所需的计算量，从而提高飞行器动态气动特性数值模拟的计算效率，有效解决数值模拟计算效率难以满足飞行器气动设计日益增长的实际应用需求的问题。
技术实现要素：
有鉴于此，本发明提供了一种飞行器动态气动特性模拟的迭代推进扰动域更新方法，用以解决现有飞行器动态气动特性数值模拟方法中存在大量无效计算的问题。
本发明提供的一种飞行器动态气动特性模拟的迭代推进扰动域更新方法，包括如下步骤：
s1：读入数据，包括飞行器流场的网格、预设计算域、边界条件和计算设置；
s2：根据来流条件或根据给定流场，对所述预设计算域中的流场进行初始化；
s3：根据流场初始化方式，建立非定常动态计算域、对流动态计算域和粘性动态计算域；
s4：将求解所需存储的数据分为两类；第一类数据是单元的固有信息，采用静态数据结构存储所述预设计算域中所有网格单元的信息，包括网格坐标以及当前时刻和上两个时刻的流场变量；第二类数据是与求解更新相关的信息，采用动态数据结构仅存储所述对流动态计算域内网格单元的信息，包括守恒量更新量和当地迭代步长；
s5：根据边界条件类型，为边界虚网格的守恒量赋值；
s6：将流动控制方程的残差项分为无粘项与粘性项两类；其中，所述无粘项包括对流通量和迭代推进法的源项，所述无粘项在对流动态计算域中计算；所述粘性项包括粘性通量和湍流模型方程源项，所述粘性项在粘性动态计算域中计算；
s7：在所述对流动态计算域中，求解守恒量更新量，更新当前时刻的流场变量；
s8：对所述对流动态计算域的所有边界单元逐个判断是否已受到无粘扰动；若是，则执行步骤s9；若否，则返回步骤s8，对所述对流动态计算域的下一边界单元进行判断；遍历所述对流动态计算域的所有边界单元后，执行步骤s10；
s9：衡量无粘扰动的传播方向，并将该边界单元的紧邻单元中位于传播方向上的单元分别纳入所述对流动态计算域和所述非定常动态计算域，返回步骤s8，对所述对流动态计算域的下一边界单元进行判断；
s10：判断所述非定常动态计算域是否有所增大；若是，则执行步骤s11和步骤s12；若否，则执行步骤s12；
s11：根据所述非定常动态计算域，重新分配所述第二类数据的存储空间；
s12：对所述对流动态计算域的所有边界单元逐个判断是否同时满足已收敛、位于可压缩流动中、位于最上游以及不再受所述对流动态计算域中其他单元的影响四个条件；若是，则执行步骤s13；若否，则返回步骤s12，对所述对流动态计算域的下一边界单元进行判断；遍历所述对流动态计算域的所有边界单元后，执行步骤s15；
s13：将该边界单元从所述对流动态计算域中移除，并判断该边界单元是否存在于所述粘性动态计算域中；若是，则执行步骤s14；若否，则返回步骤s12，对所述对流动态计算域的下一边界单元进行判断；
s14：将该边界单元从所述粘性动态计算域中移除，返回步骤s12，对所述对流动态计算域的下一边界单元进行判断；
s15：对所述粘性动态计算域的所有边界单元逐个判断是否受粘性效应主导；若是，则执行步骤s16；若否，则返回步骤s15，对所述粘性动态计算域的下一边界单元进行判断；遍历所述粘性动态计算域的所有边界单元后，执行步骤s17；
s16：将该边界单元的紧邻单元中位于所述对流动态计算域中的单元纳入所述粘性动态计算域，返回步骤s15，对所述粘性动态计算域的下一边界单元进行判断；
s17：判断内迭代是否达到收敛条件；若是，则执行步骤s18；若否，则返回步骤s5，进入内迭代下一迭代步的计算，重复执行步骤s5～s17；
s18：利用物理时间导数衡量非定常效应的影响，对所述非定常动态计算域的所有边界单元逐个判断是否不再受非定常效应主导；若是，则执行步骤s19；若否，则返回步骤s18，对所述非定常动态计算域的下一边界单元进行判断；遍历所述非定常动态计算域的所有边界单元后，执行步骤s20；
s19：将该边界单元从所述非定常动态计算域中移除，返回步骤s18，对所述非定常动态计算域的下一边界单元进行判断；
s20：判断所述非定常动态计算域是否有所缩小；若是，则执行步骤s21和步骤s22；若否，则执行步骤s22；
s21：根据所述非定常动态计算域，重新分配所述第二类数据的存储空间；
s22：更新上两个时刻的流场变量，将所述对流动态计算域重置为所述非定常动态计算域的范围，将所述粘性动态计算域重置为所述对流动态计算域中粘性效应起主导作用的区域；
s23：判断是否达到指定的末态物理时刻；若是，则执行步骤s24；若否，则返回步骤s5，进入下一物理时刻的计算，重复执行步骤s5～s23；
s24：输出结果。
在一种可能的实现方式中，在本发明提供的上述飞行器动态气动特性模拟的迭代推进扰动域更新方法中，步骤s2，根据来流条件或根据给定流场，对所述预设计算域中的流场进行初始化，具体包括：
根据来流条件初始化，将所述预设计算域中所有网格单元的守恒量赋为来流值；
根据给定流场初始化，将所述预设计算域中所有网格单元的守恒量赋为给定流场值。
在一种可能的实现方式中，在本发明提供的上述飞行器动态气动特性模拟的迭代推进扰动域更新方法中，步骤s3，根据流场初始化方式，建立非定常动态计算域、对流动态计算域和粘性动态计算域，具体包括：
根据来流条件初始化，非定常动态计算域和对流动态计算域取壁面的若干层相邻单元作为初始单元，粘性动态计算域取紧邻壁面的1层单元作为初始单元；
根据给定流场初始化，对流动态计算域的初始单元为给定流场中流动特性与来流条件不符的已受扰单元，满足：
||w-w∞||/||δw(1)||max＞εa,c(1)
其中，w表示守恒量；w∞表示来流条件的守恒量；||δw(1)||max表示对流动态计算域所有单元的守恒量更新量在第1物理时刻的第1迭代步的最大值；εa,c表示对流新增阈值，取10-4≤εa,c≤10-6；非定常动态计算域的初始单元与对流动态计算域的初始单元一致；粘性动态计算域的初始单元为对流动态计算域中的粘性效应主导单元，满足：
其中，ψ表示粘性效应衡量参数，为粘性扰动与无粘扰动的质量流量之比；
在一种可能的实现方式中，在本发明提供的上述飞行器动态气动特性模拟的迭代推进扰动域更新方法中，步骤s6，将流动控制方程的残差项分为无粘项与粘性项两类；其中，所述无粘项包括对流通量和迭代推进法的源项，所述无粘项在对流动态计算域中计算；所述粘性项包括粘性通量和湍流模型方程源项，所述粘性项在粘性动态计算域中计算，具体包括：
流动控制方程表示为：
其中，τ表示虚拟时间；δt表示物理时间步长；fc表示对流通量，fv表示粘性通量，qt表示湍流模型方程源项，
在一种可能的实现方式中，在本发明提供的上述飞行器动态气动特性模拟的迭代推进扰动域更新方法中，步骤s8，对所述对流动态计算域的所有边界单元逐个判断是否已受到无粘扰动，具体包括：
对所述对流动态计算域的所有边界单元逐个判断是否满足：
||δw||/||δw(1)||max＞εa,c(5)
其中，||δw||表示对流动态计算域的边界单元在当前迭代步中守恒量更新量的模值；
步骤s9，衡量无粘扰动的传播方向，并将该边界单元的紧邻单元中位于传播方向上的单元纳入所述对流动态计算域和所述非定常动态计算域，具体包括：
无粘扰动相对于流动以声速传播，沿无粘扰动的传播方向，无粘扰动的传播速度为正；令q表示单位方向向量，则无粘扰动沿q方向传播表示为：
u·q+a＞0(6)
其中，u表示流动速度矢量，a表示声速；q取对流动态计算域边界单元的格心指向边界单元某一格点的单位向量，若q满足式(6)，则将与对流动态计算域边界单元共享该格点的紧邻单元纳入对流动态计算域和非定常动态计算域。
在一种可能的实现方式中，在本发明提供的上述飞行器动态气动特性模拟的迭代推进扰动域更新方法中，步骤s12，对所述对流动态计算域的所有边界单元逐个判断是否同时满足已收敛、位于可压缩流动中、位于最上游以及不再受所述对流动态计算域中其他单元的影响四个条件；若是，则执行步骤s13；若否，则返回步骤s12，对所述对流动态计算域的下一边界单元进行判断；具体包括如下步骤：
s121：对所述对流动态计算域的所有边界单元逐个判断是否满足：
||δw||/||δw(1)||max＜εd(7)
其中，εd表示删除阈值，取εd≤10-7；若是，则执行步骤s122；若否，则返回步骤s121，对所述对流动态计算域的下一边界单元进行判断；
s122：判断该边界单元的马赫数是否大于0.3；若是，则执行步骤s123；若否，则返回步骤s121，对所述对流动态计算域的下一边界单元进行判断；
s123：判断该边界单元是否满足：
其中，q取对流动态计算域边界单元的格心指向边界单元紧邻单元格心的单位向量，θd表示上游单元容差角，取5°≤θd≤10°；若是，则执行步骤s124；若否，则返回步骤s121，对所述对流动态计算域的下一边界单元进行判断；
s124：对于超声速无粘流动，无需再判断，直接执行步骤s13；
对于亚声速流动和粘性流动，判断该边界单元是否满足：
||δw||+||δ(δw)||＜εd(9)
其中，||δ(δw)||表示边界单元守恒量更新量受边界单元在对流动态计算域中紧邻单元影响的模值；δt表示迭代步长，ccfl表示时间推进格式的cfl数，δri表示对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响；
对于亚声速无粘单元，对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响δri表示为：
其中，δfc表示对流通量的变化量；下标i+1表示对流动态计算域的边界单元沿正i方向的紧邻单元，下标i-1表示对流动态计算域的边界单元沿负i方向的紧邻单元，δsi+1/2表示对流动态计算域的边界单元沿正i方向单元面的面积；δsi-1/2表示对流动态计算域的边界单元沿负i方向单元面的面积；
对于粘性单元，对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响δri表示为：
若是，则执行步骤s13；若否，则返回步骤s121，对所述对流动态计算域的下一边界单元进行判断。
在一种可能的实现方式中，在本发明提供的上述飞行器动态气动特性模拟的迭代推进扰动域更新方法中，步骤s15，对所述粘性动态计算域的所有边界单元逐个判断是否受粘性效应主导，具体包括：
对所述粘性动态计算域的所有边界单元逐个判断是否满足：
其中，φ表示粘性新增阈值的缩比因子，||δw||max表示对流动态计算域所有单元的守恒量更新量在当前迭代步的最大值，(||δw||max)min表示当前时刻中内迭代的第1迭代步至当前迭代步中单步守恒量更新量最大值的最小值。
在一种可能的实现方式中，在本发明提供的上述飞行器动态气动特性模拟的迭代推进扰动域更新方法中，步骤s18中，对所述非定常动态计算域的所有边界单元逐个判断是否不再受非定常效应主导，具体包括：
对所述非定常动态计算域的所有边界单元逐个判断是否满足：
本发明提供的上述飞行器动态气动特性模拟的迭代推进扰动域更新方法，残差估计仅在对流动态计算域中计算残差的无粘项，仅在粘性动态计算域中计算残差的粘性项；时间积分也仅在对流动态计算域内执行，因此，本发明提供的上述迭代推进扰动域更新方法通过采用非定常、对流、粘性三类动态计算域，可高效避免现有方法在飞行器流场网格的预设计算域中全局更新所有网格单元所导致的无效计算，从而可有效提升飞行器动态气动特性数值模拟的计算效率，解决数值模拟计算效率难以满足飞行器气动设计日益增长的实际应用需求的问题。将所需存储的数据分为所有单元的固有信息和与求解更新相关的信息两类；其中，对于第一类数据，以静态数据结构存储预设计算域中所有网格单元的信息；对于第二类数据，以动态数据结构仅存储对流动态计算域中网格单元的信息，在每个迭代步和每个物理时刻及时根据非定常动态域调整第二类数据的存储空间，因此，本发明提供的上述迭代推进扰动域更新方法可有效降低飞行器动态气动特性数值模拟的内存需求。根据扰动传播特性，仅分别将受扰单元和粘性效应主导单元加入对流动态计算域与粘性动态计算域中；利用求解收敛特征，在不影响收敛速率的前提下及时将已收敛单元从求解更新中移除；通过衡量非定常效应的影响，及时将定常单元从求解更新中移除；因此，本发明提供的上述迭代推进扰动域更新方法可始终保持仅对求解未收敛的非定常单元求解，仅在粘性效应主导单元中考虑粘性，有效避免了现有方法在飞行器流场网格的预设计算域中全局更新所有网格单元所导致的无效计算，从而解决了数值模拟计算效率难以满足飞行器气动设计日益增长的实际应用需求的问题。
附图说明
图1为本发明实施例1提供的迭代推进扰动域更新方法的流程图；
图2为利用本发明实施例1提供的迭代推进扰动域更新方法求解跨声速naca0012翼型振荡问题的流场与非定常、对流两类动态计算域的演化过程图；
图3为利用本发明实施例1提供的迭代推进扰动域更新方法求解跨声速naca0012翼型振荡问题的非定常、对流两类动态计算域网格量变化曲线与收敛曲线；
图4为利用本发明实施例1提供的迭代推进扰动域更新方法求解超声速有翼导弹俯仰问题的流场与非定常、对流两类动态计算域的演化过程图；
图5为利用本发明实施例1提供的迭代推进扰动域更新方法求解超声速有翼导弹俯仰问题的非定常、对流两类动态计算域网格量变化曲线与收敛曲线。
具体实施方式
下面将结合本发明实施方式中的附图，对本发明实施方式中的技术方案进行清楚、完整的描述，显然，所描述的实施方式仅仅是作为例示，并非用于限制本发明。
飞行器动态气动特性模拟的迭代推进求解过程具有如下特征：(1)流场的扰动始于控制方程无法满足之处；(2)携带着间断信息的扰动会随时间以有限速度逐渐传向周围流场，尚未受到扰动的区域仍保持初始状态；(3)可压缩流动中，上游流场的收敛不会晚于下游流场；(4)无论是否分离，粘性效应仅在有限区域起主导作用；(5)非定常效应主导区域随流场演化逐渐变化。为充分利用上述5点特征，避免现有数值模拟方法中的无效计算，本发明提供了一种飞行器动态气动特性模拟的迭代推进扰动域更新方法，通过仅对求解未收敛的非定常单元求解，仅在粘性效应主导单元中考虑粘性，有效减少内迭代中单个迭代步的计算量和达到收敛态所需的迭代步数，从而提高飞行器动态气动特性数值模拟的效率。
本发明提供的一种飞行器动态气动特性模拟的迭代推进扰动域更新方法，包括如下步骤：
s1：读入数据，包括飞行器流场的网格、预设计算域、边界条件和计算设置；
s2：根据来流条件或根据给定流场，对预设计算域中的流场进行初始化；
s3：根据流场初始化方式，建立非定常动态计算域、对流动态计算域和粘性动态计算域；
s4：将求解所需存储的数据分为两类；第一类数据是单元的固有信息，采用静态数据结构存储预设计算域中所有网格单元的信息，包括网格坐标以及当前时刻和上两个时刻的流场变量；第二类数据是与求解更新相关的信息，采用动态数据结构仅存储对流动态计算域内网格单元的信息，包括守恒量更新量和当地迭代步长；
s5：根据边界条件类型，为边界虚网格的守恒量赋值；
s6：将流动控制方程的残差项分为无粘项与粘性项两类；其中，无粘项包括对流通量和迭代推进法的源项，无粘项在对流动态计算域中计算；粘性项包括粘性通量和湍流模型方程源项，粘性项在粘性动态计算域中计算；
s7：在对流动态计算域中，求解守恒量更新量，更新当前时刻的流场变量；
s8：对对流动态计算域的所有边界单元逐个判断是否已受到无粘扰动；若是，则执行步骤s9；若否，则返回步骤s8，对对流动态计算域的下一边界单元进行判断；遍历对流动态计算域的所有边界单元后，执行步骤s10；
s9：衡量无粘扰动的传播方向，并将该边界单元的紧邻单元中位于传播方向上的单元分别纳入对流动态计算域和非定常动态计算域，返回步骤s8，对对流动态计算域的下一边界单元进行判断；
s10：判断非定常动态计算域是否有所增大；若是，则执行步骤s11和步骤s12；若否，则执行步骤s12；
s11：根据非定常动态计算域，重新分配第二类数据的存储空间；
s12：对对流动态计算域的所有边界单元逐个判断是否同时满足已收敛、位于可压缩流动中、位于最上游以及不再受对流动态计算域中其他单元的影响四个条件；若是，则执行步骤s13；若否，则返回步骤s12，对对流动态计算域的下一边界单元进行判断；遍历对流动态计算域的所有边界单元后，执行步骤s15；
s13：将该边界单元从对流动态计算域中移除，并判断该边界单元是否存在于粘性动态计算域中；若是，则执行步骤s14；若否，则返回步骤s12，对对流动态计算域的下一边界单元进行判断；
s14：将该边界单元从粘性动态计算域中移除，返回步骤s12，对对流动态计算域的下一边界单元进行判断；
s15：对粘性动态计算域的所有边界单元逐个判断是否受粘性效应主导；若是，则执行步骤s16；若否，则返回步骤s15，对粘性动态计算域的下一边界单元进行判断；遍历粘性动态计算域的所有边界单元后，执行步骤s17；
s16：将该边界单元的紧邻单元中位于对流动态计算域中的单元纳入粘性动态计算域，返回步骤s15，对粘性动态计算域的下一边界单元进行判断；
s17：判断内迭代是否达到收敛条件；若是，则执行步骤s18；若否，则返回步骤s5，进入内迭代下一迭代步的计算，重复执行步骤s5～s17；
s18：利用物理时间导数衡量非定常效应的影响，对非定常动态计算域的所有边界单元逐个判断是否不再受非定常效应主导；若是，则执行步骤s19；若否，则返回步骤s18，对非定常动态计算域的下一边界单元进行判断；遍历非定常动态计算域的所有边界单元后，执行步骤s20；
s19：将该边界单元从非定常动态计算域中移除，返回步骤s18，对非定常动态计算域的下一边界单元进行判断；
s20：判断非定常动态计算域是否有所缩小；若是，则执行步骤s21和步骤s22；若否，则执行步骤s22；
s21：根据非定常动态计算域，重新分配第二类数据的存储空间；
s22：更新上两个时刻的流场变量，将对流动态计算域重置为非定常动态计算域的范围，将粘性动态计算域重置为对流动态计算域中粘性效应起主导作用的区域；
s23：判断是否达到指定的末态物理时刻；若是，则执行步骤s24；若否，则返回步骤s5，进入下一物理时刻的计算，重复执行步骤s5～s23；
s24：输出结果。
下面通过一个具体的实施例对本发明提供的上述迭代推进扰动域更新方法的具体实施进行详细说明。
实施例1：如图1所示，本发明提供的上述迭代推进扰动域更新方法具体包括如下步骤：
第一步：数据读入。
读入数据包括飞行器流场的网格、预设计算域、边界条件和计算设置。
第二步：流场初始化。
流场初始化分为根据来流条件初始化和根据给定流场初始化两种方式。当根据来流条件初始化时，将预设计算域中所有网格单元的守恒量赋为来流值；当根据给定流场初始化时，将预设计算域中所有网格单元的守恒量赋为给定流场值。
第三步：建立非定常、对流、粘性三类动态计算域。
(1)当根据来流条件初始化时，非定常、对流、粘性三类动态计算域均根据壁面边界建立。携带着流场变化信息的扰动从物体表面产生，受粘性主导的单元也为紧邻物体表面，因此，非定常、对流两类动态计算域取壁面的若干层相邻单元作为初始单元，例如可以取5～10层壁面相邻网格单元，而粘性动态计算域仅取紧邻壁面的1层单元作为初始单元。
(2)当根据给定流场初始化时，非定常、对流、粘性三类动态计算域根据给定流场的流动特性建立。非定常、对流两类动态计算域的初始单元均为给定流场中流动特性与来流条件不符的已受扰单元，粘性动态计算域的初始单元则为对流动态计算域中的粘性效应主导单元。
对流动态计算域的初始单元为给定流场中的已受扰单元，这些单元的流动特性与来流条件不符，满足：
||w-w∞||/||δw(1)||max＞εa,c(1)
其中，w∞表示来流条件的守恒量；||δw(1)||max表示对流动态计算域所有单元的守恒量更新量在第1物理时刻的第1迭代步的最大值；εa,c表示对流新增阈值，取10-4≤εa,c≤10-6。
非定常动态计算域的初始单元与对流动态计算域的初始单元一致。
粘性动态计算域的初始单元为对流动态计算域中的粘性效应主导单元。定义粘性效应衡量参数ψ，令其表示粘性扰动与无粘扰动的质量流量之比。粘性效应主导单元的粘性效应衡量参数应处于较大量级，故粘性动态计算域中的单元应满足：
其中，ψ表示粘性效应衡量参数，为粘性扰动与无粘扰动的质量流量之比；
第四步：分配存储空间。
将求解所需存储的数据分为两类；第一类数据是单元的固有信息，采用静态数据结构存储预设计算域中所有网格单元的信息，包括网格坐标以及当前时刻和上两个时刻的流场变量；第二类数据是与求解更新相关的信息，采用动态数据结构仅存储对流动态计算域内网格单元的信息，包括守恒量更新量和当地迭代步长。
对于第一类数据，采用静态数据结构存储预设计算域中所有网格单元的信息。例如，令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表示粘性通量，qt表示湍流模型方程源项，
第七步：时间积分。
在对流动态计算域中，采用时间推进格式计算式(4)的左端项；所采用的时间积分方法与现有方法一致。
现有方法中，最为耗时的残差估计、时间积分两步需对静态的预设计算域内所有网格单元执行，约占99％的总求解时间。相比之下，本发明提供的上述迭代推进扰动域更新方法，第六步的残差估计仅在对流动态计算域内计算残差的无粘项，仅在粘性动态计算域内计算残差的粘性项；第七步的时间积分也仅在对流动态计算域内执行。因此，本发明提供的上述迭代推进扰动域更新方法可有效提升飞行器动态气动特性数值模拟的计算效率。
第八步：增大非定常、对流两类动态计算域。
对对流动态计算域的所有边界单元逐个判断是否已受到无粘扰动，即对对流动态计算域的所有边界单元逐个判断是否满足：
||δw||/||δw(1)||max＞εa,c(5)
其中，||δw||表示对流动态计算域的边界单元在当前迭代步中守恒量更新量的模值；若对流动态计算域的边界单元满足式(5)，则衡量无粘扰动的传播方向，并将该边界单元的紧邻单元中位于传播方向上的单元纳入对流动态计算域和非定常动态计算域，返回第八步，对对流动态计算域的下一边界单元进行判断；若对流动态计算域的边界单元不满足式(5)，则返回第八步，对对流动态计算域的下一边界单元进行判断；遍历对流动态计算域的所有边界单元后，执行第九步。
无粘扰动相对于流动以声速传播，沿无粘扰动的传播方向，无粘扰动的传播速度为正；令q表示单位方向向量，则无粘扰动沿q方向传播表示为：
u·q+a＞0(6)
其中，u表示流动速度矢量，a表示声速；q取对流动态计算域边界单元的格心指向边界单元某一格点的单位向量，若q满足式(6)，则将与对流动态计算域边界单元共享该格点的紧邻单元纳入对流动态计算域和非定常动态计算域。非定常动态计算域应始终包含对流动态计算域。
第九步：重新分配第二类数据的存储空间。
判断非定常动态计算域是否有所增大；若是，则根据非定常动态计算域，重新分配第二类数据的存储空间后，执行第十步；若否，则执行第十步。
第十步：缩小对流、粘性两类动态计算域。
为了在不影响收敛速率的前提下，尽可能地减少单个迭代步的计算量，可以缩小对流动态计算域。对对流动态计算域的所有边界单元逐个判断以下条件实现对流动态计算域的缩小：
(1)对对流动态计算域的所有边界单元逐个判断是否已收敛，即对对流动态计算域的所有边界单元逐个判断是否满足：
||δw||/||w(1)||max＜εd(7)
其中，εd表示删除阈值，取εd≤10-7；若该边界单元满足式(7)，则继续判断(2)；若否，则对对流动态计算域的下一边界单元进行判断；
(2)判断该边界单元是否位于可压缩流动中，即判断该边界单元的马赫数是否大于0.3；若是，则继续判断(3)；若否，则返回判断(1)，对对流动态计算域的下一边界单元进行判断；
(3)判断该边界单元是否位于最上游，即判断该边界单元是否满足：
其中，q取对流动态计算域边界单元的格心指向边界单元紧邻单元格心的单位向量，θd表示上游单元容差角，取5°≤θd≤10°；若该边界单元满足式(8)，则继续判断(4)；若否，则返回判断(1)，对对流动态计算域的下一边界单元进行判断；
(4)判断该边界单元是否不再受对流动态计算域中其他单元的影响，对于超声速无粘流动，则无需判断、自动满足；对于亚声速流动和粘性流动，判断该边界单元是否满足：
||δw||+||δ(δw)||＜εd(9)
其中，||δ(δw)||表示边界单元守恒量更新量受边界单元在对流动态计算域中紧邻单元影响的模值；δt表示迭代步长，ccfl表示时间推进格式的cfl数，δri表示对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响；对于亚声速无粘单元，对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响δri表示为：
其中，δfc表示对流通量的变化量；下标i+1表示对流动态计算域的边界单元沿正i方向的紧邻单元，下标i-1表示对流动态计算域的边界单元沿负i方向的紧邻单元，δsi+1/2表示对流动态计算域的边界单元沿正i方向单元面的面积；δsi-1/2表示对流动态计算域的边界单元沿负i方向单元面的面积；对于粘性单元，对流动态计算域中紧邻单元对边界单元的残差项沿i方向的影响δri表示为：
若该边界单元满足式(9)和式(10)，则该边界单元已同时满足已收敛、位于可压缩流动中、位于最上游以及不再受所述对流动态计算域中其他单元的影响四个条件，则将该边界单元从对流动态计算域中移除，若该边界单元同时存在于粘性动态计算域中，为保证对流动态计算域始终包含粘性动态计算域，则将该边界单元从粘性动态计算域中移除，返回第十步，对对流动态计算域的下一边界单元进行判断；若该边界单元不满足式(9)和式(10)，则返回第十步，对对流动态计算域的下一边界单元进行判断。遍历对流动态计算域的所有边界单元后，执行第十一步。
第十一步：增大粘性动态计算域。
对粘性动态计算域的所有边界单元逐个判断是否受粘性效应主导，即对粘性动态计算域的所有边界单元逐个判断是否满足：
其中，φ表示粘性新增阈值的缩比因子，||δw||max表示对流动态计算域所有单元的守恒量更新量在当前迭代步的最大值，(||δw||max)min表示第1迭代步至当前迭代步中单步守恒量更新量最大值的最小值；若粘性动态计算域的边界单元满足式(13)和(14)，则将该边界单元的紧邻单元中位于对流动态计算域中的单元纳入粘性动态计算域，返回第十一步，对粘性动态计算域的下一边界单元进行判断；若粘性动态计算域的边界单元不满足式(13)和(14)，则返回第十一步，对粘性动态计算域的下一边界单元进行判断；遍历粘性动态计算域的所有边界单元后，执行第十二步。流动控制方程的粘性项呈椭圆型，因此，可以在不增大对流动态计算域的前提下，将已受粘性扰动单元的所有紧邻单元均加入粘性动态计算域。
第十二步：判断内迭代是否达到收敛条件；若是，则执行第十三步；若否，则返回第五步，进入内迭代下一迭代步的计算，重复执行第五步至第十二步。
第十三步：缩小非定常动态计算域。
利用物理时间导数衡量非定常效应的影响，对非定常动态计算域的所有边界单元逐个判断是否不再受非定常效应主导，即对非定常动态计算域的所有边界单元逐个判断是否满足：
若非定常动态计算域的边界单元满足式(15)，则将该边界单元将从非定常动态计算域中移除，返回第十三步，对非定常动态计算域的下一边界单元进行判断；若非定常动态计算域的边界单元不满足式(15)，则返回第十三步，对非定常动态计算域的下一边界单元进行判断；遍历非定常动态计算域的所有边界单元后，执行第十四步。
本发明在第八步、第十一步中，根据扰动传播特性，仅分别将受扰单元和粘性效应主导单元加入对流动态计算域与粘性动态计算域中；在第十步中，利用求解收敛特征，在不影响收敛速率的前提下及时将已收敛单元从求解更新中移除；在第十三步中，通过衡量非定常效应的影响，及时将定常单元从求解更新中移除；因此，本发明提供的上述迭代推进扰动域更新方法可始终保持仅对求解未收敛的非定常单元求解，仅在粘性效应主导单元中考虑粘性，可有效避免现有非定常数值模拟方法求解中的无效计算。
第十四步：重新分配第二类数据的存储空间。
判断非定常动态计算域是否有所缩小；若是，则根据非定常动态计算域，重新分配第二类数据的存储空间后，执行第十五步；若否，则执行第十五步。
第十五步：重置对流、粘性两类动态计算域。
更新上两个时刻的流场变量；将对流动态计算域重置为非定常动态计算域的范围，将粘性动态计算域重置为对流动态计算域中的粘性效应主导区域。
第十六步：判断是否达到指定的末态物理时刻；若是，则执行第十七步；若否，则返回第五步，进入下一物理时刻的计算，重复执行第五步至第十六步。
第十七步：输出结果。
下面分别对马赫数0.755的naca0012翼型俯仰振荡问题和马赫数1.58的有翼导弹俯仰运动问题进行模拟。
(1)对马赫数0.755的naca0012翼型俯仰振荡问题进行模拟。
图2示意了本发明提供的上述迭代推进扰动域更新方法在不同物理时刻所对应的非定常动态计算域与内迭代中的最小对流动态计算域。图2中，ηu表示非定常动态计算域与预设计算域的网格量之比，ηc,min表示内迭代最小对流动态计算域与预设计算域的网格量之比。流场以来流条件初始化，扰动从壁面边界产生，非定常动态计算域的初始单元取壁面的10层相邻单元。在第1个物理时刻，非定常动态计算域便从壁面扩展至整个预设计算域，并在此后一直保持在此范围。而在所有时刻的内迭代中，对流动态计算域都根据非定常动态计算域重置，并随着求解收敛逐渐从远场向壁面、从上游向下游收缩，最后仅剩近壁区及其下游区域。
图3给出了本发明提供的上述迭代推进扰动域更新方法在最后6个物理时刻的非定常、对流两类动态计算域网格量变化曲线，并与现有全局更新方法的收敛曲线进行对比。图3展示的对流动态计算域网格量变化曲线表明，在一次内迭代中，对流动态计算域会单调缩小，ηc从1.0降至0.44，即在内迭代结束时，仅有44％的预设计算域单元参与了求解。图2中，现有全局更新方法与本发明提供的上述迭代推进扰动域更新方法的收敛曲线的对比表明，本发明提供的上述迭代推进扰动域更新方法不仅不会对非定常数值模拟的收敛速率产生不利影响，而且可减少内迭代达到收敛态所需的求解步数。相比于现有方法，本发明提供的上述迭代推进扰动域更新方法可节省26.5％的时间。
(2)对马赫数1.58的有翼导弹俯仰运动问题进行模拟。
图4给出了本发明提供的上述迭代推进扰动域更新方法在不同物理时刻所对应的非定常动态计算域与内迭代中的最小对流动态计算域。对于超声速流动，非定常动态计算域不再如亚声速情况一般始终保持在整个预设计算域范围，而是随着导弹的俯仰运动也呈现出周期性变化。第1个物理时刻开始时，非定常动态计算域依据壁面建立，仅将紧邻壁面边界的10层单元作为初始单元；而在第1个物理时刻结束时，非定常动态计算域的ηu已达0.999。在前半个周期，非定常动态计算域一直保持着较大范围；之后，非定常动态计算域则每半个周期重复一次缩小与增大。当导弹回复至初始状态时，非定常动态计算域便会开始缩小。非定常动态计算域的缩小过程先急后缓，历经小半个周期，越接近迎角最值则缩小得越慢。相比之下，从最小范围增至最大范围的增大过程则十分迅速，仅在一个时刻内完成。
图5给出了本发明提供的上述迭代推进扰动域更新方法在最后若干个物理时刻的非定常、对流两类动态计算域网格量变化曲线，并与现有全局更新方法的收敛曲线进行对比。对比结果表明，本发明提供的上述迭代推进扰动域更新方法与现有的全局更新方法一样，在每个内迭代步都实现了求解收敛。对于超声速三维问题，本发明的加速效果愈发明显；当迭代更新扰动域更新方法已结束求解时，现有方法仍有16个时刻尚待计算。得益于缩小非定常、对流两类动态计算域所获得的计算量节省与迭代步数节省，本发明提供的上述迭代推进扰动域更新方法相比于现有方法可获得38.9％的时间节省。
本发明提供的上述飞行器动态气动特性模拟的迭代推进扰动域更新方法，残差估计仅在对流动态计算域中计算残差的无粘项，仅在粘性动态计算域中计算残差的粘性项；时间积分也仅在对流动态计算域内执行，因此，本发明提供的上述迭代推进扰动域更新方法通过采用非定常、对流、粘性三类动态计算域，可高效避免现有方法在飞行器流场网格的预设计算域中全局更新所有网格单元所导致的无效计算，从而可有效提升飞行器动态气动特性数值模拟的计算效率，解决数值模拟计算效率难以满足飞行器气动设计日益增长的实际应用需求的问题。将所需存储的数据分为所有单元的固有信息和与求解更新相关的信息两类；其中，对于第一类数据，以静态数据结构存储预设计算域中所有网格单元的信息；对于第二类数据，以动态数据结构仅存储对流动态计算域中网格单元的信息，在每个迭代步和每个物理时刻及时根据非定常动态域调整第二类数据的存储空间，因此，本发明提供的上述迭代推进扰动域更新方法可有效降低飞行器动态气动特性数值模拟的内存需求。根据扰动传播特性，仅分别将受扰单元和粘性效应主导单元加入对流动态计算域与粘性动态计算域中；利用求解收敛特征，在不影响收敛速率的前提下及时将已收敛单元从求解更新中移除；通过衡量非定常效应的影响，及时将定常单元从求解更新中移除；因此，本发明提供的上述迭代推进扰动域更新方法可始终保持仅对求解未收敛的非定常单元求解，仅在粘性效应主导单元中考虑粘性，有效避免了现有方法在飞行器流场网格的预设计算域中全局更新所有网格单元所导致的无效计算，从而解决了数值模拟计算效率难以满足飞行器气动设计日益增长的实际应用需求的问题。
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。
- 还没有人留言评论。精彩留言会获得点赞！