Supporting data and methods for the multi-scale modelling of freeze-drying of microparticles in packed-beds.

A multi-scale approach can be used to simulate the drying behavior of microparticles in packed-bed. Data outcomes from discrete element method (DEM) and computational fluid dynamics (CFD) simulations can be used to estimate some relevant product characteristics, such as the porosity, tortuosity, voids in the bed and permeability which are required by the multi scale model. Data from DEM simulations are presented, with a particular focus on the influence of the model parameters, packing characteristics and inhomogeneities (wall effect and particles segregation); computational costs and scala bility are also considered. Data on the properties of packings as modeled at the macroscale are presented with regard to the thermal conductivity of gases in the Knudsen regime and effective properties of packed-beds modeled as a pseudo-homogeneous medium. A mathematical model of the freeze-drying of single microparticles and its outcomes are first presented. Data outcomes from the mathematical model at the macroscale concerning the drying behavior of microparticles in a tray and in a vial are then presented and can be used for process design. Some further data, with detailed interpretation and discussion of the presented data, can be found in the related research data article, "A multi-scale computational framework for modelling the freeze-drying of microparticles in packed-beds" (Capozzi et al., 2019).

permeability, m2

parameter expressing the dependence of , from radiation and the contact between vial bottom and tray surface, J s−1 m−2 K−1

parameter expressing the pressure dependence of , J s−1 m−2 K−1 Pa−1

parameter expressing the pressure dependence of , Pa−1

molar concentration, mol m−3

specific heat capacity, J K−1 kg−1

binary diffusion coefficient, m2 s−1

effective binary diffusion coefficient, m2 s−1

Knudsen diffusion coefficient, m2 s−1

transport coefficients of DGM equation for molar flux of species-i, m2 s−1

particle diameter, m

pore diameter within the bed, m

pore diameter within the particle, m

emissivity for radiation heat exchange, -

restitution coefficient, -

generic view factor, -

vector of total force contact acting on particle i, N

vector of normal force contact acting on particle i, N

vector of tangential force contact acting on particle i, N

view factor between the lower heating shelf and the vial side, -

view factor between upper shelf to the top boundary of the packed-bed, -

view factor between the chamber wall and the vial side, -

view factor between chamber wall to the top boundary of the packed-bed, -

equivalent shear modulus, N m−2

gravity vector, m s−2

position of the sublimation interface, m

position of the sublimation interface within the particle, m

enthalpy of sublimation, J kg−1

momentum of inertia of particle i, kg m2

heat flux at product bottom, J m−2 s−1

heat flux at product upper surface, J m−2 s−1

heat flux at vial side, J m−2 s−1

heat flux at vial top, J m−2 s−1

heat transfer coefficient due to direct conduction from the shelf to the glass at the points of contact, J s−1 m−2 K−1

heat transfer coefficient between the technical fluid and the shelf, J s−1 m−2 K−1

overall heat transfer coefficient between the heating fluid and the product at the bottom of the vial, J s−1 m−2 K−1

overall heat transfer coefficient between the heating shelf and vial bottom, J s−1 m−2 K−1

normal elastic constant, N m−1

tangential elastic constant, N m−1

constant effective distance between the bottom of the vial and the shelf, m

thickness of the packed-bed, m

molecular weight of the i-component, kg kmol−1

equivalent mass, kg

mass of particle i, kg

molar flux, mol m−2 s−1

normal vector

pressure, Pa

chamber pressure, Pa

partial pressure of component i, Pa

vapor equilibrium pressure, Pa

equivalent radius, m

vial inner radius, m

ideal gas constant, J mol−1 K−1

radius of particle i, m

radial coordinate, m

vector position of particle i, m

normal stiffness, N m−1

tangential stiffness, N m−1

thickness of the vial wall, m

temperature, K

time, s

average temperature, K

vector of total torque acting on particle i, N m

tangential vector

Herz time, s

Rayleigh time, s

maximum velocity, m s−1

normal component of relative velocity, m s−1

velocity of sublimation front, m s−1

tangential component of relative velocity, m s−1

friction coefficient, -

equivalent Young׳s modulus, N m−2

Young׳s modulus of particle i, N m−2

molar fraction of component i

axial coordinate, m

energy accommodation coefficient, -

constant defined in Eq. (11)

temperature jump distance

particle dried layer

particle frozen layer

heat capacity ratio

normal viscoelastic constant, N s m−1

tangential viscoelastic constant, N s m−1

normal displacement, m

tangential displacement, m

bed porosity

particle porosity

total porosity

thermal conductivity, W m−1 K−1

effective thermal conductivity, W m−1 K−1

thermal conductivity of gases in rarefied conditions, W m−1 K−1

thermal conductivity of gases at atmospheric pressure, W m−1 K−1

free molecular heat conductivity at 0 °C, J s−1 m−1 K−1

viscosity of gas mixture, kg m−1 s−1

ratio of the gas to solid atomic masses, -

Poisson ratio of particle i, -

density, kg m−3

Stefan-Boltzmann constant, W m−2 K−4

particle tortuosity, -

frozen fraction function defined in (Eq. (12) [1])

volume fraction, -

layer of completely dried microparticles

layer of completely or partially frozen/dried microparticles

vial wall

angular velocity vector, rad s−1

referring respectively to ΩI, ΩII and ΩIII

dried layer

frozen layer

gas

glass vial

saturation index

ice

inert gas

interface

particle

solid

water

referring respectively to Γ and Γ

Knudsen number

