COMSOL 6.4 - Drying of a Potato Sample

Drying of porous media is an important process in the food and paper industry among others. Many physical effects must be considered: fluid flow, heat transfer with phase change, and transport of participating liquids and gases. All of these effects are strongly coupled, and predefined interfaces can be used to model these effects in an hygroscopic porous medium with COMSOL Multiphysics.

In this model, the liquid and gaseous phases are assumed to be in equilibrium inside a potato sample modeled as a porous medium, and the changing water saturation is computed over time, to model the heat and moisture transport by a two-phase flow.

Model Definition

This model describes a laminar dry airflow through a potato sample, modeled as a porous medium containing moist air and liquid water. The geometry and principle is shown in Figure 1.

Figure 1: Geometry and principle of the model.

Heat and Moisture Transport in Air

In the moist air domain, it is assumed that no liquid water is present. This means that all the moisture leaves the porous medium under vapor state.

The transport and thermodynamic properties of air with water vapor can be described using mixture laws, based on the amount of water vapor and dry air. This is done automatically in the Moist Air feature and subfeature in the Heat Transfer interface. See the Heat Transfer Module User’s Guide for the governing equations. As the input term for water vapor, the relative humidity

from the moisture transport equation is used in the Heat Transfer and Fluid Flow interfaces.

The compressible Navier–Stokes equations are solved to compute the velocity and pressure fields of the free flow for the convective transport of heat and moisture.

Another relevant effect for the moisture transport is the binary diffusion of vapor and air. The variations of the density of moist air induced by vapor transport are relatively small in the current configuration, so the diluted species formulation is used.

Two phase flow in The Potato Sample

The basic principle of modeling two phase flow in porous media is similar for many applications. First, to account for different phases, saturation variables are used that fulfill the following constraint:

(1)

where the index g is used for the gaseous phase (which in this model is moist air) and the index l is used for the liquid phase (which in this model is water).

Single-phase flow in porous media is described by the Brinkman equations. With an additional liquid phase, capillary effects also arise and the liquid flow is driven by a pressure gradient and the capillary pressure pc = pg − pl. How to deal with the latter depends on the application: sometimes the capillary pressure can be neglected, sometimes it is the driving effect and different approaches exist.

In this model, the capillary effects are treated by an additional diffusion term in the transport equation. The Brinkman equations are used to calculate the flow field ug and pressure distribution pg of moist air in the porous medium. Therefore the porosity εp must take into account that only a fraction of the void space is occupied by the gas phase.

The liquid-phase velocity is small compared to the moist air velocity and such that Darcy’s law is defined in terms of the gas-phase pressure gradient to calculate the water velocity ul according to

(2)

where κl and μl are the permeability and viscosity of the liquid phase.

Finally, due to the dimensions of the porous medium, the gravity effects on transport are neglected in both phases.

Moisture Transport in The Potato Sample

Following the ideas about mechanistic formulations from (Ref. 1), and by summing the mass conservation equations for vapor and liquid water, a single equation can be written to calculate the time-dependent evolution of moisture content, w(ϕw), in the porous medium:

(3)

The total moisture content accounts for the vapor and liquid water in the pores of the medium:

with ρl the liquid water density, ρg the moist air density, and ωv the vapor mass fraction defined as

The remaining terms in Equation 3 account for:

•

Convection of vapor due to total pressure gradient, with the flow field ug computed by the Brinkman equation.

•

Binary diffusion of water vapor and dry air in the gaseous phase, with gw defined as

A common correlation for the effective diffusivity Deff in an unsaturated medium is the Millington and Quirk equation:

with the vapor-air diffusivity Dva = 2.6⋅10−5 m2/s.

•

Convection of liquid water due to total pressure gradient, with the flow field ul computed by the Darcy’s Law in Equation 2. The permeability for the liquid phase κl depends on the overall permeability of the porous matrix κ and a relative permeability (see Permeability) κrl.

•

Capillary transport of liquid water, with glc defined as

with the diffusion coefficient Dw defined as (Ref. 2)

with ρs the solid phase density.

When summing the conservation equations to obtain Equation 3, the condensation source term in the mass conservation equation for liquid water cancels out with the evaporation source term in the mass conservation equation for vapor. However, the evaporation source can be expressed as follows:

Heat Transfer

Inside the porous domain the overall velocity field for liquid and gaseous phase contributes to the heat convection term. It is possible to account for the liquid water and moist air phases, by using the liquid saturation calculated by the moisture transport equation.

The mixture law is applicable in the gaseous phase, and averaged thermal properties (ρCp)eff and keff are defined by taking into account the properties of the porous matrix (ρs, Cp,s, ks), moist air (ρg, Cp,g, kg), and liquid water (ρl, Cp,l, kl):

(4)

Then the overall velocity can be expressed as the average of moist air and liquid water velocity, which is

(5)

The diffusive flux of enthalpy, due to vapor diffusion in air, and the capillary flux, due to the presence of the liquid water in the pores, are included in the following heat source term:

(6)

The heat of evaporation is inserted as a source term in the heat transfer equation according to

