A Framework and Numerical Solution of the Drying Process in Porous Media by Using a Continuous Model

e modelling and numerical simulation of the drying process in porous media are discussed in this work with the objective of presenting the drying problem as the system of governing equations, which is ready to be solved by many of the now widely available control-volume-based numerical tools. By reviewing the connection between the transport equations at the pore level and their up-scaled ones at the continuum level and then by transforming these equations into a format that can be solved by the control volume method, we would like to present an easy-to-use framework for studying the drying process in porous media. In order to take into account the microstructure of porous media in the format of pore-size distribution, the concept of bundle of capillaries is used to derive the needed transport parameters. Some numerical examples are presented to demonstrate the use of the presented formulas.


Introduction
e drying process plays an important role in many different industries, for example, in chemicals, pharmaceuticals, and agriculture.e drying process is one of the most complex problems that one finds in process engineering because not only heat and mass transfer takes place simultaneously in the course of the process but also because other phenomena may play a significant role.Although drying processes have been studied experimentally and theoretically for decades, simulating the coupling of heat and mass transfer and other phenomena in drying is still a challenging problem.Many researches were carried out to build suitable models and simulate numerically the drying process in engineering, and among recent works are those of Sekki and Karvinen [1], Antonov et al. [2], Azmir et al. [3], Ramos et al. [4], and Wu et al. [5].
Besides theoretical developments, numerical methods were applied successfully to simulate the drying process of porous media at the macroscopic scale and at the microscopic scale [6].At the microscopic scale, the drying of porous media can be modelled as a network of pores, and the motion of the liquid-gas interface is modelled at the pore level, for example, in the work of Laurindo and Prat [7], Prat [8], Segura and Toledo [9], Metzger et al. [10], and Hirschmann and Tsotsas [11].By using this approach [12], which we will refer to as the discrete approach, the microscopic structure and therefore the transport properties of the porous medium can be modelled with better accuracy.However, the problem becomes very large, and solving the system of equations of coupled heat-mass transfer becomes in many cases impractical, in particular when dealing with systems of large geometrical dimension.In such cases, the use of the continuous approach is more relevant (see for example [13,14] or [15]).e continuous approach is built on the assumption that the porous medium can be described as a fictitious continuum by using effective coefficients for heat and mass transfer.For many applications, by using the continuous approach, the drying characteristics of porous media can be simulated with a very good accuracy.However, one of the challenges in using the continuous approach is how to determine the transport parameters [16].Note that since the continuous approach is built on a fictitious continuum, this is also called the continuum approach.A continuous model for the drying process is therefore also called a continuum model.
In developing a drying model at the macroscopic scale, Whitaker [17] used the volume averaging technique to derive a system of macroscopic transport equations from a set of basic transport laws at the microscopic level.In Whitaker's work [18], a porous medium is assumed to be equivalent to a continuum.A set of conservation equations for mass, energy, and momentum are introduced using average state variables.e continuous model developed by Whitaker is considered as rigorous and the most advanced continuous model today.e theory of Whitaker was later applied to different porous media, for example, by Perré [19], Ouelhazi et al. [20], Boukadida and Nasrallah [21], Boukadida et al. [22], Ferguson [23], Perré and Turner [24], Truscott [25], and Truscott and Turner [26].Numerical techniques were developed to simulate the drying process using the derived average conservation equations.Among others, Perré and Turner [24] employed the control volume method to solve the problem.e advantage of this numerical method is that it ensures the conservation of mass and enthalpy through the boundaries of elements.
Despite the fact that the derivation of the governing equations of the drying problem at the continuum level can be found elsewhere [18], some effort was made to put these equations into a format that can be solved numerically [24], and there is the need to put all available knowledge in one framework, which is easy to understand and ready to be solved by many of the now-available numerical tools (as example of such tools, see [27,28] or [29]).Such a framework will not only offer us the tool to solve the drying problem quickly but will also allow us to modify the governing equations in order to reflect the different phenomena that are not yet taken into consideration.In this work, by following the previous foundation laid out by Whitaker [18], Perré and Turner [24], and others, we will revisit the continuous approach starting with the transport equations at the pore scale (microlevel).We will briefly review the set of transport equations at the macro level.After that, we will introduce the transport parameters as a function of the material microstructure, and finally, some numerical solution will be presented and discussed.