Reynolds number

Prandtl number

computational fluid dynamics

discrete element method

dusty-gas model

Quality by Design

representative elementary volume

trehalose, mannitol, dextran (10 kDa) and dextran (150 kDa) mixture

Specifications table

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.dib.2018.12.061.

The following is the Supplementary material related to this article Video 1..

Video 1

Demonstrative video showing the generation of a packing of microgranules.

Value of the data

Data show how multi-scale modeling can be used as a support to better understand the freeze-drying behavior of micro-particles in vials and trays.

The open source code LIGGGHTS is used for creating packing of microparticles.

The open source code OpenFOAM is used for studying packing properties.

Mathematical models are implemented in COMSOL Multiphysics for studying the drying behavior of a single granule or granules in trays and vials.

Data can be used for choosing optimal process conditions and developing robust freeze-drying cycles in the framework of Quality by Design and continuous freeze-drying context [2].

Data

Data concerning the multi-scale modeling of microparticles in vials and trays are presented, which were generated by Discrete element method (DEM), computational fluid dynamics (CFD) and mathematical modeling at the macro-scale.

The computational cost and scalability of the DEM simulations (performed using the open source code LIGGGTHS) are presented. The effect of the time-step used in the simulations, the REV size, and the packing inhomogeneities are also documented.

The effect of mesh refinement in CFD simulations (carried out using the open source code OpenFOAM) is presented.

The simulations at the macroscale were carried out using COMSOL Multiphysics. Data presented in the following sections concern freeze-drying of a single micro-granule or packed-beds of micro-granules in a vial (see also [1]) or a tray. They refer to the (i) model parameters for describing micro-granules as a pseudo-homogeneous medium, (ii) heat transferred to the product during freeze-drying, (iii) outcomes of the models in terms of drying duration, maximum product temperature, vapor flux and the position of the sublimation interface. Table 1 summarizes the available data.

Table 1

Synopsis of simulation conditions adopted, variables considered, and data presented.

	Investigation	Variables investigated	Data available
DEM simulations (LIGGGTHS)	Contact model	–	– Model description – Model parameters
	Computational cost and scalability	– Integration time - step – Number of cores	– Runtime – Speedup – Timesteps/h
	Particles extraction	– Dimension of REV	– Porosity – Number of particles
	Packing inhomogeneities (wall effects, particles segregation)	– Position in the packings – Number of particles – Particle polydispersity	– Porosity – Mean pore diameter – Particle diameter
CFD simulations (OpenFOAM)	Mesh refinement	– Number of cells	– Model description – Permeability coefficient
Macro-scale modeling (COMSOL Multiphysics)	Mathematical formulation of drying behavior of a single micro granule	– Particle diameter – Temperature	– Model description – Drying duration – Particle temperature – Vapor flux – Position of the sublimation interface
	Mathematical formulation of drying behavior of micro granules in vials or trays	– Container (vial or tray) – Shelf temperature – Chamber pressure – Particle dimension – Particle polydispersity	– Model description – Drying duration – Product temperature – Frozen fraction – resistance to water vapor – Thermal conductivity
	Parameters evaluation of the macro-scale model	– Thermal conductivity of gases in the Knudsen regime – Heat transfer through a porous medium – Effective properties of a particle (effective density, effective heat capacity, effective thermal conductivity) – Effective properties of packed-bed (effective density, effective heat capacity, effective thermal conductivity)

Experimental design, materials, and methods

Set-up of DEM simulations

DEM simulations solve Newton׳s second law of motion for translation and rotation for any particle i at any time t:where m is the mass of particle i, r its position, F the total force acting on it. I is the moment of inertia, ω the angular velocity, and t the total torque.

Contact model

The Hertz-Mindlin contact algorithm has been used to simulate the falling of micro-particle and, so, the generation of random packings [3]. The basis behind the soft spheres model is that it allows two particles to deform during a collision by means of an overlap. The overlap then allows the calculation of the frictional, plastic and elastic forces resulting from this collision; the magnitude of these forces depends on the size of the deformation or overlap, see Fig. 1.

Fig. 1

Schematic of Hertz-Mindlin contact model between two particles.

The Hertz-Mindlin model describes the total force on each particle after a collision between particle i and particle j as follows [4]:where kn and kt are the elasticity constants, γn and γt are viscoelastic damping constants, vn,r and vt,r are the normal and tangential component of relative velocity, δn is the normal displacement and δt is the tangential displacement vector between the two particles. The first term in this equation governs the normal force and the second one accounts for the tangential forces. The tangential displacement vector satisfies the Coulomb frictional limit:

The associated parameters of kn and kt are evaluated from the elastic theory:where Y*, G* and R* are the equivalent Young׳s modulus, the equivalent shear modulus and the equivalent radius of the two contacting bodies.

The values of the viscoelastic parameters γn and γt are:where β, Sn and St are evaluated from the contact collision theory as follows:where e is the coefficient of restitution.

The values of Y*, G*, R* and m* are:where ν is the Poisson ratio, and m∗ is the equivalent mass of the two bodies in contact.

Hertz-Mindlin contact model requires five key parameters, i.e., particle radius R, particle mass m, Young׳s modulus Y, shear modulus G, Poisson׳s ratio ν and the coefficient of restitution e. These parameters are provided in Table 2. It should also be noted that this contact model was also used to describe both particle-particle and particle-wall collisions. The demonstrative video shows generation of a packing of microgranules.

Table 2

Parameters used in DEM simulations for particle-particle and wall-particle interactions.

