Dynamic Parameters Estimation Method: Advanced Manometric

Oct 9, 2008 - The proposed approach is able to estimate not only the temperature profile ..... of the versatility of MATLABʼs high-level programming ...
2 downloads 0 Views 195KB Size
Download PDF
Ind. Eng. Chem. Res. 2008, 47, 8445–8457

8445

Dynamic Parameters Estimation Method: Advanced Manometric Temperature Measurement Approach for Freeze-Drying Monitoring of Pharmaceutical Solutions Salvatore A. Velardi, Valeria Rasetto, and Antonello A. Barresi* Dipartimento di Scienza dei Materiali e Ingegneria Chimica, Politecnico di Torino, C.so Duca degli Abruzzi 24, Torino, 10129, Italy

The aim of this paper is to describe and validate the dynamic parameters estimation (DPE) method, a noninvasive in-line monitoring technique for the freeze-drying process of pharmaceuticals that only requires pressure measurement; simulations and experimental data of pharmaceutical solutions in vials are used to this purpose. This method has to be considered as an improvement of analogous techniques, such as the manometric temperature measurement (MTM) method, which are also based on the pressure rise analysis concept. The proposed approach is able to estimate not only the temperature profile inside the frozen product at any time during the pressure rise test, but also the position of the moving interface of sublimation, the external heat transfer coefficient, and the effective diffusivity coefficient in the dried product, which are necessary for a predictive model-based control algorithm. Fairly good estimations of the product temperature and of the transport parameters, as well as of the position of the interface of sublimation, which is directly related to the state of progression of the primary drying phase, are obtained along almost all the primary drying phase. Introduction Freeze-drying is a widespread technology that, even if quite expensive, is particularly suitable for those heat-sensitive materials such as pharmaceuticals that could be damaged by the high temperature required by traditional drying treatments. In this process, water, or another solvent, is separated by freezing the wet substance (or solution) and exposing it to a very low solvent partial pressure in the drying chamber, thus causing the direct sublimation of the ice at low temperature to vapor, leaving a porous and friable structure. Primary drying is the first step of the operation where most of the solvent is removed by sublimation, followed by secondary drying, in which desorption of the bounded water occurs.1-3 Monitoring of the primary drying phase of the lyophilization process is of paramount importance since the process should be carried out at a controlled temperature value, thus avoiding irreversible damage to the product. The position of the moving front of sublimation is another variable that would be useful to monitor in addition to product temperature, because it directly provides the state of progression of the primary drying, thus allowing the determination of the end point of the sublimation phase. These aspects are particularly interesting in the phase of process development, when the operating conditions that minimize the drying time with respect to the quality constrains on the product must be identified. To be able to measure, or to estimate, these variables is then essential in order to develop an efficient control system for the process.4 Although the insertion of thin thermocouples (or thermoresistances) in some of the vials is a widely used method to measure the product temperature during the operation, especially in laboratory-scale apparatus, it is well-known that this invasive procedure can modify the elementary phenomena of ice crystal growth and nucleation so that the monitored vials may not represent the whole batch. In fact, vials with a thermocouple * To whom correspondence should be addressed. Tel.: +39-0110904658. Fax: +39-011-0904699. E-mail: [email protected].

inserted usually have less supercooling, larger ice crystals, and faster ice sublimation and, therefore, a shorter primary drying time. Differences in the drying behavior between the vials with thermocouples inserted and the rest of the batch could be a serious problem in the step of recipe development, but the bias can become even more relevant in the clean environment of sterile production, making more problematic the scale-up.5 In addition, thermocouples are usually placed in front-row vials which typically end primary drying far earlier than the interior vials due to radiation heat transfer effects, and the thermocouple itself can alter the heat flux to the product.6,7 Thus, thermocoupleequipped vials are not generally representative both because of heat transfer effects and because of freezing-bias-induced mass transfer effects. Moreover, the use of probes inserted in the product can introduce sterility problems when dealing with several pharmaceutical products, and is not generally applicable in production plants, especially if equipped with automatic loading and unloading. Concerning the monitoring of the sublimating interface, it must be considered that the thermocouples only give the information on the temperature at the point where the tip is inserted: strong variations in the slope of the temperature curve vs time are used as an indication of the passage of the ice front and thus of the end of primary drying, if the sensor is close to the bottom of the vial; however, the variability between different vials monitored with thermocouples is very large, especially toward the end of the primary drying step, for the reasons discussed before and for the uncertainty related to the exact position of the measuring junction. The pressure rise method, also known as the manometric temperature measurement (MTM) method, is an alternative monitoring approach based on the fast closing of the valve separating the vacuum chamber from the condenser, for a short time period on the order of 5-20 s. Due to the closure of the valve the pressure inside the chamber increases, for the accumulation of the water vapor sublimated from the product, at first rapidly and then more slowly when the chamber pressure

10.1021/ie7017433 CCC: $40.75  2008 American Chemical Society Published on Web 10/09/2008

8446 Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008

pressure rise tests; in particular, the new algorithm is able to provide better information on the state of the product than analogous models proposed in the literature, giving consistent results throughout about all the primary drying phase. The final goal is the development of a tool that can be used for in-line monitoring and/or control of the primary drying phase. State of the Art

Figure 1. Examples of experimental pressure rise tests run at different times (measured from vacuum) during primary drying. Upper: 5% w/w BSA solution (pc ) 10 Pa, NvA/Vc ) 1.4, Tshelf ) 268 K, type A vials). Bottom: 10% w/w sucrose solution (pc ) 5 Pa, NvA/Vc ) 0.46, Tshelf ) 258 K, type A vials).

approaches the equilibrium value with the ice surface which, in its turn, slowly increases during the test. Some examples of the evolution of the chamber pressure during the pressure rise test are shown in Figure 1 for two different products freezedried at different shelf temperatures and chamber pressures. Chamber pressure data can be collected during the pressure rise test and related to the sublimating interface temperature by means of mathematical models. If the analysis of the pressure rise curve is performed at different times during primary drying, the evolution of the product temperature can be monitored all along the process. Furthermore, the pressure rise gives also information on the entity of the sublimation flow rate, which is directly related to the slope of the pressure profile at the beginning of the pressure rise run; in the figures shown it is evident that the sublimation rate initially increases and then decreases toward the end of the primary drying. Thus, the method can potentially provide some information concerning the passage from the primary drying phase to secondary drying, characterized by the strong decrease of the vapor flux flowing from the vials into the chamber. Finally, being a model-based estimation method, the MTM approach can supply additional information on the parameters of the process, such as the heat and mass transfer resistances;7-10 of course, the reliability of these estimated parameters is strongly dependent on the accuracy of the model adopted. The pressure rise method could be also used during the secondary drying phase in order to determine the drying rate and eventually, by integration, the residual moisture content in the product, provided that the actual moisture content at the end of primary drying is known; of course, in this case the mechanisms involved are different, and the shape of the curve is almost linear. However, the vapor flow that characterizes the secondary drying stage is very low, thus making the application of this method more difficult, even if feasible in principle. Most of the proposed MTM methods provide good results at the beginning of the drying phase, while poorer predictions are given in the last part of the process.7,11 The objective of the study presented in this paper is to develop and validate an alternative dynamic mathematical model obtained considering the elementary physical mechanisms involved during the

