Computational Fluid Dynamics data for improving freeze-dryers design

Computational Fluid Dynamics (CFD) can be used to simulate different parts of an industrial freeze-drying equipment and to properly design them; in particular data concerning the freeze-dryer chamber and the duct connecting the chamber with the condenser, with the valves and vanes eventually present are given here, and can be used to understand the behavior of the apparatus allowing an improved design. Pilot and large scale freeze-drying chambers have been considered; data of a detailed simulation of a complete pilot scale apparatus, including duct and condenser, are included. Data on conductance of an empty duct with different L/D ratio, on disk valves with different geometry, and on mushroom valve are presented. Velocity, pressure, temperature and composition fields are reported on selected planes for chambers and valves. Results of dynamic simulations are also presented, to evaluate possible performance of monitoring devices in the chamber. Some further data, with detailed interpretation and discussion of the presented data can be found in the related research article by Barresi et al. [1] and Marchisio et al. [2].


Specifications
Chemical engineering More specific subject area Fluid dynamics, Freeze-drying, Pharmaceutical technology Type of data  Correlations for pressure drop and maximum flow for mushroom valves, as a function of the valve opening distance are given.
Data for critical mass flow are shown in term of mass flux, to allow validation with experimental data.

Data
Data concerning the influence of chamber geometry, number of shelves and clearance between shelves, position of the duct leading to the condenser, number and position of the inert gas injection nozzles on fluid dynamics, water mass fraction and pressure distribution in the chamber, are first presented. Data on fluid dynamics of ducts and valves, and of a pilot scale complete apparatus, are then shown. The effect of the boundary conditions selected (including the slip or no-slip assumption), depending on the Knudsen number, and of actual entrance conditions is also documented. Graphs for estimation of critical mass flow in empty ducts and in presence of butterfly and mushroom valves are given. Design charts and correlations for equivalent length of valves developed from data here presented are given and discussed in Ref. [2]. Table 1 summarizes the simulation conditions adopted, the variables considered and the data made available.

Freeze-dryer chamber
Two different equipment sizes (a small-scale apparatus and an industrial-scale apparatus) have been investigated (see Fig. 1). For the small-scale apparatus, three different positions for the duct have been considered: in the center of the rear wall, on one side of the rear wall and on the bottom of the chamber (details of the geometries are shown together with the fluid dynamic data). The laboratoryscale freeze-dryer is constituted by four shelves (450 Â 455 mm 2 ) for a total chamber volume of 0.2 m 3 . The large-scale freeze-dryer contained 14-17 shelves of 1500 Â 1800 mm 2 for a total volume of 10.3 m 3 (with a 200 mm lateral channel) and a 800 mm diameter horizontal duct (duct length/duct diameter¼2).
Three different distances between the product and the upper shelf have been considered for the small-scale apparatus (with 4 shelves) and four for the large-scale apparatus. In the industrial chamber the number of usable shelves has been varied from 14 to 17, but the position of the first and of the last shelf was not varied. In Table 2 the values of the distance between the shelves, h, and of the clearance between the product and the upper shelf, r, considered for the simulations for both the small-and large-scale apparatus have been reported.  As concerns the butterfly valve simulations, just a short piece of straight duct 830 mm long and with a diameter of 700 mm has been considered, with a butterfly valve in open position. Two different valve shapes have been analyzed: in a first series of simulations a simplified geometry has been considered, represented by a round disk (685 mm in diameter) with uniform thickness assumed equal to 40 mm. Then the butterfly valve (again represented in open position) has been reproduced more faithfully to its real geometry (see Fig. 1). The overall diameter of the valve is still 685 mm, but the valve is represented as the union of three rings of different diameter: an outer ring about 80 mm wide (from radius 342.5 mm to radius 264.4 mm), presenting a thickness of 12 mm at the valve outer border, that increases linearly up to 16 mm; an intermediate ring with a starting thickness of 16 mm (at the radius of 264.4 mm) and a final thickness of 52 mm (at the radius 235 mm), with a non-linear profile; and a central "spherical cap" with a radius of 235 mm, with a starting thickness of 52 mm (at radius 235 mm) and reaching a maximum thickness of 125 mm in its center. These thickness values refer to the whole valve but it is also necessary to consider the intersection of the valve with the pipe used for its movement, presenting a diameter of 125 mm.
The mushroom valve is modelled as a flat disk (diameter¼ 750 mm; thickness ¼ 80 mm). A standard horizontal condenser, characterized by a standard toro-spherical bottom (DIN 28011) and a conical inlet, equipped with a mushroom valve, has been considered. The smallest section is at the intersection of the entrance cone with the vessel bottom, and has a diameter D¼700 mm; this geometry is compatible with the industrial freeze dryer considered. Only the first part of the apparatus has been considered, where there is no ice formation, and thus no material sink.

Complete pilot-scale freeze-dryer
Some simulations have been carried out on the whole pilot-scale apparatus (LyoBeta ™ from Telstar, Terrassa, Spain), including the real chamber described previously, the duct and also the condenser (see Fig. 1). A simple flat disk butterfly valve has been considered in the duct, which ends with a bend in the cylindrical condenser. Two different chamber-duct connections have been compared, a sharp entrance and a rounded one, corresponding to the actual realization (details of the meshes will be shown with the data).

