基于线化矩阵混合求解模式的高阶CFD隐式时间推进方法

摘要

本发明公开了基于线化矩阵混合求解模式的高阶CFD隐式时间推进方法，属于计算流体力学技术领域。现有的方法通用性非常差，效率低，针对控制方程中的不同项，计算代码都需要重新编写，不利于编程实现，从而导致流体仿真的工作量大，效率低，无法快速有效的得到流体仿真的计算结果和流场信息。本发明通过全微分和链式法则，将依赖项对自由度的导数直接转换为依赖项对变量的导数与变量对自由度导数的乘积，相比传统方法，本发明能够大大减少计算量，提高计算效率，也更具有通用性，利于编程实现，同时兼顾线化矩阵的准确度，能够快速得到流体仿真的结果和流场信息。

权利要求

1.基于线化矩阵混合求解模式的高阶CFD隐式时间推进方法，应用于NACA0012二维翼型中，其特征在于，其包括如下步骤：

步骤一：建立基于非结构混合网格的高精度间断迦辽金计算框架；用于流体动力学的数值仿真计算，其包括Navier-Stokes控制方程、基函数、空间离散，控制方程采用Navier-Stokes方程；

步骤二：建立以Newton/GMRES为基础的隐式迭代方法；

步骤三：通过全微分和链式法则，将依赖项对自由度的导数直接转换为依赖项对变量的导数与变量对自由度导数的乘积，或者转化为依赖项对变量梯度的导数与变量梯度对自由度的导数乘积；采用手动解析和差分扰动混合的方式求解线化矩阵，带入到控制方程，求解得到仿真计算结果；

所述步骤一，其具体包括以下步骤：

步骤101、构建微分形式下的Navier-Stokes方程上式中，t代表时间，矢量U代表方程求解的守恒变量， 表示对守恒变量求一阶导数，代表守恒变量的梯度，代表向量微分算子， 代表对流通量， 代表粘性通量；在后续的公式中，分别用 简化表示 简化表示

步骤102、将流体的计算区域剖分成若干不重合的离散单元κ，在不同的离散单元内，流场中的守恒变量U能够用基函数的线性组合来表示，选择单项式等级基作为基函数，单元内的积分点根据网格单元的类型来确定；

步骤103、在方程两边同时乘以基函数φi并进行积分，得到弱形式下的控制方程(2)；

κ和 代表离散单元和其环量，dV代表离散单元的体积，dS代表离散单元的面积；

其中， 代表无粘通量，采用Riemann通量求解方式， 和 代表粘性通量和提升算子，能够采用内罚方法进行计算； 表示基函数的梯度。

2.如权利要求1所述的基于线化矩阵混合求解模式的高阶CFD隐式时间推进方法，其特征在于，所述步骤二，其具体包括以下步骤：

步骤201、采用Euler后差方法对时间导数进行离散，从而得到控制方程的全离散形式；

上式中，Un代表守恒变量U在第n迭代步下的值，Δt代表迭代步之间的时间间隔；

步骤202、式(3)等效为关于守恒变量U的非线性系统R(U)＝0，采用Newton迭代方法对非线性系统进行迭代求解；

通过Newton迭代，非线性系统转化为一系列线性系统的逼近，式(4)中，R(U)为公式(3)中的左侧部分，即包括时间积分、体积分、面积分和提升算子积分项， 又被称作线化矩阵，Un+1代表第n+1步的值，ΔUn代表两步之间的差。

3.如权利要求2所述的基于线化矩阵混合求解模式的高阶CFD隐式时间推进方法，其特征在于，所述步骤三，其具体包括以下步骤：

步骤301、用基函数的线性组合代替守恒变量U，其具体形式为：其中，φj代表基函数，j为基函数的下标，uj(t)代表对应基函数的系数，又被称作自由度，是一个关于时间的函数，后续的公式中用uj简化表示uj(t)；N代表总的基函数的个数，是一个与计算阶数直接相关的量；

步骤302、将线化矩阵中的守恒变量用(6)式来代替，分别能够得到时间积分、体积分、面积分和提升算子积分项的线化矩阵；

步骤303、计算时间积分项的线化矩阵；

步骤304、利用链式法则和全微分，对体积分项进行分解；

上式中，通量对守恒变量或者守恒变量的梯度的导数能够通过差分微小扰动的方式求出，守恒变量对自由度的偏导数能够转换成基函数，而守恒变量的梯度对自由度的偏导数能够转化为基函数的梯度；

步骤305、同样利用链式法则和全微分，对面积分项和提升算子项进行分解，并计算对应项的线化矩阵；

步骤306、将得到的线化矩阵带入步骤202，更写迭代，得到最终计算结果。

说明书

基于线化矩阵混合求解模式的高阶CFD隐式时间推进方法

技术领域

