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

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

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


【發(fā)明內容】

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

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