вход по аккаунту


Modeling of freezing step during vial freeze-drying of pharmaceuticalsЧinfluence of nucleation temperature on primary drying rate.

код для вставкиСкачать
Asia-Pac. J. Chem. Eng. 2011; 6: 288–293
Published online 19 January 2010 in Wiley Online Library
( DOI:10.1002/apj.424
Short Communication
Modeling of freezing step during vial freeze-drying
of pharmaceuticals – influence of nucleation temperature
on primary drying rate
Kyuya Nakagawa,1 Aurelie Hottot,2 Severine Vessot2 and Julien Andrieu2 *
Department of Mechanical and System Engineering, University of Hyogo, 2167 Shosha, Himeji, Hyogo 671 2280, Japan
Laboratoire Automatique etGénie des Procédés, LAGEP.UMR Q 5007 CNRS UCBL ESCPE Bât. 308 G, 43, Bd du 11 Novembre 1918, 69622
Villeurbanne, Cédex, France
Received 28 September 2009; Accepted 24 November 2009
ABSTRACT: To optimize the primary drying step of freeze-drying processes, it is very important to control
morphological parameters of the frozen matrix. A mathematical model that simulates temperature profiles during
freezing process of standard pharmaceutical formulations was set up and the ice crystal mean sizes were semi-empirically
estimated from the simulated thermal profiles. Then, water vapor mass transfer kinetics during sublimation step was
estimated from ice phase morphological parameters. All these numerical data were compared with experimental data
and a quite good agreement was observed, which confirmed the adequacy of the present model calculations. It was
confirmed that, for a given formulation, the mass transfer parameters during freeze-drying were strongly dependent on
morphological textural parameters, and consequently, on the nucleation temperatures that fix the ice phase morphology.
The influence of freezing rates was also predicted from the simulations, proving that an increase of cooling rates led
to slower primary drying rates.  2010 Curtin University of Technology and John Wiley & Sons, Ltd.
KEYWORDS: freezing; freeze-drying; ice crystal size; nucleation temperature
Ice nucleation process is well known to be a spontaneous and stochastic phenomenon. Due to the heterogeneous and stochastic nature of nucleation phenomena,
morphologies of frozen products vary from one vial to
another all through their location on the plate of the
sublimation chamber. These structural heterogeneities
mainly result in large distributions of primary drying
rates during freeze-drying.[1 – 2] Especially, for small
scale frozen system, it is considered that the nucleation
temperature is a key factor for the optimization of the
principal quality factors of the freeze-dried matrix. As
a matter of fact, the undercooling degree of the solution determines the number of nuclei and, consequently,
greatly influences the ice crystals morphology in the
frozen sample.[3 – 4]
In previous studies, we have investigated an ultrasound system that allows the ice crystal nucleation
*Correspondence to: Julien Andrieu, Laboratoire Automatique
etGénie des Procédés. LAGEP.UMR Q 5007 CNRS UCBL ESCPE
Bât. 308 G, 43, Bd du 11 Novembre 1918, 69622 Villeurbanne,
Cédex, France. E-mail:
 2010 Curtin University of Technology and John Wiley & Sons, Ltd.
Curtin University is a trademark of Curtin University of Technology
control during the freeze-drying process of model formulations in commercial glass vials.[5] As results, we
have observed that the ice crystal sizes were closely
dependent on nucleation temperatures, and that the
water vapor permeabilities of the freeze-dried cake during the sublimation step increased with the nucleation
In this article, a mathematical model of the freezing step was proposed. It allowed the simulations of
the temperature profiles of standard pharmaceutical formulations by using a commercial finite element code
taking account of actual vial geometry and of the freezedrying operating conditions. From these temperatures
profiles, a semi-empirical model was set up to estimate
the mean ice crystal sizes and, consequently, the water
vapor-dried layer permeabilities for an optimization of
the whole freeze-drying cycle.
Ten percentage (w/w) mannitol solutions prepared from
distilled water and mannitol powder (Fluka Chemie
AG) were used as model pharmaceutical formulations.
A heat exchanger combined with ultrasound transducer
Asia-Pacific Journal of Chemical Engineering
(SODEVA, France) allowed the control of the cooling
and of the nucleation processes. The nucleation of the
sample at the selected temperature was realized by ultrasound wave propagation during 1 s through the system
cooled at −1 ◦ C/min (−0.72 ◦ C/min at vial bottom) and
tubing glass vials (Verretubex, France) of 3 ml (vial
diameter d = 12 mm) were only used in this work.
Sample solution (0.75 ml) was filled in each vial (solution height h = 8.0 mm). Frozen samples were sublimated (primary drying period) with a laboratory pilot
freeze-dryer (USIFROID SMH45, France) under the
following standard sublimation conditions for therapeutic proteins: shelf temperature at −40 ◦ C, and chamber
total pressure at 10 Pa.
moderately. The freezing modeling was based again on
the heat conduction equation recalled below, by adding
two source terms, namely the heat generation due to the
ice nucleation latent heat, Qn , and the heat generation
due to the ice crystallization latent heat, Qc .
As well, it was assumed that the rate of ice nucleation
was proportional to the supercooling degree Tf − T ∗
[K], so that the positive source term Qn was expressed
as follows[6] :
Qn = Hf ki (Tf − T ∗ )
where T ∗ represents the temperature in the supercooled
liquid, and Tf , the freezing front temperature. This term
became zero as soon as the supercooled water has been
totally crystallized. Here, ki (kg/m3 s K) represents the
nucleation rate constant estimated from freezing front
velocity v (m/s).
Besides, the positive source term, Qc , corresponding
to the latent heat generated by the ice crystallization
was expressed by the classical equation:
To simulate these freezing temperature profiles, our
modeling was divided into two periods, namely the
cooling step and then the freezing step as explained
below. The 2D axi-symmetric model taking into account
the real geometry of the vial was written under FEMLAB 3.0 (COMSOL).
Firstly, cooling step was simulated to obtain the initial
condition for the temperature profile throughout the vial
bulk for the following freezing step. Thus, the cooling
step from t = 0 up to the nucleation time, noted tn , was
simulated with the classical conduction heat equation,
= ∇(k ∇T )
Qc = Hf
(ρXice )
As illustrated in Fig.1, during the ice crystal growth
period, the sample could be divided in three zones
corresponding to three different physical states, namely,
a solid zone, a mushy zone and a liquid zone. The
mushy zone is a suspension of ice in the undercooled
solution (ice fraction, Xice ) which was defined by the
temperature difference (T ) between the freezing front
Then, during the freezing step, the nucleation started at
the fixed time tn corresponding to nucleation temperature, Tn , and subsequently, ice crystal growth progressed
Sample height [m]
= ∇(k ∇T ) + Qn + Qc
Ambient temperature
280 K
Heat transfer by air contact
h1 = 6 W/m2K
Liquid zone
Mushy zone
Freezing front
Ice zone
Tf + ∆T
Controlled temperature
at -0.72 K/min
Heat transfer via
vial glass
h2 = 11 W/m2K
8 mm
4 mm
Temperature [K]
12 mm
f = 0.2 mm
Figure 1. Model of the physical formulation states through the sample.
 2010 Curtin University of Technology and John Wiley & Sons, Ltd.
Asia-Pac. J. Chem. Eng. 2011; 6: 288–293
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
and the liquid zone. ThisT parameter was a fitting
parameter, the value of which was identified from
experimental data. Moreover, we supposed that the
ice crystal growth rate was entirely controlled by the
heat transfer and that the ice fraction, Xice , varied
linearly with the temperature in the domain of phase
change. Thus, this term could be introduced in the
accumulation term be using an apparent heat capacity
value, Cp−apparent , defined as follows.[7]
: T > Tf + T
 p−liquid
 Cp−liquid + Cp−ice
+ H
T Xice :
Tf ≤ T ≤ Tf + T
Cp−ice : T < Tf
Estimation of mean ice crystal sizes
As a result of model calculation, we could obtain an
history of temperature for each freezing procedure as
a function of nucleation temperature Tn . Quantitative
estimation of ice crystal mean size, L∗ , was achieved
by using a relationship between the freezing front
rate, R, and the temperature gradient in frozen zone,
G. For each position, the time corresponding to the
freezing front passing across was determined at the
equilibrium freezing temperature (Tf ). The R and G
values were estimated from temperature profiles at
different locations from the vial bottom. Based on
literature reviews and due to our ice phase morphology
type, we adopted the following relationship type to
interpret our ice morphology data[8 – 10] :
The variations of apparent thermal conductivity and
of bulk density values as a function of temperature
were represented by a linear law all along the domain
between the liquid and the solid zone. The numerical
values that were used in this freezing model calculation
are listed in Table 1.
L∗ = aR −0.5 G −0.5
Freezing model calculation
Ice crystal size estimation
Form factor estimation
kliquid (W/m K)
kice (W/m K)
Cp−liquid (J/kg K)
Cp−ice (J/kg K)
Tf (K)
Hf (J/kg)
Xice (−)
T (K)
a (mK0.5 s−0.5 )
ε/τ (−)
333 500
Tn = -2 °C
Tn = -6 °C
Tn = -10 °C
Tn = -14 °C
L* [mm]
Position in Vial (height) [m]
where a = 12 (mK0.5 s−0.5 ) is an empirical constant
that was identified by experimental data (cf. Table 1).
Then, we could obtain the different ice crystal values
corresponding to each R and G values set. During
freezing step, these values distributed all along the
heat transfer direction and, moreover, they depended on
nucleation temperature values. Thus, one can observe on
the plots of Fig. 2(a) that the ice crystal sizes generally
increase as a function of local vial height position, but
that this size becomes smaller at the upper layer of the
Then, mean ice crystal sizes values were calculated
by averaging local ice crystal sizes throughout the whole
sample bulk; the calculated mean values, noted L∗ , are
plotted as a function of nucleation temperatures and
Table 1. Numerical values used in the calculation.
L* [mm]
Tn [°C]
Figure 2. Ice crystal mean sizes (Frozen 10% mannitol solution), (A) distribution in vial vertical direction;
(B) nucleation temperature dependency.
 2010 Curtin University of Technology and John Wiley & Sons, Ltd.
Asia-Pac. J. Chem. Eng. 2011; 6: 288–293
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
compared with the experimental values obtained from
image analysis in Fig. 2(b). Although some deviations
were observed around Tn = −2 ◦ C, most of other data
are in quite good agreement. This agreement was also
observed with Bovin Serum Albumin formulations.[11]
Permeability during freeze-drying
Dried zone permeability, noted K , could be estimated
from mean ice crystals diameters values by assuming
that the dried cake texture is represented by a bundle
of capillary tubes (diameter d p ). In our case, we could
estimate the value of Knudsen number around Kn = 4
and, consequently, from the molecular diffusion theory
in Knudsen regime, the dried layer permeability, K , by
the following equations:
Kmodel =
Dk τ
where represents the total flow contribution factor,
and the Knudsen diffusivity Dk is expressed by:
1 8Rg Ts
Dk =
· dp
3 π Mw
and the total flow contribution factor, , is equal to:
1 + dp λ
where λ represents the mean free path, given by:
kb Ts
λ= √
2π dm2 P
For this estimation, we assumed that the ice crystal
sizes corresponded to the pore diameter (d p ) in the
Eqn (8). In Eqn (7) the form factor ε/τ = 0.225 is an
empirical constant that was identified from experimental
sublimation rate data (cf. Table 1). For determining the
sublimation rates, the freeze-drying runs were stopped
after 3–4 h of sublimation that corresponded to 30–40
wt% of ice sublimated and, then, the weight loss was
Otherwise, the experimental water vapor permeability
values were estimated by the following equation
Kexp =
edry Rg Ts
Ps − Pc Mw
which assumes that the water vapor pressure decreases
linearly through the dried layer thickness, noted edry ,
and where Ps and Pc represents, respectively, the
water vapor pressure at the sublimation front and in
the sublimation chamber. Then, these experimental
 2010 Curtin University of Technology and John Wiley & Sons, Ltd.
K [10-3 m2/s]
Tn [°C]
Figure 3.
Dried layer permeability as a function of
nucleation temperature.
Kexp values are plotted as a function of nucleation
temperature in Fig. 3. Thus, we observed that the mean
Kmodel values plot represented pretty well the tendency
of variation of the corresponding experimental values
during the primary drying step.
Moreover, as previously discussed, frozen samples
presented always ice crystal sizes distributions depending on their thermal history, namely their cooling rate
and their nucleation temperature. Then, sublimation
kinetics at each sublimation time can be estimated by
freezing conditions by using Eqs (7)–(9) and the ice
crystal mean values in Fig. 2(b). For these calculations,
the value of the form factor ε/τ was still identified by
fitting several values of Kmodel and Kexp as mentioned
before Thus, we have calculated sublimation rates at
certain vial heights by using Eqn (11) with Kmodel value
deduced from ice crystal size distributions (cf. Fig. 2(a))
for two nucleation temperatures, namely −2 ◦ C and at
around −6 ◦ C, and the results are plotted in Fig. 4. It
was observed, firstly, that the experimental sublimation
rates decrease slowly with the progress of sublimation
and, secondly, that a large discrepancy between Kexp
and Kmodel values existed at the beginning and during
the first half of sublimation period with, nevertheless,
a better agreement in the second half of sublimation
period. This discrepancy could be due to a crust effect
at the top of freeze-dried layer, the effect resulting from
the cryoconcentration of the solute (mannitol) at the top
of frozen sample.
Consequently, it was quantitatively confirmed that
nucleation temperatures are key parameters that determined the ice morphology and, consequently, the sublimation times during freeze-drying processes. However,
Asia-Pac. J. Chem. Eng. 2011; 6: 288–293
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
Tn = -2 °C
Model Tn = -2 °C
Tn = -6 °C
Cooling Rate 0.2 K/min
Model Tn = -6 °C
K [10-3 m2/s]
Drying Rate [10-3kg/m2s]
0.72 K/min
5.0 K/min
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Amount Dried [kg-sublimed water/kg-initial water]
Figure 4. Drying rate kinetics during primary drying step.
cooling rate is the most experimentally accessible freezing operating parameter that controls the ice crystals
morphology. Thus, it is also worthwhile to know the
direct influence of this cooling rate on ice crystal sizes
and freeze-dried layer permeability. Figure 5 shows the
estimated dried layer permeability values simulated for
different cooling rates for mannitol samples by using the
same values for the empirical constant reported above.
These simulations clearly showed that an increase of
cooling rates led to smaller permeability values of the
dried layer. Moreover, it is noteworthy that this tendency becomes much more significant by increasing the
nucleation temperature because the permeability dependency on the nucleation temperatures becomes much
higher by decreasing the cooling rate. For example,
when samples were prepared at cooling rate equal to
−5 K/min, which corresponded to a fast cooling rate,
the variation of nucleation temperature did not result in
too large variations of permeability. The physical explanation of this behavior certainly relies on the complex
relationships that exist between cooling rates, internal
thermal gradients, crystals growth and nucleation rates
that took place all along the freezing step in the undercooled solution inside the glass vials of small size.
A physical modeling of standard pharmaceutical formulations freezing processes was proposed. Ice crystals
mean sizes were estimated from the obtained temperature profiles and they were found in fine agreement with
experimental values. We also estimated the ice crystal
sizes distribution in the axial position of the vial along
the heat flux direction. Based on molecular diffusion
 2010 Curtin University of Technology and John Wiley & Sons, Ltd.
Tn [°C]
Figure 5. Simulation of influence of cooling rates on dried
layer permeability.
theory, dried layer permeability values were estimated
from simulated ice crystal sizes and we still observed
a satisfactory agreement between experimentally determined mean permeability values and calculated ones.
Sublimation kinetics was also evaluated from estimated
permeability values. A discrepancy between experimental sublimation rates and simulated ones was observed
at the beginning of sublimation period possibly due to
a crust effect induced by the solute cryoconcentration
at the top of frozen sample. Moreover, permeability
dependency on cooling rate was also predicted, that is
to say, when the samples were frozen at faster cooling rates, the variation of nucleation temperature does
not result in too large distribution of freeze-dried cake
permeability values. Thus, it was experimentally and
theoretically confirmed that the mass transfer parameters during freeze-drying are strongly dependent on
morphological parameters of the ice crystals phase and,
consequently, on the nucleation temperatures.
a Empirical parameter in Eqn (6)
Cp Heat capacity (J/kg K)
Dk Knudsen molecular diffusion coefficient (m2 /s)
dm Diameter of water molecule (m)
dp Pore diameter (m)
d p Ice crystal mean diameter (m)
edry Dried layer thickness (m)
G Temperature gradient in frozen zone (K/m)
Hf Latent heat of crystallization (J/kg)
Asia-Pac. J. Chem. Eng. 2011; 6: 288–293
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
Convective heat transfer coefficient (W/m2 K)
Dried layer permeability (m2 /s)
Thermal conductivity (W/m K)
Boltzmann constant (J/K)
Nucleation rate constant (kg/m3 s K)
Knudsen number (−)
Ice crystal mean size (m)
Molecular weight of water (kg/kmol)
Water vapor sublimation rate (kg/m2 s)
Total pressure (Pa)
Chamber total pressure (Pa)
Sublimation front total pressure (Pa)
Latent heat of crystallization (W)
Latent heat of nucleation (W)
Freezing front rate (m/s)
Perfect gas constant (J/kmol K)
Thickness of undercooled zone (m)
Temperature (K)
Time (s)
Temperature in supercooled liquid (K)
Homogeneous undercooled temperature (K)
Freezing front temperature (K)
Nucleation temperature (◦ C)
Sublimation front temperature (K)
Ice fraction
Freezing front velocity (m/s)
Greek letters
ε Porosity of dried layer (−)
λ Water molecules mean free path (m)
ρ Density (kg/m3 )
τ Tortuosity factor (−)
Total flow contribution (−)
[1] M. Kochs, C.H. Korber, I. Heschel, B. Nunner. The influence
of the freezing process on vapour transport during sublimation
in vacuum-freeze-drying. Int. J. Heat Mass Transfer, 1991; 34,
[2] M. Kochs, C.H. Korber, I. Heschel, B. Nunner. The influence
of the freezing process on vapour transport during sublimation
in vacuum-freeze-drying of macroscopic samples. Int. J. Heat
Mass Transfer, 1993; 36, 1727–1738.
[3] J.A. Searles, J.F. Carpenter, T.W. Randolph. The ice nucleation temperature determines the primary drying rate of
lyophilization for samples frozen on a temperature-controlled
shelf. J. Pharm. Sci., 2001; 90, 860–871.
[4] J.A. Searles, J.F. Carpenter, T.W. Randolph. Annealing to
optimize the primary drying rate, reduce freezing-induced
drying rate heterogeneity, and determine Tg’ in pharmaceutical
lyophilization. J. Pharm. Sci., 2001; 90, 872–887.
[5] K. Nakagawa, A. Hottot, S. Vessot, J. Andrieu. Influence of
controlled nucleation by ultrasounds on ice morphology of
frozen formulations for pharmaceutical proteins freeze-drying.
Chem. Eng. Process, 2006; 45, 783–791.
[6] F.J.F. Qin, J.C. Zhao, A.B. Russel, X.D. Chen, J.J. Chen,
L. Robertson. Simulation and experiment of the unsteady heat
transport in the onset time of nucleation and crystallization of
ice from the subcooled solution. Int. J. Heat Mass Transfer,
2003; 46, 3221–3231.
[7] V.J. Lunardini. Finite Difference Method for Freezing and
Thawing in Heat Transfer in Cold Climates’, Van Nostrand
Reinhold Company: New York 1981.
[8] J.L. Bomben, C.J. King. Heat and mass transport in the
freezing of apple tissue. J. Food Technol., 1982; 17, 615–632.
[9] D.S. Reid. Cryomicroscope studies of the freezing of model
solutions of cryobiological interest. Cryobiology, 1984; 21,
[10] B. Woinet, J. Andrieu, M. Laurent, S.G. Min. Experimental
and theoretical study of model food freezing. Part II:
characterization and modeling of the ice crystal size. J. Food
Eng., 1998; 35, 395–407.
[11] K. Nakagawa, A. Hottot, S. Vessot, J. Andrieu. Modelling of
freezing step during freeze-drying of drugs in vials. AIChE J.,
2007; 53, 1362–1372.
apparent Apparent value
liquid Liquid zone
ice Frozen zone
 2010 Curtin University of Technology and John Wiley & Sons, Ltd.
Asia-Pac. J. Chem. Eng. 2011; 6: 288–293
DOI: 10.1002/apj
Без категории
Размер файла
124 Кб
vial, nucleation, pharmaceuticalsчinfluence, step, temperature, freezing, modeling, rate, primary, drying, freeze
Пожаловаться на содержимое документа