Abstract

The modeling and numerical simulation of drying in porous media is discussed in this work by revisiting the different models of moisture migration during the drying process of porous media as well as their restrictions and applications. Among the models and theories, we consider those are ranging from simple ones like the diffusion theory to more complex ones like the receding front theory, the model of Philip and de Vries, Luikov’s theory, Krischer’s theory, and finally Whitaker’s model, in which all mass, heat transport, and phase change (evaporation) are taken into account. The review of drying models as such serves as the basis for the development of a framework for numerical simulation. In order to demonstrate this, the system of equations governing the drying process in porous media resulting from Whitaker’s model is presented and used in our numerical implementation. A numerical simulation of drying is presented and discussed to show the capability of the implementation.

1. Introduction

The modeling and numerical simulation of drying in porous media is a topic of great interest in process engineering and has attracted the attention of research institutions for decades. The analysis of absorption or desorption of liquid in general, and of water in particular, finds its application not only in chemical engineering, food processing, and pharmaceuticals but also in electronic packaging, which is becoming more and more important in the area of the Internet of Things [18]. Over the years, many works were undertaken to understand the mass and heat transport phenomena during drying. In this work, we will revisit the different models of moisture migration during the drying of porous media. The restrictions and applications of these models will be reviewed and discussed.

Drying models were developed since the beginning of the twentieth century. Different methods and numerical solutions were presented and applied successfully for different porous media. A review of these methods until the 1980s was given by Fortes and Okos [9] and Bories [10]. Reviews of empirical, analytical, and numerical methods of drying until the 1990s can be found in the work of Tsotsas [11]. In what follows, we will examine the most relevant theories and models for drying of porous media and discuss the most suitable model that can today be used for numerical simulation. Amongst others, we will first consider the so-called diffusion theory, the receding front theory, the model of Philip and de Vries, Luikov’s theory, and Krischer’s theory. In particular, we will discuss Whitaker’s model for drying and discuss the framework based on this model for numerical simulation of the drying process. Finally, we will present a numerical simulation using our implementation of Whitaker’s model and discuss the sample drying characteristics in our numerical study. For details of the implementation of the model, refer [12].

2. The Different Theories of Drying in Porous Media

2.1. Diffusion Theory

Lewis and Sherwood are known as pioneers in developing mathematical drying models by applying the Fourier equation of heat conduction to the drying of solids. In this equation, temperature and thermal diffusivity were replaced by moisture and moisture diffusivity, respectively. Starting from the idea of Lewis, Sherwood [13] provided solutions of the diffusion equation. Sherwood showed that the moisture transport involves two independent processes: the evaporation of moisture at the solid surface and the internal diffusion of liquid to the surface. The following simple diffusion model, in which the diffusivity of liquid is constant, was used to calculate the moisture distribution in a solid during drying and compared with experimental data of some materials (e.g., slabs of wood, clay, and soap):where is vaguely defined as the moisture content, represents time, and can be considered as an effective diffusion coefficient and is determined experimentally.

At the end of the 1930s, Ceaglske and Hougen [14] and Hougen et al. [15] pointed out that the moisture distribution cannot be calculated correctly only from (1). It was noted that the moisture movement in a solid during drying is not only due to diffusion but also due to other mechanisms such as gravity, external pressure, capillarity, convection, and vaporization-condensation when a temperature gradient is applied. Experimental data were collected for different kinds of porous media (clay, paper pulp, sand, lead shot, porous brick, and wood) and compared with the numerical results obtained by Sherwood’s model to show the limitation of his model and prove their criticism.

By using the capillary theory to describe the drying of granular materials (such as coarse, medium, and fine sand) and based on the collected experimental data, Ceaglske and Hougen [14] suggested that the effective diffusion coefficient should be considered as varying during drying and proposed the following diffusion model:

The effective diffusion coefficient is now considered as a function of moisture content, temperature, material type, and drying history. In solving the diffusion equation, this parameter is usually taken as constant or in the form of linear, exponential, or polynomial functions of moisture content. One example of these functions was given by Suzuki and Maeda [16]:where , , and are the constant factors. Suzuki and Maeda also presented an approximation method to describe the moisture distribution within drying porous materials in which the effective diffusion coefficient is expressed as an exponential function of moisture content:where and are the constant factors as in (3). This is a common function used to describe the effective diffusion coefficient. Suzuki and Maeda then solved the nonlinear diffusion problem by using dimensionless variables for the diffusion coefficient, moisture content, time, and space. It was shown that the steady-flux model (pseudosteady state) is fairly accurate when it is employed in low moisture content drying. However, it was also found that the diffusion model leads to a noticeable error when it is applied to high moisture content drying. The reason is that, in this case, capillary pumping has strong effect, and this mechanism must be taken into account.

In order to determine the effective diffusion coefficient, in the past decades, a large number of works were carried out. For example, the works of Saravacos and Raouzeos [17], Jaros et al. [18], Mourad et al. [19], Ribeiro et al. [20], Li and Kobayashi [21], and Srikiatden and Roberts [22] focused on materials in the food industry and agriculture. Koponen [23] dealt with wood, Ketelaars et al. [24] studied clay, and Pel et al. [25] considered fired-clay brick, sand-lime brick, and gypsum. In developing the diffusion model, heat transfer was included in examining the problem, for example, in the work of Thijssen and Coumans [26]. A shortcut method was proposed to calculate drying rates in nonisothermal drying of particles and hollow spheres. In their work, heat transfer was taken into account in the calculation of the evaporation rate and the sample temperature. Shrinking and nonshrinking materials were considered. The method proposed by Thijssen and Coumans is based on the numerical solution of the diffusion equation with variable diffusion coefficient (as a function of moisture content) and based on the result (isothermal drying rate versus average moisture content) of drying experiments with a slab at different temperatures. By using this method, the information obtained from the isothermal drying experiments was applied to other geometries and nonisothermal conditions. Dimensionless variables (moisture content, time, space coordinate, and diffusion coefficient) were used in the study. Instead of solving numerically the nonlinear diffusion equation, a calculation procedure was applied in a step-by-step manner where the conditions for each new step were obtained from the experimental drying curves.

In order to measure moisture profiles during drying, some advanced methods such as scanning neutron radiography or nuclear magnetic resonance imaging (MRI) were used in a number of studies. Among others are the works of Blackband and Mansfield [27] on solid blocks of nylon, Schrader and Litchfield [28] on food gel, Pel et al. [29] on brick and kaolin clay, McDonald et al. [30] on sandstone and rock, and Koptyug et al. [31] on alumina pellets. By means of these advanced methods, the measured moisture profiles can be used directly to determine the effective diffusion coefficient. In the work of Schrader and Litchfield [28], MRI was used to measure moisture profiles in a cylinder of food gel during drying at room temperature. The measured profiles were then compared with the numerical results calculated by the diffusion model, and in this way, the effective diffusion coefficient was computed. As an example, Figure 1 shows the variation of the obtained effective diffusion coefficient at t = 30, 45, and 60 minutes of drying. It was pointed out that the diffusion model seems not to be a good method to predict the interior moisture profile of food gel. Koptyug et al. [31] used MRI to study the diffusion of water in alumina pellets and showed that MRI can provide good information on the real-time variation of liquid content in the course of drying of porous solids. MRI is a good experimental method because it is able to measure moisture content at any point within a complex material. Additionally, it provides a quick, accurate, nondestructive method and therefore allows the evaluation of various drying models.