Parameter	Value	Unit	Refs.
Density of particles	926.7	kg m−3
Young׳s modulus of ice particles	93.3 × 108	N m−2	[5]
Young׳s modulus of wall	65.0 × 109	N m−2	[5]
Poisson ratio of ice particles	0.325	–	[5]
Poisson ratio of wall	0.210	–	[5]
Friction coefficient particle-particle	0.33	–	[6]
Friction coefficient particle-wall	0.50	–	[6]
Restitution coefficient particle-particle	0.4	–	[7]
Restitution coefficient particle-wall	0.5	–	[7]

Computational cost and scalability

DEM simulations were performed on the Galileo supercomputer located in CINECA and equipped with 516 nodes containing two Intel Haswell 8-core processors each, with a clock of 2.40 GHz and a RAM of 128 GB/node [8].

LIGGGHTS is an open-source DEM code derived from the molecular dynamics code LAMMPS and used for simulating granular materials. LIGGGHTS supports the addition of mesh geometries and includes granular models for modeling particle-particle and particle-walls collisions. LIGGGHTS can be run as a serial or in a parallel environment through MPI. Dynamic MPI domain decomposition was used to mitigate the load-imbalance of DEM simulations.

The speedup of a parallel implementation reads,

Table 3 reports the test cases used for evaluating the simulation speedup. Table 4 shows the test cases used for evaluating the computational costs related to the integration timestep adopted for the simulations.

Table 3

Test case: 100 k monodisperse (St. Deviation = 0) falling particles of 30 µm as diameter; total simulated time = 1.1 s, no. of cores = 4.

Test	Integration time step	No. of time steps	Runtime	Timesteps/h
1	250	4,400,000	19.9	241,540
2	200	5,500,000	25.3	237,287
3	100	11,000,000	47.0	255,234
4	50	22,000,000	96.0	250,002
5	25	44,000,000	186.8	256,894
6	10	110,000,000	644.8	186,094

Table 4

Test case: 100,000 monodisperse (St. Deviation = 0) falling particles of 30 µm as diameter; integration time step = 25 ns, total simulated time = 1.1 s (48·106 timesteps).

Test case	No. of cores	Runtime	Timesteps/h
1	1	257.0	186,747
2	2	252.4	190,207
3	4	186.8	256,894
4	8	89.5	536,471
5	16	43.4	1,105,375
6	32	34.1	1,407,494
7	48	21.9	2,194,135
8	64	16.7	2,867,110
9	128	12.2	3,926,460
10	512	3.2	15,204,827

Data in Fig. 2a and b refer to simulation runtime and speedup relative to serial runs varying the number of cores from 1 (serial) to 512. The simulations refer to 100,000 monodisperse falling particles of 30 µm and a timestep of 25 ns. Increasing the number of cores from 1 to 4 did not have a dramatic impact on the runtime, and equivalently on the speedup, because of the geometry of the system. On the other hand, from 4 to 512 cores, the scalability was very strong.

Fig. 2

(a) Simulation runtime. (b) Speedup relative to serial run by varying the number of cores used in the parallelization. Dashed line refers to the ideal speedup. (c) Simulation runtime as a function of the timestep used.

Data in Fig. 2c show that the computational cost strongly depends on the integration time step used in the simulation and the number of falling particles.

Selection of the time-step

For DEM simulations, time step should be small enough to avoid instability and provide a reliable and stable set of particle spatial location data during the simulation as a whole. A rule of thumb states that the maximum acceptable time step must be set to a value between 0.1 and 0.3 of the minimum value of the semi-empirical parameters, Hertz time and Rayleigh time [4].

The Hertz time can be evaluated from:

The Rayleigh time reads:

Each simulation has been performed using a time-step equal to 10% of the Rayleigh time; to give you an example, for monodisperse microparticles of 10 µm, the time-step was 8.62 ns, whereas for microparticles of 90 µm, was 77.60 ns. In the case of polydisperse microparticles, the time-step was usually lower, because it was chosen considering the smallest particles of the distribution.

Selection of the REV size

REV needs to be chosen sufficiently extended to catch global properties and avoid local fluctuations. If the REV was small, porosity, which was chosen as the reference property, fluctuated for different REVs located one near each other. On the other hand, if the REV size is sufficiently large, the fluctuations in the proximity of a certain point in the bed become negligible.

Data shown in Fig. 3 refer to the case of frozen microparticles consisting of a water mixture 35% w/w of trehalose, mannitol, dextran (10 kDa) and dextran (150 kDa) (TMDD) in the ratio 3:3:3:1 atomized at 48 kHz (experimental data [9] and simulation results [1])

Fig. 3

35% w/w water-TMDD (3:3:3:1) atomized at 48 kHz. (a) Number of particles as a function of REV size and (b) porosity value as a function of the number of particles. REV was chosen in the vicinity of the center of the packing, with small deviations of its center coordinates randomly chosen within the range ±0.1 mm; each dot represents a value of these neighbors REVs.

Data shown in Fig. 4 refer to the packings of monodisperse and polydisperse microparticles of 50 µm as mean diameter. The data show that REV volume should contain at least 1000 particles to avoid local fluctuations of the porosity in the proximity of a certain position in the bed, which corresponds to a REV of 0.06 mm3.

Fig. 4

