CN108563843B - 定常可压缩流动的扰动区域更新方法 - Google Patents

Abstract

本发明公开的定常可压缩流动的扰动区域更新方法，属于计算流体力学技术领域。所述方法包括输入文件载入、流场初始化、动态计算域建立、边界虚网格更新、动态计算域更新、残差估计、时间积分、判断是否收敛、结果输出的步骤。本发明适用于任何基于时间推进方法的可压缩流动的求解，可与任意空间离散格式、时间推进格式和现有的其他时间推进加速技术相结合，并适用于不同的初始化方式以及复杂的、不连续的壁面情况。本发明能够大幅提高计算效率；对于亚声速流动问题，本发明能够在求解接近定常状态的后期合理地去除冗余的远场单元，有效地缩小计算域；对于超声速流动问题，本发明还能够降低求解过程的存储需求。

Description

定常可压缩流动的扰动区域更新方法

技术领域

本发明属于计算流体力学技术领域，具体涉及一种定常可压缩流动的扰动区域更新方法，涉及数值模拟收敛速率加速技术的原理与实现。

背景技术

含有时间导数项的Euler/N-S方程的数学特性对亚声速、超声速流动都始终保持双曲型，因此，当前计算流体力学中普遍采用时间推进法求解定常可压缩流动问题。时间推进法以任意条件为起始，沿时间对非定常控制方程进行积分，并以时间趋近无穷大时获得的流场为定常状态的流场。

然而，基于时间推进法的流场求解往往极其耗时，高效的收敛加速技术一直是计算流体力学的研究热点。时间推进法的求解过程可主要分为输入文件载入、流场初始化、边界虚网格更新、残差估计、时间积分、结果输出等步骤。其中，残差估计与时间积分两步由于需在每个迭代步中对所有网格单元进行操作，因而成为整个求解中最为耗时的部分。

残差估计与时间积分两步的计算量与网格量呈正比，因此，降低所述计算量最直观的方法就是降低网格量。目前，这一策略只能通过网格生成来实现，如合理地设计网格的拓扑与分布或自适应加密技术。虽然网格生成技术已经能够静态或动态地控制网格的分布与变形，但仍无法在计算中灵活地改变计算区域的范围。

现有的求解方法均采用在整个流场求解过程中始终不变的静态计算域。这种静态计算域的求解方式具有以下两点弊端：

1)增加无效计算。对于未受到扰动的网格单元和求解已达收敛条件的网格单元，其流动参数的更新对整个流场的求解并没有实际的有利贡献。

2)引入数值误差。对流场特性不应改变的未受扰动网格单元进行反复的迭代会在求解中引入更多的舍入误差或截断误差。

发明内容

从本质上来说，定常流动的时间推进求解过程实际上就是粗略地模拟从给定的初始流场发展至定常状态的时变过程。时间推进求解过程具有以下四点特征：1)定常控制方程无法满足之处会产生扰动；2)这些携带着间断信息的扰动会随着时间推进逐渐传向周围流动；3)未受这些扰动影响的区域总是保持着初始的状态；4)上游流动总是先于下游流动收敛。

因此，为了提高定常可压缩流动的数值模拟效率，本发明利用上述的时间推进法的四点特征，提出了一种定常可压缩流动的扰动区域更新方法。该方法合理利用了时间推进求解过程中，扰动以有限速度传播的特性；通过衡量扰动波的传播速度和方向，能够在流场中高效地甄别出有必要更新的区域；并采用动态数据结构存储，实现仅对受扰区域中未收敛网格单元进行更新的策略，从而达到了降低计算量的目的。

本发明提出的一种定常可压缩流动的扰动区域更新方法是基于时间推进法的特征而提出的：一方面，在时间推进过程中，扰动始于流动特性不满足于定常控制方程之处，并以有限波速向周围流动传播；另一方面，对流速度的影响使得上游流动总是先于下游流动收敛。因此，本发明提供的扰动区域更新方法使用动态计算域，仅对求解中受到扰动影响且尚未收敛的网格单元进行更新。其求解过程中，计算域随时间推进变化并遵循如下规则：首先，从扰动产生处建立；其次，随扰动的传播向给定计算域扩张；最后，随着求解收敛逐渐收缩并向下游移动。

