[发明专利]一种垃圾填埋场失稳流滑灾害数值仿真方法

摘要

本发明涉及一种移动半隐式粒子法（简称MPS）的垃圾填埋失稳流滑灾害的数值仿真方法。针对原始MPS方法存在严重的压力震荡，从三个方面对源程序进行改进，包括：核函数、Poisson压力方程以及自由表面判别。为刻画流滑灾害体的非牛顿流变特性，创新性的引入宾汉姆非牛顿流体本构关系。建立等效粘度概念，将宾汉姆流体本构转化为变粘度的牛顿流体本构，并编程实现。通过平衡压力震荡和引入宾汉姆非牛顿流体本构关系发明了改进的MPS方法。该发明能够精确仿真垃圾填埋场失稳流滑的动力过程。

主权项

一种垃圾填埋场失稳流滑灾害数值仿真方法，其特征在于具体步骤如下：1)确定垃圾填埋场的研究区域及物理参数，所述物理参数包括填埋场尺寸、垃圾体密度及强度参数；2)根据步骤1)确定的垃圾填埋场的研究区域计算粒子的初始分布；3)根据步骤2)得到的粒子的初始分布考虑粘性力项和重力项以计算粒子的临时速度；所述粘性力项采用下式(1)宾汉姆流体本构模型计算得到；$\left\{\begin{array}{cc}\mathrm{\&tau;}=\mathrm{\&eta;}\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}+{\mathrm{\&tau;}}_{min}& \mathrm{\&tau;}\&gt;{\mathrm{\&tau;}}_{\min }\\ \stackrel{\&CenterDot;}{\mathrm{\&gamma;}}=0& \mathrm{\&tau;}\&le;{\mathrm{\&tau;}}_{\min }\end{array}\right.---\left(1\right)$其中：τ为剪应力，η是粘度系数，为应变速率，τ_min为屈服剪切强度；为考虑垃圾土的特性，将摩尔‑库伦准则引入上述宾汉姆流体本构模型以界定垃圾填埋场失稳流滑的初始状态，则得到屈服剪切强度τ_min如下式(2)：τ_min＝c+P tanφ  (2)c表示垃圾土粘聚力，P表示垃圾土正应力，φ表示垃圾土内摩擦角；同时为简化数值计算，将等效粘度引入η′(如下式(3))引入宾汉姆流体本构模型得到修正后的宾汉姆流体本构模型如下式(4)；${\mathrm{\&eta;}}^{\&prime;}=\mathrm{\&eta;}+\frac{c+Ptan\mathrm{\&phi;}}{\left|\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}\right|}(1-e^{-m\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}})---\left(3\right)$$\mathrm{\&tau;}={\mathrm{\&eta;}}^{\&prime;}\stackrel{\&CenterDot;}{\mathrm{\&gamma;}}---\left(4\right)$其中m为压力增长系数；4)根据步骤3)得到的粒子的临时速度计算粒子的临时位置；5)根据步骤4)得到的粒子的临时位置计算粒子的临时粒子数密度n_i，公式如下式(5)：$n_i=\underset{j\&NotEqual;i}{\&Sigma;}w\left(\right|r_j-r_i\left|\right)---\left(5\right)$其中：rj和ri代表粒子j和i的坐标，w(r)表示核函数，r表示粒子间距离；由于传统MPS采用的核函数，当r＝0时，取值为无穷大，会导致压力出现震荡；为平衡压力震荡，采用下式(6)所述的核函数代替传统MPS采用的核函数；$w\left(r\right)=\left\{\begin{array}{cc}{(1-\frac{r}{r_e})}^3{(1+\frac{r}{r_e})}^3& (0\&le;r\&lt;r_e)\\ 0& (r\&GreaterEqual;r_e)\end{array}\right.---\left(6\right)$r_e为影响域半径；6)根据步骤5)得到的粒子的临时粒子数密度建立Poisson压力方程；采用包含一个主项和两个修正项的Poisson压力方程(下式(7))以平衡压力震荡：$\frac{1}{\mathrm{\&rho;}}{\&dtri;}^2{P}_{k+1}=\frac{n^{*}-2n^k+n^{k-1}}{n^0{\mathrm{\&Delta;t}}^2}+\frac{a}{\mathrm{\&Delta;}t}\frac{n^k-n^{k-1}}{n^0\mathrm{\&Delta;}t}+\frac{b}{{\mathrm{\&Delta;t}}^2}\frac{n^k-n^0}{n^0}---\left(7\right)$其中：ρ表示垃圾体密度，▽表示梯度算子，P_k+1表示第(k+1)步的压力，n*表示临时粒子数密度，n^k表示第k步的粒子数密度，n^k‑1表示第(k‑1)步的粒子数密度，n⁰表示初始粒子数密度，△t表示时间增量，a和b为误差项修正系数；7)根据步骤6)求解Poisson压力方程得到压力场分布；8)采用自由表面判别公式(8)确定自由表面粒子，并将自由表面粒子的压力置为零，得到修改后的压力场分布；n^*<βn⁰     (8)n*表示临时粒子数密度，n⁰表示初始粒子数密度，β为系数；9)根据步骤8)得到的修改后的压力场分布，计算压力梯度修正粒子速度及位置；10)根据步骤9)，判断是否达到预设计算时间，若达到即停止计算，否则重复步骤3)‑9)，直至达到预设计算时间。