Chemical Engineering Science 63 (2008) 5218 -- 5228 Contents lists available at ScienceDirect Chemical Engineering Science journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / c e s Influence of heating mode on drying behavior of capillary porous media: Pore scale modeling V.K. Surasani, T. Metzger ∗ , E. Tsotsas Thermal Process Engineering, Otto-von-Guericke-University, Universitaetsplatz 2, P.O. 4120, 39016 Magdeburg, Germany A R T I C L E I N F O Article history: Received 4 February 2008 Received in revised form 1 July 2008 Accepted 5 July 2008 Available online 12 July 2008 Keywords: Pore networks Phase distributions Temperature fields Bimodal pore size distribution Invasion percolation Imbibition A B S T R A C T Invasion percolation (IP) rules under non-isothermal conditions are applied to model the pore-scale events occurring during drying of capillary porous media, namely displacement of immiscible phases and cluster formation. A saturated two-dimensional network with a bimodal pore size distribution is dried by applying two different heat transfer boundary conditions; one corresponds to convective drying and the other to less resistive contact drying. Simulated macroscopic drying behavior is presented in conjunction with freely evolved microscopic temperature fields and phase distributions for both heating modes. Convective heating exhibits similar invasion patterns as those in isothermal simulations; both are dominated by the spatial distribution of pore radii. However, in contact heating, temperature dependency of surface tension produces significantly different invasion patterns. © 2008 Elsevier Ltd. All rights reserved. 1. Introduction Drying of capillary porous media has been of significant scientific interest for many researchers since it depicts a complex, coupled, and multiphase process with a wide range of applications in industry. Because of its high cost in energy, drying is an operation with a high potential for optimization with respect to energy savings. In addition, drying plays an important role in weathering of monuments, buildings, and rocks. For many years, it has been studied experimentally for measuring drying kinetics on the macro-scale and the most advanced techniques like magnetic resonance imaging (MRI) (Guillot et al., 1989; Pel et al., 1996), and neutron scattering (Pel et al., 1993) are used to measure the spatial liquid distributions in porous media. However, in most of the cases, theoretically, drying is a remarkably little understood process. Numerous theoretical models that have been proposed in the past are based on rigorous heat, mass, and momentum balance equations on a fictitious continuum. The methods, which are mostly used in such a continuum approach, are homogenization (Sanchez-Palencia, 1980) or volume averaging techniques (Quintard and Whitaker, 1993). Even though much knowledge has been gained through this macroscopic approach to drying, and validated with experiments, it fails in relating the macroscopic phenomena to the microscopic, i.e. pore-scale, events during the drying process. These events refer to the immiscible displacement of the defending phase (wetting fluid) by the invading phase (non-wetting ∗ Corresponding author. Tel.: +49 391 67 11362; fax: +49 391 67 11160. E-mail address: thomas.metzger@ovgu.de (T. Metzger). 0009-2509/$ - see front matter © 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2008.07.011 fluid), and formation of dry and wet patches (clusters) in the pore space caused by interfacial effects on the phase boundaries. Discrete network models based on fractal concepts and statistical physics such as percolation theory, and their applications to drying of porous media, mainly reviewed by Prat (1993) and Yiotis et al. (2001), are particularly suited for studying the situations where the continuum approach fails (Laurindo and Prat, 1998). A network model is a comprehensive model which incorporates the description of pore space and the physics of the pore-scale events (Berkowitz and Ewing, 1998). The basic understanding of pore network modeling commences with the geometric representation of the void space in the porous medium. Most of the literature work has adopted a network representation in the form of sites (pores) connected with bonds (throats) (Nowicki et al., 1992; Prat, 1993; Yiotis et al., 2001). Prat (1993) showed that the immiscible displacements in the drying process can be treated as a special case of invasion percolation (IP). Menisci in the pore space will respond to capillary, viscous, and buoyancy forces, and the interplay between these forces has been studied in literature. In the absence of viscous forces, gravity or thermal gradients, phase distributions (IP patterns) completely depend on the spatial disorder of pore size (Laurindo and Prat, 1998). In the presence of gravity, the phase distributions are patterns of IP in a stabilizing gradient (IPSG) or in a destabilizing gradient (IPDG) depending on the direction of the gravity vector to the main direction of vapor flow (Laurindo and Prat, 1998). Similarly, IPSG patterns can be observed in the presence of viscous forces (Yiotis et al., 2001; Metzger et al., 2007a). Additionally, the influence of pore structure on drying kinetics has been investigated using pore network modeling V.K. Surasani et al. / Chemical Engineering Science 63 (2008) 5218 -- 5228 (Metzger et al., 2007b) for various two- and three-dimensional networks. For a comprehensive record on the use of pore network models in drying refer to the reviews by Prat (2002) and Metzger et al. (2007c). In the above-mentioned literature, invasion is considered at slow drying rates. Heat transfer aspects of drying are neglected when studying the physical and structural influences on IP patterns. The influence of temperature gradients—which were imposed onto the drying product—has been investigated by Plourde and Prat (2003), Huinink et al. (2002). However, product temperature did not evolve freely as corresponding to real drying boundary conditions. Recently, Surasani et al. (2008a) considered heat transfer in a pore network model that assumes dominant capillary forces and neglects both viscous forces and gravity. This network model was the first to describe the free evolution of temperature fields in convective drying. The present work applies an improved model version to the contact heating mode and simulates the evolution of invasion patterns and temperature fields, which are then compared with those of the convective heating mode (Surasani et al., 2008a). The concept of condensation in clusters, partially implemented in the algorithm, is discussed in this article. ∞ 5219 ∞ Fig. 1. Partially dried pore network with boundary layer and several liquid clusters, definitions of pore and throat conditions and magnified control volume i (CV). 2. Pore network model Pore network models are unique in accounting both for a geometrical description of the pore space and for the physics of the pore-scale events, i.e. transport mechanisms at the pore scale. In this section, we will present the data structure model, i.e. the geometry of the porous network, boundary layer formation, and pore and throat conditions. Further, the model equations that represent quasi-stationary mass transfer relations in each phase and transient heat transfer are presented. The coupling relations between heat and mass transfer are also discussed in this section. 2.1. Data structure model Void space in the porous medium is envisaged as a network of cylindrical throats connected by pore nodes. In general, the connectivity of pores, i.e. co-ordination number, will determine the topology of void space in two- and three-dimensional domain. Even though the implemented algorithm is independent of connectivity, we stick to a two-dimensional square network (see Fig. 1) due to computational limitations of the present algorithm Surasani et al. (2008a). A pictorial representation of the solid network and the gas-side boundary layer flow conditions is given in Fig. 1. Radii of the throats are distributed randomly. The pores are volume less; each pore is associated with a control volume (CV). To describe the first drying period with its constant drying rate, boundary layer is modeled in a discretized way, by extending the network into the gas phase with additional nodes (see Fig. 1); their number is given by mass transfer coefficient and vapor diffusivity (Metzger et al., 2005; Irawan et al., 2005; Yiotis et al., 2006): =  = NBL L.  (1) The mass transfer coefficient  can be obtained from an empirical correlation as Sh = LT = f (Re, Sc),  (2) with Reynolds number Re = uLT / g and Schmidt number Sc = / g . Here,  is vapor diffusivity, u bulk air velocity, LT network length, and g kinematic viscosity of the gas. However, Eq. (1) is derived for a plate at uniform partial pressure which is not necessarily fulfilled for drying porous media (see Masmoudi and Prat, 1991). Heat transfer between bulk air and network surface is described by a heat transfer coefficient , which can be obtained from the analogy of heat and mass transfer as = g Pr  Sc  1/3 , (3) where g is gas thermal conductivity and Pr Prandtl number. Immiscible displacements upon evaporation of liquid will cause evolving IP patterns in the network. One such pattern is shown in Fig. 1 which represents the continuous vapor phase and disconnected liquid patches (clusters) in the solid network. Pore and throat conditions (Surasani et al., 2008a) (Fig. 1) will determine the transport within and between the phases. At the phase boundary, where the menisci exist, liquid flow, phase change, and vapor diffusion are simultaneous. In interfacial throats (full or partially filled), the phase boundary is either a moving meniscus or a stationary meniscus as evaporated liquid diffuses through air and, at the same time, liquid flows from or to the meniscus depending upon boundary liquid pressures. In gas throats, diffusion takes place due to vapor pressure gradients; and in liquid throats, capillary flow takes place due to capillary pressure gradients in the cluster. Pore saturations are given as the mean of the saturations of connected throats. Equilibrium vapor pressure Pv∗ (T) exists in partially filled pores, giving boundary conditions to the linear system of vapor transport. The missing boundary condition comes from partial vapor pressure (Pv,∞ ) of convective air. 2.2. Drying model Initially, the network is saturated completely with water and has temperature Tin =20 ◦ C, whereas drying air has zero absolute humidity and temperature T∞ . Since we only investigate situations well below the boiling temperature, gradients in total pressure P are neglected; constant total pressure is applied in the whole gas phase (boundary layer and network). The number of boundary layer nodes formed due to convective air is assumed to be constant. Transport mechanisms within and in between the phases are given by pore and throat conditions. Vapor diffusion through the air is assumed to be quasi-stationary and the displacements of the immiscible phases act only upon capillary forces, i.e. viscous effects are neglected in both phases. In each CV local thermal equilibrium is applied and heat transfer in the network is only by conduction. 5220 V.K. Surasani et al. / Chemical Engineering Science 63 (2008) 5218 -- 5228 2.2.1. Vapor transport The vapor flow rate Ṁv,ij between two gas pores is given by Stefan's law. From the law of conservation and the assumption of quasi-stationary vapor transport, total mass flow rates from any gas pore are equal to zero: 4  Ṁv,ij = − j=1 4  j=1 Aij  PM̃v L R̃T ln  P − Pv,i P − Pv,j  = 0. (4) Here, Aij is the cross-sectional area of the throat in the solid network (Aij = rij2 ). Since we neglect the overlap of throats at pore nodes, all throats in the network are assumed as equal in length (Lij = L). Vapor transfer in the boundary is assumed as molecular diffusion. Thus, Eq. (4) is equally valid to represent the vapor transport in the boundary layer; in this case Aij in Eq. (4) is the area of cross section between the two pores in the boundary layer (Aij = LW). The first boundary condition to Eq. (4) is equilibrium vapor pressure in partially filled gas pores (Pv,j = Pv∗ (T)). As second boundary condition, vapor pressures in bulk fluid, are applied to the gas pores at the edge of the boundary layer (Pv,j = Pv,∞ ). The system of equations for vapor transport (Eq. (4)) for gas pores both in network and boundary layer can be put in linear form to be solved for partial vapor pressures (Pv,i ) using the conjugate gradient method (CGM). The partially filled pores should be treated differently since they are the boundaries of the clusters/single throats and the vapor flow rates will play an important role in calculating the cluster evaporation rates. The right-hand side of Fig. 1 represents partially filled pore i (CVi ) with neighbor pores named as 1, 2, 3 at unknown vapor pressures (Pv,j ) and the remaining pore is at liquid pressure (Pl,4 ). Unknown vapor pressures at pores 1, 2, 3 have to found by solving the above system of equations (Eq. (4)). Then, the phase change rate at the partially filled throat (Ṁev,i4 ) can be found by the sum of the vapor flow rates in the gas throats (Ṁv,i1 , Ṁv,i2 , Ṁv,i3 ). In general, for partially saturated pores, the vapor balance can be written as 4−lm  j=1 Ṁv,ij = lm  Ṁev,ij , (5) j=1 with evaporation or condensation rates Ṁev,ij in the pore neighboring meniscus throats (lm). (For several menisci, the phase change rate is distributed according to cross-sectional areas). 2.2.2. Cluster labeling Before we proceed into discrete liquid transport modeling (Section 2.2.3), it is worth to recall the definition of a cluster and explain how to set up data structures for the formed clusters, needed for the discrete description of non-viscous capillary pumping. Cluster by definition is a region of single phase (either gas or liquid) in porous media. Here, we deal only with liquid phase clusters; for the assumed network geometry, a cluster is a group of liquid throats connected by liquid pores. For better understanding see the visualization of clusters 1 and 2 in Fig. 1. Every cluster contains (saturated or unsaturated) meniscus throats, and may or may not contain saturated throats without a meniscus. In every cluster, meniscus throats, where mass transfer in vapor and liquid phase is coupled, must be maintained at capillary equilibrium. So it is necessary to know the exact geometric structure of a cluster to solve the discrete invasion process (Section 2.2.3) due to phase change and to monitor the disintegration process due to invasion. The process of obtaining the structural data of every cluster combined with the method of effective storing of these data for computation of liquid transport is called cluster labeling. To this purpose, the Hoshen–Kopelman algorithm (well documented in literature, Al-Futaisi and Patzek, 2003; Metzger et al., 2006a, b) can be used for cluster labeling. 2.2.3. Liquid transport In the present model, viscous forces are neglected in both liquid and gas phases since throats are of relatively large size; additionally, throats are cylindrical, so that film flow is negligible. Literature work with polygonal throat cross sections shows that film flow in the corners may be important: Yiotis et al. (2004) analyzed different film regimes and also predicted that thermal gradients should usually have little influence on the extension of the film region, but that for large temperature gradients, films ought to be shortened; Prat (2007) mainly investigated the role of pore shape and contact angle. The assumption of negligible viscous forces implies that, in each cluster, liquid can be pumped from the throat with the highest liquid pressure to all other meniscus throats, resulting in one throat with a moving meniscus (MM) and stationary menisci (SM) in all remaining throats. In each cluster, the MM throat will empty first since this throat alone contributes to the total evaporation of the cluster through capillary pumping. The criterion for selecting MM throats is given by capillary pressures Pc,ij , i.e. in each cluster, the throat with highest liquid pressure (6) Pl,ij = P − Pc,ij is selected as MM throat. Capillary pressures can be given by the ratio of temperature dependent surface tension and throat radius: Pc,ij = 2 (Tij ) rij . (7) Using Eq. (7), the MM throat is chosen as the throat with the lowest capillary pressure. Fig. 1 shows an example of two clusters (labeled 1 and 2) in each of which the MM throat is partially filled. Since all mass loss due to evaporation is attributed to MM throats, the time steps to empty these throats can be obtained in each cluster (nc) from the saturation of the MM throat SMM , liquid density l , total throat volume VMM,nc and total evaporation rate of the cluster as S MM tm,nc =  mt j=1 l VMM,nc Ṁev,ij,nc . (8) 2.2.4. Condensation Consideration of heat transfer in a pore network drying model poses difficulties because condensation effects have to be included. In a non-isothermal system evolved temperature gradients over the network instigate equilibrium vapor pressure gradients at phase boundaries. The above phenomenon causes the liquid to evaporate from hotter spots in the network and condense at colder spots. At the same time, capillary effects let the continuous liquid phase split up into disconnected clusters. For a cluster or a single throat, the net phase change rate is calculated as the sum of evaporation/condensation rates on the phase boundaries. Fig. 2 represents the definitions of an evaporating and a condensing cluster. In the  sample cluster (Fig. 2a), if the net phase change rate ( Ṁev,ij ) is positive (=net evaporation), normal invasion rules apply in updating the phase distributions. In this case, capillary pumping nullifies the effect of local condensation. If the net phase change rate is negative, net condensation takes place, and the liquid cluster (Fig. 2b) should grow. For the growth of a condensing cluster, the rules for liquid transport as presented above do not apply. To maintain capillary equilibrium it is reasonable to fill the lowest capillary pressure throat (MM) until it is completely full because, in 5221 V.K. Surasani et al. / Chemical Engineering Science 63 (2008) 5218 -- 5228 filled throats is equally distributed to the neighboring CVs. The dynamic energy balance equation over the CV is given by 4 lm j=1 j=1   dTi Qij − hv,i (Ti ) Ṁev,ij , Vi ( cp )i =− dt > (10) where hv,i (Ti ) is phase change enthalpy (as positive value). The mass flow rate Ṁev,ij is positive for evaporation and negative for condensation. The meniscus position is used to allocate the heat sink/source to either CV. The heat exchange rates Q̇ij between CVs i and j are given by Fourier's law of heat conduction as < Q̇ij = Acv,ij ij Ti − Tj L . (11) The exchange area for heat conduction is Acv,ij = L2 because both void space and solid can contribute. For simplicity, solid and liquid phase resistances are considered as parallel (Strumillo and Kudra, 1986), and contribution of gas is neglected. The effective thermal conductivities ij and CV heat capacities Vi ( cp )i are functions of throat saturation (Surasani et al., 2008a). Two kinds of heat transfer boundary conditions are applied for the convective drying mode; first is heat transfer through the boundary layer as < Q̇ij = Acv,ij , (Ti − T∞ ), Fig. 2. The concept of: (a) evaporating and (b) condensing clusters; (c) condensing clusters (saturated). the capillary dominated regime, capillary pressure at every meniscus in a cluster must be constant: in condensing meniscus throats, the flattening of the meniscus by vapor deposition is compensated by capillary flow to the moving meniscus which has a prescribed curvature. Eventually, when the MM throat is full, all menisci in the cluster will become flat (i.e. having infinite radius of curvature); this situation of a completely saturated condensing cluster is represented in Fig. 2c. The time step to fill the MM throat in each condensing cluster (nc) can be obtained from the ratio of mass of liquid required to fill the MM throat to the total condensation rate of the cluster. tm,nc = − (1 − SMM ) l VMM,nc . mt Ṁ j=1 ev,ij,nc (9) If the MM throat of a condensing cluster is full but it is next to an evaporating cluster, the condensation rate is subtracted from the evaporation rate of this neighboring cluster. Further condensation is not taken into account in the present algorithm, i.e. empty throats are not filled with condensed liquid (imbibition). 2.2.5. Transient heat transfer In many processes, during transport in porous materials, solid, liquid, and gas phases coexist in a chosen scale of CV. In each CV, temperatures of gas (Tg,i ), liquid (Tl,i ), and solid (Ts,i ) may be different. But in most of the drying processes of capillary porous media are characterized by relatively low convective transport in liquid and gas phases so that only small temperature differences exist between the phases in a CV (Whitaker, 1977); therefore, we assume local thermal equilibrium (Ts,i = Tl,i = Tg,i = Ti ). In the same spirit, heat conduction between the CVs (mainly by the solid phase) is assumed dominant over convective heat transfer. In other processes, convection plays a more important role and has been modeled in literature (Satik and Yortsos, 1996; Lu and Yortsos, 2005). Fig. 1 shows how each CV represents one pore and contains four half throats, i.e. each throat participates in two CVs. Liquid in partially (12) and second, heat transfer through insulated sides of the pore network is zero (Q̇ij = 0). In the case of contact heating, a heating plate (at temperature Tcon ) is applied to the bottom edge of the network with perfect contact (numerical implementation is done by taking Eq. (11) and replacing Tj by Tcon ). The explicit scheme for solving Eq. (10) can only be stable if the thermal time step fulfills the condition ( c )i L2 . p tt <  j ij (13) 2.3. Coupling of heat and mass transfer Heat and mass transfer coupling plays an important role in the drying process, particularly at the phase boundaries (menisci) in the network. At an evaporating meniscus, the vapor takes up the enthalpy of evaporation and liquid temperature is reduced (heat sink); likewise there are heat sources at places of condensation. This coupling of mass transfer with heat transfer is already implemented by including the phase change rates (Ṁev ) in the energy balance equation (Eq. (10)). Further coupling is introduced by the temperature dependency ∗ (T ) and local surface tension of local equilibrium vapor pressure Pv,i i (Tij ). In a non-uniform temperature field, equilibrium vapor pressure varies spatially. This has two major effects on vapor diffusion: first, it causes higher drying rates at hot spots and, second, it causes back diffusion in the network with condensation at colder spots. Equilibrium vapor pressures have to be updated after each time step by use of an Antoine equation. At the same time, surface tension varies inversely with temperature. By Eq. (7), a surface tension gradient may have an effect on IP patterns depending on magnitude and orientation of the temperature gradient (Plourde and Prat, 2003). 2.4. Drying algorithm A flow sheet representation of the drying algorithm is given in Fig. 3. First, the data structures describing the solid network and boundary layer must be formed from network size; additionally, initial and boundary conditions are set (step 1) and the liquid cluster needs to be identified using the cluster labeling algorithm (step 2). The linear system for vapor transport (Eq. (4)) 5222 V.K. Surasani et al. / Chemical Engineering Science 63 (2008) 5218 -- 5228 is solved and evaporation rates of each cluster are computed as the sum of the evaporation rates on their respective menisci using Eq. (5) (steps 3 and 4). Then, from the capillary pressures Pc,ij in meniscus throats, the MM throat is identified for each cluster and the time step tmin is chosen as a minimum from the values of Eqs. (8), (9), and (13) (step 5). The mass loss due to evaporation is computed for the time step tmin and phase distributions (pore and throat saturations) are updated (step 6). Next, the new temperature field is computed by use of Eq. (10) for same Δ Fig. 3. Algorithm of pore-scale drying model. ∞ ∞ time step; and, finally, temperatures as well as temperature dependent variables are updated (step 7). The procedure will be repeated until all liquid is removed from the network. 3. Drying simulations 3.1. Heating modes and parameters We are interested in demonstrating the influence of heat transfer by changing the thermal boundary conditions. To this purpose two different drying methods are used; the first one is classical convective heating, and the second is less resistive contact heating. The pictorial representations of the two modes are shown in Fig. 4. In both cases, the top surface of the solid network is open to vapor diffusion through the boundary layer and all the remaining sides are impermeable to mass transfer. A boundary layer forms parallel to the open surface of the solid network due to air convection. It is supplied at velocity such that the boundary layer contains seven nodes ( = 10.4 mm/s). Dry air has zero humidity (Pv,∞ = 0). In convective heating mode, its temperature is 80 ◦ C and all three faces of the network are insulated (Fig. 4a). In contact heating mode, drying air is at 20 ◦ C and the heat transfer is due to perfect thermal contact of the bottom edge of the solid network (Fig. 4b) with a heating plate at 80 ◦ C; again the left and right sides are insulated. One of the major features of the pore network model is to account for the geometry of void space, and its major outcome is the formation of liquid clusters as a consequence of spatial variation in pore size distribution. Thus, radii of the throats are chosen randomly according to a normal number density function and the length of every throat is 350 m. The results presented in this section have been obtained for networks of size 51 × 51 with a bimodal throat size distribution. Every fifth throat orthogonal to the open surface of the network is chosen from a radius distribution of 100 ± 5 m (macrothroats) and all remaining throats from a distribution of 40 ± 2 m (micro-throats). Applied throat radius distributions for bimodal pore network are large enough to justify the assumption of neglecting viscous effects, which are dominant in network of small radius distributions (Metzger et al., 2007a). ∞ ∞ Fig. 4. Heating modes used in simulations. (a) Convective heating; (b) contact heating. V.K. Surasani et al. / Chemical Engineering Science 63 (2008) 5218 -- 5228 7 non−isothermal isothermal 6 Mv [mg/h] 5 ⋅ 4 3 2 1 0 0 0.2 0.4 0.6 0.8 1 S [−] Fig. 5. Isothermal mode and convective heating mode drying rate curves for bi-modal network. Many technical porous media exhibit bimodal pore size distribution and both micro- and macro-throats need to be continuous in reality. However, in the two-dimensional pore network such a topology is strictly not possible, but real drying behavior can be simulated by the above-described network (Metzger et al., 2007b). Three-dimensional simulations with more realistic spatial distributions of macro-throats will be presented in Surasani et al. (2008b). Initially, the network is saturated with water and has a uniform temperature of 20 ◦ C. The heat transfer parameters of the solid are chosen as for glass: ( cp )s = 1.7 × 106 J/(m3 K) and s = 1 W/(m K). Simulation results of both convective heating and contact heating modes are presented in the following. 3.2. Results for convective heating mode Fig. 5 shows the drying kinetics as evaporation rate Ṁv versus network saturation S. The isothermal curve is plotted for uniform temperature of 20 ◦ C to show the effect of heat transfer. Both the isothermal and the non-isothermal curve start at the same rate. The initial warming-up to constant drying conditions can be clearly seen in the non-isothermal process. The constant rate period extends for both curves to a saturation of approximately 0.65. The order of emptying of throats is determined by capillary pressure at the menisci (Eqs. (6) and (7)). In the isothermal model, invasion of air depends on throat radius alone, since temperature is constant throughout the network. In the convective heating mode, the order of empting may also be influenced by surface tension (Eq. (7)). At the start of the drying process, macro-throats empty since the effect of temperature on capillary pressure via surface tension is negligible (see below discussion). Macro-throats at high liquid pressures (low capillary pressures) can pump liquid to the micro-throats at the network surface which stays almost completely wet down to relatively low overall saturations. As a consequence, a first drying period can be observed (the few dry surface throats are additionally masked by lateral vapor transfer in the boundary layer, see Metzger et al., 2005). The evaporation rate drops drastically when micro-throats on the surface start to dry out. The evolution of IP patterns and temperature fields, corresponding to the isothermal mode and convective heating mode drying simulations in Fig. 5, are shown in Fig. 6 at different overall saturations of the network. Fig. 7 illustrates the variation of temperature at the 5223 top of the network (Tmax ), temperature at the bottom of the network (Tmin ), and network saturation S with time. The situations in convective heating at saturation S equal to 0.8 and 0.65 are located within the constant-rate period (see Figs. 5 and 7), in which the network has almost uniform temperature and only macro-throats will empty. Near wet bulb temperatures are attained at the end of the first drying period. These attained temperatures are higher than wet bulb temperature (25 ◦ C) due to the warming up of dry spots formed after emptying of macro-surface throats. Surface throats empty in the approximate saturation range of 0.65–0.48. Short quasi-constant rate periods can be observed in this falling rate period in isothermal drying when—due to random radius distribution—throat saturation near the surface stays unchanged while internal throats are invaded. In the case of convective heating, the surface heats up and evaporation rate increases, until another near surface throat empties. Then, drying rate drops to a lower value in both model versions. These quasi-constant rate periods can be seen in the drying curve (Fig. 5) between saturations 0.65 and 0.48. In the falling rate period, larger temperature gradients develop and network temperature starts to rise drastically. At the same time, drying rate drops, as can be clearly seen from the sudden change of slope on the saturation curve in Fig. 7. After all the macro- and surface throats have been emptied, the drying front recedes (see transition from S = 0.48 to 0.10) towards the bottom of the network. This increases the resistance to mass transfer; however, the resulting rise in temperature causes a higher driving force for vapor diffusion. Due to the latter effect, evaporation rates are much larger than those of the isothermal model in this region (Fig. 5). The evolved temperature profiles shown in Fig. 7 reach a uniform temperature (of drying agent) in post drying period. Fig. 6 also presents the comparison of simulated phase evolutions obtained from isothermal and convective heating mode. Very similar IP patterns are observed in the two cases, the reason for which shall now be discussed. The influence of temperature on surface tension of water can be represented by the equation (N/m) = −1.704 × 10−4 T (◦ C) + 7.636 × 10−2 , (14) so that the relative change in surface tension per 1 K temperature variation is approximately 0.2%. The maximal temperature difference over the network is about 6 K (see temperature fields in Fig. 6). This corresponds to a surface tension variation among the meniscus throats of maximally 1.2%. The relative standard deviation of throat radius is 5% for both modes of the distribution (100±5 and 40±2 m). Only macro-throats will empty at the start of the drying process. As long as a cluster contains macro-throats, they are so much larger in radius that a temperature gradient cannot lead to preferential emptying of micro-throats (S = 0.8). When a complete macro-channel is invaded, this leads to splitting up of the liquid phase into clusters (2D artifact), and micro-throats start emptying. Since the temperature gradient is very small in this period, the invasion order is essentially the same as for the isothermal case (S  0.65). When surface micro-throats start emptying, temperatures start to rise. And if the relative variation of competing meniscus radii is lower than that of surface tension, we expect preferential invasion of throats at higher local temperatures (see Fig. 6 at S = 0.48). However, in convective heating mode, evolved temperature gradient is low, and surface tension will not play a big role for the given network disorder (i.e. standard deviation of radii). This implies IP patterns are similar to that of isothermal simulations. However, if the disorder of the network is very small then the effect of surface tension gradients can become more important (see also Plourde and Prat, 2003). 3.3. Results for contact heating mode To depict the influence of heat transfer on macroscopic behavior, simulated drying rates for convective and contact heating modes are 5224 V.K. Surasani et al. / Chemical Engineering Science 63 (2008) 5218 -- 5228 Fig. 6. Evolution of IP patterns and temperature fields for convective heating mode and comparison with simulated IP patterns of isothermal mode (network saturations are indicated; temperatures values are in ◦ C). compared with isothermal drying rates (Metzger et al., 2007b) at room temperature in Fig. 8. Temperature dependency of equilibrium vapor pressure directly affects evaporation rates leading to overall higher drying rates for both non-isothermal simulations (Fig. 8). All three simulations show a first drying period because macrothroats empty first keeping the surface wet down to low network V.K. Surasani et al. / Chemical Engineering Science 63 (2008) 5218 -- 5228 80 1 0.8 S 60 0.6 S [−] T [K] Tmin Tmax 0.4 40 0.2 0 20 0 2 4 6 8 10 t [h] Fig. 7. Temporal evolution of temperatures at (Tmax ) top and bottom (Tmin ) of the network and its overall saturation S. 10 isothermal contact heating convective heating ⋅ Mv [mg/h] 8 6 4 2 5225 takes place at the bottom of the network (see Fig. 9 at S = 0.60) due to the high gradient in temperature. A 40 K temperature difference at S = 0.60 corresponds to a variation of approximately 9% in capillary pressures through surface tension gradients in comparison to the relative standard deviation in one mode of the radius distribution of 5%. Therefore, a lower surface tension at the bottom will allow to empty throats which would not have been emptied for their relatively small size. Considering the network saturations 0.60–0.40, capillary fingering occurs in the entire network, i.e., forming of numerous clusters. This complex behavior is determined by the random dominance of radius distribution over surface tensions or vice versa. The small disconnected clusters evaporate not only at the top of the network, but also clusters at the network bottom can be eroded because of the significantly higher equilibrium vapor pressure in that region. As a result, two drying fronts will form which progress to the middle of the network (see Fig. 9 from S=0.40 to 0.10). With increasing mass transfer resistance for the upper evaporation front, drying rates decrease and evaporative cooling becomes less important; consequently, temperature gradients become smaller and smaller. One of the major consequences of heat and mass transfer coupling of equilibrium vapor pressures (in pore network modeling) is diffusion of vapor from hot to cold neighboring liquid–gas interfaces and subsequent condensation. This process replaces non-wetting fluid by wetting fluid and does not follow IP rules as proposed in liquid transport. As explained in Section 2.2.4 only part of the condensing clusters (Fig. 2b) are chosen for accommodating condensed vapor. Fig. 10 shows the total condensation rate of the network ṀC for both heating modes as a function of network saturation S. The actual rate as computed from vapor balances—and used for heat balances—is plotted along with the neglected rate which could not be accommodated (because imbibition is not modeled). During the initial phase of drying, while only macro-throats are emptying, capillary effects in a spanning liquid cluster nullify condensation, so that no error is made. When the liquid phase is discontinuous, condensation cannot fully be accounted for, and the error is higher in contact mode than in convective mode. In the case of convective heating, with stabilized drying front and low temperature the neglected condensation rates appear to be negligible (see also Surasani et al., 2008a). In contact heating mode, with high temperature gradient, two drying fronts and many liquid clusters, the error in condensation is relatively high. 0 0 0.2 0.4 0.6 0.8 1 S [−] Fig. 8. Drying curves for isothermal case and both considered heating modes. saturations. Additionally, a constant surface temperature is approached which corresponds to wet-bulb temperature for pure convective heating; for contact heating, it is higher and results from steady-state heat conduction from network bottom to top of boundary layer and from the cooling by surface evaporation. The order of throat emptying depends on the temperature field by capillary pressure (Eq. (7)); this determines the position of the evaporation front which in turn controls drying rates. The evolutions of IP patterns and temperature fields for contact heating mode are shown in Fig. 9 and, for comparison, the results for the convective heating mode are also presented. In contact drying, again only macrothroats empty at the start (see Fig. 9 at S = 0.80) because of the high difference in radii for the two modes in throat size distribution. In this period, maximum temperature gradient develops in contact heating mode. Therefore, the invasion order of macro-throats differs from that of convective heating mode (Fig. 9). After breakthrough (BT), i.e., air reaching the network bottom, and when clusters without macro-throats are created, invasion also 3.4. Evolution of phase regions Fig. 11 represents the evolution of phase regions of the bimodal network for both cases of convective and contact heating. The phase evolutions are plotted in terms of dimensionless drying front position in the network versus network saturation S. The regions of evolved phases, i.e., gas, liquid, and liquid–gas interface, are named in the figure. The upper boundary in the liquid–gas interface indicates the dimensionless height above which only gas exists. Above this boundary, only vapor diffusion takes place. One can see the disconnection point (DP) in both figures. It indicates the point where all the surface throats dry out and the overall drying rates fall drastically. In convective heating, the lower boundary of liquid–gas interface indicates the dimensionless height below which only liquid exists. Only liquid flow due to capillary pumping takes place in this region. The point, at which air invades the bottom of the network, is called breakthrough point (BT). In contact heating, the liquid–gas interface lower boundary progresses in a similar way as in convective heating until BT. Later on, however, it moves upwards to give way to a gas region by empting of bottom throats. Unlike monomodal pore network drying simulations (Surasani et al., 2007) in which BT occurs after main cluster disconnection (MCD), bimodal simulations exhibit BT before DP (both MCD and DP represent the point where the surface dries out). This is due to the emptying of 5226 V.K. Surasani et al. / Chemical Engineering Science 63 (2008) 5218 -- 5228 Fig. 9. Influence of heating modes on IP patterns and temperature fields (network saturations are indicated; temperatures values are in ◦ C). macro-throats at the beginning of drying process in the bi-modal network, whereas in the case of mono-modal surface throats will empty first. 4. Conclusions A recently proposed non-isothermal pore network drying model (Surasani et al., 2008a), has been modified to compare convective and contact heating modes for the drying of a two-dimensional square pore network with a bimodal radius distribution. As already discussed in literature (Plourde and Prat, 2003; Huinink et al., 2002; Surasani et al., 2008a), the coupling between heat and mass transfer can have influences on both drying rates and evolving phase patterns in the porous medium. The first major expected effect is the enhancement of evaporation rate because the warming up of the porous body during drying significantly increases equilibrium vapor pressures at the evaporation front. Second, the temperature dependency of surface tension is expected to alter the order of pore emptying, leading to different phase patterns than under isothermal conditions. A bimodal network has been chosen because many real porous media exhibit micro- and macro-pores. A major feature of such a pore network is that capillary pumping empties macro-pores first while the micro-pores stay saturated with liquid. This ensures a period of almost constant drying rate because the network surface remains nearly saturated. In the two-dimensional representation with several macro-channels (see Fig. 8), micro-porous network-spanning liquid clusters are formed in this period (Metzger et al., 2007b). In both non-isothermal drying modes (convective and contact), the macrochannels still empty preferentially because surface tension gradients have negligible effect on capillary pressure in comparison to the large difference in pore radii between micro- and macro-pores. In the convective heating mode, evolving temperature gradients in the liquid–gas region are small throughout the drying process 5227 V.K. Surasani et al. / Chemical Engineering Science 63 (2008) 5218 -- 5228 1 0.8 ζ [−] 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1 0.6 0.8 1 S [−] 1 0.8 ζ [−] 0.6 0.4 0.2 0 0 0.2 0.4 S [−] Fig. 11. Evolution of phase regions in both (a) convective and (b) contact heating modes. Fig. 10. Actual and neglected condensation rates for (a) convective and (b) contact heating. so that only minor surface tension gradients develop. As a consequence, invasion order is almost not affected in both periods, i.e., during emptying of macro-throats and during drying out of microporous clusters. However, the temperature level gradually rises as the convective heating overrules evaporative cooling; this is because vapor diffusion resistance increases when the evaporation front recedes into the porous medium. With rising temperature, equilibrium vapor pressure at the phase boundary also increases leading to significantly higher drying rates than for the isothermal case. A similar, but even more pronounced effect on overall drying rates is observed for the contact heating mode, since higher temperatures are reached for this less resistive heat transfer. At the same time, temperature gradients are much higher, so that surface tension variation plays a vital role in this case. Already when macro-throats are emptied, the order of invasion is modified. And after microporous clusters have been formed, the direction of temperature gradient leads to preferential emptying of throats at the bottom of the network. We observe what in literature is called IP in a destabilizing gradient (IPSG) (Plourde and Prat, 2003). A second effect occurs when the micro-porous regions have split up into small clusters. Then, the high variation in local equilibrium vapor pressure at the menisci leads to fast evaporation of clusters near the network bottom. The vapor is evacuated through the partially saturated center region of the network by the empty macro-channels. This second mechanism continues the destabilized emptying of throats. In the same time, near surface clusters also evaporate, due to the easy vapor diffusion to the dry convective air. This means that two drying fronts move towards the core of the network until they meet and the network is completely dry. Very similar behavior has been reported by Huinink et al. (2002). As described above, the present model does not fully account for condensation, and the error has been quantified (Fig. 10). The main target of future work is to couple the invasion rules (for gas) with new imbibition rules (for liquid) so that fully saturated condensing clusters may grow into adjacent pores with free void volume. To this purpose, rules for the choice of the wetted pore throat have to be set up and details of reorganizing capillary equilibrium and cluster coalescence will be developed. 5228 V.K. Surasani et al. / Chemical Engineering Science 63 (2008) 5218 -- 5228 Notation A cp g H L M̃ Ṁ N P Q̇ r R̃ S t T u V W h area of cross section, m2 specific heat capacity, J/kg/K mass transfer conductance, kg/s enthalpy content, J length of throat, m molar mass, kg/kmol mass flow rate, kg/s number of pores/clusters, dimensionless pressure, Pa heat flow rate, W throat radius, m ideal gas constant, J/kmol/K saturation, dimensionless time, s temperature, K velocity, m/s volume, m3 network thickness, m specific phase change enthalpy, J/kg Greek letters       heat transfer coefficient, W/m2 /K mass transfer coefficient, m/s surface tension, N/m vapor diffusivity, m2 /s boundary layer thickness, m dimensionless height, dimensionless thermal conductivity, W/m/K kinematic viscosity, m2 /s density, kg/m3 Subscripts/superscripts BL c C con ev g i or j ij l lm nc s T v ∗ boundary layer capillary condensation contact evaporation gas pore or node throat connecting pores i and j liquid number of menisci cluster label solid total network vapor equilibrium Acknowledgments This work was financed by the German Research Foundation (DFG) in the frame of Graduate School 828 “Micro-MacroInteractions in Structured Media and Particle Systems” and received additional support from Max-Buchner-Forschungsstiftung. The first author would like to express his special thanks to Dr. M. Prat for the fruitful discussion. References Al-Futaisi, A., Patzek, T.W., 2003. Extension of Hoshen–Kopelman algorithm to nonlattice environments. Physica A 321, 665–678. Berkowitz, B., Ewing, R.P., 1998. Percolation theory and network modeling applications in soil physics. Surveys in Geophysics 19, 23–72. Guillot, G., Trokiner, A., Darrasse, L., Saint-Jalmes, H., 1989. Drying of a porous rock monitored by NMR imaging. Journal of Physics D—Applied Physics 22, 1646–1649. Huinink, H., Pel, L., Michels, M.A.J., Prat, M., 2002. Drying processes in the presence of temperature gradients—pore-scale modeling. European Physical Journal E 9, 487–498. Irawan, A., Metzger, T., Tsotsas, E., 2005. Pore network modeling of drying: combination with boundary layer model to capture the first drying period. In: Proceedings of Seventh World Congress of Chemical of Chemical Engineering, Glasgow, Scotland. Laurindo, J.B., Prat, M., 1998. Modeling of drying in capillary-porous media: a discrete approach. Drying Technology 16 (9 and 10), 1769–1787. Lu, C., Yortsos, Y.C., 2005. Dynamics of forward filtration combustion at the porenetwork level. A.I.Ch.E. Journal 51 (4), 1279–1296. Masmoudi, W., Prat, M., 1991. Heat and mass transfer between a porous media and a parallel external flow, application to drying of porous materials. International Journal of Heat and Mass Transfer 34 (8), 1975–1981. Metzger, T., Irawan, A., Tsotsas, E., 2005. Discrete modeling of drying kinetics of porous media. In: Proceedings of Third Nordic Drying Conference, Karlstad, Sweden, 2005. Metzger, T., Irawan, A., Tsotsas, E., 2006a. Remarks on the paper “Extension of Hoshen–Kopelman algorithm to non-lattice environment” by Al-Futaisi, A., Patzek, T.W., 2003. Physica A 321, 665–678. Metzger, T., Irawan, A., Tsotsas, E., 2006b. Remarks on the paper “Extension of Hoshen–Kopelman algorithm to non-lattice environment” by Al-Futaisi, A., Patzek, T.W., 2003. Physica A 363, 558–560. Metzger, T., Irawan, A., Tsotsas, E., 2007a. Isothermal drying of pore networks influence of friction for different pore structures. Drying Technology 25, 49–57. Metzger, T., Irawan, A., Tsotsas, E., 2007b. Influence of pore structure on drying kinetics: a pore network study. A.I.Ch.E. Journal 53, 3029–3041. Metzger, T., Tsotsas, E., Prat, T., 2007c. Chapter 2: Pore-network models: a powerful tool to study drying at the pore level and understand the influence of structure on drying kinetics. In: Tsotsas, E., Mujumdar, M. (Eds.), Modern Drying Technology Computational Tools at Different Scales, vol. 1. Wiley, New York, pp. 57–102. Nowicki, S.C., Davis, H.T., Scriven, L.E., 1992. Microscopic determination of transport parameters in drying porous media. Drying Technology 10 (4), 925–946. Pel, L., Ketelaars, A.A.J., Adan, O.C.G., van Well, A.A., 1993. Determination of moisture diffusivity in porous media using scanning neutron radiography. International Journal of Heat and Mass Transfer 36, 1261–1267. Pel, L., Brockem, H., Kopinga, K., 1996. Determination of moisture diffusivity in porous media using moisture concentration profiles. International Journal of Heat and Mass Transfer 39, 1273–1280. Plourde, F., Prat, M., 2003. Pore network simulations of drying of capillary porous media. Influence of thermal gradients. International Journal of Heat and Mass Transfer 46, 1293–1307. Prat, M., 1993. Percolation model of drying under isothermal conditions in porous media. International Journal of Multiphase Flow 19 (4), 641–704. Prat, M., 2002. Recent advances in pore-scale models for drying of porous media. Chemical Engineering Journal 86, 153–164. Prat, M., 2007. On the influence of pore shape, contact angle and film flows on drying of capillary porous media. International Journal of Heat and Mass Transfer 50, 1455–1468. Quintard, M., Whitaker, S., 1993. Transport in ordered and disordered porous media: volume-averaged equations, closure problems and comparison with experiments. Chemical Engineering Science 41 (14), 2537–2564. Sanchez-Palencia, E., 1980. Non-homogeneous media and vibration theory. Lecture Notes in Physics, vol. 127. Springer, Berlin. Satik, C., Yortsos, Y.C., 1996. A pore-network study of bubble growth in porous media driven by heat transfer. Journal of Heat Transfer 118, 455–462. Strumillo, C., Kudra, T., 1986. Drying principles applications and design. In: Hughes, R. (Ed.), Topics in Chemical Engineering, vol. 3. Gordon and Beach Science Publishers. Surasani, V.K., Metzger, T., Tsotsas, E., 2007. A non-isothermal pore network drying modal: influence of gravity. In: Proceedings of European Congress of Chemical Engineering (ECCE-6), Copenhagen, 2007. Surasani, V.K., Metzger, T., Tsotsas, E., 2008a. Consideration of heat transfer in pore network modelling of convective drying. International Journal of Heat and Mass Transfer 51, 2506–2518. Surasani, V.K., Metzger, T., Tsotsas, E., 2008b. Drying simulations of various 3D pore structures by a non-isothermal pore network model. Submitted to 16th International Drying Symposium, Hyderabad, India. Yiotis, A.G., Stubos, A.K., Boudouvis, A.G., Yortsos, Y.C., 2001. Pore network model of drying of single-component liquids in porous media. Advances in Water Resources 24, 439–460. Yiotis, A.G., Boudouvis, A.G., Stubos, A.K., Tsimpanogiannis, I.N., Yortsos, Y.C., 2004. Effect of liquid films on the drying of porous media. A.I.Ch.E. Journal 50 (11), 2721–2737. Yiotis, A.G., Tsimpanogiannis, I.N., Stubos, A.K., Yortsos, Y.C., 2006. Pore-network study of characteristic periods in the drying of porous media. Journal of Colloid and Interface Science 297, 738–748. Whitaker, S., 1977. Simultaneous heat, mass, and momentum transfer in porous media: a theory of drying. Advances in Heat Transfer 13, 119–203.