Loading-unloading contact law for micro-crystalline cellulose particles under large deformations

TL;DR

This paper introduces a semi-empirical contact law for micro-crystalline cellulose particles that models large deformations, capturing elastic-plastic behavior and strain-hardening effects through experimental calibration.

Contribution

It develops a generalized contact law for irregular particles incorporating shape, elastic, and plastic properties, validated against experimental force-displacement data.

Findings

Good agreement between model predictions and experimental data.

Plastic properties follow a log-normal distribution.

The model captures strain-hardening behavior during compression.

Abstract

A semi-empirical mechanistic contact law for micro-crystalline cellulose (Avicel PH-200) particles is proposed and characterized experimentally using force-displacement curves obtained from diametrical compression of single particles. The concepts of a shape factor and a master contact law are introduced first for elastic ellipsoidal particles, and subsequently generalized to plastic irregular particles. The proposed loading-unloading contact law is a function of three characteristic diameters (lengths of the principal axes of an approximated ellipsoid), a geometric parameter associated with the loading condition, three plastic and one elastic material properties. The force-displacement curves obtained using a micro-compression tester exhibit an apparent strain-hardening at distinctly different strain values, which is captured by the shape factor function and its geometric parameter. The three plastic material properties are log-normal distributions estimated from the loading experimental curves, while the elastic property is estimated from the unloading experimental curves. The study shows a very good agreement between predictions of the calibrated loading-unloading contact law and the experimental values.

Figures

Tables

Table 1. Table 1: Obtained characteristic dimensions of the studied Avicel PH-200 particles

Particle No.	$D_1\left(\mu m\right)$	$D_2\left(\mu m\right)$	$D_3\left(\mu m\right)$
$1$	$213.54$	$196.07$	$186.33$
$2$	$155.32$	$217.51$	$136.87$
$3$	$168.74$	$171.64$	$125.15$
$4$	$97.56$	$100.98$	$69.09$
$5$	$238.47$	$249.24$	$237.81$
$6$	$93.8$	$76.46$	$52.85$
$7$	$160.44$	$167.7$	$151.3$
$8$	$114.1$	$114.54$	$80.98$
$9$	$168.31$	$144.54$	$131.27$
$10$	$133.8$	$117.12$	$103.38$
$11$	$122.69$	$131.73$	$109.91$
$12$	$121.04$	$205.93$	$133.08$
$13$	$185.1$	$206.98$	$146.91$
$14$	$180.35$	$173.88$	$139.81$
$15$	$173.26$	$157.3$	$147.79$
$16$	$296.34$	$267.76$	$219.78$
$17$	$219.5$	$281.94$	$225.55$
$18$	$194.97$	$167.95$	$151.23$
$19$	$105.57$	$109.25$	$86.46$
$20$	$154.12$	$182.82$	$166.13$

Table 2. Table 2: Plastic coefficients log ⁡ ( C i ) ∼ 𝒩 ​ ( μ C i , σ C i 2 ) similar-to subscript 𝐶 𝑖 𝒩 subscript 𝜇 subscript 𝐶 𝑖 subscript superscript 𝜎 2 subscript 𝐶 𝑖 \log(C_{i})\sim\mathcal{N}(\mu_{C_{i}},\sigma^{2}_{C_{i}}) of the stress-strain relationship of the master contact law for micro-crystalline cellulose, Eqn. ( 16 ).

	${\mu }_{C_i}$	${\sigma }_{C_i}$
$\log \left(C_1\right)$	$1.9027$	$0.2985$
$\log \left(C_3\right)$	$1.1015$	$4.8914\times {10}^{-9}$
$\log \left(C_5\right)$	$0.9784$	$0.0797$

$a$	geometric factor for plastic ellipsoidal particles (-)
$c_{max}^2$	material parameter for similarity contact law (-)
$C_1$	plastic material property 1 ($\mathrm{MPa}$)
$C_2$	plastic material property 2 (${\mathrm{MPa}}^{-3}$)
$C_3$	plastic material property 2 (${\mathrm{MPa}}^{-5}$)
$D_1$	major diameter of the ellipsoidal particle when looked through MCT’s microscope ($\mathrm{mm}$)
$D_2$	minor diameter of the ellipsoidal particle when looked through MCT’s microscope ($\mathrm{mm}$)
$D_3$	vertical diameter of the ellipsoidal particle captured by MCT’s side-view camera ($\mathrm{mm}$)
$E$	Young’s modulus ($\mathrm{MPa}$)
$E[\cdot ]$	expected value (-)
$F$	force acting on the particle ($N$)
$R_e$	effective radius of the ellipsoidal particle ($\mathrm{mm}$)
$m$	Power-law hardening exponent (-)
$N$	total number of tested particles (-)
$s_{ϵ}$	standard deviation of a set of strain values (-)
$\#{𝒮}_{\sigma }$	number of experiments for which a given stress $\sigma$ is reached (-)
$\gamma$	particle deformation at the contact ($\mathrm{mm}$)
$ϵ$	effective particle strain (-)
$\overline{ϵ}$	mean of the strain values for a set of particles (-)
${ϵ}_p$	internal strain variable for unloading contact law (-)
${\eta }_e$	elastic shape factor (-)
${\eta }_p$	plastic shape factor (-)
$\kappa$	particle reference strength ($\mathrm{MPa}$)
$\nu$	Poisson’s ratio (-)
$\sigma$	effective particle stress ($\mathrm{MPa}$)
${\sigma }^m$	maximum stress value in the master contact law ($\mathrm{MPa}$)
${\sigma }_y$	uniaxial yield stress ($\mathrm{MPa}$)
${\sigma }_p$	internal stress variable for unloading contact law ($\mathrm{MPa}$)

