粉体学是药剂学的理论基础, 对制剂处方设计、制备、质量控制、包装都有重要指导意义。固体制剂是各向异性的粉体物料(主药和药用辅料) 通过搅拌、混合、黏合或荷载等外力作用形成的一种颗粒聚集体, 包括散剂、颗粒剂和片剂等常用剂型^[1]。药用辅料是片剂制备必不可少的部分, 粉末流动性差会限制其应用, 因此良好的流动性对固体制剂具有重要影响。现有研究大多通过宏观现象描述粉体流动性, 休止角是粉体堆积自由表面与水平的夹角, 是一种快速、直观表征散体粉体流动能力的方法^[2]。粉体休止角的物理检测方法包括提升缸法、漏斗注入法、剪切盒法等^[3]。然而, 不同方法得到的数据不同, 且同一种方法测得的数据也不一样。休止角大小和测量方法的差异与粉体堆积过程的细观组构及演化规律有关。

离散单元法(discrete element method, DEM) 是一种用于分析单个颗粒在每一时间点上的运动规律与力学特征的数值模拟方法。通过DEM可获得一些难以通过物理实验检测到的微观信息^[4-6], 如每个颗粒的运动轨迹、速度和受力情况。因此, 采用DEM数值模拟, 以力链为核心的颗粒跨尺度分析获得物理实验方法难以获得的颗粒微观行为变化, 将有利于实现对粉体颗粒流动的系统认知。在DEM模拟过程中材料接触参数是决定数值模拟结果可靠性的关键。

DEM参数包括材料的本征参数(如泊松比、剪切模量、实体密度等) 和材料间接触参数(如碰撞恢复系数、静摩擦因数及滚动摩擦因数)。颗粒-颗粒静摩擦系数可通过斜面实验测定^[7], 直剪仪测定颗粒-不锈钢的滚动摩擦系数, 三轴压缩试验确定刚度比、有效模量等物性参数。然而, 这些实验仪器只适用于岩土力学的大颗粒, 药用辅料的粒径大部分都在微米级, 物性参数难以通过物理实验测定, 故通过仿真模拟反求为获取物性参数提供了有效的解决方法。

本文以微晶纤维素(microcrystalline cellulose, MCC) 粉体为研究对象, 利用PFC软件进行二维仿真模拟。首先, 以提升缸法休止角(θ) 为响应值, 通过Plackett-Burman、最陡爬坡及Box-Behnken试验, 确定最佳细观参数组合, 然后以漏斗注入法休止角(α) 和剪切盒法休止角(φ) 进行稳健性考察; 最后在确定可靠的细观参数基础上, 从颗粒堆积过程中力链的演变规律来分析不同测量方法以及同一测量方法人为因素出现的差异, 为正确解释休止角测量提供参考。这种新技术的探索将为固体制剂研究提供一种新思路, 可能取得的突破将为制药工业的数字化和智能化提供理论依据和数据支撑。

药品与试剂  MCC (SH200, 批号200301, 安徽山河药用辅料股份有限公司惠赠), 经测定, 50~100、100~200、200~300 μm粒径所占的体积百分数分别为10%、60%、30%, 平均粒径为180 μm。

仪器  激光粒度仪(Malvern2000, 英国马尔文公司); 休止角测定仪(漏斗注入法, HYL-105, 丹东百特仪器有限公司; 提升缸和剪切盒, 委托重庆力耘喷嘴有限公司设计开发)。

PFC软件  PFC 5.0是基于DEM框架, 由计算机引擎和图形用户界面构成的细观分析软件。通过构建离散单元代码允许离散颗粒产生位移和旋转, 随着计算过程自动识别新的接触从而模拟颗粒的运动和相互作用。PFC集成了二维及三维, 在二维软件里将颗粒视为圆盘, 三维里表现为球体; 墙体在二维视为线段, 三维表现为三角形, 并可构成任意复杂的多边形墙体。此软件还包含C++插件选项, 用户可根据自身需求自定义接触模型, 完善软件功能。