Governing Equations at the Microscopic Scale (Pore Scale)
We consider here a rigid porous medium V with external boundary zV in which the matrix is made of some solid material and a system of interconnected voids.e voids are also called here "pores".ese pores are connected as a network of pores (voids).At the microscopic level (or pore level), we consider a small part dV of the porous medium with three phases: solid, liquid, and gas (Figure 1).e solid phase is denoted by s, the liquid phase (water) is denoted by w, and the gas phase is denoted by g. e gas phase has two components (species): air (denoted by a) and vapour (denoted by v).In drying analysis, one of the primary objectives is to compute the distribution of moisture content, temperature, and internal gaseous pressure within the porous medium during the drying process.At the pore level, the (local) moisture content, the temperature, and the gaseous pressure at each point can be determined using suitable laws of physics such as the conservation of mass, the conservation of linear momentum, and the conservation of energy of each phase: solid, liquid, and gas.ese conservation laws will be presented in the following, according to the work of Whitaker [18].
Note that we do not consider the particular shape of each pore (void) in the pore network.As presented later, the property of the pore network will be represented by the total amount of void volume in terms of porosity.e equations presented in this section are valid for each phase (liquid, air, and vapour) inside dV.By transforming the transport equations presented in this section to the macroscopic level (continuum level, Section 3), the influence of the pore shape, the number of pores, the size distribution of the pores, and how they are connected are represented by some effective  International Journal of Chemical Engineering transport properties, which can be determined by experimental measurement.
2.1.Mass Conservation at the Pore Scale.By assuming that no chemical reaction happens during drying, the total mass conservation equation of each phase can be written as and the mass conservation equation of each component in each phase is In the above equations, the first terms on the left-hand side are the change of mass due to accumulation, the second terms express the change of mass due to convection, v is the mass average velocity, and ρ is the total mass density of the phase under consideration.
where N denotes the total number of components, ρ i is the density, and v i is the velocity of component i.

Conservation of Linear Momentum at the Pore Scale.
By neglecting the contribution of body forces such as the gravitation force, the conservation of linear momentum for each phase can be written as where T is the stress tensor.e conservation of angular momentum requires this tensor to be symmetric: (5)

Energy Conservation at the Pore Scale.
e conservation of energy for each phase can be presented in the following format: where h is the enthalpy per unit mass, q is the conductive heat flux vector, τ is the viscous stress tensor, the term τ : ∇v is the viscous dissipation, P is the pressure, DP/Dt is the compression work, and Φ is the source or sink of electromagnetic energy.e conductive heat flux vector q is computed by Fourier's law: where λ represents the thermal conductivity.If the contribution of Φ is neglected and if in addition, we assume that, for liquid and gas phases, the viscous dissipation and the compression work are negligible: and then the energy equation is reduced to It is also assumed that the enthalpy is independent of pressure and that all heat capacities are constant so that the following relationship holds, where c p is the specific heat capacity and T R is the reference temperature (chosen here to be T R � 273.15 K).After having the above conservation equations, we can apply them to each phase of the three-phase system under consideration.Note that we are still considering these equations at the pore level.

Solid Phase.
e solid phase is considered to be rigid and fixed in space with zero velocity (movement of the solid phase is not considered): is assumption means that, for the solid phase, we need to study only the conservation of energy equation (6), which now becomes By making use of equation (7) and taking into account the assumption (10), we have the following conservation of energy for the solid phase: 2.5.Liquid Phase.For the liquid phase, which contains water as the only component, the mass conservation equation is e use of equations (7) and (10) allows us to write the energy conservation of the liquid phase in the following format: 2.6.Gas Phase.e gas phase is more complicated than the solid and the liquid phase since it contains two components: air and vapour.e mass conservation of the gas phase can be written in the following form: International Journal of Chemical Engineering By writing the component velocity v i in terms of the mass average velocity v g and the diffusion velocity u i , We can write the mass conservation of air and vapour in the following form: Furthermore, by expressing the diffusive flux ρ i u i as where δ v,a is the binary molecular diffusion coefficient for vapour and air, we finally have Note that for a multicomponent phase, the appropriate form of the energy equation ( 9) can be written in the following format: where h i is the enthalpy per unit mass of component i, and the mass average enthalpy h is defined in a similar way as the mass average of velocity: By using equation (21), for the gas phase, we have where In addition to the above conservation equations, ideal gas laws are assumed for partial and total gas pressures: where i stands for a, v, or g,  R is the ideal gas constant, and  M i stands for molar mass of air, vapour, or gas.e constraint for the partial and total gas pressure is simply Besides the set of equations listed above, the boundary conditions connecting the transport equations for the three separate phases need to be specified.Details can be found, for example, in the work of Whitaker [18].

Numeric Ready Continuum Equations of the
Drying Process e macroscopic transport equations can be obtained from the microscopic ones presented in Section 2 by using the volume averaging method [17].
is set of macroscopic equations is summarized as follows [14].

Conservation Equation for Water.
e mass conservation equation for water in the both liquid and gas phase can be written as follows: where ε w � V w /V and ε g � V g /V are the volume fractions of the liquid and gas phase with respect to the total volume V of the porous medium and D eff is the effective diffusivity tensor.

Conservation Equation for Air.
e second equation is the mass conservation equation for air in the gas phase:

Conservation Equation of Energy.
e third equation is the energy conservation equation for the whole system: where ε s � V s /V is the volume fraction and of the solid phase and λ eff is the effective thermal conductivity tensor.Note that, since the sum of the solid, liquid water, and gas is the total volume of the porous medium V s + V w + V g � V, the sum of all volume fractions is unity: In the above equations, the velocities are computed by using 4 International Journal of Chemical Engineering where K w and K g are the absolute permeabilities of liquid and gas phases, k w and k g are the corresponding relative permeabilities, g is the gravitational acceleration, and χ is the height.In this work, for simplicity, we will not consider the gravitational effect, and therefore, the terms with g∇χ drop out.In order to compute the gas pressure, we make use of equation ( 26): where the vapour pressure will be computed by the sorption isotherm, which can be determined by an experiment.e liquid water pressure is computed by the following equation: where P c is the capillary pressure.Also in the above equations, the enthalpies of the solid phase (h s ), liquid water (h w ), air (h a ), and vapour (h v ) are computed using (10) as where Δh vw is the latent heat of vaporization.For the computation of these values, we use here e value of c s can be given directly or indirectly depending on the ease of use and on the particular porous medium under consideration.For example, in the case of light concrete, instead of using directly c s , we can use an averaged heat capacity ρC p defined as follows: In solving the above system of conservation equations, we can use as main variables the temperature T, the volume fraction of the liquid water ε w , and the gas density ρ g at each point in our porous medium.However, since it is convenient to compute different quantities using moisture content X, we will use here as main variables the temperature T, the moisture content X, and the gas density ρ g at each point in our porous medium.e formula for moisture content X is written as By using moisture content, the volume fraction of liquid water ε w can be computed by equation (36) as ε w � ((ε s ρ s )/ρ w )X. e volume fraction of gas is where the volume fraction of the solid phase ε s can be computed by the porosity of the medium under consideration.e vapour pressure will be computed by the sorption isotherm and can be determined experimentally as the function of moisture content and temperature RT. e gas density is ρ a � ρ g − ρ v .Note that for simplicity, the density of water is considered here as constant ρ w � 1000 kg/m 3 .
Besides the above governing equations, the boundary conditions for mass and heat transfer at the external drying surfaces zV of the porous medium must be specified.Here, we assume 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.For the mass flux J w •  n at zV, we have and for the energy flux J e •  n at zV, where J w and J e are the fluxes of water and heat, respectively; P v,∞ and T ∞ are the vapour pressure and the temperature of bulk drying air;  n is the outward-pointing normal vector at the boundary surface; and β and α are mass and heat transfer coefficients.Additionally, the gas pressure at the external drying surfaces is fixed at the pressure of the bulk drying air: Sorption isotherm, capillary pressure, ideal gas laws, and enthalpy-temperature relations will help express all variables as functions of three state variables.
As summary, in solving the drying problem using the continuous approach (continuum approach), we will solve the 3 conservation equations ( 27)-( 29) with the corresponding boundary and initial conditions.e main variables to be solved are the temperature T, the moisture content X, and the gas density ρ g .e boundary condition for heat and mass transfer will be temperature and relative humidity of drying air.ese conditions will be used with the help of equations ( 38)- (40) in terms of drying air temperature (T ∞ ), vapour pressure in drying air (P v,∞ ), and total pressure of drying air (P ∞ ).e initial condition will be the temperature T, moisture content X, and gas pressure (P g ) at each point inside the porous medium to be dried at the beginning of the drying process.By using these values, the initial values of our main variables T, X, and ρ g can be computed.