More recently, Guillard et al. [32] used the diffusion model to predict the moisture distribution in multicomponent heterogeneous food where components of high/low water activity are placed adjacent to one another. The calculation of moisture distributions compared well with the experimental results, and it was proposed that, in this special application, the model could be useful. Efremov [33] used another approach for the description of the drying kinetics of porous materials. The approach is based on the analytical solution of the diffusion equation (for one-dimensional isotropic diffusion) with a flux-type boundary condition in the form of mass flux. In this work, the drying kinetics (dimensionless moisture content versus time and drying rate) were determined by applying the Laplace transformation to express the mass flux. Porto and Lisbôa [34] developed a three-dimensional model based on the diffusion model with the constant diffusion coefficient in order to describe the drying process of a parallelepipedic oil shale particle. Lim et al. [35] introduced an equation derived from a diffusion equation in which the diffusion coefficient is a function of space. Akpinar and Dincer used the diffusion model to investigate moisture transfer in a slab of potato [36] and in eggplant slices [37]. The works dealt with drying processes at different air temperatures and flow velocities. The influence of boundary conditions on drying process was investigated. The model is limited to the one-dimensional problem of an infinite slab. In their works, the thermophysical properties of the drying material are taken as constant, and the effect of heat transfer on the moisture loss is neglected.

The diffusion model can lead to an incorrect prediction and misinterpretation of the moisture distribution or of the drying behavior due to the fact that only moisture transport is considered and that the physical meaning of the diffusion coefficient is either lost (in the case of a constant) or it becomes a lumped parameter of all simultaneous effects (in the case of a variable). However, this model is still used as a simple way to describe drying in some cases.

2.2. Receding Front Theory

Different versions of the so-called receding front model were developed in order to obtain a better understanding and describe the influence of other mechanisms (capillarity, gravity, or external forces in gradients of pressure and temperature) on the motion of water during drying. According to this model, at the critical point (when the falling rate period starts), an evaporation front arises and gradually moves into the interior of the body [11]. The moving evaporation front divides the system into two zones: the wet and the dry zones as shown in Figure 2. For a hygroscopic material, the dry zone is called the sorption zone due to the adsorptive nature of moisture retention. In the dry zone, the free water content is zero and the main mechanism of moisture transfer is vapor flow. However, in this region, the movement of adsorbed water may also play an important role [38]. During drying, the position of the receding evaporation front varies with time.

The receding front model was first developed in the 1960s. A review of the development of the receding front model can be found in the work of Tsotsas [11]. The simplest version of the receding front model is a model where saturation S is 1 in the wet region and 0 in the dry region. In what follows, the model presented by Chen and Schmidt [38] is shown as one example. According to Chen and Schmidt, the set of one-dimensional equations describing the coupled heat and mass transfer can be written as (subscripts 1 and 2 denoting the wet and dry zones)

Wet zone ():where is the liquid transfer coefficient and is the specific heat capacity of water. The term is the moisture content of free water and is the effective thermal conductivity, which is calculated bywhere is the thermal conductivity of liquid, is the evaporation enthalpy, is the saturation vapor pressure, and is the vapor transfer coefficient. The vapor transfer coefficient covers the contribution of both convective and diffusive flows:where is the ratio of air and vapor diffusion coefficient, is the relative permeability of the gas phase, is the dynamic viscosity, and is the vapor diffusion coefficient.

Dry or sorption zone ( ):where is the specific heat capacity of vapor, is the thermal conductivity of vapor, is the adsorbed water content, and is the adsorbed water transfer coefficient. For a nonhygroscopic material, is zero and is negligible. denotes the molar mass of vapor, and is the partial vapor pressure.

In addition to the above equations, the mass and heat transfer at the moving boundary must fulfill the following conditions:

Sorption isotherm is applied in the model, and the surface boundary conditions are needed. For more details, refer [38].

A drawback of the receding front approach is that the diffusion equation is used instead of more fundamental concepts like capillary pressure, liquid pressure gradients, and permeability to describe the capillary activity in the wet zone. In addition, heat transfer is only described by an effective thermal conductivity. Difficulties appear in determining the boundary of the moving evaporation front and the coefficients for heat and mass transfer, which are the functions of dry and wet zones.

2.3. Drying Model of Philip and De Vries

Philip and De Vries [39] and De Vries [4041] extended the previous treatment of the diffusion equations by including effects of capillary flow and vapor transport. In their work, the thermal energy equation was also incorporated into the set of the governing equations to describe the drying process. This set of equations was treated under the combination of moisture and temperature gradients. The obtained system consists of diffusion-like equations whose coefficients must be determined by experiment. The model is briefly presented below.