休止角测定

提升缸法θ (图 1A-1) 空心缸(37.8 mm × 90 mm, d × h) 垂直放置于直径D为80 mm的圆盘台正中心, 将一定质量粉体填充于空心缸内, 缓慢提起缸, 粉体堆满于正下方的圆盘台上形成圆锥体堆, 测量并记录粉堆高度H₁, 按公式θ = arctan (2H₁/D) 计算该休止角θ, 重复3次, 求平均值。

Figure 1

Figure 1

Figure 1 Physical and simulation model of the angle of repose. A: Lifting cylinder test; B: Funnel flowing test; C: Shear box test; 1: Physical model; 2: Simulation model

漏斗注入法α (图 1B-1) 加入足量待测样品于漏斗中, 使其流下并逐渐堆积于正下方的样品圆盘台, 形成对称的圆锥体堆, 且平台四周有粉体落下时停止加料, 自带量角器测量并记录休止角α, 重复3次, 求平均值。

剪切盒法φ (图 1C-1) 将一定质量的样品倒入(44 mm × 42 mm × 47 mm, l × w × h) 无顶盖的盒中, 撤掉活动挡板, 样品在重力作用下流出, 待样品静止形成斜坡, 测定其高度H₃和宽度W₁, 按φ = arctan (H₃/W₁) 计算该休止角, 重复3次, 求平均值。

休止角仿真模拟

离散元基本原理  DEM是由Cundall等^[8]于20世纪70年代末首次提出, 其基本原理建立在牛顿第二运动定律上, 结合力与位移和相应的本构模型来描述介质的复杂力学行为, 通过时步迭代的方法求解不连续体的运动过程。假设颗粒i和颗粒j之间存在法向和切向的刚度K_n、K_s, 法向和切向的阻尼系数η_n、η_s及摩擦系数μ, 颗粒单元的接触力通过单元之间的“重叠”而产生, 通过力-位移定律将接触力与相对位移联系起来, 实现动态模拟过程。

抗滚动阻尼线性模型  接触本构模型体现了颗粒间力与力矩的相互作用, 表征颗粒单元之间的拉应力或黏滞力, 因此, 针对不同的颗粒性质、运动状态选择合适的接触本构模型。本研究中的MCC粒径较小(50~300 μm), 表面粗糙, 考虑到软件自带的抗滚动阻尼线性模型适用于非规则材料^[9], 为了使模拟更接近真实试验, 故选择抗滚动阻尼线性模型。抗滚动阻尼线性模型的主要原理即通过增加滚动阻尼系数, 增加颗粒转动时的内力矩从而降低转动能力^[10], 实现颗粒堆积, 同时滚动摩擦系数大小也可用来模拟颗粒形状^[11]。本研究基于PFC 2D软件, 为降低计算规模, 选取模具轴截平面进行二维仿真^{[4, 12]}。抗滚动阻力线性模型的力-位移定律将接触力和力矩定为[公式(1~3)]:

式中, F_c、F^l、F^d分别为颗粒间的接触力、线性接触力、阻尼力(N); M_c、$ {M}^{\mathrm{r}} $分别为接触力力矩、抗滚动阻尼力矩(N·m)。

式中, K_s为切向刚度(N·m^-1); R^*为有效接触半径(m), $ \frac{1}{{R}^{\mathrm{*}}}=\frac{1}{R\left(A\right)}+\frac{1}{R\left(B\right)} $; Δθ_r为相对旋转增量。

