2024年4月12日发(作者：)

沈阳农业大学学报，2023，54(4)：447-455

ShenyangAgriculturalUniversityJournalof

http：//

DOI:10.3969/.1000-1700.2023.04.008

任德志，张露籍，宫元娟，等.生物炭基肥离散元模型参数标定与试验[J].沈阳农业大学学报，2023，54（4）：447-455.

生物炭基肥离散元模型参数标定与试验

任德志

a,b

，张露籍

b

，宫元娟

b

，孟

摘

军

a

，赫天一

a

（沈阳农业大学a.农学院/国家生物炭研究院/农业农村部生物炭与土壤改良重点实验室，b.工程学院，沈阳110161）

要：生物炭基肥料是近年来新兴的一种农业投入品，然而生物炭基肥料还田施用的机械化程度尚属较低水平。为促进生物

炭基肥料机械化还田的理论和实践发展，寻求最优生物炭基肥离散元模型接触参数组合，指导生物炭及炭基肥料还田机械的设

计和试制，提出了一种真实试验与仿真试验结合的方法，对生物炭基肥颗粒数学模型开展研究，验证数值方法的适用性和准确

性。以生物炭颗粒堆积角的实际测量值为目标值，通过物料特性研究了创建法向力、切向力和附着力接触模型，利用离散元仿真

软件ESSSRockyDEM对生物炭基肥颗粒离散元参数进行标定，设计Placket-Burman试验获取对堆积角有显著影响的参数，通过

最陡爬坡试验得到显著因素参数区间，最后通过Box-Behnken试验优化参数组合。结果表明：生物炭基肥料颗粒粒径占比2

mm为35%，3mm颗粒占比为42%，4mm颗粒占比为23%，颗粒平均含水率为5%，测得炭颗粒堆积角实际值为23.65°；对炭基肥

颗粒接触参数有显著影响的因素为炭—炭滚动摩擦系数、炭颗粒表面能JKR以及炭—钢滚动摩擦系数，最优化参数组合为炭—炭

证了仿真参数及试验的可行性和适用性。研究结果为生物炭基肥料离散元模型参数选取提供了新的参考。

关键词：生物炭基肥；Rocky；离散元模型；参数标定

中图分类号：S3

文献标识码：A

文章编号：1000-1700（2023）04-0447-09

开放科学（资源服务）标识码（OSID):

滚动摩擦系数0.36、炭颗粒表面能0.30、炭-钢滚动摩擦系数0.20，仿真堆积角试验值与真实试验值的平均相对误差仅为0.72%，验

CalibrationandTestingofDiscreteElementModelParameters

forBiochar-basedFertilizer

RENDe-zhi

a,b

,ZHANGLu-ji

b

,GONGYuan-juan

b

,MENGJun

a

,HETian-yi

a

（eofAgronomy/NationalBiocharInstitute/KeyLaboratoryofBiocharandSoilAmelioration,

MinistryofAgricultureandRuralAffairs,eofEngineering,ShenyangAgriculturalUniversity,Shenyang110161,China;）

Abstract：Biochar-basedfertilizerisanemergingagriculturalinputinrecentyears,however,themechanizationofbiochar-based

rtopromotethetheoreticalandpracticaldevelopmentof

mechanizedbiochar-basedfertilizerfieldreturn,andtoseektheoptimalcontactparameterscombinationofbiochar-basedfertilizer

methodcombiningrealtestsandsimulationtestsisproposedtocarryoutresearchonthemathematicalmodelofbiochar-based

discreteelementmodeltoguidethedesignandtrialproductionofbiochar-andcarbon-basedfertilizerfieldreturnmachinery;a

fertilizerpartualmeasuredvaluesofbiocharpellet

andadhesionforcecontactmodels,andthediscreteelementparametersofbiochar-basedfertilizerpelletswerecalibratedusingthe

reposeanglewereusedasthetargetvalue,andthematerialcharacteristicswerestudiedtocreatethenormalforce,tangentialforce

cketBurmanexperimentwasdesignedtoobtainsignificantfactorson

reposeangle,therangesofsignificantfactorswereobtainedthroughthesteepestclimbingexperiment,andfinallytheparameter

ultsshowedthat:theproportionsofparticlesizeofbiochar-based

valueofthereposeangleofbiocharparticleswas23.65°.Thefactorsthatsignificantlyinfluencedthecontactparametersofbiochar-

fertilizerwere35%of2mm,42%of3mm,23%of4mm,andtheaveragemoisturecontentoftheparticleswas5%,andtheactual

basedfertilizerparticleswere:biochar-biocharrollingfrictioncoefficient,biocharparticlesurfaceenergyJKR,andbiochar-steel

