CN114841095A - 一种不可压缩流动的扰动域推进方法 - Google Patents

Abstract

本发明涉及计算流体力学技术领域，提供了一种不可压缩流动的扰动域推进方法，该方法包括：基于来流速度，建立初始不可压缩流动的计算域和初始流场；基于初始流场，通过SIMPLE方法，进行第一步迭代，获得更新后的速度流场和压力流场；将第一步迭代计算出的当地网格残差作为数值扰动，通过网格单元的位置和距离物体壁面的距离判断，进行扰动域的推进；在扰动域推进后进行下一步迭代循环，获得收敛的解。本发明通过建立随收敛进度变化的动态扰动域，减少了单次迭代过程中的无效计算，避免了大量计算资源的浪费，提高了不可压缩流动数值模拟求解的计算效率。

Description

一种不可压缩流动的扰动域推进方法

技术领域

本发明涉及计算流体力学技术领域，尤其涉及一种不可压缩流动的扰动域推进方法。

背景技术

计算流体力学是计算科学的一个分支，专注于求解流体力学控制方程。科学家和工程师经常通过计算流体力学方法得到流体运动的定量信息，比如在特定条件下流场内不同位置的速度、压力、密度或者温度等参数的值。当速度比较低时，上述流体力学控制方程为不可压缩的纳维-斯托克斯方程，需要求解的物理参数主要包括速度和压力，而求解该方程的主流算法是以SIMPLE（Semi-Implicit Method for Pressure-Linked Equations）系列算法为代表的压力与速度耦合方法。

尽管计算机硬件性能发展迅速，但上述对于航空航天等领域真实飞行器在飞行包线内的性能计算方法，仍需要耗费大量的计算资源和超长周期，难以满足实际工程应用的需求。为了提高数值计算的效率，研究者开展了代数多重网格、并行计算等多方面的研究，在一定程度上提高了计算效率。但随着计算流体力学应用范围越来越广，应用层次逐渐加深，在工程实践中对计算流体力学方法的效率仍有着更高的要求。

扰动域推进方法提出了一种加速求解的新思路，针对可压缩流动的数值模拟求解，传统方法预设固定计算域存在大量无效计算，而扰动域推进方法着眼于计算域在时间推进求解过程中的自适应变化，以期减少无效计算量，提高单个迭代步的收敛速率。然而，这一方法针对的是可压缩流动，不适用于不可压缩纳维-斯托克斯方程的求解。其原因在于对于不可压缩流场，声速趋于无穷大，流场内任一点的扰动可以传遍整个流场，原扰动域推进方法根据声速和当地速度作为扰动传播的判断依据对不可压缩流动而言是失效的。因此，需要针对不可压缩流动的方程性质和求解方法进行深入研究，提出不可压缩流动的扰动域推进方法。

发明内容

有鉴于此，本发明提供了一种不可压缩流动的扰动域推进方法，以解决现有技术中无法提高单次迭代过程中的无效计算以及不可压缩流动数值模拟求解的计算效率的问题。

本发明提供了一种不可压缩流动的扰动域推进方法，包括：

S1基于来流速度，建立初始不可压缩流动的扰动域和初始流场；

S2基于所述初始流场，通过SIMPLE方法，进行第一步迭代，获得更新后的速度流场和压力流场；

S3将由更新后的速度流场和压力流场计算出的当地网格残差作为数值扰动，通过网格单元的位置和距离物体壁面的距离判断，进行扰动域推进；

S4在所述扰动域推进后进行下一步迭代循环，获得收敛的解。

进一步地，所述S1中所述不可压缩流动，包括二维不可压缩流动和三维不可压缩流动。

进一步地，所述不可压缩流动为二维不可压缩流动时，所述S1包括：

根据来流速度，建立20倍飞行器长度的初始网格，并将全部所述初始网格作为初始扰动域；

将不可压缩纳维-斯托克斯方程组在每一网格上解耦分解和离散表示为动量离散方程和压力修正方程；

通过假定一个速度分布和一个压力场，计算出所述动量离散方程中的系数及常数项。

进一步地，所述不可压缩纳维-斯托克斯方程组、所述动量离散方程和所述压力修正方程，表达式分别如下：

其中，x，y为二维的两个方向，u、v为x和y两个方向上的速度，t为时间，P为压力，ρ为密度，μ为动力粘性系数，a、b、a'、b'均为不可压缩纳维-斯托克斯方程组在当前时刻离散得到的方程系数，其中，a包括a _w 、a _N 、a _E 、a _S 、a _P ；b包括b _P ；a'包括a' _P 、a' _w 、a' _N 、a' _E 、a' _S ；b'包括b' _P ；ϕ表示由速度u和v构成的矩阵，包括ϕ _P 、ϕ _W 、ϕ _E 、ϕ _S 、ϕ _N ；p'为修正压力，包括p' _P 、p' _E 、 p' _W 、p' _N 、p' _S ；下标P指当前网格，W、E、S、N指与当前网格相邻的左、右、下、上四个网格。

