Vacuum Contact Drying of Frozen Millimetric Beads in a Rotary Drum—Modeling and DEM Simulation

1. Introduction

The discrete element method (DEM) was initially devoted to simulate particulate solid mechanics, especially the flow of bulk granular materials. Nowadays this method is growingly extended to multi-physics and applied to processes like contact heating, fluidized-bed drying, wet granulation, etc. In this study vacuum contact drying of small frozen spherical beads in a rotating drum with a hot wall was considered. The context is an innovative “dynamic” (or “active”) freeze-drying process for the pharmaceutical industry. The classical “static” freeze-drying in vials is a very time, energy and equipment surface consuming process. Less extensive or continuous processes have been prospected for, for instance the “dynamic” one where the product is dried as a revolving bed of free-flowing frozen spherical particles called “micro-pellets” (or “lyo-beads”).

However, freeze-drying of micro-pellets is technologically very challenging. The first issue is that the initial liquid product has to be changed into frozen solid particles. For this purpose, calibrated size droplets are sprayed (atomized) into a countercurrent stream of cold nitrogen in a so-called “prilling” tower. The frozen droplets are then collected at the bottom of the tower and conveyed in the drum of an unconventional rotary freeze-dryer. The second challenge is operating under vacuum in a rotary drum, because of the huge sealing (interlock) problems for vapor evacuation and hot fluid supply. That is why this promising technology is still not used in the industry at a large scale. In the latest paper dealing with this technology, Peng et al. [1] presented a rotary freeze-dryer design and DEM simulations of the particles flow during operation and discharge.

In other, more or less conventional spray freeze-drying processes, the frozen particles after prilling are transferred onto the shelves of a classical vacuum chamber, or in a fluidized bed vessel at atmospheric pressure, or in a stirred vessel. Good review articles on existing technical solutions and on latest trends for spray freeze-drying are available [2,3].

In pharmaceutical applications, it is extremely important that every granule or micro-pellet in the processed bed reaches the required final water content and/or does not overpass a given temperature limit. Thus, the process simulations must provide the physical state of each individual particle. The DEM is perfectly suited for this purpose as it allows tracking each solid particle all over the drying process, i.e. basically its position and velocity but by extension of the model also its temperature and water content. Particles water content and drying rate assessment by DEM simulation has received some attention in research, especially for drying of wet particles (crop grains, wood chips, pharmaceutical granules) in a hot air stream with the numerical coupling of CFD and DEM approaches [4,5]. Papers concerning specifically DEM simulation of freeze-drying of a stirred bed of frozen particles, including the particles water content tracking and drying rate evaluation, are still missing. The present paper intends to contribute to fill this gap.

2. Modeling

Modeling of drying in stirred beds had started with the continuous medium approach developed by Schlünder and Mollekopf [6], the comparison and possible transition from the continuous to the discrete approaches for heat transfer in the bed was described by Tsotsas [7]. However, the present model is centered on a single milli-bead freeze-drying and the physical ground (thermodynamics, heat and mass transfer with phase change) can be found in textbooks [8,9]. The mechanical part of the DEM model is very classical and will not be presented.

2.1. Principle and Main Assumptions

While using the DEM approach, the bed of processed material is not considered as a whole (as an equivalent continuous medium) but as an assembly of individual solid particles, where each particle has its own spatial trajectory and its own physical state history. Balance and state equations are usually written for a single particle (here milli-bead), labeled i, in contact with a given number of other beads (or a wall) labelled j. First, the classical thermal DEM equations will be recalled, starting with a single bead heat balance equation, without considering any phase change within the material. For a single bead of a given mass m_i and specific heat capacity c_i, this equation writes:

$${\mathrm{m}}_{\mathrm{i}}{\mathrm{c}}_{\mathrm{i}}{\displaystyle \frac{\mathrm{d}{\mathrm{T}}_{\mathrm{i}}}{\mathrm{d}\mathrm{t}}}=\sum _{\mathrm{j}}{\dot{\mathrm{Q}}}_{\mathrm{i}\mathrm{j}}$$

(1)

The contact heat flow rate between the bead i and the surrounding beads j (or the wall) writes:

