[发明专利]一种欧拉方程的间断伽辽金有限元数值求解方法

权利要求书

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，λ₁λ₂…λ_r0为其整特征值，为对应的特征向量，则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≤LN)的中间变量；

则场值对系统方程展开并带入近似：

其中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按递减顺序排列的奇异值，且有σ₁≥σ₂≥...≥σ_r0；令U＝(φ₁,φ₂,...,φ_N)，分别为N×N和L×L的酉矩阵，其中φ_i(i＝1,2,...,N)为AA^T的特征向量；同样的为A^TA的特征向量；此时由降维的特征系统C＝A^TA得到POD基。