进一步地，所述S2，包括：

依次在二维方向进行迭代求解动量离散方程，获得二维方向上的新的速度分布；

求解所述压力修正方程，获得修正压力；

通过所述修正压力，获得修正速度；

利用所述修正压力和修正速度，获得更新后的所述压力流场和所述速度流场。

进一步地，所述S3，包括：

根据计算条件和计算需求，设置全局收敛参数；

将所述更新后的压力流场和速度流场，通过不可压缩纳维-斯托克斯方程组，计算求得所述当地网格的残差；

判断所述当地网格的残差对当地网格单元格心与物体壁面的距离的影响是否小于预先设定的删除阈值，若不小于，则不作为；若小于，则进一步判断所述当地网格单元格心与物体壁面的距离是否位于所述删除阈值的扰动域的最外围，若是，删除所有处于所述扰动域中的与所述网格单元紧邻的网格单元；若否，不删除；

循环判断整个扰动域，完成所述扰动域的推进。

进一步地，所述当地网格的残差对当地网格单元格心与物体壁面的距离的影响与所述预先设定的删除阈值的大小关系判断表达式如下：

其中，ɛ _loc 为当地网格的残差，ε _d 为删除阈值，d _loc 为网格单元格心距物体壁面的距离，l为所计算问题的特征尺度。

进一步地，与所述网格单元紧邻的网格单元距离壁面中心的距离，满足小于等于所述网格单元距离壁面中心的距离，表达式如下：

其中，d _adj 为紧邻的网格单元距离壁面中心的距离。

进一步地，所述S4，包括：

统计所有所述当地网格残差，并取所述当地网格残差的最大值；

判断所述当地网格残差的最大值是否小于所述全局收敛参数，若是，则输出结果完成计算；若否，跳转至S2。

进一步地，所述收敛的解，包括：所述更新后的速度流场和压力流场。

本发明与现有技术相比存在的有益效果是：

本发明通过考虑不可压缩流动扰动域中残差数值扰动与物面距离成反比的特点，通过建立随收敛进度变化的动态扰动域，有效减少单次迭代过程中的无效计算，有效提高不可压缩流动数值模拟求解的计算效率。

附图说明

为了更清楚地说明本发明中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1是本发明提供的一种不可压缩流动的扰动域推进方法的流程图；

图2是本发明提供的建立初始不可压缩流动的扰动域和初始流场的流程图；

图3是本发明提供的获得更新后的速度流场和压力流场的流程图；

图4是本发明提供的进行扰动域的推进的流程图。

具体实施方式

以下描述中，为了说明而不是为了限定，提出了诸如特定系统结构、技术之类的具体细节，以便透彻理解本发明实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本发明。在其它情况中，省略对众所周知的系统、装置、电路以及方法的详细说明，以免不必要的细节妨碍本发明的描述。

下面将结合附图详细说明根据本发明的一种不可压缩流动的扰动域推进方法。

图1是本发明提供的一种不可压缩流动的扰动域推进方法的流程图。

如图1所示，该扰动域推进方法包括：

S1，基于来流速度，建立初始不可压缩流动的计算域和初始流场。

S1中不可压缩流动，包括二维不可压缩流动和三维不可压缩流动。

图2是本发明提供的建立初始不可压缩流动的扰动域和初始流场的流程图。

如图2所示，不可压缩流动为二维不可压缩流动时，S1包括：

S11，根据来流速度，建立20倍飞行器长度的初始网格，并将全部初始网格作为初始扰动域。

参考具体计算算例的物理参数，建立20倍特征尺度的初始网格，并将全部网格作为初始扰动域，其中，物理参数指流场初始值来流速度。

S12，将不可压缩纳维-斯托克斯方程组在每一网格上解耦分解和离散表示为动量离散方程和压力修正方程；

S13，通过假定一个速度分布和一个压力场，计算出动量离散方程中的系数及常数项。

假定一个速度分布，分别记为u ⁰ ，v ⁰ ，假定一个压力流场P ^*；以计算动量离散方程中的系数及常数项。

不可压缩纳维-斯托克斯方程组、动量离散方程和压力修正方程，表达式分别如下：

（1）

（2）

（3）

