Calibration of contact parameters for bamboo powder particles based on DEM

毛竹采自安徽宣城，竹龄4~5年生，胸径为10.2 cm，密度为680 kg·m⁻³，静曲强度为156.97 MPa，抗压强度为66.24 MPa，取竹中部位，包含竹青和竹黄。将样品研磨成粉，干燥后，采用不同规格标准筛分得到100目竹粉样品。

采用漏斗法测定竹粉颗粒堆积角，参照相关标准和文献[5, 14−16]，自制堆积角测定装置如图1所示，短颈漏斗颈口内径10 mm，底盘培养皿直径100 mm，漏斗颈口距底盘培养皿90 mm。测量时将竹粉缓慢倒入并轻微搅动，避免颗粒在漏斗壁上积聚；待料堆高度不再增加时停止添加竹粉。在相同位置使用相机拍摄料堆的正视图像，重复10次并记录。

Figure 1.  Stacking angle measuring device

采用离散元软件EDEM模拟毛竹粉末料堆形成过程，建立毛竹粉末颗粒模型和堆积角模型，并选择适当的接触模型：①颗粒模型。竹粉实际为纤维状，且不同竹粉颗粒长宽比差异较大，测定物理堆积角时采用统一粒径的竹粉是为了降低颗粒模型的复杂程度；在对比实际粒径和放大粒径形成的堆积角后，将竹粉颗粒的粒径(宽)放大至1 mm ^[17]。建立的竹粉颗粒模型如图2A所示，此颗粒模型的长宽比为2∶1，粒径分布设置为fixed。②堆积角模型。依照物理堆积角测定装置建立几何模型并导入EDEM软件，在漏斗上方建立总数为25 000的颗粒工厂，生成速率1 250个·s⁻¹，静置4 s，总仿真时间设为24 s。待竹粉颗粒自由落下形成稳定的颗粒堆时，截取+X、−X、+Y、−Y共4个方向的料堆图像。建立的竹粉颗粒堆积角模型如图2B所示。③接触模型。毛竹粉末与加工设备之间的颗粒运动特征是密集和缓慢的，加上生物质材料的纤维性和联锁效应，干燥后的样品仍具有轻微的黏附性^[18−19]。选用EDEM软件提供的JKR Cohesion接触模型，该模型考虑了范德华力在接触区内的影响，也允许构建强黏合系统^{[16, 18]}。

Figure 2.  Bamboo powder particle model (A) and powder stacking angle model (B)

利用图像处理软件对物理实验和仿真实验获得的堆积角图像进行处理，调用灰度函数rgb2gray、二值函数im2bw、形态学运算函数strel (square和disk)得到内部填充、边缘柔和的堆积角图像；在GUI界面利用canny算子提取料堆边缘，并利用最小二乘法进行一阶线性拟合得到斜率(精确到0.001)；最后将斜率(正切值)转换为角度，堆积角处理过程如图3所示。实测物理堆积角为49.29°。

Figure 3.  Image processing of stacking angle

进行离散元模拟研究所需的仿真参数主要有颗粒密度、泊松比、恢复系数、剪切模量、静摩擦系数、动摩擦系数和表面能等^{[7, 20]}。壁面材料采用不锈钢板，其密度为7 850 kg·m⁻³、泊松比为0.3、剪切模量为7×10¹⁰ Pa，而竹粉的接触参数难以通过物理实验测量，需通过模拟实验进行标定^{[7, 14]}。结合预模拟实验、相关农林生物质材料参数标定文献[6, 11, 21−24]，本研究中竹粉和壁面材料设置仿真参数见表1。影响堆积效果的因素很多，先进行Plackett-Burman (P-BD)实验筛选显著因素，再进行爬坡实验缩小显著因素的参数范围，然后进行响应面实验获得堆积角和DEM参数之间的二次多项式回归模型，最后使用预测得到标定参数的最优组合进行验证。

① P-BD实验。采用Minitab软件9因子表进行P-BD实验设计，总计13组(包括1个中心点)，实验方案及结果见表2。