收稿日期：2023-05-17

基金项目：辽宁省教育厅重点项目（LSNZD201707）；现代农业产业技术体系建设专项项目（CARS-01-51）；辽宁省科学技术计划项目（2021-

MS-232）；国家自然科学基金项目（31901449）

第一作者：任德志（1988-），男，博士，讲师，从事生物炭还田炭-土耦合力学机制及本构关系研究，E-mail:************.cn

通信作者：赫天一（1987-），男，博士，讲师，从事生物炭农业应用研究，E-mail:*****************.cn

Copyright©博看网. All Rights Reserved.

imizedparameterswere:biochar-biocharrollingfrictioncoefficientof0.36,biocharparticlesurface

energyof0.30,ragerelativeerrorbetweenthesimulatedreposeangletest

studyprovidesanewreferencefortheselectionofparametersofthediscreteelementmodelofbiochar-basedfertilizer.

Keywords：biochar-basedfertilizer;Rocky;discreteelementmodel;parametercalibration

-448-

沈阳农业大学学报第54卷

valueandtherealtestvaluewasonly0.72%,whichverifi

生物质炭化还田是秸秆等农林废弃物资源化利用的可行途径，在土壤改良和作物栽培领域备受关注

[1-2]

。

目前生产中生物炭主要用于加工生物炭基肥料，施用生物炭基肥料可在不增加额外作业流程的情况下，实现

生物炭的还田，并发挥其改良土壤、促进作物生产的积极作用

[3]

。当前关于生物炭还田的研究多为探究炭土混

合后对土壤及各类农作物生长的影响

[4-5]

，而对于以何种方式高效还田及如何提高生物炭还田的机械化程度鲜

有研究。离散元法（DiscreteElementMethod）是一种用于预测散装固体行为的数值计算模拟技术，离散元模型

已被广泛用于农业领域，它可以帮助缩短和改善农业机械及其部件的设计过程

[6-9]

。利用离散元法进行仿真试

验，是反映生物炭炭基肥料还田过程的可行途径，而进行离散元仿真试验前，需对生物炭基肥料固体颗粒进行

参数标定。

国内外已有相关生物质材料方面的参数标定理论与试验方面的研究，侯俊铭等

[10]

以蓖麻休止角为响应值，

设计中心组合试验得到了蓖麻接触参数的最优组合；PACHÓN-MORALES等

[11]

利用DEM校准框架，以木制纤

维素生物质的休止角、堆积密度和附着率为响应值，得到了生物质颗粒流动和与设备交互作用的离散元仿真

模型；王韦韦等

[12]

利用EDEM软件得到了玉米秸秆粉料致密成型离散元仿真模型最优化参数。然而值得注意

的是，生物炭基肥料作为一种备受关注新兴农业投入品，尚未有离散元法参数标定和实验方面的研究。

本研究在前人研究的基础上，提出一种真实试验与仿真试验结合的方法测定炭基肥颗粒的物理特性，以

物理试验为参照建立生物炭炭基肥颗粒离散元模型并标定仿真参数，旨在为实现生物炭科学高效还田、仿真

生物炭基肥料还田过程和设计生物炭基肥料专用还田机械提供新的理论依据。

1

1.1

材料与方法

试验材料与物料特性

试验材料为沈阳农业大学国家生物炭研究院监制的生物炭基肥料，该颗粒材料平均含水率为5%，百粒质

量5.26g，球形度92%。从炭颗粒试验样品中随机抓取百粒，测定记录颗粒粒径，重复8次，炭基肥颗粒堆俯视

图和粒径分布如图1，2mm粒径占比35%；3mm颗粒占比42%；4mm颗粒占比23%。

粒径

Diameterofparticle

占比

Percentage

35%

2mm3mm4mm

42%23%

a.颗粒堆俯视图

Topviewofthepileofparticles

b.粒径分布示意图

Diagramofparticlesize

图1

Figure1

生物炭基肥颗粒及粒径分布

Biochar-basedfertilizerparticles

Copyright©博看网. All Rights Reserved.

第4期

1.2堆积角的测量

任德志等：生物炭基肥离散元模型参数标定与试验-449-

堆积角是表征颗粒物料流动、摩擦等特性的宏观参数，是仿真试验中的关键响应指标。参考球形颗粒肥

料的参数标定研究

[13-16]

，选择运用圆筒提升法测量颗粒堆积角（图2）。

堆积角θ

r

d

a.颗粒堆积角

Stackingangleofparticles

h

b.物理特性万能试验机

Physicalpropertiesuniversaltestingmachine

图2

Figure2

堆积角试验

Stackingangletest

