[发明专利]定常可压缩流动的扰动区域更新方法

权利要求书

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)的马赫角。