一種基於界帶有限元和拉格朗日坐標的流體仿真方法

2023-11-11 01:01:27 1

一種基於界帶有限元和拉格朗日坐標的流體仿真方法
【專利摘要】本發明提出了一種基於界帶有限元和拉格朗日坐標的不可壓縮流體仿真分析方法，是將二維不可壓縮流體的計算域Ω按照傳統有限元網格剖分成Ne個單元，每個單元為Ωi；構造單元Ωi的位移插值場，根據位移插值場構造流體的動力微分方程，求解該動力微分方程得到流體的各種物理參數，從而進行流體的運動分析；其特徵在於用界帶有限單元法構造位移插值場；並基於拉格朗日坐標描述法得到流體的動力微分方程。本發明將拉格朗日坐標方法與界帶有限元方法結合來解決不可壓縮流體的運動仿真問題，目的是利用界帶有限元精度高和拉格朗日坐標下邊界處理方便，通用性好的優勢，提高分析的計算效率和精度。
【專利說明】一種基於界帶有限元和拉格朗日坐標的流體仿真方法

【技術領域】
[0001] 本發明涉及流體仿真分析技術，構建了一種基於界帶有限元和拉格朗日坐標的流 體仿真方法。

【背景技術】
[0002] 流體的仿真分析廣泛應用於各個不同工程領域，如水利工程、航空飛行、高速列車 等。流體的仿真分析，從理論上看主要有歐拉坐標和拉格朗日坐標兩種，從數值分析手段來 說，主要有限差分法和有限單元法兩種。其中有限差分方法使用規則的矩形網格，因此對於 在不規則體中的流體仿真分析需要加密網格，這導致計算量大幅增加。有限單元法在固體 結構仿真分析中得到廣泛運用，該方法的特點是適用於不規則問題，但是需要變分原理。
[0003] 目前，流體的仿真分析主流是基於歐拉坐標，在歐拉坐標系可以導出不可壓縮流 體的運動微分方程，即著名的納維-斯託克斯方程，對於該方程的計算最常用的方法是有 限差分法。但是基於歐拉坐標的流體仿真分析在處理自由面流體問題時，分析困難，計算格 式複雜。由於自由面隨時間不斷變化，其計算區域也不再規則，因此用有限差分分析困難， 而且基於歐拉坐標的流體方程，很難建立變分原理，從而導致有限元分析也很困難。目前僅 能對某些特殊的流體運動問題，建立變分，使用有限元分析。在分析不可壓縮流體時，採用 普通有限元，還存在體積閉鎖等問題，影響計算精度。
[0004] 如果運用拉格朗日坐標，則可以很容易得到流體的動能表達式和勢能表達式，並 可以利用哈密爾頓變分原理導出拉格朗日坐標系的流體動力微分方程。根據動能表達式和 勢能表達式，還可以利用有限元進行分析。但是在拉格朗日坐標下，以流體中質點的位移為 基本未知量，以位移為基本未知量建立的傳統有限單元不能完全滿足不可壓縮條件，計算 精度不好。實際上在固體不可壓縮材料的仿真分析中，有學者提出界帶有限元方法，通過引 入流函數來代替位移作為未知量，並利用界帶單元進行分析，所構造的單元完全滿足不可 壓縮條件。界帶單元可以在傳統有限元網格上進行分析，因此可以利用現有有限元網格剖 分技術，而且具有計算精度高、所需的自由度小等優點。目前還沒有關於流體仿真分析中採 用界帶有限元方法的應用。


【發明內容】

