Comparison of Knudsen Diffusion and the Dusty Gas Approach for the Modeling of the Freeze-Drying Process of Bulk Food Products

Abstract

Freeze-drying is generally used to achieve high quality products and preserve thermal sensitive components; however, it is also considered as a high energy and costly process. Modeling of the process can help to optimize the process to reduce these drawbacks. In this work, a mathematical model is presented to predict the heat and mass transfer behavior for freeze-drying of porous frozen food particles during freeze-drying to optimize the process. For the mass transfer, a comparison between Knudsen diffusion and the more complex dusty-gas approach is performed. Simulation results of a single particle are validated by experiments of single-layer drying to extend the usage of this model from a single particle to a particle bed. For the moisture transfer, adaption parameters are introduced and evaluated. A comparison shows a good agreement of the model with experimental results. The results furthermore suggest a strong correlation of the drying kinetics with pore size and particle porosity. An increase in the pore diameter strongly improves the overall mass transfer rates and hence is a suitable parameter for an effective increase of the drying rates in freeze-drying.

1. Introduction

Freeze-drying is used to dry sensitive products in the food and pharma industry [1] and is seen in comparison to other drying processes as a rather energy and cost intensive technology [2]. Especially in the food industry, it needs to be ensured that in addition to the quality, the drying times are held as short as possible to be competitive with other drying technologies. Freeze-drying consists of three steps: freezing, primary and secondary drying [3,4]. For liquid food materials, e.g., soluble coffee or tea, a part of the water crystallizes and the solution is concentrated during the freezing step. The remaining water is called non-freezable or bound water and remains in the glassy matrix. During primary drying, the crystalline water is sublimated, while during secondary drying the bound water needs to be desorbed from the amorphous food matrix. Both drying steps follow different drying mechanisms and for both, the dried layer forms a resistance for the heat and mass transfer within the drying particle. The understanding of these limiting effects allows for a prediction of the drying times and therefore the possibility to predict the drying kinetics by modeling, which is a powerful and cost-effective tool for process optimization.

Several model approaches were used in the past to describe the freeze-drying process. These mainly focused on drying in vials used predominantly in pharma industry or on tray freeze-drying used for the drying of bulk materials. The approaches range from single particle models to complex pore network models of particle beds on trays.

Many authors studied the freeze-drying process in packed beds: Liapis and Bruttini [5] considered also the freeze-drying of particles in packed beds, where they extended the application to spray-freeze drying which received broad interest [6,7]. The authors used the dusty gas model for the mass transfer, already applied by Song and Yeom [8] and Sheehan and Liapis [9]. Here, a mixture of inert gases and vapor was considered during the drying approach. In the work of Liapis and Bruttini [10], this model was further extended and validated by experimental data derived by Song et al. [7] and Her et al. [11]. The authors found out that the gas flux of inert gases is negligible compared to the vapor flux, during primary drying. The recent approaches of Warning et al. [12] considered the importance of a sublimation width instead of a single sublimation front, adding more precision, but also more complexity to the model. Gruber et al. [13] and Vorhauer-Huget et al. [14] confirmed experimentally that the sublimation front is not homogeneous, but affected by local changes in heat and mass transfer. The complexity applies also to the recent progresses of pore network models as applied by Vorhauer et al. [15]. The pore network model allows the coupling of pores with various sizes, allowing a more precise prediction of the vapor flux and its dependency on pore geometry and arrangement.

The influence of structure parameters on the heat and mass transfer rates were only insufficiently investigated and validated and therefore need further investigation. As discussed in previous publications [3,16], changes in the internal porous structure can have a significant impact on the drying time. Modeling approaches in this direction were already applied and vary strongly in complexity. In general, complex models are able to achieve higher precision, mostly because of a higher degree of freedom. These advantages usually come at the cost of increased computing power and simulation time. Furthermore, in each model different experimental results were used. A comparison between a simple Knudsen diffusion approach and an approach with higher complexity such as the dusty gas model can help to judge the experimental results according to existing models, which was done in this work. For the validation, the structural parameters such as porosity and pore size were varied to observe the change in drying time. Since the focus lies on the change of internal porous structure, the internal mass transfer limitation is of high interest, which is modeled and validated on a single-particle scale.

2. Methodology

2.1. Model Structure

