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

Authors: Hong Thai Vu; Evangelos Tsotsas
Journal: International Journal of Chemical Engineering
 (2019)

---

## Abstract

The 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.


---

## Body

## 1. Introduction

The drying process plays an important role in many different industries, for example, in chemicals, pharmaceuticals, and agriculture. The 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]). The 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. The continuous model developed by Whitaker is considered as rigorous and the most advanced continuous model today. The 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. The 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.


## 2. Governing Equations at the Microscopic Scale (Pore Scale)

We consider here a rigid porous mediumV with external boundary ∂V in which the matrix is made of some solid material and a system of interconnected voids. The voids are also called here “pores”. These 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). The solid phase is denoted by s, the liquid phase (water) is denoted by w, and the gas phase is denoted by g. The 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. These conservation laws will be presented in the following, according to the work of Whitaker [18].Figure 1
Modelling the drying process in porous media: pore scale (the microscopic scale) and macroscopic scale.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. The equations presented in this section are valid for each phase (liquid, air, and vapour) insidedV. 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 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(1)∂ρ∂t+∇⋅ρv=0,and the mass conservation equation of each component in each phase is(2)∂ρi∂t+∇⋅ρivi=0,fori=1,2,…,N.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.(3)v=∑i=1Nρiρvi,ρ=∑i=1Nρi,where N denotes the total number of components, ρi is the density, and vi is the velocity of component i.


### 2.2. 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(4)ρDvDt=∇⋅T,where T is the stress tensor. The conservation of angular momentum requires this tensor to be symmetric:(5)T=TT.


### 2.3. Energy Conservation at the Pore Scale

The conservation of energy for each phase can be presented in the following format:(6)∂∂tρh+∇⋅ρhv=−∇⋅q+DPDt+τ:∇v+Φ,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. The conductive heat flux vector q is computed by Fourier’s law:(7)q=−λ∇T,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:(8)τ:∇v=0,DPDt=0,and then the energy equation is reduced to(9)∂∂tρh+∇⋅ρhv=−∇⋅q.It is also assumed that the enthalpy is independent of pressure and that all heat capacities are constant so that the following relationship(10)h=cpT−TR,holds, where cp is the specific heat capacity and TR is the reference temperature (chosen here to be TR = 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.


### 2.4. Solid Phase

The solid phase is considered to be rigid and fixed in space with zero velocity (movement of the solid phase is not considered):(11)vs=0.This assumption means that, for the solid phase, we need to study only the conservation of energy equation (6), which now becomes(12)ρs∂hs∂t=−∇⋅qs.By making use of equation (7) and taking into account the assumption (10), we have the following conservation of energy for the solid phase:(13)ρscp,s∂Ts∂t=λs∇2Ts.


### 2.5. Liquid Phase

For the liquid phase, which contains water as the only component, the mass conservation equation is(14)∂ρw∂t+∇⋅ρwvw=0.The use of equations (7) and (10) allows us to write the energy conservation of the liquid phase in the following format:(15)ρwcp,w∂Tw∂t+vw⋅∇Tw=λw∇2Tw.


### 2.6. Gas Phase

The gas phase is more complicated than the solid and the liquid phase since it contains two components: air and vapour. The mass conservation of the gas phase can be written in the following form:(16)∂ρg∂t+∇⋅ρgvg=0.By writing the component velocityvi in terms of the mass average velocity vg and the diffusion velocity ui,(17)vi=vg+ui,i=a,v.We can write the mass conservation of air and vapour in the following form:(18)∂ρi∂t+∇⋅ρivg=−∇⋅ρiui,i=a,v.Furthermore, by expressing the diffusive fluxρiui as(19)ρiui=−ρgδv,a∇ρiρg,where δv,a is the binary molecular diffusion coefficient for vapour and air, we finally have(20)∂ρi∂t+∇⋅ρivg=∇⋅ρgδv,a∇ρiρg,i=a,v.Note that for a multicomponent phase, the appropriate form of the energy equation (9) can be written in the following format:(21)∂∂t∑i=1Nρihi+∇⋅∑i=1Nρihivi=−∇⋅q,where hi 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:(22)h=∑i=1Nρiρhi.By using equation (21), for the gas phase, we have(23)ρgcp,g∂Tg∂t+vg⋅∇Tg=λg∇2Tg−∇⋅ρahaua+ρvhvuv,where(24)cp,g=ρacp,a+ρvcp,vρg.In addition to the above conservation equations, ideal gas laws are assumed for partial and total gas pressures:(25)Pi=ρiR˜TM˜i,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. The constraint for the partial and total gas pressure is simply(26)Pa+Pv=Pg.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].


