CN114880971B - 一种计算流体力学软件采用的隐式方法 - Google Patents

Abstract

本申请涉及一种计算流体力学软件采用的隐式方法，在物面边界附近的网格单元上采用点LU‑SGS方法，其对角阵上采用谱半径标量近似矩阵，保证了对角占优，使计算更稳定；在除物面边界附近外的绝大多数内部网格单元上采用块LU‑SGS方法，其对角阵上采用雅克比矩阵，更逼近原始离散NS方程，可以大幅提升计算效率，同时还可以克服现有技术中当模拟对象的几何外形越来越复杂，物面边界附近网格质量下降，导致CFD计算稳定性下降，甚至计算失败的问题。

Description

一种计算流体力学软件采用的隐式方法

技术领域

本申请涉及计算机技术领域，特别是涉及一种计算流体力学软件采用的隐式方法。

背景技术

计算流体力学（Computational Fluid Dynamics，CFD）软件，是一种利用计算机进行数值计算，模拟流体流动时的各种相关物理现象，包括流动、热传导、声场等的计算机软件。计算流体力学软件广泛应用于航空航天设计、汽车设计、生物医学工业、化工处理工业、涡轮机设计、半导体设计等诸多工程领域。我国在大型工业级CFD软件的研发上相比美欧国家还存在较明显差距，为解决CFD软件的“卡脖子”难题，国家数值风洞(NNW，NationalNumerical Wind tunnel)工程在2018年启动，专注于以空气动力模拟为核心的CFD软件开发与应用。

CFD软件采用的CFD方法是对描述流场的控制方程，即NS方程，用计算数学的方法将其离散到一系列网格节点上求其离散数值解的一种方法，网格节点包含着流体的相关物理量。上述离散过程包括空间离散和时间离散两个过程。

CFD计算中常用的时间离散方法又包括显式方法和隐式方法。上述显式方法是指离散后的每一个方程只含一个未知数，只需初始时刻方程中所需单元点的值，则可以利用初值进行迭代推进求解，从而得到不同时刻的解。上述隐式方法是指离散后的每一个方程包含多个未知数，不能单独求解，需将方程中所有单元的表达式联立起来进行时间推进求解。显式方法具有简单易实现、内存占用小的优点，但其计算时间步长小、计算稳定性不好；隐式方法可采用较大的时间步长，计算稳定性好。相比显式方法，隐式方法在工程和CFD软件中得到更成熟的应用。上下对称高斯赛德尔方法（Lower-Upper Symmetric Gauss–Seidel，LU-SGS）就是一种隐式方法，主要包括点LU-SGS方法和块LU-SGS方法。二者的区别在于，点LU-SGS方法在对角阵中采用谱半径标量近似矩阵，计算稳定性好，但计算效率低；块LU-SGS方法在对角阵中采用雅克比矩阵，计算效率更高。

随着高性能计算和数值方法的发展，工程应用对CFD数值计算效率和稳定性都提出了更高的要求。目前，大部分商用CFD软件都为国外软件。在国内数字化经济转型升级的背景下，CFD软件的应用得到迅速普及，用户对自主可控CFD软件的计算效率提出了更高的要求。特别是，随着所模拟对象的几何外形越来越复杂，用户对几何外形细节保真度的要求也越来越高，网格生成难度不断增加，物面边界附近网格质量下降，如果仍然采用块LU-SGS方法，极易产生病态的雅克比矩阵，导致CFD计算不稳定，甚至计算失败。

发明内容

基于此，有必要针对上述技术问题，提供一种能够同时提高计算稳定性和计算效率的计算流体力学软件采用的隐式方法。

本发明一个实施例的技术方案如下：

一种计算流体力学软件采用的隐式方法，包括以下步骤：

步骤102，读取模拟对象的网格信息，更新所有网格单元的数据，进行网格初始化，生成网格的壁面距离；

步骤104，判断所述壁面距离是否大于计算域阈值；

