竹材适应性强、生长期短、分布广，相比木材等其他农林生物质材料，竹材纤维长径比高，硬度和韧性更好，在家具、包装、建筑、人造板等行业占据越来越重要的市场地位^[1−3]。生物质材料的处理、运输、存储和应用依赖于适当的加工设备，而准确的接触参数是实现生物质材料机械化加工的关键，也是相关设备研发和加工参数优化的需要^[4−6]。但生物质材料的物理、力学和化学性质受与颗粒特性相关的形状、弹性、黏性、纤维性等影响，传统的研究方法难以准确获得加工设备与生物质材料粉末之间的接触参数^[4−5]。离散元法(discrete element method，DEM)是设计、建模和研究颗粒系统常用的方法，能有效模拟生物质材料采集、筛选、搅拌、粉碎和运输等过程中的颗粒行为^[7−9]。接触参数标定是颗粒系统离散元仿真计算研究的前提，国内外学者对接触参数的离散元研究主要涉及土壤、肥料、秸秆等材料。ADAJAR等^[10]根据环剪实验和直接剪切测量结果对不同作物秸秆摩擦系数、颗粒的法向刚度和剪切刚度进行了标定。ZHAO等^[11]以物理实验测定的椰糠实际堆积角为响应值，采用离散元方法获得了椰糠接触参数的最佳组合。XIA等^[12]提出了一种层析成像的DEM方法，并建立了碾磨松木颗粒的近似模型，证明了DEM模型适用颗粒状林业生物质材料的可行性。 PACHON-MORALES等^[13]分别用多球体法和Hertz-Mindlin with JKR接触模型模拟了杨木粉末颗粒的细长形状和内聚行为。但尚未发现竹材接触参数离散元研究的相关报道。

本研究对竹粉颗粒的接触参数进行标定，以毛竹Phyllostachys edulis为研究对象，基于离散元模拟仿真和实验设计(design of experiment，DOE)方法，通过Plackett-Burman (P-BD)实验、爬坡实验和响应面实验建立了毛竹粉末接触参数与堆积角的回归模型，并预测最优接触参数组合，最后验证了模拟接触参数的可靠性和准确性，以期为竹粉加工过程的模拟实验研究和相关设备的研发优化提供有效的接触参数。

1.   材料与方法

1.1   材料与装置

1.1.1   材料

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

1.1.2   堆积角测定装置

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

图 1  堆积角测定装置

Figure 1  Stacking angle measuring device

下载: 全尺寸图片 幻灯片

1.2   仿真模型

采用离散元软件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]}。

图 2  竹粉颗粒模型(A)和堆积角模型(B)

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

下载: 全尺寸图片 幻灯片

1.3   堆积角图像处理

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

图 3  堆积角图像处理

Figure 3  Image processing of stacking angle

下载: 全尺寸图片 幻灯片

1.4   参数标定

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

表 1  竹粉堆积角测定仿真参数

Table 1  Simulation parameters for measuring the stacking angle of bamboo powder

参数符号	参数名称	取值范围	参数水平
			−1	1
μ	竹粉泊松比	0.20~0.50	0.20	0.40
E	竹粉剪切模量/GPa	1~15	2	4
X1	竹粉-竹粉恢复系数	0.15~0.30	0.25	0.50
X2	竹粉-竹粉静摩擦系数	0.20~0.80	0.30	0.60
X3	竹粉-竹粉滚动摩擦系数	0~0.40	0.15	0.30
X4	竹粉-不锈钢板恢复系数	0.10~0.80	0.25	0.50
X5	竹粉-不锈钢板静摩擦系数	0.30~0.70	0.30	0.60
X6	竹粉-不锈钢板滚动摩擦系数	0~0.40	0.15	0.30
J	表面能 /(J·m−2)	0.01~0.10	0.04	0.08

下载: 导出CSV 

| 显示表格

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

表 2  竹粉堆积角P-BD实验设计及结果

Table 2  Plackett-Burman test design and results of bamboo powder stacking angle

实验号	μ	E	X1	X2	X3	X4	X5	X6	J	模拟堆积角/(°)
1	−1	−1	1	1	1	−1	1	1	−1	51.54
2	0	0	0	0	0	0	0	0	0	47.17
3	−1	−1	−1	1	1	1	−1	1	1	54.48
4	−1	−1	−1	−1	−1	−1	−1	−1	−1	42.63
5	1	−1	1	−1	−1	−1	1	1	1	41.68
6	−1	1	1	−1	1	−1	−1	−1	1	46.47
7	1	1	1	−1	1	1	−1	1	−1	37.22
8	−1	1	−1	−1	−1	1	1	1	−1	40.56
9	1	1	−1	1	1	−1	1	−1	−1	57.91
10	1	−1	1	1	−1	1	−1	−1	−1	23.71
11	−1	1	1	1	−1	1	1	−1	1	39.00
12	1	1	−1	1	−1	−1	−1	1	1	37.79
13	1	−1	−1	−1	1	1	1	−1	1	61.32
说明：符号所表示参数名称见表1。

