中药浸膏粉离散元模拟参数标定方法研究

中药浸膏粉离散元模拟参数标定方法研究

中药浸膏粉离散元模拟参数标定方法研究


中药浸膏粉离散元模拟参数标定方法研究

摘  要目的  采用离散元法(DEM)模拟中药浸膏粉休止角,确定颗粒之间、颗粒与几何体之间的接触参数。方法  以生甘草浸膏粉、独活浸膏粉、微晶纤维素(MCC)和乙基纤维素(EC4种粉末为研究对象,在Hertz-Mindlin with JKR Cohesion接触模型和颗粒缩放的基础上,通过Plackett-Burman设计筛选出对休止角模拟测定影响显著的接触参数,进而通过最陡爬坡设计确定关键接触参数的最佳区域。根据Box-Behnken试验结果,建立接触参数和模拟休止角之间的回归模型,优选最佳接触参数值并验证。结果  筛选得到的关键接触参数为颗粒-颗粒滚动摩擦系数、颗粒-颗粒恢复系数和颗粒-不锈钢恢复系数,所建回归模型对休止角的标定范围为33.30°43.64°4种粉末休止角模拟测定值与实验值的相对误差绝对值均小于2.0%,表明所建标定方法准确可靠。结论  证明了通过宏观物性参数间接标定中药颗粒体系DEM微观力学参数的可行性,为混合、输送等中药制药过程的精确模拟奠定了基础。

离散元法(discrete element methodDEM)是解决非连续介质问题的数值模拟方法。DEM将求解空间离散为单元,并模拟所有单元在任意时刻的位移、速度、加速度、角速度等物理量,进而预测离散群体行为,被广泛应用于岩土、建筑、农业和医药等领域颗粒体系的过程模拟和装备设计等[1-4]DEM运行原理为(1)在单个时间步长内使用接触模型计算相邻颗粒或者几何体对每个颗粒施加的作用力;(2)根据牛顿第2定律计算颗粒的速度和旋转速度;(3)根据时间步长计算颗粒的新位置。这一过程在每个时间步长中应用于每个颗粒,并重复至模拟结束[5]。近年来,高性能计算的发展极大提高了DEM计算效率,可支持百万级和千万级单元的动态仿真[6],使之更接近工程实际。DEM耦合计算流体力学(CFD)、多体动力学(MBD)和有限元法(FEM)可实现复杂颗粒体系的多学科仿真[7-9]

在制药工程领域,DEM已应用于粉末混合、输送、制粒和压缩等口服固体制剂的制备过程。Park[10]使用DEM5种具有不同桨叶角度设计的连续混合器中双峰分布颗粒的混合过程进行了模拟,提出轴向和横向方差指数以评估颗粒混合均匀度,并指导不同给料条件下的混合器螺杆桨叶角度设计。Mazor[11]以微晶纤维素(MCCPH 101为研究对象,采用DEM模拟干法制粒进料区的颗粒流动,采用FEM模拟MCC在辊压区的形变,DEMFEM联用可以更好的预测螺杆转速和辊轮转速对薄片密度和均匀性的影响。

中药浸膏粉成分复杂,具有易吸湿、黏性大、流动性差等特点[12],在DEM模拟中缺乏相应仿真参数的报道。DEM仿真参数主要指颗粒与颗粒之间、以及颗粒与设备几何体之间的接触参数,如恢复系数、静摩擦系数、滚动摩擦系数等。通常来说,接触参数难以由实验测定,一般通过改变参数值并进行模拟,将模拟结果与真实实验结果进行对比,直到模拟结果与实验结果一致[13]。为了提高DEM在中药制药领域的适用性,本实验以粉末休止角(α)的测定过程为例,对接触参数进行标定,研究结果可为中药粉体混合和输送等制剂过程精细化设计奠定基础。

1  材料

生甘草浸膏粉(批号180421-20180101094)、独活浸膏粉(批号180423-445800-18)均由北京康仁堂药业有限公司提供,通过水提、浓缩和喷雾干燥工艺制备。MCC Vivapur® type 200,批号5620030813,上海昌为医药辅料技术有限公司;乙基纤维素N100 Pharm,批号44547Ashland Inc.

使用BT-2001激光粒度分布仪(丹东百特仪器有限公司)对上述4种粉末的粒径分布进行测定,累积粒径分布达10%50%90%所对应的粒径值D10D50D90见表1

中药浸膏粉离散元模拟参数标定方法研究 

2  模拟方法

2.1  接触模型

接触模型是在离散元模拟时颗粒之间相互作用力计算的依据,涉及到了颗粒的表面黏连、接触、变形等。接触模型可根据对研究对象特性的理解进行选择。中药粉末具有一定黏性,因此在本实验中选用Hertz-Mindlin with Johnson-Kendall-RobertsJKRCohesion接触模型[14],该模型为凝聚力接触模型,在接触区域中考虑了液桥力、范德华力等的影响,适于模拟强粘性系统,如细干颗粒或湿颗粒等。在Hertz-Mindlin with JKR Cohesion接触模型中,颗粒间的切向阻尼力可由式(1)求出。

Ft−2(5/6)1/2β(Stm*)1/2vtrel                 1

m*为等效质量,vtrel是切向相对速度,系数β和切向刚度St可由下面2式求出。

βlne/(ln2eπ2)1/2                       2

St8G*(R*α)1/2                           3

e为恢复系数,G*为等效剪切模量,R*为等效粒子半径

切向力与摩擦力μsFn有关,此处μs为静摩擦系数,滚动摩擦可以通过接触表面上的力矩来说明。

TiμrFnRiωi                            4

μr为滚动摩擦系数,Ri为质心到接触点的距离,ωi为接触点处物体的单位角速度矢量

通过选择接触模型、设置接触参数、以及设置几何体物理模型,可对颗粒体系的运动行为进行DEM模拟。

2.2  颗粒缩放

颗粒粒径和形状的设置影响模拟的计算速度,本实验采用的中药粉末由喷雾干燥工艺制得,粒径较小且接近球型,为了提高离散元仿真的效率并节约计算时间,将模拟颗粒形状设置为球形。对于由大量微小颗粒组成的粉体系统的模拟通常采用放大颗粒的方法。Feng[15]提出了颗粒缩放应当遵循3个相似原理:几何相似、运动相似、动力相似。几何相似指缩放模型必须严格遵守经典的几何相似性原理,运动相似是指模型与原型中对应点处的速度方向相同,大小成比例,并且具有同一比率;运动相似还要求对应时间间隔成比例并具有同一比率。动力相似是指模型和原型受到相同性质力的作用,并且这些力成比例并具有同一比率。本实验将表14种粉末的颗粒半径均放大为1 mm进行后续的模拟试验。

3  粉体休止角试验测定和模拟测定

3.1  休止角实验测定

休止角测定参照《欧洲药典》第102.9.16节,具体操作如下:使用BEP2粉体流动性测试仪(英国Copley公司),固定底面圆盘,圆盘直径(Φ)为100 mm,将粉末从喷嘴上方的漏斗缓慢倒入,喷嘴直径为10 mm,喷嘴下端距离底面圆盘的高度为75 mm。粉末在底面圆盘形成锥体,当锥体高度不再增加时,读取锥体高度(h)和底面圆盘的半径(r),平行测量3次,求平均值,休止角α°)的计算公式如下。

αarctan(h/r)                            5

3.2  休止角模拟测定

构建BEP2粉末流动性测试仪漏斗和底座部分的几何模型,如图1所示。使用离散元模拟软件EDEM2018DEM SolutionsLtd.Edinburgh)进行研究。放大颗粒半径设置为r1 mm,颗粒生成方式为Dynamic(即动态生成),生成的总量为100 g,产生的颗粒数量为2 000/s。待100 g的颗粒全部堆积在底面圆盘上并形成堆积的锥体后,采用软件后处理中记录颗粒位置的功能,记录颗粒堆积最高位置随时间变化的趋势,当两个时刻颗粒的高度变化≤0.2%,认为颗粒的堆积相对静止,记录此时锥体的高度h并根据式(5)计算休止角。