试验所用装置和材料包括:物理特性万能试验机、内径40mm套筒、游标卡尺和载物平板。试验过程为：圆

筒置于平板中心位置，测量物理特性的886粒炭基肥充满于圆筒内；试验机夹具固定圆筒后，以0.1m·s

-1

的速

度将圆筒提起，使生物炭形成颗粒堆，静止5min，剔除滚出颗粒堆外的颗粒，利用游标卡尺测量颗粒堆底部最

大直径d和颗粒堆最大高度h，生物炭炭基肥颗粒堆积角为：

θr=arctan

()

2

·h

d

（1）

试验重复5次，结果如表1，最后计算得到堆积角平均值为23.65°。

表1

Table1

堆积角试验结果

试验次数

Testnumber

3

24.12

Resultsofstackingangletest

堆积角/(°)

Stackingangle

23.85

1

23.14

2

23.54

4

23.60

5

1.3接触模型

离散元软件ESSSRockyDEM提供了专用于预测农

Table2

表2

模型名称

Modelname

接触参数模型设置

模型设置

Modelsetting

Contactparametermodelsetup

业设备中颗粒行为所需的模型，仿真模型的选取直接

影响仿真结果，模型包括重力模型、接触模型以及热模

型。炭基肥颗粒之间的作用存在粘结现象，合理的选

择接触模型可以得到与实际试验相近的结果。本研究

设置接触参数如表2。

1.3.1法向力模型DEM模型中的法向力模型需要满

足两个条件：该力是一个反作用力；法向力模型满足高

能量耗散。本研究中法向力模型选择HertzianSpring

非线性函数

[17]

，其计算公式为：

切向力模型

Tangentialforcemodel

附着力模型

Adhesiveforcemodel

法向力模型

Normalforcemodel

线性弹簧库伦极限模型

Linearspringcoulomblimitmodel

JKR附着力模型

JKRadhesionmodel

赫兹弹簧-水箱模型

Hertzianspringdashpotmodel

滚动阻力模型

Rollingresistancemodel

模式C：线性弹簧滚动极限模型

TypeC:Linearspringrollinglimit

Dashpot模型，该模型与LinearSpring-dashpot模型区别在于Hertzian模型中法向力的弹性和阻尼分量是重叠的

Copyright©博看网. All Rights Reserved.

-450-

^

沈

^

阳农业

^

F

n

=K

s

+C

ss

n

3

2

H

n

大学

^

1

4

H

n

学报第54卷

（2）

s

n

为接触法线重迭量；

s

n

为触法线重迭量时间导数。式中：F

n

为法向力值；

K

H

为刚度系数；

C

H

为阻尼系数；

根据TSUJI等

[18]

研究，HertzianSpringDashpot模型中刚度系数阻尼系数的定义与LinearSpring-dashpot模

型中的定义类似，法向阻尼系数的值可通过黏性能量耗散与非弹性碰撞的能量耗散相匹配来确定，而非弹性

碰撞的能量耗散则由恢复系数的值来确定，则刚度系数和阻尼系数为：

^

4

K

H

=

E

∗

R

∗

3

C

H

=2η

H

^

（3）

（4）

m

∗

K

H

^

R

∗

为有效或等效半径；

η

为阻尼比，

m

∗

为接触的有效质量式中:

E

∗

为降低后的杨氏模量；其值与恢复系数有关；

或等效质量。其中计算公式为：

1

ì

1

ï

m

+

m

1

ï

12

=

í

m

∗

ï

1

ï

î

m

（5）

m

1

和

m

2

分别为相互接触粒子的质量，式中：而

m

是与边界接触的粒子的质量。

较好的适配性，综上考虑选择HertzianSpringDashpot模型。

1.3.2切向力模型

HertzianSpringDashpot模型适用于球体颗粒，并且该模型适用于所有切向力模型，与3种附着力模型也有

切向力模型用于计算接触力的切向分量，包括LinearSpringCoulombLimit模型、Coulomb

Limit模型和Mindlin-Deresiewicz模型。本研究选择LinearSpringCoulombLimit模型，该模型中的切向力是弹

性摩擦力，若切向力视为纯弹性，其与时间t的关系式为：

F

τ

t

,e

=F

τ

t-Δt

-K

τ

Δs

τ

（6）

F

τ

t

,e

为模型切向力；

F

τ

t-Δt

为前一瞬间切向力的值；

K

τ

为切向刚度。

Δs

τ

为粒子的切向相对位移；式中：

式（6）中的切向刚度与粒子材料相关。在此模型中，切向力不能超过库伦极限，因此切向力的完整表达式为：

F

τ

t

,e

ttt

F