2.3.1. Liquid Water Transfer

Free liquid water movement is macroscopically described by Darcy’s law:where is the gravity potential. By expressing the term as a function of and and by substituting this function into (10), the liquid water flux can be written as a combination of three components due to the moisture gradient, temperature gradient, and gravity [42].where and are the isothermal and thermal diffusivities of water given by

2.3.2. Water Vapor Transfer

The transport of water vapor by molecular diffusion is described macroscopically by Fick’s first law and by using the assumption of a steady diffusion in a closed system between an evaporation source and a condensation sink, a commonly used expression for the vapor flux in terms of moisture and temperature gradients:where and are the isothermal and thermal diffusivities of vapor, respectively. These terms are expressed as [42]where is a function of porosity and moisture content, is the diffusion coefficient of vapor in air, is the gravitational acceleration, is the saturation vapor pressure, is the average air temperature gradient, and and are the densities of vapor and liquid.

In addition to liquid and vapor transfer, Philip and De Vries assumed that the gas pressure can be treated as constant and the gas phase momentum equation can be ignored.

2.3.3. Mass and Heat Conservation Equations

The partial differential equations of mass and energy are formulated as follows [9, 39]:where is the overall thermal mass diffusivity, is the overall isothermal mass diffusivity, denotes the thermal conductivity, and is the volumetric heat capacity of the moist porous medium. Note that convective energy terms are assumed negligible.

The theory of Philip and De Vries has become generally known and has been applied to porous media other than soil, which was chosen to investigate the heat and mass transfers by the authors. The major restriction of the theory is that it does not include the gradient of gas pressure, there is no convection contribution in heat equation, and the coefficients of the model are complicated.

2.4. Luikov’s Theory

Independent of Philip and De Vries’s work, Luikov [4345] investigated the heat and mass transfer during drying of capillary-porous bodies by employing the principles of irreversible thermodynamics. In this theory, the total moisture flux is assumed to be made up of three components: the first one is due to a gradient in moisture content, the second due to a gradient in temperature, and the last due to a gradient in the total pressure [46]:where is the total moisture flux, is the moisture diffusion coefficient, is the thermal gradient coefficient, and is the pressure gradient coefficient.

The conservation equations of Luikov’s model are written in the following form [46]:where , is the thermal conductivity of the moist body, is calculated from Darcy’s law (10), with as the filtration coefficient, and is the mass rate of evaporation per unit volume:where is a dimensionless factor characterizing resistance to vapor diffusion in a body.

By using the above conservation equations, three interdependent partial differential equations involving variables , , and can be obtained aswhere the kinetic coefficients depend not only on temperature and moisture content but also on material properties and drying conditions. For example,where is the coefficient of moisture conductivity, is the moisture capacity, is the density of the dry solid, and is the thermogradient coefficient related to the moisture content difference. For more details on the computation of these kinetic coefficients, refer Luikov [45], Irudayaraj and Wu [47], and Lewis and Ferguson [48].

It is noted that under the assumption of constant gas pressure, Luikov’s equations are similar to those proposed by Philip and De Vries. The biggest problem encountered in using Luikov’s equations is the definition of the coefficient . In practice, it is often not possible to obtain these parameters to solve the full system of equations. However, Luikov’s theory provides a well-established model in the treatment of simultaneous heat and mass transfer of the drying problem. The solution of Luikov’s partial differential equations was studied numerically by Irudayaraj and Wu [47] and Lewis and Ferguson [48]. These equations are still commonly employed today and quite often solved by the finite element method.

2.5. Krischer’s Theory

Krischer is also among the first researchers who had investigated the role of heat and mass transfer during drying of porous media. The research work of Krischer was and is still used today as a basis for much of the development in drying theory. In his work [49], which was first published in 1956, Krischer proposed a set of equations to describe the moisture transport for several geometries (plate, cylinder, and sphere). Krischer assumed that moisture transfer is controlled by the combined influence of capillary flow of liquid and diffusion of vapor.

