[发明专利]一种泥石流方程的有限元数值求解方法

权利要求书

1.一种泥石流方程的有限元数值求解方法，其特征在于，包括以下步骤：

S1、建立工程对象数字模型，并对工程对象数字模型进行计算网格划分；

S2、设置模拟计算时的材料参数、边界条件、模拟时长及时间步长；

S3、将泥石流方程应用于计算网格划分后的工程对象数字模型，并应用数学分裂的方法将所述泥石流方程分裂成扩散方程、对流方程和泥深修正方程；

S4、根据所述材料参数、边界条件以及扩散方程计算得到n+1时刻的第一平均流量过渡值；

S5、根据所述第一平均流量过渡值和对流方程计算得到n+1时刻的第二平均流量过渡值；

S6、根据所述泥深修正方程，结合所述第二平均流量过渡值计算得到n+1时刻的泥深值；

S7、根据所述泥深修正方程得到平均流量修正方程，结合所述第二平均流量过渡值及n+1时刻的泥深值计算得到n+1时刻的平均流量值；

S8、采用动边界处理计算域中的干湿节点；

S9、判断当前时间步是否达到模拟时长，若达到则停止计算，否则返回步骤S2；

所述步骤S2中的材料参数包括重力加速度g、浆体宾汉应力τ_B、泥石流体的容重γ、浆体粘性系数μ_B、泥石流固相体积浓度C_s、液相浆体体积浓度C_f、泥石流整体粘性系数μ_m、泥石流固相颗粒容重ρ_s、泥石流流体密度ρ以及泥石流液相浆体容重ρ_f；

所述步骤S3具体为：

应用数学分裂的方法将泥石流方程：

分裂成扩散方程：

对流方程：

以及泥深修正方程：

对于对流方程(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为浆体的粘性系数，μ为泥石流的综合粘度，C_s为泥石流固相体积浓度，ρ_s为泥石流固相颗粒容重，ρ_f为液相的浆体容重，ρ为泥石流流体密度，为固相摩擦角，Sgn(u_i)＝±1，为二维流速向量在x,y方向上的符号式Sgn(u)＝±1,Sgn(v)＝±1。

2.根据权利要求1所述的泥石流方程的有限元数值求解方法，其特征在于，所述步骤S2中的边界条件设置为：

计算域上、下、左为不可滑移边界，右端出口为自由出流边界；初始流量为0，在半椭圆形容器下设置有泥石流材料，材料的初始高度条件为：

其中x表示x方向坐标，y表示y方向坐标，t表示时间，h(x,y,t)为材料初始高度，椭圆中心坐标(x₀-y₀)＝(5,0)，椭圆形长轴a＝4，短轴b＝2，最大初始高度

3.根据权利要求1所述的泥石流方程的有限元数值求解方法，其特征在于，所述步骤S2中的模拟时长设置为T_max＝30，时间步长设置为Δt＝0.01。

4.根据权利要求1所述的泥石流方程的有限元数值求解方法，其特征在于，所述步骤S4具体为：

对于扩散方程(4)进行化简并向后差分得到：

其中hu_i表示流体在x,y两个方向上的平均流量(hu,hv)，下标i＝1,2；根据公式(9)计算得到n+1时刻的第一平均流量过渡值