Transport Parameters as Functions of Material Microstructure
A major problem in using the continuous approach (with the governing equations presented in the previous section) is the determination of the different model parameters for a given porous material.In this section, based on the work of Metzger and Tsotsas for a bundle of capillaries [30], we present a simple and effective way to compute the needed transport properties by taking into account the characteristics of the pore structure of porous materials, namely, the pore size and its distribution.In this approach, the capillary pressure and the absolute and relative permeabilities are computed as functions of material pore size and pore-size distribution.e effective diffusivity and effective thermal conductivity are considered as dependent on porosity ψ and saturation S. ese parameters will then be used in the continuous drying model presented in the last section, which is capable of modelling the spatial and temporal evolution of moisture content, temperature, and gaseous pressure.e aim here is to provide a tool that we can use to understand, on a fundamental basis, how a variation of pore-size distribution changes the drying behaviour of porous media.However, it should be pointed out that the pore-size distribution alone is not enough to characterize the drying behaviour of a porous medium; it is only the most accessible structural information.e idea to use a bundle of capillaries to describe the pore space is not new.As an example, Krischer and Kast [31] described liquid water transport in porous media using this geometry.However, the idea was rarely used for a systematic investigation of the whole drying process.
In using the concept presented by Metzger and Tsotsas [30], we assume that a bundle of capillaries represents the void space of a porous medium, as shown in the model depicted in Figure 2. In the model, the capillary tubes are set perpendicular to the surface of the porous body and the solid phase is arranged in parallel.
ere is no lateral resistance to heat or mass transfer between solid and capillaries, making the model strictly one-dimensional.We restrict ourselves to large enough pore sizes so that for every capillary the gas-liquid phase boundary can be described by a meniscus having a capillary pressure.During drying, larger capillaries will empty first, because they have lower capillary pressure.However, this capillary pumping is subject to friction leading to a nontrivial moisture profile.
For our investigation, two types of pore-size distributions (capillary radius distributions) are used.e first one is a monomodal normally distributed pore volume where r 0 is the mean pore radius, σ 0 is the standard deviation, and C is a constant.is normal distribution is truncated at r 0 ± 2.5σ 0 .e integral of the pore-size distribution computed in this truncated range must be equal to the void volume of the sample, or in other words, the total volume fraction of liquid for S � 1 must correspond to the porosity ψ of the sample.However, in this work, this integral is set to unity for the sake of simplicity.
e second distribution type is a bimodal distribution, which consists of two monomodal models (with different mean pore radii and deviations) named here as small pore and large pore distributions.In bimodal distributions, the volume fractions of small pores and large pores and the transition region between the two kinds of pores are taken into account.Similar to the case of monomodal distribution, the integral of the pore-size distribution in this case is set to unity.Evidently, different choices are possible for the capillary radius distribution.However, the mono-and bimodal distributions are considered to be enough for a systematic investigation.
If the porous medium is partially saturated with water, the assumption of ideal lateral transfer between the capillaries implies that, for a given local free water saturation S fw , small capillaries are filled up to a maximum radius r fill such that where r min is the smallest capillary radius of the bundle.e localization of free water S fw is the key for computing effective parameters as the function of saturation.Based on equation (42), the maximum radius r fill can be computed for a given S fw .For any saturation under the irreducible value, the saturation of free water S fw is zero and the maximum radius filled by liquid r fill is set to r min .
Adsorbed water, which may play an important role in the drying of hygroscopic materials, needs to be modelled separately since it depends mainly on material properties (the influence of pore-size distribution is neglected because the condensation due to the Kelvin effect only plays an important role at very high level of relative humidity).In our work, we chose the same type of temperature-independent sorption isotherm as Perré and Turner [24] used for concrete: where φ is the relative humidity, S irr is the maximum amount of adsorbed water, S � S irr + S fw , and in the presence of free water, φ � 1.
Besides vapour pressure, the capillary pressure is also linked to saturation by a state equation given by the meniscus in the largest filled capillary: where the zero contact angle is assumed and σ is the temperature-dependent surface tension.
Let us now consider one capillary, which is fully saturated by water.On the one hand, the volumetric flow rate is calculated from Poiseuille's equation: 6 International Journal of Chemical Engineering where L is the capillary length, r is the capillary radius, η is the dynamic viscosity (temperature dependent), and P w is the water pressure.On the other hand, the mean velocity v (volumetric flow rate per total cross section of the porous medium) of the liquid can be described by the generalized Darcy law.In this calculation, we assume that gravitational effects are negligible and that velocity is small enough to neglect inertial effects.If we apply Darcy law to a fully saturated capillary (k w � 1), we obtain By comparing equations ( 45) and ( 46), the absolute permeability can be found to be By extending the above formula to the bundle of capillaries, we get where the interval [r min , r max ] is the total range of the poresize distribution.e relative permeabilities of liquid and gas phases are computed in a similar way as the absolute permeability: Note that the sum of these two quantities is unity.
In order to illustrate the influence of pore-size distribution on the effective parameters, four different cases of pore-size distribution are considered [33].e parameters of the distributions and the corresponding permeabilities are presented in Table 1 and Figure 3.For the bimodal distributions, the volume is equally distributed to the two modes.
e capillary pressure curves for these four pore-size distributions are illustrated in Figure 4. Naturally, with the decreasing saturation, the capillary pressure increases.e overall level of the capillary pressure is determined by the mean pore size r 0 of the mode, and its range of variation is determined by the standard deviation σ 0 .
e effective vapour diffusivity does not depend on the distribution of the pores but on the evaporation area.
erefore, this transport parameter is assumed to be a function of saturation, porosity, and the binary diffusion coefficient δ va [14] where ψ is the porosity, which is the ratio of the total void or pore volume (V void ) to the total volume (V total ) of the material ψ � V void /V total .e next parameter is the effective thermal conductivity.Like the effective vapour diffusivity, this transport parameter is also assumed to be independent of pore-size distribution.As heat conduction occurs in all phases in parallel, the heat flux or thermal conductivity contributions must be weighted according to the respective volume fractions of the phases.If the contribution of gas is neglected, then the effective thermal conductivity can be computed as follows [14]: (52)