## 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(1)∂ρ∂t+∇⋅ρv=0,and the mass conservation equation of each component in each phase is(2)∂ρi∂t+∇⋅ρivi=0,fori=1,2,…,N.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.(3)v=∑i=1Nρiρvi,ρ=∑i=1Nρi,where N denotes the total number of components, ρi is the density, and vi is the velocity of component i.


## 2.2. 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(4)ρDvDt=∇⋅T,where T is the stress tensor. The conservation of angular momentum requires this tensor to be symmetric:(5)T=TT.


## 2.3. Energy Conservation at the Pore Scale

The conservation of energy for each phase can be presented in the following format:(6)∂∂tρh+∇⋅ρhv=−∇⋅q+DPDt+τ:∇v+Φ,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. The conductive heat flux vector q is computed by Fourier’s law:(7)q=−λ∇T,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:(8)τ:∇v=0,DPDt=0,and then the energy equation is reduced to(9)∂∂tρh+∇⋅ρhv=−∇⋅q.It is also assumed that the enthalpy is independent of pressure and that all heat capacities are constant so that the following relationship(10)h=cpT−TR,holds, where cp is the specific heat capacity and TR is the reference temperature (chosen here to be TR = 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.


## 2.4. Solid Phase

The solid phase is considered to be rigid and fixed in space with zero velocity (movement of the solid phase is not considered):(11)vs=0.This assumption means that, for the solid phase, we need to study only the conservation of energy equation (6), which now becomes(12)ρs∂hs∂t=−∇⋅qs.By making use of equation (7) and taking into account the assumption (10), we have the following conservation of energy for the solid phase:(13)ρscp,s∂Ts∂t=λs∇2Ts.


## 2.5. Liquid Phase

For the liquid phase, which contains water as the only component, the mass conservation equation is(14)∂ρw∂t+∇⋅ρwvw=0.The use of equations (7) and (10) allows us to write the energy conservation of the liquid phase in the following format:(15)ρwcp,w∂Tw∂t+vw⋅∇Tw=λw∇2Tw.


## 2.6. Gas Phase

The gas phase is more complicated than the solid and the liquid phase since it contains two components: air and vapour. The mass conservation of the gas phase can be written in the following form:(16)∂ρg∂t+∇⋅ρgvg=0.By writing the component velocityvi in terms of the mass average velocity vg and the diffusion velocity ui,(17)vi=vg+ui,i=a,v.We can write the mass conservation of air and vapour in the following form:(18)∂ρi∂t+∇⋅ρivg=−∇⋅ρiui,i=a,v.Furthermore, by expressing the diffusive fluxρiui as(19)ρiui=−ρgδv,a∇ρiρg,where δv,a is the binary molecular diffusion coefficient for vapour and air, we finally have(20)∂ρi∂t+∇⋅ρivg=∇⋅ρgδv,a∇ρiρg,i=a,v.Note that for a multicomponent phase, the appropriate form of the energy equation (9) can be written in the following format:(21)∂∂t∑i=1Nρihi+∇⋅∑i=1Nρihivi=−∇⋅q,where hi 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:(22)h=∑i=1Nρiρhi.By using equation (21), for the gas phase, we have(23)ρgcp,g∂Tg∂t+vg⋅∇Tg=λg∇2Tg−∇⋅ρahaua+ρvhvuv,where(24)cp,g=ρacp,a+ρvcp,vρg.In addition to the above conservation equations, ideal gas laws are assumed for partial and total gas pressures:(25)Pi=ρiR˜TM˜i,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. The constraint for the partial and total gas pressure is simply(26)Pa+Pv=Pg.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].