In Krischer’s model, the liquid flux is calculated from Darcy’s law (10), and the vapor flux is written aswhere (>1) is called the diffusion resistance factor and describes the decrease of the vapor flow in the considered material in comparison with that in a stagnant gas.

By using the theory of Krischer, Berger and Pei [50] included the sorption isotherm (empirically obtained) into the model as a coupling equation among liquid, vapor, and heat transfer. Based on Krischer’s theory, Berger and Pei introduced two balanced equations for heat and mass (the gas pressure was taken as constant). In this model all phenomenological coefficients (e.g., liquid conductivity, vapor diffusivity, and thermal conductivity) are taken as constant. Heat transfer is assumed to take place only by conduction through the solid skeleton. The overall mass and heat balance equations proposed by Berger and Pei [50] are expressed in terms of moisture content and vapor pressure as follows:where is the vapor diffusivity; is the liquid conductivity; , and are the volume fraction of solid, water, and gas, respectively.

The main difficulties encountered in using Krischer’s model to predict the drying rate are the assumption of surface boundary conditions (Krischer postulated that “at the surface of the drying material, the corresponding equilibrium values of the dependent variables were reached instantaneously at the beginning of the drying process”) and the application of the sorption isotherm for the whole range of moisture content [50]. Even though sorption isotherm is taken into account, the approach of Berger and Pei does not offer much innovation over Luikov’s and Philip and De Vries’ models [46]. In addition, as for the previous models, experimental tests are needed to ensure its validity.

3. Whitaker’s Model and the Framework for the Numerical Simulation

In the late 1970s and early 1980s, Whitaker [51, 52] presented a set of equations to describe the simultaneous heat, mass, and momentum transfer in porous media. Based on the traditional conservation laws, the model proposed by Whitaker, an important milestone in the development of drying theory, incorporated all mechanisms for heat and mass transfer: liquid flow due to capillary forces, vapor and gas flow due to convection and diffusion, and internal evaporation of moisture and heat transfer by convection, diffusion, and conduction. By using the volume averaging method, the macroscopic differential equations were defined in terms of average field quantities.

Based on the model of Whitaker, a system of governing equations can be built to represent the drying process, in which the most important equations are the conservation equation for water in the liquid and the gas phases, the conservation equation for air in the gas phase, and the conservation equation of energy of the whole porous system under consideration. These equations can be solved numerically and form the framework for the numerical simulation of the drying process, in which simultaneous mass and heat transfer together with phase changes (vaporization) can be taken into account. We will briefly discuss these equations here.

The first equation governing the drying process is the conservation equation for water in both liquid and gas phases. This equation can be written aswhere , , and are the mass density of the liquid, vapor, and gas phases; and are the volume fractions of the liquid and gas phases; and are the velocities of the liquid and gas phases, and is the effective diffusivity tensor. The conservation equation for water states that the mass of water at each point inside the porous body is conserved.

The second equation governing the drying process is the conservation equation for air in the gas phase:which also states that the mass of air at each point inside the porous body is conserved.

The third equation is the conservation equation of energy, which can be formulated as

Note that in the above system of equations, the velocity of water can be computed using the equation of motion for the liquid phase:and the velocity of air can be computed using the equation of motion for the gas phase:where the dynamic viscosity of water and of gas are temperature dependent, is the absolute permeability tensor, and and denote the relative permeability tensors of liquid and gas, respectively.

In addition to the above governing equations, the conditions for mass and heat transfer at the external drying surfaces of the porous medium must be specified. For example, it can be assumed that, at the external drying surfaces, the fluxes of mass and heat are described for convective drying by the boundary layer theory with Stefan correction:where and are the fluxes of water and heat respectively, and are vapor pressure and temperature of bulk drying air, is the outward-pointing normal vector at the boundary surface, and and are mass and heat transfer coefficients. Additionally, we can assume that the gas pressure at the external drying surfaces is fixed at the pressure of the bulk drying air:

The advantage of Whitaker’s model is that it offers a very good representation of the physical phenomena occurring in porous media during drying. However, the problem encountered in using Whitaker’s model is the difficulty in determining its complicated transport coefficients, such as the effective diffusivity and permeabilities, which depend strongly on the material properties and structure. These parameters are either functions of moisture content or temperature or both moisture content and temperature. In addition, solving the coupled equations of heat and mass transfer, which are strongly nonlinear, requires very complicated numerical methods.

The theory of Whitaker was further developed and applied in drying analysis of various porous media, for example, in the drying analysis of sand by Whitaker and Chou [53], Hadley [54], Oliveira and Fernandes [55], and Puiggali et al. [56], analysis of glass beads by Quintard and Puiggali [57] and Kaviany and Mittal [58], analysis of sandstone by Wei et al. [59], analysis of porous insulators by Tien and Vafai [60], analysis of brick by Nasrallah and Perré [61], analysis of cellular materials by Crapiste et al. [62], analysis of wood by Spolek and Plumb [63], Michel et al. [64], Perré [65], and Lartigue et al. [66]. In these works, the model was usually quite successfully matched against experimental data.

Whitaker and Chou [53] simplified the theory to obtain two nonlinear equations for the distribution of saturation and temperature. In this work, the gas pressure was assumed as constant, the gas momentum equation was ignored, and a quasisteady state was applied. It is interesting to note that there is a resemblance of these two equations to the equations proposed by Luikov and by Philip and De Vries [46]. In this simplified case, the comparison between theory and experiment was made by Hougen et al. [15]. The important conclusion is that the gas phase momentum equation must be included in solving the comprehensive set of equations. Crapiste et al. [62] applied Whitaker’s model to investigate the drying of cellular materials. To validate the model, a comparison of one-dimensional drying to the experimental drying of apple and potato was presented and a good agreement was found. Wei et al. [59] applied Whitaker’s model to the drying of a cylinder of sandstone subjected to convective heating. The obtained partial differential equations in one dimension were solved by a three-point, two-level implicit finite difference method. The calculated results were compared with experimental results and showed a quite good agreement.

Ferguson [67] focused on a two-dimensional problem of the high-temperature drying of spruce. The numerical results highlighted the advantage of the discretization technique (control volume finite element method) in solving the problem with structured and unstructured meshes. A numerical investigation was conducted by Boukadida et al. [68] to study the convective drying of a slab of clay-brick. The work analyzed the influence of the properties of the surrounding drying agent (temperature, gas pressure, and vapor concentration) as well as the initial medium conditions (temperature and moisture content) on the drying process by considering several configurations. However, the full investigation of the effect of the boundary layer on the coupled heat and mass transfer still requires further work, as concluded by the authors. Silva [69], based on Whitaker’s theory, presented a general model to describe the momentum, heat, and mass transfer in drying problems with moving boundary. By using the volume averaging method, a set of equations for multiphase systems was applied to porous media. Numerical results showed a good agreement with the experimental data of kaolin drying.

One of the most significant advances in developing Whitaker’s theory as well as in modelling the drying of porous media comes from the works of Nasrallah and Perré [61] and Perré et al. [70]. In their work, the drying of two quite different porous media—clay-brick and softwood—was investigated. The most important advance in the work of Perré is the consideration of bound water [7173]. By considering bound water, the driving potential for bound water migration was assumed to be proportional to the gradient in the bound moisture content. For the case of wood, Perré and his colleagues introduced two equations to calculate the transport of this kind of bound waterwhere is the rate of bound water evaporation and the subscripts b and c denote bound water and cellulose matter, respectively. For wood, the diffusion coefficient of bound water is calculated in m2/s from the following equation [72]:where is moisture content of bound water and is temperature (Kelvin).

In his work, Perré solved the one-dimensional problem of drying with three state variables (temperature, pressure, and moisture content). The control volume method was applied to solve the nonlinear partial differential equations. The mathematical schemes for equidistant and nonequidistant meshes were discussed [61]. The authors also investigated the sensitivity upon model parameters by numerically varying the effective diffusivity, effective thermal conductivity, intrinsic and relative permeabilities, and external drying conditions (heat and mass transfer coefficients).