A shrinking core approach is used to capture the main heat and mass transfer mechanisms of the single-particle freeze-drying process. Therefore, the model considers three distinct regions: a frozen core (index fc), a porous dried layer (index dl) and an infinitesimal particle surface region (index p). The shrinking core model is schematically presented in Figure 1. When the particle dries, the frozen core region shrinks and its boundary moves towards the center of the particle. The water vapor is transported through the pores of the growing dried layer region from the frozen core to the surface. The sublimation process is divided into two different stages: The primary stage considers a pore diffusion transport of the water vapor through the porous layer. For the second stage, a simplified approach of a first-order reaction kinetic is assumed to remove the remaining ice from the particle. The one-dimensional structure of the model allows only for a resolution of the heat and mass transfer along the radial direction. The radiative heat transfer from the bottom and top shelf as well as radiation from the environment to the particle surface are considered. Due to the low pressure environment, convective heat transfer at the particle surface is neglected. From the surface, heat is transferred by conduction through the porous dried layer to the sublimation interface of the frozen core region.

The main assumptions and simplifications of the model are listed below:

• Only the vapor transport within the porous structure of the particle is investigated, while the vapor transport within the drying chamber is neglected.
• The material is homogeneously distributed within each respective model region.
• The pores all have the same pore diameter and are connected to the outside of the particle. Closed pores are neglected.
• Material properties are independent of the temperature.
• Particles have a spherical shape.
• The influence of the sample holder is neglected and only radiative heat transfer is considered.
• Only radiative heat transfer in the drying chamber is considered. The particle is modeled as a gray body.

The frozen core and the dried regions are assumed to be ideally mixed, neglecting possible gradients in composition or temperature. The transient behavior of the mass and enthalpy is modeled by a CSTR (continuously stirred tank reactor) approach and are shown below:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{m}_{fc}}{dt}}& =-{\dot{m}}_{subl}-{\dot{m}}_{fc\to dl},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{H}_{fc}}{dt}}& =-{\dot{H}}_{subl}-{\dot{H}}_{fc\to dl}+{\dot{Q}}_{dl\to fc},\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{m}_{dl}}{dt}}& ={\dot{m}}_{fc\to dl},\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{H}_{dl}}{dt}}& ={\dot{H}}_{fc\to dl}-{\dot{Q}}_{dl\to fc}+{\dot{Q}}_{p\to dl}.\hfill \end{array}$$

(4)

m and H are the mass and the enthalpy of the respective region, $\dot{m}$ the mass flow, $\dot{H}$ the enthalpy flow and $\dot{Q}$ the heat flows between the regions. The mass flow and the enthalpy flow of sublimation, ${\dot{m}}_{subl}$ and ${\dot{H}}_{subl}$, respectively, are directly removed from the frozen core. The sublimation results in a shrinking of the frozen core region. While the vapor moves through the pores of the dried layer to the surfaces, the dried matter is removed from the frozen core region and is added to dried layer region, which is expressed by the exchange flows ${\dot{m}}_{fc\to dl}$ and ${\dot{H}}_{fc\to dl}$. The temperature of the frozen core and the dried layer region are calculated by:

$$\begin{array}{cc}{T}_{fc}\hfill & =T_0+\frac{h_{fc}-h_{fc,0}}{c_{\mathrm{p},fc}},\hfill \end{array}$$

(5)

$$\begin{array}{cc}{T}_{dl}\hfill & =T_0+\frac{h_{dl}-h_{dl,0}}{c_{\mathrm{p},dl}},\hfill \end{array}$$

(6)

with the specific enthalpy $h=H/m$, the specific heat capacity $c_{\mathrm{p}}$ and the reference temperature $T_0$.

The infinitesimal particle surface is modeled by a quasi-stationary enthalpy balance, which takes into account the incoming radiative heat transfer and the outgoing conductive heat transfer from the surface into the dried layer as follows:

$${\dot{Q}}_{rad}-{\dot{Q}}_{p\to dl}=ϵA_p\sigma (T_{\infty }^4-T_p^4)-k_{dl}·4\pi ·(T_p-T_{dl})·{\displaystyle \frac{r_p·r_{dl}}{r_p-r_{dl}}}=0,$$

(7)

with the emissivity $ϵ$, the Stefan–Boltzmann constant $\sigma$, the thermal conductivity of the dried layer $k_{dl}$ and the radius of the respective region r. The implicit equation is solved internally to calculate the temperature of the particle surface $T_p$.

