企业官网建设流程全解析

简介：该PDF文档围绕离散元方法在旋转筒二元颗粒混合模拟中的具体应用展开，适合机械、化工、制药及建材等领域从事颗粒设备设计与工艺优化的工程师、研究者，以及学习EDEM软件和离散元建模的高年级本科生或研究生。资源为单份PDF，共1个文件，压缩包仅32KB，内容紧凑、便于快速阅读。目前已有129人学习，具有一定参考价值。文档从离散元接触模型、EDEM软件模拟流程、单元空间划分到蒙特卡罗统计拟合均作了系统说明，并重点给出了抄板样式、长度、角度、数量及填充率等参数对混合效果的定量影响。关键结论包括：抄板长度与筒半径比约0.4、抄板与筒壁间距与筒半径比约0.2、装填量约40%、抄板数量4~5个、抄板角度约75度时混合效率最优。这份资料兼具理论与工程参数依据，可直接用于同类混合设备的模拟方案设计与优化参考。

1. 离散元方法在颗粒混合模拟中的定位

颗粒混合看起来简单：把几种粉体倒进容器，搅一搅，混合均匀就完事。但实际生产中，片剂含量均匀度超标、电池正极浆料分散不良、3D打印粉体回料结块，这些问题的根源往往都是混合器内部颗粒运动状态不明。传统实验只能测混合终点的取样结果，取不到过程信息，更看不到颗粒尺度上的偏析、团聚和破碎。离散元方法（Discrete Element Method，简称DEM）的价值就在这里——它把每个颗粒当作独立刚体，逐个计算碰撞和运动，从机理层面复现混合过程。

这篇博文要讲的，是基于DEM做颗粒混合数值模拟分析时，从建模、标定到后处理的完整路径。面向的读者是已经接触过仿真工具、但不想停留在软件操作层面的工程师。你会在文中看到接触模型怎么选、时间步长怎么定、混合度用什么指标量化，也会看到最常见的数值发散和颗粒穿透问题在哪里发生。这些内容不依赖某个特定商业软件，思路在LIGGGHTS、EDEM、Rocky里都适用，自己写代码的也能参考。

2. 离散元方法的计算框架与接触力模型

2.1 颗粒运动方程的离散化写法

DEM的核心是把颗粒运动拆成平动和转动两部分。对第 i 个颗粒，平动方程是：

m_i * dv_i/dt = m_i * g + Σ F_c + F_cohesive

转动方程是：

I_i * dω_i/dt = Σ (r_i × F_c + M_r)

这里的 F_c 是颗粒间或颗粒与壁面间的接触力，F_cohesive 是额外施加的粘聚力（比如液体桥力或范德华力），M_r 是滚动摩擦产生的阻力矩。数值求解时常用中心差分格式，也就是速度 Verlet 或 leapfrog 积分。时间步长内假设力不变，更新速度，再更新位置。颗粒混合模拟的时间尺度通常在秒级，而单步临界时间步长在微秒级，一个搅拌混合案例动辄计算数千万步，所以计算效率直接决定了模拟方案的可行性。

接触力的计算是DEM计算量的主要部分。每个时间步都需要遍历所有可能接触的颗粒对，用邻居列表（Verlet List）配合网格搜索加速。常见的做法是把计算域划分成规则网格，颗粒半径大于网格边长，一个颗粒只可能和本网格及相邻网格内的颗粒接触，这样就把 O(N²) 的全局配对降到了 O(N) 量级。N 是颗粒总数。颗粒数超过十万以后，网格尺寸和更新频率的选择会明显影响速度，具体后面讲参数时一起说。

2.2 Hertz-Mindlin 接触模型的本构细节

颗粒混合中应用最多的接触模型是 Hertz-Mindlin 无滑移模型。法向力按 Hertz 接触理论计算，切向力按 Mindlin 增量形式累加。法向力表达式为：

F_n = (4/3) * E_eff * sqrt(R_eff) * δ_n^(3/2)

其中 E_eff 是等效弹性模量，由两个接触颗粒的材料参数决定：

1/E_eff = (1-ν_1²)/E_1 + (1-ν_2²)/E_2

R_eff 是等效接触半径，δ_n 是法向重叠量。注意：这个公式假定接触区域远小于颗粒尺寸，所以对细粉（粒径小于10微米）需要特殊处理，比如引入 JKR 或 DMT 粘聚接触模型。法向阻尼力是：

F_n_damp = -2 * sqrt(5/6) * β * sqrt(S_n * m_eff) * v_n_rel