Early works12,13 investigated the transient pressure response measured during a pressure rise test by means of graphical methods. Oetjen proposed and patented a method, called BTM, according to which the interface temperature can be estimated at the pressure at which the first derivative of the pressure rise curve reaches the maximum.14,15 More recent published works propose different mathematical models in order to compute the temperature of products undergoing lyophilization through pressure rise data.6,8,16,17 Milton et al.8 who proposed the name “manometric temperature measurement” for this technique, modeled the transient pressure response under the assumption that four mechanisms contribute to the pressure rise: (1) direct sublimation of ice through the dried product layer at a constant temperature (2) increase in the ice temperature due to continuous heating of the frozen matrix during the measurement (3) increase in the temperature at the sublimation interface when a stationary temperature profile is obtained in the frozen layer (4) leaks in the chamber The values of the vapor pressure over ice, of the product resistance, and of the heat transfer coefficient at the vial bottom, are determined with a regression analysis, by fitting the measured pressure rise response to the pressure values calculated through a model consisting of the sum of the previous terms. The thermal gradient across the ice, the temperature at the vial bottom, and the thickness of the frozen product must be known to perform the regression analysis. Finally, the interface temperature is calculated from the vapor pressure over ice using thermodynamic vapor-ice equilibrium data. Several approximations are made in developing the model, which can be potential sources of errors: the four contributions are considered purely additive; the thermal gradient across the frozen layer is assumed constant and the interface temperature is considered equal to the mean product temperature; finally, the values of the thickness and of the thermal gradient are needed but they are not exactly known. The overall heat transfer coefficient is estimated, but according to the authors its value is not very reliable. The same approach was proposed successively with some small modifications by Pikal et al.18 and Tang et al.19 A modification of the previous model was proposed by Obert.17 It considers two additional physical phenomena: the desorption of the bound water during the primary drying and the thermal inertia of the glass wall of the vial. The overall heat transfer coefficient is expressed by adopting the heat and mass transfer steady-state hypothesis, while a nonlinear regression analysis is carried out in order to estimate the vapor pressure at the interface, the mass transfer resistance in the dried product, and the desorption rate. Then, the interface temperature is found through the thermodynamic equilibrium data of water vapor-ice. Again, the temperature at the bottom of the vial and the thickness of the frozen layer should be known in order to implement the model, but they are only guessed in the proposed procedure. In a more rigorous approach, Chouvenc6 proposed the pressure rise analysis (PRA) method, where the arbitrary

Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008 8447

hypothesis of additivity of the four mechanisms contributing to the total pressure increase is abandoned and a new model is set up, based on elementary heat and mass balances; a simple macroscopic heat balance for the frozen product is adopted, assuming that the temperature increase at the interface is the same as the mean product temperature rise, and the thermal capacity of the portion of the vial glass in contact with the frozen product is included. The sublimating interface temperature at the beginning of the pressure rise and the mass transfer resistance of the porous layer are estimated by regression methods. The thickness of the frozen layer is a needed parameter of the model. The authors propose to estimate it by considering constant sublimation flow between two subsequent PRA tests. In this way, the value of the thickness at a generic run is retrieved by subtracting from the total initial mass of product the sum of the mass sublimated up to the time of the current PRA run. The application of the PRA method to various systems gave interesting results, confirming that the contribution of desorption was generally negligible. Some of the simplifications introduced in developing the previous model anyway are questionable. The assumption that the temperature increase at the interface is the same as the mean product temperature rise implies a constant temperature gradient along the ice during the pressure rise. This can be reasonable toward the end of the primary drying, where the thickness of the ice is small, but not in principle in the first part of the drying cycle, when the frozen layer thickness is larger, and accumulation effects prevent the temperature gradient from being constant. For the same motivations the assumption of constant temperature at the vial bottom during the pressure rise appears not to be realistic, since the vial bottom is continuously heated during the process. Another questionable point concerns the determination of the frozen layer thickness, mainly because it requires a relatively high frequency of PRA tests to give a good accuracy. The assumption of constant sublimation flow between two subsequent pressure rise tests introduces an error in the determination of the ice thickness which is summed to the error accumulated up to the previous test. Thus, the quality of the estimation is expected to decrease toward the end of the primary drying phase. A more complex mathematical model was developed by Liapis and Sadikoglu;16 it implements a dynamic pressure rise method as a remote sensing procedure for determining the following quantities at different times during the primary drying stage: (1) temperature of the moving interface between the dried layer and the frozen layer of the product (2) temperature close to the upper surface of the dried layer of the product (3) temperature of the bottom surface of the frozen layer of the product (4) temperature profile of the frozen layer of the product Furthermore, the authors present an expression suitable for calculating the position of the moving front by knowing the temperature of the heating shelf. Few approximations are introduced in the proposed model with respect to the unsteadystate mathematical model presented by Sadikoglu and Liapis,20 from which it is derived and that was found to describe satisfactorily the experimental dynamic behavior of the primary and secondary drying stages of bulk solution freeze-drying of pharmaceuticals in trays. The model is able to produce the values of the interface temperature and position, as well as the temperature profiles in the frozen layer, starting from pressure rise data. Nevertheless,

many parameters are requested to perform the analysis: diffusivity and permeability of the porous layer; the vial bottom heat transfer coefficient; the temperature and the partial pressure at the top of the vial. The quantity and the type of the parameters and of the measurements required to implement the model makes its practical in-line application a complicated task, even if feasible in theory; as a matter of fact, this model was validated only through simulations. As a last comment, it can be evidenced that all the previous works consider an almost instantaneous closing of the valve, a condition that is generally verified except in very large apparatus; in any case the model can be adapted to the case of slow valve closing.2,21 The Dynamic Parameters Estimation Method (DPE) A detailed one-dimensional unsteady-state model for heat transfer in the frozen layer is proposed in order to take into account the different dynamics of the temperature at the interface and at the vial bottom. The resulting partial differential equation (eq 1) describes conduction and accumulation in the frozen layer during the pressure rise test. Of course, the use of a one-dimensional model does not allow taking into account radial effects, and in particular the curvature of the sublimating interface, but considering the results by Pikal et al.,22 Sheehan and Liapis,23 and Bru¨lls and Rasmuson,24 it seems adequate for the purpose of this work, which aims for monitoring and control. In fact, the use of a very sophisticated model that takes into account all the phenomena involved in the process would require the solution of a large set of equations involving many parameters whose values are poorly or not known at all, and would be computationally very expensive; it must be remembered that DPE must work on-line, as a part of an automatic control loop, and since a best-fitting procedure is required, the model must be sufficiently simple to allow fast and robust computations. The initial condition for the temperature (eq 2) is written under the assumption that the system is in pseudostationary conditions during primary drying operation, i.e., before each pressure rise test is run. Considering initial pseudostationary conditions corresponds to assuming a linear temperature profile in the frozen layer at the beginning of the test. Concerning the boundary conditions, the heat flux at the bottom of the vial (eq 4) is determined by the energy coming from the heating plate (and includes the contributions due to radiation, contact, and conduction in the gas layer in the vial bottom gap), while at the interface it is assumed to be equal to the heat flux due to sublimation (eq 3); as the energy balance refers to the frozen layer, radiation from the upper shelf is not taken into account, as it affects the dried layer, and may be relevant only at the early beginning, when anyway it is largely limited by the shielding effect of the stopper, or in case a significant fraction of heat is transferred to the product through this mechanism. A reference system solidal with the moving interface is considered. In this approach, both radiation from the vial side and conduction along the vial sidewall are not explicitly considered, but it has been already proven that the effect of the vial wall, both for the parallel heat conduction and for the radiating contribution from the chamber wall, can be accounted for in a one-dimensional model by using an effective heat transfer coefficient Kv.25 It must be evidenced that only small radiating heat fluxes from the sidewall can be accounted for (as shown also by Hottot et al.7), and their effect will be to slightly increase the effective heat transfer contribution with respect to the value of the same heat transfer coefficient that

8448 Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008 Table 1. Values of the Variables and Parameters Used in Pressure Rise Simulations for a Generic Pharmaceutical Product parameter cp,frozen Fleak kfrozen k1 Kv L Mw pc,0 pin,0 R Vc/NvA Ffrozen ∆Hs Tshelf

value

unit

1930 0.018 2.56 2.9 × 10-3 19.6 8 × 10-3 18 10.0 0.0 8314 1.434 1030 2687.4 × 103 233.15 + 11.67t; for 0 3 exp[-20.9470( 273.156 / T - 1) 3.56654 ln( 273.156 / T ) + 2.0189(1 - 273.156 / T ) - 5.0983] Pa

pw,i

-1

J kg K-1 Pa s-1 W m-1 K-1 m2 s-1 W m-2 K-1 m kg kmol-1 Pa Pa J kmol-1 K-1 m kg m-3 J kg-1 K K Pa

would be adopted in a more detailed model that includes the radial dimension (or at least the wall); DPE, as all the other MTM approaches proposed up to now, will fail if the main heat source is by radiation from the top or side walls, because this will modify significantly the temperature profile in the product by a mechanism that is not included in the model itself. It has been shown, by means of numerical simulations that will be discussed in the following, that the contribution of the vial wall to the dynamics of the system during the pressure rise test (PRT) is negligible; for these reasons, the thermal inertia of the glass vial has not been included in the balance. Thus, heat transfer in the frozen layer is described by the following equations: kfrozen ∂2T ∂T ) ∂t Ffrozencp,frozen ∂z2 T|t)0 ) Ti,0 +

kfrozen

z kfrozen