中药浸膏粉离散元模拟参数标定方法研究

4  休止角模拟过程实验方法设计

使用离散元模拟需要设置颗粒、几何体的材料参数,以及颗粒-颗粒、颗粒-几何体的接触参数。颗粒材料参数包括:颗粒粒径、颗粒密度、颗粒形状、泊松比、剪切模量。几何体物理参数包括密度、剪切模量、泊松比等,由于使用的休止角测量实验工具为不锈钢材质,因此直接使用不锈钢的模拟参数,见表2EDEM软件接触参数包括颗粒-颗粒恢复系数、颗粒-颗粒静摩擦系数、颗粒-颗粒滚动摩擦系数、颗粒-不锈钢恢复系数、颗粒-不锈钢静摩擦系数、颗粒-不锈钢滚动摩擦系数以及反映物料粘性的JKR表面能。

中药浸膏粉离散元模拟参数标定方法研究

结合文献对药物粉末与不锈钢材料离散元模拟的参数设置[16-18]以及本课题组建立的中药原辅料物性数据库[19],设置模拟参数和(或)参数变化范围如表2所示。其中粉末泊松比0.3,固体密度为1 500 kg/m3,剪切模量100 MPa。表14种粉末与几何体的接触参数与模拟材料的密度、形状、粒径等性质有关,文献资料尚无报道,且并非每个接触参数对模拟结果都有显著性影响,因此,由标定方法获得[20]。首先采用Plackett-Burman试验从7个接触参数中筛选出对休止角模拟影响较大的参数,之后进行最陡爬坡试验使得模拟休止角结果进入最佳区域,最后通过Box-Behnken试验优化确定模拟休止角的最佳接触参数。