Equations

$$F=\frac{4E\sqrt{R_e}}{3(1-{\nu }^2)}{\left[{\eta }_e(D_1,D_2)\gamma \right]}^{3/2}$$

$$\frac{1}{{\eta }_e(D_1,D_2)}=1-{\left[{\left(\frac{D_1}{D_2}\right)}^{0.1368}-1\right]}^{1.531}$$

$$\sigma =\frac{F}{\pi R_e^2}$$

$$ϵ=\frac{\gamma }{2R_e}$$

$$\sigma :=\sigma ({\eta }_eϵ)=\frac{8\sqrt{2}E}{3\pi (1-{\nu }^2)}{\left[{\eta }_e(D_1,D_2)ϵ\right]}^{3/2}$$

$$\sigma :=\sigma ({\eta }_pϵ)$$

$${\eta }_p:={\eta }_p(D_1/D_3,D_2/D_3)\text{mbox}suchthat{\eta }_p(1,1)=1$$

$${\sigma }_i(ϵ)=\sigma ({\eta }_{p,i}ϵ)⟹\sigma (ϵ)={\sigma }_i(ϵ/{\eta }_{p,i})$$

$$ϵ(\sigma )={\eta }_{p,i}{ϵ}_i(\sigma )$$

$$ϵ(\sigma )\sim \mathcal{N}(\overline{ϵ}(\sigma ),s_{ϵ}^2(\sigma ))$$

$$\overline{ϵ}(\sigma )$$

$$s_{ϵ}^2(\sigma )$$

$${\eta }_p(D_1/D_3,D_2/D_3):=1+a\left(1-\frac{D_3^2}{D_1{D}_2}\right)$$