The conductive heat flow between the dried layer and the frozen core is calculated by:

$${\dot{Q}}_{dl\to fc}=k_{dl}·4\pi ·(T_{dl}-T_{fc})·{\displaystyle \frac{r_{dl}·r_{fc}}{r_{dl}-r_{fc}}}.$$

(8)

The mass flow of dry material from the frozen core to the dried layer and the respective enthalpy flow can be described as a function of the vapor mass flow, which is released from the sublimation process as follows:

$$\begin{array}{cc}\hfill {\dot{m}}_{fc\to dl}& ={\dot{m}}_{subl}·{\displaystyle \frac{w_w}{1-w_w}},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\dot{H}}_{fc\to dl}& ={\dot{m}}_{fc\to dl}·c_{\mathrm{p},dl}·(T_{fc}-T_0).\hfill \end{array}$$

(10)

where $w_w$ is the water weight fraction. The sublimation enthalpy flow is calculated using the equation

$${\dot{H}}_{subl}={\dot{m}}_{subl}·\Delta h_{subl,w},$$

(11)

with the specific enthalpy of sublimation $\Delta h_{subl,w}$.

Calculation of the sublimation mass flow ${\dot{\mathit{m}}}_{\mathit{subl}}$

The total sublimation mass flow ${\dot{m}}_{subl}$ depends on the vapor transport through a single pore and the number of pores which connect the frozen core with the environment. It is calculated by

$${\dot{m}}_{subl}={\epsilon }_p·4\pi r_{fc}^2·{\dot{m}}_{pore}^{″},$$

(12)

where ${\dot{m}}_{pore}^{″}$ is the vapor mass flux through the porous layer. For the calculation of the vapor mass flux, two different models are applied that are subsequently shown. The Knudsen number $Kn$ is used as a measure for the characteristic movement of particles in confined regions and can be calculated by

$$Kn={\displaystyle \frac{\lambda }{d_{pore}}}={\displaystyle \frac{k_{B}T}{\sqrt{2}\pi {\varsigma }^{2}p{d}_{pore}}},$$

(13)

where $\lambda$ is the mean free path, $k_B$ is the Boltzmann constant and $\varsigma$ is the molecular diameter. A Knudsen number above unity describes a mass transfer in pores, which is mainly influenced by molecular wall interactions. A Knudsen number 0.01 < Kn < 1 is in the transition regime between viscous flow and Knudsen flow and for Knudsen numbers Kn≪ 1 a purely viscous flow can be assumed [17].

The Knudsen number is determined for all experiments. The values range between 23 < $Kn$ < 84 for the validation experiments, considering a pore size between 2.88 and 10.49 $\mathsf{\mu }$$\mathrm{m}$ at $-30$ ${}^{\circ }\mathrm{C}$ at a chamber pressure of $0.4$ $\mathrm{mbar}$.

Mass transfer of vapor through the open pores by Knudsen Diffusion

For the simple Knudsen flow approach in case of a transitional flow regime, the sublimation mass flow rate through the porous layer is calculated using Equations (14) and (15):

$$\begin{array}{cc}\hfill {\dot{m}}_{pore,Kn}^{″}& =-{\displaystyle \frac{D_{Kn}}{RT}}{\displaystyle \frac{\partial p_w}{\partial r}},\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill D_{Kn}& ={\displaystyle \frac{4}{3}}·d_{pore}{\displaystyle \frac{{ϵ}_p}{{\tau }_p}}·\sqrt{{\displaystyle \frac{RT}{2\pi M}}},\hfill \end{array}$$

(15)

where $D_{Kn}$ is the Knudsen diffusion coefficient, R is the universal gas constant, M is the molecular weight and $\frac{\partial p_w}{\partial r}$ the gradient of the partial pressure in radial direction.

Mass transfer of vapor through the open pores by Dusty gas Model