5  结果与分析

5.1  Plackett-Burman试验

Plackett-Burman设计是一种2水平(−1+1水平)的试验设计方法,它试图用最少试验次数筛选并确定对结果影响比较显著的因素。本实验Plackett-Burman设计以粉末休止角为响应值,对影响模拟且不易确定的7个接触参数颗粒-颗粒恢复系数(particle-particle restitutioncoefficientA)、颗粒-颗粒静摩擦系数(particle-particle static friction coefficientB)、颗粒-颗粒滚动摩擦系数(particle-particle rolling friction coefficientC)、颗粒-不锈钢恢复系数(particle-stainless steelrestitution coefficientD)、颗粒-不锈钢静摩擦系数(particle- stainless steelstatic friction coefficientE)、颗粒-不锈钢滚动摩擦系数(particle-stainless steelrolling friction coefficientF)和JKR表面能(Johnson Kendall Roberts surfaceenergyG)进行筛选。Plackett-Burman试验因素水平表如表3所示,采用Design Expert 8.0.6软件(Stat-Ease, IncUnited States)生成13次试验,使用“3.2”项方法进行休止角模拟测定,试验设计方案及结果见表3。颗粒全部下落后,堆积高度随时间的变化曲线见图2。图2中,曲线起点表示在8.2 s时颗粒全部生成,曲线终点表示模拟的休止角测定总时间为20 s

中药浸膏粉离散元模拟参数标定方法研究

中药浸膏粉离散元模拟参数标定方法研究

对表3模拟的结果进行方差分析,结果如表4所示。可见,颗粒-颗粒滚动摩擦系数、颗粒-颗粒恢复系数、颗粒-不锈钢恢复系数的P值小于0.01,对粉末休止角的模拟结果影响极其显著;颗粒-颗粒静摩擦系数的P0.05,对粉末休止角的模拟结果影响显著;其余3个因素的P值>0.05,对粉末休止角的模拟结果影响较小。在后续的最陡爬坡试验以及Box-Behnken实验中选择3个对休止角的模拟影响极其显著(P0.01)的因素;根据已有的研究以及作者对于药物粉末的理解[21-23],将其余参数的数值固定为:颗粒-颗粒静摩擦系数0.6、颗粒-不锈钢静摩擦系数0.4、颗粒-不锈钢滚动摩擦系数0.15JKR表面能0.01来进行后续的试验设计。

中药浸膏粉离散元模拟参数标定方法研究

5.2  最陡爬坡试验

Plackett-Burman试验后,根据筛选出的对模拟影响显著的因素,进行最陡爬坡试验,以便能够快速进入最佳的响应区域。本实验颗粒-颗粒恢复系数(A)爬坡步长设定为0.2,颗粒-颗粒滚动摩擦系数(C)爬坡步长设定为0.15,颗粒-不锈钢恢复系数(D)爬坡步长设设定为0.2,最陡爬坡试验设计如表5所示。采用“3.1”项的方法进行休止角实验测定,生甘草浸膏粉、独活浸膏粉、MCC和乙基纤维素的休止角实测值分别为41.18°43.57°38.23°34.91°。根据EP 10.0休止角大小和粉体流动性的关系,生甘草浸膏粉和独活浸膏粉的休止角在41°45°,属于流动性尚可、可处理的范畴;MCC休止角在36°40°,表明流动性适当、无需辅助的范畴;乙基纤维素休止角在31°35°,说明其流动性良好。

采用“3.2”项的方法进行休止角模拟测定,结果如表5所示。分别计算4种粉末休止角和模拟休止角的相对误差=(实验休止角模拟休止角)/实验休止角,然后计算4种粉末分别在每个实验条件下的平均相对误差,结果13号实验的平均相对误差分别为9.81%8.39%19.92%。可见在12号水平的平均相对误差较小,且有3种粉末的实验测定休止角落在了1号水平和2号水平的范围内,因此在后续的Box-Behnken试验中以1号水平为低水平,2号水平为高水平进行休止角模拟测定实验。