$${\dot{\mathrm{Q}}}_{\mathrm{i}\mathrm{j}}={\mathrm{K}}_{\mathrm{i}\mathrm{j}}\left({\mathrm{T}}_{\mathrm{j}}-{\mathrm{T}}_{\mathrm{i}}\right)$$

(2)

If all the thermal conductances (K) between pairs of beads are known, and if the two equations are solved versus time simultaneously for all beads in the bed, the individual temperature histories of the beads can be obtained.

Nevertheless, drying implies a water phase change within the material, here ice sublimation within the initially frozen milli-beads. In order to take into account this phenomenon, the entire drying process at the bead scale will be split into three steps:

1) Initial heating-up step

Bead’s temperature increases from its freezing value to the sublimation saturation (solid-vapor thermodynamical equilibrium) value at the given vessel pressure. During this step there is no water loss from the bead and its water (ice) content remains constant.

2) Main sublimation step

Bead’s temperature remains constant at its saturation (thermodynamical) value. The unbounded solid water (ice) within the bead passes directly into vapor and leaves the bead at a steady rate. The bead’s ice content decreases.

3) Final desorption step

There is no more ice within the bead, there is some bounded liquid water which desorbs and evaporates. The bead’s water content decreases and tends to the residual water content at a given vessel pressure and wall temperature. The bead’s temperature increases and tends to the hot wall temperature.

Moreover, the temperature and water content gradients within the milli-beads will be neglected, with the sublimation taking place uniformly in the interior of the bead. Accordingly, the drying rate will be supposed to be controlled entirely by the heat transfer rate at the milli-bead’s surface. The dry basis water content (X) of the milli-bead will be used in the equations below, and that’s why all the thermophysical properties of the milli-bead’s material will be expressed also on a dry basis (per kg of dry matter).

2.2. Balance Equations

2.2.1. Step 1: Heating-Up

The heat balance equation of a single bead (i) during this step, when there is no phase change, is identical to equation (1). Considering that all beads are identical and have the same dry mass m_d, the equation writes:

$$m_{\mathrm{d}}{\mathrm{c}}_{\mathrm{d}\mathrm{b}}{\displaystyle \frac{\mathrm{d}{\mathrm{T}}_{\mathrm{i}}}{\mathrm{d}\mathrm{t}}}=\sum _{\mathrm{j}}{\dot{\mathrm{Q}}}_{\mathrm{i}\mathrm{j}}$$

(3)

where c_db is the specific heat capacity (dry basis) of the frozen material given by:

$$c_{db}=c_{dry}+c_{ice}{X}_{ini}$$

(4)

The water balance of a single bead in this step have just to represent the fact that there is no water removal and that the water content is constant:

$${\displaystyle \frac{\mathrm{d}{\mathrm{X}}_{\mathrm{i}}}{\mathrm{d}\mathrm{t}}}=0$$

(5)

2.2.2. Step 2: Sublimation

During this step, the heat balance equation of a single bead is quite different and plays also the role of the water balance equation since the two are strictly coupled. It is based on the fact that all of the heat flow rate provided from outside (by the contact with surrounding beads or the wall) is consumed by sublimation (and sublimation only) inside the bead:

$${-\mathrm{m}}_{\mathrm{d}}{\displaystyle \frac{\mathrm{d}{\mathrm{X}}_{\mathrm{i}}}{\mathrm{d}\mathrm{t}}}∆h_{subl}=\sum _{\mathrm{j}}{\dot{\mathrm{Q}}}_{\mathrm{i}\mathrm{j}}$$

(6)

On the left side of this equation, the amount of released water is multiplied by the phase change enthalpy gap (Δh_subl) to give the amount of heat consumed by sublimation. This equation determines the water flow rate out of the bead (drying rate) and by integration over time can provide the water content history of the single bead. The bead’s temperature remains constant during this step and equals the thermodynamical equilibrium value for sublimation at the given pressure in the drum. This condition is simply expressed by:

$${\displaystyle \frac{\mathrm{d}{\mathrm{T}}_{\mathrm{i}}}{\mathrm{d}\mathrm{t}}}=0$$

(7)

2.2.3. Step 3: Desorption

