Freeze-Drying

Freeze-drying, or lyophilization, is a dehydration process that is often used to preserve a perishable food or material. It is also frequently used as a way to remove water from goods, making them lighter and easier to transport. Lyophilization is widely recognized as an important technique in many industries. It is used, for example, in the pharmaceutical industry for the preservation of antibiotics, in the manufacturing of semiconductor ceramics, and by preservation initiatives, such as in the restoration of water-damaged documents. And of course, the process is widely used in the food industry, to preserve tasty snacks that can last up to 30 years. The wet substance is frozen and ice (or some other frozen solvent) is removed through sublimation in the presence of a high vacuum.

Freeze-drying relies on sublimation, where a frozen liquid undergoes phase change directly from a frozen state and into a gaseous state. As can be shown in a phase diagram, at very low pressures and temperatures, a solid can pass directly into the gaseous stage without passing through the intermediary liquid phase.

This model demonstrates the process of ice sublimation in a vial under vacuum-chamber conditions, a test case for many freeze-drying setups.

Model Definition

In this example, you model the process of sublimation of ice through the pores of a porous medium in a vial. Figure 1 depicts the model geometry cut in half.

Figure 1: Vial model geometry.

The geometry is axisymmetric but the boundary condition presented in Equation 1 below is not. For this reason, the model must be built in 3D. Because y = 0 is a symmetry plane, only half of the structure is represented to save computational cost.

Mass Transfer

The frozen product (here skim milk) initially fills 98% of the vial.The remaining gap contains the dried product and the vapor generated by sublimation, and expands during the sublimation process.The inert gas is neglected in the simulation.

The vapor flow rate Qm through the pores of the dried product is defined by the stationary Darcy’s law:

(1)

where ρv is the vapor density, μv is the vapor dynamic viscosity, pv is the vapor partial pressure, and κ is the permeability of the dried product.

The top of the vial is close to vacuum, maintained at a vapor pressure p0 of 24 Pa.

At the sublimation interface, the mass flux N0 is:

(2)

where ρice is the ice density, ε is the porosity of the product, and Vs is the sublimation interface velocity.

Heat Transfer

You solve the heat transfer equation for porous media, without convection in the frozen domain,

(3)

and with convection in the dried domain

(4)

Here, (ρCp)eff,dr and keff,dr refer to the equivalent volumetric heat capacity and thermal conductivity of the dried region, (ρCp)eff,fr and keff,fr refer to the equivalent volumetric heat capacity and thermal conductivity of the frozen region, Cp,v is the specific heat capacity of vapor, and u is the velocity obtained from Darcy’s law.

In the vapor domain, the density is calculated using the molar mass of vapor, Mv, and the ideal gas law:

(5)

where R is the universal gas constant.

Convective exchanges with the surrounding air occur all around the vial. At the bottom of the vial, convective exchanges with the shelf occur as well, and the heat transfer coefficient is lower on the center than on the exterior of the vial due to an air gap.

Ice-Vapor Interface

Assuming thermodynamic equilibrium at the phase change interface, the sublimation front temperature Ts is defined from the vapor pressure pv at this boundary, using the Clausius–Clapeyron relation:

(6)

where Ls is the sublimation latent heat.

Simultaneous heat and mass balances at the sublimation interface lead to the Stefan condition for the interface velocity Vs:

(7)

where Qs is the jump in the normal heat flux at the interface.

Results

The process of primary drying takes about 10 hours to transform most of the initial amount of ice. The frozen product has the proportions shown in Figure 2.

Figure 2: Temperature and heat flux at the end of the drying time period.

You can see from this visualization that the sublimation front has a concave curved shape that is lower around the outside edges of the vial wall. As the ice becomes a vapor, the mass of the solid decreases and therefore the mesh in the frozen domain squeezes as well.

Notes About the COMSOL Implementation

To set up the application in COMSOL Multiphysics, use the Deformed Geometry interface to track the ice surface, and then compute coupled mass and heat balances on the moving mesh. Figure 3 shows the initial mesh and Figure 4 shows the deformed mesh at the end of the simulation.

Figure 3: Initial swept mesh.

Figure 4: Deformed mesh at the end of the simulation.

In Equation 7, Qs corresponds to the jump in the normal heat flux at the interface. This quantity can be precisely evaluated through the Lagrange multiplier for temperature, T_lm. This variable is computed in the feature, which sets a fixed temperature weak constraint on the sublimation front.

References