Microparticles of () monodisperse and () polydisperse microparticles (σ = 5 µm) with a mean diameter of 50 µm. (a) Number of particles as a function of the REV size and (b) porosity as a function of the number of particles. REV was chosen in the vicinity of the center of the packing, with small deviations of its center coordinates randomly chosen within the range ±0.01 cm; each dot represents a value of these neighbors REVs.

Fig. 5 show a packing of 100,000 polydisperse micro-particles and REVs with different volume (from 10 × 10−3 mm3 to 1 × 10−1 mm3).

Fig. 5

Packing of polydisperse microparticles (σ = 5 µm) with a mean diameter of 50 µm and extracted REVs of different volumes.

Set-up of CFD simulations for estimating packing properties

In this section, the set-up of CFD simulations used in [1] for estimating packing properties is presented.

Simulation set-up

The fluid flow within the packed-bed was simulated by solving the continuity and Navier–Stokes equations, imposing a pressure drop in the z-direction between the inlet and outlet sections. The no-slip boundary condition was imposed at the particle surface and the symmetry condition at the remaining sides of the REV. The pressure drop across the REV was set sufficiently low to guarantee that Re < 0.1, and so Stokes regime occurs. In that conditions, the inertial term is negligible, and the estimated permeability corresponds to the true permeability of the medium.

Mesh refinement

The independence of the solution from the mesh was obtained by successive refinement of the grid within the REV. We performed this analysis by varying the number of cells from about 20,000 to 35,000,000. To give you an example, Fig. 6 shows the permeability as a function of the number of cells within the computational domain in the case of monodisperse microparticles of 50 µm. In this case, mesh independency of the solution was reached using 20,000,000 cells.

Fig. 6

Value of permeability coefficient as obtained from successive mesh refinement of the REV, in the case of monodisperse microparticles of 50 µm as diameter.

Modeling of freeze-drying of single microparticles

The heat and mass transfer in a single particle during the primary drying step have been described using a one-dimensional, axisymmetric, unsteady state model (Fig. 7). The computational domain is divided into two subdomains, namely, the dried layer Γ and the frozen layer Γ, divided by the sublimation interface H(t).

Fig. 7

Schematic of the computational domain for the freeze-drying of a single microparticle.

Heat transfer within the particle

The heat balance equation in the dried layer reads,where ρp,d, c,and κp,d are respectively the effective density, heat capacity and thermal conductivity of the dried particle (see Section 2.5.3), and M = MwNw + MinNin is the total mass flux.

The heat balance equation in the frozen layer is expressed as follows,where ρp,f, c,and κp,f are respectively the effective density, heat capacity and thermal conductivity of the frozen particle (see 2.5.3).

The heat transfer initial conditions are,

The boundary condition at the surface of the particle can be written as:

Mass transfer within the particle

The conservation of species equations on a molar basis for water vapor and inert gas in the dried layer Γ can be expressed as:

The frozen core in the particle is considered compact, so that no mass flux occurs within that subdomain.

Rarified conditions and the micrometer magnitude of pores in the particle lead to the transition regime (0.1 < Kn < 10). In this case, the expression of flux of the i-component can be obtained from the Dusty-Gas Model as follows:where and (i = w,in) are the transport coefficients depending on concentration and pressure gradients respectively, and can be written as:

The effective binary diffusion Deff is defined as:where εp and τp are the particle porosity and tortuosity, respectively. The effective Knudsen diffusivity for the i-component is:where dP∗ refers to the pore diameter within the particle.

The initial conditions for mass transfer are:

The boundary conditions on the particle surface read:

Mass and heat balance at the moving interface

A moving interface divides the dried layer Γ and the frozen layer Γ. The mass balance across the interface gives:

The heat balance across the moving boundary reads:

Combining Eqs. (31), (32) the velocity of the interface is:where .

Modeling of freeze-drying of a packed-bed of micro-particles in vials or trays

In this section, model details of freeze-drying of micro-particles in a vial or in a tray are reported; the schematics of these models are shown in Fig. 8.

Fig. 8

Schematic illustration of the freeze-drying process of a particle-based material within a vial and in a tray.

The model of packed-beds of microparticles in a vial is discussed in [1].

This model involves the following assumptions:

The sublimation interface between ΩI and ΩII is assumed to be sharp and regular;

There is no mass neither energy accumulation at the sublimation interface;

Water vapor and ice are in thermodynamical equilibrium at the interface;

The layer ΩI and ΩI are considered as a pseudo-homogeneous medium:

In the layer ΩII, temperature and vapor pressure gradients within the particles are neglected, whereas temperature and vapor pressure gradients within the bed are considered;

In the layer ΩII, vapor within the pores in the bed is in thermodynamical equilibrium with the ice in the particles, and the mass transfer resistance of the particles itself is neglected;

Water vapor and inert gas behave as ideal gas;

Water desorption during primary drying is neglected, and only ice sublimation is considered.

Heat transfer

The energy balance in the subdomain ΩI and ΩII is given by:where M = MwNw + MinNin is the total mass flux.

In the case of microparticles in vial, the energy balance in the subdomain ΩIII reads:

The heat transfer initial conditions are:

The boundary conditions are:and:

In the case of micro-particles in a vial:

Mass transfer

Mass transfer equation in the dried subdomain ΩI reads:and in the frozen subdomain is:where:

The concentration of water vapor in ΩII is determined from the local thermodynamic equilibrium between particles and the vapor within the interstitial voids of the bed:where is determined from Marti–Mauersberger correlation [10]. For φS > 0, the variation of vapor concentration within the bed reads:

Initial and boundary conditions for mass transfer are given as:and:

The velocity of the interface is given combining mass and heat balance at the moving interface:

Parameters of the macro-scale model

Thermal conductivity of gases in the Knudsen regime

The thermal conductivity of an unconfined gas is independent of its pressure, whereas strongly depends on pressure if the gas is confined in small gaps (Smoluchowski effect). The thermal conductivity of gases in rarefied conditions were calculated according to the Kaganer model:where is the thermal conductivity at atmospheric pressure, Kn is the adimensional Knudsen number and β∗ is defined as follows:

α is the so-called accommodation coefficient which considers the effectiveness in the energy transfer in the molecule-wall collision, γ is the heat capacity ratio and Pr is the adimensional Prandtl number. β* is the temperature jump distance, derived from the Smoluchowsky׳s temperature jump condition applied to the Fourier׳s heat conduction equation [11]:where Tgas and Tsol are the temperature of the gas and the solid surface respectively.

The thermal accommodation coefficient accounts for the effectiveness of the energy and momentum exchanged because of the interaction between the gas molecules and the solid surface. The accommodation coefficient was calculated using the Baule formula [12] as modified by Goodman [13]:where µ∗ is the ratio of the gas to solid atomic masses.

The thermal conductivity calculated for water vapor at 250 K as a function of its pressure in confined gaps of different size is shown in Fig. 9.

Fig. 9

κKngas to κ°gas ratio as a function of pressure. Gap size is reported over the lines.

Heat transfer through a porous medium

The mathematical model at the macroscale describes the frozen and the dried product as pseudo-homogeneous media, which require the knowledge of effective properties.

As shown in Fig. 10, the effective thermal conductivity can be depicted as a function of many heat transfer mechanisms [14], [15]:where κsol and κgas are respectively the thermal conductivity of the solid and gas, κconv is due to the convective contribution and κrad due to the scattering at interfaces and grain boundaries; these last two terms are considered negligible in the presented simulations.

Fig. 10

Schematic of the heat transfer mechanism in a porous medium.

Table 5 shows the main models for estimating the effective thermal conductivity through a porous medium, and Fig. 11 reports the normalized effective thermal conductivity values calculated according to those models [16].

Table 5

Summary of the main models for estimating the effective thermal conductivity through a porous medium.

Model	Schematic	Effective thermal conductivity
Parallel model		κ=φ_solκ_sol+φ_gasκ_gas
Maxwell-Eucken (continuous solid phase, dispersed gas phase)		κ=κ_sol2κ_sol+κ_gas−2φ_gas(κ_sol−κ_gas)2κ_sol+κ_gas+φ_gas(κ_sol−κ_gas)
Effective Medium theory model		φ_gasκ_gas−κκ_gas+2κ+φ_solκ_sol−κκ_sol+2κ=0
Maxwell-Eucken (continuous gas phase, dispersed solid phase)		κ=κ_gas2κ_gas+κ_sol-2φ_sol(κ_gas-κ_sol)2κ_gas+κ_sol+φ_sol(κ_gas-κ_sol)
Series model		1κ=φ_solκ_sol+φ_gasκ_gas

Fig. 11

Normalized effective thermal conductivity of the porous medium as a function of the gas fraction ϕgas according to the () parallel model, () Maxwell-Eucken (continuous solid phase, dispersed gas phase), () EMT, () Maxwell-Eucken (continuous gas phase, dispersed solid phase) and () series model. In (a) κsol/κgas = 10 and in (b) κsol/κgas = 100.

Particle properties

The value of physical properties used in the simulation of freeze-drying at the macroscale is shown in Table 6. The effective properties of the dried layer are considered as an average of solid and gas properties, using the particle porosity εp as a weight. Effective density ρp,d, heat capacity c and thermal conductivity κp,d of the dried particle can be expressed as follows:where ρgas is the density of gas mixture in the particle pores, c its specific heat capacity and κp,d its thermal conductivity. ρsol, c and κsol are respectively the density, the heat capacity and the thermal conductivity of the solid.

Table 6

Values of physical properties and parameters.

Parameter	Value	Unit
ρsol	1514	kg m−3
ρice	920	kg m−3
ρgl	2600	kg m−3
κsol	0.2014	W m−1 K−1
κice	2.56	W m−1 K−1
κgl	1.0014	W m−1 K−1
cpsol	1383	J kg−1 K−1
cpice	2100	J kg−1 K−1
cpgas	1617	J kg−1 K−1
cpgl	840	J kg−1 K−1
Mw	18	kg kmol−1
Min	28	kg kmol−1
∆Hs	2.84	MJ kg−1
Rg	8.314	J kmol−1 K−1
σB	5.67 × 10−8	W m−2 K−4
Rgl	12.0	mm
sgl	1.2	mm
Kc	6.354	W m−2 K−1
l_V	3.80 × 10−4	m

The effective properties of the frozen layer are considered as an average of solid and ice properties, using particle porosity εp as a weight. The effective density ρp,f, heat capacity c and thermal conductivity κp,f of the frozen particle can be expressed as follows:where ρice, c and κice are respectively the density, the heat capacity and the thermal conductivity of the ice.

Packed-bed properties

The packed bed can be considered a bidisperse porous medium, as it is characterized by the particle porosity εp and bed porosity εb. The dried subdomain ΩI consists of a bed of completely dried particles, thus the total porosity is:

The volume-averaged properties of the dried subdomain ΩI can be determined from the following expressions:where ρp,d, c, κp,d are respectively the effective density, heat capacity and thermal conductivity of dried particles and ρgas, c and κgas the properties of gas flowing through the packed-bed.