步骤106，若所述壁面距离大于所述计算域阈值，则判定为内部区域，并采用块LU-SGS方法对NS方程进行时间离散；

步骤108，若所述壁面距离小于或等于所述计算域阈值，则判定为物面边界附近区域，并采用点LU-SGS方法对NS方程进行时间离散；

步骤110，更新流场变量的数值解，得到下一迭代时间步上的通量值。

在其中一个实施例中，还包括：在上述步骤106中，所述采用块LU-SGS方法对NS方程进行时间离散，具体包括以下步骤：

步骤202，将网格单元中心的通量采用插值方法进行分裂，计算谱半径；

步骤204，块LU-SGS初始化，对后续扫描过程中需使用的中间变量赋零值；

步骤206，将NS方程离散化后的离散方程的左端矩阵分解成三个矩阵之和，即前扫矩阵L、矩阵D及后扫矩阵U；

步骤208，利用上述前扫矩阵L、矩阵D以及后扫矩阵U，依次执行向前扫描、通讯中间量、向后扫描的操作。

在其中一个实施例中，在上述步骤206中将离散方程的左端矩阵分解成三个矩阵之和中，具体包括以下步骤：

步骤302，所述前扫矩阵L取小于当前网格单元序号的相邻网格单元的无粘通量系数雅可比矩阵与该网格单元无粘通量谱半径之差，再加上该网格单元粘性通量系数雅可比矩阵与该网格单元粘性通量谱半径之差；

步骤304：所述矩阵D取当前网格单元的无粘通量系数雅可比矩阵与该网格单元粘性通量系数雅可比矩阵的之差，再加上该网格单元粘性通量的谱半径与粘性通量谱半径的差值；

步骤306：所述后扫矩阵U取大于当前网格单元序号的相邻网格单元的无粘通量系数雅可比矩阵与该网格单元无粘通量谱半径之差，再加上该网格单元粘性通量系数雅可比矩阵与该网格单元粘性通量谱半径之差。

在其中一个实施例中，上述步骤208中利用上述前扫矩阵L、矩阵D以及后扫矩阵U，依次执行向前扫描、通讯中间量、向后扫描的操作，具体包括以下步骤：

步骤402，L步向前扫描运算，即计算某时间迭代步的残差值并作为所述离散方程的右端项，执行L块矩阵运算，进行全场逐点向前扫描得到中间预测值；

步骤404：U步向后扫描运算，即将所述中间预测值作为所述离散方程的右端项，执行U块矩阵运算，进行全场向后扫描得到通量增量值。

在其中一个实施例中，在上述步骤108中，所述采用点LU-SGS方法对NS方程进行时间离散，具体包括以下步骤：

步骤502，将网格单元中心的通量采用插值方法进行分裂，计算谱半径；

步骤504，点LU-SGS初始化，对后续扫描过程中需使用的中间变量赋初值；

步骤506，将NS方程离散化后的离散方程左端矩阵分解成三个矩阵之和，即前扫矩阵L’、对角矩阵D’及后扫矩阵U’；

步骤508，利用上述前扫矩阵L’、对角矩阵D’以及后扫矩阵U’，依次执行向前扫描、通讯中间量、向后扫描的操作。

在其中一个实施例中，在上述步骤506中将离散方程左端矩阵分解成三个矩阵之和中，具体包括以下步骤：

步骤602：对角矩阵D’取当前网格单元无粘通量与粘性通量的谱半径之和，即当前网格单元无粘通量与粘性通量所组成合通量的系数雅可比矩阵的谱半径；

步骤604：前扫矩阵L’取小于当前网格单元序号的相邻网格单元的无粘通量系数雅可比矩阵与该网格单元无粘通量的谱半径之和，再加上该网格单元粘性通量的谱半径；

步骤606：后扫矩阵U’取大于当前网格单元序号的相邻网格单元的无粘通量系数雅可比矩阵与该网格单元无粘通量谱半径之差，再减去该网格单元粘性通量的谱半径。