τ

=min

|

F

τ,e

|

,μF

n

（7）

t

|

F

τ,e

|

F

n

t

为随时间变化的接触法向力；

μ

为摩擦系数，式中：若在接触处无滑动选择静摩擦系数

μ

s

，存在滑动选择动摩

()

擦系数

μ

d

。

当切向力

|

F

τ

t

,e

|

的大小第1次超过

μ

s

F

n

t

的极限，则认为在接触上发生了滑动，一旦该力低于

μ

s

F

n

t

的值，则认

为在接触处无滑动。

siewicz模型则会禁用滚动阻力模型。综上考虑选择适用于所有模型的LinearSpringCoulombLimit模型。

1.3.3附着力模型

在模拟计算炭基肥颗粒中引用CoulombLimit模型不适用于颗粒少精度高的仿真试验；选择Mindlin-Dere⁃

上述模型中描述的正常接触模型适合于模拟非粘性和干燥的颗粒材料，实际上，绝大多

数颗粒状材料都会含有不同的水分，所以会产生粒子与粒子、粒子与其他表明黏附的现象。为了捕捉这种行

为，必须补充附着力模型，以准确预测其流动特性。需要注意的是，颗粒之间的粘附力是应力的一个函数，

Rocky中的附着力模型通过将黏附力和接触力缩放来捕捉这一现象。

F

n,adh

=8π⋅ΓE

∗

a

3

ian理论中所所预测的略大，这是表面能加入到系统总能量中的结果。选择JKR模型后，法向力的附着力公式为：

本研究中选择的附着力模型为JOHNSON等

[19]

提出的JKR模型，在JKR模型中，粒子间的接触面积比Hertz⁃

（8）

E

∗

为有效的杨氏模量；

a

为粒子与粒子或与表式中：F

n,adh

为法向附着力；该系数的值与材料有关；

Γ

为模型参数，

明之间的接触半径。

的函数表达，即：

JKR模型会影响法向和切向力模型，HertzianSpringDashpot模型中的弹性部分需要修正，并以接触半径

a

Copyright©博看网. All Rights Reserved.

第4期

任德志等：生物炭基肥离散元模型参数标定与试验

s

n

=

a

2π⋅Γa

-

R

∗

E

∗

2

()

1

2

-451-

（9）

式中:

R

∗

为有效粒子半径。

要校准并且可以用于模拟球形颗粒。在Rocky软件中JKR模型需要与HertzianSpring-dashpot模型一起使用。

1.3.4滚动阻力模型

JKR模型拥有夯实的理论基础，被广泛应用于附着力模型中，由于表面能可以通过实验测量，该模型不需

滚动阻力是在建模中引入与粒子滚动运动相反的力矩时使用的通用名称，通常将该力

矩作为一种实用的方法来表示滚动球体上的非球形效应或其他类型粒子上的表面不规则性效应。本研究选

择的滚动阻力模型为TypeC:LinearSpringRollingLimit

[20]

。

阻力行为，其中滚动刚度为：

LinearSpringRollingLimit是一个弹塑性模型，该模型将黏性阻尼项替代为滚动刚度值以更好的模拟滚动

K

r

=R

2

r

K

τ

（10）

K

r

为滚动刚度；

K

τ

为切向刚度；

R

r

为滚动半径，式中：定义为在给定时间内连接粒子中心点和接触点的向量，若

滚动阻力为纯弹性力，则滚动阻力为：

M

r

t

,e

为模型滚动阻力；

M

r

t-Δt

为前一瞬间滚动阻力值；

K

r

为滚动刚度；

ω

rel

为相对角速度矢量，式中：定义为两个接

Δt

为仿真中模拟时间步长。触的粒子的角速度之差；

M

r

t

,e

=M

r

t-Δt

-K

r

ω

rel

Δt

（11）

作用，所以滚动阻力模型选择可配合附着力模型的TypeC:LinearSpringRollingLimit模型。

1.4离散元仿真模型

Rocky软件中提供两种滚动阻力模型，另一模型TypeA只适用于光滑滚动，由于生物炭颗粒之间存在黏连

Rocky软件材料本构提供了用于描述农业颗粒行为所需的数学模型，通过Rocky建立生物炭粒子。基于已

完成的颗粒物理特性试验得到生物炭颗粒球形度为92%，根据离散元法理论

[21]

，当颗粒非球形度小于10%时，

设置球体颗粒对仿真分析结果影响极小，所以建立Biochar粒子，形状类型设为球体颗粒，粒径分布参考物理特

性试验设置为:4mm直径的粒子数目占总颗粒数25%、3mm直径的粒子数目占总颗粒数40%、2mm直径的粒

