CN111241728B - 一种欧拉方程的间断伽辽金有限元数值求解方法 - Google Patents

Abstract

本发明属于计算流体力学数值求解领域，涉及一种欧拉方程的间断伽辽金有限元数值求解方法，基于正交分解模型降阶。本发明将间断伽辽金有限元(DG)方法与正交分解(POD)方法结合，通过构造POD基来对求解变量做近似，以此来降低模型维数，利用POD模型降阶把多维的物理过程进行低维的近似描述，降低模型维数以使得DG算法的计算过程中，在维护可接受的精度的同时使得计算模型的自由度减少以降低计算成本，减少时间和内存消耗。最终本发明既保持了有限元和有限体积的优点又克服了它们的不足。

Description

一种欧拉方程的间断伽辽金有限元数值求解方法

技术领域

本发明属于计算流体力学数值求解领域，涉及一种欧拉方程的间断伽辽金有限元数值求解方法，基于正交分解模型降阶。

背景技术

计算流体力学(Computational fluid dynamics，CFD)是现代流体力学的一个重要学科分支，它是计算机科学、流体力学和计算数学以及计算物理相结合的产物。在用CFD解决实际问题的过程中，首先要将实际问题描述成相应的数学模型，一般情况下是给定初值及边界条件的偏微分方程系统，常解的模型为Euler/Navier-Stokes(N-S)方程。在CFD领域，数值方法已经是其求解计算最重要的方式之一，包括有限元法、有限差分法、有限体积法或间断有限元方法等。

间断有限元方法，也称间断伽辽金(Discontinuous Galerkin，DG)方法，是利用完全间断的分片多项式空间作为近似解和试验函数空间的有限元方法，目前该方法已经成熟应用于计算流体力学中。DG方法不仅具有局部守恒、稳定等特性，并且可以通过提高插值函数的阶数来提高精度，所以容易处理复杂的几何形状以及有悬挂节点的不规则网格，并且在不同单元里具有不同度的逼近多项式。这样间断有限元方法既保持了有限元和有限体积的优点又克服了它们的不足。

尽管DG方法有诸多优点，但其有个主要的缺点，尤其对于稳定问题特别敏感，即DG方法在每个单元需要求解的变量数增加且随着精度的提高变量数是非线性的增长，对系统做数值模拟时会导出一个很大的偏微分方程组，导致计算时的自由度太多，需要大量计算时间和内存容量，使得对真实时间估计不可控制,因此极大的限制了其实际应用。而目前在计算流体力学数值求解领域并没有关于优化DG方法缺陷的技术手段。

发明内容

针对上述存在问题或不足，为解决DG方法在求解Euler方程过程中自由度太多而造成的大量计算时间和内存消耗问题，提高计算效率，本发明提供了一种欧拉方程的间断伽辽金有限元数值求解方法，即基于正交分解(Proper Orthogonal Decomposition，POD)模型降阶的间断伽辽金(DG)有限元数值方法。

具体技术方案包括以下步骤：

A、将目标Euler方程用间断Galerkin有限元算法进行离散，建立一个广义系统；

B、根据步骤A建立的广义系统求解出瞬时解的集合U，作为样本进行选取来构造瞬像矩阵A；

C、对步骤B中得到的瞬像矩阵A进行奇异值分解得到互相关矩阵的特征系统C＝A^TA；

D、求出步骤C中所得特征系统C＝A^TA的特征值和对应的特征向量，根据误差界从中选择POD基向量；

E、用步骤D所得的POD基向量对要求解的状态向量做近似得到

将近似的状态向量带入步骤A中的广义系统；

F、对步骤E带入近似的状态向量后的广义系统做Galerkin投影，得到降阶模型并求解。

本发明将DG方法与POD方法结合，通过构造POD基来对求解变量做近似，以此来降低模型维数，利用POD模型降阶把多维的物理过程进行低维的近似描述，降低模型维数以使得DG算法的计算过程中，在维护可接受的精度的同时使得计算模型的自由度减少以降低计算成本，减少时间和内存消耗。最终本发明既保持了有限元和有限体积的优点又克服了它们的不足。

