Pharmaceutical Tableting Process

Powder compaction is a widely adopted manufacturing process in the ceramic, automotive, and pharmaceutical industries due to its high flexibility, high material utilization, and better control over quality.

The Capped Drucker–Prager (DPC) model is popular for modeling the compaction processes of pharmaceutical powders since it is relatively easy to characterize the material parameters from experimental data.

The model is inspired by the example presented in Ref. 1, where microcrystalline cellulose (MCC) powder is compacted, and the constitutive material properties are obtained from experiments. Friction between the metal powder and the compaction tools is taken into account.

The material properties are considered density dependent, and since the formulation of the DPC model presented in Ref. 1 differs from the one in COMSOL Multiphysics, a material property mapping is applied before using these parameters directly.

Model Definition

The geometry of the workpiece (pharmaceutical powder), punches, and die are shown in Figure 1. The actual compaction process needs two punches: a fixed bottom punch and a moving top punch. Because the bottom punch and die are fixed and rigid, they are not explicitly modeled. The top punch is modeled with a moving rigid material. Due to axial symmetry, the size of the model can be reduced.

Figure 1: Geometry of the workpiece (pharmaceutical powder), punches, and die.

Material Property Mapping

A series of experiments to calibrate elastic and plastic material properties were performed for several specimens (Ref. 1) so that the samples were compacted to different final densities. During the compaction process, the powder density changes, affecting the material properties. Therefore, all material properties are expressed in terms of the relative density of the powder.

Young’s modulus and Poisson’s ratio are given as functions of the relative density in Ref. 1. But since the variation of Poisson’s ratio with respect to changes in relative density is small, a constant Poisson’s ratio of 0.16 is used instead.

The DPC model formulation presented in Ref. 1 is different than the one in COMSOL Multiphysics. Hence, a material property mapping is needed before directly using these parameters in COMSOL Multiphysics.

Capped Drucker–Prager

In Ref. 1, the Drucker–Prager yield function Fc and the plastic potential Qc are defined as

where

is the von Mises equivalent stress, p = I1/3 is the hydrostatic pressure, β is the angle of internal friction, and d is the cohesion. Here, I1 is the first stress invariant and J2 is the second deviatoric stress invariant.

The elliptic cap function Fcap and plastic potential Qcap are

where qa = d + patan β, and f = l + α − α/cos β. The eccentricity of the ellipse,

, is given by

The Drucker–Prager cone is connected with the elliptic cap using a transition surface that serves as a smooth transition between the cone and cap surfaces. By using a transition surface, it is possible to control the variables pa, pb, and R independently.

The hardening law is given by

where εpvol is the volumetric plastic strain, and A and B are the material parameters calibrated from experimental data.

The relation between the current relative density RD and the volumetric plastic strain is given by

where RD0 is the initial relative density.

Formulation in COMSOL Multiphysics

The formulation of the DPC model implemented in COMSOL Multiphysics is related to the material parameters used in Ref. 1.

The Drucker–Prager yield function Fc and the plastic potential Qc are

where α and k are the Drucker–Prager parameters. The relation between the Drucker–Prager parameters and the parameters given in Ref. 1 is

and

The elliptic cap function Fcap and plastic potential Qcap are then

where Ja is the ordinate in the

axis at I1 = Ia, and the eccentricity of the ellipse, R, is related to the eccentricity given in Ref. 1 by

In COMSOL Multiphysics, no transition zone is required between the Drucker–Prager cone and the elliptic cap, since there is a unique and smooth transition between the two surfaces. Hence, the variables Ia and Ja are determined from the values of the parameters α, k, R, and Ib.

The hardening law is given by

where pb0 is the initial location of the cap, Kiso is the isotropic hardening modulus, and εpvol,max is the maximum volumetric plastic strain.

It is clear that the hardening law given in Ref. 1 and in COMSOL are different. The parameters pb0, Kiso, and εpvol,max are chosen so that they match the results given in Ref. 1. The initial location of the cap, pb0, is defined as zero since loose powder undergoes negligible initial elastic loading.

For the small plastic strains model, the relation between the current relative density RD and the volumetric plastic strain is given by

For the large plastic strains model, the relation between the current relative density RD and the volumetric plastic strain is given by

where Jp is the plastic volume ratio. In this example, large strain plasticity is used.

