CN112818574B - 模拟泥石流起动形成、流动发展和再次淤积的数值方法 - Google Patents

Abstract

模拟泥石流起动形成、流动发展和再次淤积的数值方法，包括将各物质点上的质量和动量通过高阶B样条基函数映射到相应的背景网格，计算背景网格节点处的质量和动量，分别计算背景网格节点上的内力

和外力

引入人工阻尼力，在背景网格节点上动量方程的求解，更新背景网格节点的加速度和速度，更新后的背景网格节点的速度和加速度分别映射回相应的物质点，更新各物质点的速度和位置，输出所需物质点的相关物质信息；本发明引入人工阻尼力模拟泥石流运动过程中颗粒间的相互摩擦，引入固液转化本构模型准确模拟泥石流的瞬时运动性状态，利用数值模拟结果，对泥石流灾害的预测、评估和防治，避免或降低泥石流对各工程领域灾害的预测和防治具有重要意义。

Description

模拟泥石流起动形成、流动发展和再次淤积的数值方法

技术领域

本发明属于地质灾害技术领域，具体涉及一种模拟泥石流起动形成、流动发展和再次淤积运动全过程的数值模拟方法。

背景技术

泥石流是一种掺混大量碎屑物和液态水的多相流体，具有巨大的运动动能和极强的破坏力，其运动全过程包括起动形成、流动发展和再次淤积等不同阶段。泥石流起动形成阶段是坡体由非连续变形转化为连续变形、物理状态由固态转化为液态，具有明显的固-液相变现象，在流动发展和再次淤积阶段，随着剪切速率的变化，泥石流在宏观上表现出类固体或类液态的力学行为。

目前的模拟方法采用的基于网格类的数值方法，如有限元法和有限体积法等，受限于网格畸变和复杂的网格重构过程，难以对泥石流的大变形破坏过程进行有效模拟；基于非连续介质力学的离散元法(DEM)等方法，虽然在滑坡和泥石流等领域得到了一定的应用，但这类方法计算量较大，限制了其计算规模；基于连续介质力学的无网格粒子类方法，对该类问题的模拟具有天然的优势，如光滑粒子流体动力学方法(SPH)在滑坡、泥石流等各类土体颗粒物流动问题中得到了广泛应用，但SPH需要通过耗时的邻域粒子搜索算法来实现控制方程的建立和求导，不利于三维及大规模问题的高效求解。

物质点法(MPM)是一种结合了拉格朗日描述和欧拉描述的无网格粒子类方法，其通过拉格朗日型的物质点离散求解域，基于欧拉背景网格实现各物质点之间的相互作用和联系，其控制方程的建立不依赖于网格，在有效避免网格畸变的同时可以大大提高计算效率，已在各类大变形问题中展现了强大的求解能力和求解优势。然而，传统物质点法由于采用线性插值形函数，其插值形函数导数在网格边界处不连续，因此当物质点跨越网格边界时将引起“网格穿越误差”而造成速度和应力求解精度的降低；同时，传统物质点法难以统一有效描述泥石流运动过程中的固液转化相变现象。因此，建立一种准确模拟泥石流从起动、发展到再次淤积运动全过程以揭示其起动机理，明确其演化规律，对泥石流灾害的预测和防治具有重要的现实意义。

发明内容

针对现有技术的不足，本发明提供一种基于高阶B样条基函数和固液转化本构模型的高阶B样条物质点法，通过人工阻尼力模拟运动过程中土体颗粒间的相互摩擦，能够准确、瞬时动态模拟泥石流起动形成、流动发展和再次淤积运动的全过程的数值方法。

为实现上述目的，本发明采用以下技术方案：

模拟泥石流起动形成、流动发展和再次淤积的数值方法，该方法包括以下步骤：

步骤1、将各物质点上的质量和动量通过高阶B样条基函数映射到相应的背景网格，计算背景网格节点处的质量和动量：

式中：上标k表示计算时间步长；下标J和p分别为背景网格节点和物质点上的变量；S为高阶B样条基函数；m为质量；x和

分别为位置和位移矢量；

在物质点法中，求解域由一系列拉格朗日型物质点粒子离散，物质点粒子携带整个求解域的物理信息，包括应力、应变、质量、密度等，物质点粒子在整个计算过程中随求解域的变形而运动；同时，通过欧拉背景网格和插值形函数实现各物质之间的相互作用和联系；

