[發明專利]一種細長柔性絲狀顆粒模型的構建及其數值模擬方法

在審

• 申請號: 201911305748.1
• 申請日: 2019-12-18
• 公開（公告）號: CN111062132A
• 公開（公告）日: 2020-04-24
• 發明（設計）人: 張二強;趙敏;潘廣樂;宋祖國;王永紅;史建新;趙漢文;王宗英
• 申請（專利權）人: 陜西中煙工業有限責任公司
• 主分類號: G06F30/20
• 分類號: G06F30/20;G06F111/10
• 代理公司: 北京科億知識產權代理事務所(普通合伙) 11350
• 代理人: 宋秀珍;翟小梅

1.一種細長柔性絲狀顆粒模型的構建及其數值模擬方法，其特征在于：包括以下步驟；

1)依據細長柔性絲狀顆粒的物理結構特性，基于DEM模型，構建細長柔性絲狀顆粒模型，添加細長柔性絲狀顆粒模型約束力方程，建立細長柔性絲狀顆粒模型的運動方程；

2)通過CAD繪圖軟件畫出三維實體模型，并使用Pointwise依據三維實體模型劃分出模擬區域的計算網格節點；

3)利用OpenFoam對所述網格建立基本物理模型：基本控制方程，即連續性方程、動量方程、能量方程，并建立壁面碰撞模型；

4)定義細長柔性絲狀顆粒的物性參數，所述的物性參數包括：密度ρ、厚度h、寬度d、長度l、滾動摩擦系數μ_r；

5)將構建的柔性絲狀顆粒模型與開源的計算流體力學數值模擬軟件OpenFoam耦合，實現柔性絲狀顆粒模型離散模擬的并行模擬計算；

6)定義通道進出口邊界條件及壁面條件；

7)定義初始條件：柔性絲狀顆粒的數量，柔性絲狀顆粒模型的初始位置、初始速度，定義氣相的初始速度及密度；

8)對步驟3)的基本控制方程進行離散化，并采用步驟6)和步驟7)中定義的邊界條件和初始條件進行封閉和求解；

9)對整個計算區域初始化，設定時間步長和模擬時長，對計算區域內代數方程組進行反復迭代，直到滿足設定的模擬時長，完成柔性絲狀顆粒流態化過程的數值模擬，并運用時間步保存機制對計算結果進行保存；

10)對計算進行后處理分析和顯示。

2.根據權利要求1所述的一種細長柔性絲狀顆粒模型的構建及其數值模擬方法，其特征在于：上述步驟1)中，在DEM模型中，顆粒被假設成剛體，DEM模型處理顆粒間碰撞時，顆粒與顆粒的相互作用允許發生微小重疊，顆粒只能與相鄰的顆粒發生作用，單個顆粒所受作用可以由所有顆粒對其作用的加和；DEM模型中顆粒的碰撞作用力采用彈簧-阻尼器-滑動摩擦器模型進行處理計算；

采用Hooke模型計算顆粒間作用力，單個球型顆粒間的作用力的計算方法運動方程及其運動如下：

碰撞作用力由法向力和切向力兩部分組成，假設顆粒i質心到顆粒j質心的方向為法方向，則顆粒i受到的法向力F_n和切向力F_t為：

F_n＝k_nδn_ij-γ_nv_nij (1)

F_t＝k_tδt_ij-γ_tv_tij，max(|F_t|)＝|μF_n| (2)

k_n法向上的剛性系數，γ_n法向為阻尼系數，k_t法向上的剛性系數，γ_t切向為阻尼系數，δn_ij為顆粒間的重疊量，δt_ij切向力產生的顆粒位移，μ為顆粒的滾動摩擦系數；v_nij是碰撞顆粒相對速度在法向上的分量，v_tij是碰撞顆粒之間的滑移速度；

Hooke模型中，δt_ij忽略不計，則k_n，γ_n，k_t，γ_t的計算如下：

k_n＝k_t (6)

γ_n＝γ_t (7)

這4個系數根據相互碰撞的兩個顆粒的屬性與接觸狀態計算：

其中R為顆粒半徑，m為顆粒質量，e為系數，Y為楊氏模量，v為泊松比；

在DEM模型的基礎上，通過在分組的顆粒間彼此添加一個作用勢函數，通過添加作用勢函數將顆粒捆綁在一起，使之形成柔性絲狀顆粒；勢函數形式如下：

Ψ_ij(r_ij)＝k(r_ij-r_o) (19)

其中r_ij表示i和j兩個原子之間的距離；Ψ_ij是對勢，r_o為截斷半徑，這里相當于硬球半徑，k為作用系數；

具體捆綁方法：通過對顆粒間添加一個作用勢函數，顆粒在運動的過程中，當相鄰的顆粒間超過了設定的截斷半徑，顆粒間添加的勢函數就會起作用，對兩者施加作用力，將兩者束縛在一起，當相鄰的顆粒間小于設定的截斷半徑，顆粒在DEM模型受力的基礎上，額外添加勢函數同樣對兩者施加作用力，將兩者排斥開來，使顆粒間距恢復到設定的截斷半徑；通過捆綁方法，最終形成柔性絲狀煙絲顆粒模型；

除DEM模型以及顆粒間的作用勢函數外，顆粒平動和轉動的運動方程遵循牛頓第二定律，分別為：

其中m和I表示顆粒的質量和轉動慣量；v和ω表示顆粒的平動和轉動速度；t為時間；g表示重力加速度；R表示從顆粒質心指向碰撞接觸點的矢量，其值為顆粒半徑；

對顆粒位置與速度方程的求解過程，采用了蛙跳算法進行迭代更新，蛙跳算法用t時刻的位置與t-Δt/2時刻的速度，更新位置與速度，用更新t時刻的位置計算得到的作用力F(t)，公式如下：