β 由恢复系数 e 决定，S_n 是法向刚度，m_eff 是等效质量。恢复系数是标定过程中最容易测、也最容易测错的参数——很多初学者直接从文献里抄值，但实际仓库物料的回弹特性和文献数据可能差异很大，后面会专门讲标定。

切向力用增量形式更新：

F_t = F_t_prev + k_t * Δδ_t + c_t * Δv_t

并且受到库仑摩擦极限约束：F_t ≤ μ * F_n。滚动摩擦用定向恒定力矩模型或 Elasto-Plastic 模型实现，推荐后者，它在低应力下表现为弹性，高应力下塑性耗能，更符合颗粒混合中颗粒缓慢滚动的情况。

2.3 时间步长的理论估计与安全系数

时间步长决定了模拟的稳定性。DEM 中常用 Rayleigh 波在颗粒表面的传播时间来估计临界步长：

Δt_Rayleigh = π * R_min * sqrt(ρ / G) / (0.1631 * ν + 0.8766)

R_min 是颗粒最小半径，ρ 是颗粒密度，G 是剪切模量，ν 是泊松比。实际取Δt = 0.2 * Δt_Rayleigh 作为默认值。但颗粒混合中还有另一个约束：接触刚度越大，力传递越快，需要的步长越小。如果颗粒弹性模量很高（比如玻璃珠），Rayleigh 步长可能被压得极低，这时有两种处理思路：一是给颗粒设定一个"虚拟刚度"，只要不显著影响碰撞动力学行为（比如最大重叠量控制在颗粒半径的0.1%以内），允许比真实材料低几个数量级；二是改用显式动态松弛配合局部时间步长，但这实现复杂，多数场景不需要。

一个实用经验：设完模型后先跑100步，输出最大重叠量和系统动能，如果最大重叠量超过颗粒半径的0.5%，调小时间步长或降低颗粒刚度。如果动能呈指数增长，几乎可以确定是时间步长过大或者背后有接触穿透未被检测到。

3. 颗粒混合模拟的三步建模流程与实际参数设定

3.1 从几何模型到颗粒填充的完整链路

做颗粒混合模拟，第一个要解决的是"初始堆怎么放"。常用做法是重力沉积法：在混合器几何内按随机位置生成颗粒，让它们在重力作用下自由下落并堆积到平衡态。几何模型用 STL 或网格文件导入，颗粒工厂（Particle Factory）按指定的生成速率和初始速度在设定区域产生颗粒。关键点是生成速率不能太快，否则颗粒瞬间堆积导致大面积接触重叠，后续要花大量的迭代步来释放这些应变能。

一个更省时间的替代方案是分层生成：先让颗粒在混合器外的管道里沉积成柱，然后把柱体移至混合器内部。这个方法的优势是颗粒之间的初始接触相互独立，形成稳定的自然堆积结构，但需要额外的位移脚本。对双锥混合机这类回转设备，更常见的做法是直接在水平位置生成颗粒，然后开始转动混合器，让颗粒在运动过程中自行重新分布。这个方法会在前几百步产生短暂的能量尖峰，但数值稳定性通常可控。

3.2 接触参数的标定方法：斜面法与休止角对比

接触参数是影响混合模拟结果最敏感的因素，也不需要盲目追求精确测量。对于工程级模拟，通常需要标定以下六个参数。

最常用也最经济的标定实验是休止角实验：让颗粒从漏斗自由落体堆积成锥，测量堆积角。然后调整 μ_s 和 μ_r，使DEM模拟得到的堆积角与实验一致。注意休止角对 μ_s 的敏感性远高于 μ_r，而混合过程中的流动行为同时受两个摩擦系数影响，所以单独靠休止角标定不够。建议增加一个转鼓实验：颗粒在半充满的水平圆筒内转动，观察表面流动模式和拍动频率，用转鼓实验的结果同时验证 μ_s 和 μ_r。

有一个常见误用：直接把文献中的玻璃珠参数套用到活性药物成分（API）颗粒。API颗粒通常存在内聚力，休止角会比玻璃珠大得多。碰到这种情况，需要在接触模型中引入 JKR 粘聚力或简化线性内聚模型，内聚强度可以通过剪切盒实验估算。

3.3 一个可直接修改运行的DEM编程框架

不用商业软件，用Python也能跑通一个小规模的颗粒混合DEM，适合验证模型和算法。下面是一个极简的二维混合器DEM实现框架，核心逻辑是接触检测、力计算、位置更新三步循环。