Numerical Solution Using the Control Volume Method
In principle, the finite element method, the finite difference method, or the control volume method can be employed to solve the governing equations of the drying problem presented in the previous sections.Many works were carried out trying to find the best technique for the simulation of the drying process.In the quest for a quicker, more accurate, and less expensive solution, even mixtures of different methods appeared, for example, the so-called control volume finite element method [23].In many textbooks on numerical methods for heat and mass transfer, the control volume method is praised for its accuracy in solving problems involving conservation of heat and mass [34].e method was Figure 2: Partially saturated bundle of capillaries: relationship between free water saturation S fw and maximum radius filled r fill [32].
International Journal of Chemical Engineering applied in drying simulations by, for example, Perré and Turner [24], Truscott [25], Truscott and Turner [26], Hadley [35], Nasrallah and Perre [36], Perre and Degiovanni [37], and Turner and Ferguson [38,39].e basic idea of the control volume method is simple.In this method, the calculation domain is divided into a number of nonoverlapping control volumes, each of which is associated with a grid point or node [34].e system of differential equations is then integrated over each control volume.Piecewise profiles expressing the variation of variables and related quantities are used to evaluate the required integrals.For each control volume, the result is a discrete version of the differential equations involving the variables related to the central node of this control volume and to the nodes connected to it.e most important feature of the control volume method that makes it different from other methods is that the requirement of conservation of the basic physical quantities such as mass and energy will be satisfied at any discrete level: across a control volume element, over a group of control volume elements, or over the whole calculation domain.In this section, we will present and discuss the results of some example problems.
In the first two examples, the drying process is considered with some given transport properties.In the third example, the transport properties are computed by using the microscale properties of the porous material, as presented in the last section.e purpose here is to showcase the validity of the above system of governing equations and the effect that the microstructure (in this case, pore size and pore-size distribution) has on the drying behaviour of porous media.In the first and in the second examples, a spherical particle is considered.In the third and fourth examples, we consider an infinite plate with a thickness of 200 mm.In practical application, the infinite plate can be used for the case the thickness of a real plate is small in comparison with its length and width.e first example is presented to demonstrate that the current approach can capture important effects which lead to large differences between isothermal and nonisothermal drying.e second example shows that simplification in models like the diffusion model and receding front model can lead to large differences in drying kinetics.is example also shows that the current approach can capture the different drying characteristics of the two simplified models.e third and the fourth examples demonstrate the capability of our approach in capturing changes of transport parameters (and correspondingly changes in the drying kinetics) due to changes in the microstructure of a porous medium (pore size and pore-size distribution).e fourth example also demonstrates the accuracy of the current approach when 8 International Journal of Chemical Engineering we compare the simulation results with those of the discrete approach.
In the first example, we examine here the drying of a sphere of light concrete with radius R � 2.5 mm. e initial temperature of the sample is T 0 � 20 °C, the initial moisture content is X 0 � 1, and the initial pressure is P 0 � 1 bar for the whole sample.
e material properties of light concrete used in this example are given as follows [24]: porosity ψ � 0.8, solid density ρ s � 2500 kg•m −3 , and heat capacity: e fully saturated material has a moisture content of X sat � 1.6.e sorption isotherm is where X irr � 0.07 is the irreducible moisture content.e saturation vapour pressure P * v (T) is given as a function of temperature T ( °C) by Antoine's equation: e capillary pressure is computed from where the surface tension σ is a function of temperature T ( °C): and X fw is the moisture content of free water: X fw � X − X irr .e absolute permeability is taken as constant: K � 2 × 10 −13 m 2 .e relative permeabilities for liquid k w and for gas k g phases are calculated from the following relationship: where S fw � (X − X irr )/(X sat − X irr ) and X sat is the saturated moisture content.e effective diffusivity is calculated from where k g is the relative permeability of gas and δ va is the binary diffusivity of vapour in air [31]: where T R � 273.15K and P R � 101325 Pa are reference temperature and pressure, respectively.e effective thermal conductivity has contributions from both the solid and the liquid phase (the contribution of the gas phase is neglected).It is computed as Note that in order to solve a drying problem, there are 2 approaches that we can use, depending on the characteristics of the accompanying heat transfer process.In the first approach, only mass transfer is considered.Here the temperature is assumed to be constant.We will call this isothermal approach.A typical example of this approach in solving the drying problem is the use of the diffusion model.By using the diffusion model, only mass transfer is considered.Despite the fact that the diffusion coefficient can be formulated as a function of temperature, the heat transfer is completely neglected, and the temperature of the porous medium is considered as unchanged during the drying process.e solution using this approach is accurate enough and becomes acceptable when the isothermal condition can be actually satisfied; for example, when the transfer of heat takes place quick enough, the temperature of the whole system remains constant.In the second approach, the heat and mass transfer processes are considered as coupled.Despite the fact that this system is often difficult to solve, the corresponding solution is closer to reality.e second approach must be used in the case isothermal condition is violated; for example, when the transfer of heat takes place too slowly, the constant temperature of the system under consideration cannot be guaranteed.We will call this the nonisothermal approach.Evidently, if the solutions obtained by the 2 approaches are very similar, the preferred approach is the isothermal one since this is the simpler approach.On the contrary, if the 2 solutions are very different, the nonisothermal approach must be used.
As we will see in the following, the drying problem considered in our first example does not satisfy the isothermal condition.
e temperature of the sphere is not constant, and the drying process is therefore nonisothermal.We will call this the nonisothermal drying.By solving this nonisothermal process, we would like to demonstrate that the approach presented in Section 3 above is capable of dealing with fully coupled mass and heat transfer processes.In addition, we would like to show that the presented approach can also handle the case of isothermal drying.In order to do this, we consider another case with the same drying conditions as stated above, but with different heat transfer properties, α � 6000 W/m 2 /K (heat transfer coefficient) and λ s � 6000 W/m/K (thermal conductivity of sample's solid phase, see Section 2).In this case, the effective thermal conductivity is dominated by thermal conductivity of the sample's solid phase.Under these conditions, the heat transfer process will happen much faster.As we will show later, under these conditions, the temperature of the sphere International Journal of Chemical Engineering is almost constant, and the drying process is thus isothermal.We will call this the isothermal drying.Finally, by comparing the result of the 2 cases (isothermal and nonisothermal drying), we would like to point out that not only the temperatures of the system are different but also the drying rates are significantly different.
We would like to emphasize that, in our simulation, there is no modelling difference between calculating the sample moisture content and temperature between the case of isothermal drying and the case nonisothermal drying.Here, the same system of equations of fully coupled heat and mass transfer is solved (see Section 3).e only the difference is the condition for heat transfer inside the sample and at the interface between the sample and the drying air.With relatively small effective thermal conductivity and small heat transfer coefficient, the result of our simulation shows that the temperature inside the sample is not constant and is not equal to the temperature of the drying air.With large effective thermal conductivity and large heat transfer coefficient, the result of our simulation shows that the temperature inside the sample is constant and is equal to the temperature of the drying air.
Since the boundary conditions applied to the sample are symmetric, the drying problem of the sphere can be solved by the control volume method in one dimension.e numerical results are presented in Figures 5-7.
e result presented in Figure 5 shows that, in the first drying period of the nonisothermal case, a cooling to temperature of 13.17 °C takes place, and the temperature rises back to the initial value T 0 � 20 °C in the second drying period.For the isothermal case (Figure 6), the temperature of the sample stays almost constant during the whole drying process: a cooling to temperature of 19.95 °C takes place in the first drying period, and the temperature rises back to the initial value T 0 � 20 °C in the second drying period.e 2 constant temperatures in the first drying period are very similar to the definition of the wet bulb temperature (if the supply of moisture were infinite, the temperature would stay at these 2 temperatures for ever).Due to this reason, we call these 2 temperatures imaginary wet bulb temperatures.In the isothermal case, the imaginary wet bulb temperature is then 19.95 °C, while this quantity is 13.17 °C in the nonisothermal case.
From the results presented in Figure 7, the evaporation rate in the first period and the critical moisture content for the isothermal case are _ m v,I � 0.1310 g•m −2 •s −1 and X cr � 0.1538.ese values are 0.0394 g•m −2 •s −1 and 0.1320 for the nonisothermal case.In the isothermal case, the initial drying rate is higher compared to the nonisothermal case.
is is due to the fact that the temperature in the first drying period of the isothermal case is significantly higher (19.95 °C) than that of the nonisothermal case (13.17 °C).Due to this reason, the effective diffusivity in the isothermal case is larger than that in the nonisothermal case (because this is a function of temperature, see equations ( 59) and ( 60)).Consequently, the total drying process is significantly longer in the nonisothermal case (290.9 minutes compared to 99.2 minutes).e results in this example show a clear difference between the two cases of isothermal and nonisothermal drying.
It is evident that, in many cases, mass transfer must be considered together with heat transfer to obtain realistic results.In these cases, the framework presented above is a suitable choice.
In the second example, we consider the isothermal drying of the same spherical particle of light concrete above.We assume that heat transfer is quick enough so that the temperature of the sample does not change during the drying process.e objective here is to demonstrate that even in the case of isothermal drying, simplification in modelling of mass transfer alone may lead to significant differences in the analysis of a drying process.We will consider here three cases.In the first case, the problem will be solved by using the continuous model (using the framework presented in Section 3 above, which is also called continuum model) and by using the control volume method as in the first example.However, in this case, the conservation of energy ( 29) is not considered, so that the isothermal condition is satisfied.In the second and third cases, the problem will be solved using the diffusion model, and the   International Journal of Chemical Engineering receding front model presented hereafter.e isothermal drying is set at T ∞ � 20 °C, the drying air has zero moisture content φ � 0, the pressure of drying air is set to 1 bar, the mass transfer coefficient is β � 0.015 m•s −1 , and the initial moisture content is X 0 � 1 as in the first example.
In using the diffusion model [31], the drying problem for a sphere can be expressed in the following format: where r is the radial coordinate and the diffusion coefficient δ * takes the constant value of 2.6 • 10 −5 m 2 •s −1 [31].
In addition to the diffusion equation above, we need to specify the boundary conditions for a sphere at radius r � 0 and r � R. At the centre of the sphere, due to symmetry, we have At the boundary of sphere, the mass flux must satisfy the boundary condition (38).e mass flux at the boundary r � R is computed as follows [31]: where ρ s is the density of the solid and _ m v is the evaporation rate.Since this mass flux must satisfy the boundary condition given by the boundary condition (38) in which the vapour pressure P v is given by the sorption isotherm defined above.As summary, the diffusion model is presented by the conservation of mass (62) and the 2 boundary conditions (63) and (65).In order to get the numerical solution, the system of 3 equations ( 62), (63), and (65) is solved by using the PDE solver pdepe in MATLAB.
In using the receding front model [40], in the case of an infinite plate, the porous medium is divided into 2 zones: wet zone and dry zone (Figure 8).We assume that there is only vapour in the dry zone, and in the wet zone, only liquid exists.e mount of adsorbed water is assumed here to be negligible.Note that the diffusion takes place in the dry zone.As mentioned above, we assume that heat transfer is quick so that the temperature of the sample is constant.Due to this reason, the flux of heat is not considered.For a plate, the relationship between the dry-wet front position and the momentary drying rate is where β is the mass transfer coefficient at the surface and z � s is the position of the front at time t. e diffusion coefficient δ * is computed as an effective diffusivity, as in equation ( 59), with saturation S � 0. As we can see in equation ( 66), the mass transfer resistance is obtained by addition: 1/β is the resistance in the boundary layer and s/δ * is the resistance in the dry zone.For a spherical geometry, the diffusion process takes place through a shell with inner radius R -s and outer radius R. is means the resistance in the dry zone can be computed as R • s/[δ * (R − s)], and the time needed to evaporate the amount of water contained in a shell is as follows: where A is the surface area of the sphere, ΔV is the volume of the shell, and ψ is the porosity.
In using the receding front model, for every given value of s we can compute the drying rate _ m v by using equation (66).By using this drying rate _ m v and by computing A and ΔV as function of s, the corresponding time t can be computed by equation (67).By using the computed drying rate _ m v , the total moisture inside the sample can be computed.In addition, by using the total moisture of the sample, the averaged moisture content X av is obtained.As a consequence, by varying s from 0 to R, we have a set of three values presenting the drying kinetics of the drying process at each value of s: (X av , t, _ m v ).e different solutions of the drying problem by using the three models are presented and compared in Figures 9-11.As we can see from Figure 9, all three models capture the initial value of evaporation rate _ m v,I � 0.262 g•m −2 •s −1 since for all three models, the mass transfer is initially only controlled by the boundary layer.From Figure 10, it can be seen that the receding front model has no first drying period.e first drying period of the diffusion model is longer (X cr � 0.0722) than in the case of the continuous model (X cr � 0.1687).Because the drying rate and moisture content computed by the three models are different (Figures 9 and 10), the drying times are also different.In order to investigate this difference further,  e results show that the drying time is a linear function of the sample size for the case of the diffusion model since in this considered case we have almost only outer resistance (due to small value of mass transfer coefficient β).If outside resistance is negligible (when mass transfer coefficient β goes to infinity), we should obtain drying time as a quadratic function of the sample size.It is found that the drying time increases more than linearly with sample size in the case of the continuous and receding front model.
As the numerical results of this example demonstrate, even in the case of isothermal drying, the framework presented here is a suitable choice.
In the third example, we consider the drying of an infinite plate with a thickness of 200 mm.e system is symmetric, and one-dimensional control volume elements can be used.In this example, the porosity is taken as ψ � 0.5, and the solid phase has thermal conductivity λ s � 1 W/m/K and volumetric heat capacity (ρc) s � 2•10 6 J/m 3 /K.e initial saturation and initial temperature are S 0 � 0.9 and T 0 � 20 °C.e drying air has pressure P ∞ � 1 bar, zero moisture content φ � 0, and temperature T ∞ � 80 °C.
e heat and mass transfer coefficients for the boundary layer are given as α � 95 W•m −2 •K −1 and β � 0.1 m•s −1 , respectively.ree cases with two types of pore-size distribution (mono-modal and bimodal) with different mean pore radii r 0 and different standard deviations σ 0 are investigated.
In this example, the capillary pressure (P c ), the absolute permeability (K), and the relative permeability of liquid (k w ) Figure 10: Evaporation rates for isothermal drying of a spherical particle using different drying models.and gas (k g ) phases are computed by using formulas (44) and ( 48)-(50), respectively.In order to show changes in transport parameters due to different pore-size distributions, information about the pore-size distributions and the corresponding absolute permeabilities is presented in Table 2 and Figure 12.In our simulation, the total void volume in each case of the 3 pore-size distributions is controlled with the help of formula (41).e total pore volume computed with the help of this formula must correspond to the porosity ψ � 0.5 of the sample (see also the comment after formula (41) above).Note that the capillary pressure (P c ) and the relative permeability of liquid (k ) and gas (k g ) are also computed as function of different pore-size distributions and are implemented directly in our calculation but are not shown here to save space.In order to see how changes in pore size and pore-size distributions lead to changes in relative permeability of liquid (k w ) and gas (k g ) phases and in capillary pressure (P c ), refer Figures 3 and 4 and Section 4 above.
e results of our simulation are presented in Figures 13  and 14.By comparing the three cases, it can be seen that, for the small pores with a narrow distribution, the first drying period is short.e large pore case with a broad distribution has a first drying period.e large pores in a bimodal distribution can significantly prolong the first drying period, but the second drying period is entirely determined by the small pores.It is observed that, in the case 2 with very large pores (1000 ± 100 nm), the drying process takes place faster than in the other two cases.In our simulation, the time to remove all the free water from the samples is 37.5 hours for case 1 (100 ± nm), 21.3 hours for case 2 (1000 ± 100 nm), and 27.3 hours for case 3 (100 ± 10; 200 ± 20 nm).For more discussion about the influence of pore-size distribution on drying kinetics, the readers are referred to [15].
In the fourth example, we consider the drying of the same plate as in the third example.However, we assume here isothermal condition.e simulation result obtained by our simulation is compared with the solution obtained by a completely different approach, namely, the "discrete approach" [8,9].e purpose here is to verify the accuracy of the simulation approach presented in this work.
In our simulations, convective drying by a flow of absolutely dry air at T ∞ � 20 °C and atmospheric pressure is applied.For T ∞ � 20 °C, the imaginary wet bulb temperature is 19.3 °C.e amount of adsorbed water is set to S irr � 1% to approach the case of no sorption in the discrete models.Formulas for monopore-size distribution presented in Section 4 are used to compute effective parameters of the drying model with three different distributions (r 0 ± σ 0 ): 100 ± 5 nm, 100 ± 10 nm, and 100 ± 25 nm.Similar to the third example, the capillary pressure (P c ), the absolute permeability (K), and the relative permeability of liquid (k w ) and gas (k g ) phases are computed by using formulas (44) and ( 48)-(50), respectively.More details on the simulation using the discrete approach can be found in [41].
e numerical results are presented in Figures 15 and 16 by plotting the normalized drying rate versus saturation and by plotting the saturation profiles at the end of the first drying period.Note that, in Figure 15, the normalized evaporation rate _ v is dimensionless and is defined as the ratio of the evaporation rate ( _ m v ) to the evaporation rate of the first drying period ( _ m v,I ): _ v � _ m v / _ m v,I .From the numerical results, it is observed that the continuous approach leads to a slightly longer first drying period and slightly flatter moisture profiles than the discrete approach.However, the overall behaviour of the drying curves demonstrates that the simulation using the continuous approach agrees well with the simulation using the discrete approach.e nonsmooth curves of the discrete approach are due to assumption that partially filled throats have no resistance to vapour diffusion.

Conclusion
In this work, we revisit the continuous approach for studying the drying problem of porous media.A ready-to-use system  International Journal of Chemical Engineering of governing equations is presented.
is system of equations can be directly solved by many numerical tools that are now available.By presenting the connection between the governing equations at the pore level and the ones at the continuum level, the framework presented here offers the possibility to modify the system of governing equations in the case of need, for example to reflect the different phenomena that are not yet taken into consideration.By using the concept of bundle of capillaries, we show that this concept can be used to investigate the Simulation results using the continuous approach (solid lines) and discrete approach (dashed lines) for an infinite plate with a thickness of 200 mm: normalized evaporation rates.
14 International Journal of Chemical Engineering influence of the microstructure on drying behaviour of porous media.

Figure 1 :
Figure 1: Modelling the drying process in porous media: pore scale (the microscopic scale) and macroscopic scale.

Figure 3 :Figure 4 :
Figure 3: Relative permeabilities for four cases of mono-and bimodal distributions.

Figure 5 :
Figure 5: Nonisothermal drying: temporal evolution of temperature (centre vs. surface temperature of the sphere particle).

Figure 6 :
Figure 6: Isothermal drying: temporal evolution of temperature of the sphere particle (the temperature of the whole sample is the same).

Figure 9 :
Figure 9: Evolution of average moisture content for isothermal drying of a spherical particle using different drying models.

FigureFigure 8 :
Figure Influence of sample size: drying time versus sample size.