Modeling the Continuous Freeze-Drying Process for Making Biopharmaceuticals

The majority of lyophilization processes in the (bio)pharmaceutical industry have been operated in batch mode for many decades. The conventional configuration typically comprises a large number of vials located on the cooling/heating shelf (Figure 1A). The bottom shelf, whose temperature can be manipulated, is used to cool the vials during the freezing step and heat the vials during the drying steps. In this case, the freezing and drying processes take place in the same space but at different times.

Several drawbacks associated with batch lyophilization are extensively discussed in Ref. [6, 7]. Typical issues of any batch process include batch‐to‐batch variability, quality control, process downtime, and operational flexibility. Another disadvantage that is specific to the configuration illustrated in Figure 1A is non‐uniformity in heat transfer. For example, vials located on the side and corner of the shelf are typically affected by thermal radiation from the chamber walls, whereas the effect is almost negligible for center vials. Variability in heat transfer conditions could result in different crystal structures and drying times, hence variation in the final product quality. Consequently, process control and optimization are not straightforward. These known issues motivate the development of a continuous lyophilization concept that can be practically implemented in an industrial scale.

A comprehensive review of continuous lyophilization technologies is given in Ref. [6]. To summarize, the current continuous lyophilization technologies can be divided into four main categories. The first category employs equipment that continuously moves a bulk product through the process, e.g., a conveyor^[29^] or a revolving plate with holes.^[30^] Another type of technology relies on the concept of spray freeze drying, in which a liquid product is sprayed through nozzles to create small liquid droplets with a high surface‐to‐volume ratio, greatly improving heat, and mass transfer in the process.^[8, 28^] Both aforementioned approaches have been proposed for lyophilization of bulk material. For lyophilization of unit doses, a well‐known technique uses spin freezing to accelerate the freezing process and improve the uniformity of heat transfer in a product.^[5^]

Finally, the current state‐of‐the‐art continuous lyophilization technology for unit doses has been recently proposed by Ref. [7]. This technology employs the suspended‐vial configuration as illustrated in Figure 1B, where all vials are suspended and continuously move through the freezing and drying chambers, without any contact between the vials and shelf. Unlike conventional lyophilizers, the freezing and drying processes occur in different locations/chambers of the equipment. This suspended‐vial technology has several benefits over the existing ones. First, heat transfer in the equipment is uniform because every vial experiences the exact same heat transfer conditions along the way. In addition, there is a dedicated chamber for controlling the ice nucleation process to help reduce variation in the crystal structure caused by the stochastic nature of ice nucleation. Second, the process does not involve any complicated motions (e.g., high speed flow in spray freeze drying), resulting in a simpler design and more reliable quality control. Third, there is no contact between the vials and cooling/heating shelf, minimizing the risk of having fine particles that could lead to contamination, which is crucial in (bio)pharmaceutical applications. Finally, with the automatic filling and load‐lock systems, the process is fully continuous.

The mechanistic model developed in this work is based on this suspended‐vial configuration as illustrated in Figure 1B, the current state‐of‐the‐art continuous lyophilization technology for unit doses. Detailed modeling strategies and model equations are discussed in the next section (Section 3), while information and discussion related to the equipment design can be found in Ref. [7].

Mechanistic modeling for conventional/batch lyophilization has been studied for many decades, with various models that have been proposed and used in both academia and industry. Despite the fact that the models are for batch processes, modeling strategies underpinning those models can be used as a basis for the development of a continuous lyophilization model for suspended vials. While some past reviews briefly discuss the literature on mechanistic modeling of lyophilization,^[4, 19^] no article systematically summarizes the available modeling strategies and comprehensively discusses pros and cons of each model variation. Hence, this section first provides a systematic overview of the mechanistic modeling strategies for lyophilization, discusses pros and cons of each model variation in detail, and finally concludes with the optimal modeling strategies used in this work. The published models for continuous lyophilization mentioned in Section 1 are based on the spin and spray freeze drying, which are completely different from the suspended‐vial configuration, and so those modeling strategies are not discussed here.