The dusty gas model was shown to predict the mass transfer in catalysts [18,19,20], and was recently used for the modeling of spray freeze-dried coffee particles on trays [5]. The mass flux ${\dot{m}}_{pore,Du}^{"}$ is provided by the following equation:

$${\dot{m}}_{pore,Du}^{"}=-{\displaystyle \frac{M_w}{R_g·T_d}}·\left(k_1+k_2·p_{w,d}\right){\displaystyle \frac{\partial p_{w,d}}{\partial r}},$$

(16)

where $k_1$ and $k_2$ are constants describing the bulk diffusivity and the self-diffusivity, respectively. These constants can be calculated by the following equations, where Equation (17) describes the Knudsen flow term and Equation (18) describes the viscous flow as follows:

$$\begin{array}{cc}\hfill k_1& ={\displaystyle \frac{C_2·D_{w,in}^0·K_w}{C_2·D_{w,in}^0+K_{mx}·P}},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill k_2& ={\displaystyle \frac{K_w·K_{in}}{C_2·D_{w,in}^0+K_{mx}·P}}+{\displaystyle \frac{C_{01}}{{\mu }_{mx}}}.\hfill \end{array}$$

(18)

$K_w$ is the Knudsen diffusivity of water and $K_{in}$ is the Knudsen diffusivity of the inert gas, which are calculated by Equations (19) and (20). The Knudsen diffusivity of the binary mixture $K_{mx}$ can be calculated by its molar fractions $y_w$ and $y_{in}$, as shown in Equation (21):

$$\begin{array}{cc}\hfill K_w& =C_1·{\left({\displaystyle \frac{R_g·T_I}{M_w}}\right)}^{0.5},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill K_{in}& =C_1·{\left({\displaystyle \frac{R_g·T_I}{M_{in}}}\right)}^{0.5},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill K_{mx}& =y_w·K_w+y_{in}·K_{in}.\hfill \end{array}$$

(21)

The viscosity of the vapor phase, which is needed to calculate the viscous flow term, is calculated as function of the temperature:

$${\mu }_{mx}=18.4858·{10}^{-7}·\left({\displaystyle \frac{T_I^{1.5}}{T_I+650}}\right).$$

(22)

$C_{01}$, $C_1$ and $C_2$ are constants for the Darcy flow permeability, the Knudsen flow permeability and the bulk diffusivity. Liapis and Bruttini [10] propose to calculate these constants as a function of structural parameters such as the porosity ${\epsilon }_p$, the tortuosity ${\tau }_p$ of the particle and the pore diameter $d_{pore}$:

$$\begin{array}{cc}\hfill C_1& ={\displaystyle \frac{{\epsilon }_p}{{\tau }_p}}·{\displaystyle \frac{48.5·d_{pore}}{R_g^{0.5}}},\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill C_2& ={\displaystyle \frac{{\epsilon }_p}{{\tau }_p}}.\hfill \end{array}$$

(24)

The expression for the Darcy flow permeability is adjusted in comparison to the literature used, since in this work no particle bed but a single particle is investigated. Therefore, the equation is taken from Kast [17] who gives an expression for the Darcy flow permeability in pores:

$$C_{01}={\displaystyle \frac{{\epsilon }_p}{{\tau }_p}}·{\displaystyle \frac{d_{pore}^2}{32}}.$$

(25)

The tortuosity ${\tau }_p$ is difficult to measure experimentally and is therefore assumed as ${\tau }_p$ = $\sqrt{2}$, as it was done by other authors for soluble coffee [5,10]. The vapor pressure at the interface as a function of temperature was already proposed by Marti and Mauersberger [21]. It is calculated by:

$$p_w^{sat}={10}^{\frac{-2663.5}{T_{fc}}+12.537}.$$

(26)

Mass transfer of vapor for second drying stage modeled by first-order reaction kinetic

The secondary drying is modeled by a first-order kinetic expression as already shown in previous attempts [5,10]. Therefore, the vapor flow is assumed to decrease with proceeding drying progress, which is characterized by the desorption constant $k_{des}$:

$${\dot{m}}_{des}=-k_{des}·c_{sw}.$$

(27)

2.2. Validation Experiments

For the validation of the model, freeze-drying experiments are carried out. For the experiments, the setup as described in [22] is used to structure the material. In addition to the described trials, the samples are gassed with nitrogen using a different overrun. The overrun is defined as the amount of gas $V_g$ in comparison to the amount of liquid $V_l$:

$$OR={\displaystyle \frac{V_g}{V_l}}={\displaystyle \frac{{\rho }_l-{\rho }_f}{{\rho }_f}}={\displaystyle \frac{V_{tot}·{\rho }_l-m_f}{m_f}},$$

(28)

where ${\rho }_l$ and ${\rho }_f$ are the liquid and foam density, respectively, $V_{tot}$ is the total volume and $m_f$ corresponds to the foam mass. To change the internal porous structure, two parameters are changed: (1) the amount of overrun is varied between 50%, 100% and 150% to generate foams with varying porosity and (2) the freezing method is changed. For freezing, the slurry is distributed in a 10 $\mathrm{m}$$\mathrm{m}$ layer on a circular aluminum tray with a diameter of 80 $\mathrm{m}$$\mathrm{m}$. Three methods are applied to achieve a change in ice crystal size and structure: Airblast freezing is used inside of a $-60$ ${}^{\circ }\mathrm{C}$ freezer, using a fan (3212 JH4, EBM-papst, Germany), which is placed 120 $\mathrm{m}$$\mathrm{m}$ above the sample. The second way of freezing is carried out by placing the sample on a steel plate inside the freezer at $-60$ ${}^{\circ }\mathrm{C}$, using conductive heat transfer on the bottom. For the third method, the sample is placed within of a polystyrene box with a wall thickness of 20 $\mathrm{m}$$\mathrm{m}$ into the freezer to reduce the heat transfer and allow for lower freezing rates. The scraping speed inside the SSHE was held constant at 300 rpm and the temperature of the coolant was at $-25$ ${}^{\circ }\mathrm{C}$ for all trials. After the freezing is completed, the frozen cake is milled and fractionated by sieving at $-60$ ${}^{\circ }\mathrm{C}$. Only particles between $2.24$ $\mathrm{m}$$\mathrm{m}$ and $2.80$ $\mathrm{m}$$\mathrm{m}$ are considered for the freeze-drying process. For the drying, the experimental setup is used as shown in [22]. A total of 1 $\mathrm{g}$ of the frozen granules is placed on a polymer mesh, with a mesh size of $0.5$ $\mathrm{m}$$\mathrm{m}$, hanging from a microbalance (Martin Christ Gefriertrocknungsanlagen, Osterode am Harz, Germany) to receive heat by radiation from the bottom and the top. The reduction in moisture is monitored over time with a resolution in time of 1 s and resolution in weight of 1 $\mathrm{m}$$\mathrm{g}$. The drying profile used is summarized in

Table 1. Drying profile for single-layer drying to determine the freeze-drying kinetics of instant coffee.

Step	T${}_{\mathit{start}}$	T${}_{\mathit{end}}$	P${}_{\mathit{chamber}}$	t
	(${}^{\circ }$C)	(${}^{\circ }$C)	(mbar)	(min)
1	−30	−30	1013	1
2	−30	46	0.4	38
3	46	46	0.4	90
4	46	25	0.4	20

.

For the analysis of the dried coffee particles, the Sauter mean diameter is calculated from the particle size distribution (PSD) to confirm the homogeneity in particle size of all samples. The PSD is measured (CAMSIZER XT, Retsch GmbH, Haan, Germany) using the particle size based on the measured area $x_{area}$. Particles smaller than 1 $\mathrm{m}$$\mathrm{m}$ are neglected for the Sauter mean diameter calculation, to avoid variance due to abrasion. The pore size distribution and the porosity of the samples are measured by mercury porosimetry (Pascal 140 + 440, Thermo Fisher Scientific Inc., Waltham, MA, USA). A total of $0.15$ $\mathrm{g}$ of the sample are analyzed at pressures between $0.1$ $\mathrm{k}$$\mathrm{Pa}$ and 400 $\mathrm{M}$$\mathrm{Pa}$ with an accuracy of 0.2%. The pressure ranges allow for the analysis of pore diameters of 3.8 to 116 $\mathsf{\mu }$$\mathrm{m}$ and 0.0036 to 15 $\mathsf{\mu }$$\mathrm{m}$ for the different units, respectively. For the determination of the closed pore porosity, a helium pycnometer (AccuPyc 1330, Micromeritics, Norcross, GA, USA) is used. For the precise method to calculate the pore diameter and the closed porosity, refer to [22]. Images of selected samples are taken by scanning electron microscopy to verify the obtained results for the pore structure (Gemini 1530, Carl Zeiss Meditec, Jena, Germany). The particles are sputter-coated with a 40 $\mathrm{n}$$\mathrm{m}$ gold layer to avoid interferences during measurements. Magnifications of 25, 100, 400 and 1000 times are performed with an electron beam of 5 $\mathrm{k}$$\mathrm{V}$.

3. Results

3.1. Experimental Results

To validate the model, experimental data are generated. To compare the drying kinetics and to define the transition between the end of the primary drying and the beginning of secondary drying, the time, which is needed to remove 80% of the water, is referred to as primary drying time. The residual water is removed by desorption during secondary drying. Since primary drying and secondary drying are not strictly separated, a part of the bound water is desorbing during the primary drying step as already discussed by other authors [4,5,8,9].

The structural parameters are summarized in

Table 2. Drying time and mean evaporation rate of coffee granules as a function of Sauter mean diameter d${}_{32,p}$ of the particle size distribution and structural parameters such as mean pore size d${}_{50,pore}$ and open porosity $ϵ$.

${\mathit{d}}_{32,\mathit{p}}$	${\mathit{d}}_{50,\mathit{pore}}$	${\mathit{ϵ}}_{\mathit{p}}$	${\mathit{t}}_{80}$	${\dot{\mathit{m}}}_{80,\mathit{mean}}$
(mm)	(µm)	( - )	(min)	(kg/h)
2.65	2.88	0.743	44.4 ± 2.2	2.36 ± 0.20
2.67	7.20	0.771	39.4 ± 0.1	3.03 ± 0.00
2.84	9.49	0.759	39.1 ± 0.5	3.05 ± 0.04
2.60	5.34	0.692	44.4 ± 1.2	3.56 ± 0.10
2.50	7.71	0.796	37.4 ± 1.6	2.47 ± 0.11

. For the characterization of the particle diameter, the Sauter mean diameter is chosen, considering the surface area with respect to the volume, since freeze-dried coffee particles are not spherical. The primary drying time $t_{80}$ and the averaged evaporation rate ${\dot{m}}_{80,mean}$ are used as measures for the description of the drying kinetics, where the latter is averaged over 15 s to overcome noice suppression.

In Figure 2, the internal pore structure of the coffee granules is shown for various interconnecting pore sizes, measured by mercury porosimetry. Here, the SEM images proof the range of interconnecting pore diameters, which are connecting the gas pores. The analysis of the images with the software ImageJ shows that gas bubble size can be held constant for all experiments at ($23.1\pm 4.5$) $\mathsf{\mu }$$\mathrm{m}$.

The results show that an increasing pore diameter leads to decreasing primary drying time. These results go along with the findings of Searles et al. [16], Hottot et al. [3] and our previous findings [22]. Mass transfer limitations seem to become more significant at small pore diameters of 3 $\mathsf{\mu }$$\mathrm{m}$ or lower, but have less strong effects at higher pore diameters of 7 to 10 $\mathsf{\mu }$$\mathrm{m}$.

For the increased porosity, a decreasing trend in drying time can be observed, but the interpretation of experimental data need to be treated with care. At constant sample mass and particle diameter, but increasing overrun, the number of particles is increasing. As only a single layer is observed, the particle surface area available for radiation is increasing, which improves the heat transfer as well. In consequence, an increase of overrun leads to a reduction of the total water amount per particle and increases the drying rates. Therefore, the mean evaporation rate per particle is more meaningful compared to the drying time and shows a decreasing rate with increasing overrun, since heat transfer limitations become dominant.

Due to increased robustness the experiments are conducted using a single layer of material instead of single particle investigations. However, for a single layer of particles, not the complete surface area is available for the heat transfer due to the slightly overlapping of the particles. These effects turn out to be more significant for samples with increasing porosity and are corrected by a fitting parameter f. As shown in previous publications [23,24], vials drying on the edges of the plate dry faster compared to the center vials, since a larger part of the vial surface is exposed to a radiative source. Additionally, the shelves are not the only sources for radiative heat transfer, but also the other surfaces inside the freeze-dryer and their reflection of heat radiated by the front window. Therefore, besides the shelf temperature also an average environmental temperature needs to be defined to adapt the model to experimental conditions. These overlapping effects and the heat source are relevant for the heat transfer and are considered in the model approach as fitting parameters for the radiative heat transfer:

$${\dot{Q}}_{rad}=fϵA_p\sigma ((g·T_{sh}^4+(1-g)T_{amb}^4)-T_p^4),$$

(29)

where f is the fitting parameter to consider the reduction in surface area for the heat transfer and g gives a ratio between heat transferred by the shelves in comparison to the heat transferred by other sources, where $T_{amb}$ is the ambient temperature of 20 ${}^{\circ }\mathrm{C}$.

3.2. Model Results

3.2.1. Validation

Both model approaches are validated with the experimental data. All samples for the validation process have an initial coffee concentration of 50% w/w and are foamed with nitrogen. The samples shown in Section 3.1 differ in their porosity and pore size. In Figure 3, the model results are compared with experimental data, varying the pore size diameters of the coffee granules between $2.88$ $\mathsf{\mu }$$\mathrm{m}$ and $9.49$ $\mathsf{\mu }$$\mathrm{m}$, keeping the applied overrun constant at 100%.

The total porosity (open and closed pore porosity) is comparable for all samples with 0.79 ± 0.02. The closed pore porosity for all samples ranges between 2 and 5%. The data show that the model can predict the loss of moisture over time, showing a slightly curved profile during primary drying, which is expected as shelf temperature is ramped up. The model and the experiment show that narrow pores of $2.88$ $\mathsf{\mu }$$\mathrm{m}$ decrease the drying rates significantly. Nevertheless, this limitation seems to decrease significantly for pore sizes between $7.2$ $\mathsf{\mu }$$\mathrm{m}$ and $9.5$ $\mathsf{\mu }$$\mathrm{m}$, where the drying process gets limited rather by heat than mass transfer.

For the evaporation rate, the regions of primary and secondary drying can be observed in model and experiment. During primary drying a constant increase of the evaporation rate is observed, since shelf temperature increases in time. When shelf temperature is constant, a decreased rate is shown with time due to progressing heat and mass transfer limitations. During secondary drying a decay function can be observed, which approximates to an equilibrium value. While in the model primary drying and secondary drying are strictly separated, it can be observed that for the experiments the two mechanisms are overlapping, smoothening the curve in the transition region. Faster drying samples with usually increased pore diameter show a higher evaporation rate during primary drying (green and blue). During a later time-step between 2000 and 2500 s, the water content is decreasing significantly and the red curve shows the highest evaporation rate since primary drying is still not finished. The influence of the porosity on the drying kinetics are shown in Figure 4. The pore size is held constant in a range of ($6.5\pm 1.18$) $\mathsf{\mu }$$\mathrm{m}$.

It shows that also for the varying porosity, the model fits the experimental data. With an increasing porosity, the drying time is decreasing and the vapor flow characterized by the normalized evaporation rate is increasing. The reason here is that for the experiments, the number of particles are increasing with increasing porosity and constant particle size and sample weight. Therefore, the heat transfer is increased. In the model, the volume of a particle is held constant and thus the surface area per particle stays constant. The amount of water is decreasing per particle. The normalization of the evaporation rate makes both approaches comparable and shows good trends for the primary drying, where low porosity shows the lowest normalized evaporation rate and high porosity shows the highest, as observed in the experiments. For secondary drying, all modeling approaches are in good agreement with the experimental data. In summary, the model shows similar trends in comparison to the experimental data for the Knudsen and the dusty gas approach, where only small deviations can be observed.

Table 3. Coefficient of determination of the sample moisture evolution over time during freeze-drying, using the Knudsen approach $R_{Kn}^2$ and the dusty gas approach $R_{Du}^2$ in comparison to experimental data.

${\mathit{d}}_{32,\mathit{p}}$	${\mathit{d}}_{50,\mathit{pore}}$	${\mathit{ϵ}}_{\mathit{p}}$	${\mathit{R}}_{\mathit{Kn}}^2$	${\mathit{R}}_{\mathit{Du}}^2$
(mm)	(µm)	( - )	( - )	( - )
2.65	2.88	0.743	0.995649	0.995648
2.67	7.20	0.771	0.995247	0.995255
2.84	9.49	0.759	0.996340	0.996339
2.60	5.34	0.692	0.996624	0.996621
2.50	7.71	0.796	0.994986	0.994981

shows coeffcient of determination for both model approaches. A deviation both models is observed in the sixth decimal place. Both models show that an increasing pore diameter can decrease the limitations in mass transfer during primary drying and lead therefore to reduced drying times. An increase in porosity also shows a reduction in the drying time. However, this also leads to a reduction of absolute water content per volume of the particle on top of the improved pore space for vapor transport. The benefits and drawbacks of both parameters are discussed in detail in the sensitivity analysis.

3.2.2. Sensitivity Analysis

For the sensitivity analysis, the model results are extrapolated and the influence of pore diameter and the porosity of the particle on the overall drying time and average sublimation rate is investigated. The pore diameter is varied from 0 to 20 $\mathsf{\mu }$$\mathrm{m}$ and the porosity is varied from 0.6 to 0.9. In Figure 5, the primary drying time as a function of both structural parameters is shown.

The vapor flow through the porous dried layer is modeled with Knudsen diffusion and with the dusty gas model. For both simulations, the same fitting parameters are used. Both approaches come to similar results in prediction of the drying time, with marginal differences. An increase of the pore diameter or an increase of porosity decreases the primary drying time. As mentioned before a change in porosity shows several overlapping effects: An increase of porosity (1) decreases the amount of water which needs to be evaporated for a single particle, (2) increases the pore volume for the mass transfer and (3) reduces the conductivity of the dried layer and inhibits the heat transfer inside the particle towards the sublimation front. These effects are coupled and therefore complicate an interpretation for the optimal configuration.

To understand the impact of the porosity on the heat and mass transfer mechanisms, the sublimation flow of water vapor needs to be considered to remove the influence of the decreasing water content, as shown in Figure 6. Here, only the dependency of the drying kinetics on pore volume and conductivity is considered.

Again both models show consistent results. It can be observed here that an increasing porosity of the sample reduces the vapor flow significantly. The optimum for both models is shown with 3.3 mgh⁻¹ at a pore size of 20 $\mathsf{\mu }$$\mathrm{m}$ and an accessible porosity of 0.61 which corresponds to a sample where no gas was injected. With this configuration the optimal amount of vapor can be removed from the system during primary drying. It shows that the influence of porosity is increasing with increasing pore diameter since the freeze-drying process gets increasingly limited by the heat transfer. On the other hand, it shows that an increasing pore diameter has higher influences at low porosities of 60% in comparison to high porosities where the system is anyways limited by the heat transfer. At porosities of almost $0.9$ a change in pore diameter from 10 to 20 $\mathsf{\mu }$$\mathrm{m}$ does not significantly change the results, while a change from 1 to 10 $\mathsf{\mu }$$\mathrm{m}$ has an impact. The highest limitations are observed at small pore diameters of 1 $\mathsf{\mu }$$\mathrm{m}$ and high porosities of almost $0.9$.

4. Conclusions

In this work, a mathematical model is presented to predict the freeze-drying time for bulky food materials, validated with experimental data. A Knudsen flow and a dusty gas model were used and compared for the mass transfer through the dried layer. Both model approaches are able to fit the experimental results in terms of moisture loss and sublimation rate over time during the freeze-drying process. The model shows room for improvement in the prediction of the transition regime between primary and secondary drying, where experiments have shown that a significant amount of water is already desorbed, while primary drying is still not completed. Nevertheless, this fact does not show any influence on the correct prediction of the drying kinetics of the model in comparison to the experimental approach. As it can be shown in the sensitivity analysis of this study, the pore size might be the appropriate tool to decrease the mass transfer limitations during freeze-drying. An increase of porosity also shows a reduction of mass transfer resistances, but on the expense of reduced heat transfer through the dried layer of the particle. Therefore, an increase of porosity leads to increasing drying time in the observed range of parameters at the given 50% w/w coffee concentration. Furthermore, a change in porosity affects the density and fragility of the product which can impact the product quality.

Both model approaches, the Knudsen flow and the dusty gas approach, are comparable tools to model the freeze-drying of food granules in a range of pore sizes of 1 to 20 $\mathsf{\mu }$$\mathrm{m}$ and porosities in the range of 0.6 to 0.9, which are typical for freeze-dried coffee products. In terms of complexity, an exclusive use of Knudsen flow is shown to be sufficient in this pore size range, shown to have a comparable coefficient of determination, where the dusty gas model may be of higher interest for bigger pore diameters or higher pressures, where other transport mechanisms such as viscous flow are dominating.