由于本发明提供的扰动区域更新方法的动态计算域是未知且时变的，为了实现动态计算域的高效存储，本发明提供的一种定常可压缩流动的扰动区域更新方法有别于以往的计算流体力学程序，同时采用了静态数组与动态数据结构。扰动区域更新方法的数据结构具有如下特征：

1.所有网格单元都需存储的变量仍旧采用静态数组存储；所述变量如网格单元坐标、网格单元格心守恒变量。

2.与流场更新相关的变量采用双向环形链表存储；所述变量如残差、之前时间步的解。

3.双向环形链表中的每一项对应于网格的一行；

4.双向环形链表中各项按照网格的行标号从小到大的顺序排列，相邻项间的行标号可以不连续；

5.双向、环形的链表设计有助于提高动态计算域更新的效率。

本发明提供的一种定常可压缩流动的扰动区域更新方法的步骤如下：

第一步：输入文件载入；

第二步：流场初始化；

流场初始化包括根据远场条件初始化和根据给定流场初始化两种方式。

第三步：动态计算域建立；

对于根据远场条件初始化的方式，动态计算域将根据壁面条件建立。对于根据给定流场初始化的方式，动态计算域的建立需要衡量给定流场的残差。这种情况中，动态计算域的建立可分为两个子步骤：

(3.1)筛选流动特性与来流不同的网格单元，建立初始动态计算域；

(3.2)依据残差细化步骤(3.1)所建立的初始动态计算域。

第四步：边界虚网格更新；

根据边界条件，为边界虚网格的守恒量赋值。

第五步：动态计算域更新；

动态计算域的更新过程可概括为以下四个子步骤，这些操作必须紧密有序地完成才能确保扰动区域更新方法的正确与效率：

(5.1)更新现有行；

这一步要遍历所有的现有行，包括新增网格单元与删除网格单元两类操作。通过衡量动态计算域边界处网格单元扰动波的传播方向与大小，判断是否需要扩宽当前动态计算域的现有行。若扰动能够传出，则扩大现有行的范围；若不能传出，则通过衡量收敛性与网格单元间的相互影响，判断是否可以缩小现有行的范围。

(5.2)新增行；

若动态计算域没有包括所有预设网格中的行，通过衡量动态计算域边界处网格单元扰动波的传播方向和大小，判断是否需要在动态计算域中新增行。若扰动能够传入尚不包含在当前动态计算域内的行，则在链表中新增项用于存储新增行的信息。

(5.3)移除链表冗余项；

若链表中出现不包含网格单元的空项，将该空项从链表中移除。

(5.4)再分配存储空间。

若动态计算域有所改变，则重新分配链表的存储空间。

第六步：残差估计；

估计动态计算域内每个网格单元的对流通量、扩散通量和源项。

第七步：时间积分；

使用显式、隐式时间积分方法确定动态计算域内网格单元的守恒量修正量。

第八步：判断是否收敛；

判断动态计算域内所有网格单元的守恒量修正量的模值是否均小于所要求的收敛阈值。若不满足，则返回第四步；若满足，则退出迭代计算，继续下一步。

第九步：结果输出。

本发明的优点在于：

1.本发明提出的扰动区域更新方法适用于任何基于时间推进方法的可压缩流动的求解，可与任意空间离散格式、时间推进格式和现有的其他时间推进加速技术相结合，并适用于不同的初始化方式以及复杂的、不连续的壁面情况。

2.本发明提出的扰动区域更新方法能够大幅提高计算效率。计算效率的提高得益于两方面：一方面，该方法只计算已受扰动和未收敛的网格单元，从而大幅降低了单个迭代步的求解量；另一方面，该方法合理地去除了上游已收敛的网格单元，避免其在定常状态附近产生小幅振荡，使得求解更容易收敛至定常状态，从而减小了达到相同收敛阈值的总迭代步数。