|

∆Hs

for t > t0, 0 e z e Lfrozen (1)

( )(

)

k1Mw pw,i,0 - pw,c,0 RTi,0 L - Lfrozen for 0 e z e Lfrozen (2)

( )(

k1Mw pw,i - pw,c ∂T ) ∆Hs ∂z z)0 RTi L - Lfrozen

kfrozen

|

∂T ) Kv(Tshelf - TB) ∂z z)Lfrozen

)

[

for t g t0

Lfrozen Tshelf - Ti,0 k1Mw pw,i,0 - pw,c,0 kfrozen ∆Hs RTi,0 L - Lfrozen

( )(

)

Nv

m ˙ w,tot )

(4)

]

∑ m˙

where the subscript j refers to the jth vial. Applying the ideal gas law and rewriting the mass flow rates as functions of the pressure driving force between the interface and the chamber, it follows, considering instantaneous closing of the valve:

( )

N

[

v MwVc dpw,c k1,jMw pw,i(Ti,j) - pw,c ) Aj RTc dt RTi,j L - Lfrozen,j j)1

∑

(5)

By this way, Kv is a function of the initial value of the interface temperature, Ti,0, of the effective diffusivity of water vapor in the porous matrix, k1, and of the thickness of the frozen product,

]

(7)

where Aj is the actual surface of sublimation of the jth vial. Equation 7 can be simplified introducing a correction coefficient γ which accounts for the heterogeneity of the vial batch, whose expression is given by Nv

γ)

∑ j)1

[

Aj

] )

k1,jMw pw,i(Ti,j) - pw,c RTi,j L - Lfrozen,j

( )(

k1Mw pw,i(Ti) - pw,c NvA RTi L - Lfrozen

(8)

Of course γ ) 1 in the case of a uniform batch. In eq 8 A is the cross section area of the vial (that is, the geometric surface of the product, independently of the real surface from which sublimation occurs). Substituting the expression for γ in eq 7, we obtain

( )

( )(

MwVc dpw,c k1Mw pw,i(Ti) - pw,c ) γNvA RTc dt RTi L - Lfrozen

)

for t g t0 (9)

Finally, the total pressure is computed through eqs 10 and 11, taking into account a constant leakage in the chamber: pc ) pw + pin ) pw + Fleakt + pin,0 pw|t)0 ) pc,0 - pin,0

-1

(6)

w,j

j)1

for t g t0 (3)

where T ) T(z,t), Ti ) T(t)|z)0, and TB ) T(t)|z)Lfrozen. Thermodynamic equilibrium is assumed at the sublimating front, corresponding to the axial position z ) 0; thus pw,i ) pw,i(Ti). The interpolation function used in the simulations is given in Table 1. The heat fluxes at z ) 0 and at z ) Lfrozen, corresponding to the internal bottom of the vial, are generally not equal during the pressure rise run, because of accumulation in the frozen layer, except at the beginning because of the pseudostationary hypothesis. Thanks to this assumption, the expression for the heat transfer coefficient, Kv, can be derived by equating the boundary conditions (3) and (4), both taken at t ) t0. Thus Kv )

Lfrozen. It can be noted that Kv actually is dependent on the pressure over the shelf,22,24 which varies during the PRT; nonetheless, a constant value for Kv can be used, since this dependence is not relevant during the PRT, due to the thermal inertia of the system, as has been evidenced by Chouvenc21 and confirmed by our simulations. Vapor flow in the dried product takes place in the Knudsen regime, and neglecting the pressure gradient of the inert gas, it has been modeled in a simple way in eqs 2 and 3 as dependent on the effective mass diffusivity and on the dried cake thickness, and proportional to the difference of pressure between the sublimating interface and the chamber. In this approach the resistance of the stopper and additional resistances over the shelf are neglected as well. However, these resistances can be easily taken into account in the previous equations, and alternatively a global mass transfer resistance can be adopted, as discussed in the following (see eqs 16 and 17). The previous equations are completed with the equation providing the dynamics of the water vapor pressure rise in the chamber, which consists of the material balance for the vapor flowing into the chamber environment. The total vapor flow rate is given by the sum of the flow rates from each vial:

for t g t0

for t ) t0

(10) (11)

If the value of the chamber temperature Tc is not available, it can be substituted with the product temperature at the interface, usually committing a small error. Of course, the actual thickness of the frozen layer is needed to perform the calculation. Rather than calculating Lfrozen a posteriori as in Chouvenc,6 in our model the ice thickness is

Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008 8449

determined through a material balance written across the moving interface, which is solved contemporaneously to the previous equations. For a single vial the water vapor flow rate at the interface is equal to the difference between the rate of disappearance of the frozen mass and the rate of formation of the dried mass, according to the following equation: m ˙ w,j ) FfrozenAj

dLfrozen,j dLfrozen,j - FdriedAj dt dt

(12)

After the introduction of the correction coefficient γ, the material balance at the interface can be integrated in time between the previous pressure rise test and the actual one, obtaining (-1) Lfrozen ) Lfrozen

Mw R∆F

∫

t0

t0(-1)

(

γ

)

k1 pw,i - pw,c dt Ti L - Lfrozen

(13)

where ∆F ) Ffrozen - Fdried, and the superscript “(-1)” refers to quantities calculated or measured in the previous pressure rise run. A simplified mathematical model for primary drying could be employed in order to solve eq 13, such as the model proposed by Velardi and Barresi25 which comprises two differential equations giving the time evolution of Lfrozen and Ti. These equations can be integrated from t (-1) to t0 with initial 0 conditions known (since determined after the application of the DPE method at time t 0(-1)), obtaining the value of Lfrozen to be used in the current DPE test. Indeed, as the application of this tool is for in-line monitoring and control, and it is thus preferable to reduce the computational demand, it seemed preferable to use a quadrature formula to carry out the integration; here the integral on the right-hand side was solved applying the trapezoidal rule of integration: (-1) Lfrozen ) Lfrozen

[( (

)

Mw k1 pw,i,0 - pw,c,0 γ + R∆F Ti,0 L - Lfrozen γ(-1)

)]

(-1) (-1) k(-1) pw,i,0 - pw,c,0 t0 - t (-1) 1 0 (14) (-1) (-1) 2 T L-L i,0

frozen

The use of the composite form of the trapezoidal rule is not possible in this case since the function to be integrated according to eq 13 is only known in the extremes of the interval t 0(-1) and t0 (i.e., at the times when pressure rise tests are performed). At the first pressure rise test, when no data are available from previous runs, eq 13 can be integrated from the time when vacuum was made in the chamber, where the interface temperature can be assumed equal to the shelf temperature and the frozen layer thickness equal to the initial thickness of the product. It may be worth reminding that, once the Lfrozen value has been estimated, eq 2 immediately gives the initial temperature profile in the frozen layer, including the value at the bottom. The spatial domain of the frozen layer was discretized in order to transform the partial differential equation (eq 1) in a DAE (differential algebraic equations) system; the orthogonal collocation method was employed to obtain the values of T(z,t) in the nodes of the spatial grid. For further details on the discretization method, refer to Velardi and Barresi.25 At each pressure rise test the discretized system of equations is integrated in time; the integration is performed in the internal loop of a curvilinear regression analysis, where the cost function to minimize in a least-squares sense is the difference between the simulated values of the chamber pressure and the actual values measured during the pressure rise. In the simplest approach, the parameters to be estimated are the initial interface temperature Ti,0 and the vapor effective diffusivity inside the dried

Figure 2. Steps of the optimization procedure in the DPE algorithm.

product, k1, while the correction factor γ can be set equal to 1. In order to obtain a better fitting of the experimental data, γ can be estimated too, and in this case, of course, it is variable with process time. The Levenberg-Marquardt method was used in order to perform the minimization of the cost function. The steps of the optimization procedure, sketched in Figure 2, are the following: (1) initial guess of Ti,0 and k1 (and γ) (2) determination of Lfrozen from eq 14 using the values measured and computed in the pressure rise test run at time t ) t 0(-1) (3) determination of Kv from eq 5 (4) determination of the initial temperature profile in the frozen mass, from eq 2 (5) integration of the discretized DAE system in the time interval (t0,tf), where tf - t0 is the time length of the pressure rise test (6) iteration of points 1-5 and evaluation of Ti,0 and k1 (and γ) that best fit the simulated chamber pressure pc to the measured data, pc,meas, in order to solve the nonlinear-least-squares problem: min