[0001] 本发明涉及基于线化矩阵混合求解模式的高阶CFD隐式时间推进方法，属于计算流体力学技术领域。

背景技术

[0002] 由于在色散和耗散方面的巨大优势，高阶间断迦辽金方法在计算流体力学(英文缩写：CFD)领域中得到了广泛的研究和应用。考虑到稳定性条件等因素的影响，隐式时间推进方法在高阶间断迦辽金方法的应用最为普遍。而在隐式时间推进的过程中，最核心的内容就是求解每个未知量的线化矩阵。如何能够高效准确的求解线化矩阵，是隐式计算方法中的重要研究内容。

[0003] 以往求解线化矩阵主要有两种方法，手动解析方法和自动微分方法。

[0004] 手工解析求出线化矩阵，这种方法是最准确的，然而对于Navier-Stokes方程中的很多部分，如Riemann通量、边界条件、提升算子等部分，很难通过解析的方式求出其关于自由度的导数，而且该方法计算复杂，公式推导十分繁琐，适用性也一般，效率较低，无法快速有效的得到流体仿真的计算结果和流场信息。

[0005] 自动微分方法依托相应的PDE求解库，针对不同的求导问题，编写相应的代码计算导数项。该方法计算量小，但是代码通用性非常差，效率低。针对控制方程中的不同项，计算代码都需要重新编写，不利于编程实现，从而导致流体仿真的工作量大、效率低，同样无法快速有效的得到流体仿真的计算结果和流场信息。

发明内容

[0006] 针对现有技术的缺陷，本发明的目的在于提供一种通过全微分和链式法则，将依赖项对自由度的导数直接转换为依赖项对变量的导数与变量对自由度导数的乘积，或者转化为依赖项对变量梯度的导数与变量梯度对自由度的导数乘积；能够大大减少计算量，提高计算效率，利于编程实现，同时兼顾线化矩阵的准确度，能够快速得到流体仿真的结果和流场信息的基于线化矩阵混合求解模式的高阶CFD隐式时间推进方法。

[0007] 为实现上述目的，本发明的技术方案为：

[0008] 基于线化矩阵混合求解模式的高阶CFD隐式时间推进方法，其包括如下步骤：

[0009] 步骤一：建立基于非结构混合网格的高精度间断迦辽金计算框架；用于流体动力学的数值仿真计算，其包括Navier-Stokes控制方程、基函数、空间离散，控制方程采用Navier-Stokes方程；

[0010] 步骤二：建立以Newton/GMRES为基础的隐式迭代方法；

[0011] 步骤三：通过全微分和链式法则，将依赖项对自由度的导数直接转换为依赖项对变量的导数与变量对自由度导数的乘积，或者转化为依赖项对变量梯度的导数与变量梯度对自由度的导数乘积；采用手动解析和差分扰动混合的方式求解线化矩阵，带入到控制方程，求解得到仿真计算结果。

[0012] 为了能够简化线化矩阵的求解，提高高精度间断迦辽金方法隐式迭代的效率，本发明提出了基于线化矩阵混合求解模式的高阶CFD隐式时间推进方法，该方法通过全微分和链式法则，将依赖项对自由度的导数直接转换为依赖项对变量的导数与变量对自由度导数的乘积，或者转化为依赖项对变量梯度的导数与变量梯度对自由度的导数乘积。变量(导数)对其自由度的导数，能够通过解析的方式得到。对于Riemann通量或者边界积分，能够采用差分微小扰动的方式进行求解，相比传统方法，本发明能够大大减少计算量，提高计算效率，也更具有通用性，利于编程实现，同时兼顾线化矩阵的准确度，保证隐式迭代有高效的收敛性，能够快速得到流体仿真的结果和流场信息。

[0013] 作为优选技术措施，

[0014] 所述步骤一，其具体包括以下步骤：

[0015] 步骤101、构建微分形式下的Navier-Stokes方程

[0016]

[0017] 上式中，t代表时间，矢量U代表方程求解的守恒变量， 表示对守恒变量求一阶导数， 代表守恒变量的梯度，代表向量微分算子， 代表对流通量， 代表粘性通量；在后续的公式中，分别用 简化表示 简化表示

[0018] 步骤102、将流体的计算区域剖分成若干不重合的离散单元κ，在不同的离散单元内，流场中的守恒变量U能够用基函数的线性组合来表示，选择单项式等级基作为基函数，单元内的积分点根据网格单元的类型来确定；

[0019] 步骤103、在方程两边同时乘以基函数φi并进行积分，得到弱形式下的控制方程(2)；

[0020]

[0021] κ和 代表离散单元和其环量，dV代表离散单元的体积，dS代表离散单元的面积；

[0022] 其中， 代表无粘通量，采用Riemann通量求解方式， 和 代表粘性通量和提升算子，能够采用内罚方法进行计算； 表示基函数的梯度。