附图说明

图1是本发明的流程图。

图2是实施例步骤B-D求POD基的流程示意图。

图3为实施例三维球DG与POD降维后的压力场图。

图4为实施例三维球降维前后表面场压力系数变化对比。

图5为实施例三维球降维前后表面场相对误差变化。

具体实施方式

下面结合附图来详细说明本发明的技术方案。

本发明属于计算流体力学数值求解领域，涉及一种欧拉方程的间断伽辽金有限元数值求解方法，基于正交分解模型降阶，参照图1，包括以下步骤：

A.将目标Euler方程用间断Galerkin有限元算法进行离散，建立一个广义系统；

三维非定常可压缩欧拉/纳维-斯托克斯(Euler/Navier-Stokes)方程组的守恒形式为：

当ivis＝0时为Euler方程，ivis＝1时为N-S方程。

以三维欧拉(Euler)方程为例，则ivis＝0。式中U为守恒变量，

为无粘通量张量，其具体形式如下：

式中ρ,p分别为密度和压力；x,y,z分别为笛卡尔坐标系下的坐标分量；u,v,w分别为笛卡尔坐标系下的速度分量；e,h分别为内能和焓；E,H分别为总能和总焓。

对于完全气体，有：

其中R表示气体常数，对于空气，一般取287.053焦耳/千克开。c_p,c_v,γ分别为定压比热、定容比热和比热比。

用间断Galerkin有限元算法进行离散，化简后得如下的广义系统：

M为质量矩阵，其元素为m_ij，只与单元坐标及单元类型相关；u_i为单元自由度，通过时间推进求解，R(u)为右端项，u＝(ρ,ρu,ρv,ρw,ρE)^T，t∈[0,T_F]，T_F为总的求解时间。

B.根据步骤A建立的广义系统求解出瞬时解的集合U，作为样本进行选取来构造瞬像矩阵A；

瞬像矩阵的构造对于求解POD基很关键，一般在一个很小的时间区间[0,T₀](T₀远小于T_F)上进行选取，将这个区间不同时刻的瞬时解当做样本。现将时间区间[0,T₀]上分成N个相等子区间，即时间步长

则可得到样本矩阵：

其中n为单元个数，u_i(t_nl)(i＝1,2,...,4n,l＝1,2,...,N)为第i个点上第nl时刻的解。接下来从样本矩阵里选取L个时间点的解构造瞬像矩阵如下：

对于L的确定，我们可以通过样本数据的后验误差来进行判断：

L确定的具体步骤如下：

①设定一个L的初始值L₀，L₀ ²＝O(N)，构造出瞬像矩阵A后计算出样本区间的降维值

(如采用后面的步骤C-F方法)，并与原模型的值做误差

记为error(0)；

②重复①中的步骤求出L₀+1对应的误差error(1)，若error(1)≤error(0)，则循环步骤②直至error(k)＞error(k-1)时停止计算，L₀+k-1即为所求的L。

C.对步骤B中得到的瞬像矩阵A进行奇异值分解得到互相关矩阵的特征系统C＝A^TA，具体如下：

对瞬像矩阵做奇异值(SVD)分解

r为瞬像矩阵的秩，Σ_r×r＝diag(σ₁,σ₂,...,σ_r)，其中σ_i(i＝1,2,...,r)为A按递减顺序排列的奇异值，且有σ₁≥σ₂≥...≥σ_r＞0。令U＝(φ₁,φ₂,...,φ_N)，

分别为N×N和L×L的酉矩阵，其中φ_i(i＝1,2,...,N)为AA^T的特征向量；同样的

为A^TA的特征向量。此时POD基可由降维的特征系统C＝A^TA得到。

D.求出步骤C中所得特征系统C＝A^TA的特征值和对应的特征向量，根据误差界从中选择POD基向量；

一个互相关矩阵的特征系统C＝A^TA，λ₁＞λ₂＞…＞λ_r＞0为其整特征值，

为对应的特征向量，则POD基可定义为：

或者定义为