Figure 2 summarizes various modeling strategies for the key steps of lyophilization, namely freezing, primary drying, and secondary drying. Published models mostly focus on primary drying because this step is recognized as the most time‐consuming and expensive step, thus the first target for process optimization and improvement. A variety of modeling strategies are available, ranging from the simplest lumped capacity model to the high‐fidelity multidimensional model. Models for secondary drying have also been well developed because the final product quality (i.e., moisture content) is governed by this step. The mechanistic modeling of freezing receives less attention in the lyophilization literature due to its complicated behavior, including supercooling and stochastic ice nucleation. As a result, the number of available models for freezing is more limited compared to that of the drying models.

For the design proposed by Ref. [7], suspended vials are cooled by using the cryogenic gas at a controlled temperature and flow rate. Besides, a dedicated chamber is added to control the nucleation temperature using the vacuum‐induced surface freezing (VISF) technology.^[38^] With this setup, all surfaces of a vial experience a similar heat transfer condition, which is different from batch lyophilization where only the bottom surface is cooled by the cooling shelf.

Our model (Figure 3A) for freezing of suspended vials is divided into five main steps: 1) preconditioning, 2) VISF, 3) nucleation, 4) solidification or ice formation, and 5) final cooling.

In conventional lyophilization, vials are placed on the heating shelf, and thus heat transfer at the bottom surface is driven by conduction (at the point of contact), convection (gas), and radiation, whereas only radiation dominates at the side and top surfaces.^[18^] For continuous lyophilization, suspended vials are heated by the below heating shelf without any contact between the vials and shelf. As a result, heat transfer at the bottom surface is driven only by thermal radiation and natural convection.

The model for primary drying is formulated in the rectangular coordinate system with one spatial dimension (z) and time (t) (Figure 3B). Define the primary drying step to start at t₀ and complete at t_d1. If the primary drying model is simulated consecutively after the freezing step, then set t₀ = t_f4. Otherwise, t₀ should be set to 0 for a standalone primary drying simulation.

The governing equations for the primary drying step consist of the Equation (1) energy balance in the frozen region and Equation (2) mass balance at the sublimation front. By assuming that the supplied heat is used in the frozen region only, the energy balance for the frozen region can be described by the partial differential equation (PDE) (32)ρfCp,f∂T∂t=kf∂2T∂z2+QradVf,S<z<Hwhere T(z, t) is the temperature, S(t) is the sublimation front/interface position, k is the thermal conductivity, ρ is the density, C_p is the heat capacity, V is the volume, H is the height of the product, and the subscript 'f' denotes the frozen region. We refer to Appendix B for the calculations of some relevant parameters. The radiative heat transfer from the sidewall Q_rad is(33)Qrad=σArFs3Tc4−T4where T_c is the chamber wall temperature and A_r = πdH is the side area of the product. We refer to Section 3.5 for the detailed derivation of Equation (33). The mass balance of water at the sublimation front gives(34)dSdt=Nwρf−ρewhere N_w is the sublimation flux and ρ_e is the effective density of the dried region above the sublimation front. The driving force for mass transfer at the sublimation interface is^[1, 18, 22^](35)Nw=pw,sat−pw,cRpwhere p_w,sat is the saturation/equilibrium pressure of water, p_w,c is the partial pressure of water in the chamber (environment), and R_p is the mass transfer resistance. The saturation pressure for sublimation is described by Ref. [22](36)pw,sat=exp−6139.9T+28.8912The variation of the mass transfer resistance can be approximated by the empirical expression^[4, 18, 22^](37)Rp=Rp0+Rp1SRp2+Swhere R_p0, R_p1, and R_p2 are the constants to be estimated from data. Physically, R_p0 represents the resistance associated with mass convection above the product, which could include the presence of a skin layer on the cake surface and the effect of the stopper,^[44^] whereas Rp1SRp2+S corresponds to the resistance associated with diffusion through the porous dried layer in the product.

Modeling the heat transfer and sublimation front in one dimension, as described in Equations () and (), implies an assumption that the sublimation front is flat. This assumption has been widely used for conventional lyophilization as the heat flux through the bottom of the vial is much stronger than that from the sidewall, mainly due to heat conduction. In lyophilization of suspended vials, due to the absence of conductive heat transfer at the bottom surface, the sidewall heat flux becomes more comparable to the bottom heat flux. However, the bottom heat flux still dominates for two key reasons. First, the bottom shelf remains the primary heat source in the system, which directly interacts with the bottom surface of the vial through thermal radiation and convection. Second, since most vials are surrounded by others at similar temperatures, the driving force for sidewall heat transfer is significantly weaker than that for bottom heat transfer. A systematic analysis for comparing different heat transfer modes can be done using relevant dimensionless numbers, e.g., in Appendix . 32 34 A