The subdomain ΩII consists of a bed of particles that can be completely frozen, completely or partially dried, or frozen with condensed ice over the particle surface. Those situations are described using the frozen fraction function φS (see Eq. (12) in Ref. [1]), and the effective particle properties are functions of the frozen fraction function φS. The total porosity in ΩII can be written as:

The volume-averaged properties of subdomain ΩII can be determined from the following expressions:where ρp, c, κp are respectively the effective density, heat capacity and thermal conductivity of dried particles and ρgas, c and κgas the properties of gas flowing through the particle bed. The value of ρp, c and κp depends on the function φS as follows:where ρp,f, ρp,d, c, c, κp,f and κp,d are calculated as reported in 2.5.3.

Heat transfer at the vial bottom

Heat transfer between the heating shelf and product is due to the heat conduction by direct contact between shelf and vial, the thermal radiation and the conduction through the rarefied gas in the gap between the vial and the shelf. According to the literature [17], the heat transfer coefficient between the shelf and the vial varies with temperature and in particular with chamber pressure, as follows:where:

The terms Kc and are respectively the heat transfer coefficient due to the direct contact between the shelf and the vial and the effective separation distance of the vial, which are independent of both temperature and pressure and can be estimated experimentally for the individual type of vial. The overall heat transfer coefficient can be expressed as follows:where is the heat transfer coefficient between heating fluid and shelf and is the glass vial resistance.

Heat transfer by radiation

Radiative heat transfer during freeze-drying involves the radiation from the heating shelves and chamber walls to the vials in the chamber. Heat flux due to radiation depends on the temperature of the product and radiant surfaces, and the view factors between these surfaces and each vial. The view factors used in this work were: Fup = 0.86, Fwt = 0.06, Fws = 0.0, Flp = 0.0203 [18].

Data analysis

Packing inhomogeneity

Packings of particles can show inhomogeneities due to the wall effects but also to particles segregation; in Fig. 12 porosity has been mapped within the entire packed-bed.

Fig. 12

Value of porosity within the bed for a packed-bed of polydisperse microparticles of 30 µm and σ = 5 µm.

Wall effects

Data in Fig. 13 show the typical wall effects that characterize the packing in a confined container. These data refer to the simulation of packings of microparticles in a container with a diameter of 0.4 mm for particles of 10 µm, 1.2 mm for 30 µm, 2 mm for 50 µm, 2.8 mm for 70 µm and 3.6 mm for 90 µm.

Fig. 13

Porosity value in the axial direction across the diameter of the packed bed for different particle sizes in a monodisperse packed-bed ( 30 µm, 50 µm, 70 µm, 90 µm).

On the other hand, the freeze-drying of microparticles is usually done in vials with a diameter of 5 to 20 mm; in that context, the wall effects are negligible as the ratio of vial diameter to particle diameter is usually higher than 100.

Particles segregation

Particle segregation can occur because of (a) sifting, (b) fluidization, (c) entrainment of particles in the airstream, (d) flotation due to vibration and (e) agglomeration [19]. The data in Figs. 14 and 15 show the particle segregation in a packed bed of granules of 50 µm, monodisperse and polydisperse (σ = 5 µm) respectively; the packings were obtained by simulating 100,000 particles falling into a cylinder of 2 mm as diameter.

Fig. 14

Packed bed of 100,000 monodisperse particles of 50 µm as diameter. (a) Porosity and (b) mean pore diameter along packing height.

Fig. 15

Packed bed of 100,000 polydisperse particles of 50 µm (σ = 5 µm) as diameter. (a) Porosity and (b) mean pore diameter along packing height.

In Fig. 16 the data of porosity along the packing height are compared for packings obtained by simulating the falling of 100,000, 200,000 and 500,000 monodisperse particles of 50 µm and polydisperse particles (σ = 5 µm).

Fig. 16

Porosity (a) 100k (b) 200k and (c) 500k (,,) monodisperse particles of 50 µm as diameter and (,,) polydisperse particles (σ = 5 µm).

Freeze-drying of a single microparticle

Results concerning the behavior of a single particle during freeze-drying are reported for an aqueous solution of mannitol (21.5% w/w) used as model product. Primary drying was performed at 10 Pa and has been supposed that the particle is completely irradiated from a surface at 253 K. Particles from 10 to 100 µm and having pore diameters from 0.5 to 10 µm have been considered. As shown in Fig. 17a, drying time for a single micro-particle depends on both particle diameter and pore diameter within the particle. In the case of a particle of 10 µm, drying time is below one minute, no matter the dimension of pores within it. On the other hand, in the case of bigger particles, the dimension of pores is relevant. In fact, particles of 100 µm and having pores of 10 µm showed drying time of 78 min, whereas the same particles with pores of 2 µm were completely dried in about 380 min. The data in Fig. 17b refer to the drying time of a particle with a diameter of 50 µm and having pores of 5 µm; drying time is plotted as a function of the temperature of the radiant surface.

Fig. 17

(a) Drying time of a single micro-particle as a function of particle diameter Dp and for different : 0.5 µm (), 2.0 µm (), 5.0 µm () and 10.0 µm (). Primary drying was carried out at 10 Pa and supplying heat by radiation from a surface at 253 K. (b) Drying time of a single micro-particle (Dp = 50 µm,  = 5 µm) as a function of the temperature of the radiative surface.