1. A. Alonso and others, “Time-scale modeling and optimal control of freeze-drying,” Journal of Food Engineering, vol. 111, pp. 655–666, 2012.

2. W.J. Mascarenhas, H.U. Akay, and M.J. Pikal, “A computational model for finite element analysis of the freeze-drying process,” Computer Methods on Applied Mechanics and Engineering, vol. 148, pp. 105–124, 1997.

3. Hriberšek and others, “Lyophilization model of mannitol water solution in a laboratory scale lyophilizer”, Journal of Drug Delivery Science and Technology, vol. 45, pp. 28–38, 2018.

Heat_Transfer_Module/Phase_Change/freeze_drying

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

In the tree, select .

Click .

In the tree, select .

Click .

Click

In the tree, select .

Click

For this model, some parameters are needed. Start by loading all of them from a text file.

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file freeze_drying_parameters.txt.

Add empty materials that will be linked later to the porous materials.

Vapor

In the window, right-click and choose .

In the window for , type Vapor in the text field.

Ice

Right-click and choose .

In the window for , type Ice in the text field.

Product (Skim Milk)

Right-click and choose .

In the window for , type Product (Skim Milk) in the text field.

Geometry 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Cylinder 1 (cyl1)

In the toolbar, click

In the window for , locate the section.

In the text field, type R0.

In the text field, type Z0.

Click to expand the section. In the table, enter the following settings:


Layer name	Thickness (mm)
Layer 1	Zi

Clear the check box.

Select the check box.

Cylinder 2 (cyl2)

In the toolbar, click

In the window for , locate the section.

In the text field, type R0.

In the text field, type Z0.

Locate the section. In the table, enter the following settings:


Layer name	Thickness (mm)
Layer 1	R0/2

Block 1 (blk1)

In the toolbar, click

In the window for , locate the section.

In the text field, type R0*5.

In the text field, type R0*2.

In the text field, type Z0.

Locate the section. In the text field, type -R0*2.

In the text field, type -R0*2.

Difference 1 (dif1)

In the toolbar, click

and choose .

In the window for , locate the section.

Click

In the dialog box, type cyl1 cyl2 in the text field.

Click .

In the window for , locate the section.

Find the subsection. Select the

toggle button.

Select the object only.

Click

Click the

button in the toolbar.

These steps create a geometry similar to that in Figure 1.

Definitions

Interface

In the window, expand the node.

Right-click and choose .

Define the sublimation temperature at the phase change interface according to the Clausius-Clapeyron equation.

In the window for , type Interface in the text field.

Locate the section. In the table, enter the following settings:


Name	Expression	Unit	Description
T_s	2.19e-3DelHs1[kg/J]/(28.89-log(p1[1/Pa]))1[K]	K	Interface temperature

Integration 1 (intop1)

In the toolbar, click

and choose .

In the window for , locate the section.

Click

In the dialog box, type 1 3 5 in the text field.

Click .