The PDE represented by Equation (32) requires two boundary conditions. Heat transfer at the bottom surface of the frozen product follows Newton's law of cooling, (38)−kf∂T∂z=hb(T−Tb),z=Hwhere T_b is the bottom shelf temperature and h_b is the overall heat transfer coefficient that combines the effects of thermal radiation and natural convection. At the top surface, the energy balance associated with sublimation and thermal radiation from the upper surface of the chamber is (39)NwΔHsub=kf∂T∂z+σFs1Tu4−T4,z=Swhere ΔH_sub is the heat of sublimation.

The initial conditions for Equations (32) and (34) are (40)T(z,t0)=T0,0≤z≤H(41)S(t0)=0For consecutive simulation with the freezing step, set T₀ = T(t_f4). Otherwise, T₀ can be set arbitrarily for a standalone primary drying simulation.

The primary drying model is simulated until the interface position is equal to the height of the product, i.e., S = H, indicating that there is no frozen material left, which marks the end of the primary drying step at t_d1. In some cases, there might be an additional heating period at the end of primary drying to adjust the temperature and ensure complete sublimation before starting the secondary drying step. However, the model contains only a simple heat equation, and so it is not detailed here.

The model for secondary drying is formulated in the rectangular coordinate system with one spatial dimension (z) and time (t) (Figure 3C), which is consistent with the primary drying model. The secondary drying step is defined to start at t₀ and complete at t_d2. If the secondary drying model is simulated consecutively after the primary drying step, then set t₀ = t_d2. Otherwise, t₀ should be set to 0 for a standalone secondary drying simulation.

The governing equations for the secondary drying step comprise the 1) energy balance in the dried region and 2) desorption kinetics. The energy balance of the dried product is (42)ρeCp,e∂T∂t=ke∂2T∂z2+ρdΔHdes∂cw∂t+QradVe,0≤z≤Hwhere T(z, t) is the product temperature, c_w(z, t) is the concentration of bound water (aka moisture content, residual moisture, residual water), ρ_d is the density of the dried region (solid and vacuum), ΔH_des is the heat of desorption, Q_rad is as defined in Equation (33), and the other parameters are as defined in Equation (32), with the subscript 'e' denoting the effective properties considering both solid and gas in the pores. The desorption kinetics of bound water is described by (43)∂cw∂t=kd(cw∗−cw)where cw∗ is the equilibrium concentration of bound water and k_d is the rate constant for desorption that exhibits Arrhenius temperature dependence^[4, 15, 16^](44)kd=fae−Ea/RTwhere f_a is the frequency factor (aka collision frequency), E_a is the activation energy, and R is the gas constant. The above desorption kinetics, Equation (43), is known as the linear driving force model, one of the simplest adsorption/desorption models that can accurately predict the dynamics of bound water and has been widely used in the literature.^[4, 15, 19, 23, 45^] Further simplification that is relatively common can be done by setting cw∗=0. This simplification produces insignificant error as shown in Ref. [16] and eliminates the need for equilibrium data and detailed knowledge about the solid structure.^[4^]

The governing PDE, Equation (42), requires two boundary conditions. The bottom surface of the dried product is heated by the heating shelf, which follows Newton's law of cooling (45)−ke∂T∂z=hb(T−Tb),z=Hwhere the value of h_b is approximated to have the same value as that used in primary drying, Equation (38). Heat transfer at the top surface is mainly thermal radiation, resulting in the boundary condition (46)−ke∂T∂z=σFs1Tu4−T4,z=0

The initial conditions for Equations (42) and (43) are (47)T(z,t0)=T0,0≤z≤H(48)cw(z,t0)=cw,0,0≤z≤HFor consecutive simulation with the primary drying step, set T₀ = T(z,t_d1). Otherwise, T₀ can be set arbitrarily for a standalone secondary drying simulation.

The secondary drying model should be simulated until the concentration of bound water is below the target value, denoted as c_w,∞, which marks the end of the secondary drying step at t_d2.

Convection and radiation are the two important modes of heat transfer in suspended, vials as discussed extensively in the previous sections. During the freezing step, there could be a combination of natural/forced convection and radiation, depending on how the cryogenic gas is fed to the system. In primary and secondary drying, only natural convection and radiation are important.