步骤2、分别计算背景网格节点上的内力f^int和外力f^ext：

式中：▽S表示高阶B样条基函数的梯度；σ为柯西应力张量；

和b分别表示面力和体力矢量；h表示面力边界厚度；V为体积；

步骤3、引入人工阻尼力，在背景网格节点上进行动量方程的求解，更新背景网格节点的加速度和速度：

式中：

和

分别为速度和加速度矢量；f^damp为人工粘性力；

其中，为了消除颗粒物流动过程中的非物理数值振荡，引入一种线性人工阻尼力为：

f^damp＝-ζ_d·m_p·v_p (7)

式中：m_p和v_p分别表示物质点p的质量和速度；ζ_d为人工阻尼系数，求解公式为：

其中：ξ为无量纲的阻尼系数,h为背景网格尺寸，对一维、二维和三维问题，分别取h＝dx,

和

E和ρ分别为物质点的弹性模量和密度；

人工阻尼力f^damp的大小与物质点的运动速度和声速成正比，且方向与运动速度相反，即：人工阻尼力可以看作是颗粒材料在运动过程中，颗粒与颗粒之间的相互摩擦力，其中，无量纲阻尼系数ξ通过模型试验确定，相较人工粘性力具有更明确的物理意义；

步骤4、将更新后的背景网格节点的速度和加速度分别映射回相应的物质点，更新各物质点的速度和位置：

步骤5、将更新后的物质点速度再次映射回相应的背景网格节点，得背景网格节点速度

步骤6、利用背景网格节点速度

计算各物质点的变形率张量和应变增量：

步骤7、更新物质点的密度；

步骤8、基于固液转化本构模型更新各物质点的应力：

(ρ_p)^k+1＝(ρ_p)^k/[1+(Δε_p)^k+1] (13)

a)若物质点的密度小于密度阈值时则视为流体状态，其中密度阈值一般取1800kg/m³～2000kg/m³，即为泥水混合物的密度，步骤8采用基于非牛顿流体本构模型更新其应力，其更新过程如下：

基于引入广义Cross非牛顿流体本构模型以描述泥石流的非牛顿流体流动，其等效粘度系数写为：

式中：μ₀和μ_∞分别表示流体在低剪切率和高剪切率下的粘度系数；K和m为常系数；利用等效粘度系数μ_eff，可将广义Cross模型转化为牛顿流体模型形式：

得总的应力张量为：σ_ij＝-Pδ_ij+τ_ij，其中，基于人工状态方程求解流体的静水压力P：

式中：c为人工声速；ρ₀和ρ分别为流体的初始密度和当前密度；r为常数，一般取为7；

b)若物质点密度大于等于密度阈值则视为固体状态，步骤8采用基于Drucker-Prager弹塑性本构模型更新其应力，其更新过程如下：

式中：I₁＝σ_kk为应力张量的第一不变量；J₂＝s_ijs_ij/2为应力偏量的第二不变量；q_φ为摩擦系数，反映了压力对屈服极限的影响程度；k_φ为纯剪切状态下时的屈服应力；q_φ和k_φ可由材料的内聚力c和摩擦角φ确定；通过选取参数：

和

则Drucker-Prager屈服面可在π平面上外接Mohr-Coulomb屈服面，在Drucker-Prager模型中，剪切势函数采用非关联塑性流动法则，即：ψ^s＝τ+q_ψσ_m

式中：材料常数q_ψ由剪胀角ψ确定；q_ψ与ψ之间的关系和q_φ与φ之间的关系相同；如果q_ψ＝q_φ，则为关联塑性流动；如果q_ψ＝0剪胀角为零，材料满足塑性不可压缩条件；

步骤9、根据需要选择是否重复步骤2至步骤8输出所需物质点的相关物质信息或结束计算。