(7)

where Lv (J/mol) is the latent heat of evaporation.

Permeability

The permeability of the porous matrix κ defines the absolute permeability. When two phases are present, the permeability of each phase depends also on the saturation. This is defined by the relative permeabilities κrl and κrg for liquid and gaseous phase respectively, so that κl = κκrl and κg = κκrg. The determination of relative permeability curves is often done empirically or experimentally and the form strongly depends on the porous material properties and the liquids themselves. The functions that are used in this model (Ref. 2) are defined such that they are always positive:

The variable sli is the irreducible liquid phase saturation, describing the saturation of the liquid phase that will remain inside the porous medium.

Results and Discussion

Inside the potato sample, the relative humidity decreases from about 100% everywhere at the beginning of the simulation to about 10% at the end, which is the ambient air condition (Figure 2 and Figure 3).

Figure 2: Relative humidity after 5000 s.

At the beginning of the simulation (time 5000 s), the vapor concentration is close to its saturation value in the moist air enclosed in the pores. It goes down to the ambient concentration toward the end of the simulation.

Figure 3: Relative humidity after 40,000 s.

Because moisture is composed of vapor and liquid water in the porous medium, a good quantifier of the drying process is the liquid saturation. Figure 4 shows the drying front in the porous medium at time 5000 s, starting from the top-left corner, which is facing the dry air flow.

Figure 4: Liquid saturation after 5000 s.

The temperature plot (Figure 5) shows the corresponding cooling in the potato sample at time 5000 s, due to evaporation of the liquid water close to the drying front.

Figure 5: Temperature field after 5000 s.

Finally, the evolution of the moisture content in the potato sample over time is shown below.

Figure 6: Moisture content in the potato sample over time.

Notes About the COMSOL Implementation

Using a proper mesh size is important to resolve the steep gradients at the interface boundaries. Therefore the default mesh is refined.

To get good convergence of the time dependent behavior, first solve the stationary flow equations only. This solution will then be used for the time dependent study step. This approximation neglects the evaporation mass source in the fluid flow computation.

References

1. A.K. Datta, “Porous media approaches to studying simultaneous heat and mass transfer in food processes. I: Problem formulations,” J. Food Eng., vol. 80, 2007.

2. A.K. Datta, “Porous media approaches to studying simultaneous heat and mass transfer in food processes. II: Property data and representative results,” J. Food Eng., vol. 80, 2007.

3. A. Halder, A. Dhall, and A.K. Datta, “Modeling Transport in Porous Media with Phase Change: Applications to Food Processing,” J. Heat Transfer, vol. 133, no. 3, 2011.

Heat_Transfer_Module/Phase_Change/potato_drying

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, 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 potato_drying_parameters.txt.

Define the relative permeabilities for water vapor and liquid water as functions of liquid water saturation. The advantage of using functions is that they can be plotted immediately.

Definitions

Relative Permeability, Moist Air

In the window, expand the > node.

Right-click and choose > .

In the window for , type Relative Permeability, Moist Air in the text field.

In the text field, type kappa_rma.

Locate the section. In the text field, type S_l.

Find the subsection. In the table, enter the following settings:


Start	End	Function
0	1/1.1	1-1.1*S_l
1/1.1	1	eps

To force the saturation to be always positive eps is used as minimum value.

Click

Relative Permeability, Liquid Phase

In the toolbar, click

In the window for , type Relative Permeability, Liquid Phase in the text field.

In the text field, type kappa_rl.

Locate the section. In the text field, type S_l.

Find the subsection. In the table, enter the following settings:


Start	End	Function
0	S_il	eps
S_il	1	((S_l-S_il)/(1-S_il))^3

Click

Define the sorption isotherm by interpolation of tabulated data.

Sorption Isotherm

In the toolbar, click

In the window for , type Sorption Isotherm in the text field.

Locate the section. In the text field, type wc_int.

Click

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

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


Argument	Unit
t	1

In the table, enter the following settings:


Function	Unit
wc_int	kg/m^3

Click

Now, define the geometry.

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.15.

In the text field, type 0.05.

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.04.

In the text field, type 0.005.

Locate the section. In the text field, type 0.04.

Fillet 1 (fil1)

In the toolbar, click

On the object , select Points 3 and 4 only.

It might be easier to select the points by using the window. To open this window, in the toolbar click and choose . (If you are running the cross-platform desktop, you find in the main menu.)

In the window for , locate the section.

In the text field, type 2e-3.

In the toolbar, click

Click the

button in the toolbar.

Definitions

Define the ambient conditions used later in the domain and boundary conditions.

Ambient Properties 1 (ampr1)

In the toolbar, click

and choose .

In the window for , locate the section.

In the Tamb text field, type T0.

In the ϕamb text field, type phi_0.

Heat Transfer in Moist Air (ht)

Define the domain and boundary conditions for each interface. Start with the interface, by setting the initial and boundary conditions, and by adding a domain condition.

Initial Values 1

Use the ambient temperature defined previously as input.

In the window, under > click .