Convection can be simply modeled using Newton's law of cooling. Heat transfer coefficients for convection can be estimated from correlations that entail dimensionless groups such as the Nusselt number, Prandtl number, Reynolds number (forced convection), and Grasholf number (natural convection). Nevertheless, it is more common and accurate to estimate the heat transfer coefficient from data, e.g., using the techniques suggested in Ref. [4].

Thermal radiation is significantly more complicated in terms of modeling. The theories and model equations described below are mainly based on in Ref. [31] and; Ref. [39] the former discusses general thermal radiation theories, while the latter provides a detailed study on the modeling of thermal radiation in lyophilization. Here we summarize only the key elements needed for modeling thermal radiation in the suspended‐vial configuration; more detailed discussion can be found in the aforementioned references. Consider radiation exchange between two diffuse, gray surfaces of finite size denoted as surfaces 1 and 2, respectively. In general, the net radiant energy receiving by surface 1 can be written as (49)Qrad=σA1FT24−T14where A is the surface area, F is the transfer factor, and the subscripts '1' and '2' denote surfaces 1 and 2, respectively. The transfer factor is dependent on the geometry and material properties (emissivity) of both surfaces. For any two diffuse, gray surfaces that form an enclosure, (50)Qrad=σT24−T141−ε1ε1A1+1A1F1−2+1−ε2ε2A2where F₁₋₂ is the view/shape factor. If surface 1 is surrounded by surface 2 (i.e., F₁₋₂ = 1) and surface 2 is much larger than surface 1, Equation (50) can be simplified as(51)Qrad=ε1σA1T24−T14In Equation (51), the transfer factor F=ε1. Next, consider the case where a single vial (surface 1) is surrounded by the chamber walls (surface 2). In this case, Equation (51) becomes (52)Qrad=εglσArTc4−T4where ε_gl is the emissivity of the glass vial. Here the transfer factor F=εgl. The amount of radiative heat for the single‐vial case represented by Equation (52) defines the upper bound on the radiative heat that one vial can receive from the chamber walls. This is because, when there are multiple vials, the view factor F₁₋₂ is less than 1, and so the radiative heat is shared among the vials. As a result, the transfer factor F for the multiple‐vial case must be less than ε_gl. Therefore, Equation (52) needs to be modified to account for interactions between vials.

Typical lyophilization of unit doses always consists of a large number of vials, and so thermal radiation exchange exists not only between vials and chamber walls but also between those multiple vials. In batch lyophilization, vials are arranged as an array consisting of many rows, which results in significant differences in heat transfer conditions between the outer and inner vials. In such cases, the outer vials with higher temperature transfer significant heat to the inner vials while simultaneously exchanging heat with the chamber walls. A rigorous way of modeling this complicated radiation is the radiation network approach, which is comprehensively discussed and demonstrated in the context of lyophilization in Ref. [39]. In the suspended‐vial configuration, however, vials are typically aligned in a few rows, e.g., a single row in the equipment proposed by Ref. [7]. As a result, all vials experience nearly the same heat transfer heat condition, resulting in negligible heat transfer between vials. Therefore, instead of using the radiation network approach, we modify Equation (52) to follow Equation (49) as(53)Qrad=σArFTc4−T4where F must be less than ε_gl as described in the previous paragraph. This transfer factor F is best to be estimated from experimental data but can also be approximated mechanistically, e.g., using Equation (50). In our model, the transfer factor Fs3 appears in Equations (6), (32), and (42). If the chamber design and vial configuration are identical for both freezing and drying steps, the same value of Fs3 can be used for all equations. A similar analysis can also be done for Fs1 in Equations (39) and (46).

When there exists both convection and radiation, it is more convenient to linearize and rewrite the radiation part in the form of Newton's law of cooling; at the temperature of about 300 K, the error caused by linearization is about 0.1% for the temperature difference of 20 K and 2% for the temperature difference of 100 K (see Section 1.3.2 of Ref. [31]), which is tiny. The linearized equation can then be combined with the convection part, with the corresponding heat transfer coefficient taking into account of both convection and radiation, e.g., Equations (8), (38), and (45). This approximation helps simplify the equations and also facilitates the use of an overall heat transfer coefficient (e.g., Equation (29)) as both convection and radiation parts are written in the form of Newton's law of cooling.