## 3. Numeric Ready Continuum Equations of the Drying Process

The macroscopic transport equations can be obtained from the microscopic ones presented in Section2 by using the volume averaging method [17]. This set of macroscopic equations is summarized as follows [14].


### 3.1. Conservation Equation for Water

The mass conservation equation for water in the both liquid and gas phase can be written as follows:(27)∂∂tρwεw+εgρv+∇⋅ρwvw+ρvvg=∇⋅ρgDeff⋅∇ρvρg,where εw=Vw/V and εg=Vg/V are the volume fractions of the liquid and gas phase with respect to the total volume V of the porous medium and Deff is the effective diffusivity tensor.


### 3.2. Conservation Equation for Air

The second equation is the mass conservation equation for air in the gas phase:(28)∂∂tεgρa+∇⋅ρavg=∇⋅ρgDeff⋅∇ρaρg.


### 3.3. Conservation Equation of Energy

The third equation is the energy conservation equation for the whole system:(29)∂∂tεsρshs+εwρwhw+εgρvhv+εgρaha+∇⋅ρwhwvw+ρvhv+ρahavg=∇⋅ρghaDeff⋅∇ρaρg+∇⋅ρghvDeff⋅∇ρvρg+∇⋅λeff∇T,where εs=Vs/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 Vs+Vw+Vg=V, the sum of all volume fractions is unity:(30)εs+εw+εg=1.In the above equations, the velocities are computed by using(31)vw=−Kwkwμw∇Pw−ρwg∇χ,vg=−Kgkgμg∇Pg−ρgg∇χ,where Kw and Kg are the absolute permeabilities of liquid and gas phases, kw and kg 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):(32)Pg=Pa+Pv=ρaR˜TM˜a+Pv,where the vapour pressure will be computed by the sorption isotherm, which can be determined by an experiment. The liquid water pressure is computed by the following equation:(33)Pw=Pg−Pc,where Pc is the capillary pressure.Also in the above equations, the enthalpies of the solid phase (hs), liquid water (hw), air (ha), and vapour (hv) are computed using (10) as(34)hi=ciT−TR,i=s,w,a,hv=Δhvw+cvT−TR,where Δhvw is the latent heat of vaporization. For the computation of these values, we use here cw=4185J⋅kg−1⋅K−1, cv=1874J⋅kg−1⋅K−1, and ca=1005.683J⋅kg−1⋅K−1. The value of cs 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 cs, we can use an averaged heat capacity ρCp¯ defined as follows:(35)ρCp¯=εsρscs+εwρwcw+εgρvcv+εgρaca.In solving the above system of conservation equations, we can use as main variables the temperatureT, 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. The formula for moisture content X is written as(36)X=εwρwεsρs.By using moisture content, the volume fraction of liquid waterεw can be computed by equation (36) as εw=εsρs/ρwX. The volume fraction of gas is(37)εg=1−εs−εw=1−εs−εsρsρwX,where the volume fraction of the solid phase εs can be computed by the porosity of the medium under consideration. The vapour pressure will be computed by the sorption isotherm and can be determined experimentally as the function of moisture content and temperature Pv=PvX,T. The vapour density is ρv=PvM˜v/R˜T. The gas density is ρa=ρg−ρv. Note that for simplicity, the density of water is considered here as constant ρw=1000kg/m3.Besides the above governing equations, the boundary conditions for mass and heat transfer at the external drying surfaces∂V 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 Jw⋅n^ at ∂V, we have(38)Jw⋅n^=m˙v∂V=βPgM˜vR˜TlnPg−Pv,∞Pg−Pv,and for the energy flux Je⋅n^ at ∂V,(39)Je⋅n^=q˙+Δhvm˙v∂V=αT−T∞+ΔhvβPgM˜vR˜TlnPg−Pv,∞Pg−Pv,where Jw and Je are the fluxes of water and heat, respectively; Pv,∞ 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:(40)Pg=P∞.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. The main variables to be solved are the temperature T, the moisture content X, and the gas density ρg. The boundary condition for heat and mass transfer will be temperature and relative humidity of drying air. These conditions will be used with the help of equations (38)–(40) in terms of drying air temperature (T∞), vapour pressure in drying air (Pv,∞), and total pressure of drying air (P∞). The initial condition will be the temperature T, moisture content X, and gas pressure (Pg) 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.