3.对于超声速情况，本发明提出的扰动区域更新方法只计算激波附近及波后区域，因此能够降低求解过程的存储需求。

4.对于亚声速情况，本发明提出的扰动区域更新方法在求解接近定常状态的后期能够合理地去除冗余的远场网格单元，从而有效地缩小计算域。

附图说明

图1为本发明的扰动区域更新方法的流程图；

图2为马赫数0.85、1°攻角条件下，NACA 0012翼型采用扰动区域更新方法的求解过程以及定常解示意图；

图3为NACA 0012翼型算例的动态计算域网格量随迭代数的变化，以及扰动区域更新方法和全局更新方法最大残差随迭代数的变化；

图4为马赫数2.5时，超声速圆柱绕流采用扰动区域更新方法的求解过程，以及其定常解与全局更新方法、激波装配法结果的对比；

图5为超声速圆柱绕流算例的动态计算域网格量随迭代数的变化，以及扰动区域更新方法和全局更新方法最大残差随迭代数的变化。

具体实施方式

下面结合附图和实施例对本发明做进一步说明。

本发明提供一种定常可压缩流动的扰动区域更新方法，如图1所示流程，具体包括如下步骤：

第一步：输入文件载入；

载入的文件包括网格信息文件、边界条件文件、计算条件设置文件。

第二步：流场初始化；

流场初始化包括根据远场条件初始化和根据给定流场初始化两种方式。初始化的方式将决定动态计算域的建立方式。

第三步：动态计算域建立；

对于根据远场条件初始化的方式，动态计算域将根据壁面条件建立。以根据远场条件初始化的时间推进过程相当于在均匀流动中突然放置一个物体后周围流动变化的非定常过程。因此，流场中的变化与扰动始于壁面。这种情况中，动态计算域的建立相对简单，动态计算域的初始网格单元就是所有紧邻壁面边界的那一层网格单元。

对于根据给定流场初始化的方式，动态计算域的建立需要衡量给定流场中网格单元的残差。这种情况中，动态计算域的建立可分为两个子步骤：

(3.1)筛选流动特性与来流不同的网格单元，建立初始动态计算域；

在给定流场中，若一个网格单元的流动特性与来流相同，则该网格单元不是未受扰动就是已经收敛，而这类网格单元一定不会包含于动态计算域中。为了降低残差判断的计算量，利用这个特征，先建立只包含流动特性与来流条件不同的网格单元的初始动态计算域；其中，流动特性与来流条件不同的判据可表示为：

||W-W_∞||＞ε_d (1)

式中，W表示求解的网格单元的守恒量，W_∞表示来流条件对应的守恒量，ε_d表示收敛阈值。

(3.2)依据残差细化步骤(3.1)所建立的初始动态计算域。

离散形式的非定常流动控制方程可表示为：

式中，ΔW表示网格单元每一个时间步的守恒量修正量，R表示网格单元的残差，Δt表示网格单元的当地时间步长，Ω表示网格单元的体积。

虽然本发明流场时间推进迭代求解的收敛条件是守恒量修正量ΔW小于收敛阈值ε_d，但在步骤(3.2)中本发明采用残差判断网格单元是否参与后续的流场时间推进迭代求解。其一，由于残差R的确定不涉及时间推进过程，因此确定残差R的计算量远小于确定守恒量修正量ΔW的计算量；其二，残差R为守恒量修正量ΔW中的主导项，足以表示守恒量修正量ΔW的量级。因此本发明根据初始动态计算域中所有网格单元的残差筛选需要参与后续流场时间推进迭代求解的网格单元。为了考虑守恒量修正量ΔW与残差R之间的差异，判断中，残差R乘以一个小于1.0的缩放因子以避免未收敛网格单元的遗漏。因此，流场时间推进迭代求解尚未收敛网格单元，即仍需保留在动态计算域中的网格单元，应满足如下条件：

式中，

为残差的缩放因子。

其余不满足条件(3)的网格单元从步骤(3.1)建立的初始动态计算域中去除，得到参与后续流场求解的动态计算域。