仿真参数  抗滚动阻尼线性模型需设定的参数为有效模量、刚度比、颗粒-颗粒摩擦系数、颗粒-不锈钢摩擦系数和滚动阻尼系数, 其中MCC密度为1 600 kg·m^-3。由于PFC软件常用于岩石等大颗粒的离散元模拟, 粉末材料的参数甚少, 结合文献^{[13, 14]}, 初始孔隙率设置为0.3。对抗滚动阻尼线性模型的微观参数进行初步赋值, 有效模量为1×10⁶ kPa, 刚度比为2, 颗粒-颗粒摩擦系数为0.5, 颗粒-不锈钢摩擦系数为0.5, 滚动阻尼系数为0.5。

仿真模拟   为了提高仿真速率, 本研究将样品近似看成球(圆盘)。首先在指定尺寸的模具内按照粒径分布要求用distribute方式生成若干圆盘颗粒, 对颗粒及模具设置的细观参数, 指定重力加速度g = 9.81 m·s^-2, 采用重力沉降法, 使颗粒自由下落, 待系统动能趋于0时, 分别模拟不同方法测量休止角。其中, 提升缸法以一定速度提升空心缸, 漏斗注入法打开漏斗下端口释放颗粒, 剪切盒法撤掉侧板, 分别使颗粒流落或堆积在底盘上(图 1A-2、B-2、C-2), 颗粒动能耗散且相对误差小于0.005%时认为颗粒堆积停止, 导出图片, 采用CAD距离测定功能, 测定高度计算休止角大小。为了消除粒径效应对模拟结果的影响, 在保持样品粒径大小尺寸的前提下, 以θ测量法进行模具几何尺寸缩放, 确定缩放倍数, 并用于后续仿真模拟。

试验设计  Plackett-Burman试验设计通过考察目标值与各因子间的关系, 比较各因子2水平间的差异确定因子显著性。本文以θ为响应值, 根据多次调试, 将初步赋值设为低水平, 高水平为低水平(初步赋值) 的2倍, 从有效模量(X₁)、刚度比(X₂)、颗粒-颗粒摩擦系数(X₃)、颗粒-不锈钢摩擦系数(X₄)、滚动阻尼系数(X₅) 中筛选出对θ影响显著的参数, 再通过最陡爬坡设计, 即以结果变化的梯度方向为爬坡方向, 根据各因素影响值的大小确定变化步长, 快速有效地确定关键接触参数的最佳区域, 然后根据Box-Behnken试验建立接触参数和休止角之间的回归模型。

参数优化与验证  优化回归模型, 应用Design Expert 12软件数值求解优化功能, 根据目标值, 得到最佳参数组合, 并使用该组参数模拟提升缸法休止角, 比较模拟值与真实值的差异, 验证参数组合的准确性。

参数稳健性考察  基于优化后的最佳参数组合进行模拟, 进行漏斗注入法和剪切盒法进行仿真模拟, 比较模拟值与真实值的差异, 考察参数组合的稳健性。

颗粒细观分析  在颗粒体系内部, 由于自身重力和外部荷载作用, 颗粒物质在颗粒界面接触处产生接触力, 接触力沿着颗粒接触网络传递, 并最终形成完整的力链网络^[15]。力链作为连接宏观、细观特性之间的桥梁, 通过力链能进一步加强对颗粒流动的宏观认识, 因此, 本研究利用PFC 2D中的fish语言输出力链的接触数目、大小, 提取3种方法的休止角力链演化图, 从力链的形成、大小及数量等细观角度^[16]分析不同休止角测量方法及同一测量方法的人为因素产生差异的原因。

数据分析  采用SPSS、Auto CAD、Design-Expert软件进行试验设计和数据分析, P < 0.05表示有显著性差异。

本文分别对提升缸法模具的整体尺寸分别缩小2、2.5、3、4、6、12倍, 结果表明, 当模型缩小至原尺寸4倍及以内时, 根据SPSS回归数据分析P = 0.843 > 0.05, 休止角无显著变化; 当模型缩小至原尺寸的4倍以上时, 休止角大小随模具缩小倍数呈增大趋势, 因此, 当模型缩小4倍及以内可忽略粒径效应。本研究以θ测量法进行模具几何尺寸缩放, 确定缩放倍数用于后续仿真模拟(图 1A-2、B-2、C-2)。