Minimum 1 (minop1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Click

In the dialog box, type 6 13 20 in the text field.

Click .

Add porous materials to the materials node, but do not define their properties at this point. After setting up the physics interface, COMSOL Multiphysics automatically detects which material properties are needed.

Materials

Porous Material 1 (poromat1)

In the window, under right-click and choose .

In the window for , locate the section.

In the list, choose , , and .

Click

Select Domains 2, 4, and 6 only.

Fluid 1 (poromat1.fluid1)

Right-click and choose .

Solid 1 (poromat1.solid1)

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

In the θs text field, type 1-por_p.

Porous Material 2 (poromat2)

In the window, right-click and choose .

Select Domains 1, 3, and 5 only.

Fluid 1 (poromat2.fluid1)

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Solid 1 (poromat2.solid1)

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

In the θs text field, type 1-por_p.

Now, set up the physics interfaces. Start with the Darcy’s law interface, that simulate the vapor flow through the pores of the dried layer.

Click the

button in the toolbar.

In the dialog box, select in the tree.

In the tree, select the check box for the node .

Click .

Darcy’s Law (dl)

In the window, under click .

Select Domains 2, 4, and 6 only.

In the window for , click to expand the section.

From the list, choose .

Locate the section. In the pref text field, type 0.

Click to expand the section. Select the check box.

Fluid and Matrix Properties 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Initial Values 1

In the window, click .

In the window for , locate the section.

In the p text field, type Pc.

Pressure 1

In the toolbar, click

and choose .

Select Boundaries 7, 14, and 21 only.

In the window for , locate the section.

In the p0 text field, type Pc.

Mass Flux 1

In the toolbar, click

and choose .

Select Boundaries 6, 13, and 20 only.

In the window for , locate the section.

In the N0 text field, type -ht.pci1.vn*rho_ice*por_p.

Symmetry 1

In the toolbar, click

and choose .

Select Boundaries 4, 12, and 25 only.

Heat Transfer in Porous Media (ht)

Now, set up the Heat Transfer in Porous Media interface, in both dried and frozen layers.

Dried Layer

In the window, under right-click and choose .

In the dialog box, type Dried Layer in the text field.

Click .

Fluid 1

In the window, click .

In the window for , locate the section.

From the u list, choose .

Porous Matrix 1

In the window, click .

In the window for , locate the section.

From the list, choose .

Initial Values 1

In the window, click .

In the window for , locate the section.

In the T text field, type Ti.

Frozen Layer

In the toolbar, click

and choose .

In the window for , type Frozen Layer in the text field.

Select Domains 1, 3, and 5 only.

Porous Matrix 1

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

Ambiant Heat Flux

In the toolbar, click

and choose .

In the window for , type Ambiant Heat Flux in the text field.

Locate the section. Click the button.

In the h text field, type 3.6.

In the Text text field, type Ta.

Locate the section. Click

In the dialog box, type 2 5 7 14 21 22 23 in the text field.

Click .

Surface-to-Ambient Radiation 1

In the toolbar, click

and choose .

Select Boundaries 7, 14, and 21 only.

In the window for , locate the section.

From the ε list, choose . In the associated text field, type 0.9.

In the Tamb text field, type Ta.

Shelf Heat Flux (Center)

In the toolbar, click

and choose .

In the window for , type Shelf Heat Flux (Center) in the text field.

Locate the section. Click

In the dialog box, type 10 in the text field.

Click .

In the window for , locate the section.

Click the button.

In the h text field, type 11.

In the Text text field, type Ts.

Shelf Heat Flux (Exterior)

In the toolbar, click

and choose .

In the window for , type Shelf Heat Flux (Exterior) in the text field.

Locate the section. Click

In the dialog box, type 3 17 in the text field.

Click .

In the window for , locate the section.

Click the button.

In the h text field, type 62.3.

In the Text text field, type Ts.

Symmetry 1

In the toolbar, click

and choose .

Select Boundaries 1, 4, 9, 12, 24, and 25 only.

Phase Change Interface 1

In the toolbar, click

and choose .

In the window for , locate the section.

Click

In the dialog box, type 6 13 20 in the text field.

Click .

In the window for , locate the section.

In the Tpc text field, type T_s.

In the Ls → f text field, type DelHs.

From the list, choose .

Deformed Geometry (dg)

Finish with the Deformed Geometry interface under the Definitions node. Roller boundaries are set up along the vial and the symmetry axis. The feature in the Heat transfer interface already specifies a prescribed normal mesh velocity.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Free Deformation 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Prescribed Mesh Displacement 2

In the toolbar, click

and choose .

In the window for , locate the section.

Click

Select Boundaries 1, 2, 4, 5, 9, 12, and 22–25 only.

Locate the section. Clear the check box.

Global Definitions

The material properties can now be specified.

Vapor (mat1)

In the window, expand the node, then click .

In the window for , locate the section.

Click

In the dialog box, type pressure in the text field.

Click

In the tree, select .

Click .

In the window for , locate the section.

Click

In the dialog box, type temperature in the text field.

Click

In the tree, select .

Click .

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	M_v*pA/R_const/T	kg/m³	Basic
Dynamic viscosity	mu	mu_v	Pa·s	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	k_v	W/(m·K)	Basic
Heat capacity at constant pressure	Cp	Cp_v	J/(kg·K)	Basic
Ratio of specific heats	gamma	1.33	1	Basic

Ice (mat2)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	k_ice	W/(m·K)	Basic
Density	rho	rho_ice	kg/m³	Basic
Heat capacity at constant pressure	Cp	Cp_ice	J/(kg·K)	Basic
Ratio of specific heats	gamma	1	1	Basic

Product (Skim Milk) (mat3)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Porosity	epsilon	por_p	1	Basic
Permeability	kappa_iso ; kappaii = kappa_iso, kappaij = 0	kappa_p	m²	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	k_p	W/(m·K)	Basic
Density	rho	rho_p	kg/m³	Basic
Heat capacity at constant pressure	Cp	Cp_p	J/(kg·K)	Basic