第四步：边界虚网格更新；

根据第一步中所述的边界条件，为边界虚网格的守恒量赋值。

第五步：动态计算域更新；

动态计算域的更新过程可概括为以下四个子步骤：

(5.1)更新现有行；

这一步要遍历所有的现有行，包括新增网格单元与删除网格单元两类操作。通过衡量扰动波沿网格方向的传播方向与大小，判断扰动波是否会传出当前动态计算域的现有行。若能够传出，则沿网格行方向扩大该现有行；若不能，则通过衡量收敛性与网格单元间的相互影响，判断是否可以缩小现有行的范围。

(5.1.1)新增网格单元操作；

新增网格单元操作是依据扰动波的方向与大小，将受到扰动的网格单元加入动态计算域，主要在求解前期发挥作用。定义I、J表示网格和动态计算域的方向，I表示行方向，J表示列方向，则(I,J)对应动态计算域中的任意一个网格单元，记为网格单元(I,J)。

以I正向的情况为例，I负向采用同样的方法。为了便于表述，假设右向为网格I标号的增大方向，假设网格单元(I,J)已在之前的时间步被纳入动态计算域但其右侧网格单元还在动态计算域之外，若要确定网格单元(I,J)右侧网格单元是否需要加入动态计算域需要依次经过以下三个子步骤的判断：

(5.1.1.1)判断网格单元(I,J)的流动特性是否已经受到了扰动的影响；

(5.1.1.2)确定影响网格单元(I,J)的扰动是否能够向右传播；

(5.1.1.3)估计在下一个时间步n+1能够受到影响的右侧网格的数目。

其中，若子步骤(5.1.1.1)或(5.1.1.2)为否定结果，则说明右侧网格单元在当前时间步n无需加入动态计算域，余下子步骤也无需执行。

子步骤(5.1.1.1)可以通过衡量网格单元(I,J)在上一个时间步的守恒量修正量来完成。如果网格单元(I,J)已经受到了扰动的影响，则它的守恒量修正量应在一个不可忽视的数量级。给定一个预设的新增网格单元阈值，若超过该预设的新增网格单元阈值则认为网格单元(I,J)已被影响；若低于该预设的新增网格单元阈值则认为未受到影响。因此，代表已受影响的网格单元(I,J)的条件可表示为：

式中，ΔW_I ^(n-1)表示n-1时间步网格单元(I,J)的守恒量修正量，ε_a为预设的新增网格单元阈值。

对于子步骤(5.1.1.2)，当且仅当网格单元(I,J)沿I正向的最大信号速度大于0时，扰动才能沿I正向(即右向)传播。扰动沿任一方向传播的最大信号速度为其逆变速度与当地声速之和。因此，代表扰动沿I正向传播的条件为：

式中，S_R表示沿I正向的波速，

表示网格单元(I,J)与网格单元(I+1,J)的交界面在I正向的逆变速度，a表示网格单元(I,J)所在位置的声速，上标I表示I方向，而下标I则表示单元编号(I,J)。逆变速度

中的下标I+1/2表示网格单元(I,J)与网格单元(I+1,J)的交界面，逆变速度

中暗含的方向向量是交界面I+1/2指向网格单元(I+1,J)的单位法向量。

同理，如果扰动是沿I负向(即左向)传播，则子步骤(5.1.1.2)需判断向左传播的最大信号速度。当且仅当沿I负向的最大信号速度大于0时，扰动才能对网格单元(I,J)的左侧网格单元产生影响，即扰动沿I负向传播。因此，代表扰动沿I负向传播的条件可表示为：

此时，S_L表示沿I负向的波速，

表示网格单元(I,J)与网格单元(I-1,J)的交界面在I负向的逆变速度，下标I-1/2表示网格单元(I,J)与网格单元(I-1,J)的交界面，逆变速度

中暗含的方向向量为交界面I-1/2指向网格单元(I-1,J)的单位法向量。