Placket-Burman试验设计及结果如表 1所示。对结果进行方差分析, 可知X₅的P = 0.008 < 0.01, 对休止角影响极其显著, X₃、X₄的P = 0.025、0.039 < 0.05, 对休止角影响显著, 而X₁、X₂的P = 0.109、0.354 > 0.05, 对休止角影响极小。X₁与粉末的弹性模量有关, X₂与粉末的泊松比有关, 但对粉末类的解释甚少且对休止角影响极小, 本文的参数筛选也证明了这一点, X₁、X₂对休止角无显著影响, 故根据文献资料^[10]固定X₁为1×10⁶ kPa、X₂为2进行最陡爬坡及响应面设计。

表1(Table 1)

Table 1 Design and response value of Plackett-Burman test (X1: Effective modulus, X2: Stiffness ratio, X3: Particle-particle coefficient of friction, X4: Particle-stainless steel coefficient of friction, X5: Rolling damping coefficient) No. X1 X2 X3 X4 X5 Angle of repose θ/° 1 1 -1 -1 -1 1 34.75 2 -1 1 1 1 -1 36.44 3 -1 1 1 1 1 35.69 4 -1 -1 -1 -1 1 36.42 5 0 0 0 0 0 37.53 6 -1 -1 -1 -1 -1 33.15 7 1 -1 1 1 -1 34.00 8 1 1 1 -1 -1 33.47 9 -1 1 1 -1 1 36.87 10 -1 -1 -1 1 -1 35.56 11 1 -1 1 1 1 38.77 12 1 1 -1 -1 -1 31.84 13 1 1 -1 1 1 34.94

Table 1 Design and response value of Plackett-Burman test (X1: Effective modulus, X2: Stiffness ratio, X3: Particle-particle coefficient of friction, X4: Particle-stainless steel coefficient of friction, X5: Rolling damping coefficient)

根据Placket-Burman试验筛选出的参数进行最陡爬坡试验, 以便快速找到最接近真实值的区域。本研究最陡爬坡试验设计及结果如表 2所示(多次调试, Ⅰ、Ⅱ为调试试验, Ⅲ为正式试验), 调试试验结果的相对误差无法达到最接近真实值的效果, 故最终拟定X₅的爬坡步长定为0.1, X₃数的爬坡步长定为0.05, X₄的爬坡步长定为0.05。MCC休止角的实测值是37.29°, 结果表明14号的相对误差最小, 13~15号的相对误差先减小再增大, 因此选14号为中水平, 13号为低水平, 15号为高水平进行后续的响应面设计。

表2(Table 2)

Table 2 Design and results of steepest ascent test Serial number X3 X4 X5 Angle of repose θ/° Relative error/% Ⅰ 1 0.5 0.5 0.5 33.15 11.11 2 0.6 0.6 0.57 36.21 1.55 3 0.7 0.7 0.64 35.73 2.85 4 0.8 0.8 0.71 36.10 1.85 5 0.9 0.9 0.78 35.17 4.38 Ⅱ 6 0.5 0.5 0.5 33.15 11.11 7 0.55 0.55 0.54 34.59 5.95 8 0.6 0.6 0.58 34.63 5.85 9 0.65 0.65 0.62 36.4 1.03 10 0.7 0.7 0.66 35.33 3.94 Ⅲ 11 0.5 0.5 0.5 33.15 11.11 12 0.55 0.55 0.6 34.92 6.36 13 0.6 0.6 0.7 36.3 2.65 14 0.65 0.65 0.8 37.35 0.05 15 0.7 0.7 0.9 38.1 2.17

Table 2 Design and results of steepest ascent test