Ti,0,k1,γ

1 1 ||p - pc,meas||22 ) min Ti,0,k1,γ 2 2 c

∑ [p - (p c

c,meas)k

]2

(15)

k

Multiple solutions due to local minima cannot be excluded a priori, and a good initial guess of the parameters at point 1 is thus required. This is usually not a problem, since the values of Ti,0, k1, and γ determined in the previous DPE run can be used, as the time elapsed between two pressure rise tests is relatively short, compared to the time scale of the whole process. A multistep method based on the backward differentiation formulas was initially employed to integrate in time the discretized system of equations as the first version of the numerical code for the DPE algorithm was developed using the software MATLAB, to take advantage of the versatility of MATLAB’s high-level programming language. As a matter of fact, this method gave satisfactory results, with computing times compatible with the regression procedure of the DPE algorithm. Anyway, a FORTRAN version of the code is being set up for the in-line application of DPE, as part of a control system, where

8450 Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008 Table 2. Geometric Characteristics of the Vials and Stoppers Used in the Experimental Tests parameter type volume internal diameter wall thickness bottom thickness

vial A tubing 4 14.25 1 0.7

vial B

tubing 12 20.85 1.2 1.0 (min) 1.5 (max) contact gap (maximum) 0.4 1.1 stopper type four parallel vents igloo stopper opening 3 (trapezoidal, largest side) 8 (width)

unit mL mm mm mm mm mm

more efficient routines will be used: in particular, modern routines based on a multivalue algorithm where the polynomial that interpolates numerically the solution is computed in one point through its successive derivatives by means of the Nordsieck vector are particularly interesting, as this approach allows modifying at any step both the order of the method and the step size without reinitializing the procedure as would happen with a multistep method; it will be necessary to verify if they really will be more efficient also in the case of our DAE system, where the multivalue algorithm could be ill-conditioned. Equations 6-9 have been developed for the case of product in vials, but the method works in principle also for bulk product distributed in a tray, provided that a sufficient uniformity of the product thickness and of the heat transfer conditions is obtained. Also in this case, it is possible to introduce a corrective coefficient that accounts for nonuniformity; the surface must be ideally divided into Nv zones, and the subscript j in eqs 6-8 and eq 12 refers to each of these zones. Results and Discussion The applicability and reliability of the proposed dynamic parameters estimation method was evaluated through both simulations and experimental lyophilization cycles. In the first case, the detailed model for primary drying described and validated in Velardi and Barresi25 was used as a source of simulated experimental pressure rise data, in the typical conditions adopted during the freeze-drying of pharmaceuticals in vials; the scope of validating the DPE algorithm using simulated data as a first step, similarly to what done by Liapis and Sadikoglu,16 was the validation of the proposed approach and of the used algorithm, with particular attention to the problem of regression and possible multiple solutions. In fact, by this way it is possible to verify that the behavior of the product inside the vial, which cannot be measured experimentally with the required accuracy, is estimated correctly. In particular, the reconstruction of the introduced data makes it possible to evaluate the accuracy of the prediction of the heat and mass transfer coefficients that are affected by large uncertainty when measured experimentally, and also to verify the prediction of nonmeasurable parameters and quantities. The Velardi and Barresi25 model can supply reliable simulations of the real system, as it is based on unsteadystate heat transfer balances in the dried layer, in the frozen layer, and along the vial glass sidewall. Sublimation at the moving interface and mass transfer of a binary mixture of water vapor and inert gas through the porous matrix are modeled with the Dusty Gas Model approach. Thus, the equations of the model allow describing the time and spatial evolution of the product temperature and of the mass fluxes inside the dried material, as well as the progression of the moving front of sublimation; it can be underlined that the model was experimentally validated using the smallest vials employed also in the present work. In

Figure 3. Upper: Simulated temperature profiles at the bottom and at the interface of a generic pharmaceutical product during primary drying. A pressure rise test is carried out at each hour during the drying time. Lower: Comparison of interface temperature values predicted by the detailed model and estimated by DPE algorithm. Conditions of Table 1, type A vials. (s) Interface product temperature (Ti); (- -) bottom product temperature (TB); (b) Ti values estimated by DPE in correspondence with PRT (20 s).

order to simulate in detail the behavior of a vial during a PRT, the influence of the variation of the pressure during the test on the sublimation rate had to be taken into account; this was done by coupling the model of the vials to the model of the chamber. The detailed model was also used during the development of the DPE algorithm, in order to assess the relevance of some simplification, as the assumption of a constant Kv during the PRT, or the role of the thermal inertia of the glass vial and of heat conduction along the vial wall, and to quantify the uncertainty deriving for example from the presence of radiating heat transfer from the chamber walls. Indeed, the dynamic model used for the DPE derives from the detailed model, but as it has to work efficiently in-line, it has been conveniently simplified; in particular, it does not take into account explicitly the glass wall. A typical lyophilization cycle for pharmaceuticals was simulated using the data given in Table 1 and type A vials (see Table 2), considering a stopper with negligible resistance; the simulated primary drying time was found to be 10.26 h. Ten pressure rise experiments were run, one every hour, simulating the closure of the chamber-to-condenser valve for a time period of 20 s and with a pressure acquisition frequency of 10 Hz. The pressure evolution in the chamber during the pressure rise test was simply obtained by a mass balance to the “closed” system, as the total mass flowing out from the vials was set equal to the mass accumulation in the chamber. Figure 3(upper) shows the temperature at the bottom and at the interface during the main drying phase. The peaks observable in the temperature profiles are due to the pressure rise tests and are more evident toward the end of the drying. The last peak is very pronounced because the pressure rise virtually coincided with the end of the drying phase; here in fact, the heat provided at the bottom mainly contributes to increase the product temperature due to the very small mass of still frozen product. This phenomenon had been already evidenced, and as a consequence these methods should be used with care in the last part of primary drying, in order to avoid damages to the product: a possibility is to properly reduce the measuring time, as the temperature increase is proportional to the duration of the pressure rise test.

Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008 8451

Figure 4. Simulated pressure rise test run after 4 h of primary drying. Conditions of Table 1. Upper: Pressure rise dynamics in the chamber and at the moving interface. (b) pw,i actual values; (O) pc actual values; (s) pw,i estimated by DPE; (- -) pc fitted by DPE. Lower: Temperature profiles in the frozen product. (b) actual T(z) values at t ) 0 s; (O) actual T(z) values at t ) 4 s; (s) T(z) estimated by DPE at t ) 0 s; (- -) T(z) estimated by DPE at t ) 4 s.

Figure 4 shows as an example the results of the analysis for the pressure rise experiment run after 4 h of primary drying. Four seconds of data acquisition was sufficient in this case to perform the optimization analysis. This behavior is in agreement with the experimental results of Milton,8 where it is also observed that just a few seconds of pressure rise may provide sufficient data for the analysis of the response in the first phase of primary drying. In the figures the actual values predicted by the detailed model are compared with those estimated from the DPE algorithm by means of the fitting and optimization procedure, to confirm that the correct values are recovered and no artifact or compensation between the estimated parameters occurs; it must be remarked that it is not possible to measure directly the interface temperature, and it would be extremely difficult to measure the complete temperature profile of the solid, at least without heavily disturbing the system. Very good agreement is observed for the profiles of the pressure in the chamber and at the interface (Figure 4, upper graph). In Figure 4, lower graph, the ice temperature profiles along the frozen product are given at the beginning and at the end of the test; it can be seen that the proposed DPE model is able to predict fairly well the temperature increase in the ice, at any axial position, with only a very small discrepancy between the estimated and the true temperature values at the bottom of the product. The interface temperature Ti,0 was estimated to be 241.36 K, practically coincident with the value provided by the detailed model. The values of Ti,0 obtained by the regression analysis of the pressure rise tests throughout the primary drying step are given in Figure 3 (lower plot) and compared with the actual (simulated) ones. For the sake of clarity, the temperature peaks recorded in correspondence with the pressure rise tests are not shown in this case. The agreement is very satisfactory,

Figure 5. Upper: (b) Estimated values of the effective diffusivity k1; (0) estimated values of the effective heat transfer coefficient Kv; (- -) average j v. Lower: Estimated mass transfer effective heat transfer coefficient K resistance of the dried cake.