其中，x，y为二维的两个方向，u、v为x和y两个方向上的速度，t为时间，P为压力，ρ为密度，μ为动力粘性系数，a、b、a'、b'均为不可压缩纳维-斯托克斯方程组在当前时刻离散得到的方程系数，其中，a包括a _w 、a _N 、a _E 、a _S 、a _P ；b包括b _P ；a'包括a' _P 、a' _w 、a' _N 、a' _E 、a' _S ；b'包括b' _P ；ϕ表示由速度u和v构成的矩阵，包括ϕ _P 、ϕ _W 、ϕ _E 、ϕ _S 、ϕ _N ；p'为修正压力，包括p' _P 、p' _E 、 p' _W 、p' _N 、p' _S ；下标P指当前网格，W、E、S、N指与当前网格相邻的左、右、下、上四个网格。

S2, 基于初始流场，通过SIMPLE方法，进行第一步迭代，获得更新后的速度流场和压力流场。

图3是本发明提供的获得更新后的速度流场和压力流场的流程图。

如图3所示，S2包括：

S21，依次在二维方向进行迭代求解动量离散方程，获得二维方向上的新的速度分布。

依次在x和y方向上进行迭代求解动量离散方程式（2），获得ϕ=[u*,v*]。

其中，u*和v*是x和y方向上的新的速度。

S22，求解压力修正方程，获得修正压力。

求解压力修正方程式（3），得修正压力

。

S23，通过修正压力

，获得修正速度。

按照式（4），由修正压力

，计算得到修正速度

；

（4）

其中，u _P ’、v _P ’分别为当前网格在x和y方向上的速度修正，

为动量离散方程离散在当前网格所得到的系数，

分别代指左、右、下、上四个网格上的压力修正值。

S24，利用修正压力和修正速度，获得更新后的压力流场和速度流场。

利用修正压力和修正速度，可以得到改进后的压力场和速度场，即

，

，其中，p*表示新的压力。

S3，将由更新后的速度流场和压力流场计算出的当地网格残差作为数值扰动，通过网格单元的位置和距离物体壁面的距离判断，进行扰动域的推进。

图4是本发明提供的进行扰动域的推进的流程图。

如图4所示，S3包括：

S31，根据计算条件和计算需求，设置全局收敛参数。

在本发明中全局收敛参数设置为10^-8。

S32，将更新后的压力流场和速度流场，通过不可压缩纳维-斯托克斯方程组，计算求得当地网格的残差。

将S2计算所得更新后的压力场和速度场带入不可压缩纳维-斯托克斯方程组（1）中，计算得到不满足该方程的部分，即为当地网格的残差ɛ _loc 。

S33，判断当地网格的残差对当地网格单元格心与物体壁面的距离的影响是否小于预先设定的删除阈值，若不小于，则不作为；若小于，则进一步判断当地网格单元格心与物体壁面的距离是否位于删除阈值的扰动域的最外围，若是，删除所有处于扰动域中的与网格单元紧邻的网格单元；若否，不删除。

综合考虑当地网格与物体壁面的距离，判断当地网格的残差对与当地网格距离物体壁面距离的影响是否小于可以删除的阈值，若小于，则判断当地网格单元格心距物体壁面的距离是否位于删除阈值的扰动域的最外围，若是，删除所有处于扰动域中的与网格单元紧邻的网格单元；若否，不删除所有处于扰动域中的与网格单元紧邻的网格单元；若不小于，则无需判断当地网格单元格心距物体壁面的距离是否位于删除阈值的扰动域的最外围。

当地网格的残差对当地网格单元格心与物体壁面的距离的影响与预先设定的删除阈值的大小关系判断表达式如下：

其中，ɛ _loc 为当地网格的残差，ε _d 为删除阈值，d _loc 为网格单元格心距物体壁面的距离，l为所计算问题的特征尺度。

S34，循环判断整个扰动域，完成扰动域的推进。

与网格单元紧邻的网格单元距离壁面中心的距离，满足小于等于与网格单元距离壁面中心的距离，表达式如下：

其中，d _adj 为紧邻的网格单元距离壁面中心的距离。

S4，在扰动域推进后进行下一步迭代循环，获得收敛的解。

S4，包括：

S41，统计所有当地网格残差，并取当地网格残差的最大值。

统计S3中得到的当地网格残差ɛ _loc ，取所有网格残差中的最大值；

S42，判断当地网格残差的最大值是否小于全局收敛参数，若是，则输出结果完成计算；若否，跳转至S2。

本发明通过建立随收敛进度变化的动态扰动域，改善了传统方法对不可压缩流动的数值模拟求解全都是基于固定区域的计算域的固有模式；减少了单次迭代过程中的无效计算，避免了大量计算资源的浪费，提高了不可压缩流动数值模拟求解的计算效率。