根据最陡爬坡试验设计筛选出的因素, 对X₃、X₄、X₅进行Box-Behnken试验设计(表 3), 对自变量和响应量进行二次多项拟合, 得到回归方程(θ = 37.42 + 0.221 2X₃ + 0.342 5X₄ + 0.733 7X₅ - 0.042 5X₃X₄ - 0.08X₃X₅ + 0.592 5X₄X₅ + 0.468 3X₃² + 0.051 8X₄² - 0.011 7X₅²)。

表3(Table 3)

Table 3 Design and results of Box-Behnken test No. X3 X4 X5 Angle of repose θ/° 1 1 0 0 38.76 2 0 -1 -1 37.25 3 0 1 1 39.02 4 -1 1 1 38.1 5 -1 0 0 38.52 6 0 0 0 37.35 7 0 1 1 36.42 8 1 -1 -1 37.8 9 1 1 1 38.5 10 0 -1 -1 37.02 11 0 0 0 37.54 12 0 0 0 37.38 13 -1 -1 -1 37.23 14 -1 0 0 36.84 15 1 0 0 37.4

Table 3 Design and results of Box-Behnken test

模型方差分析结果所得拟合模型P < 0.000 1, 失拟项P_L = 0.555 3 > 0.05, 表明模型良好, 没有弯曲失拟现象发生, 试验中变异系数CV = 0.27%, 说明试验有较好可靠性。决定系数R² = 0.994; 校正决定系数R_adj² = 0.982; 预测决定系数R_pre² = 0.935; 三值都 > 0.9, 表明模型能真实地反映实际情况。

基于回归模型, 将对休止角影响不显著(P > 0.05) 的因素(X₃X₄、X₃X₅、X₅²、X₄²) 移除, 得到优化后的二次多项式回归方程(θ = 37.43+ 0.221 2X₃ + 0.342 5X₄ + 0.733 7X₅ + 0.592 5X₄X₅ + 0.468 3X₃²)。

优化后的P < 0.000 1, 失拟项P_L = 0.630 1 > 0.05, 决定系数R² = 0.989; 校正决定系数R_adj² = 0.984; 预测决定系数R_pre² = 0.968; 三值都 > 0.9, 表明优化后的模型仍具有良好可靠性。应用数值求解优化功能, 以θ平均实测值37.29° (n = 3) 为目标, 得到MCC的接触参数X₃为0.6、X₄为0.6、X₅为0.79, 得到θ的预测值为37.295°。使用该组参数进行仿真模拟, 得到模拟的θ为37.13°, 预测值与模拟值的相对误差为0.44%, 表明该参数组合可靠。

为进一步验证参数组合的稳健性, 分别采用3种检测方法进行仿真模拟, 结果表明(图 2), 仿真模拟值与真实试验值无显著差异, 说明该参数组合稳健性好。不同方法测定结果有所差异, θ、φ大于α。可能原因是漏斗注入法测量时, 粉体受到重力加速度影响, 堆积过程中将分解额外剪切应力, 所以θ、φ大于α。但由于缸体采用不锈钢304材质, 摩擦力较小, 所以θ与φ差异不显著。

Figure 2

Figure 2

Figure 2 Experimental and simulated values of angle of repose measured by different methods. θ: Lifting cylinder test; α: Funnel flowing test; φ: Shear box test

一般认为, 休止角的变化在2°之内, 符合休止角允许的误差范围(差值约为2°)^[4], 提升缸法中提升缸速度和漏斗注入法的漏斗高度为休止角测量的人为变量, 图 3为过程变量对休止角测量大小的影响, 结果表明, 提升速度对休止角的大小有一定影响(图 3A), 当提升速度小于0.25 m·s^-1时, 对休止角的测量结果影响不大, 但当速度大于0.5 m·s^-1时, 对测量结果有较大影响, 可能是因为速度较小时, 颗粒与墙体之间有充分的接触并受到摩擦缓慢堆积, 而速度大于0.5 m·s^-1时, 由于速度过快, 颗粒还未开始堆积已脱离墙体, 受惯性作用, 颗粒坍塌角度偏小。当漏斗高度在8.25~35 mm内(图 3B), 对测量结果无显著影响, 当高度为70 mm时(相当于常规高度的4倍), 对测量结果影响较大, 可能是因为高度越高, 颗粒碰撞到地面的速度越大, 颗粒间相互作用力增大使其坍塌, 角度偏小。