and a very small deviation is observed only at the end of primary drying, where the optimal PRT time is shorter. The estimated effective diffusivity in the dried layer, k1, is shown in Figure 5 together with the heat transfer coefficient Kv obtained from eq 5. Comparing the values of the two parameters estimated at different times with the corresponding values reported in Table 1 (a constant value of these parameters was assumed, equal to 2.9 × 10-3 m2 s-1 and 19.6 W m-2 K-1, respectively), very good agreement is found concerning k1 while Kv, as expected, is slightly overestimated. This can be explained if we consider that part of the energy provided by the heating shelf during drying is transferred to the product through the vial sidewall, as evidenced by Velardi and Barresi.25 Since the DPE model does not account for heat transfer along the vial sidewall, differently from the detailed model, slightly larger Kv values are predicted by DPE in order to compensate this additional heat input, as discussed in the previous section; in any case, as shown in Figure 5, this value corresponds to the j v ) 21.2 W m-2 K-1, which average effective coefficient, K would be necessary to consider adopting a simplified onedimensional model in order to obtain the same results provided j v has been calculated through eq 16 by the detailed one. K following the procedure presented by Velardi and Barresi:25

[

j v ) average K

T˜B(t) - T˜i(t) kfrozen T˜shelf(t) - T˜B(t) L˜frozen(t)

]

(16)

where the symbol “∼” refers to quantities obtained through the detailed model simulation. From the values of Ti,0, k1, and Lfrozen provided by the DPE method, it is possible to calculate the mass transfer resistance of the dried layer, RP, as follows: k1Mw 1 1 ) RP RTi L - Lfrozen

(17)

8452 Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008

The calculated mass transfer resistance RP exhibits a quasi-linear dependence upon the thickness of the dried cake, as shown in Figure 5 in the lower graph, as expected in the case of uniform cake structure, that does not change (as a consequence of collapse or microcollapse) with time. In the case considered, the resistance of the stopper was neglected; in any case, it is very easy to take it into account including the term RS in eqs 2 and 3, defining a global resistance:

(

RTi L - Lfrozen 1 1 ) ) + RS RPS RP + RS Mw k1

)

-1

(18)

Simulations with the detailed model have also evidenced that during the short time of a pressure rise test the increase in the overall heat transfer coefficient, consequent to the increase in the total pressure that favors the transfer from the shelf to the vial bottom, has a very limited effect on the interface temperature; in fact, the assumption of pseudosteady state is no longer valid during the test, and the temperature gradients are modified in the frozen product, as heat is accumulated preferentially in the lowest part of the product. The same simulations evidenced that it is very difficult to treat in a simple way the role of the wall, because depending on the time considered during the primary drying, and the axial position, the wall can be at higher or lower temperature than the product; therefore, it does not appear correct to simply include a portion of the vial glass in the energy balance, apart from the difficulty of estimating the right value of this fraction, that in fact seems variable with the process time.21 In any case the same simulations have evidenced that accumulation in the vial wall has a negligible effect during the PRT; normal conduction in the sidewall, as shown in the previous figure, is properly taken into account through an effective Kv coefficient. Using the detailed model and the same procedure described before (see eq 16), it was also confirmed that even in the presence of limited radiation from sidewalls the DPE is able to correctly estimate the interface temperature; under this condition a reasonable effective heat transfer coefficient is estimated, even if with less accuracy than in cases where no radiation is present. As the analysis of simulated cases evidenced that the developed monitoring system can be considered reliable, a series of experimental cycles was carried out in order to evaluate the performance of the DPE method applied to real cases. The experimental testing was run at the Politecnico di Torino in a pilot unit (LyoBeta-25 by Telstar, Terrassa, Spain) with a chamber volume of 0.20 m3 (0.178 m3 excluding the volume of the shelves). The freeze-drier was equipped for these tests with a set of thermocouples, an absolute capacitance manometer (MKS Instruments Baratron 626A, Wilmington, MA), a thermal conductivity pressure gauge (Inficom Pirani, Balzers, Liechtenstein), and a special analytical balance developed by Politecnico di Torino26,27 able to weigh up to 15 vials in real time during the process. The apparatus can work with either controlling the pressure with a valve on the vacuum pump or with a controlled leakage valve, and the DPE has been tested in both cases. The runs shown in the paper have been carried out with the second system, stopping the inert flow in correspondence with the PRT. The experimental runs were carried out using several products characterized by different glass transition temperatures: bovine serum albumin (BSA) (supplied by Sigma Aldrich), sucrose, and mannitol (not shown) (Ph. Eur. Analytical specification, supplied by Riedel-de Hae¨n). Different operating conditions were tested, with the number of vials ranging from about 40 to 700 and the chamber pressure varying from 5 to 20 Pa. Two

different sizes were used for the tubing vials; a stopper with higher resistance, having four parallel trapezoidal vents, and an igloo-type one were employed for the smaller and larger vials, respectively. The geometric characteristics are given in Table 2; the height of frozen product was 7.1 mm for type A vials and 9.9 mm for type B vials, corresponding to 1 and 3 mL of filling solution, respectively. The diameters considered are in the range of values generally used for freeze-drying processes; in particular, they correspond approximately to some of those considered by Pikal22 and are respectively slightly larger and slightly smaller than those used by Chouvenc21 for validation of the similar PRA method. The vials were arranged on a tray, partially stoppered and surrounded by empty vials in order to minimize the effect of radiation from the walls and of thermal conduction from the edges of the tray. Trays of different sizes were used, depending on the size of the lot. The product temperature of selected vials was monitored by placing a thin copper-constantan thermocouple in close contact with the internal part of the bottom of the vial. The liquid product was initially frozen down to -45 °C, with a cooling rate of approximately 1 °C/min; then vacuum was made in the chamber and the shelf temperature was raised to the selected temperature. Once the set point pressure selected for the primary drying phase was reached, pressure rise tests were carried out every 20 or 30 min, depending on the runs. During each pressure rise test the butterfly valve between the chamber and the condenser, after an almost instantaneous closing (τ < 0.5 s), was kept closed for less than 30 s. During this time interval the chamber pressure data were acquired at a sampling rate of 50 Hz using the Baratron capacitive manometer (accuracy 0.25% of reading, full scale range 100 Pa). Then, the collected samples were analyzed by the DPE solver. Figure 6 shows an example of comparison between the pressure rise data acquired during a test (carried out after 6 h from the beginning of the primary drying phase, run at 5 Pa with a batch of 609 vials) and the pressure values estimated by the DPE. As reported in the previous part, the DPE algorithm is based on a nonlinear regression method which determines the best fitting between the experimental pressure rise data and the simulated ones. Thus, the closer to the experimental data is the curve obtained with the DPE solver, the better is the quality of the regression analysis. As can be observed, when the two-parameter (Ti,0, k1) regression method is employed, the fitting between the experimental and the estimated pressure data is good, especially in the first seconds of the pressure rise. When the three-parameter (Ti,0, k1, γ) regression analysis is performed, the fitting quality is even better, enabling a very good approximation of all the pressure rise curve. Indeed, this could be an expected behavior, since on increasing from 2 to 3 the number of degrees of freedom of a curve-fitting problem the goodness of fit improves as well. Nonetheless, the parameter γ was introduced as a term accounting for the heterogeneity of the system, playing an important role especially in small-scale apparatus; thus its adoption has a physical foundation rather than being a numerical trick. In any case, the DPE algorithm performances were generally good even with the two-parameter fitting, so this option was selected for the following cases. In the bottom plot of Figure 6 the increase of the bottom temperature TB during the pressure rise test is shown. The vial is continuously heated from the bottom while the heat removal at the interface is reduced in the last part of the test due to the increased chamber pressure, which reduces the driving force for sublimation. This leads to an increase in the product

Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008 8453

Figure 6. Freeze-drying of 10% w/w sucrose solution (pc ) 5 Pa, NvA/Vc ) 0.46, Tshelf ) 258 K, type A vials). Upper: Comparison between measured and estimated pressure rise curve. (O) pc,meas experimental values; (- -) pc estimated by DPE with two-parameter regression; (s) pc from DPE with three-parameter regression. Lower: product temperature at the bottom. (s) Estimated with two-parameter regression; (- -) estimated with threeparameter regression; (b) temperature measured by thermocouple.

