[發明專利]一種經GPU加速的浸沒邊界?格子玻爾茲曼流固耦合模擬方法

有效

• 申請號: 201510355061.4
• 申請日: 2015-06-24
• 公開（公告）號: CN104866695B
• 公開（公告）日: 2017-10-24
• 發明（設計）人: 吳家陽;程永光;張春澤;刁偉
• 申請（專利權）人: 武漢大學
• 主分類號: G06F17/50
• 分類號: G06F17/50
• 代理公司: 武漢科皓知識產權代理事務所(特殊普通合伙)42222
• 代理人: 程欣

1.一種經GPU加速的浸沒邊界-格子玻爾茲曼流固耦合模擬方法，其特征在于：步驟如下：

一、在專業前處理軟件中建立待測結構體及其內部流體的浸沒邊界體型，并劃分邊界網格以供程序讀取；

二、程序讀取邊界體型數據，設置計算區域的邊界條件和程序模型參數；

三、根據邊界當前時刻構型和參考構型利用有限元法計算邊界施加在流體上的作用力；

此步驟中：

(1)當浸沒邊界體型為二維時，采用一維線網格離散結構體，同時考慮結構體抵抗拉伸和彎曲的能力，定義在節點上的抗拉和抗彎力密度分別按下式計算：

上式中的(F_s)_l和(F_b)_l分別表示定義在編號為l節點上的抗拉和抗彎作用力線密度；K_s和K_b分別為抗拉系數和彎曲剛度；Δα_l為編號l的節點的原始弧長；X_m+1、X_m-1和X_m分別為編號m+1、m-1和m的節點坐標；δ_ml為克羅內克符號，當m＝l時其取值為1，其余情況下取值為0；n_f為用來離散結構體的總節點數；

(2)當浸沒邊界體型為三維時，采用非結構化三角形網格離散結 構體，

結構體在編號為m節點上的作用力密度F_m的計算公式如下：

其中：W為結構的總能量，為偏微分符號，X_m為節點m的坐標；結構總能量的計算式為：

W＝∫w^sdA

式中的dA為面元面積，而w^s為結構的拉能密度，其計算式為：

其中E為彈性模量，I₁和I₂分別為第一和第二應變不變量；可以證明應變不變量能夠按下式計算：

其中λ₁和λ₂是法平面內沿一對垂直方向的拉伸比，而a、b、c的計算式為：

上式中幾何參數的定義如下：在結構體參考構型和變形后構型上選取兩塊相應的三角形面元，分別在上述三角形中任意選取兩條邊，在參考構型上其邊長分別為l₀和l₀′，而在變形后構型上其邊長分別為l和l′，和分別為上述兩條邊在參考構型和變形后構型上的夾角；

四、將邊界的作用力經由狄拉克函數擴散至周邊流體；

五、利用格子玻爾茲曼方法求解攜帶外力項的納維-斯托克斯方程；

六、同樣經由狄拉克函數將流體速度插值到邊界上得到邊界的移動速度，進而更新邊界位置；

七、重復步驟三-六，直至達到計算終點。

2.如權利要求1所述的經GPU加速的浸沒邊界-格子玻爾茲曼流固耦合模擬方法，其特征在于：所述專業前處理軟件為Gambit或ICEM或ProE。

3.如權利要求1所述的經GPU加速的浸沒邊界-格子玻爾茲曼流固耦合模擬方法，其特征在于：所述步驟四和步驟六中的作用力擴散方程為：

上式中的f(x,t)為t時刻點x處流體受到的作用力密度，F_m為編號為m節點上的作用力密度；分別針對二維和三維條件，Δs_k分別是定義在編號為k的結構體節點上的弧長和面積，而δ(x-X_m)即為狄拉克函數，一般采取如下具體形式：

為了能夠得到結構體形狀的變化過程，我們需要根據結構體網格點的速度更新結構體的構型，在已知流場速度分布條件下，結構體速度同樣將利用狄拉克函數通過從流場插值得到：

上式中的X_m和U_m分別為結構體網格點的坐標和速度，u(x,t)是t時刻 點x處流體的速度，Δh是網格步長。

4.如權利要求1所述的經GPU加速的浸沒邊界-格子玻爾茲曼流固耦合模擬方法，其特征在于：所述步驟五中，格子玻爾茲曼方程為：

式中，g_α(x,t)表示t時刻x處沿離散速度方向ξ_α的分布函數，而g_α(x+ξ_αδ_t,t+δ_t)則為t+δ_t時刻x+ξ_αδ_t處沿離散速度方向ξ_α的分布函數；是t時刻x處沿離散速度方向ξ_α的平衡態分布函數；τ是無量綱的松弛時間，與分子間的粘度系數有關；δ_t為時間步長。