In the window for , locate the section.

From the T list, choose .

Moist Porous Medium 1

In the toolbar, click

and choose .

Select Domain 2 only.

Porous Matrix 1

In the window, click .

In the window for , locate the section.

From the list, choose .

Inflow 1

In the toolbar, click

and choose .

Select Boundary 1 only.

In the window for , locate the section.

From the Tustr list, choose .

Outflow 1

In the toolbar, click

and choose .

Select Boundary 9 only.

Laminar Flow (spf)

Continue with the interface.

In the window, under click .

In the window for , locate the section.

From the list, choose .

The next step enables additional features for the flow interface to account for porous domains also.

Select the checkbox.

Porous Medium 1

In the toolbar, click

and choose .

Select Domain 2 only.

Inlet 1

In the toolbar, click

and choose .

Select Boundary 1 only.

In the window for , locate the section.

From the list, choose .

Locate the section. In the Uav text field, type u0.

Outlet 1

In the toolbar, click

and choose .

Select Boundary 9 only.

Moisture Transport in Air (mt)

Finally, set the domain, initial, and boundary conditions for the interface.

Initial Values 1

In the window, under > click .

In the window for , locate the section.

From the ϕw,0 list, choose .

Hygroscopic Porous Medium 1

In the toolbar, click

and choose .

Select Domain 2 only.

Liquid Water 1

In the window, click .

In the window for , locate the section.

From the list, choose .

In the κrl text field, type kappa_rl(mt.sl).

Moist Air 1

In the window, click .

In the window for , locate the section.

In the κrg text field, type kappa_rma(mt.sl).

Initial Values 2

In the toolbar, click

and choose .

Select Domain 2 only.

In the window for , locate the section.

In the ϕw,0 text field, type phi_1.

Inflow 1

In the toolbar, click

and choose .

Select Boundary 1 only.

In the window for , locate the section.

From the Tustr list, choose .

From the ϕw,ustr list, choose .

Outflow 1

In the toolbar, click

and choose .

Select Boundary 9 only.

Porous Interface 1

In the toolbar, click

and choose .

Select Boundaries 4, 6, 7, 10, and 11 only.

Free-Porous Interface 1

In the window, click .

In the window for , locate the section.

From the list, choose .

Materials

Next, define the hygrothermal properties for the porous medium, based on measured data for potato.

Potato

In the window, under right-click and choose > .

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

Select Domain 2 only.

Add the subnode and define its material properties.

Locate the section. Click

Solid 1 (pmat1.solid1)

In the window, click .

In the window for , locate the section.

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

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


Property	Variable	Value	Unit	Property group
Density	rho	rho_s	kg/m³	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	k_s	W/(m·K)	Basic
Heat capacity at constant pressure	Cp	cp_s	J/(kg·K)	Basic
Porosity	epsilon	1	1	Porous model

Potato (pmat1)

Set the remaining material properties in the section.

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Diffusion coefficient	D_iso ; Dii = D_iso, Dij = 0	1e-8[m^2/s]exp(-2.8+2w_c/rho_s/(1-por))	m²/s	Basic
Water content	w_c	wc_int(phi)	kg/m³	Basic
Permeability	kappa_iso ; kappaii = kappa_iso, kappaij = 0	kappa	m²	Basic
Porosity	epsilon	1	1	Porous model

As the water content is defined as a function of the relative humidity, add the relative humidity as a model input required by the material.

In the window, click .

In the window for , locate the section.

Click

In the dialog, type relative in the text field.

In the tree, select > .

Click .

Definitions (comp1)

Define an integration operator in the porous domain, in order to check the moisture content evolution.

Integration Porous Medium

In the toolbar, click

and choose .

Select Domain 2 only.

In the window for , type Integration Porous Medium in the text field.

In the text field, type intopPorous.

To resolve all effects and get a better convergence for this highly nonlinear problem, refine the mesh to resolve the interface properly.

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Click

Study 1

Add a stationary step that solves for the fluid flow only. This approximation does not account for the evaporation mass source in the porous medium for the fluid flow initialization. The solution is used for the computation of the time dependent solution of all the processes.

Step 2: Stationary

In the toolbar, click

Drag and drop above .

In the window for , locate the section.

In the column of the table, under , clear the checkboxes for and .

Step 2: Time Dependent

In the window, click .

In the window for , locate the section.

In the text field, type range(0,10,90) range(100,100,900) range(1000,1000,40000).

Reduce the solver tolerance to enhance accuracy on mass balance.

From the list, choose .

In the text field, type 0.001.

In the toolbar, click

Results

The default plots show the velocity, the pressure, the relative humidity (Figure 3), the temperature field and the temperature with fluid flow.

Relative Humidity (mt)

Change the plot time to 5000 s for the relative humidity and the temperature to obtain Figure 2 and Figure 5.

In the window, click .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Temperature (ht)

In the window, click .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Result Templates

Insert one of the plots from the to see the distribution of the liquid saturation in the porous medium at time 5000 s as in Figure 4, in order to check the evolution of the drying front.

In the toolbar, click