Computational details
Governing equations used in the CFD model can be found in the research articles [1,2]. The three-dimensional simulations carried out to obtain the data presented are based on structured computational grids of hexahedral cells, representing the geometry of the freeze-drying chamber, by using standard numerical methods; successive grid refinement was used to verify the grid independence of the solution for all the cases investigated. All these pieces of equipment, geometrically reproduced and meshed, have been modelled by means of Computational Fluid Dynamics, with the commercial CFD code Ansys Fluent that, in certain cases, has been coupled to User Defined Functions (UDF) developed on purpose.
The continuity and Navier-Stokes equations are solved by resorting to a finite-volume discretization. The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm was used to solve the pressure-velocity coupling, whereas, in order to contrast the insidious effects of numerical diffusion, the Quadratic Upwind Interpolation for Convective Kinematics (QUICK) was used, by refining an initial solution obtained with a first-order upwind interpolating scheme [3]. Very restrictive convergence criteria were used, with normalized residuals smaller than 10 À 6 , and, in some cases, small under-relaxation factors needed to be used to reach convergence.
It is important to highlight here that for the operating conditions considered, the Knudsen number, evaluated estimating the vapor mean free path from kinetic theory and using the shelf clearance as characteristic dimension, was generally verified to be sufficiently small (of the order of 0.01) to assume that the gas flows in the continuum regime. For a number of cases (in particular for the small scale apparatus with the smallest clearance, where the Kn was higher) and for the butterfly valves, the possibility of imposing slip boundary conditions on the walls (typical of the transition regime between the continuum and molecular flow) was tested. As the results evidenced that, for the geometries and the operating conditions investigated, this has little effect on the predictions of the pressure drops, and almost no effect on the flow field predictions, the main part of the simulations and in particular those for the large scale apparatus were carried out with no-slip boundary conditions. Steady-state simulations were generally carried out, considering the water vapor as a compressible fluid, whose density is evaluated according to the ideal gas law, whereas the viscosity, μ, is calculated  with the standard kinetic theory [4]: where M is the mass of the molecule (expressed in kg), k B is the Boltzmann constant (1.38066 Á 10 À 23 J K À 1 ), T is the temperature (K) and d is the hard-sphere parameter for the molecule of mass M. Flow field data have been also reported in dimensionless form, as Mach number (Ma), defined as the ratio between the gas velocity and the sound speed estimated as follows: where U is the flow velocity and a is the speed of sound of the fluid (a ¼ ffiffiffiffiffiffiffiffi ffi kRT p ). From the ideal gas theory a value k ¼c p /c v ¼4/3 is assumed for the specific heat ratio of water vapor.

Freeze-dryer chamber
The grid was made of about 300,000 (for the small-scale apparatus chamber) or 600,000 (for the large-scale apparatus) cells. Inlet boundary conditions were set for the sublimation surfaces, corresponding to the vials placed on each shelf, while a standard pressure-outlet boundary condition was used for the final section of the duct connecting the chamber to the condenser.
The layer of vials has been modelled as a slab, with a thickness of 43 mm corresponding to that of the vials partially stoppered, with an equivalent uniform vapor source on the upper side.
Preliminary simulations showed that the flow field predicted with detailed (array of vials) and simplified (slab) mass-flow-inlets is identical, mainly due to the pressure drop that the water vapor has to overcome to flow through the clearance between the shelves [5]. Therefore, the simple continuous slabs have been considered as mass-flow-inlet boundary conditions. A uniform mass flux was thus considered at the upper surface of the slab, corresponding to the actual average sublimation flux obtained over the shelf (that depends on the real sublimation rate per unit surface of product, on the vial arrangement, and vial wall thickness).
Most of the cases investigated included only water vapor as solvent, but in a few selected test cases injection of inert gases (i.e. nitrogen) at different flow rates to control the chamber pressure, was also simulated.
The CFD simulations have been run considering typical values of the operating parameters: operating pressure of 10 Pa (value imposed at the end of the portion of the duct considered in these simulations) has been generally set, but a few simulations have been carried out also at different pressures.
The chamber walls were considered at 283 K, with a shelf temperature of 258 K. In some simulations concerning the whole apparatus the case of adiabatic walls has been also considered; in any case, it must be noted that the condition at the external walls does not influence the pressure distribution over the shelves, which is strongly dependent on the clearance, pressure, and sublimation rate.
Steady-state simulations of the primary-drying phase have been carried out considering the maximum available sublimation rate of 1 kg h À 1 m À 2 with an interface product temperature of 239 K; other simulations have been also carried out at different sublimation rates (from 0.50 kg h À 1 m À 2 to 0.90 kg h À 1 m À 2 ). The temperature of the vapor was assumed equal to the temperature of the sublimating surface (T i ).
For the simulations concerning the mushroom valve section, the data refer to the maximum flow estimated for the large scale equipment: in particular the configuration with the maximum number of shelves (17), corresponding to 54 m 2 , and the conventional maximum sublimation rate (1 kg h À 1 m 2 ), leading to 0.015 kg s À 1 has been considered as reference.
Some scalar transient simulations have been also carried out to investigate the dynamic response of the small apparatus to variations of the sublimation rate on the different shelves.
As concerns the effect of radiation on the gas fluid dynamics, radiation contribution of the gas has been also taken into account by means of P1 radiation model. The P1 radiation model is the simplest case of the P-N model, which solves the radiation transfer equation (RTE) by means of the expansion of the radiation intensity into an orthogonal series of spherical harmonics.