import numpy as np class Particle: def __init__(self, x, y, r, m, vx=0.0, vy=0.0): self.x = x self.y = y self.r = r self.m = m self.vx = vx self.vy = vy self.fx = 0.0 self.fy = 0.0 def compute_contact(p1, p2, kn, kt, mu): dx = p2.x - p1.x dy = p2.y - p1.y dist = np.sqrt(dx*dx + dy*dy) overlap = p1.r + p2.r - dist if overlap <= 0: return nx = dx / dist ny = dy / dist # 法向弹性力（线弹性模型） fn = kn * overlap p1.fx += fn * nx p1.fy += fn * ny p2.fx -= fn * nx p2.fy -= fn * ny # 切向摩擦力 tx = -ny ty = nx rel_vx = p2.vx - p1.vx rel_vy = p2.vy - p1.vy vt = rel_vx * tx + rel_vy * ty ft = -kt * vt max_ft = mu * fn if ft > max_ft: ft = max_ft elif ft < -max_ft: ft = -max_ft p1.fx += ft * tx p1.fy += ft * ty p2.fx -= ft * tx p2.fy -= ft * ty

关键参数 kn 和 kt 就是法向和切向刚度。二维线弹性模型虽然比 Hertz 模型粗糙，但算法骨架完全一致。小规模测试时颗粒数控制在几千个，时间步长取1e-5秒，同时给边界加上重力即可输出颗粒轨迹。这个框架的价值在于你能直接看到每一步哪些力被计算了、哪些被忽略了。比如上面的代码忽略了阻尼项，实际运行一段时间后颗粒会持续振荡，需要加入阻尼力。

对于生产级模拟，推荐使用 LIGGGHTS 这样的开源 DEM 软件。输入脚本大致是这个结构：

# 定义单元类型和颗粒材料 units si atom_style granular boundary m m m # 定义材料参数 fix m1 all property/global youngsModulus 1e7 poissonsRatio 0.3 \ coefficientRestitution 0.5 coefficientFriction 0.4 \ coefficientRollingFriction 0.05 # 生成颗粒 region box block -0.05 0.05 -0.05 0.05 -0.05 0.05 units box create_box 1 box create_particles 1 radius 0.001 distribution uniform \ region box n 5000 # 重力与时间步 fix grav all gravity 9.81 timestep 0.0000001 run 100000

这段脚本创建5000个粒径1mm的颗粒，设定了恢复系数0.5、滑动摩擦0.4、滚动摩擦0.05，时间步长1e-7秒。核心逻辑是：fix命令定义材料属性，create_particles生成颗粒，timestep必须小于 Rayleigh 步长，run指定步数。实际混合模拟时，颗粒数通常要达到10万到百万量级，这时单核计算会慢到不可接受，需要启用多核并行：

processor 8 * * comm_style tiled

8表示在x方向划分8个子域。LIGGGHTS按空间域分配进程，颗粒在域间迁移时自动通信。颗粒混合中大量颗粒会聚集在混合器底部，均匀划分域会造成负载不均，可以用fix balance命令做动态负载均衡。

3.4 转速、填充率和叶片角的初始设定

混合操作的工艺条件直接影响模拟结果。转速用无量纲的Froude数描述更通用：

Fr = ω²R/g

ω是角速度，R是混合器半径，g是重力加速度。逆流混合器通常Fr在0.1到2之间。填充率（装填体积占混合器容积的百分比）推荐从30%到70%扫描。填充率太低，颗粒抛落高度不足，混合慢；填充率太高，颗粒层内部缺乏流动，混合不均匀。叶片角对剪切混合影响显著，30度桨叶产生强剪切，60度偏重对流。做研究时至少扫三个角度才能看出趋势。

模拟运行结束后的输出控制也值得重视。每1000步输出一次颗粒坐标，生成VTK格式供ParaView可视化。每100步记录一次系统中的动能和各颗粒坐标，用于计算混合速率。计算混合速率需要知道每条颗粒轨迹上的局部颗粒种类占比，所以输出频率不能太低，否则时间平均后信息丢失。

4. 从模拟轨迹到混合程度的定量评价

4.1 混合度的定义：Lacey混合指数及其前提

混合程度用Lacey混合指数M表示是一个经典做法。它的定义需要三个浓度方差。第一个是取样得到的实际方差S²，第二个是理论上完全分离状态的方差S₀²，第三个是理论上完全混合状态的方差Sᵣ²。

S₀² = p(1-p)

Sᵣ² = p(1-p)/N_p

p是目标组分占总体积的分数，N_p是取样区域内的颗粒数量。Lacey指数定义为：