[0005] 本發明針對現有技術不足，提出一種分析二維不可壓縮流體的仿真分析方法，該 方法將拉格朗日坐標方法與界帶有限元方法結合來解決不可壓縮流體的運動仿真問題，目 的是利用界帶有限元精度高和拉格朗日坐標下邊界處理方便，通用性好的優勢，提高分析 的計算效率和精度。
[0006] 為此，本發明提供了一種基於界帶有限元和拉格朗日坐標的流體仿真方法，將二 維不可壓縮流體的計算域Ω按照傳統有限元網格剖分成凡個單元，每個單元為Qi ;構造 單元Qi的位移插值場，根據位移插值場構造流體的動力微分方程，求解該動力微分方程得 到流體的各種物理參數，從而進行流體的運動分析；本發明用界帶有限單元法構造位移插 值場；並基於拉格朗日坐標描述法得到流體的動力微分方程；具體方法如下：
[0007] (a)將二維不可壓縮流體的計算域Ω採用傳統有限元網格剖分成凡個單元，每個 單元為Qi;
[0008] (b)在單元Ωi上建立流函數V(X，y)的插值場時，以單元Ωi為本體，將Ωi周 邊的單元視為Qi的界帶，將這些單元的節點合作進行插值，構造單元Qi上的插值函數 ψ(x,y)；
[0009] (c)將步驟（b)中所述的流函數表達式求偏導，得到流體中在坐標（X，y)處質點的 位移表達式；
[0010] (d)根據上一步的位移場表達式，得到每個單元上流體的質量矩陣Mi,並把所有單 元的質量矩陣計算出來後，通過累加得到總體質量矩陣M;
[0011] (e)根據步驟（c)的位移場表達式，得到流體的剛度矩陣：
[0012] K(Ψ) =Κχ+Κη (Ψ)
[0013] 其中K是剛度矩陣，Ψ是所有節點的流函數值所組成的向量，K1是線性剛度矩陣， Kn是非線性剛度矩陣；
[0014] (f)流體的邊界包括兩種，一種為自由面rf，一種為不可穿過邊界Γη，根據不可 穿過邊界rni所有有限元單元節點的編號，把剛度矩陣和質量矩陣中相應編號的行和列 划去，得到描述流體運動的非線性微分方程：
[0015] Μψ+Κ{ψ)ψ= 0
[0016] (g)利用非線性微分方程求解軟體求解上述非線性微分方程，得到不同時間點上 的流函數向量Ψ;
[0017] (h)根據流函數向量Ψ，得到流體中各個質點在不同時間上的位移；根據流體中 各質點在不同時間點的位移，利用有限元後處理程序，可以得到流體的動態仿真圖；
[0018] (i)在步驟（g)中，剛度矩陣包含非線性和線性兩個部分，如果非線性因素較小而 無需考慮時，可以得到線性的動力微分方程AT#=0，此時利用微分方程求解軟體求解 該線性微分方程，可以得到線性情況下的流函數；
[0019] (j)如果要分析流體的振動模態和頻率，利用特徵值求解軟體求解下面的特徵值 方程
[0020] K1F=ω2ΜΨ
[0021] 其中，ω即為水的振動頻率。
[0022] 在步驟（b)中，界帶單元上流函數的插值函數具體格式為：
[0023] ψ(x,y) =NtΨi；Nt =pT (x,y)P-1
[0024] 其中
[0025] ρτ(χ,y) = (I,χ,y,χ2,xy,y2, ···)
[0026] Pt =[ρ(X1,Y1),P(χ2,12), ···,P(χΝ,yN)]
[0027] Ψτ=(ψ(χ1；Y1),Ψ(χ2,γ2),···, Ψ(χΝ,yN))
[0028] Nt是形函數向量。（Xj,yj)是界帶單元的節點坐標,共N個,包含單元本體節點和 本體單元的周邊單元節點坐標，N是界帶單元的節點總數。pT (x，y)是插值多項式，其項數 與N相同。
[0029] 在步驟（C)中，流體中在坐標（x，y)處質點的位移為：