其中(a_n)_i为特征向量

里的元素，U_n为A的列向量，式(7)与式(8)完全等价。式(7)与式(8)此时的维数为特征系统C＝A^TA的秩，但并非最佳POD基的维数，为估计POD基的维数d，定义如下误差界：

我们要求得到的POD能量小于规定的限度ρ，即：

通常ρ＝10^-3，则我们选择最小整数d使得：

通过(9)式，得到POD基的维数d的估计值，则取前d项构成POD基。

E.用步骤D所得的POD基向量对要求解的状态向量做近似得到

将近似的状态向量带入步骤A中的广义系统；

在由步骤D获得POD基向量之后，我们扩展任意几何形状的流动解。为了说明情况，三维流场的自由度将被表示为：

统一上述守恒变量系数，作如下近似：

u＝ρ,ρu,ρv,ρw,ρE

其中

为降维系统的待求解，

为POD基组成的矩阵，

为降维后维数为d_u(d_u≤r≤L＜＜N)的中间变量。

则场值

对系统方程

展开并带入近似：

其中M为4n×4n(n为单元个数)的对角矩阵，

为4n×5的矩阵。整理后重写(9)式有：

记自由度为d＝d^ρ+d^ρu+d^ρv+d^ρw+d^ρE，其中ψ＝(ψ^ρ,ψ^ρu,ψ^ρv,ψ^ρw,ψ^ρE)为4n×d的矩阵，α(t)为d×5的矩阵：

F.对步骤E带入近似的状态向量后广义系统做Galerkin投影，得到降阶模型并求解；最后，对步骤E中的广义系统，即(11)式做Galerkin投影，则可以得到如下的降维模型：

令M_n＝ψ^T·M·ψ，为d×d的矩阵，

为d×5的矩阵，重新化简(12)式，得到最终的降维系统

由于d＜＜4n，求解系统的迭代矩阵从原来4n×4n的M矩阵降维到维数只有d×d的M_n矩阵，右端项从原来4n×5的

矩阵降维到维数只有d×5的R_n(u)矩阵，所需求解的自由度，也由原来的4n×5个降到d个，从而大大降低了求解系统的维数，减少计算时间，提高计算效率。

实施例：

以一个半径为0.5的球体以及[-20,20]×[-20,20]×[-20,20]的空气盒子为模型验证POD模型降阶的有效性，考虑一个亚声速绕流问题，其马赫数为M_∞＝0.2，在单元数为7893的简单球体绕流模型上进行计算。其场图的结果对比如图3所示：

图3-a为原DG的压力场图，图3-b为POD降维后的压力场图，由图可知两者得到的场图结果一致；图4为POD降维前后表面场压力系数变化的对比，图5为POD降维前后表面场的相对误差变化，由图可知降维前后相对误差控制在0.2％左右，并不影响计算精度。

表1则给出了计算结果

表1三维球降维前后计算结果对比

由表中可见，降维前算法自由度为31572个，POD降维后5个守恒变量的自由度在每个点上都相等，分别为1,2,17,17,1，总自由度仅38个，降维前RKDG算法时间步长为10^-5,POD降维后为5*10^-5,时间步长扩大5倍；而整体计算时间提高了2.5倍。采用该实施例，进一步说明，本发明的方法减少了计算时间，提高了计算效率。

Claims (2)

1.一种欧拉方程的间断伽辽金有限元数值求解方法，包括以下步骤：

A、将目标欧拉Euler方程用间断伽辽金Galerkin有限元算法进行离散，建立一个广义系统；

三维非定常可压缩欧拉/纳维-斯托克斯方程组的守恒形式为：

当ivis＝0时为Euler方程，ivis＝1时为N-S方程；

式中U为守恒变量，

为无粘通量张量，其具体形式如下：

式中ρ,p分别为密度和压力；x,y,z分别为笛卡尔坐标系下的坐标分量；u,v,w分别为笛卡尔坐标系下的速度分量；e,h分别为内能和焓；E,H分别为总能和总焓；

对于完全气体，有：