子数目占总颗粒数35%。在SolidWorks中建立圆筒和钢板的三维模型，导出为STL格式并导入Rocky软件中

（图3）。

为更好地对照真实试验，参考设置圆筒上升速度为0.1m·s

-1

，粒子生成数目886个，生成时间2s，总仿真时

长8s。为记录设置不同的接触参数组合下所产生的炭颗粒堆积角，在每次仿真结束后创建在t=8s处粒子的

交叉散点图，它可以观测到在某一瞬间粒子堆积在X-Y平面上的投影（图4），根据交叉散点图可确定本次接触

参数组合所得到的堆积角，计算方法同式（1）。

-0.3

Z/m

0

Y

平

面

投

影

/

m

m

Y

-

p

l

a

n

e

p

r

o

j

e

c

t

i

o

n

生物炭颗粒Biochar

140

120

100

80

60

40

20

0

X/m

0.3

0

-120

-80

-40

0

8080120

0.3

0.25

Y/m

0.13

0

-0.3

X平面投影/mm

X-planeprojection

图3

Figure3

仿真模型

Simulationmodel

图4

Figure4

交叉散点图

Cross-scatterplot

Copyright©博看网. All Rights Reserved.

-452-

沈阳农业大学学报

表3接触参数模型设置

Contactparametermodelsetup

因素

Factor

第54卷

2

2.1

试验设计与结果分析

Placket-Burman试验

由于物料的接触参数中并非所有参数都对堆积

Table3

角有显著影响，所以，本研究设计Placket-Burman试

验，筛选出关键仿真接触参数，P-B试验因素水平如

表3。

P-B试验方案及结果如表4，其中H、J、K、L为虚

拟参数，即空白列。对试验结果进行方差分析如表5，

结果表明炭-炭滚动摩擦系数（C）、炭-钢滚动摩擦系

数（F）和炭颗粒表面能（G）对生物炭颗粒的堆积角影

响显著。

2.2最陡爬坡试验

对选择的显著影响参数进行最陡爬坡试验以确

定显著性参数最优值区间，试验方案及结果如表6，结

表4

Table4

试验序号

Testnumber

1

2

3

4

5

6

7

8

9

10

11

12

13

C（炭-炭滚动摩擦系数）

C（Biochar-biocharrollingfrictioncoefficient）

E（炭-钢静摩擦系数）

E（Biochar-steelstaticfrictioncoefficient）

G（表面能）

G（SurfaceenergyJKR）

D（炭-钢恢复系数）

D（Biochar-steelrecoveryfactor）

B（炭-炭静摩擦系数）

B（Biochar-biocharstaticfrictioncoefficient）

A（炭-炭恢复系数）

A（Biochar-biocharrecoveryfactor）

0.20

0.30

0.10

0.40

0.40

0.10

0.00

-1

水平

Level

0.80

1.00

0.60

0.80

1.00

0.35

0.05

1

F（炭-钢滚动摩擦系数）

F（Biochar-steelrollingfrictioncoefficient）

Plackett-Burman试验方案及结果

试验因素

Factors

F

1

-1

-1

1

1

1

-1

1

-1

1

-1

-1

1

Plackett-Burmantestprotocolandresults

堆积角/(°)

Stackingangle

23.47

36.57

24.60

41.20

10.21

20.23

4.65

12.76

19.18

38.45

17.98

20.42

23.47

A

-1

1

-1

1

1

-1

-1

1

-1

-1

1

1

-1

B

-1

1

1

-1

1

-1

-1

1

1

1

-1

-1

-1

C

1

1

1

1

-1

-1

-1

-1

-1

1

-1

1

1

D

-1

-1

1

1

-1

1

-1

1

1

-1

-1

1

-1

E

1

-1

-1

-1

-1

-1

-1

1

1

1

1

1

1

G

-1

1

-1

1

-1

1

-1

-1

1

1

1

-1

-1

H

1

-1

1

1

1

-1

-1

-1

1

-1

1

-1

1

J

1

1

-1

-1

1

1

-1

-1

1

-1

-1

1

1

K

1

1

1

-1

-1

1

-1

1

-1

-1

1

-1

1

L

-1

-1

1

-1

1

1

-1

-1

-1

1

1

1

-1

表5

Table5

方差来源

Sourceofvariation

A

B

C

D

E

F

G

效应

Effect

1.43

2.30

16.62

1.18

0.87

3.82

12.92

Plackett-Burman试验方差分析

均方和

Meansumofsquare

6.11

15.92

828.34

4.15

2.25

43.78

500.52

影响率/%

Influencerate

0.43

1.13

58.63

0.29

0.16

3.10