Choked flow in ducts
For the duct simulations (either empty or with the butterfly valve), the inlet of the pipe is modelled as a pressure inlet (the stagnation pressure and temperature are imposed and calculated by assuming an isentropic compression), whereas the outlet section is modelled as a pressure outlet (the outlet pressure is imposed), considering for all the cases an inlet static temperature of 239 K.
Simulations are three-dimensional in space and steady-state in time. A preliminary study was carried out to establish the grid size for which no influence on the results was observed; a mesh with 80,000 nodes was necessary for the longest duct (L/D¼50). The convergence criteria for the residuals were set at 10 À 6 .
In order to determine the critical mass flow, a series of simulations has been carried out at constant inlet pressure and different values for the outlet pressure, reducing progressively the outlet pressure, to identify the value at which no further increase in the flow occurred; the results confirmed that the critical flow conditions could be reliably estimated. Then different inlet pressures have been investigated, repeating the previous procedure. It must be evidenced that certain static pressures, and Mach numbers, have been supposed on the inlet and the relative stagnation conditions have been calculated, but inlet and outlet pressure given at the beginning are considered approximated values, and are adjusted by the code: a relevant number of simulations have been lunched to get a good estimation of the critical conditions.
To investigate the relevance of the boundary conditions (adiabatic, or isothermal wall at 239 K) some simulations have been carried out; a small influence of the boundary conditions was observed (see Fig. 2); it has to be noted that for the wall the same temperature of the gas was used, but the inclusion of heat supply from the wall would make the results very dependent on the case considered. The jet flow calculations in this case underestimate the conductance of the duct, but only because null inlet velocity is considered. Due to the low pressure and short residence time, the assumption of adiabatic flow is realistic. Low pressure slip boundary conditions, that should be applied in case of transitional regime (as discussed in Ref. [1], have been also tested, but the results confirm that the assumption of viscous flow (and thus the no-slip boundary conditions)) are acceptable at least for larger ducts.
An almost flat velocity profile is considered at the inlet of the duct, as a uniform pressure is imposed. The velocity profile will develop in the duct, even if for small L/D values the profile at the exit will be still far from the parabolic one corresponding to fully developed flow. An example is shown in Fig. 3. In the short duct the profile remains flat, while in the long duct the velocity can increase, reaching sonic conditions, and a parabolic profile can develop.
In all cases the absolute (static) pressure is given, but it must be considered that two different values can be considered: the "contour pressure", that is the value measured on the contour of the duct. For the inlet it will give the maximum value at the entrance, and for the outlet the minimum value at the exit. the "surface average pressure", that is the average value calculated on the whole surface considered.
The data presented will generally refer to this pressure.

Conductance of butterfly valves
A procedure similar to that previously described for the empty duct was used for the simulation of the piece of duct containing the butterfly valve. The computational mesh used in the CFD simulations contains about 650,000 cells. Fig. 4 shows that if the inlet pressure is increased, keeping constant the outlet pressure, the mass flow increases, and keep increasing also after that choked conditions are obtained; the data refer to subcritical flow conditions, with an outlet pressure close to 4 Pa. In the example the influence of different boundary conditions is also shown; for a DN 700 valve the Knudsen number is sufficiently low to make acceptable the no-slip assumption, thus the slight effect is mainly due to the influence of temperature at the wall. Graphic correlations shown in the following have been developed for adiabatic wall and no-slip conditions.

Simulation of complete small scale equipment
The mesh used for the CFD simulations is made up of 619,401 cells in the case of the straight entrance and of 623,671 cells in the case of the rounded entrance. The incondensable gas pipe is not meshed, and its first section is modelled as the inert gas exit from the condenser.
The water vapor deposition in the condenser is modelled as a finite rate process taking place on the condenser refrigerated walls (its lateral walls); further details can be found in Ref. [6]. The two chemical species (water vapor and nitrogen) enter the chamber through uniform sources placed on the upper side of four slabs representing the layers of vials, where the inert gas represents the 5% of the overall entering mass flow. In the base case the sublimation rate considered is 1 kg h À 1 m À 2 , with an interface product temperature of 239 K, and the pressure of 4 Pa is imposed at the condenser exit (at the entrance of the incondensable gas pipe), modelled as a pressure-outlet. The low-pressureboundary slip option is activated in Fluent laminar model, considering the low pressure values and the small duct size.