## 3.1. Conservation Equation for Water

The mass conservation equation for water in the both liquid and gas phase can be written as follows:(27)∂∂tρwεw+εgρv+∇⋅ρwvw+ρvvg=∇⋅ρgDeff⋅∇ρvρg,where εw=Vw/V and εg=Vg/V are the volume fractions of the liquid and gas phase with respect to the total volume V of the porous medium and Deff is the effective diffusivity tensor.


## 3.2. Conservation Equation for Air

The second equation is the mass conservation equation for air in the gas phase:(28)∂∂tεgρa+∇⋅ρavg=∇⋅ρgDeff⋅∇ρaρg.


## 3.3. Conservation Equation of Energy

The third equation is the energy conservation equation for the whole system:(29)∂∂tεsρshs+εwρwhw+εgρvhv+εgρaha+∇⋅ρwhwvw+ρvhv+ρahavg=∇⋅ρghaDeff⋅∇ρaρg+∇⋅ρghvDeff⋅∇ρvρg+∇⋅λeff∇T,where εs=Vs/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 Vs+Vw+Vg=V, the sum of all volume fractions is unity:(30)εs+εw+εg=1.In the above equations, the velocities are computed by using(31)vw=−Kwkwμw∇Pw−ρwg∇χ,vg=−Kgkgμg∇Pg−ρgg∇χ,where Kw and Kg are the absolute permeabilities of liquid and gas phases, kw and kg 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):(32)Pg=Pa+Pv=ρaR˜TM˜a+Pv,where the vapour pressure will be computed by the sorption isotherm, which can be determined by an experiment. The liquid water pressure is computed by the following equation:(33)Pw=Pg−Pc,where Pc is the capillary pressure.Also in the above equations, the enthalpies of the solid phase (hs), liquid water (hw), air (ha), and vapour (hv) are computed using (10) as(34)hi=ciT−TR,i=s,w,a,hv=Δhvw+cvT−TR,where Δhvw is the latent heat of vaporization. For the computation of these values, we use here cw=4185J⋅kg−1⋅K−1, cv=1874J⋅kg−1⋅K−1, and ca=1005.683J⋅kg−1⋅K−1. The value of cs 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 cs, we can use an averaged heat capacity ρCp¯ defined as follows:(35)ρCp¯=εsρscs+εwρwcw+εgρvcv+εgρaca.In solving the above system of conservation equations, we can use as main variables the temperatureT, 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. The formula for moisture content X is written as(36)X=εwρwεsρs.By using moisture content, the volume fraction of liquid waterεw can be computed by equation (36) as εw=εsρs/ρwX. The volume fraction of gas is(37)εg=1−εs−εw=1−εs−εsρsρwX,where the volume fraction of the solid phase εs can be computed by the porosity of the medium under consideration. The vapour pressure will be computed by the sorption isotherm and can be determined experimentally as the function of moisture content and temperature Pv=PvX,T. The vapour density is ρv=PvM˜v/R˜T. The gas density is ρa=ρg−ρv. Note that for simplicity, the density of water is considered here as constant ρw=1000kg/m3.Besides the above governing equations, the boundary conditions for mass and heat transfer at the external drying surfaces∂V 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 Jw⋅n^ at ∂V, we have(38)Jw⋅n^=m˙v∂V=βPgM˜vR˜TlnPg−Pv,∞Pg−Pv,and for the energy flux Je⋅n^ at ∂V,(39)Je⋅n^=q˙+Δhvm˙v∂V=αT−T∞+ΔhvβPgM˜vR˜TlnPg−Pv,∞Pg−Pv,where Jw and Je are the fluxes of water and heat, respectively; Pv,∞ 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:(40)Pg=P∞.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. The main variables to be solved are the temperature T, the moisture content X, and the gas density ρg. The boundary condition for heat and mass transfer will be temperature and relative humidity of drying air. These conditions will be used with the help of equations (38)–(40) in terms of drying air temperature (T∞), vapour pressure in drying air (Pv,∞), and total pressure of drying air (P∞). The initial condition will be the temperature T, moisture content X, and gas pressure (Pg) 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.