在其中一个实施例中，上述步骤508中利用上述前扫矩阵L’、对角矩阵D’以及后扫矩阵U’，依次执行向前扫描、通讯中间量、向后扫描的操作中，具体包括以下步骤：

步骤702：L’块向前扫描运算，即计算某时间层上的残差值，执行前扫矩阵L’块算子运算，进行全场逐点向前扫描得到中间预测值；

步骤704：U’块向后扫描运算，即将所述中间预测值作为所述离散方程的右端项，执行后扫矩阵U’块算子运算，进行全场向后扫描得到通量增量值。

本发明的计算流体力学软件采用的隐式方法，在物面边界附近的网格单元上采用点LU-SGS方法，其对角阵上采用谱半径标量近似矩阵，保证了对角占优，使计算更稳定；在除物面边界附近外的绝大多数内部网格单元上采用块LU-SGS方法，其对角阵上采用雅克比矩阵，更逼近原始离散NS方程，可以大幅提升计算效率，同时还可以克服现有技术中当模拟对象的几何外形越来越复杂，物面边界附近网格质量下降，导致CFD计算稳定性下降，甚至计算失败的问题。

附图说明

图1为一个实施例中计算流体力学软件采用的隐式方法的流程示意图；

图2为一个实施例中块LU-SGS方法的流程示意图；

图3为一个实施例中块LU-SGS方法构造扫描算子的流程示意图；

图4为一个实施例中块LU-SGS方法执行扫描操作的流程示意图；

图5为一个实施例中点LU-SGS方法的流程示意图；

图6为一个实施例中点LU-SGS方法构造扫描算子的流程示意图；

图7为一个实施例中点LU-SGS方法执行扫描操作的流程示意图。

具体实施方式

为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。

CFD软件采用CFD方法用于求解NS方程获得流场变量，上述NS方程是Navier-Stokes方程的简称，其是用于描述流体力学的控制方程。CFD方法求解NS方程获得流场变量，通过空间离散和时间离散两个过程，将连续的NS方程离散为Ax=B形式的线性方程组。上述空间离散是将连续的NS方程空间项离散为线性方程组Ax=B的右端项B；上述时间离散是将连续的NS方程时间项离散为线性方程组Ax=B的左端项A，即根据前一时刻的流场变量，计算得到新时刻的流场变量。

CFD方法中的时间离散包括显式方法和隐式方法，其中显示方法是根据前一时刻的流场变量，直接计算依次得到新时刻的流场变量；而隐式方法根据前一时刻的流场变量，形成大型联立方程组，通过求解此方程组同时得到新时刻的所有流场变量。

本发明主要对现有技术中的CFD方法时间离散的隐式方法进行改进。

在一个实施例中，如图1所示，提供了一种计算流体力学软件采用的隐式方法，包括以下步骤：

步骤102，读取模拟对象的网格信息，更新所有网格单元的数据，进行网格初始化，生成网格的壁面距离；

上述模拟对象例如可以是导弹、火箭、汽车、动车、建筑物等。

步骤104，判断所述壁面距离是否大于计算域阈值；

所述计算域阈值可以设定为一个适合的参数，主要用于判断所计算的区域属于内部区域还是物面边界附近的区域；

步骤106，若所述壁面距离大于所述计算域阈值，则判定为内部区域，并采用块LU-SGS方法对NS方程进行时间离散；

所述块LU-SGS方法，其在对角阵中采用雅克比矩阵，可以提高计算效率。雅可比矩阵指函数的一阶偏导数以一定方式排列成的矩阵。

步骤108，若所述壁面距离小于或等于所述计算域阈值，则判定为物面边界附近区域，并采用点LU-SGS方法对NS方程进行时间离散；

所述点LU-SGS方法，其对角阵中采用谱半径标量近似，保证模拟对象的复杂外形时的计算稳定性。谱半径指矩阵所有的特征值绝对值集合的最大值。