Fig. 18 shows the data of temperature, vapor flux and position of the sublimation interface during drying of a microparticle of 50 µm and having pores of 5 µm at different temperature of the radiant surface (253 K, 258 K, and 263 K) and a chamber pressure of 10 Pa.

Fig. 18

(a) Temperature, (b) vapor flux and (c) position of the sublimation interface during drying of a microparticle of Dp = 50 µm and  = 5 µm. Primary drying was carried out at 10 Pa and supplying heat by radiation from a surface at () 253 K, () 258 K and () 263 K.

Freeze-drying of microparticles in packed-bed

Case studies: microparticles in a tray

The data presented in this section refer to freeze-drying of packed-beds of microparticles in a tray and have been obtained using the model described in Section 2.4. The data presented in this section were compared with the experimental data published by Song and Yeom [20]. The condition of the experimental runs are reported in Table 7. Data refer to the atomization of albumin 3% w/w atomized by using a flow rate of 25 ml/min using compressed air at 2 bar; the particles are frozen in a pool containing liquid nitrogen.

Table 7

Variables of the test cases.

	Case (i)	Case (ii)
Tshelf	263 K	247 K	Shelf temperature
Pc	15 Pa	15 Pa	Chamber pressure
Lp	7 mm	7 mm	Product height
Dp	15 µm	15 µm	Mean particle diameter
εp	0.785	0.785	Particle porosity

Data shown in Fig. 19 refer to the product temperature as measured by using T-type thermocouples placed 2 mm depth in the product and the corresponding outcome from the mathematical model. Experimental data and model outcomes are also in agreement concerning the primary drying duration: case (i) 10 h, and case (ii) 25 h.

Fig. 19

Product temperature during primary drying performed at 15 Pa and using a shelf temperature of (a) 263 K and (b) 243 K. () experimental data [20] and () model outcomes; temperature is measured at a product depth of 2 mm.

Case studies: microparticles in a vial

The data presented in this section support the results presented in [1] and refer to freeze-drying of packed-beds of microgranules in vials. These data elucidate some peculiar aspects of freeze-drying of micro-particles such as the formation of a dried layer in the bottom part of the granular packing. Moreover, these data can be also useful for the process design in the framework of Quality by Design context.

Prediction of a bottom sublimation interface

During primary drying of frozen solutions in a vial, ice sublimates creating a single sublimation front that recedes downward until the product is completely dried. In that case, the frozen layer is usually compact and there is no formation of further sublimation interfaces, except in the case of shrinkage of the product. If the product shrinks, the cake being dried breaks away from the vial walls, and a new lateral interface might form.

Contrary, a frozen particle-based material is typically a bidisperse medium, where two porosities are present, i.e., bed porosity and particle porosity. In this case, each particle has its own sublimation interface, but sublimation can only occur if the vapor partial pressure at the particle interface is lower than the vapor equilibrium pressure; if the vapor partial pressure is higher than the vapor equilibrium pressure, the deposition of ice over the particles occurs.

The model presented in [1] describes this behavior through the frozen fraction function. The model predicts the formation of a dried layer at the bottom of the vial, due to the fact the particles at the bottom receive a sufficient amount of heat for sublimation to occur. At the same time, the vapor from the bottom particles moves upwards, increasing the local partial pressure in upper locations in the bed. At that point, the vapor pressure becomes higher than the vapor equilibrium pressure, and the vapor starts depositing over particles; the frozen fraction function is higher than 1. This behavior is clearly shown in Fig. 20a. Moreover, as a bottom dried layer is formed, the thermal conductivity decreases, creating a further resistance to heat transferred to the particles being dried, see Fig. 20b.

Fig. 20

Contour plot of the (a) frozen fraction function and (b) thermal conductivity in the bed of microparticles in a vial. Data refer to microparticles of 50 µm, Tshelf = 253 K and Pc = 10 Pa.

Drying behavior

In this section, the data of temperature and frozen fraction within the packing in the vial are reported for microgranules of 10 µm (Fig. 21), 30 µm (Fig. 22), 50 µm (Fig. 23), 70 µm (Fig. 24) and 90 µm (Fig. 25); the figures also show the sublimation interface going downwards during drying. The simulations were performed using a Tshelf of 253 K and a chamber pressure of 10 Pa. Simulations refer to central vials in the batch.

Fig. 21

Contour plot of the (a) product temperature and (b) frozen fraction during primary drying of microparticles of 10 µm. The simulation was performed setting Tshelf = 253 K and Pc = 10 Pa.

Fig. 22

Contour plot of the (a) product temperature and (b) frozen fraction during primary drying of microparticles of 30 µm. The simulation was performed setting Tshelf = 253 K and Pc = 10 Pa.

Fig. 23

Contour plot of the (a) product temperature and (b) frozen fraction during primary drying of microparticles of 50 µm. The simulation was performed setting Tshelf = 253 K and Pc = 10 Pa.

Fig. 24

Contour plot of the (a) product temperature and (b) frozen fraction during primary drying of microparticles of 70 µm. The simulation was performed setting Tshelf = 253 K and Pc = 10 Pa.

Fig. 25

Contour plot of the (a) product temperature and (b) frozen fraction during primary drying of microparticles of 90 µm. The simulation was performed setting Tshelf = 253 K and Pc = 10 Pa.

Process design

The effect of shelf temperature and chamber pressure on drying time of microparticles packed in vials was evaluated for a Tshelf in the range 238–273 K. Three cases were analyzed, i.e., monodisperse microparticles having (i) 10 (ii) 50 and (iii) 90 µm as diameter, using the structural parameters shown in the previous section.