中药浸膏粉离散元模拟参数标定方法研究

5.3  Box-Behnken试验及回归模型

根据“5.2”项确定的因素范围,对颗粒-颗粒恢复系数、颗粒-颗粒滚动摩擦系数、颗粒-不锈钢恢复系数3个因素进行Box-Behnken试验设计,因素水平见表6。采用“3.2”项下方法进行休止角模拟测定,结果见表6

中药浸膏粉离散元模拟参数标定方法研究

根据表6实验结果对自变量和响应变量进行全项2次多项式回归拟合,得到模型1α26.229.88 A76.42 C12.19 D3.73 AC59.22 AD11.43 CD28.10 A237.67 C252.25 D2。模型方 差分析结果表明所得拟合模型P0.000 1,失拟项PL0.482 50.05,决定系数R20.995 4,表明回归模型拟合良好。颗粒-颗粒恢复系数(A)、颗粒-颗粒滚动摩擦系数(C)、颗粒-不锈钢恢复系数(D)的P值均小于0.01,对休止角的影响极其显著;颗粒-颗粒恢复系数颗粒-不锈钢恢复系数(AD)以及颗粒-不锈钢恢复系数的二次项(D2P0.05,对休止角的影响显著。将对休止角影响不显著(P0.05)的项(ACCDA2C2)移除,得到优化的2次多项式回归模型,即模型2α32.3328.71 A56.77 C7.15 D59.22 AD48.73 D2。模型2P0.000 1,失拟项PL0.575 30.05,决定系数R20.992 5,表明优化后的模型仍具有良好的可靠性。响应值的范围为33.30°43.64°,该范围涵盖了表14种物料实验测定休止角,表明模型2可用于表14种物料的接触参数ACD的标定。

5.4  模拟参数优化和验证

以生甘草浸膏粉接触参数ACD的标定为例,以休止角试验测定值41.18°为目标,基于“5.3”项优化后的回归模型2,应用Design Expert软件数值优化求解功能,得到生甘草浸膏粉的接触参数为颗粒-颗粒恢复系数0.32,颗粒-颗粒滚动摩擦系数0.26,颗粒-不锈钢恢复系数0.34;将该组参数代入模型2,得到休止角预测值为41.14°;使用该组参数组合对休止角进行DEM模拟,所得的休止角为41.48°,与真实值41.18°的相对误差为−0.73%。采用类似优化方法,对其他3种粉末的接触参数进行优化求解和验证,结果如表7所示。可见对4种粉末模拟所得休止角与真实测定休止角的相对误差绝对值均小于2.0%,表明优化所得接触参数可用以模拟粉末休止角。

中药浸膏粉离散元模拟参数标定方法研究

恢复系数用以衡量2个物体在碰撞后的反弹程度。弹性碰撞时恢复系数为1;完全非弹性碰撞时恢复系数为0,即2个物体黏贴在一起[24]。本实验中4种粉末的恢复系数在0.290.35,属于非弹性碰撞的范畴,颗粒碰撞时损失的系统总动能较多。滚动摩擦系数与颗粒表面性质有关,例如表面越粗糙,滚动摩擦系数越大[25]。本实验中2种中药粉末的滚动摩擦系数均大于2种纤维素类辅料的滚动摩擦系数,表明2种中药物料颗粒间的滚动摩擦和阻碍较大。

6  讨论

6.1  颗粒粒径的影响

在离散元模拟中,颗粒粒径大小及分布是影响DEM模拟效果和效率的关键因素之一。在实际应 用中一般将颗粒简单的球形化并进行放大,从而更适合应用于颗粒数目众多的制药过程。Hassanpour[26]对比了在桨叶式混合器中真实实验与不同大小粒径离散元模拟的结果,将他们的粒径分别设置2.264.527.2011.40 mm,在相同的填充水平下对应的颗粒数目分别为500 00060 00015 0007 000,采用Hertz-Mindlin model接触模型模拟实际混合10 s,耗费的计算时间分别为5805463 h,在同一时间点截取不同颗粒粒径的模拟情况的图片,结果表明不同粒径的颗粒,在混合罐中具有相同的速度。Radete[27]将颗粒粒径分别设置为0.450.982.103.30 mm,分别将他们与自身粒径相同的颗粒混合,混合后计算相对标准偏差(RSD)和莱西指数(Lacey index),结果表明3.30 mm颗粒的RSDLacey指数曲线与其他3种颗粒不同,说明3.30 mm颗粒与其他3种小粒径颗粒的混合动力学不同。本实验尝试在相同接触参数设置条件下,将颗粒半径分别设置为0.51.01.5 mm,粒径分布设置为正态分布,标准差为0.05 mm进行模拟。结果发现颗粒半径为0.5 mm时不能在圆形底座上形成堆积;颗粒半径为1.0 mm时,得到的休止角为41.93°,与未设置粒径分布时的结果41.48°存在1.08%的相对误差;当颗粒半径为1.5 mm时,会堵塞漏斗下出口,颗粒无法下落。粒径是作为离散元模拟的输入参数,虽然对接触模型中各种参数的优选是有一定影响的,但通过标定实验,可以保证放大颗粒的计算结果精度满足要求。

6.2  接触参数的设置

MCC是制药过程中的常用辅料,Weis[28]在使用离散元模拟MCCvivapur® 102JRS PharmaGermany)与α-乳糖在滚圆过程的混合行为时,采用滞后弹簧模型(Hysteretic spring model)作为接触模型,将颗粒半径设为(0.600±0.047mm,密度1 406 kg/m3,将颗粒的恢复系数设置为0.096 9,滑动摩擦系数0.90,滚动摩擦系数0.05Kwan[29]采用仅考虑弹性接触的自定义接触模型[30]模拟了MCCFMC Corporation)与α-乳糖的碾磨行为,将MCC颗粒的密度设置为了1 696kg/m3杨氏模量为8 GPa,泊松比和摩擦系数分别设置为0.30.6。可见,同样属于纤维素类,上述文献中MCC接触参数的设置与本实验所得结果并不一致。在离散元模拟中,不同的研究对象、不同的颗粒粒径放大方式和形状简化方式、以及不同的接触模型选择,将导致所使用的参数种类和具体数值的不同,在实际应用中应具体问题具体分析。