步骤110，更新流场变量的数值解，得到下一迭代时间步上的通量值。

在一个实施例中，如图2所示，在上述步骤106中，所述采用块LU-SGS方法对NS方程进行时间离散，具体包括以下步骤：

步骤202，将网格单元中心的通量采用插值方法进行分裂，计算谱半径；

具体的，计算无粘通量的雅可比矩阵，求该雅可比矩阵的所有特征值，选取特征值的最大绝对值作为无粘通量雅可比矩阵的谱半径。计算粘性通量的雅可比矩阵，求该矩阵的所有特征值，选取所有特征值的最大绝对值作为粘性通量雅可比矩阵的谱半径。所述无粘通量指NS方程中非粘性流动流过垂直于法向方向单位面积的质量、动量和能量。所述粘性通量指NS方程中粘性流动流过垂直于法向方向单位面积的质量、动量和能量。

步骤204，块LU-SGS初始化，对后续扫描过程中需使用的中间变量赋零值；

步骤206，将NS方程离散化后的离散方程的左端矩阵分解成三个矩阵之和，即前扫矩阵L、矩阵D及后扫矩阵U。

步骤208，利用上述前扫矩阵L、矩阵D以及后扫矩阵U，依次执行向前扫描、通讯中间量、向后扫描的操作。

在一个实施例中，如图3所示，在上述步骤206中将离散方程的左端矩阵分解成三个矩阵之和中，该步骤用于构造扫描算子，为保证矩阵对角占优，构造近似雅可比矩阵，使得前扫矩阵L、矩阵D及后扫矩阵U具有对角占优性，具体的，可以进一步包括以下步骤：

步骤302，所述前扫矩阵L取小于当前网格单元序号的相邻网格单元的无粘通量系数雅可比矩阵与该网格单元无粘通量谱半径之差，再加上该网格单元粘性通量系数雅可比矩阵与该网格单元粘性通量谱半径之差；

步骤304：所述矩阵D取当前网格单元的无粘通量系数雅可比矩阵与该网格单元粘性通量系数雅可比矩阵的之差，再加上该网格单元粘性通量的谱半径与粘性通量谱半径的差值；

步骤306：所述后扫矩阵U取大于当前网格单元序号的相邻网格单元的无粘通量系数雅可比矩阵与该网格单元无粘通量谱半径之差，再加上该网格单元粘性通量系数雅可比矩阵与该网格单元粘性通量谱半径之差。

在一个实施例中，如图4所示，上述步骤208中利用上述前扫矩阵L、矩阵D以及后扫矩阵U，依次执行向前扫描、通讯中间量、向后扫描的操作，具体还可以包括以下步骤：

步骤402，L步向前扫描运算，即计算某时间迭代步的残差值并作为所述离散方程的右端项，执行L块矩阵运算，进行全场逐点向前扫描得到中间预测值；

步骤404：U步向后扫描运算，即将所述中间预测值作为所述离散方程的右端项，执行U块矩阵运算，进行全场向后扫描得到通量增量值。

在一个实施例中，如图5所示，在上述步骤108中，所述采用点LU-SGS方法对NS方程进行时间离散，具体包括以下步骤：

步骤502，将网格单元中心的通量采用插值方法进行分裂，计算谱半径。

具体的，计算无粘通量的雅可比矩阵，求该雅可比矩阵的所有特征值，选取特征值的最大绝对值作为无粘通量雅可比矩阵的谱半径。计算粘性通量的雅可比矩阵，求该矩阵的所有特征值，选取所有特征值的最大绝对值作为粘性通量雅可比矩阵的谱半径。所述无粘通量指NS方程中非粘性流动流过垂直于法向方向单位面积的质量、动量和能量。所述粘性通量指NS方程中粘性流动流过垂直于法向方向单位面积的质量、动量和能量。

步骤504，点LU-SGS初始化，对后续扫描过程中需使用的中间变量赋零值。

