散体材料广泛存在于自然界之中，对其运动的研究在多个行业都有着重要的应用。由于其离散的特性，散体的运动受周围不断变化的接触物质的影响，具有非连续的特征。目前对散体运动数值方法的研究尚有不足，难以体现散体复杂的运动特征及碰撞过程。发展离散块体的数值计算方法，对散体的运动规律与堆积分布特征进行研究，对提高地质灾害的监测预警及防灾减灾能力，以及改进工业农业中散体裝置设计等方面具有重要意义。

本文针对散体的运动堆积问题，建立了适用于描述大位移小变形块体运动的时间两尺度模型，提出基于潜在接触对半平面/半空间关系的精细接触检测计算方法，实现了复杂形状块体间的接触检测。随后对不同尺寸块体在不同堆积情况下的运动规律进行了试验研究，基于材料特性发展了非线性的接触本构模型，完善了接触力的计算方法，通过对比试验和模拟结果，验证了所建立计算模型的准确性和适用性。将数值方法应用于海底滑坡的研究之中，为实际工程应用提供依据。主要研究内容和结论如下：

提出了块体运动的时间两尺度模型，通过刚体运动与弹性变形这两类参数对块体的运动过程进行刻画。基于拉格朗日能量系统进行单元动力学方程的推导，并在连续-非连续单元方法的框架内使用显式方法实现方程的求解。该模型可以有效地处理大位移和大转动问题，同时准确地描述块体的局部变形。通过分别考虑块体的刚体计算时步及变形计算时步，提高了计算效率。该模型的可靠性通过块体运动变形的多个案例得到了验证。

提出了基于潜在接触对半平面/半空间关系的精细接触检测计算方法，实现复杂形状块体间的接触计算。定义了多边形的角-边潜在接触对、多面体的角-面潜在接触对和棱-棱潜在接触对。基于潜在接触对的半平面/半空间包含关系判断接触状态，对未满足半平面/半空间包含关系的接触对的类型和数量进行统计，从而确定接触类型。该算法简化了接触检测的复杂度，在确保接触检测无遗漏的情况下提高了检测效率。数值案例验证了该方法在处理任意形状的凸多边形及凸多面体接触问题时的准确性及鲁棒性。

设计了一个用于研究块体运动特征的立方体分布试验，分析不同尺寸立方体的运动规律。将四种尺寸的立方体按照三种排列方式堆积于空心圆柱体内，统计去掉周围圆柱约束后立方体失稳运动的最终分布位置。对40次试验获得的具有统计学意义的平均结果进行分析，发现立方体的分布均匀性受其堆积方式的影响，分布位置受堆积方式和立方体尺寸的影响。该试验为相关领域的试验研究提供了参考，为离散材料的运动模型提供了验证数据。

基于材料性质，引入考虑界面压力影响的非线性接触本构模型，改进了多面体接触检测方法中接触力的计算。定义了新的势函数，采用势函数法对接触的法向力进行计算。在计算接触力时，首先计算接触界面的压力，再通过压力确定接触界面的摩擦系数，从而计算摩擦力。对立方体分布试验进行了模拟，模拟结果与试验结果之间的一致性验证了该模型在计算块体运动方面的可靠性。

对天然气水合物分解诱发的海底滑坡进行模拟，分析不同天然气水合物分解长度对滑坡运动规律的影响，并借助流体程序对海面波浪进行模拟。模拟结果表明，随着分解长度增加，滑坡规模与运动距离均增加。

Discrete materials are widely found in nature, and the study of their motion has important applications in several industries. Due to its discrete property, the motion of discrete materials is influenced by the changing surrounding contact material and has a discontinuous character. The current research on numerical methods for discrete materials is not sufficient to reflect the motion characteristics and collision processes between them. To improve the monitoring and early warning of geological hazards and disaster prevention and mitigation, as well as to improve the design of installations in industry and agriculture, it is important to develop numerical methods for discrete blocks to accurately assess the movement patterns and stacking distribution characteristics of the discrete blocks.

In this paper, for the movement problem of discrete materials, a multilevel motion model to describe the movement of blocks with large displacement and small deformation is established, and a delicate contact detection calculation method based on the half-plane/half-space relationship of potential contact pairs is proposed to realize the contact detection between blocks with complex shapes. Subsequently, an experimental study on the movement regulation of blocks of different sizes under different stacking conditions is carried out, a non-linear contact constitutive model is developed based on material properties, and the calculation of contact forces is improved. The accuracy and applicability of the established model are verified by comparing the test and simulation results. The numerical method is applied to the study of submarine landslides to provide a basis for practical engineering applications. The main research contents and conclusions are as follows:

A multilevel motion model of the block is developed and the movement process of the block is characterized by two types of parameters: rigid body motion and elastic deformation. The derivation of the element dynamics equations is based on the Lagrangian equations, which are solved in an explicit way within the framework of the continuum-discontinuum element method (CDEM). The model can effectively deal with large displacements and large rotations, while accurately describing the local deformation of the block. An accurate description of the block motion can be achieved with fewer elements, which improves computational efficiency. The simulation accuracy of this model is verified by several cases.

A delicate contact detection algorithm based on the half-plane/half-space relationship of potential contact pairs is developed to realize the contact detection between complex-shaped blocks. Vertex-edge potential contact pair is defined for polygons. Vertex-face and edge-edge potential contact pairs are defined for polyhedra. The contact state is determined based on the half-plane/half-space inclusion relationship of potential contact pairs. The type and the number of contact pairs that do not satisfy the half-plane/half-space inclusion relationship are then counted to determine the contact type. The algorithm simplifies the complexity of contact detection and improves detection efficiency while ensuring that no contact is missed. Numerical cases validate the accuracy and robustness of this algorithm in dealing with arbitrarily shaped convex polygons and convex polyhedra in contact problems.

A distribution experiment for studying the movement characteristics of cubes is designed to analyze the movement patterns of cubes of different size. Four sizes of cubes are stacked inside a hollow cylinder according to three size arrangements. The distribution of cubes after removing the surrounding cylindrical constraints is counted. The statistically significant results obtained from 40 experiments are analyzed. It is found that the distribution position of cubes is influenced by their stacking method and cube size, and the uniformity of distribution is influenced by the stacking height of cubes. The experiment provides a reference for experimental studies in related fields and validation data for the motion model of discrete materials.
    
Based on the material properties, a nonlinear contact constitutive model considering the influence of interface pressure is introduced. The calculation of contact forces in the polyhedral contact detection method is improved. A new potential function is defined, and the potential function method is used to calculate the normal force of the contact. In calculating the contact force, the pressure at the contact interface is calculated first, and then the friction coefficient at the contact interface is determined from the pressure to calculate the friction force. Simulations are carried out for the cube distribution tests, and the consistency between simulation results and test results verifies the reliability of the method in calculating block motion.

The submarine landslide induced by gas hydrate dissociation is simulated to analyze the influence of different lengths of gas hydrate dissociation on the landslide movement. The sea surface waves are simulated in combination with a fluid program. The simulation results show that the landslide scale and movement distance increase with the increase of dissociation length.