This last step starts when water content falls below à threshold value (X_des) representing the amount of water bounded to the material. In the hygroscopic domain the behavior of the material is dictated by the thermodynamical equilibriums between the solid matrix and the sorbed liquid, represented graphically by the so-called sorption-desorption isotherms or isobars. Here, the process is isobaric and thus the slope (negative) of the desorption isobar will be used to couple the thermal kinetic with the hydric kinetic of as single bead. This slope is defined by:

$$C_{des}={\left({\displaystyle \frac{\partial X}{\partial T}}\right)}_{P=cst}^{equil}$$

(8)

The above coefficient has to be identified from experimental desorption data in the considered product water content interval. If the interval is large, these data do not give a linear trend and are approximated with an average slope. The water balance equation of a single bead takes then the form:

$${\displaystyle \frac{{dX}_i}{dt}}=-\left|C_{des}\right|{\displaystyle \frac{{dT}_i}{dt}}$$

(9)

The heat balance equation of a single bead during the desorption step is much more complex. It describes the fact that the heat provided by the surroundings serves three different purposes: desorption of bounded water, vaporization of the water and warming up of the bead. During this step the bead’s temperature is rising, because the material water content is very low and thus the latent thermal effects are weak. The heat equation writes:

$${{\mathrm{m}}_{\mathrm{d}}{\mathrm{c}}_{\mathrm{d}\mathrm{b}\mathrm{i}}{\displaystyle \frac{\mathrm{d}{\mathrm{T}}_{\mathrm{i}}}{\mathrm{d}\mathrm{t}}}-\mathrm{m}}_{\mathrm{d}}∆h_{toti}{\displaystyle \frac{\mathrm{d}{\mathrm{X}}_{\mathrm{i}}}{\mathrm{d}\mathrm{t}}}=\sum _{\mathrm{j}}{\dot{\mathrm{Q}}}_{\mathrm{i}\mathrm{j}}$$

(10)

where Δh_tot is the total enthalpy total gap including desorption and vaporization:

$$∆h_{toti}=∆h_{des}\left(X_i\right)+∆h_{vap}\left(T_i\right)$$

(11)

Using equation 9, the water content variation can be eliminated and equation 10 can be rewritten with just one dependent variable, the temperature:

$${{\mathrm{m}}_{\mathrm{d}}{\mathrm{c}}_{\mathrm{d}\mathrm{b}\mathrm{i}}{\displaystyle \frac{\mathrm{d}{\mathrm{T}}_{\mathrm{i}}}{\mathrm{d}\mathrm{t}}}+\mathrm{m}}_{\mathrm{d}}∆h_{toti}\left|C_{des}\right|{\displaystyle \frac{\mathrm{d}{\mathrm{T}}_{\mathrm{i}}}{\mathrm{d}\mathrm{t}}}=\sum _{\mathrm{j}}{\dot{\mathrm{Q}}}_{\mathrm{i}\mathrm{j}}$$

(12)

In order to make this equation more compact, all the left side terms associated with the temperature variation will be grouped in an apparent specific heat capacity (dry basis) defined by:

$$c_{dbi}^{app}=c_{dry}+c_{wat}{X}_i+\Delta h_{toti}\left|C_{des}\right|$$

(13)

Using this parameter, the heat equation of a single bead writes finally:

$${\mathrm{m}}_{\mathrm{d}}{c}_{dbi}^{app}{\displaystyle \frac{\mathrm{d}{\mathrm{T}}_{\mathrm{i}}}{\mathrm{d}\mathrm{t}}}=\sum _{\mathrm{j}}{\dot{\mathrm{Q}}}_{\mathrm{i}\mathrm{j}}$$

(14)

The water and heat balance equations of all beads have to be solved simultaneously in time with the material properties being actualized at each time step.

2.3. Constitutive Equations

In order to complete the model, a constitutive equation is needed which defines the heat flow rate between a single bead and the surroundings, i.e. adjacent beads or a hot wall. As already stated, this heat flow rate writes:

$${\dot{\mathrm{Q}}}_{\mathrm{i}\mathrm{j}}={\mathrm{K}}_{\mathrm{i}\mathrm{j}}\left({\mathrm{T}}_{\mathrm{j}}-{\mathrm{T}}_{\mathrm{i}}\right)$$

(2)

Correctly evaluating the bilateral thermal conductances (K_ij) is a not trivial task. In vacuum drying processes, the pressure of the gas in the vessel is very low and heat transfer convection can be neglected. If radiation is also neglected, the bilateral conductances represent heat transfer by conduction only between two solids in direct contact. That is the assumption made here. These conductances depend on thermal conductivities (λ) of the beads but also on the number of contacts of one bead with the others during a time laps and on the mechanical momentum of the contacts. In the DEM framework, two beads in contact interpenetrate each other to a depth depending on their dynamics (trajectories and velocities) and their mechanical properties (elasticity modulus E mainly). The penetration depth (δ_ij) represents the length for heat conduction between the beads and from that the bilateral conductance can be derived.

This purely mechanistic approach, combining the Hertz’s mechanical interaction model with Fourier’s thermal conduction was widely used for DEM simulations of particulate solid contact heating processes. The equations describing this method are written below.

$$K^{ij}=2{\lambda }_{ij}\sqrt{R_{ij}{\delta }^{ij}}$$

(15)