Finally, to ensure physically reasonable heat transfer parameters in our model, we note that typical heat transfer coefficients for natural convection in air vary between 3 and 25 W·m⁻²·K⁻¹, and those for forced convection vary between 10 and 200 W·m⁻²·K⁻¹.^[31^]

This section defines the default parameter values for the model developed in Sections 3.2, 3.3, and 3.4. These parameter values are either obtained from the literature or set based on the typical values used in lyophilization. This default set of parameters assumes a complete simulation in which the models for freezing (with controlled nucleation), primary drying, and secondary drying are simulated consecutively. Table 1 lists the default model parameters; parameter values different from those reported in the table are stated explicitly in each section or case study.

The parameters that should be estimated from experimental data include heat transfer coefficients, mass transfer coefficients and cake resistance, and desorption‐related parameters (frequency factor and activation energy). Operating conditions that are not measured, e.g., wall temperatures, could also be estimated from data.

The mechanistic model presented in Section 3 is validated using the experimental data reported in Ref. [7]. Our model validation is carried out for each lyophilization step separately, with the specific parameters reported in Table 2.

For the freezing step, the simulated temperature and freezing time agree well with the experimental data (Figure 5A). The product temperature decreases from its initial value of 280 K to the nucleation temperature of about 263 K, where nucleation starts. Upon the first nucleation, the temperature instantaneously rises up from 263 K to the freezing point of about 272 K. The temperature slowly decreases during the solidification phase for about 30 min, following the freezing‐point depression. Subsequently, the temperature decreases faster and reaches equilibrium with the environment (chamber wall and cryogenic gas), with the simulated temperature approaching the equilibrium slightly faster than the actual temperature.

The model accurately predicts the evolution of the product temperature and drying time during primary drying for both low (263 K) and high (313 K) shelf temperatures (Figure 5B). The temperature rises up rapidly at the beginning and gradually increases at a slower rate toward the end of the process. This is because the driving force for sublimation is initially low, and so most of the heat input is used to increase the temperature. As time progresses, the driving force for sublimation is larger; thus, most of the heat input is used for sublimation. The maximum error in temperature prediction is about 2–3 K, while the error in drying time prediction is about 0.2 h for Case 2a and 0.4 h for Case 2b (≈6%). This error could be partly because of the uncertainty in cake resistance, which results from the stochastic nature of the freezing process. The datasets in Figure 5B, which represent the lower and upper limits of the operating boundaries, demonstrate the validity and generalizability of the model across the full range of operating conditions.

Similar to primary drying, in secondary drying, we consider two different concentration levels for the bound water (residual moisture): 0.088 and 0.075 kg water/kg solid. The simulated concentration of bound water closely aligns with the experimental data for both concentration levels (Figure 5C). The concentration decreases exponentially from its initial value to the final concentration of 0.01 kg water/kg solid, the threshold set for terminating the simulations. This concentration profile is expected given the linear driving force model is used to describe the desorption process. It is important to note that the data used in Figure 5C are based on the shelf temperature of 293 K. To rigorously estimate the values of f_a and E_a, it is recommended to obtain the concentration profiles at various shelf temperatures so that the effects of temperature on the desorption process can be quantified.

Overall, our model is able to accurately predict the time evolution of the critical process parameters for all three steps of lyophilization. The maximum deviation between the predicted and measured temperature is about 3 K, which is about 1% compared to the absolute operating temperature used in lyophilization. This error is smaller than the measurement noise of some non‐invasive temperature sensors, e.g., thermal imaging cameras^[35^] For the concentration of bound water, the maximum error is less than 0.01 kg water/kg solid, which is typically the threshold for identifying the endpoint of secondary drying^[1^] Therefore, the accuracy of our model is sufficient for general process design, optimization, and control.

We note that, in the suspended‐vial lyophilization considered in this article, every single vial moves through the process following the same trajectory as shown in Figure 1B, ensuring identical heat transfer conditions as a function of time for all vials. For batch lyophilization designs where the vials are distributed in different locations (Figure 1A), the evolution of the critical process parameters in each vial also varies. For instance, the outermost vials in a batch lyophilizer dry the fastest due to thermal radiation from the environment, while the thermal radiation effect is much weaker on the inner vials.^[39^]