The size of the die is the same as given in Ref. 1, which is 12.5 mm in height and 12 mm in diameter. The true density of the powder is taken as 1590 kg/m3 (Ref. 1), while the loose bulk density is taken as 360 kg/m3 in order to get a similar hardening and final relative density range as given in Ref. 1. These values give an initial relative density of 0.2264.

The penalty contact method with Coulomb friction (coefficient of friction equal to 0.1) is used to model the contact interaction between the powder and the die, as well as between the powder and punches (Ref. 1). Nonlocal plasticity is used to achieve mesh objectivity.

Boundary Conditions

The applied boundary conditions are:

•

The die and bottom punch are fixed.

•

The vertical displacement of the top punch is controlled by a parameter called para.

Results

Figure 2: von Mises stress at different stages of compaction and decompression.

Figure 2 shows the von Mises stress at the middle of the compaction process, at the end of the compaction process, and after decompression. The stress is at its maximum on the top periphery, and at its minimum on the bottom periphery for all stages of compaction. The higher and lower stress rings are visible at the top and bottom surfaces of the compacted mold, which is consistent with the experimental observations; see Ref. 1. The stress relaxes once the top punch is moved upward from the mold (decompression stage).

Figure 3 shows the volumetric plastic strain at the end of compaction for the pharmaceutical powder mold. There is a large variation in volumetric plastic strain from the bottom face to the top face, with the maximum plastic strain occurring in the top region.

The relative density distribution at different stages of compaction processes is shown in Figure 4. During all stages of compaction, the high-density zone is formed at the top periphery while a low-density zone is formed at the bottom periphery. Due to friction, a nonuniform density is observed at the powder mold, which is consistent with the experimental observations reported in Ref. 1. The decompression stage has negligible impact on the relative density, unlike the stress plot, as relative density is dependent on plastic deformation which is irreversible.

Figure 3: Volumetric plastic strain at the end of compaction.

Figure 4: Relative density at different stages of compaction and decompression.

Figure 5 shows the punch pressure versus axial compaction in the compaction and decompression process. The yielding starts occurring at the beginning of the compaction process. The curve matches the numerical and experimental results presented in the Ref. 1.

The relative density and average relative density during the compaction process are shown in Figure 6. The difference between them can be better explained by the volume ratios presented in the same plot. The average relative density of the powder is related to the plastic volume ratio, while the tablet’s relative density is related to the total volume ratio. The elastic deformation is small during the compaction process.

Figure 5: Axial punch pressure versus axial compaction.

Figure 6: Relative density and volume ratio during axial compaction.

Notes About the COMSOL Implementation

In the compaction process, the interaction between the workpiece and the die as well as the workpiece and the punches is modeled using a contact node. The die and bottom punch are assumed to be rigid due to their high stiffness compared to the powder mold. As the bottom punch and die are rigid and fixed, they do not need to be modeled explicitly. In the contact node, the workpiece is taken as the destination boundary. To model this situation, the option called is selected.

Reference

1. A. Baroutaji, S. Lenihan, and K. Bryan, “Combination of finite element method and Drucker-Prager Cap material model for simulation of pharmaceutical tableting process,” Material Science and Engineering Technology, vol. 48, no. 11, 2017.

Nonlinear_Structural_Materials_Module/Porous_Plasticity/pharmaceutical_tableting_process

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Geometry 1

Model parameters are available in text file.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file pharmaceutical_tableting_process_parameters.txt.

Young’s Modulus

In the toolbar, click

and choose .

In the window for , type Young's Modulus in the text field.

In the text field, type EE.

Locate the section. In the text field, type 111.96*exp(4.395*x).

Locate the section. In the table, enter the following settings:


Argument	Unit
x	1

In the text field, type MPa.

Locate the section. In the table, enter the following settings:


Argument	Lower limit	Upper limit	Unit
x	0.3	1	1

Drucker Prager Parameter k

In the toolbar, click

and choose .

In the window for , type Drucker Prager Parameter k in the text field.

In the text field, type Kd.

Locate the section. In the text field, type 0.2955*exp(4.5642*x)/sqrt(3).

Locate the section. In the table, enter the following settings:


Argument	Unit
x	1

In the text field, type MPa.

Locate the section. In the table, enter the following settings:


Argument	Lower limit	Upper limit	Unit
x	0.6	0.875	1

Drucker Prager Parameter alpha

In the toolbar, click

and choose .

In the window for , type Drucker Prager Parameter alpha in the text field.

In the text field, type Alpha.

