Generalized loading-unloading contact laws for elasto-plastic spheres with bonding strength
GGeneralized loading-unloading contact laws for elasto-plasticspheres with bonding strength
Marcial Gonzalez ∗ School of Mechanical Engineering, Purdue University,West Lafayette, IN 47907, USAMarch 21, 2018
Abstract
We present generalized loading-unloading contact laws for elasto-plastic spheres with bond-ing strength. The proposed mechanistic contact laws are continuous at the onset of unloading bymeans of a regularization term, in the spirit of a cohesive zone model, that introduces a smalland controllable error in the conditions for interparticle breakage. This continuity property is insharp contrast with the behavior of standard mechanistic loading and unloading contact theories,which exhibit a discontinuity at the onset of unloading when particles form solid bridges duringplastic deformation. The formulation depends on five material properties, namely two elasticproperties (Young’s modulus and Poisson’s ratio), two plastic properties (a plastic stiffness anda power-law hardening exponent) and one fracture mechanics property (fracture toughness), andits predictions are in agreement with detailed finite-element simulations. The numerical robust-ness and efficiency of the proposed formulation are borne out by performing three-dimensionalparticle mechanics static calculations of microstructure evolution during the three most impor-tant steps of powder die-compaction, namely during compaction, unloading, and ejection. Thesesimulations reveal the evolution, up to relative densities close to one, of microstructural features,process variables and compact mechanical attributes which are quantitatively similar to those ex-perimentally observed and in remarkable agreement with the (semi-)empirical formulae reportedin the literature.
Many physical mechanisms are required to convert a powder bed confined inside a rigid die into acompressed solid compact by the sole application of a compaction force. Typically, the initial stageof this process is characterized by rearrangement of particles that leads to the formation of a closelypacked granular system. In the subsequent stage, the porosity or the packing volume cannot be furtherreduced by particle rearrangement and therefore particles undergo brittle fracture or plastic deforma-tion, or both [8, 4]. It is indeed these dissipative and irreversible processes, in which the volume ofthe powder bed is reduced, that ultimately give rise to compact formation inside the die. Specifically,fracture and permanent deformation generate particle-to-particle contact surface and thus the oppor-tunity for bond formation. The understanding of microstructure formation and evolution during this ∗ [email protected] a r X i v : . [ c ond - m a t . s o f t ] M a r . Gonzalez Generalized loading-unloading contact laws process is therefore of paramount importance to elucidate strength formation. Particle size, shape,and roughness affect the initial stage of compaction, but it is fragmentation and plastic deformationthat dominate the synthesis of highly dense compacts [19]. For polymeric solids, fracture and plas-tic behavior are dependent on the physical form of the material; that is amorphous polymers havea tendency to ductile, elasto-plastic deformation whereas crystalline polymers exhibit brittle failureat room temperature [40]. In addition, most materials exhibit a brittle-ductile transition temperaturethat is pressure dependent and strain rate dependent—high temperature and low strain rate promoteductile, plastic behavior, while low temperature and high strain rate promote brittle fracture [40]. Itis worth noting that the formation of particle-to-particle contact surface then clearly depends on bothmaterial properties and process variables, such as compaction speed and temperature.Bonding surface area can be regarded as the effective surface area that is involved in the inter-action between particles. According to Rumpf [62], the bonding mechanisms participating in thisinteraction can be classified into five different types: (i) formation of solid bridges (driven by pro-cesses such as sintering, melting, crystallization of amorphous solids, or chemical reactions); (ii)bonding between movable liquids (caused by capillary and surface tension forces in the presenceof some moisture); (iii) non-freely movable binder bridges (resulting from adhesive and cohesiveforces in binders such as those used in wet granulation); (iv) long-range attractive forces betweensolid particles (such as van der Waals and hydrogen bonding interactions); and (v) mechanical par-ticle interlocking. For particles with low aspect ratio and low roughness, mechanical interlockingcan be neglected. Among the remaining mechanisms, there is general agreement that solid bridgeformation and attractive interfacial forces are the major contributions to strength formation. Figure 1illustrates these different mechanisms and that attractive interfacial interactions are dominant undersmall deformations whereas solid bridge formation occurs under large deformations. −−−−− − −+− −+ +++++ . + .+ . van der Waals andhydrogen bonding interactionsadhesive and cohesiveforces solid bridge sintering, melting,crystallization of amorphous solids,chemical reactions Figure 1: Different bonding mechanisms participating in the interaction between particles. Attractiveinterfacial interactions are dominant under small deformations whereas solid bridge formation occursunder large deformations.During powder compaction, regions of high plastic deformation are formed around particle con-tact interfaces. Plastic deformation dissipates energy as heat and thus locally increases the tempera-ture. This change in temperature may, in turn, promote molecular movement and, consequently, theformation of a new solid region that bridges the particles in contact. It is then the formation of aninterconnected network of solid bridges that enables strength formation in the solid compact. Thestrength of these solid bridges will vary from material to material depending on the forces that holdthe (poly)crystals and amorphous solids together [63, 16, 49, 3]. It is worth noting that this process isirreversible, and thus it is not possible to divide the compacted system into its original particles. The. Gonzalez Generalized loading-unloading contact laws system can only be separated by fracture of either solid bridges or particles, whichever is weaker. Thisirreversibility is one of the main differences between solid bridges and the attractive interfacial inter-actions. This observation also suggests to characterize the strength of the solid bridges with fracturemechanics properties and to characterize the strength of the solid compact by its tensile strength.In this work, we will restrict attention to powder blends used by the pharmaceutical industry tofabricate solid tables, the most popular dosage form in use today. Therefore, it is worth noting that thecontact surface created during compaction of these powders also allows for the formation of attractiveforces such as van der Waals and hydrogen bonding interactions [35, 33, 13, 14, 34]. These long-rangeforces are of lower energy than covalent bonding forces, and they are present in many excipients (suchas sugars, celluloses, and starches) and active pharmaceutical ingredients. It is also important to notethat some polymers can experience a change in physical form with relative humidity (e.g., lactose)and that some active ingredients can experience strain- and temperature-driven solid-state transitionsand amorphization. These additional physical mechanisms clearly increase the complexity of theanalysis and, even though necessary for specific powder blends [67], their study will be beyond thescope of this work. We will specifically restrict attention to the formation of solid bridges, as it is aphysical mechanism that dominates the synthesis of many, but not all, pharmaceutical excipients. Wewill simplify the powder morphology to spherical particles, and we will consider that these particlesare amenable to elasto-plastic deformation without brittle failure.A quantitative elucidation of strength formation requires not only the identification of the defor-mation and bonding mechanisms of interest but also of the bonding surface involved in the process.Unfortunately, it is not possible to experimentally measure the actual interfacial area that is availableduring tableting [39, 52]. However, one can assume that an upper bound for the bonding surfaceinvolved in the formation of solid bridges is the particle-to-particle contact area created during com-paction. It bears emphasis that the presence of lubricants in the formulation can diminish the bondingsurface, and thus the tablet strength [39, 12, 61]. The most common lubricant used in pharmaceuticaltablets is magnesium stearate, typically prepared as a small particle size ingredient. Before tableting,the lubricant is mixed with the particles in the formulation, partially coating their surface and thusaltering the tribochemical properties of the particle-to-particle contact area. In the context of thiswork, we will assume that the lubricant alters the fracture mechanics properties of the solid bridges.In this paper, we report three-dimensional particle mechanics static calculations that enable us topredict microstructure evolution during compaction, unloading, and ejection—that is during the threemost important steps of die-compaction of solid tablets (see Figure 2)—of elasto-plastic sphericalparticles capable of forming solid bridges. To this end, we develop and employ generalized loading-unloading contact laws for elasto-plastic spheres with bonding strength. The proposed loading-unloading contact laws are continuous at the onset of unloading by means of a regularization termthat introduces a small, controllable error in the solid bridge breakage force and the critical contactsurface. The contact laws are explicit in terms of the relative position between the particles, and theirstrain path dependency is accounted for incrementally. The resulting formulation is then numericallyrobust and efficient, and mechanistically sound.The paper is organized as follows. The generalized loading-unloading contact laws for elasto-plastic spheres with bonding strength are presented and validated in Section 3, after reviewing thestate of the art in Section 2. The particle mechanics approach used to generate a sequence of staticequilibrium configurations of granular systems at high levels of confinement is presented in Section 4.The evolution of microstructural statistical features and of macroscopic effective properties duringcompaction, unloading and ejection is investigated in Section 5. Specifically, we study the evolutionof the mechanical coordination number (number of non-zero contact forces between a particle and its. Gonzalez Generalized loading-unloading contact laws A" B" C" D" compac+on" unloading" ejec+on"
Figure 2: Transformation of the powder bed during the stages of die compaction: (A) die filling, (B)compaction, (C) unloading, and (D) tablet ejection.neighbors), punch and die-wall pressures, in-die elastic recovery, residual radial pressure and ejectionpressure, the network of contact forces and granular fabric anisotropy, bonding surface area, Young’smodulus and Poisson’s ratio of the compacted solid, and a microstructure-mediated process-structure-property-performance interrelationship. Finally, concluding remarks are collected in Section 6.
Loading contact laws for elasto-plastic spheres have been developed by Stor˚akers and co-workers[72, 73] using a rigid plastic flow formulation [32] and assuming a power-law plastic hardeningbehavior, i.e., σ = κ(cid:15) /m where κ is the plastic stiffness and m is the plastic law exponent. Thesecontact laws have been generalized to dissimilar particles [45, 46]. Specifically, for particles withthe same hardening exponent or when one particle is assumed to be rigid, an analytical self-similarsolution is derived by assuming a contact radius sufficiently small compared to the particle size andby neglecting elastic behavior. For particles with different hardening exponents, Skrinjar and Larsson[68, 69, 70] derived and verified an approximate formulae based on the self-similar solution proposedby Stor˚akers. These loading contact laws for elasto-plastic spheres are successful in simulating thedeformation of soft metals, such as bronze and aluminum [56], and pharmaceutical excipients, such asmicrocrystalline cellulose and lactose monohydrate [82, 83], and of harder materials, such as ceramicsand cemented carbides, when both elastic and plastic deformations are properly accounted for duringthe loading phase [57].Unloading contact laws for elasto-plastic spheres with bonding strength, or adhesion, have beendeveloped by Mesarovic and Johnson [47] assuming elastic perfectly-plastic behavior and usinga rigid punch decomposition [31]. Olsson and Larsson [57] have extended these laws to elasto-plastic spheres that exhibit power-law plastic hardening behavior, and have verified their validitywith detailed finite element simulations. This formulation assumes elastic behavior, approximated byHooke’s law, and Irwin’s fracture mechanics to describe the elastic recovery of the deformed spheresand the breakage of the solid bridge.We present next these loading and unloading contact laws for elasto-plastic spheres with bondingstrength. Specifically, we consider two elasto-plastic spherical particles of radii R and R , Young’smoduli E and E , Poisson’s ratios ν and ν , plastic stiffnesses κ and κ , and plastic law exponent m , that deform plastically under loading and relax elastically under unloading. For particles locatedat x and x , the relative position between them γ (see Figure 1) is give by γ = R + R − (cid:107) x − x (cid:107) . Gonzalez Generalized loading-unloading contact laws and the contact radius a is given by a = c ¯ Rγ =: a P plastic loading (cid:34) a P − (cid:18) E ( γ P − γ )3 n P a /m P (cid:19) (cid:35) + elastic (un)loading (1)where the effective radius ¯ R , the effective elastic stiffness ¯ E and the plastic law coefficient n P aregiven by ¯ R = (cid:18) R + 1 R (cid:19) − ¯ E = (cid:18) − ν E + 1 − ν E (cid:19) − n P = πk ¯ R − /m (cid:18) κ m + 1 κ m (cid:19) − /m with k = 3 × − /m , c = 1 .
43 e − . /m [71, 22, 36], and [ · ] + = max {· , } . The permanent plasticdeformation is characterized by the plastic relative position γ P and the plastic contact radius a P , whichare related by a P = 2 c ¯ Rγ P . In addition, the elasto-plastic spherical particles are capable of forming asolid bridge characterized by its fracture toughness K Ic . Therefore, the plastic and elastic (un)loadingforce is defined by P = n P a /m plastic loading n P π a /m P (cid:34) arcsin (cid:16) aa P (cid:17) − aa P (cid:114) − (cid:16) aa P (cid:17) (cid:35) − K Ic π / a / elastic (un)loading (2)This force acts in direction ( x − x ) / (cid:107) x − x (cid:107) on particle 2 and in direction ( x − x ) / (cid:107) x − x (cid:107) on particle 1.This formulae is predictive at small deformations for particles that do not form solid bridges dur-ing plastic deformation, as it is depicted in Figures 3(a) and 3(b) when compared with detailed finiteelement simulations of an elasto-plastic continuum solid [21]. However, the formulation exhibits adiscontinuity at the onset of unloading when particles form solid bridges during plastic deformationthat it is not present in detailed finite element simulations as shown in Figure 3(c). Specifically, thereis a discontinuity at a = a P , that is P ( a + P ) = (cid:40) n P a /m P = P ( a − P ) if solid bridge is not formed n P a /m P − K Ic π / a / P (cid:54) = P ( a − P ) if solid bridge is formedwhere P ( a − P ) and P ( a + P ) correspond to the contact force right before and after unloading, respec-tively.The loading-unloading contact laws proposed in this work, and described next in Section 3, arecontinuous at the onset of unloading by means of a regularization term that introduces a small, con-trollable error in the solid bridge breakage force and the critical contact surface. Finally, at moderateto large deformations, detailed finite element simulations show a dependency of the response on theloading configuration and confinement of the particles (see Figure 3(d) and [23, 24, 37, 77]). Thisbehavior is not captured by the above local contact formulation and it calls for the development. Gonzalez Generalized loading-unloading contact laws Deformation [%] - ( γ /D) C on t a c t R ad i u s - ( a / R ) Deformation [%] - ( γ /D) P r e ss u r e [ M P a ] - ( P / π R ) (a) (b) Deformation [%] - ( γ /D) P r e ss u r e - ( P / R κ ) Deformation [%] - ( γ /D) P r e ss u r e - ( P / R κ ) (c) (d) Figure 3: Loading-unloading contact laws (solid lines) and finite element results (symbols) for elasto-plastic spheres. (a) Evolution of contact radius a , and (b) of contact force P , under diametricalcompression, small deformations and no formation of solid bridges for ( (cid:35) ) a spherical particle of R = 2 mm and mechanical properties E = 80 GPa, ν = 0 . , κ = 1 . GPa, m = 2 . and for ( (cid:51) ) aspherical particle of R = 10 mm and mechanical properties E = 200 GPa, ν = 0 . , κ = 1 . GPa, m = 2 . [21]. (c) Evolution of contact force P under diametrical compression and small defor-mations for a spherical particle of R = 1 . µ m, mechanical properties E = 233 GPa, ν = 0 . , κ = 12 . GPa, m = 1 . , that forms a solid bridge with fracture properties K Ic = 0 . MPa m / [18]. (d) Evolution of contact force P under large deformations and different loading configurations,namely diametrical compression ( (cid:77) ), diametrical compression and lateral confinement ( (cid:35) ) and triax-ial compression ( (cid:79) ), for a spherical particle of R = 5 mm and mechanical properties κ = 15 . MPa, m = 2 . [29].. Gonzalez Generalized loading-unloading contact laws of plastic nonlocal contact formulations (see [25, 26] for an elastic nonlocal contact formulation).The systematic development of nonlocal contact formulations for elasto-plastic spheres with bondingstrength is a worthwhile direction of future research and, though beyond the scope of this work, it iscurrently being pursued by the author. We adopt the formulae developed by Mesarovi and co-workers for elasto-plastic spheres with power-law plastic hardening behavior [45, 46, 47], presented in the previous section, and we propose aregularization of the contact force that does not modify the evolution of contact area a which is giveby a = c ¯ Rγ =: a P plastic loading (cid:34) a P − (cid:18) E ( γ P − γ )3 n P a /m P (cid:19) (cid:35) + elastic (un)loading (3)with γ P = a P / c ¯ R . It is worth noting that a solid bridge breaks during unloading at γ = a P / c ¯ R − n P a /m P / E (see Figure 4b). The regularized contact force P is defined by P = n P a /m plastic loading n P π a /m P (cid:34) arcsin (cid:16) aa P (cid:17) − aa P (cid:114) − (cid:16) aa P (cid:17) (cid:35) − K Ic π / a / ξ B ) [ a B − a ] + (1+ ξ B ) a B − a elastic (un)loading (4)where ξ B > is the regularization parameter and a B is the radius of the bonded area, or solid bridge,which evolves as follows a B := a P if mechano-chemical conditions are favorable, i.e., when ˙ a P > a B := 0 if solid bridge is broken, i.e., when a = 0˙ a B = 0 otherwise, i.e., the size of the bonded area does not change (5)Furthermore, the fracture toughness of the solid bridge is given by K Ic = √ G ¯ E , where the dissi-pated energy G includes interfacial fracture energy ω (i.e., surface and field forces at direct contact)and plastic or other type of dissipation G p , that is K Ic = (cid:115) ( ω + ω + G p )2 E E (1 − ν ) E + (1 − ν ) E A solid bridge between two particles of the same material then reduces to K Ic = (cid:112) GE/ (1 − ν ) .It bears emphasis that the contact force is continuous at a = a P , and it is equal to zero at a = 0 ,for any ξ B > and any value of a B . Therefore, this generalized loading-unloading contact laws forelasto-plastic spheres with bonding strength are continuous at the onset of unloading after formationof a solid bridge, and at the onset of plastic loading after breakage of a solid bridge or bonded surface(see Figure 4a). This is achieved by means of a regularization term that introduces a small, con-trollable error in the solid bridge breakage force and the critical contact surface—that we will studyin Section 3.2, after introducing a set of non-dimensional parameters in Section 3.1. The proposedcontact laws are validated in Section 3.3, followed by a reinterpretation of the regularization term asa cohesive zone model and by a sensitivity analysis in Sections 3.4 and 3.5, respectively.. Gonzalez Generalized loading-unloading contact laws Deformation [%] − ( γ /2R) F o r c e − ( P / κπ R ) plastic deformationand formation ofbonding strengthelastic(un)loadingbreakage ofbonded surfacebreakage ofbonded surface Deformation [%] − ( γ /2R) C on t a c t R ad i u s − ( a / R ) plastic deformationand formation ofbonding strength breakage ofbonded surfaceelastic(un)loading (a) (b) Figure 4: Generalized loading-unloading contact laws for elasto-plastic spheres with bondingstrength. (a) Evolution of contact force P under loading and two subsequent unloading-loading cy-cles that break the solid bridge. (b) Evolution of contact radius a . A regularization parameter ξ B equalto . is used. We recast the generalized loading-unloading contact laws for elasto-plastic spheres with bondingstrength presented above using the following non-dimensional parameters: (i) elastic recovery γ/γ P and a/a P , (ii) plastic deformation a P / ¯ R , (iii) bonded surface a B /a P , (iv) contact force P /πk ¯ κ ¯ R , (v)ratio of elastic to plastic stiffness ψ = 2 ¯ E/ πk ¯ κc , (vi) ratio of bonding to stored elastic energy χ = 4 πK Ic /n P a /m P . Therefore, the non-dimensional contact radius then simplifies to aa P = plastic loading, i.e., γ/γ P = 1 , a P := (cid:112) c ¯ Rγ P (cid:20) − (cid:16) ψ (cid:16) − γγ P (cid:17) (cid:0) a P ¯ R (cid:1) − /m (cid:17) (cid:21) / elastic (un)loading, i.e., γ/γ P < (6)and the non-dimensional contact force is Pπk ¯ κ ¯ R = (cid:0) a P ¯ R (cid:1) /m plastic loading π (cid:0) a P ¯ R (cid:1) /m (cid:34) arcsin (cid:16) aa P (cid:17) − aa P (cid:114) − (cid:16) aa P (cid:17) (cid:35) − χ (cid:0) a P ¯ R (cid:1) /m (cid:16) aa P (cid:17) / ξ B ) [ a B /a P − a/a P ] + (1+ ξ B ) a B /a P − a/a P elastic (un)loading (7)It is worth noting that bonding to stored elastic energy ratio χ is π − /π ≈ . times largerthan the bonding energy to elastic energy ratio in [47]. The proposed generalized loading-unloading contact laws introduce a controllable error in the solidbridge breakage force and the critical contact surface. In order to study these errors, we first derivethe critical force P ξ BC and contact radius a ξ BC for a give regularization parameter ξ B assuming, for. Gonzalez Generalized loading-unloading contact laws simplicity, a B /a P = 1 . The critical point occurs at the maximum tensile force, that is at ∂P∂γ = 0 = 16 a P ¯ E π ( aa P ) / (cid:18) aa P (cid:19) / − π / K Ic n P a / /m P × (1 − aa P ) + ξ B (1 − aa P )(1 − aa P + ξ B ) (1 + ξ B ) − (cid:115) − (cid:18) aa P (cid:19) (8)and thus, after performing a power series expansion in ξ B , the critical contact radius or critical pull-off force is given by the solution of (cid:18) a C a P (cid:19) = 9 π K Ic n P a /m P (cid:34) − (cid:18) a C a P (cid:19) − a C /a P )(3 − a C /a P )(3 a C /a P − − a C /a P ) ξ B + O ( ξ B ) (cid:35) (9)Specifically, the critical contact radius for the regularized contact law, a ξ BC , is give by the solution ofthe following equation (cid:32) a ξ BC a P (cid:33) / − π / K Ic n P a / /m P × (1 − a ξ BC /a P ) + ξ B (1 − a ξ BC /a P )(1 − a ξ BC /a P + ξ B ) (1 + ξ B ) − (cid:118)(cid:117)(cid:117)(cid:116) − (cid:32) a ξ BC a P (cid:33) = 0 which reduces to ( a ξ BC /a P ) = χ (3 π/ (1 − ( a ξ BC /a P ) ) for ξ B = 0 (cf. [46]). Finally, the corre-sponding critical force is given by equation (4), i.e., by P ξ BC = P ( a ξ BC ) , and, naturally, the correctlimiting behavior is retained, i.e., a C = a ξ B → C and P C = P ξ B → C . Figure 5 shows the error incurred inthe critical force and the critical contact radius for different regularization parameters and bonding toelastic energies ratios. Bonding to Elastic Energies Ratio -1 C r i t i c a l F o r c e ( P c ξ B / P c ) ξ B Bonding to Elastic Energies Ratio -1 C r i t i c a l R ad i u s ( a c ξ B / a c ) ξ B (a) (b) Figure 5: Error incurred in the critical force P ξ BC and the critical contact radius a ξ BC for differentregularization parameters ξ B = { , . , . , . , . , . } and bonding to elastic energiesratio χ .We next define the optimal regulation parameter ¯ ξ B as follows ¯ ξ B := min (cid:110) ξ B > s.t. (cid:15) = | − a ξ BC /a C | , ξ B > s.t. (cid:15) = | − P ξ BC /P C | (cid:111) where (cid:15) is the maximum relative error incurred in the critical contact radius and the critical force.Figure 6 shows that a regularization parameter of ξ B = 0 . ensures a moderate error in the prediction. Gonzalez Generalized loading-unloading contact laws of the critical point over a wide range of bonding to elastic energy ratio conditions. It is worth notingthat this controllable error is the cost we pay to achieve a numerically robust and efficient as well asa mechanistically sound formulation. Bonding to Elastic Energies Ratio -1 R egu l a r i z a t i on P a r a m e t e r ( ξ B ) -2 -1 ± ± Figure 6: Optimal regularization parameter ¯ ξ B for introducing less than or equal to and errorin a ξ BC and in P ξ BC simultaneously. Thin solid lines correspond to the optimal values of ξ B that keep a ξ BC within given bounds. Thin dashed lines correspond to the optimal values of ξ B that keep P ξ BC withingiven bounds. We compare next the proposed formulation with detailed finite element simulations with Lennard-Jones stresses at the interface [17] performed by Du et al. [18]. The finite element simulationscorrespond to an elasto-plastic spherical particle with radius R = 1 . µ m, Young’s modulus E =233 GPa, Poisson’s ratio ν = 0 . , yield stress equal to 1.94 GPa and linear hardening equal to 2 % of Young’s modulus. We approximate the bi-linear elasto-plastic law by an exponential plastic lawwith plastic stiffness κ = 12 . GPa and plastic law exponent m = 1 . . The bond interaction isrepresented by the Lennard-Jones potential between two parallel surfaces, which we approximateby an equivalent fracture toughness K Ic = 0 . MPa m / . The maximum separation is γ max =28 . nm. Figure 7 shows Du’s finite element calculations and the predictions of the generalizedloading-unloading contact law for elasto-plastic spheres with bonding strength using three differentvalues of ξ B and adopting a B /a P = 1 . It is evident in the figure that the proposed loading-unloadingcontact law is continuous at the onset of unloading by means of the regularization term (cf. Figure 3c). We show next that the proposed regularization, which provides continuity in the contact force atthe onset of unloading, is in the spirit of a cohesive zone model. Specifically, we show that therelationship between interfacial separation traction ˆ T and separation displacement ˆ γ follows a typicalcohesive law curve under strain control, for ξ B > . For simplicity, we restrict attention to a B /a P = 1 but, in general, the non-dimensional bonded surface is ≤ a B /a P ≤ .It is worth noting that the separation force P B is the term in the unloading contact force thatcorresponds to separation and breakage of the solid bridge. Similarly, the separation displacement. Gonzalez Generalized loading-unloading contact laws Deformation - ( γ / γ max ) F o r c e - ( P / P m a x ) ξ B Deformation - ( γ / γ max ) C on t a c t R ad i u s - ( a / a m a x ) (a) (b) Figure 7: Detailed finite element simulations (symbols) performed by Du et al. [18] and generalizedloading-unloading contact law predictions (solid line) for R = 1 . µ m, E = 233 GPa, ν = 0 . , κ = 12 . GPa, m = 1 . , and K Ic = 0 . MPa m / —the bonding energy to elastic energy ratiois χ = 0 . . Three different values of ξ B are use , . , and . . Dimensionless force vsdeformation (a) and dimensionless contact radius vs deformation (b) correspond to a single loading-unloading cycle with γ max = 28 . nm.is zero at the onset of unloading and equal to the critical separation ∆ γ c at solid bridge breakage.Therefore, a non-dimensional separation force ˆ P and a non-dimensional separation displacement ˆ γ are defined as ˆ P = P B K Ic π / a / P = ( a/a P ) / (1 + ξ B ) [1 − a/a P ] + (1 + ξ B ) − a/a P ˆ γ = ( γ P − γ )∆ γ c with ∆ γ c = 3 n P a /m P E and thus they are related by ˆ P = (1 + ξ B ) (cid:16) − (cid:112) − ˆ γ (cid:17) (cid:0) − ˆ γ (cid:1) / ξ B − (cid:112) − ˆ γ (10)We next define a separation traction T B as T B = P B K Ic n P π / a / /m P C ( ξ B ) where the correction factor C ( ξ B ) enforces G = (cid:82) ∆ γ c T B d( γ P − γ ) and is equal to C ( ξ B ) = (cid:90) (1 + ξ B ) (cid:16) − (cid:112) − ˆ γ (cid:17) (cid:0) − ˆ γ (cid:1) / ξ B − (cid:112) − ˆ γ dˆ γ The nondimensional separation traction thus simplifies to ˆ T = T B a /m P n P K Ic = 1 C ( ξ B ) (1 + ξ B ) (cid:16) − (cid:112) − ˆ γ (cid:17) (cid:0) − ˆ γ (cid:1) / ξ B − (cid:112) − ˆ γ (11). Gonzalez Generalized loading-unloading contact laws Figure 8 shows the non-dimensional separation force ˆ P and non-dimensional separation traction ˆ T asa function of the non-dimensional separation displacement ˆ γ for different regularization parameters ξ B . The similarity between a typical cohesive traction-separation curve under stain control and thetraction-separation curves depicted in the figure is evident, for ξ B > . Nondimensional Separation N ond i m en s i ona l F o r c e ξ B Nondimensional Separation N ond i m en s i ona l T r a c t i on ξ B (a) (b) Figure 8: Non-dimensional separation force ˆ P (a) and non-dimensional separation traction ˆ T (b) asa function of the non-dimensional separation displacement ˆ γ for different regularization parameters ξ B = { , . , . , . } . The proposed loading-unloading contact law depends on five material properties, namely two elasticproperties ( E and ν ), two plastic properties ( κ , m ) and one fracture mechanics property ( K Ic ). Inorder to gain insight into the role of these parameters and the coupled mechanisms involved, weperformed a sensitivity analysis whose results are presented in Figure 9. It is interesting to note thatthe bonding surface can be controlled by changing the Young’s modulus E without changing thepeak force, that is without changing the compaction force (see Fig. 9a). Bonding surface can also bemanipulated by changing the plastic stiffness κ (Fig. 9d), plastic law exponent m (Fig. 9e) and theporosity or relative density through γ (Fig. 9c) but, in contrast, this inevitably results in a change ofthe compaction force. This is valuable insight for product and process design integration, since, asmentioned above, compact strength is directly correlated to the bonding surface created during thecompaction. The particle mechanics approach for granular systems under high confinement, developed by Gon-zalez and Cuiti˜no [26], describes each individual particle in the powder bed, and the collective rear-rangement and deformation of the particles that result in a compacted specimen. This approach hasbeen used to predict the microstructure evolution during die-compaction of elastic spherical particlesup to relative densities close to one. By employing a nonlocal contact formulation that remains predic-tive at high levels of confinement [25], this study demonstrated that the coordination number dependson the level of compressibility of the particles and thus its scaling behavior is not independent of. Gonzalez Generalized loading-unloading contact laws Deformation [%] − ( (cid:97) / (cid:97) max ) F o r c e − ( P / P m a x ) E Deformation [%] − ( (cid:97) / (cid:97) max ) F o r c e − ( P / P m a x ) R Deformation [%] − ( (cid:97) / (cid:97) max ) F o r c e − ( P / P m a x ) γ Deformation [%] − ( (cid:97) / (cid:97) max ) C on t a c t R ad i u s − ( a / a m a x ) E Deformation [%] − ( (cid:97) / (cid:97) max ) C on t a c t R ad i u s − ( a / a m a x ) R Deformation [%] − ( (cid:97) / (cid:97) max ) C on t a c t R ad i u s − ( a / a m a x ) γ (a) E ± (b) R ± (c) γ ± Deformation [%] − ( (cid:97) / (cid:97) max ) F o r c e − ( P / P m a x ) κ Deformation [%] − ( (cid:97) / (cid:97) max ) F o r c e − ( P / P m a x ) m Deformation [%] − ( (cid:97) / (cid:97) max ) F o r c e − ( P / P m a x ) K Ic Deformation [%] − ( (cid:97) / (cid:97) max ) C on t a c t R ad i u s − ( a / a m a x ) κ Deformation [%] − ( (cid:97) / (cid:97) max ) C on t a c t R ad i u s − ( a / a m a x ) m Deformation [%] − ( (cid:97) / (cid:97) max ) C on t a c t R ad i u s − ( a / a m a x ) K Ic (d) κ ± (e) m ± . (f) K Ic ± Figure 9: Sensitivity analysis. Reference values P max and a max correspond to Fig. 7 and [18], thatis R = 1 . µ m, E = 233 GPa, ν = 0 . , κ = 12 . GPa, m = 1 . , K Ic = 0 . MPa m / , γ max = 28 . nm. The bonding energy to elastic energy ratio is χ = 0 . .. Gonzalez Generalized loading-unloading contact laws material properties as previously thought. The study also revealed that distributions of contact forcesbetween particles and between particles and walls, although similar at jamming onset, are very dif-ferent at full compaction—being particle-wall forces are in remarkable agreement with experimentalmeasurements reported in the literature.In this work, we extend the particle mechanics approach to the treatment of internal variables(i.e., a P and a B ) and their equations of evolution (i.e., Equations 3 and 5) under quasi-static evolution.Therefore, an equilibrium configuration is defined by the solution of a system of nonlinear equationsthat corresponds to static equilibrium of the granular system, that is sum of all elasto-plastic contactforces acting on each particle equals zero, that is (cid:88) j ∈N i P ( a ( R i + R j − (cid:107) x i − x j (cid:107) , a P ,ij ) , a P ,ij , a B ,ij ) x i − x j (cid:107) x i − x j (cid:107) = (12)where x i and N i are the position and all the neighbors of particle i , respectively, a P ,ij = a P ,ji and a B ,ij = a B ,ji by definition, P ( a, a P , a B ) is given by equation (4), and a ( γ, a P ) is given byequation (3). A sequential strategy is proposed to treat the nonlinear problem (see Algorithm 1).The equations of static equilibrium are solved for x = ( x T , ..., x TN ) T , for given internal variables a P = ( a P , , a P , , ..., a P ,N − ,N ) T and a B = ( a B , , a B , , ..., a B ,N − ,N ) T , by employing a trust-regionmethod [10, 11] that successfully overcomes the characteristic ill-posedness of the problem (e.g.,due to metastability [48]). The basic trust-region algorithm requires the solution of a minimizationproblem to determine the step between iterations, namely the trust-region step. This minimizationproblem is of the form min { ψ ( s ) : (cid:107) s (cid:107) ≤ ∆ } , where s = n +1 x − n x is the trust-region step, ∆ is atrust-region radius and ψ is a quadratic function that represents a local model of the objective functionabout n x , that is ψ ( s ) = (cid:107) n F + n Ks (cid:107) = (cid:104) n F , n F (cid:105) + (cid:104) n Ks , n F (cid:105) + (cid:104) n Ks , n Ks (cid:105) with n F and n K the global force vector and stiffness matrix at n x —the first term in the above equa-tion is not required in the minimization problem. It is worth noting that the trust-region step is notnecessarily in the direction of a quasi-Newton step and that the trust-region radius acts as a regu-larization term that controls the growth in the size of the least squares solution observed in mostill-posed [80]. Trading accuracy for performance, Byrd [6], among others, proposed to approximatethe minimization problem by restricting the problem to a two-dimensional subspace. Furthermore,the two-dimensional subspace can be determined by a preconditioned conjugate gradient process andthe trust-region radius can be adjusted over the iterative process (see, e.g., [50]). Here we adopt theimplementation available in MATLAB R2016a Optimization Toolbox.It is worth noting that an equilibrium configuration is not obtained by artificially damped orcooled-down dynamic processes but rather by iterative solvers that follow the energy landscapearound the solution of static equilibrium. The strain path dependency of the contact law is accountedfor incrementally by updating internal variables a P and a B at the new equilibrium configuration.Specifically, new particle-to-particle contact radii a new ij are computed at the converged equilibriumconfiguration and internal variables are updated as follows: a newP ,ij ← a new ij , a newB ,ij ← a new ij if a new ij > a P ,ij plastic loading and bond formation a newB ,ij ← if a new ij = 0 solid bridge is broken a newP ,ij ← a P ,ij , a newB ,ij ← a B ,ij otherwise elastic (un)loading (13). Gonzalez Generalized loading-unloading contact laws ParticleMechanicsApproach: sequential strategy for solving the equilibrium problemand for updating history-dependent internal variables.
Require:
Initial guess for particles’ coordinates x , current state of internal variables { a P , a B } ,TOL1, trust-region radius ∆ and the simplicial complex X Error ← TOL1 n ← while Error ≥ TOL1 do / ∗ Compute global force and global stiffness. ∗ / { n F , n K } ← ContactFormulation ( n x , a P , a B , X ) / ∗ Update coordinates with the trust-region step obtained by restricting the problem to a two-dimensional subspace [6]. ∗ / n s ← argmin s (cid:104) n Ks , n F (cid:105) + (cid:104) n Ks , n Ks (cid:105) subject to (cid:107) s (cid:107) ≤ ∆ n +1 x ← n x + n s / ∗ Compute a measure of convergence. ∗ / Error ← (cid:107) n s (cid:107) n ← n + 1 end while / ∗ Update internal variables based on converged solution. ∗ / { a newP , a newB } ← UpdateInternalVariables ( n x , a P , a B ) return { n x , a newP , a newB } We next report three-dimensional particle mechanics static calculations that enable us to study mi-crostructure evolution during die-compaction up to relative densities close to one, unloading and ejec-tion of elasto-plastic spherical particles with bonding strength. We employ the generalized loading-unloading contact laws presented in Section 3 which result in a numerically robust and efficient for-mulation. The contact laws are continuous at the onset of unloading by means of a regularization term,they are explicit in terms of the relative position between the particles, and their strain path depen-dency is accounted for incrementally. Here we adopt a regularization parameter equal to ξ B = 0 . .Three different relative densities are defined and used in the study, namely (i) the maximumrelative density of the compact inside the die ρ in-diemax which occurs at the shortest gap between the twopunches (i.e., (B) in Fig. 2), (ii) the minimum relative density of the compact inside the die ρ in-diemin which occurs right after separation of the upper punch from the compact (i.e., (C) in Fig. 2), and (iii)the relative density of the tablet out of the die ρ tablet which is approximately equal to ρ in-diemin (i.e., (D) inFig. 2). Elastic Deformation Plastic Deformation Bonding & Fracture
E ν κ m K Ic Material 1 5 GPa 0.25 150 MPa 2.00 1.26 MPa m / ( ω = 150 J/m , G p = 0 )Material 2 30 GPa 0.25 900 MPa 2.00 6.19 MPa m / ( ω = 600 J/m , G p = 0 ) Table 1: Material properties.We specifically study a noncohesive frictionless [43] granular system comprised by 6,512 weight-. Gonzalez Generalized loading-unloading contact laws less spherical particles with radius R = 220 µ m, and two sets of material properties, namely (i) Mate-rial 1 with Young’s modulus E = 5 GPa, and Poisson’s ratio ν = 0 . , plastic stiffness κ = 150 MPa,plastic law exponent m = 2 . , and fracture toughness K Ic = 1 . MPa m / , and (ii) Material 2 with Young’s modulus E = 30 GPa, and Poisson’s ratio ν = 0 . , plastic stiffness κ = 900 MPa,plastic law exponent m = 2 . , and fracture toughness K Ic = 6 . MPa m / (see Table 1). Thesematerials properties do not correspond to any material in particular but rather represent lower andupper bounds for many pharmaceutical powders, including drugs and excipients (see, e.g., [44, 64]and references therein). The granular bed, which is numerically generated by means of a ballisticdeposition technique [38], is constrained by a rigid cylindrical die of diameter D = 10 mm. As-suming a sufficiently small compaction speed, we consider rate-independent material behavior andwe neglect traveling waves, or any other dynamic effect [26]. The deformation process is thereforedescribed by a sequence of static equilibrium configurations using the particle mechanics approachpresented in Section 4. In this work we employ 115 quasi-static load steps and we consider 12 unload-ing points, namely ρ in-diemax = { . , . , , . , . , . , . , . , . , . , . , . } . Figure 10 shows the compacted granular bed at ρ in-diemax = 0 . and at ρ in-diemax = 0 . . Figure 11 shows the evolution of deformations of a single particle located insidethe powder bed, where the particle deformed configuration is estimated from neighboring parti-cles’ displacements and contact radii. The similitude with the experimentally observed shape ofdie-compacted spherical granules formed from microcrystalline cellulose [51] is striking. (a) (b) Figure 10: Compacted granular bed of Material 1 at ρ in-diemax = 0 . (a) and at ρ in-diemax = 0 . (b).We investigate jamming transition, evolution of the mean mechanical coordination number (num-ber of non-zero contact forces between a particle and its neighbors) in Section 5.1, punch force anddie-wall reaction during compaction and unloading in Section 5.2, in-die elastic recovery duringunloading in Section 5.3, and residual radial pressure after unloading and ejection pressure in Sec-tion 5.4. We also investigate microstructure evolution by studying probability density function ofcontact forces as well as the anisotropic granular fabric after compaction, unloading and ejection inSection 5.5. Finally, we evaluate the evolution of bonding surface area during all stages of die com-paction in Section 5.6 and we estimate the Young’s modulus and Poisson’s ratio of the compactedsolid in Section 5.7. We close by depicting a microstructure-mediated process-structure-property-performance interrelationship of the compaction process in Section 5.8.. Gonzalez Generalized loading-unloading contact laws (a) ρ in-diemax = 0 . (b) ρ in-diemax = 0 . (c) ρ in-diemax = 0 . (d) ρ in-diemax = 0 . (e) ρ in-diemax = 0 . (f) ρ in-diemax = 0 . (g) ρ in-diemax = 0 . (h) ρ in-diemax = 0 . Figure 11: Deformed configuration of a single particle located inside the powder bed ofMaterial 1 at eight different levels of confinement ρ in-diemax = { . , . , . , . , . , . , . , . } . The mean coordination number ¯ Z evolves as a power law of the following form ¯ Z − ¯ Z c = ¯ Z ( ρ in-diemax − ρ c, ¯ Z ) θ (14)where ρ c, ¯ Z is the critical relative density, ¯ Z c is the minimal average coordination number and θ is thecritical exponent. This well-known critical-like behavior has an exponent consistent with / for dif-ferent pair-interaction contact laws, polydispercity and dimensionality of the problem [53, 54, 20]. Itis known, however, that this power law is a first order approximation to the behavior of a deformablematerial for which, as demonstrated by Gonzalez et al. [26] for elastic materials, the coordinationnumber depends on the level of compressibility, i.e., on Poisson’s ratio, of the particles and thus itsscaling behavior is not independent of material properties as previously thought. This more realisticbehavior is predicted by nonlocal contact formulations [25, 26] and it will not be the focus of this pa-per, as it was previously stated. Figure 12 shows the results obtained from the particle contact mechan-ics simulations and their best fit to equation (14). Jamming occurs at ¯ Z c = 4 . and ρ c, ¯ Z = 0 . with θ = 0 . for Material 1, and at ¯ Z c = 4 . and ρ c, ¯ Z = 0 . with θ = 0 . for Material2. The fit to numerical results is good not only near jamming but also at large relative densities. Itis worth noting that the isostatic condition for frictionless packings implies a critical coordinationnumber equal to 6 and a critical density close to 0.64. In addition, however, there exists a body ofwork that indicates that ρ c, ¯ Z depends on the protocol used for obtaining jammed configurations andon the particle-die size ratio, and that monodisperse systems are prone to crystallization [5, 7, 66, 79].Here we restrict our discussion to post-jamming behavior and to one preparation protocol. It is alsointeresting to note that the jamming transition occurs later for K Ic = 0 —cf. [27], that is ρ c, ¯ Z = 0 . with ¯ Z c ≈ , for the same preparation protocol.. Gonzalez Generalized loading-unloading contact laws Relative Density M ean C oo r d i na t i on N u m be r Relative Density M ean C oo r d i na t i on N u m be r (a) Material 1 (b) Material 2 Figure 12: Mean coordination number as a function of relative density ρ in-diemax . Solid line correspondsto the best fit of equation (14) to the mean coordination number obtained from the particle contactmechanics simulation of the granular bed (symbols). The pressures applied by the punches and the reaction at the die wall are macroscopic variablesrelevant to powder die-compaction that are effectively predicted by the particle contact mechanicssimulation. These predictions are presented in Figure 13a for solid compacts compressed at 12 differ-ent relative densities. The compaction process is dominated by plastic deformations and formation ofsolid bridges, while the unloading stage is characterized by elastic recovery and breakage of bondedsurfaces. The numerical simulation accounts for these different physical mechanisms and it success-fully predicts a residual radial stress after unloading. If there is friction between the solid compactand the die wall during the ejection stage, the residual radial stress will lead to an ejection force. Fig-ure 13b shows the evolution of the punch force during compaction, unloading and ejection (assuming,for simplicity, a friction coefficient of 1). We also note that the compaction pressure follows a powerlaw of the following form σ punch = K P ( ρ in-diemax − ρ c, ¯ Z ) β P (15)where ρ c, ¯ Z is obtained from the evolution of ¯ Z , and the coefficients K P = 210 MPa and β P = 1 . are best-fitted to the numerical results for Material 1— K P = 1 . GPa and β P = 1 . for Material2. It is interesting to note that a factor of two in κ translates into a factor of two in K P , as notedin [27] for K Ic = 0 . As mentioned above, the compaction curves shown in Figure 13 representlower and upper bounds for many pharmaceutical powders, including drugs and excipients [44, 64]—e.g., ammonium chloride’s compaction curve is similar to Material 1, and lactose monohydrate’scompaction curve to Material 2 [61]. The in-die elastic recovery is further investigated and two alternative definitions found in the literatureare proposed, namely the elastic recovery in terms of relative density (cid:15) ρ and in terms of tablet height (cid:15) H , that is (cid:15) ρ = ρ in-diemax − ρ in-diemin ρ in-diemax , (cid:15) H = H in-diemin − H in-diemax H in-diemax (16). Gonzalez Generalized loading-unloading contact laws Relative Density P r e ss u r e - ( F / A t [ M P a ] ) Tablet Press Step [-] P un c h F o r c e [ k N ] (a) Material 1 (b) Material 1 Relative Density P r e ss u r e - ( F / A t [ M P a ] ) Tablet Press Step [-] P un c h F o r c e [ k N ] (c) Material 2 (d) Material 2 Figure 13: Punch and die-wall pressures as a function of relative density ρ in-die (a)&(c) and punchforce (b)&(d) as a function of tablet press step during the powder die-compaction process—labelscorrespond to Fig. 2.. Gonzalez Generalized loading-unloading contact laws Naturally, these two definitions are interrelated and, adopting a linear relationship between (cid:15) ρ and rel-ative density, one obtains they are the same to fist order in (cid:15) ρ . Specifically, the following relationshipshold (cid:15) ρ = (cid:15) ρ in-diemax − ρ c,(cid:15) − ρ c,(cid:15) (17) (cid:15) H = (cid:15) ( ρ in-diemax − ρ c,(cid:15) )1 − ρ c,(cid:15) − (cid:15) ( ρ in-diemax − ρ c,(cid:15) ) = (cid:15) ρ + O ( (cid:15) ρ ) (18)where the in-die elastic recovery at full compaction (cid:15) = 3 . and the critical relative density ρ c,(cid:15) = 0 . are best-fitted to the numerical results for Material 1— (cid:15) = 4 . and ρ c,(cid:15) = 0 . for Material 2. Figure 14 shows the results obtained from the particle mechanics simulations andtheir fit to the above equations. It is worth noting that these values are in the lower range of manypharmaceutical excipients [30, 81], which highlights the ability of the proposed model to decouplethe loading and unloading parts of the compaction curve by properly choose material properties. Relative Density I n - D i e E l a s t i c R e c o v e r y [ - ] ( ρ max - ρ min )/ ρ max ( H max - H min )/ H max Relative Density I n - D i e E l a s t i c R e c o v e r y [ - ] ( ρ max - ρ min )/ ρ max ( H max - H min )/ H max (a) Material 1 (b) Material 2 Figure 14: In-die elastic recovery as a function of relative density ρ in-diemax . Solid line corresponds to thebest fit of equation (17) to the in-die elastic recovery in terms of relative density obtained from theparticle contact mechanics simulation of the granular bed (symbols). The residual radial pressure and ejection pressure are further investigated and the following relationsare proposed σ residual = σ res , ρ in-diemax ( ρ in-diemax − ρ c,e )1 − ρ c,e (19) σ ejection = µ σ res , Wρ t πD ρ in-diemax − ρ c,e − ρ c,e (20)where the residual radial radial pressure at full compaction σ res , = 9 . MPa and the critical relativedensity ρ c,e = 0 . are best-fitted to the numerical results for Material 1— σ res , = 59 . MPa and ρ c,e = 0 . for Material 2. The two equations presented above are equivalent and obtained by using. Gonzalez Generalized loading-unloading contact laws the relationship between the punch gap and in-die relative density, i.e., H in-diemax = 4 W/ ( πD ρ t ρ in-diemax ) —where ρ t is the true density of the material and W is the weight of the powder inside the die withdiameter D . Figure 15 shows the results obtained from the particle mechanics simulations and their fitto the above equations. These values are similar to those observed in many pharmaceutical excipients(see, e.g., [1, 15]). Relative Density P r e ss u r e [ M P a ] Ejection Pressure / µ Residual Radial Pressure
Relative Density P r e ss u r e [ M P a ] Ejection Pressure / µ Residual Radial Pressure (a) Material 1 (b) Material 2
Figure 15: Residual radial pressure and ejection pressure as a function of relative density ρ in-diemax . Solidlines correspond to the best fit of equations (19)-(20) to the residual radial pressure and ejectionpressure obtained from the particle contact mechanics simulation of the granular bed (symbols). In previous subsections we studied the evolution of macroscopic, effective properties of the com-paction process (i.e., punch force and die-wall reaction during compaction and unloading, in-dieelastic recovery during unloading, and residual radial pressure after unloading and ejection pressure).Next, we turn attention to the evolution of some of the microstructural features that give rise to suchmacroscopic behavior. Under increasing confinement, powders support stress by spatial rearrange-ment and deformation of particles and by the development of inhomogeneous force networks. A forcenetwork is typically characterized by the probability distribution of its inter-particle contact forces andtheir directional orientation.For simplicity of exposition, we restrict attention to the behavior of Material 1—results are similarfor Material 2. Figure 16 shows the distribution of contact forces after loading, unloading and ejectionat three different levels of compaction. It is worth noting that we show force histograms rather than,as it is typically used in the literature, probability distributions of contact forces non-dimensionalizedby their mean value. In turn, it is evident from the figure that: (i) the range of compressive forcesincreases with relative density more than the range of tensile forces does; (ii) compressive forcesreduce significantly in magnitude during unloading, while tensile forces hardly change; and (iii)distributions after unloading and after ejection are very similar and both exhibit some symmetryabout zero.Granular fabric anisotropy is a complementary aspect of stress transmission and it can be de-termined from particle mechanics descriptions of granular systems under static equilibrium. Thisanisotropy is influenced by different factors, such as particle shape, die filling protocol, and deforma-. Gonzalez Generalized loading-unloading contact laws Contact Force - [N] C oun t CompactionUnloadingEjection ρ max = 0.73 Contact Force - [N] C oun t CompactionUnloadingEjection ρ max = 0.84 Contact Force - [N] C oun t CompactionUnloadingEjection ρ max = 1.00 Figure 16: Distribution of contact forces after loading, unloading and ejection at three different levelsof compaction ρ in-diemax equal to . , . and . . Broken solid bridges are not included anddistributions are obtained from the particle contact mechanics simulation of the granular bed.tion history. We specifically study the orientation distribution function of contact normals and of themean contact force. Using spherical coordinates with azimuth and zenith angles θ and φ , we definethe contact orientation vector by n = (sin( θ ) cos( φ ) , sin( θ ) sin( φ ) , cos( θ )) and, for axial symmetry around the zenith axis or the direction of compaction [65], the sphericalharmonics spectrum of the orientation distribution function of contact normals ξ ( n ) = ξ ( θ, φ ) [9] by ξ ( θ, φ ) = 14 π (cid:16) a θ ) −
1) + a θ ) −
30 cos( θ ) + 3) (21) + a
16 (231 cos( θ ) −
315 cos( θ ) + 105 cos( θ ) − a
128 (6435 cos( θ ) − θ ) + 6930 cos( θ ) − θ ) + 35) (cid:17) with (cid:82) π (cid:82) π ξ ( θ, φ ) sin( θ )d θ d φ = 1 . Similarly, the orientation distribution function of the meancontact force is assumed to be f ( θ, φ ) f avg = 14 π (cid:18) b θ ) −
1) + b θ ) −
30 cos( θ ) + 3) (22) + b
16 (231 cos( θ ) −
315 cos( θ ) + 105 cos( θ ) − b
128 (6435 cos( θ ) − θ ) + 6930 cos( θ ) − θ ) + 35) (cid:19) with (cid:82) π (cid:82) π f ( θ, φ ) sin( θ )d θ d φ = f avg . Figures 17 and 18 show the orientation distribution functionof the mean contact force and of contact normals after compaction, unloading and ejection obtainedfrom the particle contact mechanics simulation of the granular bed at relative densities ρ in-diemax equal to . and . , respectively. It is evident from the figure that: (i) the small number of large forcesare oriented in the loading direction after compaction, while the large number of intermediate to small. Gonzalez Generalized loading-unloading contact laws forces are oriented at ± ◦ from the loading direction; (ii) the orientation distribution of contactnormals does not significantly change during unloading and ejection due to the plastic, permanentnature of the deformations; (iii) after unloading, most large, vertically oriented forces are relaxedand, after ejection, most radially oriented forces are relaxed; (iv) compressive residual forces in theejected solid compact are mostly oriented at ± ◦ from the loading direction, and there is a smallnumber of tensile residual forces that are oriented in the direction of loading; and (v) the orientationdistribution of residual mean contact forces is different for different relative densities ρ in-diemax . ρ = 0.73f avg = 1.85 ρ = 0.73f avg = 0.24 ρ = 0.73f avg = 0.08 After compaction After unloading After ejection(a) Orientation distribution function of the mean contact force f ( θ, φ ) /f avg . ρ = 0.73Z avg = 7.73 ρ = 0.73Z avg = 7.65 ρ = 0.73Z avg = 7.59 After compaction After unloading After ejection(b) Orientation distribution function of contact normals ξ ( θ, φ ) . Figure 17: Granular fabric anisotropy, adopting axial symmetry around the direction of compaction,for relative density ρ in-diemax = 0 . . Solid lines correspond to the best fit of equations (22) and(23) to the distributions obtained from the particle contact mechanics simulation of the granular bed(symbols).Table 2 shows the coefficients a i and b i determined by fitting equations (22) and (23) to the distri-butions illustrated in Figures 17-18 and obtained from the particle contact mechanics simulation of thegranular bed. It is worth noting that a eighth-order approximation is the lowest-order representationof the orientation distributions functions that captures the characteristics of directional distributionsat the desired approximation accuracy. The eighth-order expansion of ξ ( n ) takes the form ξ ( n ) = C + C ij n i n j + C ijkl n i n j n k n l + C ijklpq n i n j n k n l n p n q + C ijklpqrs n i n j n k n l n p n q n r n s where four independent coefficients a , a , a and a emerge after enforcing symmetry aboutthe zenith axis and (cid:82) π (cid:82) π ξ ( θ, φ ) sin( θ )d θ d φ = 1 —cf. equation (22). In a similar manner, fourindependent coefficients b , b , b and b characterize the eighth-order expansion of f ( n ) /f avg —cf. equation (23).. Gonzalez Generalized loading-unloading contact laws ρ = 1.00f avg = 4.04 ρ = 1.00f avg = 0.70 ρ = 1.00f avg = 0.24 After compaction After unloading After ejection(a) Orientation distribution function of the mean contact force f ( θ, φ ) /f avg . ρ = 1.00Z avg = 9.42 ρ = 1.00Z avg = 9.30 ρ = 1.00Z avg = 9.20 After compaction After unloading After ejection(b) Orientation distribution function of contact normals ξ ( θ, φ ) . Figure 18: Granular fabric anisotropy, adopting axial symmetry around the direction of compaction,for relative density ρ in-diemax = 0 . . Solid lines correspond to the best fit of equations (22) and(23) to the distributions obtained from the particle contact mechanics simulation of the granular bed(symbols). Distribution of mean contact force Distribution of contact normals ρ in-diemax = 0 . a a a a b b b b After compaction 0.9491 -0.05312 -0.006645 -0.06142 -0.6826 -0.2216 0.2513 -0.2356After unloading -0.8806 -0.09625 0.02114 -0.1343 -0.656 -0.2475 0.2506 -0.2305After ejection -0.3673 -0.1978 0.07271 -0.3248 -0.6495 -0.2555 0.2556 -0.2324 ρ in-diemax = 0 . a a a a b b b b After compaction 1.137 -0.2228 -0.005998 -0.0717 -0.8201 -0.01669 0.09396 -0.176After unloading -0.8273 -0.1439 0.1995 -0.1348 -0.776 -0.086 0.1392 -0.197After ejection -0.7724 -0.3899 0.2906 -0.2151 -0.7657 -0.1047 0.1529 -0.2025
Table 2: Material 1. The coefficients correspond to the best fit of equations (22) and (23) to thedistributions obtained from the particle contact mechanics simulation of the granular bed.. Gonzalez Generalized loading-unloading contact laws As discussed above, the quantitative elucidation of strength formation requires not only the identifica-tion of the deformation and bonding mechanisms of interest but also of the bonding surface involvedin the process. Here we assume that an upper bound for the bonding surface involved in the forma-tion of solid bridges is the particle-to-particle contact area created during compaction. Therefore, westudy the bonding surface area by defining a parameter ¯ A b that is proportional to the ratio between thetotal bonding surface and that total available surface in the powder bed, i.e., ¯ A b = ( (cid:80) a P ) / ( R N P ) with N P being the number of particles in the monodisperse bed. Specifically, we investigate the evo-lution of the bonding surface parameter ¯ A b during all stages of die compaction (see Figure 19) andwe identify that the following relationship holds ¯ A b = ¯ A b , ρ in-diemin − ρ c,b − ρ c,b (23)where the bonding surface parameter of a fully dense tablet ¯ A b , = 1 . and the critical relativedensity ρ c,b = 0 . are best-fitted to the numerical results for Material 1— ¯ A b , = 0 . and ρ c,b = 0 . for Material 2.It bears emphasis that the formation of bonding surface area is controlled by both plastic defor-mations and elastic deformations. It is known that extensive particle elasticity could cause a drasticdecrease in tablet strength, due to the breakage of solid bridges and thus reduction of bonding surfacearea. Figure 19 shows that during unloading and ejection the lost in bonding surface area is largerfor higher relative densities, i.e., for higher elastic recovery (cf. Figure 14). These trends are consis-tent with the behavior of many pharmaceutical excipients obtained from permeametry measurements[2, 52]. Finally, it is interesting to note that ρ c, ¯ Z ≈ ρ c,b . The Young’s modulus of the compacted solid, obtained from central one-third of the unloading curveusing Hooke’s law [28, 76], is given by E tablet = E (cid:20) ρ in-diemin − ρ c,E − ρ c,E (cid:21) n (24)where the Young’s modulus of a fully dense tablet E = 2 . GPa, the exponent n = 0 . andthe critical relative density ρ c,E = 0 . are best-fitted to the numerical results for Material 1— E = 9 . GPa, n = 0 . and ρ c,E = 0 . for Material 2. These values shown in Figure 20 arein agreement with those obtained for two grades on lactose, at different lubrication levels, using anultrasound transmission technique [60, 61]. Equation (24) is a semi-empirical relationship derived byPhani and Niyogi for porous solids [58]. The exponent n is regarded as a material constant dependenton particle morphology and pore geometry of the material, which is also suggested by these resultsfor which the packings are identical and the values of n are equal. Finally, it is also interesting to notethat ρ c,e ≈ ρ c,E .As point out above, the granular bed develops anisotropic mechanical properties during loading,unloading and ejection. In this section, for simplicity, we assume isotropic behavior and thus de-termine two elastic properties, i.e., Young’s modulus and Poisson’s ratio, from the unloading curve.This assumption can be relaxed and, if a transversely isotropic material is assumed, five anisotropiccontinuum properties can be determined from the loading curve [65]. The extension of this analysisto unloading and ejection stages, though beyond the scope of this work, is currently being pursued bythe author.. Gonzalez Generalized loading-unloading contact laws Tablet Press Step [-] B ond i ng S u r f a c e - ( Σ ( a P / R ) / N P ) Relative Density B ond i ng S u r f a c e - ( Σ ( a P / R ) / N P ) (a) Material 1 (b) Material 1 Tablet Press Step [-] B ond i ng S u r f a c e - ( Σ ( a P / R ) / N P ) Relative Density B ond i ng S u r f a c e - ( Σ ( a P / R ) / N P ) (c) Material 2 (d) Material 2 Figure 19: Bonding surface area. (a)&(c) bonding surface parameter ¯ A b as a function of tablet pressstep for loading, unloading and ejection; (b)&(d) bonding surface parameter as a function of relativedensity ρ in-diemin . Solid line in (b)&(d) corresponds to the best fit of equation (23) to the bonding surfaceparameter obtained from the particle contact mechanics simulation of the granular bed (symbols).. Gonzalez Generalized loading-unloading contact laws Relative Density Y oung ‘ s m odu l u s [ G P a ] Relative Density P o i ss on ‘ s r a t i o (a) Material 1 (b) Material 1 Relative Density Y oung ‘ s m odu l u s [ G P a ] Relative Density P o i ss on ‘ s r a t i o (c) Material 2 (d) Material 2 Figure 20: Young’s modulus and Poisson’s ratio. (a)&(c) Young’s modulus of the compacted solidas a function of relative density ρ in-diemin ; (b)&(d) Poisson’s ratio of the compacted solid as a functionof relative density ρ in-diemin . Solid line in (a)&(c) corresponds to the best fit of equation (24) to thebonding surface parameter obtained from the particle contact mechanics simulation of the granularbed (symbols).. Gonzalez Generalized loading-unloading contact laws In the spirit of Olson’s design framework that integrates process, structure, property, and perfor-mance [55], as well as of the Quality by Design (QbD) principles recently adopted by the U.S. Foodand Drug Administration (FDA) [41, 42], we study next the interrelationship between two processvariables (namely compaction and ejection pressures), particle-scale material properties (i.e., E , ν , κ , m and K Ic ) and a critical quality attribute of the compacted solid product (namely the tablet Young’smodulus). This relationship is derived from microstructure formation and evolution predicted withthe proposed particle mechanics approach and the generalized loading-unloading contact laws. Fig-ure 21 uses equations (15), (20), (16) and (24) to represent this interrelationship for Materials 1 and2. It is worth noting that this interrelationships can not only be used to design a material with targetedquality attributes (see, e.g., [75, 78, 59, 60, 61]) but it can also be used to control the manufacturingprocess to assure such quality attributes are achieved (see, e.g., [74]). C o m p a c t i o n P r e s s u r e [ M P a ] T ab l e t Y oung ' s M odu l u s [ % ] E j e c t i on P r e ss u r e / µ [ M P a ] Material 2Material 1
Figure 21: Interrelationship between process variables (namely compaction and ejection pressures),material properties (i.e., E , ν , κ , m and K Ic ) and critical quality attributes of the compact (namelythe tablet Young’s modulus). Materials 1 and 2 are depicted in gray and in black, respectively. We have reported three-dimensional particle mechanics static calculations that enabled us to predictmicrostructure evolution during the three most important steps of die-compaction of solid tablets,namely during compaction, unloading, and ejection [84]. Specifically, we have simulated the com-paction, inside a rigid cylindrical die, of monodisperse elasto-plastic spherical particles capable offorming solid bridges. To this end, we have developed and employed generalized loading-unloading. Gonzalez Generalized loading-unloading contact laws contact laws for elasto-plastic spheres with bonding strength. The proposed loading-unloading con-tact laws are continuous at the onset of unloading by means of a regularization term, in the spiritof a cohesive zone model, that introduces a small, controllable error in the solid bridge breakageforce and the critical contact surface. This continuity property is in sharp contrast with the behaviorof standard mechanistic loading and unloading contact theories, which exhibit a discontinuity at theonset of unloading when particles form solid bridges during plastic deformations. In addition, thesegeneralized contact laws are explicit in terms of the relative position between the particles, and areupdated incrementally to account for strain path dependency. Furthermore, the three-dimensionalparticle mechanics static calculations show that the formulation is numerically robust, efficient, andmechanistically sound.We have exemplified the effectiveness and versatility of the particle mechanics approach by study-ing two sets of material properties, which do not correspond to any material in particular but ratherrepresent lower and upper bounds for many pharmaceutical powders, including drugs and excip-ients. These simulations reveal the evolution, up to relative densities close to one, of (i) meanmechanical coordination number, (ii) punch force and die-wall reaction, (iii) in-die elastic recov-ery, (iv) ejection pressure, (v) network of contact forces and granular fabric anisotropy, (vi) bond-ing surface area, (vii) Young’s modulus and Poisson’s ratio of the compacted solid. Our resultsare quantitatively similar to those experimentally observed in many pharmaceutical formulations[1, 2, 15, 30, 44, 52, 64, 58, 60, 61, 81]. Moreover, the evolution during compaction of these processvariables (such as punch, die-wall and ejection pressures, and in-die elastic recovery) and compactattributes (such as Young’s modulus and Poisson’s ratio) is in remarkable agreement with the for-mulae reported in the literature—i.e., (semi-)empirical relationships developed over the last decades[20, 26, 27, 53, 54, 58, 60, 74]. Furthermore, these relationships have enabled the development ofmicrostructure-mediated process-structure-property-performance interrelationships for QbD productdevelopment and process control. Table 3 summarizes these findings. It is evident from the table thata small number of parameters, with well-defined physical meaning, is required to describe such evo-lution with relative density. The systematic investigation of the relationship between these parameterswith particle-level material properties (i.e., E , ν , κ , m and K Ic ), particle morphology (such as parti-cle size distribution) and process variables (such as tablet weight, dimensions and composition) is aworthwhile direction of future research—see [27] for a systematic study of particles with hardeningplastic behavior, but no elastic unloading and formation of bonding strength. It is worth noting thatthe development of these relationships and, by extension, of predictive contact mechanics formula-tions for highly confined systems, is key to better design, optimize and control many manufacturingprocesses widely used in pharmaceutical, energy, food, ceramic and metallurgical industries. Acknowledgements
The author gratefully acknowledges the support received from the National Science Foundation grantnumber CMMI-1538861, from Purdue University’s startup funds, and from the Network for Compu-tational Nanotechnology (NCN) and nanoHUB.org. The author also thanks Prof. Alberto Cuiti˜no forinteresting discussions.. Gonzalez Generalized loading-unloading contact laws Property Evolution during compaction Eqn. Parameter Material 1 Material 2Coordination number ¯ Z = ¯ Z c + ¯ Z ( ρ in-diemax − ρ c, ¯ Z ) θ (14) ¯ Z c .
366 4 . ρ c, ¯ Z . . θ . . Punch pressure σ punch = K P ( ρ in-diemax − ρ c, ¯ Z ) β P (15) K P MPa . GPa β P .
561 1 . In-die elastic recovery (cid:15) ρ = (cid:15) ρ in-diemax − ρ c,(cid:15) − ρ c,(cid:15) ≈ (cid:15) H (17) (cid:15) .
55% 4 . ρ c,(cid:15) ≈ ρ c, ¯ Z ≈ ρ c, ¯ Z Residual radial pressure σ residual = σ res , ρ in-diemax ( ρ in-diemax − ρ c,e )1 − ρ c,e (19) σ res , . MPa . MPaEjection pressure σ ejection = µ σ res , Wρ t πD ρ in-diemax − ρ c,e − ρ c,e (20) ρ c,e . . In-die minimum relative density ρ in-diemin = ρ in-diemax (1 − (cid:15) ρ ) (16)Bonding surface area ¯ A b = ¯ A b , ρ in-diemin − ρ c,b − ρ c,b (23) ¯ A b , .
030 0 . ρ c,b ≈ ρ c, ¯ Z ≈ ρ c, ¯ Z Young’s modulus of compact E tablet = E (cid:20) ρ in-diemin − ρ c,E − ρ c,E (cid:21) n (24) E . GPa . GPa n . . ρ c,E ≈ ρ c,e ≈ ρ c,e Poisson’s ratio of compact ν tablet ≈ . ≈ . Table 3: Summary of microstructure formation and evolution during compaction, unloading andejection. Material 1: E = 5 GPa, ν = 0 . , κ = 150 MPa, m = 2 , K Ic = 1 . MPa m / . Material2: E = 30 GPa, ν = 0 . , κ = 900 MPa, m = 2 , K Ic = 6 . MPa m / .. Gonzalez Generalized loading-unloading contact laws References [1] Abdel-Hamid S., and Betz G. Study of radial die-wall pressure changes during pharmaceuticalpowder compaction.
Drug development and industrial pharmacy :4, 387–395.[2] Adolfsson, A., Gustafsson, C. and Nystr¨om, C. Use of tablet tensile strength adjusted for sur-face area and mean interparticulate distance to evaluate dominating bonding mechanisms. Drugdevelopment and industrial pharmacy :6, 753–764.[3] Ahlneck, C. and Alderborn, G. Moisture adsorption and tabletting. II. The effect on tensilestrength and air permeability of the relative humidity during storage of tablets of 3 crystallinematerials. International journal of pharmaceutics :2, 143–150.[4] Alderborn G. and Nystr¨om C. Pharmaceutical powder compaction technology. Marcel DekkerInc., New York, 1996.[5] Baranau V. and Tallarek U. On the jamming phase diagram for frictionless hard-sphere packings, Soft matter , 7838.[6] Byrd R.H., Schnabel R.B., and Shultz G.A. Approximate solution of the trust region problem byminimization over two-dimensional subspaces, Mathematical Programming , 247–263.[7] Chaudhuri P., Berthier L., and Sastry S. Jamming transitions in amorphous packings of friction-less spheres occur over a continuous range of volume fractions, Physical Review Letters , 165701.[8] C¸ elik, M. Pharmaceutical Powder Compaction Technology, Second Edition, 2016.[9] Chang C. S., and Misra A. Packing structure and mechanical properties of granulates.
ASCE J.Engng Mech. Div. :5, 1077–1093.[10] Coleman T.F. and Li Y. An interior trust region approach for nonlinear minimization subject tobounds,
SIAM Journal on Optimization :2, 418–445.[11] Conn A., Gould N., and Toint P. Trust Region Methods , Society for Industrial and AppliedMathematics, 2000.[12] De Boer A.H., Bolhuis G.K., and Lerk C.F. Bonding characteristics by scanning electron mi-croscopy of powders mixed with magnesium stearate.
Powder Technology :1, 75–82.[13] Derjaguin B.V. The force between molecules. Scientific American :1, 47–53.[14] Derjaguin B.V., Abrikosova I.I., and Lifshitz E.M. Direct measurement of molecular attractionbetween solids separated by a narrow gap.
Quarterly Reviews Chemical Society :3,295–329.[15] Doelker E., and Massuelle D. Benefits of die-wall instrumentation for research and developmentin tabletting. European journal of pharmaceutics and biopharmaceutics :2, 427–444.[16] Down G.R.B., and McMullen M.I. The effect of inteparticulate friction and moisture on thecrushing strength of sodium chloride compacts. Powder Technology :2, 169–174.. Gonzalez Generalized loading-unloading contact laws [17] Du Y., Chen L., McGruer N.E., Adams G.G., Etsion I. A finite element model of loading andunloading of an asperity contact with adhesion and plasticity, Journal of Colloid and InterfaceScience :2, 522–528.[18] Du Y., Adams G.G., McGruer N.E., Etsion I. A parameter study of separation modes of adheringmicrocontacts,
Journal of Applied Physics , 064902.[19] Duberg M. and Nystr¨om C. Studies on direct compression of tablets XII. The consolidation andbonding properties of some pharmaceutical compounds and their mixtures with Avicel 105,
IntJ Pharm Tech Prod Manuf :2, 17–25.[20] Durian D.J. Foam mechanics at the bubble scale, Physical Review Letters :26, 4780–4783.[21] Etsion I., Kligerman, Y. and Y C. Unloading of an elastic-plastic loaded spherical contact, In-ternational Journal of Solids and Structures , 3716–3729.[22] Fleck N.A., Storakers B., and McMeeking R.M. The viscoplastic compaction of powders, Pro-ceedings of IUTAM Symposium Mechanics of Granular Flow and Particle Compaction
Materials Letters , 365–368.[24] Frenning, G. Towards a mechanistic contact model for elastoplastic particles at high relativedensities, Finite Elements in Analysis and Design , 56–60.[25] Gonzalez M. and Cuiti˜no A.M. A nonlocal contact formulation for confined granular systems,
Journal of the Mechanics and Physics of Solids , 333–350.[26] Gonzalez M. and Cuiti˜no A.M. Microstructure evolution of compressible granular systems un-der large deformations, Journal of the Mechanics and Physics of Solids , 44–56.[27] Gonzalez M., Poorsolhjouy P., Thomas A., Liu J., and Balakrishnan K. Statistical characteri-zation of microstructure evolution during compaction of granular systems composed of sphereswith hardening plastic behavior, arXiv preprint arXiv:1801.09156; under review Interna-tional Journal of Solids and Structures :10, 3088–3106.[29] Harthong B., J´erier J.-F., Dor´emus P., Imbault D. and Donz´e F.-V. Modeling of high-densitycompaction of granular materials by the Discrete Element Method, International Journal ofSolids and Structures , 3357–3364.[30] Haware R.V., Tho I. and Bauer-Brandl A. Evaluation of a rapid approximation method for theelastic recovery of tablets. Powder Technology :1-3, 71–77.[31] Hill R. and Storakers, B. A concise treatment of axisymmetric indentation in elasticity. In: Ea-son, G., Ogden, R.W. (Eds.), Elasticity: Mathematical Methods and Applications. Ellis Hor-wood, Chichester, pp. 199-210, 1990.. Gonzalez Generalized loading-unloading contact laws [32] Hill R. The mathematical theory of plasticity. Oxford, United Kingdom: Oxford UniversityPress; 1998.[33] Israelachvili J.N. Intermolecular and surface forces. Third Edition. Academic press, 2011.[34] Israelachvili J.N., and Tabor D. Van der Waals forces: theory and experiment. In Progress insurface and membrane science , 1–55.[35] Joesten M.D., and Schad Q. Hydrogen Bonding. Marcel Dekker Inc, 1974.[36] Johnson K.L., Contact Mechanics. Cambridge University Press, London, 1985.[37] Jonsson H., Gr˙asj¨o J., and Frenning G. Mechanical behaviour of ideal elastic-plastic particlessubjected to different triaxial loading conditions. Powder Technology
Journal of Colloid and Interface Science :1, 265–272.[39] Karehill P.G., B¨orjesson E., Glazer M., et al. Bonding surface area and bonding mechanisms-two important factors for the understanding of powder compactability.
Drug Dev Ind Pharm , 2143–2196.[40] Kinloch A.J., ed. Fracture behaviour of polymers. Springer Science & Business Media, 2013.[41] Lawrence X.Y. Pharmaceutical quality by design: product and process development, under-standing, and control. Pharmaceutical research :4, 781–791.[42] Lee S.L., O’Connor T.F., Yang X., Cruz C.N., Chatterjee S., Madurawe R.D., Moore C.M.,Lawrence X.Y., and Woodcock J. Modernizing pharmaceutical manufacturing: from batch tocontinuous production. Journal of Pharmaceutical Innovation :3, 191–199.[43] Mahmoodi F., Alderborn G., and Frenning G. Effect of lubrication on the distribution of forcebetween spherical agglomerates during compression. Powder Technology :1, 69–74.[44] Mahmoodi F., Klevan I., Nordstr¨om J., Alderborn G., and Frenning G. A comparison betweentwo powder compaction parameters of plasticity: The effective medium A parameter and theHeckel 1/K parameter.
International journal of pharmaceutics :2, 295–299.[45] Mesarovic S.D. and Fleck N.A. Spherical indentation of elastic-plastic solids.
Proc. R. Soc. A :1987, 2707–2728.[46] Mesarovic S.D. and Fleck N.A. Frictionless indentation of dissimilar elastic-plastic spheres.
International Journal of Solids and Structures , 7001–7091.[47] Mesarovic S.D. and Johnson K.L. Adhesive contact of elasticplastic spheres. Journal of theMechanics and Physics of Solids , 2009–2033.[48] Mehta A. Granular physics , Cambridge University Press, Cambridge, 2007.[49] Mitchell A.G., and Down C.O.X. Recrystallization after powder compaction.
Internationaljournal of pharmaceutics :2-3, 337–344.. Gonzalez Generalized loading-unloading contact laws [50] Mor´e J.J. and Sorensen D.C. Computing a trust region step, SIAM Journal on Scientific andStatistical Computing :3, 553–572.[51] Nordstr¨om J., Persson, A. S., Lazorova, L., Frenning, G., and Alderborn, G. The degree of com-pression of spherical granular solids controls the evolution of microstructure and bond probabil-ity during compaction. International journal of pharmaceutics :1, 3–12.[52] Nystr¨om C, Karehill PG. Studies on direct compression of tablets XVI. The use of surfacearea measurements for the evaluation of bonding surface area in compressed powders.
PowderTechnology , 201–209.[53] O’Hern C.S., Langer S.A., Liu A.J., and Nagel S.R. Random packings of frictionless particles, Physical Review Letters :7, 075507.[54] O’Hern C.S., Silbert L.E., Liu A.J., and Nagel S.R. Jamming at zero temperature and zeroapplied stress: The epitome of disorder, Physical Review E :1, 011306.[55] Olson G.B. Computational design of hierarchically structured materials, Science :5330, 1237-?1242.[56] Olsson E. and Larsson P.-L. On the Effect of Particle Size Distribution in Cold Powder Com-paction,
Journal of Applied Mechanics :5, 051017.[57] Olsson E. and Larsson P.-L. On force-displacement relations at contact between elastic-plasticadhesive bodies. Journal of the Mechanics and Physics of Solids , 1185–1201.[58] Phani K.K., and Niyogi S.K. Young’s modulus of porous brittle solids. Journal of materialsscience :1, 257–263.[59] Razavi S.M., Gonzalez M. and Cuiti˜no A.M. General and mechanistic optimal relationships fortensile strength of doubly convex tablets under diametrical compression. International journalof pharmaceutics :1, 29–37.[60] Razavi S.M., Callegari G., Drazer G. and Cuiti˜no A.M. Toward predicting tensile strength ofpharmaceutical tablets by ultrasound measurement in continuous manufacturing.
Internationaljournal of pharmaceutics :1, 83–89.[61] Razavi S.M., Gonzalez M., and Cuiti˜no A.M. Quantification of lubrication and particle sizedistribution effect on mechanical strength of tablets, arXiv preprint arXiv:1801.02577; underreview , 2018.[62] Rumpf H. Basic principles and methods of granulation: I, II.
Chemie-Ingr-Tech. , 138–144.[63] Rumpf H. The Strength of Granules and Agglomerates. Agglomeration, Interscience, New York,1962.[64] Panelli R., and Ambrozio Filho F. A study of a new phenomenological compacting equation. Powder Technology :1-3, 255–261.[65] Poorsolhjouy P., and Gonzalez M. Connecting discrete particle mechanics to continuum gran-ular micromechanics: Anisotropic continuum properties under compaction, arXiv preprintarXiv:1802.04905; under review , 2018.. Gonzalez Generalized loading-unloading contact laws [66] Schreck C.F., O’Hern C.S., and Silbert L.E. Tuning jammed frictionless disk packings fromisostatic to hyperstatic. Physical Review E , 011305.[67] Sebhatu T., Elamin A.A., and Ahlneck C. Effect of moisture sorption on tabletting characteris-tics of spray dried (15% amorphous) lactose. Pharmaceutical research :9, 1233–1238.[68] Skrinjar O., and Larsson P.-L. Cold compaction of composite powders with size ratio. ActaMaterialia :7, 1871–1884.[69] Skrinjar O. and Larsson P.-L. On the local contact behavior in regular lattices of compositepowders. Journal of Materials Processing Technology , 312–318.[70] Skrinjar O., Larsson P.-L. and Storakers, B. Local Contact Compliance Relations at Compactionof Composite Powders,
Journal of Applied Mechanics :1, 164–168.[71] Stor˚akers B., and Larsson P.-L. On Brinell and Boussinesq indentation of creeping solids. Jour-nal of the Mechanics and Physics of Solids :2, 307–332.[72] Stor˚akers B. Local Contact Behaviour of Visco-plastic Particles. Proceedings of IUTAM Sympo-sium Mechanics of Granular Flow and Particle Compaction
InternationalJournal of Solids and Structures :24, 3061–3083.[74] Su Q. Bommireddy Y., Gonzalez M., Reklaitis G., and Nagy Z.K. Variation and Risk Analysisin Tablet Press Control for Continuous Manufacturing of Solid Dosage via Direct Compaction, Proceedings of the 13th International Symposium on Process Systems Engineering PSE 2018(San Diego, CA) , July, 2018.[75] Sun C., Grant D.J.W. Influence of Crystal Structure on the Tableting Properties of SulfamerazinePolymorphs.
Pharmaceutical Research , 274-?280.[76] Swaminathan S., Hilden J., Ramey B., and Wassgren C. Modeling the formation of debossedfeatures on a pharmaceutical tablet. Journal of Pharmaceutical Innovation :3, 214–230.[77] Tsigginos, C., Strong, J., and Zavaliangos, A. On the force-displacement law of contacts be-tween spheres pressed to high relative densities. International Journal of Solids and Structures , 17–27.[78] Tye C.K., Sun C., and Amidon G.E. Evaluation of the Effects of Tableting Speed on the Re-lationships between Compaction Pressure, Tablet Tensile Strength, and Tablet Solid Fraction. Journal of Pharmaceutical Sciences Physical Review E , 031307.[80] Vicente L.N. A comparison between line searches and trust regions for nonlinear optimization, Investigao Operacional , 173–179.[81] Yohannes B., Gonzalez M., Abebe A., Sprockel O., Nikfar F., Kang S., and Cuiti˜no A.M. Therole of fine particles on compaction and tensile strength of pharmaceutical powders. PowderTechnology , 372–378.. Gonzalez Generalized loading-unloading contact laws [82] Yohannes B., Gonzalez M., Abebe A., Sprockel O., Nikfar F., Kang S., and Cuiti˜no A.M. Evo-lution of the microstructure during the process of consolidation and bonding in soft granularsolids. International journal of pharmaceutics :1, 68–77.[83] Yohannes B., Gonzalez M., Abebe A., Sprockel O., Nikfar F., Kang S., and Cuiti˜no A.M. Dis-crete particle modeling and micromechanical characterization of bilayer tablet compaction.
In-ternational Journal of Pharmaceutics529