简介：离散元仿真中，接触模型决定了颗粒间相互作用的准确性，而JKR等经典模型常因内聚力参数固定而难以模拟真实工况。粉体输送、湿颗粒干燥等场景下，内聚力随位置、时间或含水率动态变化，固定参数会导致仿真结果“趋势对、数据错”。通过EDEM二次开发，利用API在接触回调中动态调整内聚力调制系数，可让模型更贴合实际。本文从接触模型原理出发，讲解可变内聚力的驱动源分类、C++插件实现、Visual Studio编译配置及性能优化要点，为粉体工程和颗粒仿真提供了一套可落地的解决方案。 做了几年EDEM二次开发，我发现最常被问到的问题不是“模型怎么跑不起来”，而是“内聚力能不能让它变起来”。颗粒物料筛分、输送、造粒，几乎没有一个工业过程的内聚力是恒定值：料仓底部压实的粉体和顶部松散层，粘结行为完全是两回事；刚进入干燥段的湿颗粒和出口处的干颗粒，液桥力差了一个量级。这时候如果还抱着固定参数的JKR模型不放，仿出来的结果就是既像又不像，客户看完只说一句“趋势对了，数据对不上”。所以当时我们立项做2_Variable_Cohesion_API_EDEM_这个方案，本质思路很简单：不再去改颗粒参数表，而是直接通过EDEM的API在接触计算层把内聚力改成“变量”，让它随坐标、时间、场数据实时变化。这篇文章就记录了这个过程里我认为最有价值的部分——模型怎么设计、API怎么写、插件怎么调，以及那些文档里不会写的坑。
整个方案适合两类人看：一是刚接触EDEM二次开发、想搞清楚自定义接触模型到底能做什么的工程师；二是已经在做粉体或湿颗粒仿真、被固定内聚力参数折磨过，想给模型增加“动态感”的朋友。如果你只是拿EDEM跑跑重力卸料，那这篇文章可以收藏了以后再读；但只要你做的物料跟水、温度、粘结剂沾边，可变内聚力几乎是你迟早要迈过去的一道坎。
1. 到底要不要上API：先弄明白“内聚力变化”是什么变化
1.1 固定内聚力模型的尴尬位置
EDEM默认自带的那几个接触模型，比如Hertz-Mindlin加JKR，在工程里用得最多。JKR模型加了一项表面能参数，用来模拟颗粒之间因为有粘结剂或者液桥而产生的粘附力，湿颗粒、粉体、甚至部分秸秆类生物质物料都能用它。问题在于：JKR的表面能是一个常量，进了模型参数表之后，全程仿真都不变。颗粒从含水率20%的区域走到含水率5%的区域，它跟邻居颗粒之间的粘附力还是一模一样，这显然不物理。
更麻烦的是，真实工况里很多过程本身就是驱动物性变化的。典型例子是滚筒干燥：颗粒在滚筒入口处还是湿的，随着翻转往前运动，水分不断蒸发，颗粒间液桥越来越小，到最后几乎是干颗粒的低粘聚力状态。如果整个滚筒用一个固定JKR表面能，那只能取一个“平均”值，取低了进口段液桥力失真，取高了出口段根本不该粘在一起的地方全粘住了。这时候你调参数调得再精细，模型结构本身就压制了你的精度上限。
1.2 可变内聚力的四类常见驱动源
“让内聚力变起来”这句话听起来简单，但第一步要先定义清楚：它到底跟随什么量在变？我在实际项目里遇到过的需求基本可以归成下面四类，搞清楚这四类，后面写API的时候才知道该从哪个数据渠道取数。
| 驱动类型 | 物理含义 | 典型工业场景 | 实现难度 |
|---|---|---|---|
| 空间坐标驱动 | 内聚力随位置变化，例如料仓轴向不同高度含水量不同 | 料仓储存、竖式干燥器 | 较低 |
| 时间驱动 | 内聚力随仿真时间变化，模拟干燥、固化过程 | 滚筒干燥、固化炉 | 低 |
| 颗粒属性驱动 | 内聚力跟随颗粒自身属性，如半径、含水率标签 | 包衣造粒、粉体混合 | 中 |
| 外部场数据驱动 | 从温度场、水分场或实验数据表查表获取内聚力 | CFD-DEM耦合、数据驱动数字孪生 | 较高 |
为什么这个区分很关键？因为EDEM API的接触模型在每个接触处被调用时，你能拿到的信息是当前的接触点坐标、两个颗粒的坐标速度半径、重叠量、材料属性等，不一定能拿到你想要的“含水量”。如果你要的驱动量不在接触模型能访问的数据范围内，就得提前把数据通过自定义粒子属性或者全局属性挂进去，这个设计必须在写代码之前就定好，否则后面改起来非常痛苦。我那版2.0项目里，第一版就是没想好驱动源，写到一半发现接触点坐标能拿到但不知道颗粒含水量，被迫返工。
2. EDEM API自定义接触模型：它到底在哪个环节插了一脚
2.1 API不是用来“改软件”的，是用来“挂数据”的
很多第一次接触EDEM二次开发的人有一个误解，以为API能帮他们改求解器、改接触判定逻辑。实际上EDEM API给使用者开放的是插件化扩展点，它允许你自定义一些EDEM主程序在执行到特定环节时调用的外部逻辑。对颗粒仿真最有用的三类扩展点就是：自定义接触模型（Contact Model）、自定义粒子体力模型（Particle Body Force）、自定义仿真域逻辑。
我们的Variable Cohesion模型属于第一类——接触模型。它的调用时机非常明确：每当地质求解器检测到两个颗粒（或颗粒与几何体）发生接触，并且你选用的接触模型被激活时，API会回调你写好的函数，让你去计算这个接触产生的力。计算完了，把力返回给求解器，求解器把它跟重力、流体曳力等其他受力汇总，更新粒子的运动状态。所以你可以把自定义接触模型理解成一个“临时代办”：求解器把当前接触的现场信息交给你，你把这笔账算清楚，账目写回给求解器。
2.2 接触回调里能拿到什么数据
我以EDEM API 2022之后比较通用的接触模型接口形态来说（具体命名以你本地安装的API头文件为准，不同小版本会有改动）：接触回调函数会拿到一个接触数据对象，通过它基本上能访问到这几类信息：
- 颗粒基础数据：位置、速度、半径、角速度、质量；
- 接触几何数据：法向重叠量、切向重叠量、接触法向量、接触面积；
- 材料与属性数据：通过参数引用拿到你在EDEM界面上配置的模型参数；
- 颗粒附加属性：如果你通过API或耦合模块给颗粒挂了自定义属性，这里也能读。
这些数据里，位置和属性字段是我们做可变内聚力最常用的两条路径。比如做“随坐标变化”的内聚力，直接在回调里读颗粒坐标，取两个接触颗粒位置的均值当作接触点坐标，再映射到内聚力系数；做“随颗粒属性变化”，就给颗粒定义一个自定义属性（比如含水率标签），仿真过程中提前更新这个属性，接触模型里只负责查表。
2.3 性能问题：回调函数是热路径
有一件事必须在设计阶段就明白：接触模型回调是仿真里的热路径，一个时间步内可能被调用成千上万次。EDEM自带的高性能模型都是C++写的，而你写的自定义模型如果逻辑太胖，很短时间内就能拖垮整个仿真速度。我见过有人为了算一个“随压力变化的内聚力”，在回调里写了一个两层循环去搜索附近所有邻居，结果原本一晚上的仿真跑了两天还没到一半。
所以设计可变内聚力函数时，我给自己定了一条规矩：回调里的代码只做“查表、插值、一次性乘法”，所有需要预计算的数据（比如坐标-内聚力映射网格、数据表插值系数）都在模型初始化时算好，不要在接触回调里现场算复杂函数。这条规矩在后面几次大仿真里帮我保住了不少时间。
3. 从需求到公式：可变内聚力模型的设计细节
3.1 用一个“调制系数”去乘基准内聚力
可变内聚力模型的设计不一定要从零发明一套全新的接触本构。工程上更聪明的做法是：选取一个经过验证的基准内聚力模型（最常见的是JKR的表面能，或者简化的粘性力模型），然后在这个基准内聚力的基础上乘一个动态调制系数k。最终内聚力等于基准内聚力乘以k，k的范围通常取0到1，代表当前状态下内聚力相对基准状态保留的比例。
这个设计的最大好处是解耦了“物性基准”和“变化规律”。基准内聚力还是可以通过实验标定获得，比如通过剪切测试确定物料在某个状态下的粘结强度；而变化规律完全由k的表达式决定，改起来也很方便。如果实验数据只给了你“从湿到干内聚力衰减40%”，那你只需要把k设计成对应映射就行，不用去动底层接触力学公式，比直接改模型方程稳得多。
3.2 坐标驱动型k怎么定
以我们项目里最常用的“沿轴向变化的料仓内聚力”为例，假设实验测得了料仓高度方向上“底部压实区内聚力强、顶部松散区内聚力弱”的分布规律，并且可以用一个归一化高度z*来描述，那么k可以定义为：
k = k_min + (k_max - k_min) * f(z*)
其中f(z*)可以是一段分段线性曲线，也可以是一段样条。实际代码里我更推荐把曲线离散成一张查找表，比如每5%高度存一个k插值节点，接触回调里直接用线性插值。这样做的好处是：以后实验数据更新了，只需要改表格数据，不用重新编译插件。
同理，时间驱动型模型就是把归一化高度替换为t/T，颗粒属性驱动型则把粒子半径映射到k。万变不离其宗，核心思路是“基准力不变，调制系数变”。
3.3 力怎么叠加回去
拿到调制系数k之后，下一步是把它作用到接触力上。对JKR这类引入粘附项的法向接触模型，内聚力通常表现为法向的吸引力。在自定义接触模型里，实现方式有两种：
- 如果你的基准模型是EDEM内置模型，你可以在API里先调用内置模型得到基准法向力，自己再加一项
k放大的附加拉力； - 如果你完全自己写接触模型，那就把基准内聚力计算公式写好，乘上
k后直接写入法向力输出。
我的建议是：能调内置模型就先调内置模型，只在外面包一层调制，这样代码量小、稳定性高，也不容易因为自己写的接触力学方程缺失某个阻尼项导致仿真发散。只有当你对接触力学非常熟、并且内置模型确实无法满足需求时，才考虑完全自写。
4. C++代码实现：一个能跑的可变内聚力接触模型
4.1 工程规划和代码骨架
以下代码以EDEM API常见接口形态为例，目的是展示完整逻辑链路，具体类名和方法名请以你本地的EDEMPluginDefines.h头文件为准。不管接口名称怎么变，下面的逻辑骨架基本是通用的：注册模型、读取参数、计算调制系数、输出法向力。
#include "EDEMPluginDefines.h" #include <cmath> #include <vector> using namespace EDEM; // 查询表节点，用于坐标-调制系数插值 struct CoeffNode { double normalizedPos; double k; }; class VariableCohesionModel : public ContactModel { public: VariableCohesionModel() {} virtual ~VariableCohesionModel() {} // 初始化，在模型被加载时调用 bool initialize(PluginDataEDEM* pkg) override { // 读取在EDEM界面中配置的参数，属性名要与插件描述文件里一致 baseCohesion = pkg->getDouble("BaseCohesion"); kMin = pkg->getDouble("Kmin"); kMax = pkg->getDouble("Kmax"); // 生成坐标-调制系数查询表 // 这里用10个节点做线性插值，实际项目可以按实验数据加密 for (int i = 0; i <= 10; ++i) { CoeffNode node; node.normalizedPos = i / 10.0; node.k = kMin + (kMax - kMin) * (1.0 - node.normalizedPos); lookupTable.push_back(node); } return true; } // 接触回调，核心计算逻辑 bool calculateContact(ContactData* cdata) override { // 1. 拿到两个颗粒的位置 const ParticleData* p1 = cdata->getParticle1(); const ParticleData* p2 = cdata->getParticle2(); if (!p1 || !p2) return false; double px = 0.5 * (p1->getPosition().x + p2->getPosition().x); double pz = 0.5 * (p1->getPosition().z + p2->getPosition().z); double normZ = pz / domainHeight; // domainHeight 在初始化时从环境读取 // 2. 限制范围并查表得到调制系数 if (normZ < 0.0) normZ = 0.0; if (normZ > 1.0) normZ = 1.0; double k = interpolateK(normZ); // 3. 基准法向力基础上叠加可变内聚力 double normalForce = cdata->getNormalForce(); XYZ normal = cdata->getNormal(); double cohesionForce = k * baseCohesion; // 根据接触状态判断方向，内聚力表现为沿法线方向吸引 normalForce -= cohesionForce; cdata->setNormalForce(normalForce); // 如有需要，也可以对切向力做类似调制 return true; } private: double interpolateK(double px) { if (px <= lookupTable.front().normalizedPos) return lookupTable.front().k; if (px >= lookupTable.back().normalizedPos) return lookupTable.back().k; for (size_t i = 1; i < lookupTable.size(); ++i) { if (px <= lookupTable[i].normalizedPos) { double dx = lookupTable[i].normalizedPos - lookupTable[i-1].normalizedPos; double ratio = (px - lookupTable[i-1].normalizedPos) / dx; return lookupTable[i-1].k + ratio * (lookupTable[i].k - lookupTable[i-1].k); } } return kMax; } double baseCohesion = 0.0; double kMin = 0.0; double kMax = 1.0; double domainHeight = 1.0; std::vector<CoeffNode> lookupTable; }; // 导出工厂函数，EDEM通过它识别并加载这个接触模型 extern "C" __declspec(dllexport) ContactModel* createContactModel() { return new VariableCohesionModel(); }这段代码的逻辑很直白：初始化时从EDEM界面读三个参数（基准内聚力、最小/最大调制系数），并预生成一张基于归一化高度的查询表；接触回调里读取颗粒坐标、归一化、插值得到k，再叠加到法向力上。整个过程没有任何复杂计算，性能上不会成为瓶颈。
4.2 代码里容易写错的两处
第一处是法向力的方向。内聚力本质是吸引力，它应该沿着法线方向往回拉。但不同EDEM API版本对法线正方向的定义可能不一样，有些版本法线是从粒子1指向粒子2，有些是反过来。我在第一版就踩过这个坑，引力方向写反了，结果颗粒之间不但不粘，反而互相排斥，整个料仓塌成蘑菇云。建议你在接入自己版本API时，先用一个只有两个颗粒的极简单点测试跑一下，打印法向力方向，确认正负号再跑大仿真。
第二处是查询表的边界处理。如果颗粒跑到料仓范围之外，normZ会小于0或者大于1。代码里我用了一个简单截断，把越界值强行拉回0/1。虽然看起来粗糙，但在大多数实际工况里，越界的颗粒只可能是初始生成瞬间的偶然溢出，截断处理完全够用。如果你的是开放域仿真，建议在初始化里就把domainHeight取得足够大，避免截断成为新的误差源。
5. 编译、加载与参数配置：让EDEM认你的插件
5.1 Visual Studio工程配置
EDEM的自定义插件本质是一个Windows DLL，开发环境我建议直接用Visual Studio 2019或2022，语言选C++。工程创建时选“动态链接库”，然后做三件事：第一，在工程属性里把EDEM安装目录下的API头文件路径加进附加包含目录；第二，把API对应的库文件路径加进附加库目录，并在附加依赖项里填上对应的lib文件名；第三，最重要的是把平台改成x64，EDEM本体是64位程序，你编译出一个32位DLL它根本不会加载。
这三步做完基本就能编译。但有一个细节很隐蔽：EDEM插件的C++运行时库要和EDEM本体保持一致。我在调试时遇到过一次插件编译纯属没问题，一加载就报错，最后发现是工程默认使用了/MT静态运行时，而EDEM本体用的是/MD动态运行时，两者混用导致底层内存分配器不一致。解决办法很直接：在C/C++代码生成选项里把运行库改成“多线程DLL(/MD)”，重新编译就正常了。
5.2 插件加载和模型挂接
DLL编译好后，进入EDEM界面，在工具栏里找到插件管理器，把你编译出的DLL路径添加进去。添加成功之后，新建一个仿真，在接触模型下拉框里就能看到你注册的模型名称。我这里再强调一个容易忽略的点：接触模型必须同时配置到“颗粒-颗粒”和“颗粒-几何体”两个接触对上，如果你只配了颗粒-颗粒，颗粒跟挡板、壁面接触时还是会走默认模型，整个仿真的内聚力效果就变成“半吊子”。
模型选上之后，你需要在参数栏里填写初始化时读取的那几个参数：BaseCohesion、Kmin、Kmax。这里的参数名必须和代码里getDouble的字符串完全一致，大小写都不能错。我第一次联调时把Kmin打成了kMin，结果模型加载后参数栏显示正常，但初始化时拿不到值，跑出来内聚力全是0，排查了半天才发现是名字不匹配。
5.3 最小验证用例
插件挂好之后，不要直接跑大场景。我强烈建议先搭一个最小验证用例：一个方形料斗，装50个颗粒，设定好自定义接触模型后在底部开一个出口，观察颗粒是否因为内聚力不同而表现出明显的粘结差异。这个用例花费时间不超过十分钟，却能一次性检验三件事：插件能不能被加载、参数能不能被正确读取、内聚力叠加方向是否正确。
我用这个用例发现过很多隐藏问题。比如有一次k的插值表写反了，底部k最小、顶部k最大，结果料斗出口处的颗粒一点都不粘，直接哗啦啦漏光。如果没有这个最小用例，直接跑工业级大场景，这种错误光靠肉眼观察结果几乎不可能定位。
6. 调试实战：加载失败、力跳变与性能瓶颈
6.1 插件加载失败的四步定位法
EDEM自定义插件调试中最常见的报错就是“插件加载失败”或者“模型不可用”。这个报错的原因其实很集中，按下面顺序排查能省下大量时间：
| 排查顺序 | 检查项 | 解决方法 |
|---|---|---|
| 1 | 编译架构是否为x64 | 工程平台改为x64，重新编译 |
| 2 | 运行时库是否用的是/MD | 工程属性里代码生成→运行库改为多线程DLL |
| 3 | 导出函数名字是否被C++改名 | 导出函数前加上extern "C"，防止命名粉碎 |
| 4 | 依赖的EDAM库版本和本机是否一致 | 确认lib文件来自当前EDEM安装目录，不要跨版本 |
第三步是最容易被忽视的。C++默认会对导出的函数做符号修饰，EDEM按固定名字去查找导出函数时找不到，就会直接判定加载失败。加一行extern "C" __declspec(dllexport)是最省事的解决办法。另外注意，EDEM不同大版本之间API接口可能不兼容，你按EDEM 2021的API编译的DLL，拿到EDEM 2023上加载，大概率要报错。所以开发时最好用和仿真一致的EDEM版本。
6.2 力跳变和颗粒爆炸：时间步与光滑化
可变内聚力模型跑起来的典型翻车现场是：颗粒在某个位置突然被粘住，下一瞬间又弹飞，整个模拟很快就发散。这通常有两个原因。
第一个原因是时间步太长。接触力在法向方向是高度非线性的，当内聚力在空间上发生突变的区域，颗粒在几步之内经历的力变化非常剧烈。解决方法是把EDEM的固定时间步适当调小，经验法则是让相邻时间步内颗粒移动距离小于最小颗粒半径的1/20。代价是计算量上升，但稳定性比什么都重要。
第二个原因是k的插值曲线有尖角。比如你在查询表里让k从1.0瞬间跳到0.3，这个阶跃映射到接触力上就是一个非常大的跳变，颗粒当然会抖动。解决办法是给插值曲线加过渡段，不要用阶跃点，用一段斜率有限的线性过渡或者平滑样条过渡。我在实际项目里的经验是：过渡段至少占整个坐标范围的5%，跳变力就能被有效缓冲。
6.3 性能瓶颈：查询表预计算比什么都管用
最后说性能。可变内聚力模型如果比EDEM内置模型慢一个数量级，那再好的物理效果也落不了地。我自己的经验是，性能优化核心在于“把能提前算的全提前算”。你在初始化阶段把查询表、插值系数、甚至多项式拟合系数全部算好，接触回调里只做查表和乘法，这样单次回调的额外耗时能控制在微秒级，跟内置模型的差距就非常小。
有一个反面教训：我见过有人为了模拟“含水率随位置非线性变化的内聚力”，在接触回调里直接调用了一个复杂的指数对数混合函数，结果同一个仿真，用内置JKR模型跑两小时，用他的自定义模型跑了两天。后来改成查表插值，时间从两天压回三个小时，物理结果几乎一模一样。如果你的可变内聚力公式特别复杂，请一定记住：先计算成表，再插值查表。这是这个项目里最想强调的实战经验之一。
这个项目做完之后再回头看，2_Variable_Cohesion_API_EDEM_的价值不只是“让内聚力可变”这一个功能点，而是给了整个仿真流程一个更灵活的数据接入方式：坐标、时间、颗粒属性都能变成驱动接触行为的变量，后续接实验数据、接CFD-DEM耦合结果、接数字孪生场的路径都被打通了。你不需要迷信API有多神秘，它就是求解器在关键节点上留出来的把手，攥住它，把你想加的计算逻辑挂上去，EDEM就能按你的规则干活。如果你正准备做类似的可变接触模型，我的建议是：第一版一定要选最简单的驱动方式和最直接的力叠加逻辑，先把链路跑通，再做复杂化。等你自己把加载、调试、性能这几个关卡都闯过去以后，再回头看这个功能，会发现它其实没有你想的那么难。
本文还有配套的精品资源，点击获取