对于子步骤(5.1.1.3)，假设网格单元(I,J)与其附近网格单元尺度相近，则扰动波在当前时间步内沿I方向传播的网格数M可估计为：

式中，ΔS^I表示交界面I-1/2与交界面I+1/2的平均面积；

表示向下取整。分式中，分子(|U^I|+a)Δt表示以最大信号速度在推进时间步长内扰动沿I方向能够传播的距离，分母Ω/ΔS^I表示网格单元沿I方向的尺度。

对于式(5)至(7)，类似的表达也同样可沿J方向推出。

(5.1.2)删除网格单元操作；

删除网格单元操作是从上游开始将已达到收敛要求的网格单元从动态计算域中移除，为求解后期的主导操作。以网格单元(I,J)为例，假设网格单元(I,J)紧邻动态计算域的边界，若要判断网格单元(I,J)是否从动态计算域中删除需同时满足以下三个条件：

(a)网格单元(I,J)的求解已经收敛；

(b)网格单元(I,J)位于最上游；

(c)动态计算域中的其他网格单元不再对网格单元(I,J)产生影响。

以上三个条件都满足才能删除网格单元(I,J)。

对于条件(a)，与新增网格单元操作类似，删除网格单元操作也可以通过衡量网格单元(I,J)在上一个时间步的守恒量修正量来实现。网格单元(I,J)的求解已经收敛意味着它的守恒量修正量应该低于给定的收敛阈值。因此，代表网格单元(I,J)已经收敛的条件可表示为：

条件(b)可以通过判断网格单元的逆变速度来完成。若逆变速度指向动态计算域内，则该紧邻动态计算域边界的网格单元位于最上游。因此，代表位于动态计算域左边界的网格单元(I,J)位于最上游的条件可表示为：

同理，代表位于动态计算域右边界的网格单元(I,J)位于最上游的条件可表示为：

对于动态计算域上、下边界，应分别判断J正向和负向的逆变速度。

条件(c)分为亚声速和超声速两种情况。对于当地马赫数小于1的亚声速流动，当网格单元(I,J)的毗邻网格单元都已收敛，且毗邻网格单元各自对网格单元(I,J)守恒量修正量的影响都小于收敛阈值时，该具有亚声速流动的网格单元(I,J)被视为不再受到毗邻网格单元的影响。因此，亚声速流动中，网格单元(I,J)受其毗邻网格单元影响可忽略的条件可表示为：

式中，Δ(ΔW)表示毗邻网格单元对网格单元(I,J)的守恒量修正量的影响。

而对于当地马赫数大于1的超声速流动，当一个网格单元的依赖域内没有未收敛网格单元，则该网格单元被视为不再受到动态计算域中其他网格单元的影响。以网格单元几何中心的坐标表示网格单元的位置，若网格单元(I,J)指向未收敛相邻网格单元的位置向量与网格单元(I,J)速度矢量的反方向间的夹角大于当地马赫角，那么该相邻网格单元即位于网格单元(I,J)的影响域以外。因此，超声速流动中，网格单元(I,J)不再受动态计算域中其他网格单元影响的条件可表示为：

式中，(x,y)和(u,v)分别表示网格单元(I,J)毗邻网格单元的几何中心坐标和速度分量，(x_I,y_I)和(u_I,v_I)表示网格单元(I,J)的几何中心坐标和速度分量，μ_I表示网格单元(I,J)的马赫角。

(5.2)新增行；

若动态计算域没有包括所有预设网格中的行，则通过使用新增网格单元操作(5.1.1)衡量扰动波沿网格列方向的传播方向与大小，判断扰动波沿网格列方向是否会传出当前动态计算域。若能够传出，则在链表中新增一项并存储新增行的信息。

(5.3)移除链表冗余项；

若链表中出现不包含网格单元的空项，将该空项从链表中移除。

(5.4)再分配存储空间。

若动态计算域有所改变，则重新分配链表的存储空间。

第六步：残差估计；

估计动态计算域内每个网格单元的对流通量、扩散通量和源项，三者之和即为所述的残差。

第七步：时间积分；