One of the important aspects of continuous manufacturing is to ensure that the process is operated smoothly, optimally, and safely without human intervention, to minimize downtime and maximize production. In the context of continuous lyophilization, our model can be used to study the evolution of critical process parameters – e.g., temperature, moisture content – throughout the process at various operating conditions, which can help guide process design and optimization.

This work not only presents the first mechanistic model for continuous lyophilization of suspended vials, but also is one of the very few studies that develops a complete model incorporating all three steps of lyophilization, including freezing, primary drying, and secondary drying. This is important specifically for continuous manufacturing where the entire process should be designed and optimized simultaneously.

With the default model parameters in Table 1, a complete continuous lyophilization cycle can be simulated, including 1) freezing with controlled nucleation via VISF, 2) primary drying, and 3) secondary drying (Figure 6). The time evolution of the product temperature (Figure 6A), mass of ice (Figure 6B), and concentration of bound water (Figure 6C) can be obtained, given the operating pressure (Figure 6D) and temperature (Figure 6E). This information could be useful for various purposes. For example, it can guide process and equipment design concerning the chamber sizing and velocity profile of each vial, ensuring that the target residence time is met for each chamber. Simulation results can also be used to help identify pressure and temperature profiles to achieve some specific objectives, e.g., minimization of drying time. It is important to note that the entire simulation shown in Figure 6 can be computed in less than 1 s on a normal laptop, which makes the model very practical to be implemented for any purposes.

As discussed in Section , one of the key considerations during the drying steps is to ensure that the product temperature at any location does not exceed the upper limit. Consequently, spatiotemporal data are needed. Spatial variation of temperature and concentration in the product are dependent on several factors. For instance, temperature gradients are larger when the sample thickness and heat transfer coefficient increase (see the analysis in Appendix ). 3.1 A

To demonstrate, we set the heat transfer coefficient and sample thickness to 30 Wm⁻²·K⁻¹ and 0.02 m,^[17^] respectively. These values are noticeably higher than the default values in Table 1, and so the resulting gradients are more significant. Examples of spatiotemporal data obtained from our model are shown in Figures 7 and 8. During primary drying, the temperature gradient is about 6 K at the beginning and gradually decreases as the sublimation front recedes (Figure 7). A similar behavior can also be observed for the concentration of bound water in secondary drying, except that there is no sublimation front (aka moving interface).

A key benefit of mechanistic models is that they can sometimes provide important insights about the process when the experimental data are limited. This section employs the model and simulation to study the freezing step, which usually receives little attention compared to the drying steps. Spontaneous (aka uncontrolled) ice nucleation was already investigated in detail by Ref. [12, 13]. However, uncontrolled nucleation is not ideal for continuous lyophilization as described in Ref. [7]. Therefore, we focus on controlled nucleation, with the VISF technique as used in Ref. [7]. The input and parameter data are given in Tables 3 and 4.

Several experimental studies on VISF have been conduced in the literature. Here, we use our model and simulation results to provide insights into the process, and use that understanding to explain the previous experimental data. To clearly understand VISF, we first consider cases where that the nucleation temperature is fixed at 260 K; i.e., the effect of stochastic ice nucleation is excluded from the model (Figure 9A,B). In VISF, the key idea is to reduce the system pressure, typically from the atmospheric pressure (10⁵ Pa) to the VISF pressure, to evaporate a small amount of liquid, resulting in a fast decrease in the product temperature and hence nucleation. The product temperature decreases faster with a lower VISF pressure (Figure 9A). For example, when VISF starts at 0.25 h, the temperature drops from the initial value of 268 K to the nucleation temperature of 260 K almost instantaneously for the VISF pressure of 100 Pa, while it takes about 10 min when the VISF pressure is 10⁴ Pa. Consequently, the first nucleation occurs earlier at a lower VISF pressure (dash line in Figure 9B). During VISF, a small amount of liquid evaporates, which increases with a decrease in the VISF pressure (solid line in Figure 9B). This is expected because of a higher mass transfer driving force for evaporation at lower VISF pressure. In general, if VISF is carried out properly, the amount of liquid evaporating is small. In this case, the original liquid mass is about 2.9×10⁻³ kg, and so the mass loss of about 3×10⁻⁵ is almost negligible.