本发明的有益效果：1)本发明在现有物质点法基础上，引入高阶B样条插值形函数以消除物质点法的网格穿越误差，引入人工阻尼力模拟泥石流运动过程中颗粒间的相互摩擦，引入固液转化本构模型能够准确模拟泥石流的瞬时不同运动性状态；2)本发明能够统一模拟泥石流的起动、加速流动、减速流动和再次淤积的整个运动过程，对分析泥石流的动力学行为具有重要意义；3)本发明利用数值模拟结果明确其演化规律，指导泥石流灾害的预测、评估和防治，避免或降低泥石流对各工程领域灾害的预测和防治具有重要的现实意义。

附图说明

图1为本发明流程示意图；

图2为本发明同一网格划分下，物质点法线性插值形函数和高阶次B样条基函数比较；

图3-1为本发明实施例1土体颗粒物自由坍塌模型初始构形；

图3-2为本发明实施例1土体颗粒物自由坍塌模型最终塌落构形；

图4为本发明实施例1无量纲坍塌前缘位置X*和动能K*随无量纲时间T*的变化规律示意图；

图5为本发明实施例1不同时刻下物质点法模拟所得土体坍塌速度云图；

图6为本发明实施例2膨润土沿斜面滑动模型正视图；

图7为本发明实施例2无量纲坍塌前缘位置X*和动能K*随无量纲时间T*的变化规律示意图；

图8为本发明实施例2不同时刻下物质点法模拟所得土体滑坡速度云图。

具体实施方式

下面结合附图和具体实施列对本发明作进一步详细说明。

实施例1

图1为本发明的算法流程图；图2比较了相同网格划分情况下物质点法线性插值形函数和高阶次B样条基函数，可以看出：相较传统物质点法的线性插值形函数，本发明给出了更多数目和更为光滑的高阶次B样条基函数。

如图3-1、图3-2所示，在长、宽、高分别为l₀×b₀×h₀＝0.2m×0.1m×0.14m的立方体无粘砂土，初始时刻静止于容器左下角，基于本发明模拟其在自重作用下的自由坍塌。

分别定义无量纲时间T^*、无量纲坍塌前缘位置X^*和无量纲的系统动能K^*为：

X^*＝l_s/h₀和

式中：t为模拟时间；g＝9.81m/s²为重力加速度；h₀为土体颗粒物的初始高度；l_s表示t时刻坍塌体的前缘位置；m和K分别为系统的总质量和总动能。

如图3-1、图3-2比较了试验和本发明得到的无量纲坍塌前缘位置随无量纲时间的变化规律，同时给出了本发明得到的无量纲系统动能的变化规律。可以看出本发明所得坍塌无量纲前缘位置的变化规律与试验结果一致，说明本发明可以准确模拟泥石流的起动形成、流动发展和再次淤积的整个运动过程，验证了固液转化本构模型的有效性。图4可以看出：坍塌过程中可以分为起动形成、加速滑动、减速滑动和再次淤积四个阶段。在开始阶段，由于自重作用，土体开始坍塌，进而不断加速，约在T^*＝2时刻，系统动能达到极值；随着颗粒物的运动，系统构形不断改变，自重与人工阻尼力逐渐达到平衡，系统逐渐减速进入平面滑动，并最终停止运动，系统动能保持为零，坍塌前缘位置保持不变，泥石流再次淤积。本发明所得各个阶段的系统速度云图如图5所示，可以看出本发明所得速度云图较为光滑，说明本发明采用的高阶次B样条基函数(如图2所示)消除了物质点法的网格穿越误差。

实施例2

本实施例为膨润土沿斜面的自由滑动模型正视图，如图6所示，其中：h₀＝0.0472m，l₀＝0.1888m，θ＝15°，z方向厚度为0.34m。在本发明中，通过分解重力加速度来模拟土体沿斜面的下滑，即：分别沿着斜面和垂直于斜面建立直角坐标系，将重力加速度g沿斜面分解成g_x＝gsinθ和g_y＝gsinθ。

如图7给出了无量纲化的滑动的前缘位置和动能随无量纲时间的变化规律，可以看出：本发明所得结果与试验结果一致，说明本发明能够准确模拟泥石流的起动形成、流动发展和再次淤积的整个运动过程，再次验证了固液转化本构模型的有效性。可以看出：滑动过程中可以分为起动形成、加速滑动、减速滑动和再次淤积四个阶段。在初始滑动阶段中，由于自重作用，系统动能不断增加，约在T^*＝6时刻，达到极值；随后，系统进入减速滑动，系统动能逐渐减小；最终，土体在斜面上缓慢滑动，自重与人工阻尼力逐渐达到平衡，系统动能几乎保持为常数，泥石流沿斜坡缓慢移动。图8给出了模拟所得不同时刻系统的速度云图，可以看出本发明所得速度云图较为光滑，说明本发明采用的高阶次B样条基函数(如图2所示)消除了物质点法的网格穿越误差。

