一种土体大变形流动的计算模拟方法

摘要

本发明公开了一种土体大变形流动的计算模拟方法，该方法基于分子动力学（简称MD）方法，以牛顿运动方程为控制方程，引入Hertzian型摩擦公式和斯多克公式分别描述固‑固颗粒间相互作用和固‑液颗粒间相互作用，运用Verlet算法求解积分过程，建立了固‑液耦合的土体大变形流动计算模型。本发明可以呈现土体颗粒复杂流动的基本动力学特征及其运动规律，有利于进一步完善和提高我国防治土体大变形流动灾害的研究水平，同时也为分子动力学计算方法在土体大变形流动灾害中的广泛应用和推广奠定一定基础。

主权项

一种土体大变形流动的计算模拟方法，其特征在于：所述土体大变形流动的计算模拟方法包括以下步骤：1)在分子动力学数值计算模拟软件LAMMPS中输入粒子信息文件：基于计算体系特征，首先设置计算体系的维数、单位、粒子类型、粒子半径、边界条件以及积分算法，创建模拟单元；然后设置计算体系基本信息，所述计算体系基本信息包括但不限于系综、输出的时间间隔、步长以及运行步数；所述计算体系是要进行数值计算模拟的对象，所述要进行数值计算模拟的对象是土体；所述计算体系特征是要进行数值计算模拟的对象的物理力学性质及所处边界条件；2)初始设置模块启动，分子动力学数值计算模拟软件LAMMPS对计算体系的所有粒子赋予初始位置以及初始速度，进行初始化；输入模型参数、粒子初始密度及粒子所受外力；所述粒子是组成要进行数值计算模拟的对象的分子、原子以及离子；3)模拟控制模块启动，质点邻近搜索：根据粒子影响半径，运用Verlet邻近搜索法对计算粒子影响范围内所有粒子进行搜索，逐一判断周围粒子与计算粒子之间的间距与影响半径之间的关系，确定周围粒子与计算粒子之间的间距小于影响半径的周围粒子作为计算粒子的邻近粒子；直至确定每个粒子的所有邻近粒子；所述计算粒子是正在被进行计算的粒子；所述周围粒子是计算粒子周围的粒子；所述影响半径是依据要进行数值计算模拟的对象的颗粒粒径大小设置的距离；所述影响半径≤L/2，其中L是模拟单元立方体盒子的边长；4)计算每个粒子的所有邻近粒子对这些粒子所产生的作用力、瞬时加速度以及瞬时速度：根据步骤2)输入的模型参数、粒子初始密度及粒子所受外力计算步骤3)所得到的每个粒子的所有邻近粒子对这些粒子所产生的作用力，然后通过矢量求和求得每个粒子上所受合力，由牛顿运动方程求得瞬时加速度，采用Verlet积分算法对牛顿运动方程求解，得到每个粒子的瞬时速度和位置；所述每个粒子的所有邻近粒子对这些粒子所产生的作用力的计算方式是采用Hertzian型摩擦公式和斯多克定律计算得到的，所述Hertzian型摩擦公式是：$F_{hz}\left(t\right)=\sqrt{\delta }\sqrt{\frac{R_i{R}_j}{R_i+R_j}}[\left(k_n{\mathrm{\delta n}}_{ij}\right(t)-m_{eff}{\gamma }_n{v}_n(t\left)\right)-\left(k_t{\mathrm{\Delta s}}_t\right(t)+m_{eff}{\gamma }_t{v}_t(t\left)\right)]$其中：δ是第i个粒子和第j个粒子之间重叠的距离；R_i、R_j分别是第i个粒子以及第j个粒子的半径；k_n、k_t分别是第i个粒子和第j个粒子之间的法向弹性系数以及切向弹性系数；所述k_t＝2/7k_n；γ_n、γ_t分别是第i个粒子和第j个粒子之间的法向阻尼系数以及切向阻尼系数；所述γ_t＝1/2γ_n；m_eff是第i个粒子和第j个粒子之间的有效质量，所述m_eff＝m_im_j/(m_i+m_j)；m_i，m_j分别为第i个粒子和第j个粒子的质量；Δs_t是第i个粒子和第j个粒子之间的切向位移矢量；n_ij是连接第i个粒子和第j个粒子中心的单位矢量；v_n、v_t分别是第i个粒子和第j个粒子相对速度的法向分量以及切向分量；t是时间；i是第i个粒子；j是除第i个粒子外的所有粒子；所述斯多克定律的计算公式是：f＝‑3πηdv其中：η是水分子的粘度；d是第i个粒子的粒径；v是第i个粒子的运动速度；所述每个粒子的瞬时速度的计算公式是：$v_i\left(t\right)=\frac{{\mathrm{dr}}_i}{dt}=\frac{1}{2\left(\Delta t\right)}\left[r_i\left(t+\Delta t\right)-r_i\left(t-\Delta t\right)\right]+O\left({\left(\Delta t\right)}^2\right)$其中：i是第i个粒子，即正在被进行计算的粒子；t是时间；v_i(t)是第i个粒子在t时刻的速度；r_i(t)是第i个粒子在t时刻的位置；Δt是时间间隔；r_i(t‑Δt))是第i个粒子在t‑Δt时刻的位置；r_i(t+Δt)是第i个粒子在t+Δt时刻的位置；O((Δt)²)是Δt的二阶无穷小；5)重复步骤4)，对步骤4)所得到的每个粒子的所有邻近粒子对这些粒子所产生的作用力、瞬时加速度以及瞬时速度进行更新；6)实时计算结果输出模块启动，统计所需物理量，求出模拟系统中所有粒子在相空间的轨迹后,应用统计物理原理得到计算体系的宏观特性和结构特点；7)判断模拟时间是否大于设置的总时间，若模拟时间大于设置的总时间时，步骤5)的更新过程结束，程序终止；所模拟时间小于设置的总时间时，重复步骤5)直至模拟时间大于设置的总时间；所述模拟时间是步骤1)中步长与运行步数的乘积；所述设置的总时间是步骤1)中步长与运行步数总和的乘积。