步骤506，将NS方程离散化后的离散方程左端矩阵分解成三个矩阵之和，即前扫矩阵L’、对角矩阵D’及后扫矩阵U’。

具体的，该步骤用于构造扫描算子，主要包括无粘通量作线化处理展开，省略展开后的二阶及高阶项，构造不包含块矩阵的前扫矩阵L’、后扫矩阵U’。

步骤508，利用上述前扫矩阵L’、对角矩阵D’以及后扫矩阵U’，依次执行向前扫描、通讯中间量、向后扫描的操作。

在一个实施例中，如图6所示，在上述步骤506中将离散方程左端矩阵分解成三个矩阵之和中，为保证矩阵对角占优，构造近似雅可比矩阵，使得前扫矩阵L’、对角矩阵D’及后扫矩阵U’具有对角占优性，且对角矩阵D’是标量对角矩阵。具体的，可以进一步包括以下步骤：

步骤602：对角矩阵D’取当前网格单元无粘通量与粘性通量的谱半径之和，即当前网格单元无粘通量与粘性通量所组成合通量的系数雅可比矩阵的谱半径。

步骤604：前扫矩阵L’取小于当前网格单元序号的相邻网格单元的无粘通量系数雅可比矩阵与该网格单元无粘通量的谱半径之和，再加上该网格单元粘性通量的谱半径。

步骤606：后扫矩阵U’取大于当前网格单元序号的相邻网格单元的无粘通量系数雅可比矩阵与该网格单元无粘通量谱半径之差，再减去该网格单元粘性通量的谱半径。

在一个实施例中，如图7所示，上述步骤510中利用上述前扫矩阵L’、对角矩阵D’以及后扫矩阵U’，依次执行向前扫描、通讯中间量、向后扫描的操作中，因为LU-SGS方法在三维情况下无条件稳定，且对角矩阵D’为简单的标量矩阵，其求逆过程为直接除以对角阵D’的对角元素，只需进行前扫矩阵L’、后扫矩阵U’两次扫描，两步计算中均不含块矩阵求逆，可提高求解过程的计算效率。具体地，该步骤可以进一步包括：

步骤702：L’块向前扫描运算，即计算某时间层上的残差值，执行前扫矩阵L’块算子运算，进行全场逐点向前扫描得到中间预测值。

步骤704：U’块向后扫描运算，即将所述中间预测值作为所述离散方程的右端项，执行后扫矩阵U’块算子运算，进行全场向后扫描得到通量增量值。

应该理解的是，虽然图1-7的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1-7中的至少一部分步骤可以包括多个步骤或者多个阶段，这些步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤中的步骤或者阶段的至少一部分轮流或者交替地执行。

本发明的计算流体力学软件采用的隐式方法，在物面边界附近关键单元面上采用点LU-SGS方法，其对角阵上采用谱半径标量近似矩阵，保证了对角占优，使计算更稳定，在除物面边界附近的其他单元面上采用块LU-SGS方法，其对角阵上采用雅克比矩阵，更逼近原始离散NS方程，大幅提升计算效率，同时还可以克服现有技术中当模拟外形复杂，物面边界附近网格质量会下降，导致CFD计算稳定性下降，甚至计算失败的问题。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和易失性存储器中的至少一种。非易失性存储器可包括只读存储器（Read-Only Memory，ROM）、磁带、软盘、闪存或光存储器等。易失性存储器可包括随机存取存储器（Random Access Memory，RAM）或外部高速缓冲存储器。作为说明而非局限，RAM可以是多种形式，比如静态随机存取存储器（Static Random Access Memory，SRAM）或动态随机存取存储器（Dynamic Random Access Memory，DRAM）等。

以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。

Claims (3)

1.一种计算流体力学软件采用的隐式方法，所述方法包括：

步骤102，读取模拟对象的网格信息，更新所有网格单元的数据，进行网格初始化，生成网格的壁面距离；

步骤104，判断所述壁面距离是否大于计算域阈值；