temperature during the DPE test, as already evidenced from the simulations, and the longer is the DPE run, the larger is the temperature increase. Considering the product temperature variation during the test is very important, as the pressure rise time should be long enough to allow the acquisition of a sufficient number of pressure samples, but as short as possible in order to avoid potentially damaging temperature overshoots. The filled circle at t ) 0 is the average value of the temperatures measured by the thermocouples placed in contact with the bottom of the monitored vials. No data are available giving the experimental temperature profile during the pressure rise test, as the maximum acquisition rate allowed by the apparatus was limited to one sample per minute. Nevertheless, the estimated temperature at the beginning of the test (in the case of both two- and three-parameter regressions) is in good agreement with the experimental data, within a standard deviation interval, and the difference between the two predictions, even if amplified by the scale of the plot, is always very small, below the accuracy of the thermocouples generally employed. Since the direct measurement of the interface temperature is not feasible in practice, some comparisons are proposed between the values of the temperature TB estimated by the DPE solver and the product temperatures monitored along the primary drying by the thermocouples inserted at the bottom of sampled vials. Figure 7 (upper graph) shows the results for a run with BSA solution, where six thermocouples inserted in different vials were used to monitor the product temperature; the average values and the standard deviation bars (only in correspondence with the DPE estimation) are plotted. The agreement is very good, with the temperatures being well estimated within a standard deviation interval, except at the very end where, the behavior of the batch is less uniform as evidenced by the increased variance in the values of the temperature measured

Figure 7. Freeze-drying of 5% w/w BSA solution (pc ) 10 Pa, NvA/Vc ) 1.4, Tshelf ) 268 K, type A vials). Upper: (s) Mean temperature at the vial bottom measured by thermocouples and error bars giving standard deviation; (b) TB estimated by DPE; (0) sublimation flux estimated by DPE. Lower: (b) Estimated values of the effective diffusivity k1; (0) estimated values of the effective heat transfer coefficient Kv. Time is given from beginning of primary drying.

Figure 8. Comparison between experimental and estimated temperatures at the bottom of the product for a cycle run at 5 Pa with 43 vials (type A) (10% w/w sucrose solution). ( · · · ) TB measured by thermocouples placed in four different vials; (- -) heating/refrigerant fluid temperature; (b) TB estimated by DPE; (s) Pirani to Baratron pressure ratio. The start of the primary drying, corresponding to pressure reduction in the chamber, is indicated by the arrow.

by the thermocouples. The problem of batch uniformity will be discussed later in more detail. Similarly, the comparison given in Figure 8 for a sucrose solution evidences good agreement between the experimental and estimated temperature data almost up to the last part of the primary drying step. Indeed, as discussed in the Introduction, it is well-known that the thermocouples are not very reliable when used as a monitoring technique for the determination of the end of the primary drying phase, as a consequence of the uncertainty in the position of the thermocouple inside the vial and of the variability exhibited by the thermocouple-monitored

8454 Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008

Figure 9. Comparison between the ice thickness determined from mass measurements through eq 19 (s) and estimated by the DPE solver (O) for a cycle run at 5 Pa with 43 vials (type A), 15 of which are weighed by the balance (10% w/w sucrose solution).

samples. The data shown in Figure 8 confirm this, evidencing that, relying on the thermocouple response, the duration of primary drying can be largely underpredicted. Thus, in order to confirm that the DPE tool is effective almost up to the end of the primary drying, its duration was evaluated independently monitoring the ratio between the pressure measured by the Pirani gauge and the Baratron sensor. This signal can be used to detect the end of the ice sublimation step, following the suggestions by Nail,28 and previous work29 has shown that this indication is reliable and consistent with that of other devices commonly employed. The special weighing device operating in-line in the freezedrier chamber was thus employed in order to show the capability of the DPE method to estimate easily the position of the moving interface between the still frozen and the already dried product, and thus to monitor the progress of primary drying. In fact, as evidenced by previous authors, the MTM methods are generally quite robust in estimating the interface temperature, but its position is often predicted with lower accuracy. The mass of the vials measured by the balance throughout the cycle can be converted in ice thickness according to the following linear relation, assuming a constant porosity of the dried cake throughout the sublimation step: Lfrozen(t) )

[

]

m(t) - m(tfinal) L m(tinitial) - m(tfinal)

(19)

where Lfrozen(t) is the ice thickness at time t, L is the initial thickness of the product, m(t) is the mass weighed by the balance, tinitial is the time at the beginning of the sublimation phase, and tfinal is the end time of primary drying; the residual weight at the end of the primary drying can be calculated a priori, possibly considering the residual moisture in the dried product, if known, for better accuracy. The ice thickness values calculated with eq 19 were compared with those computed by the DPE method. The results reported in Figure 9 show that the agreement between the experimental and estimated data is fairly good. The obtained results suggest that the DPE method provides consistent results all over the primary drying phase. Nevertheless, in correspondence with the end of the sublimation phase, identified by the balance and by the pressure sensors, the estimated temperature exhibits a strong decrease, as outlined already in Barresi.30 A similar behavior has been reported also by Oetjen14 for the BTM approach and by Tang et al.11 using MTM. A possible explanation for the phenomenon is that not all the vials contribute to the pressure rise in the final part of the process as a consequence of the nonuniformity of the batch.

Figure 10. Freeze-drying of 10% w/w sucrose solution: example of DPE estimations obtained in two cycles run with vials of different dimensions (see Table 2). Left: Type A vials, placed on a medium-size rectangular tray and not shielded (pc ) 10 Pa, NvA/Vc ) 0.14, Tshelf ) 263 K; total primary drying time 16 h 35 min). Right: Type B vials, placed on a smaller circular tray and shielded by empty vials (pc ) 10 Pa, NvA/Vc ) 0.082, Tshelf ) 253 K; total primary drying time 17 h 17 min). (b) Moving front temperature; (2) global mass transfer resistance; (9) ice thickness estimated by the DPE solver.

In fact, at the end of the sublimation phase, mostly vials positioned in the center zone of the tray contribute to the pressure increase, because they are protected from the edge effect that affects the vials in peripheral positions; these dry faster, and close to the end of the step may have no ice remaining. Thus, the reduction of the area of sublimation is compensated by the DPE model with a drop in the estimated product temperature. Besides, the prediction of the transfer coefficients becomes less reliable because of the phenomenon described above. Figure 7 (lower graph) shows the values of k1 estimated along the primary drying and the values of Kv calculated through eq 5. As can be seen, DPE provides almost constant values in the central part and up to the final part, when they decrease in correspondence with the decrease of the sublimation rate (Figure 7, upper graph), estimated by the initial slope of the pressure rise curves. At this time, in fact, the batch has partially lost its uniformity because some vials have completed sublimation, as can be also inferred by the temperature increase recorded by some thermocouples that indicates the loss of thermal contact between the sensor and the ice. In Figure 7 k1 has been calculated taking into account the contribution of the stopper to the total mass transfer resistance between the sublimating interface and the chamber, according to eq 18; the value has been estimated experimentally as described below. Figure 10 shows the evolution of Ti, RPS, and Lfrozen/L during the primary drying for two cycles run with two types of vials of different diameters (and with different stopper shapes). At the beginning of the sublimation phase, when Lfrozen ) L, the global resistance is given by the stopper; thus an estimation of RS can be obtained by the value of RPS in correspondence with time zero. As expected, type B vials show a significantly lower RS, because the single large opening of the stopper used in this case causes a lower pressure drop than the four smaller channels present in the stopper used for type A vials. RPS increases with the drying progress, at first steeply in the first part of the drying phase, during which the temperature of the product is rising, and then almost linearly and with a lower slope in the zone

Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008 8455

leads to the estimation of an effective Kv higher than the real Kv evaluated for a two-dimensional system, as evidenced in Figure 5; the smaller the diameter of the vial, the percentage larger the contribution of the wall to heat transfer, and thus the increase in the effective heat transfer coefficient. Conclusions

Figure 11. Overall maximum heat transfer coefficient estimated by DPE at different chamber pressures. (b) Type A vials on medium-size rectangular tray; (0) type B vials on a smaller circular tray. Freeze-drying of 10% w/w sucrose solution, shielded vials.