Next, we include the effect of stochastic ice nucleation into our VISF model for a more realistic analysis (Figures 9CD), with a Monte Carlo simulation of 10⁴ runs for proper statistics. Our simulation result show that VISF can significantly reduce the variation in the product temperature, compared to the uncontrolled nucleation case (Figure 9C). A decrease in the VISF pressure also reduces the degree of variation. Besides the product temperature, The variation in the nucleation time (defined as the time when the first nucleation occurs) also decreases with the VISF pressure, with the first nucleation occurring earlier (Figure 9D). These results indicate that the nucleation process is well controlled with VISF (i.e., less variation in the nucleation time and temperature), agreeing with general experimental observation in the literature.

Finally, for model validation and implementation in real‐world applications, the model can be fitted to experimental data obtained from the system of interest for more accurate results; the two important parameters are the mass transfer coefficient and heat transfer coefficient associated with evaporation. For example, by using the data from Ref. [7] the temperature profile simulated by our model agrees well with the experimental data (Figure 9E).

The results presented in this section assume minimal vial‐to‐vial thermal interactions because the vials are arranged in a single row, in which all the vials experience the same heat transfer conditions throughout the process, as described in Section 3.5. This condition is not true in conventional lyophilization where a large number of vials are arranged in a hexagonal or rectangular array. In such cases, vial‐to‐vial interactions could significantly affect the distribution of nucleation times and freezing rates. For example, when one vial nucleates, an increase in the temperature of that vial could slightly heat the adjacent vials, which subsequently delays the first nucleation of those vials. The effects of vial‐to‐vial interactions on freezing have been discussed in Refs. [12, 13, 48]. In the context of modeling, these interactions can be incorporated straightforwardly by coupling our freezing model with the radiation network approach described in Section 3.5 and.Ref. [39]. In any case, such effects become much less important in controlled nucleation.

In summary, this section demonstrates how the model and simulation results elucidate the role of VISF and its operating conditions in promoting controlled nucleation and influencing the nucleation process, improving the uniformity of the product. Furthermore, we show that the model can explain real data relatively well. Hence, our model can be used to understand and guide the design and optimization of the VISF method for controlling the nucleation process in lyophilization.

Typical process modeling focuses on behaviors of the system during its normal operation where the process is operated steadily under the desired conditions. A better process model should be able to simulate important abnormal conditions, e.g., equipment failure, that could occur due to various unexpected scenarios, which is critical in continuous manufacturing as these abnormal operations could affect the reliability, availability, and maintainability (RAM) as well as the safety of the plant/process.

In lyophilization, one of the crucial design considerations is associated with a condenser. During normal operation, the maximum capacity of a condenser must be higher than the rate of vapor production via sublimation, and choked flow should be avoided.^[1^] This operational constraint is to ensure that there will not be vapor accumulation in the equipment, which could subsequently lead to in a pressure increase. This scenario may not lead to safety‐related issues because the operating pressure in lyophilization is low and the amount of ice/water is generally not high enough to create overpressure. Nevertheless, an increase in the total pressure reduces the driving force for sublimation (see Equation (35)), which then prolongs the primary drying step. As a result, the final product quality could be affected. This section explores how to incorporate those dynamic behaviors into our model.

Most primary drying models, including our model, assume that the system pressure is well controlled, which is a reasonable assumption. To incorporate the condenser dynamics into the model, we consider cases where the total condenser capacity is not enough to condense the water vapor produced via sublimation during primary drying. In such cases, the mass balance for the water vapor in a chamber is^[22^](71)dpw,cdt=(jw−jw,max)RT¯VcMwwhere j_w is the total mass flow rate of water vapor resulting from sublimation, j_{w, max} is the maximum condenser capacity, V_c is the chamber volume, and T_avg is the average temperature in the chamber assumed to be constant for simplification. The mass flow rate of water vapor can be calculated by (72)jw=nvialAzNwwhere n_vial is the number of vials in the chamber. By coupling Equations (71) and (72) with the primary drying model developed in Section 3.3, the effects of condenser failure or choked flow can be quantified.

With the specific parameters given in Table 5, the total pressure increases from the normal operating value of 3 to 20 Pa (Figure 10A). This pressure increase results from the total vapor flow exceeding the maximum condenser capacity, thus vapor accumulation in the chamber. When the system pressure increases, the driving force for sublimation (also the sublimation flux) decreases. After about 1 h, the system pressure becomes constant, indicating the equalization of the sublimation flux and condenser capacity. In this abnormal operation, the drying time and product temperature are higher than for the normal operation (Figures 10BC).