M = (S₀² - S²) / (S₀² - Sᵣ²)

M=0表示完全未混合，M=1表示完全混合。实际过程中M从0迅速上升到0.9左右，随后接近渐近线。混合速率通常把M对时间曲线拟合成指数衰减模型：M(t) = M_max - A*exp(-kt)，k 就是混合速率常数，也是不同工况对比的核心指标。

计算Lacey指数的前提是取样区域大小要合理。取样区域太小，颗粒数量少，Sᵣ²大；取样区域太大，空间分辨率不足。工程惯例是把混合器划分成若干立方体单元，每个单元内包含至少50到100个颗粒，这样统计噪声可控。颗粒粒径分布越宽，取样区域的边长至少要高于最大粒径的5倍，否则边界效应会明显污染统计结果。

4.2 用Python处理DEM轨迹数据并计算Lacey指数

下面是处理DEM输出数据的示例脚本，使用Python处理LIGGGHTS的dump文件、读入颗粒坐标和类型，然后划分网格并计算Lacey指数。

import numpy as np def compute_lacey(data, p_label, grid_size, box_limits): """ data: numpy数组,每行 [x, y, z, type] p_label: 目标组分类型(0或1) grid_size: 取样单元边长 box_limits: 计算域边界 [[xmin,xmax],[ymin,ymax],[zmin,zmax]] """ p_total = 0 variances = [] xmin, xmax = box_limits[0] ymin, ymax = box_limits[1] zmin, zmax = box_limits[2] # 各取样子单元内的颗粒数量 nx = int((xmax - xmin) / grid_size) ny = int((ymax - ymin) / grid_size) nz = int((zmax - zmin) / grid_size) concentrations = [] for i in range(nx): for j in range(ny): for k in range(nz): x_range = (xmin + i*grid_size, xmin + (i+1)*grid_size) y_range = (ymin + j*grid_size, ymin + (j+1)*grid_size) z_range = (zmin + k*grid_size, zmin + (k+1)*grid_size) mask = ((data[:,0] >= x_range[0]) & (data[:,0] < x_range[1]) & (data[:,1] >= y_range[0]) & (data[:,1] < y_range[1]) & (data[:,2] >= z_range[0]) & (data[:,2] < z_range[1])) cell_data = data[mask] if len(cell_data) < 10: # 忽略颗粒过少的子单元 continue n_target = np.sum(cell_data[:,3] == p_label) p_local = n_target / len(cell_data) concentrations.append(p_local) p_overall = np.mean(concentrations) S2 = np.var(concentrations) S0_2 = p_overall * (1 - p_overall) N_avg = np.mean([len(data[(data[:,0] >= xmin) & (data[:,0] < xmax) & (data[:,1] >= ymin) & (data[:,1] < ymax) & (data[:,2] >= zmin) & (data[:,2] < zmax)]) / len(concentrations)]) Sr2 = p_overall * (1 - p_overall) / N_avg lacey = (S0_2 - S2) / (S0_2 - Sr2) return lacey

代码逻辑分三步：第一步按 grid_size 划分空间网格，第二步统计每个网格中目标组分的局部浓度，第三步按 Lacey 公式计算方差比。注意len(cell_data) < 10这个阈值：取样单元内颗粒太少时局部浓度的泊松噪声会主导方差，导致计算结果不可信。这个阈值通常取50，用户可以根据颗粒总数调。

4.3 颗粒尺度的偏析指数：径向分布函数与配位数

Lacey指数反映整体均匀程度，但颗粒尺度的微观混匀状态需要用径向分布函数（RDF）检查。RDF 描述的是距某个颗粒距离 r 附近出现颗粒的概率密度相对均匀分布时的比值。混合模拟中对比不同种类颗粒RDF的差异，可以揭示偏析趋势。

计算RDF时，对每一对颗粒计算距离，按直方图统计落入各距离区间的对数。需要特别注意：RDF对有限尺寸效应很敏感，系统边界附近的颗粒对距离分布会产生偏差。解决办法是施加周期性边界条件：如果颗粒i在边界附近，在镜像域中复制颗粒i的邻居，然后计算距离。这个方法在LIGGGHTS中正好对应boundary p p p的设置。

配位数（每个颗粒平均接触数量）是另一个微观指标。配位数高的区域通常对应压实的颗粒堆积区，直接影响混合器内的应力分布。用Python统计配位数时，只需要遍历dump文件中坐标间距小于两个颗粒半径之和的颗粒对。这个操作和DEM主循环里的邻居搜索一样，暴力循环在颗粒数较多时不可行，需要网格加速。