[0023] 作为优选技术措施，

[0024] 所述步骤二，其具体包括以下步骤：

[0025] 步骤201、采用Euler后插方法对时间导数进行离散，从而得到控制方程的全离散形式；

[0026]

[0027] 上式中，Un代表守恒变量U在第n迭代步下的值，Δt代表迭代步之间的时间间隔；

[0028] 步骤202、式(3)等效为关于守恒变量U的非线性系统R(U)＝0，采用Newton迭代方法对非线性系统进行迭代求解；

[0029]

[0030] 通过Newton迭代，非线性系统转化为一系列线性系统的逼近，式(4)中，R(U)为公式(3)中的左侧部分，即

[0031]

[0032] 包括时间积分、体积分、面积分和提升算子积分项， 又被称作线化矩阵，上式中Un代表第n步的守恒量的值，Un+1代表第n+1步的值，ΔUn代表两步之间的差。

[0033] 作为优选技术措施，

[0034] 所述步骤三，其具体包括以下步骤：

[0035] 步骤301、用基函数的线性组合代替守恒变量U，其具体形式为：

[0036]

[0037] 其中，φj代表基函数，j为基函数的下标，uj(t)代表对应基函数的系数，又被称作自由度，是一个关于时间的函数，后续的公式中用uj简化表示uj(t)；N代表总的基函数的个数，是一个与计算阶数直接相关的量；

[0038] 步骤302、将线化矩阵中的守恒变量用(6)式来代替，分别能够得到时间积分、体积分、面积分和提升算子积分项的线化矩阵；

[0039] 步骤303、计算时间积分项的线化矩阵；

[0040]

[0041] 步骤304、利用链式法则和全微分，对体积分项进行分解；

[0042]

[0043] 上式中，通量对守恒变量或者守恒变量的梯度的导数能够通过差分微小扰动的方式求出，守恒变量对自由度的偏导数能够转换成基函数，而守恒变量的梯度对自由度的偏导数能够转化为基函数的梯度；

[0044] 步骤305、同样利用链式法则和全微分，对面积分项和提升算子项进行分解，并计算对应项的线化矩阵；

[0045] 步骤305、将得到的线化矩阵带入步骤202，更写迭代，得到最终计算结果。

[0046] 本发明解除了线化矩阵对自由度的依赖，不用对每个变量的自由度进行扰动，只需要对守恒变量或者守恒变量的梯度进行相应的扰动就能求解线化矩阵，减少计算量，同时兼顾线化矩阵的准确度，保证隐式迭代有高效的收敛性，能够快速得到流体仿真的结果和流场信息。相比其他方法，采用本发明方法所需的计算时间有效减少，计算效率得到了提升。

[0047] 与现有技术相比，本发明具有以下有益效果：

[0048] 1、本发明将依赖项对自由度的导数转换为对变量的导数与变量对自由度导数的乘积，以及对变量梯度的导数与变量梯度对自由度的导数乘积，解决了Riemann通量积分和边界积分项的导数计算问题，具备良好的通用性。

[0049] 2、本发明解除了线化矩阵对自由度的依赖，不用对每个变量的自由度进行扰动，只需要对守恒变量或者守恒变量的梯度进行相应的扰动就能求解线化矩阵，减少计算量，提高计算效率，也更具有通用性，利于编程实现，同时兼顾线化矩阵的准确度，保证隐式迭代有高效的收敛性，能够快速得到流体仿真的结果和流场信息。

附图说明

[0050] 图1为本发明方法的流程图。

[0051] 图2为本发明采用Newton/GMRES隐式时间推进的流程图。

[0052] 图3为采用本发明方法计算得到的不同精度下二维圆柱绕流问题的压力分布等值线示图。

[0053] 图4为计算二维圆柱绕流问题时本发明与其他方法在计算量、CPU计算时间上的对比示图。

[0054] 图5为采用本发明计算NACA0012二维翼型问题得到的Mach数和流场分布情况示图。

[0055] 图6为计算NACA0012二维翼型问题时本发明的收敛曲线与传统方法的对比示图。

具体实施方式

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

[0057] 相反，本发明涵盖任何由权利要求定义的在本发明的精髓和范围上做的替代、修改、等效方法以及方案。进一步，为了使公众对本发明有更好的了解，在下文对本发明的细节描述中，详尽描述了一些特定的细节部分。对本领域技术人员来说没有这些细节部分的描述也可以完全理解本发明。

[0058] 如图1-2所示，基于线化矩阵混合求解模式的高阶CFD隐式时间推进方法，包括以下步骤：

[0059] 步骤1、建立基于非结构混合网格求解Navier-Stokes方程的高精度间断迦辽金框架，主要用来进行空气动力学的数值仿真计算，包括Navier-Stokes控制方程、基函数、空间离散。