② 爬坡实验。根据P-BD实验分析结果设计爬坡实验。显著因素效应值为正，则水平逐次增加；效应值为负，则水平逐渐减少^{[11, 25−26]}。一次爬坡实验设置步长为0.1，二次爬坡实验设置步长为0.05，分别进行5组实验，爬坡实验方案及结果见表3。

③ 响应面实验。采用Design Expert软件中的Box-Behnken方法(B-BD)进行响应面实验。选择与实际堆积角误差最小的一组参数作为响应面实验设计的水平值，根据参数个数在Design Expert软件中生成响应面实验顺序表(加3个中心点)，响应面实验方案及结果见表4。

如表5所示：模型的决定系数(R²)为0.989 3，表明该模型可用于解释因素对响应值的影响。X₁、X₃、X₅对竹粉的堆积角有显著影响(P＜0.05)，影响堆积角的显著因素从大到小依次为X₃、 X₁、 X₅，其他因素影响较小。原因为干燥后的竹粉黏附性较低，因此竹粉与竹粉之间运动状态主要取决于运动方式和碰撞力度，采用多球体方法构建的颗粒模型之间多为滚动；非球形颗粒的不规则形状减小了其流动性，因此影响竹粉与不锈钢板之间运动状态的主要原因是静摩擦。

从表3可以看出：一次爬坡实验前3组堆积角相对误差分别为21.04%、4.50%、17.57%，呈先减小后增大的趋势；后2组实验的堆积图像无效。原因为恢复系数过小，静摩擦系数过大，颗粒不易运动，造成无效堆叠，堆积效果如图4。缩短步长后二次爬坡实验5组实验堆积角均有效，且根据相对误差可知最优区间在3组附近，2、4组之间，则选择2、3、4组实验的参数值作为响应面实验的各水平值。

Figure 4.  Effective (A) and invalid (B) stacking angles image

根据表4的数据，建立了竹粉堆积角和显著因素的二次回归模型：$ {\theta }_{1} $ =50.54−2.19 X₁+3.69 X₃+3.59 X₅+1.17 X₁X₃+0.68 X₁X₅−0.087 X₃X₅−0.37 X₁²−1.44 X₃²−0.42 X₅²。从表6可以看出：模型P＜0.000 1，表明该模型用于描述参数与响应值之间的关系极显著，可以根据模型预测堆积角；变异系数(C_V)为1.18%，表明该测试具有很高的可靠性；R²为0.993 6，修正决定系数($R_{{\rm{adj}}}^2 $)为0.982 2，这些数值表明回归方程的拟合度较好(一般认为2个系数需大于80%)，回归方程可以用来代替实际测试来分析结果；精度(A_P)为30.541，表明该模型具有良好的精度。

从表6可以看出：线性项X₁、X₃、X₅的P均小于0.000 1，与P-BD实验结果完全一致，在模拟实验参考范围内对堆积角影响极显著；交互项X₁X₃的P为0.010 4，对堆积角有显著影响，其余交互项X₁X₅、X₃X₅对堆积角影响很小(P＜0.05)；二次项X₃X₃的P为0.005 2，对堆积角有显著影响，X₁X₁、X₅X₅对堆积角影响很小(P＜0.05)。

当 X₅为0.40时，利用Design-Expert软件得到X₁X₃的响应曲面，如图5所示。在X₁X₃作用下竹粉堆积角整体表现为增大趋势，原因为竹粉-竹粉滚动摩擦系数的增大增加了竹粉颗粒间接触的摩擦阻力，使颗粒堆不易滑散、趋于稳定；而竹粉-竹粉恢复系数的增大增加了颗粒间的碰撞反弹程度，使颗粒堆易飞散，不易成形。但竹粉-竹粉滚动摩擦系数的响应曲面曲线较竹粉-竹粉恢复系数的响应曲面曲线更陡峭，这表明竹粉-竹粉滚动摩擦系数对堆积角的影响比竹粉-竹粉恢复系数的影响更大，这也与P-BD实验得到的结论一致。

Figure 5.  Response surface of interaction term X₁X₃