【權利要求】
1. 一種基於界帶有限元和拉格朗日坐標的流體仿真方法，將二維不可壓縮流體的計算 域Ω按照傳統有限元網格剖分成凡個單元，每個單元為Qi ;構造單元Qi的位移插值場， 根據位移插值場構造流體的動力微分方程，求解該動力微分方程得到流體的各種物理參 數，從而進行流體的運動分析；其特徵在於用界帶有限單元法構造位移插值場；並基於拉 格朗日坐標描述法得到流體的動力微分方程；具體方法如下： (a) 將二維不可壓縮流體的計算域Ω採用傳統有限元網格剖分成凡個單元，每個單元 為Qi; (b) 在單元Ωi上建立流函數Ψ(X，y)的插值場時，以單元Ωi為本體，將Ωi周邊的單 元視為Qi的界帶，將這些單元的節點合作進行插值，構造單元Qi上的插值函數V(X，y); (c) 將步驟（b)中所述的流函數表達式求偏導，得到流體中在坐標（X，y)處質點的位移 表達式； (d) 根據上一步的位移場表達式，得到每個單元上流體的質量矩陣Mi,並把所有單元的 質量矩陣計算出來後，通過累加得到總體質量矩陣M; (e) 根據步驟（c)的位移場表達式，得到流體的剛度矩陣： κ(ψ) =Kx+Kn(￥) 其中K是剛度矩陣，Ψ是所有節點的流函數值所組成的向量，K1是線性剛度矩陣，1是 非線性剛度矩陣； (f) 流體的邊界包括兩種，一種為自由面rf，一種為不可穿過邊界Γη，根據不可穿過 邊界1\上所有有限元單元節點的編號，把剛度矩陣和質量矩陣中相應編號的行和列划去， 得到描述流體運動的非線性微分方程： Μψ + Κ(ψ)ψ = <) (g) 利用非線性微分方程求解軟體求解上述非線性微分方程，得到不同時間點上的流 函數向量Ψ; (h) 根據流函數向量Ψ，得到流體中各個質點在不同時間上的位移；根據流體中各質 點在不同時間點的位移，利用有限元後處理程序，可以得到流體的動態仿真圖； (i) 在步驟（g)中，剛度矩陣包含非線性和線性兩個部分，如果非線性因素較小而無需 考慮時，可以得到線性的動力微分方程+ = 此時利用微分方程求解軟體求解該線 性微分方程，可以得到線性情況下的流函數； (j) 如果要分析流體的振動模態和頻率，利用特徵值求解軟體求解下面的特徵值方程 K1ψ=ω2Mψ 其中，ω即為水的振動頻率。
2. 根據權利要求1所述一種基於界帶有限元和拉格朗日坐標的流體仿真方法，其特徵 在於步驟（b)中，界帶單元上流函數的插值函數具體格式為： Ψ (x, y) = Nt Ψ ^ Nt = pT (x, y) P-1 其中 ρτ (X，y) = (1，X，y，X2, xy，y2,…） Pt =[p(χι,Υι),P(χ2,y2), *··,P(?,yN)] Ψτ =(ψ(χ1；Y1),Ψ(x2,J2),···,Ψ(xN,yN)) Nt是形函數向量；(χ」，y」）是界帶單元的節點坐標，共N個，包含單元本體節點和本體單 元的周邊單元節點坐標，N是界帶單元的節點總數；ρτ (X，y)是插值多項式，其項數與N相 同。
3. 根據權利要求1所述一種基於界帶有限元和拉格朗日坐標的流體仿真方法，其特徵 在於步驟（c)中,流體中在坐標（x，y)處質點的位移為： ii(x,y) = = /VI// -,V(x,j) = - = -I//- 其中Nx和Ny分別表示形函數向量N對坐標求偏導。
4. 根據權利要求1所述一種基於界帶有限元和拉格朗日坐標的流體仿真方法，其特徵 在於步驟（d)中，對單元進行積分以計算質量矩陣，積分方案選擇數值積分，也即 M1 =^ηρ{χη,yn){NxN： +NyN]) H-I 其中，（xn，yn)為單元內部的數值積分點，Wn為積分點的權函數，η是積分點數目。
5. 根據權利要求1所述一種基於界帶有限元和拉格朗日坐標的流體仿真方法，其特徵 在於步驟（e)中，根據步驟（c)的位移場表達式，得到流體自由面上rf的位移u和ν，在重 力作用下，其剛度矩陣的計算表達式為： κ(ψ) =Kx+Kn(￥) K=JrfJrn(V) =Jrf/7(.T，.v)giVTA^[；V? 其中K是剛度矩陣，Ψ是所有節點的流函數值所組成的向量，K1是線性剛度矩陣，&是 非線性剛度矩陣，P是流體的密度，g是重力加速度，rf是流體的自由面。
6. 根據權利要求1所述一種基於界帶有限元和拉格朗日坐標的流體仿真方法，其特徵 在於在步驟（a)中，對流體進行有限元網格剖分，可以採用現有有限元軟體，如Ansys軟體。
7. 根據權利要求1所述一種基於界帶有限元和拉格朗日坐標的流體仿真方法，其特徵 在於在步驟（g)和（i)中，求解微分方程的軟體可以採用現有的軟體包，如大連理工大學開 發的SIPESC軟體，SIPESC軟體是大連理工大學工程力學系開發的工程計算分析軟體平臺； 其功能包括集成開發環境、面向系統集成的活動流程圖定製、工程資料庫管理系統、開放式 結構有限元分析系統、集成優化計算系統等，其中有限元分析系統中集成了方程求解模塊、 有限元後處理模塊等，可以對本發明中微分方程進行求解；但不僅限於SIPESC軟體，只要 採用步驟（g)和（i)中的微分方程進行分析，均在本發明權利要求範圍之內。
8. 根據權利要求1所述一種基於界帶有限元和拉格朗日坐標的流體仿真方法，其特徵 在於步驟（h)中，後處理顯示流體的動態運動的軟體，可以採用現有的後處理軟體，如大連 理工大學開發的SIPESC.POST、SIPESC.Jifex軟體，SIPESC.POST和SIPESC.Jifex均屬於 SIPESC軟體中有限元系統模塊，其功能為有限元結果後處理，可視化動畫顯示等。
9. 根據權利要求1所述一種基於界帶有限元和拉格朗日坐標的流體仿真方法，其特徵 在於在步驟（j)中，求解特徵值方程，可以採用現有軟體，如大連理工大學開發的SIPESC軟 件，但不僅限於SIPESC軟體，只要採用步驟（j)中的特徵值方程進行分析，均在本發明權利 要求範圍之內。
【文檔編號】G06F17/50GK104317985SQ201410483986
【公開日】2015年1月28日 申請日期:2014年9月19日 優先權日:2014年9月19日
【發明者】吳鋒, 徐小明, 陳飆松, 鍾萬勰 申請人:大連理工大學