[發明專利]一種歐拉方程的間斷伽遼金有限元數值求解方法

有效

• 申請號: 202010005524.5
• 申請日: 2020-01-03
• 公開（公告）號: CN111241728B
• 公開（公告）日: 2023-05-05
• 發明（設計）人: 徐立;朱蘭;楊中海;李斌
• 申請（專利權）人: 電子科技大學
• 主分類號: G06F30/23
• 分類號: G06F30/23
• 代理公司: 電子科技大學專利中心 51203
• 代理人: 閆樹平

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基。