[0060] 其具体包括以下步骤：

[0061] 步骤101、构建微分形式下的Navier-Stokes方程

[0062]

[0063] 上式中，t代表时间，矢量U代表方程求解的守恒变量， 表示对守恒变量求一阶导数， 代表守恒变量的梯度， 代表向量微分算子， 代表对流通量， 代表粘性通量。在后续的公式中，分别用 简化表示 简化表示

[0064] 步骤102、将流体的计算区域剖分成若干不重合的离散单元κ，在不同的离散单元内，流场中的守恒变量U可以用基函数的线性组合来表示，选择单项式等级基作为基函数，单元内的积分点根据网格单元的类型来确定。

[0065] 步骤103、在方程两边同时乘以基函数φi并进行积分，得到弱形式下的控制方程(2)。

[0066]

[0067] κ和 代表离散单元和其环量，dV代表离散单元的体积，dS代表离散单元的面积。

[0068] 其中， 代表无粘通量，采用Riemann通量求解方式， 和 代表粘性通量和提升算子，

[0069] 可以采用内罚方法进行计算。 表示基函数的梯度。

[0070] 步骤2、构建基于Newton/GMRES的隐式迭代框架。

[0071] 其具体包括以下步骤：

[0072] 步骤201、采用Euler后插方法对时间导数进行离散，从而得到控制方程的全离散形式。

[0073]

[0074] 上式中，Un代表守恒变量U在第n迭代步下的值，Δt代表迭代步之间的时间间隔。

[0075] 步骤202、式(3)等效为关于守恒变量U的非线性系统R(U)＝0，采用Newton迭代方法对非线性系统进行迭代求解。

[0076]

[0077] 通过Newton迭代，非线性系统转化为一系列线性系统的逼近，式(4)中，R(U)为公式(3)中的左侧部分，即

[0078]

[0079] 包含时间积分、体积分、面积分和提升算子积分项， 又被称作线化矩阵，上式中Un代表第n步的守恒量的值，Un+1代表第n+1步的值，ΔUn代表两步之间的差。

[0080] 步骤3：采用混合求解方法计算线化矩阵，带入到控制方程，求解得到仿真计算结果。

[0081] 其具体包括以下步骤：

[0082] 步骤301、用基函数的线性组合代替守恒变量U，其具体形式为：

[0083]

[0084] 其中，φj代表基函数，j为基函数的下标，uj(t)代表对应基函数的系数，又被称作自由度，是一个关于时间的函数，后续的公式中用uj简化表示uj(t)。N代表总的基函数的个数，是一个与计算阶数直接相关的量。

[0085] 步骤302、将线化矩阵中的守恒变量用(6)式来代替，分别可以得到时间积分、体积分、面积分和提升算子积分项的线化矩阵。

[0086] 步骤303、计算时间积分项的线化矩阵。

[0087]

[0088] 步骤304、利用链式法则和全微分，对体积分项进行分解。

[0089]

[0090] 上式中，通量对守恒变量或者守恒变量的梯度的导数能够通过差分微小扰动的方式求出，守恒变量对自由度的偏导数能够转换成基函数，而守恒变量的梯度对自由度的偏导数可以转化为基函数的梯度。

[0091] 步骤305、同样利用链式法则和全微分，对面积分项和提升算子项进行分解，并计算对应项的线化矩阵。

[0092] 步骤305、将得到的线化矩阵带入步骤202，更写迭代，得到最终计算结果。

[0093] 如图3-6所示，应用本发明计算圆柱问题和NACA0012翼型问题的具体实施例：

[0094] 图4为计算二维圆柱绕流问题时，本发明与其他方法在计算量(迭代步数和CPU时间)的对比，纵坐标为计算残差。收敛到相同残差时，本发明在CPU计算时间上具有明显优势。

[0095] 图6为计算NACA0012二维翼型问题时，本发明的收敛曲线与传统方法的对比。可以看到收敛到相同残差的情况下，本发明的计算时间上明显少于传统方法。

[0096] 结论：在计算圆柱问题和NACA0012翼型问题时，相比其他方法，采用本发明方法所需的计算时间有效减少，计算效率得到了提升。

[0097] 如表1所示，采用本发明的方法与手动解析方法和自动微分方法在不同精度下计算量的对比，从表1数据可以看出，本发明在求导计算量方面具备明显优势，且求导次数不会随着精度的提高而增加。

[0098] 表1三维情况下本发明与传统方法线化矩阵计算量的对比

[0099]   一阶精度 二阶精度 三阶精度 四阶精度本发明求导次数 5 5 5 5
手动解析方法求导次数 5 5*4 5*10 5*20
自动微分方法求导次数 5 5*4 5*10 5*20

[0100] 以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。