35.43

VarianceanalysisofPlackett-Burmantest

表6

Table6

试验序号

Testnumber

1

2

3

4

5

6

试验因素

Testfactor

G

0.00

0.10

0.20

0.30

0.40

0.50

最陡爬坡试验方案及结果

Steepestclimbtestprotocolandresults

堆积角/(°)

Stackingangle

17.65

19.33

22.14

27.31

30.44

33.21

相对误差/%

Relativeerror

25.35

18.26

6.38

15.49

28.71

40.42

C

0.10

0.20

0.30

0.40

0.50

0.60

F

0.10

0.15

0.20

0.25

0.30

0.35

Copyright©博看网. All Rights Reserved.

第4期

任德志等：生物炭基肥离散元模型参数标定与试验-453-

果表明:随着炭-炭滚动摩擦系数（C）、炭-钢滚动摩擦系数（F）和炭颗粒表面能（G）的值逐渐增大，仿真堆积角

逐渐增大，相对误差先减小再增大，其中第3组试验的相对误差最小，则选择第3组显著性参数值为中心水平，

建立回归模型求解参数组合最优值。

2.3Box-Behnken响应面试验

根据最陡爬坡试验结果，选择第3组试验因素值为中心水平值，第2组、第4组试验因素值为-1和1水平

值，设计Box-Behnken响应面试验，试验方案和结果如表7。

表7

Table7

试验序号

Testnumber

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

Box-Behnken试验方案及结果

试验因素

Testfactor

G

1.00

0.00

0.00

0.00

1.00

-1.00

0.00

-1.00

0.00

1.00

0.00

0.00

-1.00

1.00

-1.00

0.00

0.00

Box-Behnkentestprotocolandresults

堆积角/(°)

Stackingangle

16.40

20.12

19.25

20.45

23.14

12.30

19.66

20.64

19.13

21.50

12.56

10.34

22.93

22.36

10.47

19.84

21.35

C

-1.00

0.00

0.00

0.00

0.00

0.00

0.00

1.00

1.00

0.00

-1.00

-1.00

0.00

1.00

-1.00

0.00

1.00

F

0.00

0.00

0.00

0.00

-1.00

-1.00

0.00

0.00

-1.00

1.00

1.00

-1.00

1.00

0.00

0.00

0.00

1.00

通过DesignExpert试验分析软件对试验结果进行显著性检验，得到关于堆积角的方差分析如表8。通过检

验结果可发现，模型R

1

的显著值p<0.001，并且失拟项p>0.05，说明本次多因素试验的回归方程检验达到高度显

著，并且拟合度较好，因此可以对本次试验数据结果进行进一步的分析和优化。去除不显著项后得到的多项

式回归方程为：

R

1

=19.86+4.21A+2.13B+1.68C-1.05AB-3.07BC-3.26A

2

+0.86B

2

-0.76C

2

（12）

B

为炭颗粒表面能

C

为炭-钢滚动摩擦系数式中:

A

为炭-炭滚动摩擦系数（C）编码值；（G）编码值；（F）编码值。

将DesignExpert中得到的多项式回归方程导入绘图软件Origin的矩阵方程中，通过3D映射曲面图绘制出

三因素对堆积角的响应面图如图5。通过方差分析和响应面图可知，炭-炭滚动摩擦系数（C）和炭颗粒表面能

（G）单独对堆积角影响显著；炭颗粒表面能（G）和炭-钢滚动摩擦系数（F）的交互影响对堆积角影响显著，并且

当炭颗粒表面能（G）增大时，堆积角越趋近于真实试验值23.65°。综上可知炭颗粒表面能（G）对堆积角的影响

最为关键。

通过Design-Expert软件中试验数据参数优化模块，设置边界条件为堆积角趋近于目标值23.65°，得到最

优参数组合为:炭-炭滚动摩擦系数0.36，炭颗粒表面能0.30，炭-钢滚动摩擦系数0.20。

2.4验证试验

为验证所得最优化接触参数值的准确性，以炭-炭滚动摩擦系数0.36、炭颗粒表面能0.30、炭-钢滚动摩擦

系数0.20、其他接触参数取中间水平值（炭-炭恢复系数0.5、炭-炭静摩擦系数0.65、炭-钢恢复系数0.6、炭-钢

静摩擦系数0.7，其余设置不变进行离散元仿真试验，测得堆积角为23.48°，与实际试验测定堆积角23.65°的相

对误差为0.72%。仿真试验图与实际试验图对比如图6，两者无明显差异。

Copyright©博看网. All Rights Reserved.

-454-

沈

Table8

阳

表8

农业大学学报第54卷

Box-Behnken试验方差分析