Figure 3

Figure 3

Figure 3 Effect of lifting velocity (A) and funnel height (B) on angle of repose

休止角的形成可用力链结构进行解释, 在重力作用下, 粉体颗粒间相互挤压形成与压力一致方向的力链, 传递重力, 球形颗粒在接触点处轻微变形, 可支撑与轴线方向有一定夹角的剪切应力, 当剪切应力超过一定数值时, 则粉堆内部强力链结构断裂, 粉堆表面上力链末端颗粒首先发生滑动, 同时新力链不断形成, 直到剪切应力小于该临界数值, 休止角保持恒定。

为进一步揭示休止角测量方法导致测量结果差异的内在原因, 将模拟时间平均分成10个阶段, 分别提取每一段的力链总数, 结果发现(图 4), 提升缸在提升过程和剪切盒在撤板过程, 颗粒体系在外载荷下发生剧烈运动, 颗粒体系内部的力链不断断裂重组, 因此力链数波动较大, 后期颗粒堆积最终达到平衡, 力链趋于稳定; 漏斗注入法无外载荷作用, 颗粒在初期尚未生成力链, 随颗粒不断堆积, 颗粒体系内部的力链数目不断增加, 且变化较快, 因为重力加速度的作用, 颗粒撞击后以一定速度向周围崩塌, 导致休止角偏小。

Figure 4

Figure 4

Figure 4 Variation law of the number of internal force chains in the particle system during the simulation of the angle of repose

同时, 为进一步解释人为因素导致测量结果的影响原因, 本研究分别提取θ法和α法过程变量在4个不同时间点的力链网络变化情况(图 5)。图 5A、B为提升缸速度分别为0.05和0.5 m·s^-1颗粒力链网络变化, 从4个时间点看, 提升速度主要影响起始状态力链的分布。1/4时间时(图 5A-1、B-1), 当速度较慢时, 初始位置的力链还没来得及分散, 而0.5 m·s^-1的休止角雏形已成, 强力链主要集中在下方且网络空隙较小, 0.05 m·s^-1下方的力链方向性较0.5 m·s^-1更剧烈; 2/4时间点时(图 5A-2、B-2), 随颗粒的持续放出, 颗粒沿力链方向滚落, 堆积表面的力链网络空隙逐渐减小, 存在局部应力集中的现象; 3/4及4/4时(图 5A-3、A-4、B-3、B-4), 力链趋于稳定并最终达到平衡状态, 力链数分别为4 866、4 420, 初始状态力链数目、大小分布一致, 证明了提升速度越大, 最终力链数越小, 颗粒越少, 角度越小。图 5C、D是漏斗高度分别为17.5和70 mm颗粒力链网络变化, 1/4时间点时(图 5C-1、D-1), 7.5 mm的休止角底部生成了较强的力链且聚集, 呈现方向性, 而70 mm的底部生成不稳定的弱力链分散在地面, 2/4时(图 5C-2、D-2), 70 mm的力链方向较平缓, 颗粒体沿方向滚落, 故二者相比, 70 mm的角度偏小; 3/4到4/4 (图 5C-3、C-4、D-3、D-4), 随颗粒体的堆积, 堆积体表面的力链网络空隙逐渐减小, 且17.5 mm的表面顶端力链较凸出, 而70 mm的整个堆积体表面较平滑, 休止角明显偏小, 可能是由于漏斗高度过高, 颗粒体受自身重力作用自由下落, 与颗粒、底盘发生碰撞速度过大而分散, 堆积困难, 进而导致堆积体表面的颗粒坍塌, 角度偏小。

