[發(fā)明專利]一種泥石流方程的有限元數(shù)值求解方法

有效

• 申請?zhí)枺?/td>: 201810313622.8
• 申請日: 2018-04-10
• 公開（公告）號: CN108509741B
• 公開（公告）日: 2022-07-29
• 發(fā)明（設計）人: 王大國;譚國榮;鄭雪平;韓超超
• 申請（專利權）人: 西南科技大學
• 主分類號: G06F30/23
• 分類號: G06F30/23;G06F30/28;G06F111/10;G06F113/08;G06F119/14
• 代理公司: 成都正德明志知識產權代理有限公司 51360
• 代理人: 鮑利蕊

1.一種泥石流方程的有限元數(shù)值求解方法，其特征在于，包括以下步驟：

S1、建立工程對象數(shù)字模型，并對工程對象數(shù)字模型進行計算網格劃分；

S2、設置模擬計算時的材料參數(shù)、邊界條件、模擬時長及時間步長；

S3、將泥石流方程應用于計算網格劃分后的工程對象數(shù)字模型，并應用數(shù)學分裂的方法將所述泥石流方程分裂成擴散方程、對流方程和泥深修正方程；

S4、根據所述材料參數(shù)、邊界條件以及擴散方程計算得到n+1時刻的第一平均流量過渡值；

S5、根據所述第一平均流量過渡值和對流方程計算得到n+1時刻的第二平均流量過渡值；

S6、根據所述泥深修正方程，結合所述第二平均流量過渡值計算得到n+1時刻的泥深值；

S7、根據所述泥深修正方程得到平均流量修正方程，結合所述第二平均流量過渡值及n+1時刻的泥深值計算得到n+1時刻的平均流量值；

S8、采用動邊界處理計算域中的干濕節(jié)點；

S9、判斷當前時間步是否達到模擬時長，若達到則停止計算，否則返回步驟S2；

所述步驟S2中的材料參數(shù)包括重力加速度g、漿體賓漢應力τ_B、泥石流體的容重γ、漿體粘性系數(shù)μ_B、泥石流固相體積濃度C_s、液相漿體體積濃度C_f、泥石流整體粘性系數(shù)μ_m、泥石流固相顆粒容重ρ_s、泥石流流體密度ρ以及泥石流液相漿體容重ρ_f；

所述步驟S3具體為：

應用數(shù)學分裂的方法將泥石流方程：

分裂成擴散方程：

對流方程：

以及泥深修正方程：

對于對流方程(5)，由于泥石流的不可壓縮性，即則有：

式(2)-(8)中，u_i、u_j均表示流體在x,y兩個方向上的垂線平均速度(u,v)，x_i表示流體在x,y兩個方向上的坐標，x_j表示特征線上的位置分量，下標i＝1,2，j＝1,2，且i≠j，g為重力加速度，h為泥深，t為時間，Z₀為原始地形高程，τ_B漿體賓漢應力，γ為泥石流體的容重，μ_B為漿體的粘性系數(shù)，μ為泥石流的綜合粘度，C_s為泥石流固相體積濃度，ρ_s為泥石流固相顆粒容重，ρ_f為液相的漿體容重，ρ為泥石流流體密度，為固相摩擦角，Sgn(u_i)＝±1，為二維流速向量在x,y方向上的符號式Sgn(u)＝±1,Sgn(v)＝±1。

2.根據權利要求1所述的泥石流方程的有限元數(shù)值求解方法，其特征在于，所述步驟S2中的邊界條件設置為：

計算域上、下、左為不可滑移邊界，右端出口為自由出流邊界；初始流量為0，在半橢圓形容器下設置有泥石流材料，材料的初始高度條件為：

其中x表示x方向坐標，y表示y方向坐標，t表示時間，h(x,y,t)為材料初始高度，橢圓中心坐標(x₀-y₀)＝(5,0)，橢圓形長軸a＝4，短軸b＝2，最大初始高度

3.根據權利要求1所述的泥石流方程的有限元數(shù)值求解方法，其特征在于，所述步驟S2中的模擬時長設置為T_max＝30，時間步長設置為Δt＝0.01。

4.根據權利要求1所述的泥石流方程的有限元數(shù)值求解方法，其特征在于，所述步驟S4具體為：

對于擴散方程(4)進行化簡并向后差分得到：

其中hu_i表示流體在x,y兩個方向上的平均流量(hu,hv)，下標i＝1,2；根據公式(9)計算得到n+1時刻的第一平均流量過渡值