上述所有可选技术方案，可以采用任意结合形成本申请的可选实施例，在此不再一一赘述。

应理解，上述实施例中各步骤的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本发明实施例的实施过程构成任何限定。

以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围，均应包含在本发明的保护范围之内。

Claims (10)

1.一种不可压缩流动的扰动域推进方法，其特征在于，包括：

S1基于来流速度，建立初始不可压缩流动的扰动域和初始流场；

S2基于所述初始流场，通过SIMPLE方法，进行第一步迭代，获得更新后的速度流场和压力流场；

S3将由更新后的速度流场和压力流场计算出的当地网格残差作为数值扰动，通过网格单元的位置和距离物体壁面的距离判断，进行扰动域的推进；

S4在所述扰动域推进后进行下一步迭代循环，获得收敛的解。

2.根据权利要求1所述的扰动域推进方法，其特征在于，所述S1中所述不可压缩流动，包括二维不可压缩流动和三维不可压缩流动。

3.根据权利要求2所述的扰动域推进方法，其特征在于，所述不可压缩流动为二维不可压缩流动时，所述S1包括：

根据来流速度，建立20倍飞行器长度的初始网格，并将全部所述初始网格作为初始扰动域；

将不可压缩纳维-斯托克斯方程组在每一网格上解耦分解和离散表示为动量离散方程和压力修正方程；

通过假定一个速度分布和一个压力场，计算出所述动量离散方程中的系数及常数项。

4.根据权利要求3所述的扰动域推进方法，其特征在于，所述不可压缩纳维-斯托克斯方程组、所述动量离散方程和所述压力修正方程，表达式分别如下：

其中，x，y为二维的两个方向，u、v为x和y两个方向上的速度，t为时间，P为压力，ρ为密度，μ为动力粘性系数，a、b、a'、b'均为不可压缩纳维-斯托克斯方程组在当前时刻离散得到的方程系数，其中，a包括a _w 、a _N 、a _E 、a _S 、a _P ；b包括b _P ；a'包括a' _P 、a' _w 、a' _N 、a' _E 、a' _S ；b'包括b' _P ； ϕ表示由速度u和v构成的矩阵，包括ϕ _P 、ϕ _W 、ϕ _E 、ϕ _S 、ϕ _N ；p'为修正压力，包括p' _P 、p' _E 、p' _W 、 p' _N 、p' _S ；下标P指当前网格，W、E、S、N指与当前网格相邻的左、右、下、上四个网格。

5.根据权利要求4所述的扰动域推进方法，其特征在于，所述S2，包括：

依次在二维方向进行迭代求解动量离散方程，获得二维方向上的新的速度分布；

求解所述压力修正方程，获得修正压力；

通过所述修正压力，获得修正速度；

利用所述修正压力和修正速度，获得更新后的所述压力流场和所述速度流场。

6.根据权利要求5所述的扰动域推进方法，其特征在于，所述S3，包括：

根据计算条件和计算需求，设置全局收敛参数；

将所述更新后的压力流场和速度流场，通过不可压缩纳维-斯托克斯方程组，计算求得所述当地网格的残差；

判断所述当地网格的残差对当地网格单元格心与物体壁面的距离的影响是否小于预先设定的删除阈值，若不小于，则不作为；若小于，则进一步判断所述当地网格单元格心与物体壁面的距离是否位于所述删除阈值的扰动域的最外围，若是，删除所有处于所述扰动域中的与所述网格单元紧邻的网格单元；若否，不删除；

循环判断整个扰动域，完成所述扰动域的推进。

7.根据权利要求6所述的扰动域推进方法，其特征在于，所述当地网格的残差对当地网格单元格心与物体壁面的距离的影响与所述预先设定的删除阈值的大小关系判断表达式如下：

其中，ɛ _loc 为当地网格的残差，ε _d 为删除阈值，d _loc 为网格单元格心距物体壁面的距离，l为所计算问题的特征尺度。

8.根据权利要求7所述的扰动域推进方法，其特征在于，与所述网格单元紧邻的网格单元距离壁面中心的距离，满足小于等于所述网格单元距离壁面中心的距离，表达式如下：

其中，d _adj 为紧邻的网格单元距离壁面中心的距离。

9.根据权利要求8所述的扰动域推进方法，其特征在于，所述S4，包括：

统计所有所述当地网格残差，并取所述当地网格残差的最大值；

判断所述当地网格残差的最大值是否小于所述全局收敛参数，若是，则输出结果完成计算；若否，跳转至S2。

10.根据权利要求1所述的扰动域推进方法，其特征在于，所述收敛的解，包括：所述更新后的速度流场和压力流场。

Images