COMSOL 6.4 - 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.

This model is inspired by the example presented in Ref. 1, where material properties for a microcrystalline cellulose (MCC) powder are obtained from experiments. A specimen of the powder is then compacted. Friction between the metal powder and the compaction tools is taken into account.

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 Properties

A series of experiments were performed in Ref. 1 for several specimens compacted to different final densities to calibrate elastic and plastic material properties of the MCC powder. During the compaction process, the powder density changes, affecting the material properties. Therefore, many 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.

It is clear that the hardening law given in Ref. 1 and the exponential law used in COMSOL Multiphysics are different. The parameters pc0, Kc, and εpvol,max are therefore chosen so to match the results given in Ref. 1. The initial location of the cap, pc0, is set to zero since loose powder undergoes negligible initial elastic loading.

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.

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.

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.

Results

Figure 2 to Figure 4 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 during decompression.

Figure 2: Distribution of von Mises stress at the middle of the compaction process.

Figure 3: Distribution of von Mises stress at the end of the compaction process.

Figure 4: Distribution of von Mises stress after decompression.

Figure 5 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.

Figure 5: Volumetric plastic strain at the end of the compaction process.

The relative density distribution at different stages of compaction processes is shown in Figure 6 to Figure 8. 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 6: Relative density at the middle of the compaction process.

Figure 7: Relative density at the end of the compaction process.

Figure 8: Relative density at different after decompression.

Figure 9 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 10. 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 9: Axial punch pressure versus axial compaction.

Figure 10: 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.

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:


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	x	0.3	1	0

Yield Stress

In the toolbar, click

and choose > .

In the window for , type Yield Stress in the text field.

In the text field, type sigmay.

Locate the section. In the text field, type 0.2955*exp(4.5642*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:


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	x	0.6	0.875	0

Yield Function Parameter

In the toolbar, click

and choose > .

In the window for , type Yield Function Parameter in the text field.

In the text field, type a1.

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

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:


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	x	0.6	0.875	0

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 R0*1.1.

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*1.25.

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

In the text field, type -H0/8.

Click

Form Union (fin)

In the window, under > click .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

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

Definitions

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 checkbox.

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.715H0x
1	2	0.715H0-0.05H0*(x-1)

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

In the text field, type m.

The die is considered as rigid and fixed, hence there is no to include them in the 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 .

Locate the section. From the list, choose .

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.

In the window for , locate the section.

From the ρ list, choose .

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).

Gravity 1

In the toolbar, click

and choose .

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(solid.lemm1.popl1.rhorel)	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	0.16	1	Young’s modulus and Poisson’s ratio
Density	rho	Rhobulk	kg/m³	Basic
Initial yield stress	sigmags	sigmay(solid.lemm1.popl1.rhorel)	Pa	Poroplastic material model
Initial void volume fraction	f0	1-Rhorel0	1	Poroplastic material model
Yield function parameter	a1yield	a1(solid.lemm1.popl1.rhorel)	1	Pressure-dependent plasticity
Initial pressure limit	pc0	0	Pa	Cap and cutoff
Initial ellipse centroid	pcc0	-0.25[MPa]	Pa	Cap and cutoff
Hardening modulus	Kc	37.5[MPa]	Pa	Cap and cutoff
Maximum plastic volumetric strain	epvolmax	-log(Hf/H0)	1	Cap and cutoff