$$\sqrt{{\delta }^{ij}}={\displaystyle \frac{1}{\sqrt{{\mathrm{R}}_{\mathrm{i}\mathrm{j}}}}}{\left[{\displaystyle \frac{3F_n^{ij}{R}_{ij}}{4E_{ij}}}\right]}^{{\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}($$

(16)

$${\lambda }_{ij}={\displaystyle \frac{{\lambda }_i{\lambda }_j}{{\lambda }_i+{\lambda }_j}}$$

(17)

$$R_{ij}={\displaystyle \frac{R_i{R}_j}{R_i+R_j}}$$

(18)

$${\displaystyle \frac{1}{E_{ij}}}={\displaystyle \frac{1-{\nu }_i^2}{E_i}}+{\displaystyle \frac{1-{\nu }_j^2}{E_j}}$$

(19)

The symbols read:

• - λ - thermal conductivity of the bead’s material,
• - R - bead’s radius,
• - F_nij - normal force between bead i and bead j, provided by the mechanical part of the DEM model,
• - E - elasticity modulus of the bead’s material,
• - ν - Poisson’s ratio of the bead’s material.

However, another, more phenomenological method exists, where the conductances are derived from the equivalent thermal conductivity of a static bed of solid particles [7]. The proper evaluation of the bead’s material thermal conductivity is also an issue, because a frozen (or wet) bead is composed of dry material and ice (water). Conductivity is not an additive property and it depends not only on volume fractions of components but also on components’ arrangement within the volume. The simple hypothesis adopted here is that the dry matter and water are arranged in a parallel way and so that the thermal conductivity varies linearly with ice (water) content. Additionally, the contribution of water vapor replacing progressively the ice during the process was neglected.

The thermal conductivity of the bead i is given by:

$${\lambda }_i=\left(1-{\epsilon }_{dry}\right){\lambda }_{dry}+{\epsilon }_{dry}{\displaystyle \frac{X_i}{X_{ini}}}{\lambda }_{ice}$$

(20)

ε_dry is the porosity of the dry bead’s material, X_ini is the initial (maximum) water content of the material when it is completely saturated with water.

2.4. Numerical Solver and Input Data

The simulations of milli-beads stirring and drying were realized with the 2017 version of the commercial software EDEM (DEM Solutions, Edinburgh, UK) but in a customized way. The bead’s temperature and water content were implemented as two additional dependent state variables using the functionality called “custom property” and modifying a preexisting EDEM routine called “calculate force”. Once the new variables were properly declared, and their variation over a time step defined in the routine, the core EDEM time integration loop handled their actualization. As concerns the mechanical interactions between beads, the default EDEM contact model (Hertz-Mindlin model) was used. The water and heat balance equations of all beads were solved simultaneously in time with the material properties being actualized at each time step. The bead’s thermal conductivity changes with its water content and thus evolves during the drying process. This effect was implemented in our simulations. However, the evolution of the bead’s elastic modulus was not implemented. The energy (enthalpy) needed for water desorption was neglected in regard of that for water vaporization.

The test product was a bed of several hundred, millimetric size, spherical beads made of microcrystalline cellulose, water saturated and frozen prior to processing. The test dryer was a small scale horizontal rotary drum with a hot peripheral wall and without baffles. The filling ratio of the drum was 12 %. The bed of milli-beads was in the rolling mode (regime) flow. These settings for filling ratio and flow regime correspond to industrial practice for tumble drying of sensitive pharmaceutical materials. The parameters values used for simulations are given in Table 1 and Table 2 below.

3. Results and Discussion

The overall (averaged over all beads) and individual (one singular bead) water content and temperature profiles were simulated. The overall drying rate and heating rate curves were derived from the average water content time profile and from the average temperature profile respectively. The results are shown on Figure 1 and Figure 2.

Statistically, the "T_min" or "X_min" and "T_max" or "X_max" values delimit the confidence interval at 100 % for temperature and water content respectively, i.e. every single bead value lies within this interval. "S_min" and "S_max" values are the lower and upper limits of a classical confidence interval at 68 % estimated based on standard deviation. The following analysis of the curves focus on the thermal and hydric dispersion within the processed bed.

According to Figure 1, all the milli-beads did not reach nearly the same temperature until the very end of the processing time. The average trends were typical of the freeze-drying process, with a long plateau corresponding to the sublimation stage and a slow rise for the desorption stage. As concerns the singular temperature profiles, during the heating-up stage the temperature in the bed was rather uniform and was quite uniform during the most part of the sublimation stage. But during the desorption step huge temperature differences between individual beads showed up. That means that the average bed temperature rise is not a good indicator for the end of the sublimation stage for the entire bed but that the heating rate drop could be a better one. The heating rate in an industrial process can be evaluated by monitoring the heating fluid outlet temperature.

According to Figure 2, all the milli-beads reached nearly the same water content after 2500 s of processing and so sooner than reaching the same temperature. The average trends were typical of the freeze-drying process, with a constant water removal rate for the sublimation stage and a falling rate for the desorption stage. It was interesting to note that the drying rate curve exhibited here clearly two falling rate periods: the first one corresponding to the increasing number of totally dry breads and the second one corresponding roughly to desorption onset.

As concerns the singular water content profiles, already from the start important water content differences developed in the bed, with peak differences at the end of the sublimation stage. That means that neither the average bed water content, nor the average drying rate, are good indicators of the end of the sublimation stage. That means also that for an industrial process, the measured vapor pressure drop in the chamber will not be a reliable sign for the end of primary drying. An additional time will be needed to ensure the hydric homogeneity for all beads.

5. Conclusions

In this study vacuum contact drying of small spherical beads in a rotating drum with a hot wall was considered. Both the mechanical and thermal interactions between particles were described according to the Hertz contact model which stipulates a small interpenetration of two adjacent particles. As concerns specifically heat transfer, only direct particle to particle conduction was considered, what may be acceptable in a vacuum process. The drying rate was supposed to be entirely controlled by heat transfer, either particle to particle or particle to wall.

The overall (averaged over all beads) and individual (one singular bead) water content and temperature profiles were simulated. The overall drying rate and heating rate were also evaluated. By extracting and plotting the time evolution of the minimal, maximal and average beads water content and temperature, the hydric and thermal dispersions within the bed were globally assessed. This kind of insight is very important because the final product homogeneity is one of the most difficult issues the pharmaceutical powders manufacturing is confronted to. For given operating conditions, DEM simulations can provide the processing time required to obtain a product which is uniformly dry, i.e. a batch where not a single particle has a water content above the limit fixed by manufacturing standards. Furthermore, the question was addressed whether global process data like average heating rate and drying rate can be used to identify the end of the sublimation stage (primary drying) for a rotating bed of small beads. The end of primary drying is a critical information for freeze-drying process control and operation.

The proposed model gave interesting insights in the physical state of the material during the process but was not yet validated against experimental data. Some data obtained on an industrial pilot scale freeze-dryer exist and will be used, if permission to disclose them is granted, in a forthcoming paper.