以上所述仅为本发明的实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.模拟泥石流起动形成、流动发展和再次淤积的数值方法，其特征在于，该方法包括以下步骤：

步骤1、将各物质点上的质量和动量通过高阶B样条基函数映射到相应的背景网格，计算背景网格节点处的质量和动量：

式中：上标k表示计算时间步；下标J和p分别为背景网格节点和物质点上的变量；S为高阶B样条基函数；m为质量；x和

分别为位置和速度矢量；

步骤2、分别计算背景网格节点上的内力f^int和外力f^ext：

式中：

表示高阶B样条基函数的梯度；σ为柯西应力张量；

和b分别表示面力和体力矢量；h表示面力边界厚度；V为体积；

步骤3、引入人工阻尼力，在背景网格节点上进行动量方程的求解，更新背景网格节点的加速度和速度：

式中：

和

分别为速度和加速度矢量；f^damp为人工粘性力；

其中，所述人工阻尼力f^damp为：

f^damp＝-ζ_d·m_p·v_p (7)

式中：m_p和v_p分别表示物质点p的质量和速度；ζ_d为人工阻尼系数，求解公式为：

其中：ξ为无量纲的阻尼系数,h为背景网格尺寸，对一维、二维和三维问题，分别取h＝dx,

和

E和ρ分别为物质点的弹性模量和密度；

步骤4、将更新后的背景网格节点的速度和加速度分别映射回相应的物质点，更新各物质点的速度和位置：

步骤5、将更新后的物质点速度再次映射回相应的背景网格节点，得背景网格节点速度

步骤6、利用背景网格节点速度

计算各物质点的变形率张量和应变增量：

其中：

和Δε_p分别为变形率张量和应变增量；

步骤7、更新物质点的密度：

(ρ_p)^k+1＝(ρ_p)^k/[1+(Δε_p)^k+1] (13)

步骤8、基于固液转化本构模型更新各物质点的应力；

a)若物质点的密度小于密度阈值时则视为流体状态，其中密度阈值一般取1800kg/m³～2000kg/m³，即为泥水混合物的密度，步骤8采用基于非牛顿流体本构模型更新其应力，其更新过程如下：

基于引入广义Cross非牛顿流体本构模型以描述泥石流的非牛顿流体流动，其等效粘度系数写为：

式中：μ₀和μ_∞分别表示流体在低剪切率和高剪切率下的粘度系数；K和m为常系数；利用等效粘度系数μ_eff，可将广义Cross模型转化为牛顿流体模型形式：

得总的应力张量为：σ_ij＝-Pδ_ij+τ_ij，其中，基于人工状态方程求解流体的静水压力P：

式中：c为人工声速；ρ₀和ρ分别为流体的初始密度和当前密度；r为常数，一般取为7；

b)若物质点密度大于等于密度阈值则视为固体状态，步骤8采用基于Drucker-Prager弹塑性本构模型更新其应力，其更新过程如下：

式中：I₁＝σ_kk为应力张量的第一不变量；J₂＝s_ijs_ij/2为应力偏量的第二不变量；q_φ为摩擦系数，反映了压力对屈服极限的影响程度；k_φ为纯剪切状态下时的屈服应力；q_φ和k_φ可由材料的内聚力c和摩擦角φ确定；通过选取参数：

和

则Drucker-Prager屈服面可在π平面上外接Mohr-Coulomb屈服面，在Drucker-Prager模型中，剪切势函数采用非关联塑性流动法则，即：ψ^s＝τ+q_ψσ_m

式中：材料常数q_ψ由剪胀角ψ确定；q_ψ与ψ之间的关系和q_φ与φ之间的关系相同；如果q_ψ＝q_φ，则为关联塑性流动；如果q_ψ＝0剪胀角为零，材料满足塑性不可压缩条件；

步骤9、输出所需物质点的相关物质信息。

Images