使用显式、隐式时间推进方法确定动态计算域内网格单元的守恒量修正量。

第八步：判断是否收敛；

判断动态计算域内所有网格单元守恒量修正量的模值是否均小于所要求的收敛阈值。若不满足，则返回第四步；若满足，则退出迭代计算，继续下一步。

第九步：结果输出。

实施例1：对NACA 0012翼型在来流马赫数0.85、1°攻角条件下进行模拟，流场以远场边界条件进行初始化。网格采用C型拓扑；为了避免有限外边界对结果的影响，给定计算域约为弦长的20倍。图2示意了稀疏的网格、动态计算域的变化以及最终的模拟结果。图2所示的收敛结果中，翼型上下表面的激波均能够清晰地捕捉。本发明提出的扰动区域更新方法所得翼型表面升、阻力系数与传统的全局更新方程结果的相对误差不足万分之一。

图3给出了动态计算域的网格单元量随迭代的变化。动态计算域依据壁面边界条件建立。由于没有激波的阻挡，亚声速流动的动态计算域会随着扰动的传播逐渐扩展至整个给定计算域，如图2中40％总迭代步的情况所示。随着时间推进，上游的网格单元开始收敛，动态计算域也逐渐向下游移动；而在向下游移动的同时，动态计算域的外边界也逐渐向壁面靠近，如图2中95％总迭代步的情况所示。动态计算域的缩减过程表明，随着求解越来越接近定常状态，过大的计算域也会变得越来越冗余。

相比于传统的全局更新方法，扰动区域更新方法节省了27.7％的计算时间。对比图3中两者的最大残差变化曲线发现，在大部分时间内两者最大残差的下降趋势基本一致，而在求解的末阶段，扰动区域更新方法的最大残差迅速下降，这也使得其迭代总步数减少了18.4％。这得益于扰动区域更新方法的动态计算域去除了上游已收敛网格单元，有效避免了上游已收敛网格单元的小幅振荡提高问题的非线性，从而使得求解更容易收敛至定常状态。

实施例2：对来流马赫数为2.5的圆柱绕流进行模拟，流场以远场边界条件进行初始化。图4展示了稀疏网格、显/隐式时间格式的动态计算域变化，以及显/隐式时间格式扰动区域更新方法与激波装配法及隐式全局更新法的数值模拟对比。图4右侧子图的对比表明，本发明提出的扰动区域更新方法与其他方法所得流场的激波位置、流场特性都能够很好的吻合；且无论是显式时间推进还是隐式时间推进，扰动区域更新方法都能够获得与全局更新同等精度的数值结果。

图5给出了动态计算域的网格单元量随迭代的变化。动态计算域依据壁面边界条件建立。在求解初期，如2％总迭代步的情况，新增网格单元操作使得动态计算域随着扰动的传播逐渐增大，并包含整个弓形激波波后流场。在求解后期，如95％总迭代步的情况，删除网格单元操作从上游开始逐渐减少动态计算域中的网格单元，其中被移除的网格单元都已收敛。在求解结束时，动态计算域消失，整个预设计算域中的流场都达到了定常状态。

相比于全局更新法，扰动区域更新方法对于显式、隐式时间推进格式分别节省了55.0％和59.3％的计算时间，最大存储需求分别降低了14.8％和10.4％。此外，扰动区域更新方法不仅降低了单个迭代步的计算量，还降低了迭代的总次数；显式的迭代次数降低了39.8％，隐式降低了36.4％，如图5所示。

综上所述，本发明为计算流体力学中采用时间推进法求解的可压缩流动，提供了一种切实可行、并能显著加快收敛速率、降低存储需求的加速技术。

Claims (2)

1.定常可压缩流动的扰动区域更新方法，其特征在于：所述方法包括如下步骤：

第一步：输入文件载入；

第二步：流场初始化；

流场初始化包括根据远场条件初始化和根据给定流场初始化两种方式；

第三步：动态计算域建立；

对于根据远场条件初始化的方式，动态计算域将根据壁面条件建立；