步骤106，若所述壁面距离大于所述计算域阈值，则判定为内部区域，并采用块LU-SGS方法对NS方程进行时间离散；

所述步骤106包括：步骤202，将网格单元中心的通量采用插值方法进行分裂，计算谱半径；步骤204，块LU-SGS初始化，对后续扫描过程中需使用的中间变量赋零值；步骤206，将NS方程离散化后的离散方程的左端矩阵分解成三个矩阵之和，即前扫矩阵L、矩阵D及后扫矩阵U；步骤208，利用上述前扫矩阵L、矩阵D以及后扫矩阵U，依次执行向前扫描、通讯中间量、向后扫描的操作，得到中间预测值以全场通量值；

所述步骤206包括：步骤302，所述前扫矩阵L取值为：“小于当前网格单元序号的相邻网格单元”的无粘通量系数雅可比矩阵与该网格单元无粘通量谱半径之差，再加上该网格单元粘性通量系数雅可比矩阵与该网格单元粘性通量谱半径之差；步骤304：所述矩阵D取值为：当前网格单元的无粘通量系数雅可比矩阵与该网格单元粘性通量系数雅可比矩阵的之差，再加上该网格单元粘性通量的谱半径与粘性通量谱半径的差值；步骤306：所述后扫矩阵U取值为：“大于当前网格单元序号的相邻网格单元”的无粘通量系数雅可比矩阵与该网格单元无粘通量谱半径之差，再加上该网格单元粘性通量系数雅可比矩阵与该网格单元粘性通量谱半径之差；

步骤108，若所述壁面距离小于或等于所述计算域阈值，则判定为物面边界附近区域，并采用点LU-SGS方法对NS方程进行时间离散；

所述步骤108包括：步骤502，将网格单元中心的通量采用插值方法进行分裂，计算谱半径；步骤504，点LU-SGS初始化，对后续扫描过程中需使用的中间变量赋初值；步骤506，将NS方程离散化后的离散方程左端矩阵分解成三个矩阵之和，即前扫矩阵L’、对角矩阵D’及后扫矩阵U’；步骤508，利用上述前扫矩阵L’、对角矩阵D’以及后扫矩阵U’，依次执行向前扫描、通讯中间量、向后扫描的操作，得到中间预测值以全场通量值；

所述步骤506包括：步骤602：对角矩阵D’取值为：当前网格单元无粘通量与粘性通量的谱半径之和，即当前网格单元无粘通量与粘性通量所组成合通量的系数雅可比矩阵的谱半径；步骤604：前扫矩阵L’取值为：“小于当前网格单元序号的相邻网格单元”的无粘通量系数雅可比矩阵与该网格单元无粘通量的谱半径之和，再加上该网格单元粘性通量的谱半径；步骤606：后扫矩阵U’取值为：“大于当前网格单元序号的相邻网格单元”的无粘通量系数雅可比矩阵与该网格单元无粘通量谱半径之差，再减去该网格单元粘性通量的谱半径；

步骤110，更新流场变量的数值解，得到下一迭代时间步上的通量值。

2.根据权利要求1所述的方法，其特征在于，所述步骤208包括：

步骤402，L步向前扫描运算，即计算某时间迭代步的残差值并作为所述离散方程的右端项，执行L块矩阵运算，进行全场逐点向前扫描得到中间预测值；

步骤404：U步向后扫描运算，即将所述中间预测值作为所述离散方程的右端项，执行U块矩阵运算，进行全场向后扫描得到通量增量值。

3.根据权利要求1所述的方法，其特征在于，所述步骤508包括：

步骤702：L’块向前扫描运算，即计算某时间层上的残差值，执行前扫矩阵L’块算子运算，进行全场逐点向前扫描得到中间预测值；

步骤704：U’块向后扫描运算，即将所述中间预测值作为所述离散方程的右端项，执行后扫矩阵U’块算子运算，进行全场向后扫描得到通量增量值。

Images