[發明專利]一種基于LBM-DDA耦合的堆積體滲流侵蝕破壞模擬計算方法

有效

• 申請號: 202110578213.2
• 申請日: 2021-05-26
• 公開（公告）號: CN113268874B
• 公開（公告）日: 2022-03-22
• 發明（設計）人: 黃達;孟秋杰;宋宜祥;岑奪豐;鐘助
• 申請（專利權）人: 河北工業大學
• 主分類號: G06F30/20
• 分類號: G06F30/20;G06F113/08;G06F119/14
• 代理公司: 天津翰林知識產權代理事務所(普通合伙) 12210
• 代理人: 王瑞

1.一種基于LBM-DDA耦合的堆積體滲流侵蝕破壞模擬計算方法，其特征在于，該方法包括以下步驟：

步驟1、將流場劃分為歐拉網格，并設定歐拉網格中的初始參數，得到LBM流場；根據所要模擬的堆積體塊體的級配，編程生成代表堆積體塊體幾何形狀的塊體，并隨機選擇產生位置，塊體邊界區域內的網格設為拉格朗日網格，定義塊體的接觸參數，得到DDA塊體；進而得到LBM-DDA流固耦合模型；

步驟2、在LBM-DDA流固耦合模型中加入插值格式的邊界條件方程，將DDA塊體邊界對歐拉網格的作用離散到整個LBM流場中，得到DDA塊體對LBM流場的單向作用；

步驟3、完成DDA塊體對LBM流場的單向作用計算后，進行LBM流場的歐拉網格之間的碰撞和遷移，完成LBM流場的計算；

步驟4、根據LBM流場，應用基于伽利略不變的動量交換法計算LBM流場作用于DDA塊體的流體曳力，再將流體曳力作為載荷施加到相應的DDA塊體上；

步驟5、DDA塊體的接觸搜索：判定DDA塊體間的相對狀態；相對狀態包括接觸狀態和張開狀態，接觸狀態包括鎖定狀態和滑動狀態；

若判定結果為接觸狀態，則設置接觸彈簧形成的接觸子矩陣，并對DDA塊體間的接觸面上的粘聚力進行計算；若判定結果為張開狀態，則不設置接觸子矩陣，DDA塊體間不存在接觸面上的粘聚力；在粘聚力的計算中，如果DDA塊體間出現過滑動狀態或張開狀態，則此DDA塊體間的接觸面上的粘聚力消失，在后續的時間步中此DDA塊體間永遠不存在粘聚力；

步驟6、DDA塊體的動力學計算：根據步驟5的判定結果，當DDA塊體間的相對狀態為接觸狀態時，計算DDA塊體的單元剛度矩陣和外力矩陣，得到DDA塊體產生相對滑動或碰撞時的接觸力；然后將該接觸力與步驟4得到的流體曳力進行合成，得到合力以及合力矩；再根據牛頓第二定律求出DDA塊體的加速度，再對時間積分進而求得DDA塊體的速度和位移，從而不斷累加，實現DDA塊體的大位移計算；

步驟7、首先，判斷是否到達最大計算時間，若沒有達到則返回步驟2，若達到，首先通過輸出DDA塊體的坐標，獲取滲流侵蝕破壞后的LBM-DDA流固耦合模型的孔隙率以及模型受到滲流侵蝕破壞后損失的DDA塊體的質量和直徑分布；然后，通過LBM流場得到該模型在受到滲流侵蝕破壞后的滲透系數，進而得到滲流侵蝕破壞后的LBM-DDA流固耦合模型的計算結果。

2.根據權利要求1所述的基于LBM-DDA耦合的堆積體滲流侵蝕破壞模擬計算方法，其特征在于，步驟1中，歐拉網格的初始參數包括分布密度函數和速度邊界條件。

3.根據權利要求2所述的基于LBM-DDA耦合的堆積體滲流侵蝕破壞模擬計算方法，其特征在于，步驟1中，歐拉網格中的分布密度函數由平衡分布密度函數確定，平衡分布密度函數如式(1)所示：

式(1)中，f_i^eq為平衡分布密度函數；ρ為歐拉網格的密度，f_i為分布密度函數；ω_i為權重；e_i為離散速度矢量；v＝(v_x,v_y)代表出入口流速。

4.根據權利要求1所述的基于LBM-DDA耦合的堆積體滲流侵蝕破壞模擬計算方法，其特征在于，步驟2中，插值格式的邊界條件方程如式(2)所示：

式(2)中，f_i'和f′_inv(i)為完成碰撞之后的分布密度函數，下標inv(i)和i代表其速度方向相反；q為某個歐拉網格到DDA塊體邊界的直線距離；δf_i'為DDA塊體表面運動引起的速度矯正：

式(3)中，c_s為網格速度，u_ω為DDA塊體邊界處的速度。

5.根據權利要求1所述的基于LBM-DDA耦合的堆積體滲流侵蝕破壞模擬計算方法，其特征在于，步驟3中，碰撞和遷移過程是：對每個歐拉網格中的分布密度函數在每個時間步長內通過碰撞進行重新分配，再將重新分配后的分布密度函數沿離散速度方向遷移到最近的歐拉網格上。