6.3  仿真计算时间

计算机硬件配置是影响DEM计算时间的关键因素。本研究考察了两台不同配置的计算机用于离散元模拟:1.使用IntelRXeonRCPU E5-2620v42.10GHz服务器,但没有装配图形处理单元(Graphics ProcessingUnitGPU);2.使用AMD Ryzen 2400G 3.6GHz作为处理器并装配有AMD Radeon RX 580 2048SP作为GPU的个人电脑。配置1在颗粒数目较少时模拟速度较快,但是在颗粒数目增多以后速度变慢,对表3和表6每次实验进行模拟时耗时约20 h。配置2对表3和表6中每次实验进行模拟耗时约为13h。可见,选择具有高性能GPU有助于进行颗粒数量较多的模拟。

7  总结和展望

本实验在Hertz-Mindlin with JKR Cohesion接触模型和颗粒缩放的基础上,采用离散元方法模拟生甘草、独活、MCC和乙基纤维素4种粉末的休止角,并通过试验设计的方法筛选出对休止角模拟测定影响显著的接触参数,建立了接触参数与模拟休止角之间的回归模型,基于回归模型对每种粉末的接触参数进行了标定和验证。4种粉末休止角模拟结果与真实测量值的相对误差绝对值均小于2.0%,证明了所建标定方法的可靠性。本研究所建回归模型对休止角的标定范围为33.30°43.64°,在此范围内可通过建立的休止角预测方程实现关键接触参数的快速标定,节约了时间和成本。如需改变标定范围,可结合目标物料的流动性指标,改变试验设计的参数范围。

本实验通过宏观物性参数间接标定了4种代表性物料的DEM微观力学参数,可为这些物料在混合、输送等制药过程的精确模拟奠定基础。通常来说不是所有的接触参数都对粉末的宏观行为有影响,由于接触模型选择的不同以及模拟的工况不同,需要标定的接触参数也会发生变化,例如休止角标定试验通常是对粉末流动性的标定,然而在模拟颗粒的压缩时则需要使用单轴压缩试验进行标定。应用多种不同的模拟标定试验有助于获得更准确更稳健的接触参数[31],是制药过程DEM参数标定未来的研究方向。

参考文献(略) 

来  源: 石辰风,杨茂蕊,唐正馨,王  欣,徐  冰,乔延江.中药浸膏粉离散元模拟参数标定方法研究 [J]. 中草药, 2020, 51(24):6205-6212.

中药浸膏粉离散元模拟参数标定方法研究

中药浸膏粉离散元模拟参数标定方法研究

中药浸膏粉离散元模拟参数标定方法研究

生物医学科研方法

单细胞lncRNA分析,有它就够!

2020-12-31 3:15:13

生物医学科研方法

淡豆豉炮制过程中产纤溶酶微生物的筛选和鉴定

2020-12-31 3:45:51