## 4. 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. The effective diffusivity and effective thermal conductivity are considered as dependent on porosity ψ and saturation S. These 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. The 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. The 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. There 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.Figure 2
Partially saturated bundle of capillaries: relationship between free water saturationSfw and maximum radius filled rfill [32].For our investigation, two types of pore-size distributions (capillary radius distributions) are used. The first one is a monomodal normally distributed pore volume(41)dVdr=C2πσ0e−1/2r−r0/σ02,where r0 is the mean pore radius, σ0 is the standard deviation, and C is a constant. This normal distribution is truncated at r0±2.5σ0. The 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.The second distribution type is a bimodal distribution, which consists of two monomodal models (with different mean pore radii and deviations) named here assmall 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 saturationSfw, small capillaries are filled up to a maximum radius rfill such that(42)Sfw=∫rminrfilldVdrdr,where rmin is the smallest capillary radius of the bundle.The localization of free waterSfw is the key for computing effective parameters as the function of saturation. Based on equation (42), the maximum radius rfill can be computed for a given Sfw. For any saturation under the irreducible value, the saturation of free water Sfw is zero and the maximum radius filled by liquid rfill is set to rmin.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:(43)φ=SSirr⋅2−SSirr,forS≤Sirr,where φ is the relative humidity, Sirr is the maximum amount of adsorbed water, S=Sirr+Sfw, 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:(44)PcSfw,T=2⋅σTrfillSfw,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:(45)V˙=18η⋅ΔPwL⋅π⋅r4,where L is the capillary length, r is the capillary radius, η is the dynamic viscosity (temperature dependent), and Pw 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 (kw=1), we obtain(46)v¯=Kη·ΔPwL.By comparing equations (45) and (46), the absolute permeability can be found to be(47)K=18r2.By extending the above formula to the bundle of capillaries, we get(48)K=18∫rminrmaxr2dVdrdr,where the interval rmin,rmax is the total range of the pore-size distribution.The relative permeabilities of liquid and gas phases are computed in a similar way as the absolute permeability:(49)kwSfw=18K∫rminrfillr2dVdrdr,(50)kgSfw=18K∫rfillrmaxr2dVdrdr.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]. The 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. The capillary pressure curves for these four pore-size distributions are illustrated in Figure 4. Naturally, with the decreasing saturation, the capillary pressure increases. The overall level of the capillary pressure is determined by the mean pore size r0 of the mode, and its range of variation is determined by the standard deviation σ0.Table 1
Pore-size distributions and absolute permeabilities.
r
0
±
σ
0 (nm)
K (10−15 m2)
Case 1
100 ± 5
1.238
Monomodal
Case 2
1000 ± 100
124.6
Case 3
100 ± 10; 200 ± 20
3.119
Bimodal
Case 4
100 ± 10; 2000 ± 200
245.1Figure 3
Relative permeabilities for four cases of mono- and bimodal distributions.
(a)
(b)Figure 4
Capillary pressure for four cases of mono- and bimodal distributions.The effective vapour diffusivity does not depend on the distribution of the pores but on the evaporation area. Therefore, this transport parameter is assumed to be a function of saturation, porosity, and the binary diffusion coefficientδva [14](51)DeffS,T,P=1−S⋅ψ⋅δvaT,P,where ψ is the porosity, which is the ratio of the total void or pore volume (Vvoid) to the total volume (Vtotal) of the material ψ=Vvoid/Vtotal.The 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)λeffS,T=1−ψ⋅λsT+S⋅ψ⋅λwT.


