摘要：数学有限单元法原理的扩展-离散元法（Discrete Element Method, DEM），是一种通过模拟单个离散粒子的运动和相互作用来研究物质行为的数值计算方法。尤其适用于处理颗粒流体动力学问题，如颗粒物质的堆积、流动以及与容器壁的相互作用。在离散元法的
数学有限单元法原理的扩展-离散元法（Discrete Element Method, DEM），是一种通过模拟单个离散粒子的运动和相互作用来研究物质行为的数值计算方法。尤其适用于处理颗粒流体动力学问题，如颗粒物质的堆积、流动以及与容器壁的相互作用。在离散元法的数学模型中，最基本的组成部分是描述粒子运动的动力学方程。方程可以通过牛顿第二定律来描述，即对于每一个粒子，其加速度等于作用于它的所有外力与该粒子质量的比值。也就是用：初始条件、计算外力、应用牛顿第二定律、更新粒子速度、计算粒子位置、……、循环直至收敛。然后就是运动方程问题、……、等等。来模拟其过程。下面就用Python语言给个例子来看一看。程序涉及颗粒材料堆积过程中参数变化对稳定性的影响。
先看看编程所需相关的知识要点：
物理知识要点
1.重力作用：在地球上，任何物体都会受到重力的影响。在该程序里，重力加速度 gravity = 9.81 m/s² ，每一个时间步长 dt ，颗粒在垂直方向上的速度 vy 会因重力作用而减小，公式为 vy -= gravity * dt 。
2.颗粒间碰撞：当两个颗粒间的距离小于两倍的颗粒半径时，就会发生碰撞。在碰撞过程中，需要考虑法向力和切向力。
o法向力：采用弹簧 - 阻尼模型来计算，公式为 Fn = 100 * overlap - 0.1 * np.dot(v_rel, normal) 。这里的 overlap 是颗粒的重叠量， v_rel 是相对速度， normal 是单位法向量。
o切向力：依据库仑摩擦定律来计算，公式为 Ft = -params.mu * abs(Fn) * vt_rel / abs(vt_rel) 。其中 params.mu 是摩擦系数， Fn 是法向力， vt_rel 是相对切向速度。
3.坡面边界处理：当颗粒的位置低于坡面时，就会和坡面发生碰撞。碰撞时要考虑法向反弹和切向摩擦。
o法向反弹：法向速度依据恢复系数 params.restitution 进行更新，公式为 vn_new = -params.restitution * vn 。
o切向摩擦：切向速度按照库仑摩擦定律更新，当切向速度大小大于 0 时， vt_new = vt - np.minimum(params.mu * np.linalg.norm(vn_new), vt_mag) * vt / vt_mag 。
数学知识要点
1.向量运算：在程序中，向量运算被广泛运用，例如计算相对速度、单位法向量、切向向量等。
o相对速度：v_rel = np.array([vx[j]-vx[i], vy[j]-vy[i]]) 。
o单位法向量：normal = np.array([dx, dy]) / distance 。
o切向向量：tangent = np.array([-normal[1], normal[0]]) 。
2.三角函数：在计算坡面、滑移面以及力的分解时会用到三角函数。
o坡面计算：slope_tan = np.tan(params.slope_angle) ，slope_y = np.tan(params.slope_angle) * slope_x 。
o滑移面生成：x_slip = r * np.cos(theta) ，y_slip = r * np.sin(theta) ，其中 r = 3 * np.exp(params.mu * theta) 。
3.距离计算：运用 scipy.spatial.distance.cdist 函数来计算颗粒间的距离矩阵，从而判断颗粒是否发生碰撞。
4.数值积分：通过时间步长 dt 对颗粒的速度和位置进行更新，采用的是简单的欧拉法。
o速度更新：vx += (Fx / m) * dt ，vy += (Fy / m) * dt （程序中部分力的计算已经包含了质量的影响）。
o位置更新：x += vx * dt ，y += vy * dt 。
再看看最终的程序代码：
#SlipLineLooseLoessDEM.py
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider, Button
from scipy.spatial.distance import cdist
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
matplotlib.rcParams['axes.unicode_minus'] = False
num_particles = 200
particle_radius = 0.1
dt = 0.01
gravity = 9.81
np.random.seed(42)
x = np.random.uniform(0, 5, num_particles)
y = np.random.uniform(3, 6, num_particles)
vx = np.zeros(num_particles)
vy = np.zeros(num_particles)
fig, ax = plt.subplots(figsize=(10, 6))
plt.subplots_adjust(left=0.1, right=0.75, top=0.9, bottom=0.3)
scat = ax.scatter(x, y, s=20, c='gold', edgecolor='sienna')
slope_line, = ax.plot(, , 'r--', lw=1.5)
slip_surface, = ax.plot(, , 'b-', lw=2)
ax_phi = plt.axes([0.2, 0.15, 0.6, 0.03])
phi_slider = Slider(ax_phi, '内摩擦角 φ (°)', 20, 45, valinit=30)
ax_angle = plt.axes([0.2, 0.1, 0.6, 0.03])
angle_slider = Slider(ax_angle, '坡面角度 (°)', 15, 60, valinit=35)
reset_button_ax = plt.axes([0.8, 0.15, 0.1, 0.04])
reset_button = Button(reset_button_ax, '重置')
class DEMParams:
def __init__(self):
self.mu = np.tan(np.radians(30))
self.restitution = 0.3
self.slope_angle = np.radians(35)
params = DEMParams
def update_params(val):
params.mu = np.tan(np.radians(phi_slider.val))
params.slope_angle = np.radians(angle_slider.val)
phi_slider.on_changed(update_params)
angle_slider.on_changed(update_params)
def reset_simulation(event):
global x, y, vx, vy
reset_button.on_clicked(reset_simulation)
def dem_step:
global x, y, vx, vy
positions = np.column_stack((x, y))
dist_matrix = cdist(positions, positions)
np.fill_diagonal(dist_matrix, np.inf)
collision_pairs = np.argwhere(dist_matrix
for i, j in collision_pairs:
dx = x[j] - x[i]
dy = y[j] - y[i]
distance = np.hypot(dx, dy) + 1e-8
overlap = 2*particle_radius - distance
if overlap > 0:
normal = np.array([dx, dy]) / distance
v_rel = np.array([vx[j]-vx[i], vy[j]-vy[i]])
Fn = 100 * overlap - 0.1 * np.dot(v_rel, normal)
tangent = np.array([-normal[1], normal[0]])
vt_rel = np.dot(v_rel, tangent)
if abs(vt_rel)
Ft = 0.0
else:
Ft = -params.mu * abs(Fn) * vt_rel / abs(vt_rel)
vx[i] -= (Fn*normal[0] + Ft*tangent[0]) * dt
vx[j] += (Fn*normal[0] + Ft*tangent[0]) * dt
vy[i] -= (Fn*normal[1] + Ft*tangent[1]) * dt
vy[j] += (Fn*normal[1] + Ft*tangent[1]) * dt
slope_tan = np.tan(params.slope_angle)
for i in range(num_particles):
if y[i]
normal = np.array([slope_tan, -1])
normal /= np.linalg.norm(normal)
v_rel = np.array([vx[i], vy[i]])
vn_mag = np.dot(v_rel, normal)
vn = vn_mag * normal
vt = v_rel - vn
vn_new = -params.restitution * vn
vt_mag = np.linalg.norm(vt)
if vt_mag > 0:
vt_new = vt - np.minimum(params.mu * np.linalg.norm(vn_new), vt_mag) * vt / vt_mag
else:
vt_new = vt
vx[i], vy[i] = vn_new + vt_new
vy -= gravity * dt
x += vx * dt
y += vy * dt
def update_slip_surface:
theta = np.linspace(0, np.pi/2, 50)
r = 3 * np.exp(params.mu * theta)
return r * np.cos(theta), r * np.sin(theta)
def update_frame(frame):
dem_step
scat.set_offsets(np.column_stack((x, y)))
x_slip, y_slip = update_slip_surface
slip_surface.set_data(x_slip, y_slip)
slope_x = np.linspace(0, 5, 50)
slope_y = np.tan(params.slope_angle) * slope_x
slope_line.set_data(slope_x, slope_y)
return scat, slope_line, slip_surface
ax.set_xlim(0, 5)
ax.set_ylim(0, 6)
ax.set_xlabel('Horizontal Distance (m)')
ax.set_ylabel('Height (m)')
ax.set_title('Granular Material Slip Surface Simulation')
ax.grid(True, linestyle='--', alpha=0.5)
from matplotlib.animation import FuncAnimation
ani = FuncAnimation(fig, update_frame, frames=200, interval=50, blit=True)
plt.show
最后是运行瞬间的截图。
来源：苑博教育