Figure 5

Figure 5

Figure 5 Effect of the lifting speed and the funnel height on angle of repose. A: 0.05 m·s-1; B: 0.5 m·s-1. A-1, A-2, A-3 A-4: 1/4, 2/4, 3/4, 4/4 period of the whole stacking process; B-1, B-2, B-3, B-4: 1/4, 2/4, 3/4, 4/4 period of the whole stacking process; C: 17.5 mm; D: 7 m·s-1. C-1, C-2, C-3 C-4: 1/4, 2/4, 3/4, 4/4 period of the whole stacking process; D-1, D-2, D-3, D-4: 1/4, 2/4, 3/4, 4/4 period of the whole stacking process

DEM包含多种接触模型, 不同接触模型的参数设置有所不同。石辰风等^[17]采用DEM模拟中药浸膏粉休止角, 运用Hertz-Mindlin with JKR Cohesion接触模型, 标定的MCC颗粒-颗粒的滚动摩擦系数为0.26; Mohajeri等^[18]采用自定义的接触模型得到滚动摩擦系数为0.6。本研究采用抗滚动阻尼线性模型得到的滚动摩擦系数为0.6, 与Mohajeri等^[18]研究结果一致, 但与石辰风等^[17]研究有差异, 说明接触模型的选择对模拟结果有影响。

固体制剂中的物料粒径往往在微米级, 而粒径越小计算机模拟效率越低。因此, 为了提高模拟效率, 有必要对颗粒或设备进行缩放。Mohajeri等^[18]对剪切盒模型进行缩放的倍数对模拟休止角有一定的影响; Jiang等^[13]对圆筒直径的选取应大于颗粒最大粒径的4倍以上可忽略粒径效应影响。而本研究发现圆筒直径大于颗粒粒径的30倍即可忽略粒径效应的影响, 可能是因为前人研究的是单一粒径的颗粒, 而本研究按照粒径分布生成的颗粒, 差异较大。同时粒径分布对模拟结果也存在差异, Jiang等^[13]在相同参数下, 对平均粒径35~50 μm的颗粒进行均匀分布、单粒径分布、高斯分布并模拟, 结果表明单粒径分布的颗粒的休止角小于其他两种。本文仅采用粒径均匀分布, 只有粒径分布模拟越接近真实颗粒, 标定的参数才越可靠。

休止角是最快速、最直观的检测方法, 然而测定方法不同, 结果也有所差异。张昱等^[3]采用漏斗法和剪切盒法, 测量沙粒光滑玻璃底盘堆积, 结果两种休止角无显著差异。而本文测得结果存在差异, 可能与颗粒本身的属性如粒径大小以及形态有关。

本研究采用PFC抗滚动阻尼线性模型, 模拟二维休止角。通过试验设计筛选出对休止角影响显著的因素及最佳参数组合, 并以此模拟了3种方法(提升缸法θ、漏斗注入法α、剪切盒法φ) 的休止角, 验证该参数组合具有一定的稳健性。从颗粒细观角度诠释了不同操作和测量方法对休止角存在差异性的原因, 可为物料物性参数的标定和下一步固体制剂研究提供理论参考。以上实验结果及分析还需其他物理模型(直剪实验) 进一步验证提高其适用性, 该标定方法尚且只能得出一种辅料的参数, 后续可对此进行延伸, 改进标定方法, 从而获得多种辅料的参数, 建立系统的药物辅料数据库。

作者贡献: 戴传云和孙会敏是本课题负责人, 负责提出课题思路、设计研究方案和论文修订; 谢文影是本课题主要执行人, 完成论文数据, 撰写论文初稿; 白玉菱、赵孟涛、周康明、范仁宇协助实验实施; 管天冰和任建兵参与部分实验设计。

利益冲突: 所有作者均同意投稿且不存在与本文相关的利益冲突。