With the rapid development of the computer technology, modern computers allow the simulation of drying not only in one dimension but also in two and three dimensions. Besides, numerical methods are also more efficient in obtaining accurate results and reducing the computational time. Among the advancements during the 1990s in the study of drying of porous media is the simulation of drying processes in two dimensions with unstructured meshes proposed by Perré [74] and Perré and Turner [75]. The first comprehensive three-dimensional drying model using structured meshes was introduced by Perré and Turner [76]. In this work, a homogeneous model, which employed the full set of conservation equations, was considered. A cube of light concrete (isotropic medium) and a board of wood (anisotropic medium) were chosen to investigate the influence of the number of exchange faces. Several simulation results for low- and high-temperature drying of softwood were presented and discussed. By comparing the different simulation results, the study showed that a three-dimensional model is required to describe correctly the drying behavior of porous media.

Concerning the heterogeneity of material properties, Perré [74] developed a heterogeneous drying model for wood. The variation of the material property information such as capillary pressure and absolute permeability was taken into account with the help of experiments [77, 78]. The material information of wood obtained from this work was later applied to a two-dimensional heterogeneous drying model [79]. In this work, the effects of material heterogeneity and local material direction on the heat and mass transport during drying were investigated. Two cases of low- and high-temperature drying were considered. Following this direction, more recently, Truscott [80] and Truscott and Turner [81] developed a three-dimensional heterogeneous drying model for wood. The work considered the heterogeneity of the material properties, which vary within the transverse plane with respect to the position that defines the radial and tangential directions. Two nonlinear partial equations for moisture content and temperature (pressure was assumed as constant) were solved.

4. Numerical Simulation Based on Whitaker’s Model

By solving the system (23)–(30) resulting from the model of Whitaker, the drying process in porous media can be simulated numerically taking into account complex mass and heat transport phenomena. Different numerical methods can be used to solve the above system of equations, for example, the finite element method, the finite difference method, or the control volume method.

In simple cases, the numerical solution can be obtained with relatively small computational effort. One of such examples can be found in [52], where numerical simulation was compared with experimental measurement [14] of a sand plate under isothermal drying conditions. In the work, the above system (23)–(30) was applied. The plate of sand has a thickness of 5.08 cm. The drying took place by air. The initial saturation was set to be 100%. One surface of the plate was considered impermeable. The other surface was in contact with air. Since the width and length of the plate are much larger than its thickness, the problem can be reduced to be 1-dimensional. The numerical solution was obtained with the help of the finite difference method. Note that during the experiment [14], the averaged saturation was determined at different time instants. Corresponding to these time instants, the saturation profiles were measured. The measurement delivered then the saturation profiles at different values of averaged saturation, namely, at 36%, 53%, 79%, and 89%. The numerical solution was then compared with the experimental result as shown in Figure 3. In Figure 3, the saturation profile is plotted as a function of the normalized distance from the impermeable surface of the plate. The comparison shows that the model of Whitaker delivers a reasonable result. The adequacy of Whitaker’s model in modeling the drying process of porous media can also be seen by comparing this model with other models mentioned above. From Table 1, it can be observed that Whitaker’s model offers the most complete picture of the different phenomena happening during the drying of a porous medium.

In this work, the control volume method is chosen to solve the model of Whitaker. The reason for using the control volume method lies in the fact that this method satisfies the conservation requirement of the basic physical quantities such as mass and energy at any discrete level. This means without the need to enforce this requirement by using additional constraints, the heat and mass flows across a boundary of a control volume element or over the boundary of the whole porous medium are automatically conserved. In this work, we will not discuss the details of the application of the control volume method in solving the system (23)–(30) and limit ourselves in presenting a numerical example in order to demonstrate the capability of the approach using the model of Whitaker. For more details on the approach and numerical implementation, refer [12].