利用Design-Expert软件的预测模块，以物理堆积角49.29°作为目标值，选择一组与实际物理堆积角度数据相近的最优解，即竹粉-竹粉恢复系数X₁=0.44，竹粉-竹粉滚动摩擦系数X₃=0.22，竹粉-不锈钢板静摩擦系数X₅=0.45，其余非显著参数取值与爬坡实验相同，即泊松比μ=0.33，剪切模量E=0.033 GPa，竹粉-竹粉静摩擦系数X₂=0.4，竹粉-不锈钢板恢复系数X₄=0.3，竹粉-不锈钢板滚动摩擦系数X₆=0.1，表面能J=0.01 J·m⁻²。

采用上述标定的竹粉最优参数组合进行验证，测得堆积角为48.06°，相对误差为2.50%，误差在可接受的范围内；再选择2组测定实验目标值(47.95°和48.52°)，预测得到最优解后进行验证，分别测得堆积角为46.58°和46.65°，相对误差分别为2.86％和3.85％，进一步证明了该模型在竹粉颗粒参数标定中的有效性。

农林生物质材料的含水率影响其微观结构和物理性质，进而影响材料的生产、运输、存储和应用^[24]。含水率受样品年龄、采样部位等多因素的影响，但在含水率较低的情况下，木材散碎物料的堆积特性变化不大^[27]。本研究采用竹中部位的竹段，磨粉干燥后控制含水率低于2％，所以本研究未设计含水率对模拟堆积角的影响实验。

模拟不同比例、不同粒径、不同长宽比的竹粉纤维会增加模拟实验的复杂性，因此构建4球面、长宽比2∶1的颗粒模型，并选择fixed平均分布，与物理堆积角采用统一粒径竹粉对应；当粒径＞1.5 mm时会堵塞漏斗下端出口，而真实粒径0.15 mm需要的颗粒总量会按量级增加，相应的计算时间也会增加。当表面能超过0.12 J·m⁻²时，颗粒会黏附在漏斗下端出口造成堵塞；当摩擦系数超过0.85时，颗粒落到料堆上运动状态阻碍太大造成不规则团块堆叠。根据实际情况合理简化颗粒形状和缩放粒径，在正式实验前需要做预实验选择接触模型和参数范围。

剪切模量量级、颗粒生成速率和计算机硬件配置都会影响仿真时间。在P-BD实验得出剪切模量对堆积角的影响很小的结论后，将剪切模量降低量级到10⁷。在不同配置的中央处理器(CPU)和图像处理器(GPU)上仿真，并试验颗粒生成速率5 000、2 500、1 250个·s⁻¹时的仿真时间，发现GPU比CPU的速度稳定，例如12核CPU仿真时间约46 h，而6核GPU仿真时间约22 h；减小颗粒生成速率初期仿真较慢，但是为了降低颗粒生成位置的计算时间，会缩短后期仿真时间。在颗粒数量更大时，选择GPU是更高效的。

本研究所得主要结论如下：①通过P-BD、爬坡和B-BD响应面等方法，建立了检验因素与堆积角的回归模型，该模型能够准确描述堆积角与各参数的关系。证明了竹粉-竹粉滚动摩擦系数、竹粉-竹粉恢复系数、竹粉-不锈钢板静摩擦系数的线性项，竹粉-竹粉滚动摩擦系数的二次项，竹粉-竹粉滚动摩擦系数和竹粉-竹粉恢复系数的交互项对堆积角有显著影响。②以竹粉物理堆积角为目标，通过求解回归模型获得了最优因素组合，即竹粉-竹粉恢复系数0.44，竹粉-竹粉滚动摩擦系数0.22，竹粉-不锈钢板静摩擦系数0.45。将回归模型的最优解作为DEM参数进行验证，测得堆积角的相对误差在允许范围内。③本研究在41.20°~55.27°堆积角范围内建立的回归多项式方程可以对显著因素进行快速预测，结果具有较高的可靠性。当具体的工作环境和接触模型发生变化时，线性项和二次项的相互作用会增加标定接触参数的难度，最优解需要根据生物质物料的堆积角、料堆形状和动态行为确定。

本研究得到的回归模型证明了离散元法应用于竹材粉末研究的可行性，生物质材料特性、不同响应指标检验方法、颗粒模型的准确性等是下一步的研究方向。