对于根据给定流场初始化的方式，动态计算域的建立需要衡量给定流场的残差，这种情况中，动态计算域的建立分为两个子步骤：

(3.1)筛选流动特性与来流不同的网格单元，建立初始动态计算域；

流动特性与来流不同的判据表示为：

||W-W_∞||＞ε_d

式中，W表示求解的网格单元的守恒量，W_∞表示来流条件对应的守恒量，ε_d表示收敛阈值；

(3.2)依据残差细化步骤(3.1)所建立的初始动态计算域；

仍旧保留在动态计算域中的网格单元的判据表示为：

式中，R表示网格单元的残差，

为残差的缩放因子；

其余不满足上述判据的网格单元从步骤(3.1)建立的初始动态计算域中去除，得到参与后续流场求解的动态计算域；

第四步：边界虚网格更新；

根据边界条件，为边界虚网格的守恒量赋值；

第五步：动态计算域更新；

(5.1)更新现有行；

遍历所有的现有行，包括新增网格单元与删除网格单元两类操作，通过衡量动态计算域边界处网格单元扰动波的传播方向与大小，判断是否需要扩宽当前动态计算域的现有行；若扰动能够传出，则扩大现有行的范围；若不能传出，则通过衡量收敛性与网格单元间的相互影响，判断是否可以缩小现有行的范围；

(5.2)新增行；

若动态计算域没有包括所有预设网格中的行，通过衡量动态计算域边界处网格单元扰动波的传播方向和大小，判断是否需要在动态计算域中新增行；若扰动能够传入尚不包含在当前动态计算域内的行，则在链表中新增项用于存储新增行的信息；

(5.3)移除链表冗余项；

若链表中出现不包含网格单元的空项，将该空项从链表中移除；

(5.4)再分配存储空间；

若动态计算域有所改变，则重新分配链表的存储空间；

第六步：残差估计；

第七步：时间积分；

使用显式、隐式时间积分方法确定动态计算域内网格单元的守恒量修正量；

第八步：判断是否收敛；

判断动态计算域内所有网格单元的守恒量修正量的模值是否均小于所要求的收敛阈值；若不满足，则返回第四步；若满足，则退出迭代计算，继续下一步；

第九步：结果输出。

2.根据权利要求1所述的定常可压缩流动的扰动区域更新方法，其特征在于：所述的(5.1)更新现有行具体包括：

(5.1.1)新增网格单元操作；

定义I、J表示网格和动态计算域的方向，则(I,J)对应动态计算域中的任意一个网格单元，记为网格单元(I,J)；

以I正向的情况为例，I负向采用同样的方法；为了便于表述，假设右向为网格I标号的增大方向，假设网格单元(I,J)已在之前的时间步被纳入动态计算域但其右侧网格单元还在动态计算域之外，若要确定网格单元(I,J)右侧网格单元是否需要加入动态计算域需要依次经过以下三个子步骤的判断：

(5.1.1.1)判断网格单元(I,J)的流动特性是否已经受到了扰动的影响；

(5.1.1.2)确定影响网格单元(I,J)的扰动是否能够向右传播；

(5.1.1.3)估计在下一个时间步n+1能够受到影响的右侧网格的数目；

其中，若子步骤(5.1.1.1)或(5.1.1.2)为否定结果，则说明右侧网格单元在当前时间步无需加入动态计算域，余下子步骤也无需执行；

子步骤(5.1.1.1)通过衡量网格单元(I,J)在上一个时间步的守恒量修正量来完成；给定一个预设的新增网格单元阈值，若超过该预设的新增网格单元阈值则认为网格单元(I,J)已被影响；若低于该预设的新增网格单元阈值则认为未受到影响；因此，代表已受影响的网格单元(I,J)的条件表示为：

式中，ΔW_I ^(n-1)表示n-1时间步网格单元(I,J)的守恒量修正量，ε_a为预设的新增网格单元阈值；

对于子步骤(5.1.1.2)，当且仅当网格单元(I,J)沿I正向的最大信号速度大于0时，扰动才能沿I正向传播；扰动沿任一方向传播的最大信号速度为其逆变速度与当地声速之和，因此，代表扰动沿I正向传播的条件为：