Mushroom valve fluid dynamics
The duct inlet is modelled as a mass flow inlet, imposing the inlet mass flow; the stagnation temperature must also be imposed, and has been taken equal to 239 K: it must be evidenced that differently from previous cases (pieces of duct, empty or with the butterfly valve), where the stagnation temperature to be fixed at the inlet was calculated for the different flow rates in such a way that the static temperature was always the same, here the same value of stagnation temperature is considered for the different cases, and as a consequence the static temperature may change at the inlet (the variation is of about 5 K between the lowest and the highest flow rate considered). The choice of constant stagnation temperature allows having the same static temperature in correspondence of the center of the disk, where the flow is stopped.
The outlet section, a plane positioned in the vessel behind the mushroom disk, is modelled as a pressure outlet; the outlet pressure is set equal to 4 Pa in all cases. It can be reminded that practically the total pressure is constant in a condenser, as the reduction in partial pressure of water vapor is compensated by the increase in the partial pressure of inert gas. The low pressure boundary slip conditions have been adopted; the walls of duct and vessel are modelled as adiabatic surfaces and radiation is neglected.

Analysis of the Knudsen number
The evaluation of the Knudsen number, defined in Eq. (3) is necessary in order to assess the validity of the continuum approach; this is very important in the freeze-drying chamber since the governing equations are often pushed to their limit of applicability under common operating conditions, possibly resulting in an over-estimation of pressure drops and fluid velocities [1]. . The slab that represents the volume occupied by trays or vials is shown, and the size of the free clearance r, is evidenced. Velocity magnitude (m/s) plotted on x-z plane (middle graphs) and velocity magnitude (m/s) plotted on y-z plane for the S1 case (r¼ 57 mm) for the three configurations; a mass flow rate equal to 4.5 Á 10 À 5 kg s À 1 , and outlet pressure 10 Pa, corresponding to typical values experimentally measured in lab operation has been considered in this case. If Kn is very small (Kn o0.01) the gas flow is in the continuum regime, and the standard continuity and momentum balance equations can be used, while when Kn is large (Kn4 1) the gas is in the socalled molecule regime and the BE has to be solved instead. If the mean free path of the gas molecules is neither very large nor very small, as compared to the macroscopic length-scale of the flow, the system is under the so-called transitional regime. In this case the boundary conditions for the governing equations are modified to take into account for velocity slip and temperature jump at the walls [7]. This is needed as the gas-phase velocity at solid surfaces differs from the velocity at which the wall moves and the gas temperature at the surface differs from the wall temperature.
For the cases reported in Table 2, the parameter r can be considered the characteristic dimension of the system. In order to evaluate the dependence of the Knudsen number on the clearance r it is necessary to evaluate the mean free path (λ) of the vapor molecules. From the kinetic gas theory, the mean free path can be written as in the following equation: where k B is the Boltzmann constant, T is the temperature, d is the Lennard-Jones characteristic length of the molecule and P is the absolute pressure. It can be observed that λ is proportional to the temperature and inversely proportional to the pressure so, if the temperature is maintained constant, decreasing the pressure value the mean free path increases. In Fig. 5 the Kn number depending on r has been reported for some of the cases investigated. The Knudsen number decreases if either the distance between the shelves or the pressure is increased.
The data confirm that, with the conditions considered, the regime is generally laminar, thus the simulations give reliable results, but moving to smaller free distances, at the lowest pressure, it is possible to enter in the transition regime. In particular, for very low values of r pure viscous flow (with Kn o0.01) occurs only at higher pressure; for higher clearance values the pressure can reach smaller values but also for this condition below 10 Pa the transitional regime for the vapor flow occurs.

Data analysis
3.1. Fluid dynamics in the pilot scale apparatus. Influence of clearance and duct position on flow field and pressure distribution In Fig. 6 the contour plots of the velocity magnitude are reported on two different planes for all the investigated positions of the duct; the x-z plane is positioned perfectly in the middle of the chamber, while the position of the y-z plane has been chosen in order to keep the same distance from the duct and the plane for all the configurations. For these contour plots it can be observed that the maximum of the velocity is reached near the duct; moreover, when the duct is positioned at the bottom of the chamber, stronger velocity gradients are observed. Fig. 7 shows the pressure distribution over the shelves for different duct positions and shelf clearances (the two ones with the duct on the rear wall are very similar).

3.2.
Temperature distribution in the chamber. Influence of inter-shelf clearance Fig. 8 shows the temperature distribution on two different orthogonal planes positioned in the middle of the chamber; no appreciable differences are observable for these clearances between the case with slip and no-slip boundary conditions (data not shown). In the upper graphs the walls are all at 283 K and it can be seen that the vapor in the outer clearance zone is at higher temperature and close to the wall temperature, but the gas temperature in the zone between the shelves is determined mainly by the sublimation process (temperatures are in the range 240-260 K, with the product at 239 K).
Changing the shelf clearance, the situation remains similar. In the bottom graphs a different boundary condition is adopted for the top surface (an adiabatic wall is considered).
The area-weighted average temperature calculated on the two planes shown resulted to be 266.59 and 258.71 K in the two cases (266.35 and 258.56 K respectively with slip boundary conditions).