where the interface moves at an almost constant velocity. It can be noted that, notwithstanding the lower resistance of the stopper, it consistently assumes higher values in the right-hand plot (Figure 10), as in this case the thickness of the product is higher. The strong increase at the beginning of primary drying is probably mainly due to the finer structure that is generally observed close to the upper surface, and to the formation of a thin crust with low permeability. At the end of the primary drying step the global resistance steeply increases again, but this is an expected behavior related to the poorer estimation of k1 in the very last phase of the process. In Figure 10 the estimation of the position of the frozen interface is also shown. The behavior observed for the interface temperature and RPS (and similarly for k1 and Kv, not shown), may indicate that the heterogeneity in this case is larger for the small vials batch. The overall heat transfer coefficients between the heating fluid and the vial bottom, calculated for some different cases, are shown in Figure 11; the maximum estimated Kv is plotted against the chamber pressure for two different configurations (type A and B vials on different trays). It can be observed that DPE provides Kv values in good agreement with the findings reported by Hottot et al.7 and consistent with those calculated using the correlations of Pikal et al.22 These correlations, based on the theory of conduction in low-pressure gases, evidence the dependence of the overall heat transfer coefficient on the total pressure; this behavior is correctly reproduced by the DPE method. It must be taken into account that the presence of the tray approximately halves the heat transfer coefficient; this is confirmed by our results, as for example at the pressure of 10 Pa Kv reduces from 16.7 to 9.8 W m-2 K-1 for the type B vials as a consequence of the presence of the tray. The difference between the vials of different size, at the same pressure, is mainly due to the different geometries of the vials used and possibly to the different configuration of the tray, as in varying the size of the tray (and eventually its planarity) and its surface finishing the heat transfer coefficient may change. In fact, the maximum gap at the bottom, and thus the contribution of conduction in the gas layer, which depends on chamber pressure, is different in the two cases: type A vials not only have a smaller maximum gap, but also a different profile, with a very small gap in the middle, which further reduces the resistance. Furthermore, it must be remembered, as discussed before, that the mathematical model used in the DPE algorithm does not account explicitly for the energy exchanged with the vial sidewall and by radiation from the environment, which finally

The dynamic parameters estimation method (DPE) is presented in this paper as an advanced noninvasive technique for in-line monitoring the primary drying step of the lyophilization process. The DPE is based on the well-known pressure rise analysis concept, but evidences some improvements respect to analogous methods proposed in the literature. The developed mathematical model, based on unsteady-state heat and mass balances, describes in detail the physical phenomena occurring during the pressure rise test. Thus, differently from previous methods, the DPE is able to provide the frozen product temperature at any axial position (from the vial bottom to the interface) at any time during the pressure rise test. This aspect is quite important in process monitoring, since the product must not be freeze-dried above a maximum temperature value that should not be overcome even during the short time of the test in order to avoid the damage of the product. The DPE method was validated both with modeling simulations and with experimental cycles. Data from simulations were used to evaluate the applicability of the method to typical operating conditions used for the freeze-drying of pharmaceuticals, showing its capability to estimate well different process parameters. Also, mathematical simulations were usefully employed to assess the impact of the simplifications introduced in the DPE algorithm with the aim of reducing the computational effort and thus the processing time for in-line applications. The results obtained from the experimental campaign evidenced that the DPE method allows computing consistent temperature values almost up to the end of the primary drying stage. Different products and chamber pressures were tested, obtaining reliable estimations of the transfer coefficients even in the case of relatively small vials, where the role of the vial wall may be more relevant; in these cases it must be remembered that the effective coefficients estimated by DPE may differ from those evaluated for a two-dimensional system. The thickness of the still frozen product, taken as an indication of the progression of the sublimation process, is well estimated too if compared with the data obtained through a special inline weighing device. An extensive test of the method including also large-scale industrial apparatus and a larger selection of products and experimental conditions is planned; the scope is to evidence strengths and weaknesses and to compare it with the other methods previously proposed. Some points, common to practically all methods based on the manometric temperature measurement approach, remain open and deserve further investigation. One of them is the selection of the optimal acquisition time; it has been shown that in some cases a few seconds suffice, and that prolonged times expose the product to higher temperatures, but the minimum time required to allow a correct estimation may vary, depending on chamber volume, shelf loading, operating conditions, product characteristics, and also changes during the same cycle. The minimum acquisition time can be evaluated a posteriori, but hardly estimated a priori, unless a predictive model, such as the one used in this work for validation, is adopted. An acquisition time of 20-30 s is generally a good compromise for reasonable loading, but it must be considered that too

8456 Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008

prolonged times with respect to the minimum required can reduce the accuracy of the estimation, even if to a small extent; thus a check of the results as a function of the acquisition time considered may be recommended. The other important aspect is the nonuniformity of the batch; it was shown that the introduction of a correction factor, variable with the process time, for this effect is theoretically grounded, and it can improve the quality of the prediction, even if further developments are required to allow the estimation of the batch variance. In fact, during the drying, the effect of the differences in ice structure, radiation and heat flux from shelf, just to mention the principal sources of batch variability, can vary in an unpredictable way the sublimation rate of the single vials, but this affects not only the value of the γ factor defined in eq 9, but also the average value of the parameters estimated. As a matter of fact, the DPE (like the other MTM methods) finally calculates an “effective” average value of the interface temperature and of the other parameters. Preliminary results shown in Rasetto et al.31 evidence that radiation and fluid dynamics inside the chamber can affect the results of DPE, and the relative contribution of these phenomena varies according to the equipment size; these different behaviors must be taken into account when transferring the results from small- to large-scale units. When the primary drying is at the end and a fraction of the vials is already dried, the reduction of the sublimating interface is the main cause that affects the predictions of these methods, and in that case the correction factor can be interpreted with good approximation as the percentage of vials with residual frozen product. At the current stage of development, DPE does not give a clear estimation of the end point of the main drying, and in any case the recourse to the pressure rise test at the end of the sublimation stage could lead to product damaging, because of excessive temperature increase. Thus DPE (and generally MTM methods) can be coupled with another sensor to establish the end of primary drying;29 anyway, some outputs of this method, and in particular the extrapolation of the frozen layer thickness or the relative variation of the sublimation rate, that can be accurately evaluated also in the secondary drying, can be considered to this purpose. The final objective is thus to use the DPE as a reliable tool for monitoring effectively the process, especially in combination with other devices (e.g., Baratron and Pirani sensors, weighing device) that can give complementary information. Furthermore, DPE can be effectively employed in combination with a simplified model of the primary drying in order to realize a novel control system for automatically finding the best cycle recipe.32 In this approach, DPE provides at different times during the process the estimation of the product temperature and of the other key parameters (transport coefficients, product thickness, etc.). These process variables are fed to the control system, based on an advanced predictive control logic that determines in real time the most efficient heating strategy in order to minimize the drying time with respect to the quality constraints on the product. A preliminary application of this control tool is shown by Galan et al.33 and Barresi et al.29

di Torino) for experimental data and for detailed simulations respectively is gratefully acknowledged. Notation A ) cross surface of the vial, m2 cp ) specific heat at constant pressure, J kg-1 K-1 Fleak ) leakage rate, Pa s-1 ∆Hs ) heat of sublimation, J kg-1 k ) thermal conductivity, J m-1 s-1 K-1 k1 ) effective diffusivity coefficient, m2 s-1 Kv ) overall heat transfer coefficient, J m-2 s-1 K-1 L ) total product thickness, m Lfrozen ) frozen layer thickness, m m ) mass, kg m ˙ ) mass flow rate, kg s-1 M ) molecular weight, kg kmol-1 Nv ) number of vials p ) pressure, Pa rsub ) sublimation flux, kg m-2 s-1 R ) ideal gas constant, J kmol-1 K-1 RP ) mass transfer resistance in the dried layer, m s-1 RS ) mass transfer resistance in the stopper, m s-1 RPS ) global mass transfer resistance, m s-1 t ) time, s T ) temperature, K Tshelf ) temperature of the heating/cooling fluid in the shelf, K V ) volume, m3 z ) axial coordinate, m Greek Symbols γ ) correction coefficient F ) mass density, kg m-3 Subscripts and Superscripts 0 ) value at t ) 0 (-1) ) pressure rise test before the actual one B ) vial bottom, corresponding to z ) Lfrozen c ) chamber dried ) dried product f ) end of pressure rise run frozen ) frozen product i ) sublimating interface, corresponding to z ) 0 in ) inert gas j ) generic vial k ) generic chamber pressure sample acquired during a pressure rise test meas ) measured tot ) total w ) water vapor AbbreViations BSA ) bovine serum albumin DPE ) dynamic parameters estimation MTM ) manometric temperature measurement PRA ) pressure rise analysis PRT ) pressure rise test