其中R表示气体常数，对于空气，取287.053焦耳/千克开；c_p,c_v,γ分别为定压比热、定容比热和比热比；

用间断Galerkin有限元算法进行离散，化简后得如下的广义系统：

M为质量矩阵，其元素为m_ij，只与单元坐标及单元类型相关；u_i为单元自由度，通过时间推进求解，R(u)为右端项，u＝(ρ,ρu,ρv,ρw,ρE)^T，t∈[0,T_F]，T_F为总的求解时间；

B、根据步骤A建立的广义系统求解出瞬时解的集合U，作为样本进行选取来构造瞬像矩阵A；

将时间区间[0,T₀]不同时刻的瞬时解当做样本，分成N个相等子区间，T₀远小于T_F，即时间步长

则可得到样本矩阵：

其中u_i(t_nl)为第i个点上第nl时刻的解，n为单元个数，i＝1,2,...,4n；l＝1,2,...,N；接下来从样本矩阵里选取L个时间点的解构造瞬像矩阵如下：

对于L的确定，通过样本数据的后验误差来进行判断：

L确定的具体步骤如下：

①设定一个L的初始值L₀，L₀ ²＝O(N)，构造出瞬像矩阵A后计算出样本区间的降维值

并与原模型的值做误差

记为error(0)；

②重复①中的步骤求出L₀+1对应的误差error(1)，若error(1)≤error(0)，则循环步骤②直至error(k)＞error(k-1)时停止计算，L₀+k-1即为所求的L；

C、对步骤B中得到的瞬像矩阵A进行奇异值分解得到互相关矩阵的特征系统C＝A^TA；

D、求出步骤C中所得特征系统C＝A^TA的特征值和对应的特征向量，根据误差界从中选择POD基向量；

一个互相关矩阵的特征系统C＝A^TA，λ₁>λ₂>…>λ_r>0为其整特征值，

为对应的特征向量，则POD基可定义为：

或者定义为

其中(a_n)_i为特征向量

里的元素，U_n为A的列向量，式(7)与式(8)完全等价；式(7)与式(8)此时的维数为特征系统C＝A^TA的秩，但并非最佳POD基的维数，为估计POD基的维数d，定义如下误差界：

要求得到的POD能量小于规定的限度ρ，即：

通常ρ＝10^-3，选择最小整数d使得：

通过(9)式，得到POD基的维数d的估计值，则取前d项构成POD基；

E、用步骤D所得的POD基向量对要求解的状态向量做近似得到

将近似的状态向量带入步骤A中的广义系统；

在由步骤D获得POD基向量之后，扩展任意几何形状的流动解，三维流场的自由度将被表示为：

统一守恒变量系数，作如下近似：

其中

为降维系统的待求解，

为POD基组成的矩阵，

为降维后维数为d_u(d_u≤r≤L<<N)的中间变量；

则场值

对系统方程

展开并带入近似：

其中M为4n×4n的对角矩阵，n为单元个数，

为4n×5的矩阵，整理后重写(9)式有：

记自由度为d＝d^ρ+d^ρu+d^ρv+d^ρw+d^ρE，其中ψ＝(ψ^ρ,ψ^ρu,ψ^ρv,ψ^ρw,ψ^ρE)为4n×d的矩阵，α(t)为d×5的矩阵：

F、对步骤E带入近似状态向量后的广义系统做Galerkin投影，得到降阶模型并求解。

2.如权利要求1所述欧拉方程的间断伽辽金有限元数值求解方法，其特征在于，所述步骤C具体为：

对瞬像矩阵做奇异值SVD分解

r为瞬像矩阵的秩，Σ_r×r＝diag(σ₁,σ₂,...,σ_r)，其中σ_i(i＝1,2,...,r)为A按递减顺序排列的奇异值，且有σ₁≥σ₂≥...≥σ_r>0；令U＝(φ₁,φ₂,...,φ_N)，

分别为N×N和L×L的酉矩阵，其中φ_i(i＝1,2,...,N)为AA^T的特征向量；同样的

为A^TA的特征向量；此时由降维的特征系统C＝A^TA得到POD基。

Images