Fluid dynamics in the large scale apparatus. Influence of clearance on flow field and pressure distribution
In order to give a general idea of the fluid dynamics in the chamber, the streaklines of the vapor flowing from the shelves to the duct connecting the drying chamber to the condenser can be used to visualize the fluid trajectory in the drying chamber. An example is given in Fig. 9: as it is possible to see from the velocity vectors, the vapor sublimating from the vials is forced to flow toward the edge of the shelf and, from there, it directly goes into the duct positioned on the rear wall, to be then collected in the condenser. To design purposes, it can be useful to quantitatively evaluate the flow distribution between the different sides of the shelves, to verify if the clearance between the shelf-pack and the walls is appropriate; the geometrical solution which gives a uniform flux through the lateral zones of the shelves will assure the minimum pressure drop in the chamber. In Fig. 10 the coefficients α xi and α yi as functions of the z position are shown (0 add 1 in the pedex refer respectively to the side closer to and farer from the axis origin. These coefficients represent for each shelf the fraction of vapor flowing out of the selected side and are calculated as the ratio between the mass flow rate m xi , m yi and the total mass flow rate of vapor; these values have been obtained by integrating the mass flux obtained by CFD simulations over the surfaces that delimit the shelf volume, and can be useful if the pressure drop over different shelves has to be compared. Vapor rate fraction coefficients sum to unity: It is interesting to observe that the coefficients of the vapor distribution relative to the x-direction are not affected by the number of shelves, while this effect exists, even if it is weak, in the y-direction; moreover, for the plates at the bottom of the chamber, far from the duct, an equivalent vapor distribution through the four sides of the shelves is observed (it must be remembered that one side is 20% longer than the other one); moving up, the value of α y1 increases and reaches the maximum value of 60% for the shelves close to the duct. Fig. 11. Influence of the clearance, r, on the location of maximum pressure, y max , for the shelves at the level of the duct (a), and on ΔP max , for three different shelves, numbered from bottom (b): ▲ , third shelf; ○ , sixth shelf; • , shelf close to the duct. On the lower graphs, the variation of y max (c) and ΔP max (d), with the z-coordinate for the L3 configurations at two different sublimation fluxes: Δ, 1 kg h À 1 m À 2 ; ▲ , 0.5 kg h À 1 m À 2 .
In the apparatus considered the duct is on one side, thus pressure profile over the shelves is not symmetric, and changes with the shelf position (the farer from the duct, the more symmetric they are: see also Figs. 4 and 6 in Ref. [1]). In Fig. 5 in Ref. [1] y max , the y-coordinate of the maximum pressure value for each shelf, and the corresponding maximum pressure difference over the shelf are reported for all the configurations, as a function of the vertical position of the shelf.
The dependence of these quantities on the inter-shelf clearance is explicitly shown in Fig. 11. In particular, in Fig. 11a it is shown the variation of the location of the maximum pressure in the shelves at the level of the duct. The maximum pressure drop (ΔP max ) over three different shelves is plotted as a function of clearance in Fig. 11b (at maximum sublimation flux of 1 kg h À 1 m À 2 ). ΔP max decreases when clearance r increases; it can be noted that the dependence of ΔP max on r is different for the different shelves. Very important is the effect of the sublimation rate. In Fig. 11c-d the effect of the mass flow rate on y max and ΔP max is evidenced for the case L3 (16 shelves). The major effect of the mass flow rate is to increase ΔP max , but this increase is larger for the shelves close to the duct.
In Fig. 12 (upper graph) the values of the pressure differences over the shelves and in the chamber are compared (for the L3 configuration, with 16 shelves). It can be noted that in the central zone of the shelves, pressure values higher than the reference one can be observed: the maximum value occurs in the lowest shelf, less affected by the duct, but the maximum pressure increase, for the configuration shown is in any case very small, lower than 0.5 Pa. In Fig. 12 (lower graph) The effect of the chamber reference pressure is shown.

High flow rate regime for the injected inert gas
The data presented have been obtained for a pilot scale drying chamber where the inert gas inlet has been supposed positioned in the top corner of one of the lateral walls of the chamber. Results concerning the range 1.20 Á 10 À 8 -1.20 Á 10 À 5 kg s À 1 , which covers the values normally observed during primary drying, have been presented in Fig. 7 in Ref. [1]. In the secondary drying, in particular, much higher values can be measured, as the desorption flux of the solvent has significantly lower values. Fig. 13. Inert gas mass fraction distribution in the small scale drying chamber for very high inlet flow rates (1.20 Á 10 À 4 kg s À 1 ). Outlet pressure¼ 10 Pa; water sublimation rate 1 kg h À 1 m À 2 . Velocity magnitude (upper left, m/s) and inert gas fraction (upper right, -) are shown on the plane of the inlet inert gas jet. In the bottom graphs the inert gas fraction (left, -) and the absolute pressure (right, Pa) are shown on two orthogonal middle planes. Red colour corresponds to maximum value for each case, yellow to 75%, green to 50%, light blue to 25% and dark blue to 0.   15. Evolution of interface temperature and sublimation flux, J w , during primary drying, depending on the local gas composition in the chamber (% of inert gas fraction). 1 ml of 10% w/w sucrose solution in vial (14.25 mm internal diameter); chamber pressure ¼ 10 Pa, T shelf ¼ À10°C. Contour plots (on the jet plane) of inert gas fraction and velocity magnitude for the highest inert gas bleed rate investigated (1.2 Á 10 À 4 kg s À 1 ) are shown in Fig. 13 (top graphs). In this case extremely high velocities are obtained, and the jet reaches the opposite wall generating local high pressure values ( 450 Pa).
In Fig. 13 (bottom left) it is evidenced that, in the clearance between the shelves, the inert gas mass fraction remains low, while with this configuration higher values are observed in the lateral space. In the bottom right image the pressure is shown in two cross planes, where it ranges from about 20 (blue) to 24 (red) Pa; for comparison, it can be mentioned that, without inert gas injection, a quite uniform 14 Pa value is obtained.

Influence of inert gas fraction on the heat transfer coefficient
Even if the average value of the inert gas fraction is low (but not null) during primary drying with gas-bleeding pressure control, depending on the geometric configuration, locally relative high values of inert gas fraction can occur. As the heath transfer coefficient between shelf and vials depends on the gas composition, this fact can affect the uniformity of the batch. In any case the heat transfer during secondary drying will be strongly affected, and this must be taken into account in cycle design.
An example of K v calculated for various values of total chamber pressure and of nitrogen percentage in the drying chamber is shown in Fig. 14. The case study here analyzed is the freeze-drying of a 10% w/w sucrose aqueous solution in tubing vials. Values shown in the figure have been calculated for a vial having an internal diameter of 14.25 mm, with a wall thickness of 1 mm, a bottom thickness of 0.7 mm, and the maximum gap between the bottom and the shelf of 0.4 mm; the shelf temperature is assumed to be equal to 253 K. The operating conditions used for the calculations are: P c ¼10 Pa, T shelf ¼ À10°C. It can be seen that the higher is the nitrogen percentage in the drying chamber, the lower is the value of K v .
An example of the influence of the local inert gas fraction on the drying behavior of the product in vials is shown in Fig. 15. A relatively wide range is considered for the nitrogen fraction, for a sensitivity analysis; the pure water vapor case can be representative of processes where chamber pressure is maintained acting on the vacuum pump, neglecting the leakages from the environment.

Transient chamber modelling
The details of the conditions of all the simulations carried out are given in Table 3. In the simulations D1-D6 all the shelves are considered loaded with product and then all the vapor sources are considered operating; in these cases the flow field in the drying chamber is the same obtained with the steady-state simulations. For the simulations D7 and D12 only the shelf at the top of the chamber (ss 4 ) has been considered operating; in these cases new steady-state simulations are carried out before adding the scalar because the flow field of the previous simulation is no more representative for this system.
In Fig. 16 cases D1 and D7 are compared; in these simulations the full load and the partial load (but on the top shelf, close to the virtual detectors) of the chamber are compared. In the graphs, the dynamics of the system in different positions is shown. The dynamics is quite similar when the response is detected in points 2, 3 and 4 (thus only point 2 is shown for brevity); this aspect is very important for the sensor response interpretation. In point 1 the delay for the partial-load system (D7) is higher than for the fully loaded system (please note that the time scale on the axis is different). This results can be explained comparing the flow-field of the two systems, in fact the velocity in D1 is higher than in D7. Fig. 17. Response in the four detection points for the fully loaded system, with disturbance applied to all the shelves. Response for different sublimation rates are compared: -, J w ¼ 1 kg h À 1 m À 2 (D1); ---, J w ¼0.5 kg h À 1 m À 2 (D8); À Á À, J w ¼0.1 kg h À 1 m À 2 (D9).
For the sensor response interpretation it is important that differences in the sensor response are small in full and partial-load cases; for sensor location 1 delays are slightly different (even if differences are relatively small at high sublimation rate). But it must be underlined that this consideration is true only if the load is on the upper shelf ss 4 ; in fact, in case D4 no signal can be detected in points 2, 3 and 4.
The effect of the sublimation flow rate on the dynamic response of the system has been assessed with the simulations D1, D8 and D9 for the full-load system. In Fig. 17 the results of these simulations are reported. Decreasing the sublimation rate both the response delay and the rise time increase.
Similar results were obtained with partial load (simulations D7 and D12) (see Fig. 18). The delay increases significantly at low sublimation rates, but at very low mass flow rate the effect of scalar diffusivity becomes important. In Fig. 19 the dynamic response calculated taking or not the diffusivity into account is compared. When the scalar diffusivity is accounted for, the system response is faster; Fig. 18. Response in the four detection points for the single shelf loaded system (ss 4 ,), with different sublimation rates: -, J w ¼1 kg h À 1 m À 2 (D7); ---, J w ¼0.1 kg h À 1 m À 2 (D12). moreover the response detected in the four different positions is quite similar. Further data and discussion can be found in Ref. [1].

Maximum flow in ducts: dependence on duct length, diameter and pressure
The influence of the duct length on the maximum available vapor mass flow has been investigated considering different L/D ratios: the range 1-5 is the one of largest interest, and in literature data up to now available are limited to this range, but here data are presented in an extended range from 1.2 to 50, to allow checking the validity of the equivalent length concept for valves (see Ref. [2]). The simulations have been carried out for DN 700, but the results have been generalized using the mass flux; eventually a correction factor can be applied for smaller ducts.
All the plots refer to surface average pressures.
The results obtained are shown in Fig. 20 (data at L/D¼5 and 7 have been omitted for sake of clarity). It is evident that a linear relationship exists between the inlet (chamber) pressure and the critical mass flow density.
The relationship between pressure drop and flow rate, and the maximum allowed flow in choked conditions, depends on the characteristics of the compressible gas, that is on the inert gas fraction in the water-nitrogen mixture. Here some data obtained with "jet flow" calculation are presented in Fig. 19. Response in the four detection points, without species diffusivity (simulation D9, line) and with species diffusivity (D10, symbols); fully loaded system, low sublimation rate, J w ¼ 0.1 kg h À 1 m À 2 . order to show which the influence of the different factors is. The jet flow depends on gas temperature and inlet Mach number (see [2]). Qualitatively the results for a water-nitrogen mixture are very similar to those obtained for pure water, but the total mass flux increases in presence of inert gas.
The critical mass flux increases with the initial Mach number, as can be seen in Fig. 21. The figure also shows how the critical flux increases progressively when the fraction of nitrogen gas increases; the molar flow rate is higher for pure water vapor and decreases when the inert gas fraction increases.
From the operating point of view, it is interesting to evaluate which is the influence of the inert gas on the mass flow rate of water that is possible to remove in the new conditions, as this is related directly to the sublimation rate in the drying chamber.     water mass flux as a function of the inert gas fraction in the mixture, at different Mach number, for a given chamber pressure (similar results have been observed for other pressures). Both quantities increase with the chamber pressure, but while the total mass flow rate increases with the inert gas fraction in the mixture, the water mass flow rate is reduced significantly, and more than proportionally. Fig. 25. Zoom of the contours of the water mass fraction in the condenser, for the case with 5% inert gas (left) and 1% inert gas (right). Sublimation rate¼ 1 kg h À 1 m À 2 . Fig. 26. Velocity vectors on different condenser sections, taken from the bottom to the front: (a) at the entrance; (b) at 9 cm from the entrance, towards the front; (c) at 13 cm from the entrance, towards the bottom; (d) at incondensable gas pipe entrance. Sublimation rate ¼0.4 kg h À 1 m À 2 , with 5% w/w inert gas.
In the jet flow calculations, the inlet Mach number, and thus the inlet kinetic energy of the flow, has been taken into account, as the critical flux is significantly affected by this variable; but while this is easy to estimate in case of flow through a hole, because it corresponds to the velocity in the vessel (and is generally close to zero), it is much more complicate to estimate it in a duct. But this is not really a problem, because the interest of the designer is in the reduction of the water mass flow that occurs in case inert gases are present. All the curves scale in the same way, thus the relative reduction in the water mass flux is the same independently of the inlet Mach number and chamber pressure.  If the behavior of straight ducts is analyzed, using CFD simulations, it appears that similar results are obtained (see Fig. 21). In this case the pressure at the inlet and outlet is set, thus the inlet Mach number is not a free variable, but is determined by the flow rate and the pressure drop in the duct. In   the example a relatively long duct is considered; it can be noted that in this case the jet flow calculation significantly overestimates the mass flux, even if zero inlet velocity is considered. It can be seen that in case of adiabatic flow the temperature decreases significantly, and the temperature of the wall (assumed equal to the inlet static temperature) has some influence on the developing of the velocity profiles and finally on the allowable mass flow. For a DN 700 valve the Knudsen number is sufficiently low to make acceptable the no-slip assumption, thus the slight effect is mainly due to the influence of temperature at the wall.

Flow around butterfly valves
In the right hand column, the corresponding pressure, temperature, velocity and Mach number profiles obtained for the valve with profiled disk are shown.

Inlet effects in the large-scale apparatus
As concerns the large-scale apparatus, just a straight duct, with a L/D¼2 has been considered. Fig. 23a shows a sketch of the apparatus, with the exact location of the planes where radial pressure and velocity profiles are given. Fig. 23b shows a 3D contour plot of the pressure profiles along the duct, showing the large gradients at the chamber-duct connection, and the nonuniformities still remaining in the duct (the radial pressure profile are shown in Fig. 6a in Ref. [2]). Fig. 23c shows the pressure profile along the centerline of the duct, at different value of the mass flow rate (different configurations, with variable number of loaded shelves, and different sublimation rates have been considered). The end of the inlet zone is approximately at y/D ¼1, slightly affected by the mass flow rate as evidenced in Fig. 23d. It can be noted that a new inlet occurs, due to the development of the final parabolic profile.
Also the velocity profile is very different from the ideal one in a simple empty duct: Fig. 23e shows the radial velocity profile in the same planes of Fig. 23a and b (the zone close to the wall has been omitted for sake of clarity, but the slip effect is negligible). It can be seen that the profile has two lateral maximums, with the minimum on the centerline; this is a consequence of the particular velocity profile in the chamber clearance, but it becomes even more pronounced at the entrance plane where the maximum velocity almost doubles. It can be noted that in the following sections these strong differences disappear, but the profile remains irregular and is still asymmetric also at the exit.
3.8.2. Conductance of the lab-scale freeze-dryer. Influence of the shape of the entrance The whole apparatus, including the condenser chamber, has been modelled in this case. Fig. 24a shows the details of the duct, with location of the plane where the radial pressure profiles shown in Fig. 6 (bottom graphs) of the article [2] are taken.
Two different entrances have been compared: a straight entrance and a rounded one. Two simulations have been carried out to investigate this effect, differing only in this geometric feature.
The details of the meshed geometries are shown in Fig. 24c and d, for the two cases. The effect of the entrance type on the values of pressure and axial velocity along the duct axis is shown in Fig. 24e  It must be evidenced that the profile drawn is just that one on the axis; the pressure at the wall keep decreasing, while the radial profile changes continuously (see Fig. 6 in Ref. [2]). There the radial profiles at different distances are shown: four in the chamber clearance, and 13 in the duct, each one distanced of 10 mm (for sake of clarity profiles corresponding to line 6, 8 and 9 have been omitted).
The influence of the shape of the entrance is shown in Fig. 24b: in case of sharp edges the entrance profile presents two lateral maxima, similarly to what observed in the larger apparatus, even if less marked. The rounded edge generates a flatter profile; the maximum velocity is lower, but it must be taken into account that the inlet area is slightly larger just at the inlet in this case.
The profile at the end of this section (at z¼ À0.13 m, with the origin of the axis in the chamber below the bottom shelf), just before the valve, is also shown; the velocity profile is still developing and some differences between the two cases are still evident.
The profile in the second part of the duct (and then in the condenser) can be seen in Fig. 24g; the gap in the graph is due to the presence of the butterfly valve. The complex behavior is again related to the developing pressure and velocity profile after the valve, with the change of section; the influence of the shape of the duct entrance is no longer significant.

Fluid dynamics in the condenser
In Fig. 24g the profile of the absolute pressure along the whole path, is also shown, from the duct inlet to the condenser exit (as the center of the condenser is occupied by the exit pipe of the inert gases, the pressure is taken in the middle of the annular zone, as shown in Fig. 24h); the distance covered along the path is reported on the abscissa axis.
Different symbols and colors are used to represent the profiles along the different tracts of the path: the duct, the radial tract after the condenser entrance, the tract along the condenser (line 1 condenser), parallel to the condenser axis, the second radial tract (line 2 condenser) and the line in front of the incondensable gas pipe entrance (line entrance incond. pipe). This confirms that in this case there are no significant pressure gradients in the condenser, except just at the outlet.
What changes along the condenser is the distribution of the water mass fraction, due to the icing at the wall surface. This can be clearly seen in Fig. 25, where two cases, with a different inert gas fraction, are compared; as mentioned in Refs. [1,6], the inert gas is responsible for a blanketing effect in the condenser.
CFD is also useful in this case to evidence the flow pattern in the apparatus. Fig. 26 shows, as an example, the velocity vectors of the gas mixture on different condenser sections.

Mushroom valve
Fluid dynamics data for different values of opening valve distance and mass flow rate, depending on the sublimation rate on the chamber shelves, are shown; finally a graphic correlation for the critical water flux as a function of the condenser pressure will be shown. Fig. 27 shows the contour plots of the velocity along the plane x ¼0 for different values of sublimation rate: the values observed are proportional to the mass flow rate, but the qualitative behavior is similar; in all cases higher velocities (and corresponding higher mass flows) occur in the plane y¼0 in the outer region. Fig. 28 shows for one of the valve distances the pressure distribution in the system, for different values of the sublimation rate. The strong reduction in velocity caused by the presence of the disk determines a significant increase in the absolute pressure that locally, in correspondence of the disk, can be higher than the inlet value. Fig. 29 shows the contour plots of the absolute pressure for the three configurations investigated, evidencing that the valve distance significantly affects the pressure drop, and in particular the pressure in the inlet zone: absolute pressure increases significantly in front of the valve disk, and the shorter the valve distance (that is the smaller the clearance for the passage of the vapor) the higher the pressure at the inlet; differences in the rear zone are much weaker.
To get a deeper insight into the hydrodynamics of the apparatus, and to better understand the influence of the vapor mass flow rate, the profiles of the variables of interest are shown. Fig. 30 (upper graphs) compares the axial velocity of the vapor along the axis, from the inlet of the conic duct up to the first tube bundle section, for the two extreme valve-position configurations. It is evident that while at the lowest sublimation rate, when the flow is subcritical, the velocity is significantly lower than in the other cases, the three curves corresponding to the higher sublimation rates are almost overlapping. This apparently strange behavior can be easily explained, considering the results already shown in Fig. 28, from which it can be seen that the pressure in the inlet duct increases with mass flow rate; as the inlet pressure increases practically in a linear way, then the increase in mass flow rate is compensated by the increase in the density value, leading approximately to the same velocity value. This is true not only on the axis, but in all the zone upstream the disk, while around and beyond the disk the flow field velocity is dependent on the mass flow rate (as shown in Fig. 27). Fig. 30 (bottom graphs) compares the absolute pressure axial profiles for the two configurations; it is evident that while behind the valve zone the pressure remains constant and practically equal to the exit value, for all flow rates, in the inlet it increases significantly with the mass flow rate. The increase is higher for the valve with a lower disk distance, because this causes a stronger flow restriction, but the variation in the zone before the disk (that is similar to that observed for the axial velocity) is larger in the valve with larger clearance.
It is interesting to compare the behavior of the mushroom valve with that of the butterfly valves investigated before. Fig. 31 shows that the conductance of mushroom and butterfly valves is comparable. It appears in particular that the conductance of the simple butterfly valve with flat disk in intermediate between that of the mushroom valve with smallest and largest clearance, while that of the butterfly valve with profiled disk is smaller.