配位数受时间影响强烈：混合器转动时，颗粒受到离心力，会在器壁附近形成高配位数区域；中心区域配位数较低。对比不同转速下的配位数分布，能直观看到离心混合区如何形成和扩展。

5. 常见数值病态与高级技巧：刚度过高、颗粒穿透、多时间尺度处理

5.1 颗粒穿透问题的三个直接原因

模拟中出现颗粒互相嵌入并且不弹开，是最常见的问题。这个现象通常由三个原因导致。第一个是时间步长过大：颗粒在一个时间步内移动的距离超过重叠量，接触检测错过碰撞。对粒径1mm、刚度1e7 Pa的颗粒，速度1m/s时，一个时间步1e-5秒内移动10微米，看似很小，但接近重叠量的量级，不安全；要缩小到1e-6秒。

第二个原因是接触搜索邻居列表更新频率太低。LIGGGHTS 中neigh_modify every 1 delay 0 check yes设置每个时间步更新邻居列表。颗粒混合中颗粒运动剧烈，如果每10步才更新一次，高速颗粒可能穿越邻居列表的皮肤厚度，直接落入"无接触"状态。表现就是颗粒穿透后，在下一次列表更新前已经移动到很深的位置。

第三个原因是初始重叠过大后力没有释放时间。比如颗粒工厂一次性生成大量颗粒，初始堆叠时的嵌入深度超过半径的1%。此时即使时间步长正常，接触力也会造成数值振荡。应对方法：把颗粒生成速率调低，让颗粒逐一降落并稳定后再生成下一批。

5.2 混合模拟的并行效率优化与负载均衡

颗粒混合模拟的并行效率受颗粒空间分布影响很大。混合器内的颗粒在重力作用下大部分时间沉积在底部，空间域分解后会有一个进程承担大量颗粒，其他进程空闲。启动动态负载均衡：

fix balance all rebalance 100 1.1 shift xy 100 1.2

这个命令让LIGGGHTS每100步检查一次进程间颗粒数差异，超过1.1倍则触发重新划分，直到差异小于1.2倍。注意shift只在指定方向迁移域边界，混合器中颗粒分布梯度最大的是竖直方向。另外，固定壁面区域的网格搜索开销也不小，合理做法是给壁面网格增加一层虚拟颗粒，而不是让所有进程都保留壁面网格副本。

颗粒总数超过500万以后，时间步长反而成了瓶颈。此时需要更多关注 GPU 加速的 DEM 求解器，比如在统一计算设备架构（CUDA）环境下实现接触力计算。GPU端把邻居列表构建和力计算放在kernel里，每帧传递颗粒位置和力数组。颗粒数100万时，GPU加速比单核CPU快20倍以上，但需要OpenGL传递数据的开销低于一定比例才划算。

5.3 非球形颗粒的替代方案：多球颗粒与超椭球模型

颗粒混合模拟如果只做球形颗粒，结果和真实粉体的偏差会体现在休止角偏小、孔隙率偏高、混合速率偏快。处理非球形有两条路线。一条是多球颗粒（Multi-sphere）：用一个刚体把多个球粘在一起，在LIGGGHTS中用球簇单元生成，比如药片形状用两个球簇拼圆柱。这个方法实现简单、计算快，缺点是表面有台阶状凸起，接触面积虚高。另一条是超椭球模型，在开源库里实现不多，多数商业软件原生支持，公式复杂但形态更准确。

工程上的务实建议是：球形颗粒结果和实验差异不超过15%，优先用球。如果对混合速率精度有硬性要求，至少把颗粒替换成双球或三球簇，再标定一次休止角。

5.4 验证模拟结果的可用技巧

模型建完后，第一步验证做颗粒堆积模拟：把颗粒填充到混合器中，测量堆积高度和孔隙率和实际设备作对比。如果堆积高度偏差超过5%，检查颗粒刚度和摩擦系数标定是否准确。第二步做混合进程的"时间-浓度"曲线和实验取样的对比，至少对比三个时间点。

一个经常被忽略的验证技巧是模拟中采样位置要和实验取样枪尺寸一致。实验取样枪直径20毫米，那么模拟中的采样球半径也取10毫米，统计该球范围内的组分占比。取样位置选在混合器两端和中间，模拟结果和实验的方差逐点对比。通过对比，能进一步判定摩擦系数的细微差异需要调整，还是混合动力学模型本身需要修正。

本文还有配套的精品资源，点击获取