As shown in Fig. 26a, drying time decreased sharply as the shelf temperature increased, especially for the smallest particles. For microparticles having 50 µm as diameter, drying time decreased from 56 to about 13 h as shelf temperature rose from 243 to 273 K. In the same range of temperature, for microparticles of 10 µm, drying time fell from 161 to 23 h, and from 47 to 10 h for microparticles of 90 µm. At higher temperatures, e.g., 273 K, drying time was roughly the same for the case (ii) and (iii), indicating that, at such high-temperature, the mass transfer was not the controlling mechanism in the drying process. Fig. 26b also shows the effect of shelf temperature on the maximum product temperature during primary drying. As the Tshelf increased, the maximum product temperature increased as well, and the difference in temperature among the three cases can be imputed to the differences in the product structure.

Fig. 26

Effect of shelf temperature on (a) drying time and (b) maximum product temperature in the case of monodisperse microparticles of () 10 µm, () 50 µm and () 90 µm.

The effect of chamber pressure on drying time is shown in the range 2–30 Pa, for a shelf temperature equal to 253 K. As shown in Fig. 27, case (ii) and (iii) displayed an optimal value of drying time in correspondence of chamber pressure set at 10 Pa. For lower pressures, drying time was slightly longer, e.g., at 2 Pa drying time was about 30 min longer than in the case of chamber pressure equal to 10 Pa. On the other hand, increasing the pressure over the optimal pressure of 10 Pa, simulations showed a sharp increase of drying time, e.g., 2 h longer at 15 Pa and 4 h longer at 20 Pa. That is a typical behavior in lyophilization of products in vials, where chamber pressure plays a double role; it affects the driving force for mass transfer within the dried product, but also the heat transfer coefficient for the heat transferred from the shelf to vial bottom. In fact, lowering pressure in the chamber, the pressure difference between sublimation interface and the upper surface of the product increased. In contrast, decreasing chamber pressure, the value of Kv decreased as well, and the heat flux to the bottom felt down. For instance, Kv, as referred to the central vials simulated in this work, at 2 Pa was about 11.5 W m−2 K, 17.8 W m−2 K at 10 Pa, and rose to 23.9 W m−2 K when chamber pressure was 20 Pa. In the case (i), the optimal value of pressure was 5 Pa; this is due to the fact that the packing of microparticles of 30 µm showed a mass transfer resistance much higher than that of the case (ii) and (iii). In the case (i), the optimal value of pressure was 5 Pa; this is due to the fact that the packing of microparticles of 30 µm showed a mass transfer resistance much higher than that of the case (ii) and (iii) and, although the heat supplied to the product was much lower, a further reduction of pressure was still beneficial for reducing drying duration.

Fig. 27

Effect of chamber pressure on drying time of particle-based product within vials and constituted of particles of () 10 µm, () 50 µm and () 90 µm.

Conventional vs granules in packed-bed

In this section, particle-based lyophilized products are compared with bulk products in vials in terms of resistance to vapor flow. Data in Fig. 28 show a comparison of the resistance to vapor flow for particle-based material and bulk products. The mannitol-based bulk products were produced using the conventional freezing [21], using suspended-vial freezing [22] and vacuum-induced surface freezing (VISF) [23], [24]. For the particle-based material, the frozen microparticles consisted of a water mixture 35% w/w of TMDD atomized at 48 kHz and 24 kHz [1].

Fig. 28

Resistance to water vapor as a function of product depth in the case of different freezing protocols for mannitol solutions: () VISF [23], () suspended-vial freezing [22], () conventional freezing [21], and in the case of water-TMDD particle-based products atomized at () 24 kHz and () 48 kHz [1].

In the case of mannitol solution, conventional freezing produced lyophilized products with a mean pore diameter of 30 to 50 µm; usually, high variability in the product structure is exhibited because of the intrinsic stochasticity of nucleation phenomena. In the case of suspended-vial freezing, the mean pore diameter ranges between 80 and 100 µm. On the other hand, using VISF, nucleation temperature is controlled, and consequently, the product structure; here, VISF was performed using a nucleation temperature of 253 K and a holding time of 2 h.

Subject area	Chemical and pharmaceutical engineering
More specific subject area	Discrete element method, Fluid dynamics, Freeze-drying, Mathematical modeling, Multi-scale modeling, Packed-bed, Pharmaceutical technology, Spray-freeze drying
Type of data	Model details, table, image of packing, video (mp4), text file, graph, figure
How data were acquired	DEM simulations, CFD simulations, FEM simulations, Paraview
Data format	Analyzed
Experimental factors	DEM simulation: particle dimension, polydispersity and number. CFD simulations: mesh refinement. Mathematical model at the macroscale: porosity, tortuosity, particle diameter, polydispersity, temperature, pressure.
Experimental features	Generation of packings of micro-granules, extraction of representative elementary volumes (REVs), evaluation of packing properties (porosity, tortuosity, particle diameter, mean pore diameter, permeability), evaluation of drying duration, product temperature, vapor flux, mass transfer resistance, modeling packed-beds as pseudo-homogeneous media
Data source location	Politecnico di Torino, Torino, Italy
Data accessibility	Data available with article
Related research article	L. C. Capozzi, A. Barresi, R. Pisano, A multi-scale computational framework for modelling the freeze-drying of microparticles in packed-beds, Powder Technology 343 (2019) 834–846. [1]