Locate the section. In the text field, type tan((12.628*x+56.194)[deg])/(3*sqrt(3)).

Locate the section. In the table, enter the following settings:


Argument	Unit
x	1

In the text field, type 1.

Locate the section. In the table, enter the following settings:


Argument	Lower limit	Upper limit	Unit
x	0.6	0.875	1

Hardening Function

In the toolbar, click

and choose .

In the window for , type Hardening Function in the text field.

In the text field, type Pbh.

Locate the section. In the text field, type -KIso*log(1+x/Epvolmax).

Locate the section. In the table, enter the following settings:


Argument	Unit
x	1

In the text field, type Pa.

Locate the section. In the table, enter the following settings:


Argument	Lower limit	Upper limit	Unit
x	-Epvolmax	0	1

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type R0.

In the text field, type H0.

Rectangle 2 (r2)

Right-click and choose .

In the window for , locate the section.

In the text field, type H0/10.

Locate the section. In the text field, type H0.

Rectangle 3 (r3)

Right-click and choose .

In the window for , locate the section.

In the text field, type -H0/10.

Rectangle 4 (r4)

Right-click and choose .

In the window for , locate the section.

In the text field, type R0/4.

In the text field, type H0.

Locate the section. In the text field, type R0.

In the text field, type 0.

Click

Form Union (fin)

In the window, under click .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

In subsequent steps, the die and punch are rigid domains. Hence use the toggle button in to switch the boundaries, so that the workpiece boundaries are chosen as destination boundaries.

Definitions

Contact Pair 2 (ap2)

In the window, expand the node, then click .

In the window for , click the

button.

Contact Pair 3 (ap3)

In the window, click .

In the window for , click the

button.

Add a nonlocal integration coupling operator to compute the axial force and pressure.

Integration 1 (intop1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select Boundary 7 only.

Integration 2 (intop2)

Right-click and choose .

Add a nonlocal integration coupling operator to compute the axial compaction.

In the window for , locate the section.

Click

Select Boundary 8 only.

Locate the section. Clear the check box.

Variables 1

In the window, right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
Punchforce	intop1(-solid.sz)	N	Punch force
Punchpressure	Punchforce/A0	N/m²	Punch force
Rho	PowderMass/(A0*intop2(1))	kg/m³	Current powder density

Piecewise 1 (pw1)

In the toolbar, click

In the window for , type punchDisp in the text field.

Locate the section. Find the subsection. In the table, enter the following settings:


Start	End	Function
0	1	0.75H0x
1	2	0.75H0-0.15H0*(x-1)

Locate the section. In the text field, type 1.

In the text field, type m.

Domains 3 (die) is considered as rigid and fixed, hence there is no need to consider them in physics, only a mesh is required.

Solid Mechanics (solid)

In the window, under click .

In the window for , locate the section.

Click

Select Domains 2 and 3 only.

Linear Elastic Material 1

In the window, under click .

Porous Plasticity 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Find the subsection. From the list, choose .

In the pb0 text field, type Pb0.

Click to expand the section. From the list, choose .

In the lint,m text field, type 1.6[mm].

Contact 1

In the window, under click .

Friction 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the μ text field, type 0.1.

Rigid Material 1

In the toolbar, click

and choose .

Select Domain 3 only.

Prescribed Displacement/Rotation 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the w0 text field, type -punchDisp(para).

Materials

Microcrystalline Cellulose (MCC)

In the window, under right-click and choose .

In the window for , type Microcrystalline Cellulose (MCC) in the text field.

Select Domain 2 only.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Young’s modulus	E	EE(nojac(solid.lemm1.popl1.rhorel))	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	Nu	1	Young’s modulus and Poisson’s ratio
Density	rho	Rho	kg/m³	Basic
Initial void volume fraction	f0	F0	1	Poroplastic material model
Drucker-Prager alpha coefficient	alphaDrucker	Alpha(nojac(solid.lemm1.popl1.rhorel))	1	Drucker-Prager
Drucker-Prager k coefficient	kDrucker	Kd(nojac(solid.lemm1.popl1.rhorel))	Pa	Drucker-Prager
Isotropic hardening modulus	Kiso	KIso	N/m²	Mohr-Coulomb
Maximum plastic volumetric strain	epvolmax	Epvolmax	1	Mohr-Coulomb
Ellipse aspect ratio	Rcap	Rc	1	Mohr-Coulomb