下载: 导出CSV 

| 显示表格

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

表 3  竹粉堆积角爬坡实验设计及结果

Table 3  Climbing test design and results of bamboo powder stacking angle

实验号	一次爬坡						二次爬坡
	X1	X3	X5	模拟堆积角/(°)	相对误差/%		X1	X3	X5	模拟堆积角/(°)	相对误差/%
1	0.50	0.15	0.30	38.92	21.04		0.50	0.15	0.30	39.67	19.52
2	0.40	0.25	0.40	48.49	4.50		0.45	0.20	0.35	41.53	15.75
3	0.30	0.35	0.50	57.95	17.57		0.40	0.25	0.40	49.32	5.20
4	0.20	0.45	0.60	−	−		0.35	0.30	0.45	52.82	7.17
5	0.10	0.55	0.70	−	−		0.30	0.35	0.50	57.27	16.19
说明：−表示数据无效。符号所表示参数名称见表1。

下载: 导出CSV 

| 显示表格

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

表 4  竹粉堆积角响应面实验设计及结果

Table 4  Box-Behnken test design and results of bamboo powder stacking angle

序号	X1	X3	X5	模拟堆积角/(°)
1	0(0.40)	1(0.30)	−1(0.35)	48.69
2	0	0(0.25)	0(0.40)	48.89
3	1(0.45)	0	1(0.45)	52.17
4	1	0	−1	43.18
5	−1(0.35)	0	−1	48.71
6	0	0	0	48.51
7	−1	1	0	54.98
8	0	−1(0.20)	−1	41.53
9	−1	0	1	54.97
10	0	−1	1	48.85
11	1	−1	0	41.20
12	0	0	0	50.54
13	1	1	0	51.69
14	0	1	1	55.27
15	−1	−1	0	48.11
说明：括号内数值为参数值。参数X1、X3、X5所表示名称见表1。

下载: 导出CSV 

| 显示表格

2.   结果与分析

2.1   P-BD实验结果分析

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

表 5  竹粉堆积角P-BD实验参数显著性分析

Table 5  Significance parameters analysis of bamboo powder stacking angle Plackett-Burman test

参数	效应	系数	P	显著性排序
μ	−2.506	−1.253	0.230	7
E	−2.735	−1.368	0.203	6
X1	−9.179	−4.589	0.025	2
X2	−0.904	−0.452	0.601	9
X3	13.927	6.963	0.011	1
X4	−3.623	−1.811	0.132	5
X5	8.286	4.143	0.030	3
X6	−1.297	−0.649	0.470	8
J	4.527	2.264	0.091	4
说明：参数所表示名称见表1。

下载: 导出CSV 

| 显示表格

2.2   爬坡实验结果分析

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

图 4  有效堆积图像(A)和无效堆积图像(B)

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

下载: 全尺寸图片 幻灯片

2.3   响应面实验结果分析及回归模型

根据表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  Box-Behnken实验回归模型方差分析

Table 6  ANOVA of modified model of Box-Behnken

方差来源	均方	自由度	平方和	F	P
模型	265.87	9	29.54	86.60	＜0.000 1
X1	38.22	1	38.22	112.04	＜0.000 1
X3	108.81	1	108.81	318.97	＜0.000 1
X5	103.35	1	103.35	302.96	＜0.000 1
X1X3	5.44	1	5.44	15.93	0.010 4
X1X5	1.86	1	1.86	5.45	0.066 9
X3X5	0.03	1	0.03	0.09	0.776 9
X1X1	0.49	1	0.49	1.44	0.283 2
X3X3	7.63	1	7.63	22.36	0.005 2
X5X5	0.64	1	0.64	1.88	0.228 4
残差	1.71	5	0.34
纯误差	0.015	2	0.007
合计	267.58	14
R2=0.9936，Radj2=0.9822，CV=1.18，AP=30.541
说明：参数X1、X3、X5所表示名称见表1。

下载: 导出CSV 

| 显示表格

从表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实验得到的结论一致。

图 5  交互项X₁X₃响应面

Figure 5  Response surface of interaction term X₁X₃

下载: 全尺寸图片 幻灯片

2.4   参数优化和模型验证

利用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％，进一步证明了该模型在竹粉颗粒参数标定中的有效性。

3.   讨论

农林生物质材料的含水率影响其微观结构和物理性质，进而影响材料的生产、运输、存储和应用^[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是更高效的。

4.   结论

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

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