堆积角/(°)

Stackingangle

均方

Meansquare

142.05

36.38

22.55

4.43

0.00

37.64

44.73

3.14

2.43

0.58

1.08

0.21

32.57

VarianceanalysisofBox-Behnkentest

方差来源

Sourceofvariation

模型

Model

C

G

F

CG

CF

GF

C2

G2

F2

残差

Residual

失拟

Lackoffit

误差

Error

总和

Summation

24

22

堆

积

角

/

(

°

)

S

t

a

c

k

i

n

g

a

n

g

l

e

20

平方和

Quadraticsum

293.17

142.05

36.38

22.55

4.43

0.00

37.64

44.73

3.14

2.43

4.07

3.24

0.83

297.23

24

22

堆

积

角

/

(

°

)

S

t

a

c

k

i

n

g

a

n

g

l

e

20

自由度

Degreeoffreedom

9

1

1

1

1

1

1

1

1

1

7

3

4

16

F

56.08

p

<0.0001

<0.0001

<0.0001

0.0004

0.0280

1.0000

<0.0001

<0.0001

0.0531

0.0801

0.0723

244.55

62.63

38.81

7.63

0.00

64.80

77.01

5.40

4.18

5.21

24

22

堆

积

角

/

(

°

)

S

t

a

c

k

i

n

g

a

n

g

l

e

20

18

14

12

16

18

14

12

16

18

14

12

16

0.100.400.30

0.16

0.16

）

c

i

e

n

t

0.15

F

0.350.350.25

0.18

0.18

（

f

炭

BB

i

o

f

数

0.20

F

）

0.200.20

颗

粒

0.20

o

e

n

t

c

h

a

炭

-

炭

0.30

）

i

o

c

h

a

炭

-

炭

0.30

擦

系

t

i

o

n

c

0.22

0.22

（

G

r

-

B

r

-

B

（

e

f

f

c

i

e

0.25

表

0.250.250.15

滚滚

能

数

摩

c

面

o

i

o

c

i

o

c

面

r

i

系

c

h

a

r

动

摩

擦

h

a

r

动

摩

擦

J

K

R

能

（

G

0.20

0.30

0.20

0.24

钢

滚

动

l

l

i

n

g

f

0.10

0.24

摩

擦

c

t

i

o

n

粒

表

R

r

r

o

l

l

）

i

颗

-

o

o

系系

r

K

r

动

l

l

i

n

J

炭

i

n

g

炭

t

e

e

l

数

g

f

（

g

f

r

数

钢

滚

r

o

l

l

i

n

f

r

i

c

（

C

）

S

-

i

c

t

i

C

）

-

r

t

i

o

n

l

a

炭

t

e

e

o

n

c

o

e

B

i

o

c

h

c

o

e

S

f

f

c

i

f

f

c

i

a

r

-

h

c

e

n

t

e

n

t

B

i

o

b.炭-炭滚动摩擦系数与炭-钢滚动摩擦系数的c.表面能与钢滚动摩擦系数的交互作用

a.炭-炭滚动摩擦系数与炭颗粒表面能的交互作用

交互作用

InteractionofsurfaceenergyJKRandbiochar-steel

Interactionofbiochar-biocharrollingfriction

Interactionofbiochar-biocharrollingfrictioncoeffi⁃rollingfrictioncoefficient

coefficientandsurfaceenergyJKR

cientandbiochar-steelrollingfrictioncoefficient

0.40

图5

Figure5

多因素交互作用的响应面图

Multi-factorresponsesurfacediagram

a.实际堆积角图

Actualofstackinganglediagram

图6

Figure6

试验对比图

b.仿真堆积角图

Simulationofstackinganglediagram

Testcomparisonfigure

Copyright©博看网. All Rights Reserved.

第4期

任德志等：生物炭基肥离散元模型参数标定与试验-455-

3讨论与结论

国内关于生物炭的研究基本来源于由陈温福院士带领成立的生物炭工程技术研究中心，其研制的Ⅰ级生

物炭适用于东北地区土壤，但该团队未对生物炭如何还田以及如何更好的还田做出详细的研究。关于离散元

参数标定的研究已经趋于成熟，与其他研究相比，生物炭还田的仿真过程类似于肥料颗粒入土，但两者物理特

性的差异使得需要对生物炭颗粒进行参数标定，以为后续的生物炭还田技术研究提供理论参考，其结果的准

确性直接影响后续对于生物炭还田工作过程仿真的真实性。

本研究结果表明，基于实际物理实验堆积角测定值23.65°，标定生物炭颗粒接触参数。设定离散元接触