Acknowledgment The authors would like to acknowledge E.U., for financial support in the framework of the research project LYO-PRO: Optimization and control of the freeze-drying of pharmaceutical proteins (GROWTH Project GRD1-2001-40259-RTD). Contributions from Roberto Pisano and Davide Fissore (Politecnico

Literature Cited (1) Jennings, T. A. Lyophilization: introduction and basic principles; Interpharm/CRC Press: Boca Raton, 1999. (2) Oetjen, G. W.; Haseley, P. Freeze-Drying; Wiley-VCH: New York, 2004.

Ind. Eng. Chem. Res., Vol. 47, No. 21, 2008 8457 (3) Rey, L.; May, J. C. Freeze-Drying/Lyophilization of pharmaceutical and biological products, 2nd ed.; Marcel Dekker: New York, Basel, 2004. (4) Fissore, D.; Velardi, S. A.; Barresi, A. A. In-line control of a freezedrying process in vial. Drying Technol. 2008, 26 (6), 685–694. (5) Roy, M. L.; Pikal, M. J. Process control in freeze-drying: determination of the end point of sublimation drying by an electronic moisture sensor. J. Parenter. Sci. Technol. 1989, 43, 60–66. (6) Chouvenc, P.; Vessot, S.; Andrieu, J.; Vacus, P. Optimization of the freeze-drying cycle: a new model for pressure rise analysis. Drying Technol. 2004, 22 (7), 1577–1601. (7) Hottot, A.; Vessot, S.; Andrieu, J. Determination of mass and heat transfer parameters during freeze-drying cycles of pharmaceutical products. PDA J. Pharm. Sci. Technol. 2005, 59 (2), 138–153. (8) Milton, N.; Pikal, M. J.; Roy, M. L.; Nail, S. L. Evaluation of manometric temperature measurement as a method of monitoring product temperature during lyophilization. PDA J. Pharm. Sci. Technol. 1997, 51 (1), 7–16. (9) Rambhatla, S.; Ramot, R.; Bhugra, C.; Pikal, M. J. Heat and mass transfer scale-up issues during freeze-drying: II. Control and characterization of the degree of supercooling. AAPS Pharm. Sci. Technol. 2004, 5 (4), article 58, E1-E9 (http://www.aapspharmscitech.org). (10) Tang, X. C.; Nail, S. L.; Pikal, M. J. Evaluation of manometric temperature measurement (MTM), a Process Analytical Technology tool in freeze-drying, part III: heat and mass transfer measurement. AAPS Pharm. Sci. Technol. 2006, 7 (4), article 97, E1-E7 (http://www.aapspharmscitech. org). (11) Tang, X. C.; Nail, S. L.; Pikal, M. J. Freeze-Drying Process Design by Manometric Temperature Measurement: Design of a Smart Freeze-Dryer. Pharm. Res. 2005, 22, 685–700. (12) Neumann, K. H. Determining temperature of ice. In Freeze-drying of Foods and Biologicals; Noyes, R., Ed.; Noyes Development Corp.: Park Ridge, NJ, 1968. (13) Willemer, H. Measurement of temperature, ice evaporation rates and residual moisture contents in freeze-drying. DeV. Biol. Stand. 1991, 74, 123–136. (14) Oetjen, G. W. Freeze-Drying; Wiley-VCH: New York, 1999. (15) Oetjen, G. W.; Haseley, P.; Klutsch, H.; Leineweber, M. Method for controlling a Freeze-Drying process. U.S. Patent, 6,163,979, 2000. (16) Liapis, A. I.; Sadikoglu, H. Dynamic pressure rise in the drying chamber as a remote sensing method for monitoring the temperature of the product during the primary drying stage of freeze-drying. Drying Technol. 1998, 16 (6), 1153–1171. (17) Obert, J. P. Mode´lisation, optimisation et suivi en ligne du proce´de´ de lyophilisation. Application a` l’ame´lioration de la productivite´ et de la qualite´ de bacteries lactiques lyophilise´es. Ph.D. Thesis, INRA ParisGrignon, 2001. (18) Pikal, M. J.; Tang, X.; Nail, S. L. Automated process control using manometric temperature measurement. U.S. Patent US 6,971,187 B1, 2005. (19) Tang, X. C.; Nail, S. L.; Pikal, M. J. Evaluation of manometric temperature measurement, a Process Analytical Technology tool for freezedrying: part I, product temperature measurement. AAPS Pharm. Sci. Technol 2006, 7 (1), article 14, E1-E9 (http://www.aapspharmscitech.org). (20) Sadikoglu, H.; Liapis, A. I. Mathematical modeling of the primary and secondary stages of bulk solution freeze-drying in trays: parameter

estimation and model discrimination by comparison of theoretical results with experimental data. Drying Technol. 1997, 13, 43–72. (21) Chouvenc, P.; Vessot, S.; Andrieu, J.; Vacus, P. Optimization of the freeze-drying cycle: adaptation of the Pressure Rise Analysis to noninstantaneous isolation valves. PDA J. Pharm. Sci. Technol. 2005, 5, 298– 309. (22) Pikal, M. J.; Roy, M. L.; Shah, S. Mass and heat transfer in vial freeze-drying of pharmaceuticals: role of the vial. J. Pharm. Sci. 1984, 73 (9), 1224–1237. (23) Sheehan, P.; Liapis, A. I. Modeling of the primary and secondary drying stages of the freeze-drying of pharmaceuticals products in vials: numerical results obtained from the solution of a dynamic and spatially multi-dimensional lyophilisation model for different operational policies. Biotechnol. Bioeng. 1998, 60 (6), 712–728. (24) Bru¨lls, M.; Rasmuson, A. Heat transfer in vial lyophilization. Int. J. Pharm. 2002, 246, 1–16. (25) Velardi, S. A.; Barresi, A. A. Development of simplified models for the freeze-drying process and investigation of the optimal operating conditions. Chem. Eng. Res. Des. 2008, 86, 9–22. (26) Vallan, A.; Parvis, M.; Barresi, A. A. Sistema per la misurazione in tempo reale di massa e temperatura di sostanze sottoposte a liofilizzazione. Italian Patent Application B02005A000320, 2005. (27) Vallan, A. A measurement system for lyophilization process monitoring. Proceedings of Instrumentation and Measurement Technology ConferencesIMTC, Warsaw, Poland, 13 May 2007; IEEE: Piscataway, NJ, USA, 2007, 5 pp. DOI: 10.1109/IMTC.2007.379000. (28) Armstrong, J. G. Use of the capacitance manometer gauge in vacuum freeze-drying. J. Parenter. Drug Assoc. 1980, 34, 473–483. (29) Barresi, A. A.; Pisano, R.; Fissore, D.; Rasetto, V.; Velardi, S.; Vallan, A.; Parvis, M.; Galan, M. Monitoring of the primary drying of a lyophilization process in vials. Chem. Eng. Process 2008, DOI: 10.1016/ j.cep.2008.05.004, in press. (30) Barresi, A. A.; Pisano, R.; Rasetto, V.; Velardi, S.; Vallan, A.; Parvis, M.; Galan, M. Monitoring, control and optimisation of freeze-drying process. Proceedings of the European Drying Conference AFSIA 2007, Biarritz, France, 24-25 May 2007; Cahier de l’AFSIA No. 22, pp 78-79. (31) Rasetto, V.; Marchisio, D. L.; Fissore, D.; Barresi, A. A. Modelbased monitoring of a non-uniform batch in a freeze-drying process. In 18th European Symposium on Computer Aided Process Engineering; Braunschweig, B.; Joulia, X., Eds.; Computer-Aided Chemical Engineering Series, Vol. 25; Elsevier: Amsterdam, 2008, paper FP_00210, 6 pp; CD Edition. (32) Velardi, S.; Barresi, A. Method and system for controlling a freezedrying process. European Patent application PCT /EP2007/059921 (19/ 09/2007), 2007. [Priority EP200619587 (19/09/2006), European Patent EP1903291, published 26/03/2008]. (33) Galan, M.; Velardi, S. A.; Pisano R.; Rasetto, V.; Barresi, A. A. A gentle PAT approach to in-line control of the lyophilization process. New Ventures in Freeze-Drying (IIR & Aerial). Strasbourg, France, 7-9 NoVember 2007; Refrigeration Science and Technology Proceedings; Institut International du Froid: Paris, 2007, 17 pp; CD-ROM Edition.

ReceiVed for reView December 21, 2007 ReVised manuscript receiVed August 1, 2008 Accepted August 7, 2008 IE7017433