## 5. 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]. The method was 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]. The 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]. The 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. The 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. The 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. The first example is presented to demonstrate that the current approach can capture important effects which lead to large differences between isothermal and nonisothermal drying. The second example shows that simplification in models like the diffusion model and receding front model can lead to large differences in drying kinetics. This example also shows that the current approach can capture the different drying characteristics of the two simplified models. The 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). The fourth example also demonstrates the accuracy of the current approach when 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. The initial temperature of the sample is T0 = 20°C, the initial moisture content is X0 = 1, and the initial pressure is P0 = 1 bar for the whole sample. The heat transfer coefficient is α = 14.25 W·m−2·K−1, and mass transfer coefficient is β = 0.015 m·s−1. The drying air has temperature T∞ = 20°C and relative humidity φ = 50%.The material properties of light concrete used in this example are given as follows [24]: porosity ψ = 0.8, solid density ρs=2500kg⋅m−3, and heat capacity:(53)ρCp¯=εsρs840+4185XJ⋅kg−1⋅K−1.The fully saturated material has a moisture content ofXsat=1.6. The sorption isotherm is(54)ϕX,T=PvPv∗T=1,ifX>Xirr,XXirr2−XXirr,ifX≤Xirr,where Xirr=0.07 is the irreducible moisture content. The saturation vapour pressure Pv∗T is given as a function of temperature T (°C) by Antoine’s equation:(55)Pv∗=133.32⋅e18.584−3984.2/233.426+T.The capillary pressure is computed from(56)Pc=40⋅σT⋅e8.4057×10−0.3476Xfw,where the surface tension σ is a function of temperature T (°C):(57)σT=−1.3⋅10−7⋅T2−1.58⋅10−4⋅T+0.07606N⋅m−1,and Xfw is the moisture content of free water: Xfw=X−Xirr.The absolute permeability is taken as constant:K = 2 × 10−13 m2. The relative permeabilities for liquid kw and for gas kg phases are calculated from the following relationship:(58)kw=0,X>Xirr,Sfw3,X≤Xirr,kg=1,X>Xirr,1+2Sfw−3Sfw2,X≤Xirr,where Sfw=X−Xirr/Xsat−Xirr and Xsat is the saturated moisture content.The effective diffusivity is calculated from(59)Deff=0.2⋅δva⋅kg,where kg is the relative permeability of gas and δva is the binary diffusivity of vapour in air [31]:(60)δvaT,P=2.26⋅10−5⋅TTR1.81PRPg,where TR=273.15K and PR=101325Pa are reference temperature and pressure, respectively.The 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(61)λeff=0.142+0.46XW⋅m−1⋅K−1.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 thisisothermal 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. The 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. The 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. The 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 Section3 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/m2/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 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 Section3). The 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. The numerical results are presented in Figures5–7. The 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 T0 = 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 T0 = 20°C in the second drying period. The 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.Figure 5
Nonisothermal drying: temporal evolution of temperature (centre vs. surface temperature of the sphere particle).Figure 6
Isothermal drying: temporal evolution of temperature of the sphere particle (the temperature of the whole sample is the same).Figure 7
Isothermal and nonisothermal drying: drying rate curves.From the results presented in Figure7, the evaporation rate in the first period and the critical moisture content for the isothermal case are m˙v,I=0.1310g⋅m−2⋅s−1 and Xcr=0.1538. These 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. This 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).The 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. The 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 receding front model presented hereafter. The 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 X0 = 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:(62)∂Xr,t∂t=1r2∂∂rr2δ∗∂X∂r,where r is the radial coordinate and the diffusion coefficient δ∗ takes the constant value of 2.6 · 10−5 m2·s−1 [31].In addition to the diffusion equation above, we need to specify the boundary conditions for a sphere at radiusr=0 and r=R. At the centre of the sphere, due to symmetry, we have(63)∂X∂rr=0=0.At the boundary of sphere, the mass flux must satisfy the boundary condition (38). The mass flux at the boundary r=R is computed as follows [31]:(64)m˙vr=R=−δ∗ρsdXdrr=R,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), i.e., m˙vr=R=Jw⋅n^, we have(65)−δ∗ρsdXdrr=R=βPgM˜vR˜T⋅lnPg−Pv,∞Pg−Pv,in which the vapour pressure Pv 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. The 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(66)m˙v=−1β+sδ∗−1PgM˜vR˜T⋅lnPg−Pv,∞Pg−Pv∗,where β is the mass transfer coefficient at the surface and z=s is the position of the front at time t. The 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. This 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:(67)t=ρw⋅ΔVA⋅ψ⋅m˙v,where A is the surface area of the sphere, ΔV is the volume of the shell, and ψ is the porosity.Figure 8
Receding front model for a plate.In using the receding front model, for every given value ofs 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 Xav 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: Xav,t,m˙v.The different solutions of the drying problem by using the three models are presented and compared in Figures9–11. As we can see from Figure 9, all three models capture the initial value of evaporation rate m˙v,I=0.262g⋅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. The first drying period of the diffusion model is longer (Xcr = 0.0722) than in the case of the continuous model (Xcr = 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, different radii R of the sample are examined and the corresponding drying times are computed. As result, the drying time versus radius of the sample is plotted in the logarithmic scale in Figure 11. Note that the simulation results presented in Figures 9 and 10 are obtained with radius R = 2.5 mm, whereas in Figure 11, the results are obtained with radii R = 1 mm–15 mm. The 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.Figure 9
Evolution of average moisture content for isothermal drying of a spherical particle using different drying models.Figure 10
Evaporation rates for isothermal drying of a spherical particle using different drying models.Figure 11
Influence of sample size: drying time versus sample size.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. The 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·106 J/m3/K. The initial saturation and initial temperature are S0 = 0.9 and T0 = 20°C. The drying air has pressure P∞ = 1 bar, zero moisture content φ = 0, and temperature T∞ = 80°C. The 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. Three cases with two types of pore-size distribution (mono-modal and bimodal) with different mean pore radii r0 and different standard deviations σ0 are investigated.In this example, the capillary pressure (Pc), the absolute permeability (K), and the relative permeability of liquid (kw) and gas (kg) 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). The 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).Table 2
Absolute permeability as function of pore-size distributions.
r
0
±
σ
0 (nm)
K (10−15 m2)
Case 1
100 ± 5
1.238
Monomodal
Case 2
1000 ± 100
124.6
Case 3
100 ± 10; 200 ± 20
3.119
BimodalFigure 12
Monomodal pore-size distribution ((a) 100 ± 5 nm) and bimodal pore-size distribution ((b) 100 ± 10 nm and 200 ± 20 nm).
(a)
(b)Note that the capillary pressure (Pc) and the relative permeability of liquid (kw) and gas (kg) 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 (kw) and gas (kg) phases and in capillary pressure (Pc), refer Figures 3 and 4 and Section 4 above.The results of our simulation are presented in Figures13 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. The large pore case with a broad distribution has a long first drying period. The 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 ± 5 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].Figure 13
Saturation profiles from left to right: (a) case 1 (100 ± 5 nm); (b) case 2 (1000 ± 100 nm); (c) case 3 (100 ± 10; 200 ± 20 nm).
(a)
(b)
(c)Figure 14
Drying curves for different pore-size distributions.In the fourth example, we consider the drying of the same plate as in the third example. However, we assume here isothermal condition. The simulation result obtained by our simulation is compared with the solution obtained by a completely different approach, namely, the “discrete approach” [8, 9]. The 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 atT∞ = 20°C and atmospheric pressure is applied. For T∞ = 20°C, the imaginary wet bulb temperature is 19.3°C. The amount of adsorbed water is set to Sirr = 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 (r0 ± σ0): 100 ± 5 nm, 100 ± 10 nm, and 100 ± 25 nm. Similar to the third example, the capillary pressure (Pc), the absolute permeability (K), and the relative permeability of liquid (kw) and gas (kg) 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].The numerical results are presented in Figures15 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. The nonsmooth curves of the discrete approach are due to assumption that partially filled throats have no resistance to vapour diffusion.Figure 15
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.Figure 16
Simulation results using the continuous approach (solid lines) and discrete approach (dashed lines) for an infinite plate with a thickness of 200 mm: saturation profiles at the end of the first drying period.


## 6. Conclusion

In this work, we revisit the continuous approach for studying the drying problem of porous media. A ready-to-use system of governing equations is presented. This 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 influence of the microstructure on drying behaviour of porous media.



---

*Source: 9043670-2019-04-01.xml*