$${\sum }_{i=1}^N\frac{{\int }_0^{{\sigma }_{\sigma ,i}^{\mathrm{m}}}{\left[\frac{1}{\#{\mathcal{S}}_{\sigma }}{\sum }_{j\in {\mathcal{S}}_{\sigma }}{\eta }_{p,j}{ϵ}_j(\sigma )-{\eta }_{p,i}{ϵ}_i(\sigma )\right]}^2\mathrm{d}\sigma }{{\int }_0^{{\sigma }_{\sigma ,i}^{\mathrm{m}}}{\left[\frac{1}{\#{\mathcal{S}}_{\sigma }}{\sum }_{j\in {\mathcal{S}}_{\sigma }}{\eta }_{p,j}{ϵ}_j(\sigma )\right]}^2\mathrm{d}\sigma }$$

$$\frac{1}{\#{\mathcal{S}}_{\sigma }}\sum _{i\in {\mathcal{S}}_{\sigma }}{\eta }_{p,i}{ϵ}_i(\sigma )\approx {\eta }_{p,N}{ϵ}_N(\sigma )$$

$$a:=\arg \underset{a}{\min }{\sum }_{i=1}^{N-1}\frac{{\int }_0^{{\sigma }_{\sigma ,i}^{\mathrm{m}}}{\left[{ϵ}_N(\sigma )-\frac{{\eta }_{p,i}}{{\eta }_{p,N}}{ϵ}_i(\sigma )\right]}^2\mathrm{d}\sigma }{{\int }_0^{{\sigma }_{\sigma ,i}^{\mathrm{m}}}{\left[{ϵ}_N(\sigma )\right]}^2\mathrm{d}\sigma }$$

$$\sigma (ϵ)/\mathrm{MPa}=C_1{\eta }_pϵ-{\left(C_3{\eta }_pϵ\right)}^3+{\left(C_5{\eta }_pϵ\right)}^5$$

$$\mathrm{E}⟨\sigma (ϵ)⟩/\mathrm{MPa}={\mathrm{e}}^{{\mu }_{C_1}+{\sigma }_{C_1}^2/2}{\eta }_pϵ-{\mathrm{e}}^{3{\mu }_{C_3}+9{\sigma }_{C_3}^2/2}{\left({\eta }_pϵ\right)}^3+{\mathrm{e}}^{5{\mu }_{C_5}+25{\sigma }_{C_5}^2/2}{\left({\eta }_pϵ\right)}^5$$

$$\mathrm{E}⟨\sigma (ϵ)⟩/\mathrm{MPa}=7.0094{\eta }_pϵ-27.235{\left({\eta }_pϵ\right)}^3+144.23{\left({\eta }_pϵ\right)}^5$$

$$\sigma (ϵ)=4(3^{1-1/m})(c_{max}^2{)}^{1+1/2m}\kappa {ϵ}^{1+1/2m}$$

$$\kappa ={\sigma }_y^{1-1/m}{E}^{1/m}$$

$$c_{max}^2=1.43\exp \left(\frac{-0.97}{m}\right)$$

$$\sigma (ϵ)=17.16{\sigma }_yϵ$$

$$\sigma(\epsilon)=\left\{\begin{array}[]{l}\left[C_{1}~{}\eta_{p}\epsilon-\left(C_{3}~{}\eta_{p}\epsilon\right)^{3}+\left(C_{5}~{}\eta_{p}\epsilon\right)^{5}\right]\mbox{MPa}\hskip 65.04256pt\mbox{plastic loading}\\ \frac{2}{\pi}\sigma_{\mbox{\tiny P}}\left[\arcsin\left(\phi(\epsilon,\epsilon_{\mbox{\tiny P}})\right)-\phi(\epsilon,\epsilon_{\mbox{\tiny P}})\sqrt{1-\phi(\epsilon,\epsilon_{\mbox{\tiny P}})^{2}}\right]\hskip 25.29494pt\mbox{elastic (un)loading}\end{array}\right.$$

$$\varphi (ϵ,{ϵ}_{\text{mbox}P})={\left[1-4.12\bar{E}\frac{({ϵ}_{\text{mbox}P}-ϵ{)}^2{ϵ}_{\text{mbox}P}}{{\sigma }_{\text{mbox}P}}\right]}_{+}^{1/2}$$

Full text

1 Introduction

Over the last decade, pharmaceutical industry has undergone a paradigm shift (Ierapetritou et al., 2016) towards adoption of efficient process, control and material usage strategies in manufacturing, commonly referred to as Quality by Design (QbD) (Yu et al., 2014; Etzler & Uddin, 2013), which provide guidelines for optimal control over quality of the final product by acquiring a detailed understanding of the effect of process parameters on raw material attributes through (semi) empirical models. These strategies are further complemented by Quality by Control (QbC) (Su et al., 2018b) techniques, which are more prevalent in the continuous manufacturing practices as they provide increased flexibility to the design space and efficiency in characterization processes (Su et al., 2017, 2018a).

Understanding the connection between mechanical behavior of a single particle and tabletability of the powder, i.e., the ability of the powder to gain strength under confined pressure, is essential for realizing QbD and QbC objectives (Šimek et al., 2017; Yi et al., 2018). An alternative approach utilizing empirical models has been extensively attempted, albeit with limited success, to try to capture the stress-strain relationship of a powder bed under load. Heckel and Kawakita (Denny, 2002; Šimek et al., 2017) models are amongst the ones commonly used to relate effective compaction parameters with deformation mechanisms and material properties of single particles. In general, these empirical relationships cannot be interpreted unambiguously, unless a plethora of empirical compaction equations is selectively used (Nordström et al., 2012). Similar empirical efforts have been made, with rather limited success, to relate effective strength parameters of the compact with particle material properties, e.g., by means of the Leuenberger equation (Kuentz & Leuenberger, 2000).

The challenges in achieving a systematic experimental characterization of micro-sized particles under confined conditions and large deformations are associated to the small size and irregular morphology of the particles, typical of pharmaceutical powders. Diametrical compression of single micro-particles using indentation equipment or Micro-Compression Testers (MCT) is one of the most common confinement techniques used to characterize the mechanical properties of various materials, especially polymers (Egholm et al., 2006; He et al., 2009; Tanaka et al., 2007b; Liu et al., 1998). Tanaka et al. (2007a) used a micro compression equipment to study the compaction behavior of individual 50 micron alumina granules with PVA (Poly Vinyl Alcohol) or PAA (Poly Acrylic Acid) as the binder matrix. In the context of pharmaceutical materials, Yap et al. (2008) studied the diametrical compression of single particles of various pharmaceutical excipients to relate the mechanical properties characterized at particle-scale to bulk compression properties of the materials. The technique has also been used to characterize visco-elastic properties of agarose (Yan et al., 2009) and alginate (Nguyen et al., 2009; Wang et al., 2005) microspheres. A major recent development in confinement methods for single particles has been made by Jonsson et al. (2015), who have developed a novel apparatus for triaxial testing of single particles, providing a more realistic insight into the behavior of individual particles under confinement. The apparatus has been used to investigate and characterize the yield and fracture behavior of micro-crystalline cellulose (MCC) granules (Jonsson & Frenning, 2016), highlighting the critical differences in their fracture behavior during triaxial and uniaxial compression, and correlating the triaxial response of single particles with bulk powder compression response.

In most of the mentioned works, the shape of studied particles is either spherical or assumed to be spherical, primarily for employment of existing analytical contact formulations that were developed specifically for spherical particles due to obvious geometrical simplicity. The elastic mechanical properties, typically the Young’s Modulus and Poisson’s ratio, are commonly characterized using Hertz contact law (Hertz, 1881), while the inelastic properties, typically the hardness or representative strength of the material, are characterized using various contact models developed for describing spherical indentation (Tabor, 1951; Johnson, 1985; Biwa & Storåkers, 1995) and contact of inelastic solids of revolution (Storåkers et al., 1997). It is important to note that these contact models are applicable only in the range of small deformations due to various limiting assumptions, such as independent contacts and simplified contact surface profile curvatures (Gonzalez & Cuitiño, 2012; Agarwal & Gonzalez, 2018).

In this work, we restrict attention to micro-crystalline cellulose (Avicel PH-200) particles, with average diameters ranging between $50$ and $300~{}\mu$m, under large diametrical compression. Micro-crystalline cellulose is one of the widely used binding materials in pharmaceutical powder blends for die-compaction of solid tablets (Jivraj et al., 2000). MCC particles exhibit significant plastic deformations under diametrical compression and elastic relaxation during unloading, which makes them relevant for this study (Mashadi & Newton, 1987; Mohammed et al., 2006). We address the first set of challenges described above by using a Shimadzu MCT-510 micro-compression tester equipped with a flat punch tip of $500~{}\mu$m in diameter, a top microscope, and a side camera capable of recording the deformation process. The micro-compression tester is capable of applying loading-unloading cycles under load control mode within a wide loading range of $9.8$ to $4803$ mN and a displacement range of [math] to $100~{}\mu$m, at a minimum increment of $0.001~{}\mu$m. We then characterize particle size and shape using the top microscope and the side camera, and record accurate force-displacement curves over a range of large deformations. The characterization of these irregular particles under lateral confinement clearly increases the experimental complexity and, even though it provides complementary information for building a semi-empirical mechanistic model, it will be beyond the scope of this work.

Despite the challenges described above, it is expected, as is the case for a wide range of physical responses, that a special structure exists in the response of micro-sized particles under large diametrical compression, which can be discovered and exploited to create a semi-empirical mechanistic model. Specifically, we assume that the force-displacement response exhibits a special structure known as an active subspace (see Tripathy et al. (2016) and references therein), i.e., a manifold of the stochastic space of inputs (such as particle size, shape, surface roughness, internal porosity, material properties, loading conditions and confinement) characterized by maximal variation of the multivariate response function (i.e., the contact force). We identify this low dimensional manifold based on mechanistic understanding of the problem, project onto it the high-dimensional space of inputs, and link the projection to the measured output contact force to build a semi-empirical mechanistic contact law. It is worth noting that the high-dimensional space of inputs is typically not amenable to a full experimental characterization and thus some components of the active subspace need to be estimated as part of the process of building the response function. Specifically, we estimate three plastic and one elastic material properties, and one geometric parameter associated with the loading condition through the shape factor function. Finally, we account for model inaccuracies and/or imperfections by endowing the three plastic material properties with a probability distribution, namely a log-normal distribution.

It bears emphasis that a semi-empirical mechanistic contact law for elasto-plastic particles is also relevant to three-dimensional particle mechanics calculations (Gonzalez & Cuitiño, 2012, 2016; Yohannes et al., 2016, 2017; Gonzalez et al., 2018). The particle mechanics approach enabled the prediction of microstructure evolution during the three most important steps of powder die-compaction (namely compaction, unloading, and ejection) using generalized loading-unloading contact laws for elasto-plastic spheres with bonding strength (Gonzalez, 2019). Moreover, these detailed calculations enabled the development of microstructure-mediated process-structure-property-performance interrelationships for QbD and QbC product development and process control, which depend on a small number of parameters with well-defined physical meanings. Therefore, the work presented in this paper, in combination with particle mechanics calculations, can contribute towards establishing the relationship between particle-level material properties and tablet performance.

The paper is organized as follows. Contact mechanics formulations for elastic ellipsoidal particles are reviewed in Section 2, followed by the introduction of shape factor and master contact law concepts, which are subsequently generalized to plastic irregular particles. The experimental analysis of MCC particles under diametrical compression is discussed in Section 3. Plastic shape factor for MCC particles is proposed in Section 4, followed by the description of an optimization procedure used to determine the geometric parameter associated with the loading condition. in Section 5, a master contact law for micro-crystalline cellulose particles is proposed, and log-normal distributions for the three plastic material properties are estimated from experimental data. Section 6 extends the work to loading-unloading contact laws, and Section 7 shows the comparison between predictions of the calibrated contact law and the experimental values. Finally, a summary is presented in Section 8.

2 Contact mechanics formulations for ellipsoidal particles

The contact mechanics of ellipsoidal elastic particles, characterized by three diameters $D_{1}$, $D_{2}$ and $D_{3}$ (Fig. 1), indicates that the contact force $F$ applied on a particle by two rigid plates in the direction of $D_{3}$ is given by (Zheng et al., 2013; Johnson, 1985)

[TABLE]

where $\gamma$ is the relative displacement of the plates, $E$ and $\nu$ are the Young’s modulus and the Poisson’s ratio respectivley, $R_{e}=D_{1}D_{2}/2D_{3}$ is the effective radius of the ellipsoidal particle, and $\eta_{e}$ is a shape factor dependent on $D_{1}$ and $D_{2}$ , given by

[TABLE]

By defining effective stress and strain as

[TABLE]

[TABLE]

the following stress-strain relationship $\sigma(\eta_{e}\epsilon)$ is recovered

[TABLE]

which depends on material properties ($E$ and $\nu$) and particle geometry ($D_{1}$ and $D_{2}$). It is worth noting that for spherical particles, the effective radius $R_{e}$ is equal to the radius of the sphere, the elastic shape factor $\eta_{e}$ simplifies to 1, and Hertz theory is recovered. Another noteworthy consideration from the above equation is that $\eta_{e}$ is independent of the particle height $D_{3}$. This is an artefact of the local character of the formulation, which assumes that the contacts evolve locally and independently without any interactions with other contacts on the same particle. Consequently, evolution of the contact is dependent only on the in-plane dimensions, i.e. $D_{1}$ and $D_{2}$. In the context of confined granular systems, this assumption has been found to hold only during the inital stages of compression (Mesarovic & Fleck, 2000). For linear-elastic spherical particles, the nonlocal contact formulation by Gonzalez & Cuitiño (2012), and its later extension by Agarwal & Gonzalez (2018) to correct particle contact areas, is a significant contribution towards relaxation of this limiting assumption. The formulation invokes the principle of superposition to describe the normal contact deformation and evolution of inter-particle contact area as a sum of local (i.e., Hertzian) deformation and nonlocal deformations generated by other contact forces acting on the same particle. For elasto-plastic materials, development of a closed-form analytical formulation for general contact configurations that can describe contact interactions at large deformations (Tsigginos et al., 2015) still remains an open problem. Initial progress in this regard has been made by Frenning (2013), who proposed a truncated sphere model applicable for small to moderate deformations, that utilizes the plastic incompressibility assumption to relate the average pressure in the particle due to elastic volumetric strain to the mean pressure generated at the particle contacts.

In this work, we interpret the function $\eta_{e}$ as the elastic shape factor that depends solely on particle geometry and loading condition, and the function $\sigma(\cdot)$ as the master contact behavior of the particle that depends solely on the particle deformation behavior (i.e., on its elastic behavior). Furthermore, we generalize this interpretation to plastic ellipsoidal particles as follows

[TABLE]

where the plastic master behavior depends on the deformation mechanisms of plastic particles and the shape factor is

[TABLE]

It is worth noting that the plastic shape factor depends on $D_{3}$ in addition to $D_{1}$ and $D_{2}$, which emphasizes its applicability to large deformations. In order to demonstrate this generalization, we develop an experimental and numerical procedure to identify (i) the plastic master contact behavior $\sigma(\cdot)$ for micro-crystalline cellulose (Avicel PH-200) under large plastic deformations, and (ii) the shape factor $\eta_{p}$ for irregular particles that can be approximated by ellipsoids (Fig. 2).

It is important to understand the characteristics of micro-crystalline cellulose particles and their dependency on environmental and process variables (see, e.g., Sun (2008); Inghelbrecht & Remon (1998); Westermarck et al. (1999)). These properties gather further importance in implementing QbD and QbC strategies, which require the understanding of raw material attributes for robust in-line monitoring and control of manufacturing process (Thoorens et al., 2014). Avicel has excellent compactibility, high dilution potential and superior disintegration properties, making it a prime filler material for the challenging direct compression tabletting for a wide range of active pharmaceutical ingredients (APIs) (Jivraj et al., 2000; Çelik et al., 1996; Doelker, 1993; Reier & Shangraw, 1966). Different grades of micro-crystalline cellulose exist in the pharmaceutical industry, increasing the flexibility of its usage to not only direct compression but also to dry granulation and wet granulation processes (Inghelbrecht & Remon, 1998; Westermarck et al., 1999). Some of the grades have different functional usages and not just physical variations such as to counter the weakening effect of lubricants on tablet strength (van Veen et al., 2005).

3 Diametrical compression of micro-crystalline cellulose particles

Micro-crystalline cellulose (Avicel PH-200) particles were tester under diametrical compression using a Shimadzu MCT-510 micro-compression tester (see Fig. 3) equipped with a flat punch tip of $500~{}\mu$m in diameter. Load-unload tests were performed on 20 particles at a constant stress rate of $2$ mN/s, with maximum applied stress in the range of $7-16$ MPa. The characteristic diameters $D_{1}$ and $D_{2}$ were measured from the top view using the in-built microscope, while $D_{3}$ was calculated from a side view obtained using the side camera. The obtained dimensions for all particles are listed in Table 1, while morphology of the tested particles is illustrated in Fig. 4, showing the relationship between their average diameter, i.e., $(D_{1}+D_{2}+D_{3})/3$, and their effective diameter, i.e., $D_{1}D_{2}/D_{3}$. It is interesting to note that this relationship has a linear trend, and that the ratio of the two values is about 1.5 for Avicel PH-200.

As an illustrative example, Fig. 5 shows an Avicel PH-200 particle at three increasingly higher levels of diametrical compression. It is evident from the figure that the particle undergoes large deformations without apparent brittle failure, but rather permanent plastic deformation, as confirmed from the unloading elastic curves in Fig. 6. It is also evident from the figures that surface irregularities are eliminated at low levels of diametrical compression. However, Fig. 6 illustrates that the strain required to eliminate these irregularities, which is measured from the first contact between the particle surface and the flat punch tip, is different for each particle and, therefore, is not predictable by a deterministic contact mechanics formulation such as equations (6)-(7). Therefore, in this study, we adopted a threshold stress of $0.2$ MPa to separate stochastic from deterministic deformation behaviors (i.e., to separate the initial deformation behavior dominated by uncharacterized surface irregularities from the subsequent behavior dominated by plastic deformations of an equivalent ellipsoidal particle). Specifically, we assumed a linear stress-strain response for the initial deformation behavior, with strain at $0.2$ MPa equal to the averaged value obtained from the 20 diametrical compression tests. The stress-strain response of each tested particle was then adjusted by offsetting the strain such that all deformation curves had the same averaged strain value at $0.2$ MPa. Fig. 6 shows raw and adjusted curves for three different particles. It is worth noting that the stress of $0.2$ MPa corresponds to 1-3$\%$ of the maximum stress applied during the tests. The adjusted diametrical compression curves were subsequently used to determine the shape factor $\eta_{p}(D_{1}/D_{3},D_{2}/D_{3})$ and the master contact law $\sigma(\cdot)$ for micro-crystalline cellulose particles.

4 Plastic shape factor for micro-crystalline cellulose particles

We proposed in Section 3 that the plastic shape factor for ellipsoidal plastic particles $\eta_{p}$ must depend solely on $D_{1}/D_{3}$ and $D_{2}/D_{3}$ and that it must simplify to 1 for spherical particles (Eqn. (7)). In addition, using Eqn. (6), we propose that the stress-strain relationship $\sigma_{i}(\epsilon)$ of an ellipsoidal particle $i$ is related to the master contact law $\sigma(\epsilon)$ by

[TABLE]

Equivalently, the strain-stress relationship $\epsilon(\sigma)$ of the master contact law and that of an ellipsoidal particle $i$, $\epsilon_{i}(\sigma)$, are related by

[TABLE]

Even though Eqns. (8) and (9) are equivalent, we shall use the experimentally characterized curves $\epsilon_{i}(\sigma)$ to determine the master contact law $\epsilon(\sigma)$ and the plastic shape factor $\eta_{p}$. This choice is motivated by the observation that micro-crystalline cellulose particles exhibit an apparent strain-hardening at distinctly different strain values (see Fig. 6). Hence, the domain of the function $\epsilon_{i}(\sigma)$, that is $[0,\sigma^{\mathrm{m}}_{i}]$, is quite similar for all the characterized particles—with $\sigma^{\mathrm{m}}_{i}$ being the maximum applied stress during diametrical compression of particle $i$. Moreover, Avicel PH-200 particles are inhomogeneous and not perfectly ellipsoidal; thus, a unique function (9) that holds true for all tested particles does not exist. Therefore, we propose that the master contact law is given by a normal distribution

[TABLE]

with expectation given by the sample mean $\bar{\epsilon}$ and variance given by square of the sample standard deviation $s^{2}_{\epsilon}$, i.e.

[TABLE]

for all $\sigma$ such that $\#\mathcal{S}_{\sigma}>N/5$. In the above equations $N$ is the total number of tested particles (i.e., $20$ in this study), $\mathcal{S}_{\sigma}$ is the set of experiments for which a stress of $\sigma$ has been measured or reached before unloading, and $\#\mathcal{S}_{\sigma}$ is the cardinality or the number of elements of $\mathcal{S}_{\sigma}$. In order to achieve statistical significance, a sufficient number of particles needs to be tested. For the purpose of demonstrating the proposed method, we adopt such number to be $N/5$. Next, we propose the following shape factor for micro-crystalline cellulose particles

[TABLE]

where $a$ is a geometric parameter associated with the loading condition, to be determined from the experimental strain-stress curves $\epsilon_{i}(\sigma)$. Specifically, the parameter $a$ is such that the relative standard deviation of $\epsilon$ is minimized over a range of stresses $[0,\sigma^{\mathrm{m}}_{\sigma,i}]$, i.e., $a$ minimizes

[TABLE]

with $\sigma^{\mathrm{m}}_{\sigma,i}=\min\{\sigma^{\mathrm{m}}_{i},\sigma^{\mathrm{m}}\}$, $\sigma^{\mathrm{m}}$ being the maximum stress value in the master contact law $\bar{\epsilon}(\sigma)$, i.e., in Eqn. (11). In order to simplify the optimization process, we first identify a particle with $D_{3}^{2}/D_{1}D_{2}$ close to 1 which is diametrically compressed at a high maximum stress and label it as particle $N$, and then rewrite Eqn. (14) using

[TABLE]

which leads to the following approximate minimization problem

[TABLE]

with $\sigma^{\mathrm{m}}_{\sigma,i}=\min\{\sigma^{\mathrm{m}}_{i},\sigma_{N}^{\mathrm{m}}\}$. It is worth noting that the experimental curves are linearly interpolated and numerically integrated for the purpose of solving Eqn. (15).

The proposed optimization process results in a shape factor, Eqn. (13), for the tested Avicel PH-200 particles with an optimal coefficient $a=1.43$. The master contact behavior for these particles is determined in the next section.

5 Master contact law for micro-crystalline cellulose particles

We proposed in the previous section that the master contact law is given by a normal distribution $\epsilon(\sigma)\sim\mathcal{N}(\bar{\epsilon}(\sigma),s^{2}_{\epsilon}(\sigma))$, where the sample mean $\bar{\epsilon}(\sigma)$ and the sample standard deviation squared $s_{\epsilon}^{2}(\sigma)$ are given by Eqns. (11)-(12) and the shape factor $\eta_{p}$ is given by Eqn. (13) with $a=1.43$. Therefore, the master contact law is readily available from the experimental curves and is shown in Fig. 7.

The stress-strain relationship $\sigma(\epsilon)$ of the master contact law for micro-crystalline cellulose particles is given by the inverse of $\epsilon(\sigma)$. We propose to approximate $\sigma(\epsilon)$ by

[TABLE]

with $\log(C_{i})\sim\mathcal{N}(\mu_{C_{i}},\sigma^{2}_{C_{i}})$. Therefore, the expected value of $\sigma(\epsilon)$ is given by

[TABLE]

The approximated stress-strain relationship is readily available by calibration and it is shown in Fig. 8 with coefficients $C_{i}$ given in Table 2. The expected value is thus given by

[TABLE]

The stress-strain relationship given by Eqn. (16) is in spirit of the curvature-corrected nonlocal contact formulation for elastic spherical particles (Agarwal & Gonzalez, 2018). The first term in the relationship corresponds to the initial linear deformation behavior (represented by a dashed-dotted curve in Fig. 8), which has been previously described with contact models developed for spherical indentation (Tabor, 1951; Johnson, 1985; Biwa & Storåkers, 1995) and contact of inelastic solids of revolution (Storåkers et al., 1997) in the context of small deformations. Of particular interest is the similarity solution by Biwa & Storåkers (1995), which describes the stress-strain relationship at the contact of rigid-plastic power law hardening solids. For spherical particles under diametrical compression, the relationship is given by

[TABLE]

where $m$ is the power law hardening exponent, $\kappa$ is a representative strength given by (Mesarovic & Fleck, 1999, 2000)

[TABLE]

where $\sigma_{y}$ is the uniaxial yield stress and $E$ is the Young’s modulus, and $c^{2}_{max}$ is a calibrated material parameter, well approximated by the following relationship (Storåkers & Larsson, 1994)

[TABLE]

Assuming perfectly plastic material behavior (i.e., $m\to\infty$), the stress-strain relationship reduces to

[TABLE]

and, by correlating the calibrated value of plastic coefficient $C_{1}$ with the above equation, the uniaxial yield stress is estimated as $\sigma_{y}=0.4085\,\mathrm{MPa}$. Other terms in the master contact law stress-strain relationship describe moderate-to-large deformation behavior of the particle. Following the initial linear response, a softening or reduction in stiffness is observed as a dip in the slope of the stress-strain curve. This regime of deformations has been termed as the contact interaction regime (Tsigginos et al., 2015), which starts with coalescence of two or more plastically deforming zones surrounding the contacts (Frenning, 2013; Tsigginos et al., 2015). With the onset of this regime, contacts can no longer be assumed independent. The second term of the master stress-strain relationship describes the contact response in this regime. The third term describes the large deformation strain-hardening response, referred to as geometric hardening (Sundstrom & Fischmeister, 1973) or the low compressibility regime (Tsigginos et al., 2015), during which the average contact pressure unboundedly rises due to increasing contact interactions.

6 Loading-unloading contact law for micro-crystalline cellulose particles

Unloading contact laws for elasto-plastic spheres with bonding strength, or adhesion, have been developed (Mesarovic & Fleck, 1999, 2000; Mesarovic & Johnson, 2000), assuming elastic perfectly-plastic behavior and using a rigid punch decomposition (Hill, 1990). Olsson & Larsson (2013) have extended these laws to elasto-plastic spheres that exhibit power-law plastic hardening behavior, and have verified their validity with detailed finite element simulations. This formulation assumes elastic behavior, approximated by Hooke’s law, and Irwin’s fracture mechanics to describe elastic recovery of the deformed spheres and the breakage of solid bridges. Gonzalez (2019) developed generalized loading-unloading contact laws for elasto-plastic spheres with bonding strength, which are continuous at the onset of unloading by means of a regularization term, in the spirit of a cohesive zone model. These contact laws are updated incrementally to account for strain path dependency and have been shown to be numerically robust, efficient, and mechanistically sound in three-dimensional particle mechanics static calculations. Here we endow the master contact law for micro-crystalline cellulose particles, proposed in Section 5, with elastic relaxation during unloading. Specifically, we follow Gonzalez (2019) and modify Eqn. (16) accordingly, i.e.

[TABLE]

where the internal variables are updated, i.e., $\{\epsilon_{\mbox{\tiny P}},\sigma_{\mbox{\tiny P}}\}\leftarrow\{\epsilon,\sigma(\epsilon)\}$, during plastic loading (i.e., when $\epsilon\geq\epsilon_{\mbox{\tiny P}}$). In the above equation, the function $\phi(\epsilon,\epsilon_{\mbox{\tiny P}})$ is given by

[TABLE]

for an elastic, perfectly plastic particle (i.e., for $m\rightarrow\infty$) with effective elastic modulus equal to $\bar{E}=E/(1-\nu^{2})=16$ GPa. The assumption of perfectly plastic material is in agreement with the observation made in Section 5 and the value of the effective elastic modulus is the result of a simple calibration.

7 Results and discussion

We have proposed that the experimentally characterized micro-crystalline cellulose particles deform under diametrical compression following a master contact law $\sigma(\epsilon)$, given by Eqn. (23) with plastic coefficients $C_{i}$ in Table 2, effective elastic modulus $\bar{E}=16$ GPa, and a shape factor $\eta_{p}$, given by Eqn. (13) with $a=1.43$. Figs. 9, 10 and 11 exhibit a very good agreement between predictions of the calibrated contact law and the experimental values. It is evident from the figures that the apparent plastic strain-hardening occurs at distinctly different strain values and that the proposed contact law captures its dependency on particle dimensions—i.e., on $D_{1}$, $D_{2}$ and $D_{3}$ (see Fig. 4). It is also evident that the elastic relaxation during unloading is accurately predicted with low uncertainty.

8 Summary

We have proposed a semi-empirical mechanistic contact law for micro-crystalline cellulose (Avicel PH-200) particles. This loading-unloading contact law has been characterized experimentally using diametrical compression force-displacement curves, obtained with a Shimadzu MCT-510 micro-compression tester. The irregular MCC particles have been approximated by an ellipsoid and the lengths of their principal axes have been measured using an in-built microscope and a side camera. To generalize the contact mechanics of an elastic ellipsoidal particle to an elasto-plastic irregular particle approximated by an ellipsoid, we have introduced the concepts of a shape factor and a master contact law. It is worth noting that the force-displacement curves exhibit an apparent strain-hardening at distinctly different strain values. We have postulated that this deformation mechanism depends on particle dimensions and loading configuration,which is captured by the shape factor function. The proposed loading-unloading contact law is, therefore, a function of (i) three characteristic diameters (lengths of the principal axes of an approximated ellipsoid), (ii) a geometric parameter associated with the loading condition through the shape factor function, (iii) three plastic material properties, and (iv) one effective elastic material property. The three plastic material properties are log-normal distributions estimated from the loading experimental curves, while the effective elastic property is estimated from the unloading experimental curves. It bears emphasis that the proposed loading contact law is in spirit of the curvature-corrected nonlocal contact formulation for elastic spherical particles (Agarwal & Gonzalez, 2018), with the first term correponding to the stress-strain relationship described by small-strain contact models developed for spherical indentation (Tabor, 1951; Johnson, 1985; Biwa & Storåkers, 1995) and contact of inelastic solids of revolution (Storåkers et al., 1997). Similarly, the proposed unloading contact law follows from generalized loading-unloading contact laws for elasto-plastic spheres with bonding strength (Gonzalez, 2019). The study shows a very good agreement between predictions of the calibrated loading-unloading contact law and the experimental values. However, it is recommended that a larger number of experiments be performed to estimate the five parameters of the contact law, in order to reduce the uncertainty in the model predictions.

We close by pointing out that the proposed semi-empirical mechanistic contact law is relevant to three-dimensional particle mechanics calculations (Gonzalez & Cuitiño, 2012, 2016; Yohannes et al., 2016, 2017; Gonzalez et al., 2018; Gonzalez, 2019). Therefore, the work presented in this paper, in combination with these detailed calculations, can contribute to develop microstructure-mediated process-structure-property-performance interrelationships and, thus, to establish the relationship between particle-level material properties and tablet performance. Ultimately, these relationships are needed to assist QbD and QbC product development and process control (Yi et al., 2018).

Acknowledgements

The author gratefully acknowledges the support received from the United States National Science Foundation grant number CMMI-1538861 and from the United States Food and Drug Administration grant number DHHS-FDA U01FD005535. The views expressed by authors do not necessarily reflect the official policies of the Department of Health and Human Services; nor does any mention of trade names, commercial practices, or organization imply endorsement by the United States Government.

Nomenclature