模型，通过Placket-Burman试验筛选得到对生物碳颗粒堆积角有显著影响的参数:炭-炭滚动摩擦系数、炭颗粒

表面能和炭-钢滚动摩擦系数。根据最陡爬坡试验结果确定参数的最优质区间，设计Box-Behnken响应面试验

建立堆积角与3个显著参数的回归模型。显著性检验和响应面分析表明:炭-炭滚动摩擦系数和炭颗粒表面能

对堆积角均有显著影响，并且炭颗粒表面能和炭-钢滚动摩擦系数之间存在明显交互作用。根据最优化参数

组合:炭-炭滚动摩擦系数0.36、炭颗粒表面能0.30、炭-钢滚动摩擦系数0.20，重新设置接触参数并进行验证试

验。试验表明，在此接触参数组合下仿真得到的堆积角23.48°，与实际测定值相对误差为0.72%，所标定的生

物炭颗粒离散元模型参数真实可靠。

参考文献：

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

陈温福,张伟明,孟

魏永霞,肖敬萍,王

徐国鑫,王子芳,高

程效义,孟

曾智伟,马

军,等.生物炭应用技术研究[J].中国工程科学,2011,13(2):83-89.

鹤,等.施加生物炭对黑土区坡耕地改土培肥效应的持续影响[J].农业机械学报,2021,52(3):305-314.

明,等.秸秆与生物炭还田对土壤团聚体及固碳特征的影响[J].环境科学,2018,39(1):355-362.

夏桂敏,王宇佳,王淑君,等.灌溉方式与生物炭对花生根系、磷素利用及产量的影响[J].农业机械学报,2022,53(2):316-326.

VILLETTES,COINTAULTF,PIRONE,fugalspreading:Ananalyticalmodelforthemotionoffertiliserparticles

onaspinningdisc[J].BiosystemsEngineering,2005,92(2):157-164.

旭,曹秀龙,等.离散元法在农业工程研究中的应用现状和展望[J].农业机械学报,2021,52(4):1-20.

军,黄玉威,等.生物炭对玉米根系生长和氮素吸收及产量的影响[J].沈阳农业大学学报,2016,47(2):218-223.

LEBLICQT,SMEETSB,RAMONH,eteelementapproachformodellingthecompressionofcropstems[J].Com⁃

方会敏,姬长英,FARMANALICHANDIO,等.基于离散元法的旋耕过程土壤运动行为分析[J].农业机械学报,2016,47(3):22-28.

calibrationusingbulktests[J].AdvancedPowderTechnology,2019,30(4):732-750.

putersandElectronicsinAgriculture,2016,123:80-88.

[10]侯俊铭,李金澎,姚恩超,等.典型蓖麻蒴果离散元参数标定与分析[J].沈阳农业大学学报,2019,50(5):565-575.

[11]PACHÓN-MORALESJ,DOH,COLINJ,ellingforflowofcohesivelignocellulosicbiomasspowders:Model

[12]王韦韦,蔡丹艳,谢进杰,等.玉米秸秆粉料致密成型离散元模型参数标定[J].农业机械学报,2021,52(3):127-134.

[13]温翔宇,袁洪方,王

[15]贾富国,韩燕龙,刘

[16]彭飞,王红英,方

[14]李永祥,李飞翔,徐雪萌,等.基于颗粒缩放的小麦粉离散元参数标定[J].农业工程学报,2019,35(16):320-327.

扬,等.稻谷颗粒物料堆积角模拟预测方法[J].农业工程学报,2014,30(11):254-260.

刚,等.颗粒肥料离散元仿真摩擦因数标定方法研究[J].农业机械学报,2020,51(2):115-122,142.

芳,等.基于注入截面法的颗粒饲料离散元模型参数标定[J].农业机械学报,2018,49(4):140-147.

[17]eBerührungfesterelastischerKörper(Onthecontactofelasticsolids)[J].JournalFüRDieReineUnd

[18]TSUJIY,TANAKAT,giannumericalsimulationofplugflowofcohesionlessparticlesinahorizontalpipe

[19]JOHNSONKL,KENDALK,dingsoftheRoyalSoci⁃

[20]WENSRICHCM,gfrictionasatechniqueformodellingparticleshapeinDEM[J].PowderTechnol⁃

[21]赵仕威.颗粒材料物理力学特性的离散元研究[D].广州:华南理工大学,2018.

ogy,2012,217(none):409-417.

etyA:Mathematical[J].PhysicalandEngineeringSciences,1971,324(1558):301-313.

[J].PowderTechnology,1992,71(3):239-250.

AngewandteMathematik,1896,92:156-171.

[责任编辑亓国]

Copyright©博看网. All Rights Reserved.