式中，S_R表示沿I正向的波速，

表示网格单元(I,J)与网格单元(I+1,J)的交界面在I正向的逆变速度，a_I表示网格单元(I,J)所在位置的声速，逆变速度

中的下标I+1/2表示网格单元(I,J)与网格单元(I+1,J)的交界面，逆变速度

中暗含的方向向量是交界面I+1/2指向网格单元(I+1,J)的单位法向量；

同理，如果扰动是沿I负向传播，则子步骤(5.1.1.2)需判断向左传播的最大信号速度，当且仅当沿I负向的最大信号速度大于0时，扰动才能对网格单元(I,J)的左侧网格单元产生影响，即扰动沿I负向传播，因此，代表扰动沿I负向传播的条件表示为：

此时，S_L表示沿I负向的波速，

表示网格单元(I,J)与网格单元(I-1,J)的交界面在I负向的逆变速度，下标I-1/2表示网格单元(I,J)与网格单元(I-1,J)的交界面，逆变速度

中暗含的方向向量为交界面I-1/2指向网格单元(I-1,J)的单位法向量；

对于子步骤(5.1.1.3)，假设网格单元(I,J)与其附近网格单元尺度相近，则扰动波在当前时间步内沿I方向传播的网格数M估计为：

式中，Δt表示当地时间步长，Ω表示网格单元(I,J)的体积，ΔS^I表示交界面I-1/2与交界面I+1/2的平均面积；

表示向下取整；

(5.1.2)删除网格单元操作；

以网格单元(I,J)为例，假设网格单元(I,J)紧邻动态计算域的边界，若要判断网格单元(I,J)是否从动态计算域中删除需同时满足以下三个条件：

(a)网格单元(I,J)的求解已经收敛；

(b)网格单元(I,J)位于最上游；

(c)动态计算域中的其他网格单元不再对网格单元(I,J)产生影响；

以上三个条件都满足才能删除网格单元(I,J)；

对于条件(a)，代表网格单元(I,J)已经收敛的条件表示为：

条件(b)通过判断网格单元的逆变速度来完成，若逆变速度指向动态计算域内，则该紧邻动态计算域边界的网格单元位于最上游，因此，代表位于动态计算域左边界的网格单元(I,J)位于最上游的条件表示为：

同理，代表位于动态计算域右边界的网格单元(I,J)位于最上游的条件表示为：

对于动态计算域上、下边界，应分别判断J正向和负向的逆变速度；

条件(c)分为亚声速和超声速两种情况，对于当地马赫数小于1的亚声速流动，当网格单元(I,J)的毗邻网格单元都已收敛，且毗邻网格单元各自对网格单元(I,J)守恒量修正量的影响都小于收敛阈值时，该具有亚声速流动的单元(I,J)被视为不再受到毗邻网格单元的影响；因此，亚声速流动中，网格单元(I,J)受其毗邻网格单元影响可忽略的条件表示为：

式中，Δ(ΔW)表示毗邻网格单元对网格单元(I,J)的守恒量修正量的影响；

而对于当地马赫数大于1的超声速流动，当一个网格单元的依赖域内没有未收敛网格单元，则该网格单元被视为不再受到动态计算域中其他网格单元的影响；以网格单元几何中心的坐标表示网格单元的位置，若网格单元(I,J)指向未收敛相邻网格单元的位置向量与网格单元(I,J)速度矢量的反方向间的夹角大于当地马赫角，那么该相邻网格单元即位于网格单元(I,J)的影响域以外；因此，超声速流动中，网格单元(I,J)不再受动态计算域中其他网格单元影响的条件表示为：

式中，(x,y)和(u,v)分别表示网格单元(I,J)毗邻网格单元的几何中心坐标和速度分量，(x_I,y_I)和(u_I,v_I)表示网格单元(I,J)的几何中心坐标和速度分量，μ_I表示网格单元(I,J)的马赫角。

Images