We examine here the drying of a sphere of light concrete with radius R = 2.5 mm. The temperature of the sample at the start of the drying process is T0 = 20°C, the moisture content before drying is X0 = 1, and the pressure of air inside the sample is P0 = 1 bar. In order to dry the sample, the sample is put in an oven with dried air. The temperature of the oven is T = 80°C and the vapor pressure in the oven is Pv, = 0. At the boundary of the sample, heat and mass transfer takes place during the drying process. For these two processes between the sample and the surrounding air of the oven, we assume that the heat transfer coefficient is α = 14.25 W·m−2·K−1 and mass transfer coefficient is β = 0.015 m·s−1. Since the sample is symmetric and the boundary conditions applied to the sample can also be considered as symmetric, the drying problem of our sample can be solved by the control volume method in one-dimensional space. For other details concerning transport properties, refer [12].

As numerical results, the temporal evolution of moisture, temperature, and pressure for approximately every 0.5 mm in distance along the radial direction of the sample from the center to the boundary of the sample is shown in Figures 46. In Figure 4, the dashed curve presents the average moisture content of the whole sample.

From the numerical results presented here, some important drying characteristics can be observed. Starting from a uniform initial moisture content X0 = 1 (corresponding to saturation S0 = 62.5%) over the whole sample, after a short warm-up period, the moisture content is reduced at a constant rate (constant slope of the curves in Figure 4). This is called the first drying period (or constant rate drying period). During the first drying period (approximately 17.6 minutes), free water is brought to the surface by capillary forces where heat supplied by the convectional air is used for the rapid vaporization of water. Due to this reason, the sample remains at the wet bulb temperature of the drying air (Twb = 23.81°C, Figure 5). The moisture gradient increases (relative permeability kw decreases), and our analysis shows that the moisture profiles as the function of radius appear fairly flat. Within this period, the pressure stays constant at atmospheric pressure (Figure 6). As the drying process continues, the moisture content at the surface reaches the irreducible value Xirr = 0.07 (the average moisture content of the whole sample reaches the critical moisture content Xcr = 0.1774) and the second drying period commences. In the second drying period, the dominating forces are viscous forces. Our analysis shows that a front separating the regions of adsorbed water and free water recedes from the surface into the sample. This process is finished when the moisture content everywhere in the sample is below the irreducible value, that is, when all free water of the sample has been removed. During this period, heat transfer is almost unchanged (resistance in the sample is slightly increased), but mass transfer experiences an important additional resistance. Heat is used not only to evaporate water but also to increase the temperature of the sample. Therefore, the temperature of the sample starts to rise from the wet bulb temperature. The center of the sample stays cooler than the outside (Figure 5). This is due to the fact that the evaporation of water takes place not at the surface but at a place inside the sample (front). The free water region heats up until a new equilibrium is attained (if we assume a stationary front). At the front, heat is consumed for evaporation. As we can see from Figure 6, in the second drying period, an overpressure appears due to the Stefan effect. The pressure inside the sample increases to its maximum value, whereas the pressure at the surface always stays at the atmospheric pressure (1 bar). When the receding front has passed through the whole sample, the sample becomes dry and the entire porous medium is in the hygroscopic zone. The moisture content goes down to the equilibrium value. The temperature of solid gradually approaches the dry bulb temperature of the drying air (T = 80°C, Figure 5), and the pressure falls back to the atmospheric pressure (Figure 6).

By considering the numerical example presented here, it is easy to see that the implemented simulation framework based on the model of Whitaker can be used effectively in studying the complex drying process of porous media. Especially important is the incorporation of the simultaneous heat and mass transfer processes together with the evaporation process in the simulation.

5. Conclusion

In this work, a review of the development of some drying models, their application, and their restrictions is presented. Among the most complex and modern models is the model developed by Whitaker. By using Whitaker’s model, the basic laws of mass and heat transport can be formulated at macroscopic level. The result is a system of equations governing the drying process of porous media, which can be solved numerically using modern numerical methods, in particular the control volume method. Our numerical implementation of Whitaker’s model is presented through a numerical example to show to capability of this model. The numerical example and the different characteristics of drying observed in our simulation results show that the implemented simulation framework can be used to study realistic drying processes in porous media.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.