Multiscale Modeling of Non-Isothermal Fluid Transport Involved in Drying Process of Porous Media

Abstract

To preserve the product quality as well as to reduce the logistics and storage cost, drying process is widely applied in the processing of porous material. In consideration of transport phenomena that involve a porous medium during drying, the complex morphology of the medium, and its influences on the distribution, flow, displacement of multiphase fluids are encountered. In this chapter, the recent advanced mass and energy transport models of drying processes are summarized. These models which were developed based on both pore- and continuum-scales, may provide a better fundamental understanding of non-isothermal liquid–vapor transport at both the continuum scale and the pore scale, and to pave the way for designing, operating, and optimizing drying and relevant industrial processes.

1. Introduction

In both natural systems (i.e. clays, aquifers, oil and gas reservoirs, plants, and biological tissues) and industrial systems (i.e. fuel cells, concrete, textiles, polymer composites, capillary heat pipe, and paper, etc.), porous media are often encountered. From a morphological point of view, porous media is composed of a persistent solid matrix and an interconnected void space that can be occupied by fluid phases. In several porous systems, the void space is initially filled by liquid water. To maintain the product quality, prolong the storage time as well as reduce the logistics cost, the liquid water is needed to be removed by using various drying techniques. During the drying process, the complex morphology of the media and its influences on the capillary flow, liquid – vapor phase change under thermal effects are faced. These processes are complex by themselves [1, 2]. Additionally, it should be noted that the industrial dryers consume approximately 12% of the total energy used in manufacturing processes [3]. For several industrial sectors (i.e. tissue, food, and agriculture) the energy consumption ratio of drying may reach 33%. As a complicated, important, and high energy consumption process, the dryer designing, and operation should not be done by trial and error. A better fundamental understanding of this drying process at the pore- and continuous scales (i.e. understanding the influence of porous structure on the drying behavior, drying time determination, and energy consumption prediction) may help to design better products and processes.

The drying process of porous media may be considered at different scales. At different scales, the approaches used to describe the transport processes may be quite different [4]. For instance, in a drying plant, thousands or millions of drying products can be dried simultaneously and the state of the drying agent evolves over time and space. Thus, the intra-sample distribution of process parameters should be lumped to avoid the extremely expensive computational cost. As a result, the lumped model is the suitable model that should be implemented in drying plant simulation. Coming to the sample drying process, the internal transport process is taken into account. The heat, mass, and momentum conservation equations are developed based on the first physics principles. In a control volume, the process properties are assumed to be homogenous and isotropic and the impact of micro-heterogeneity on the transport process is omitted. The macro-scale models often result in a system of partial differential equations. Both globe and local drying behavior of the sample can be well predicted by using the continuous models if the size of the domain is large enough. For a thin layer of the porous medium, where micro-heterogeneity plays a role, the continuous models are not valid. Due to the interaction between the solid skeleton and the fluid phases, the transport equations should be derived for individual pores. Since the transport phenomena are directly investigated at the pore-scale, the impact of changes at the pore or microstructure such as wettability, pore shape, pore size distribution, and pore structure on drying behavior has been effectively interpreted by pore-scale models [5, 6]. One should notes that these models are not standing individually. The results of pore scale model can be upscaled to generate the effective parameters of macro scale models which can be reduced to the lumped models. An example multiscale modeling approach for wood particle drying process is presented in Figure 1.

Figure 1.

Multiscale modeling approach for porous particle drying process.

In this chapter, a review of the recent advanced model in a selection of drying phenomena involved in porous media is illustrated. We discussed the state-of-the-art numerical methods as complementary ways to get more insight. The future challenges and the hint at solutions to accommodate for them are also given.

2. Pore-scale models

In 1950s, a new modeling approach to describe the transport process in the porous medium was proposed by Fatt. In his model, the void space of the porous media was modeled by a network of cylindrical tubes. In individual tubes, the mass transfer equation has been written and the discrete mass conservation equations are written for each pore node. This model has been named as pore network model. Afterward, several pore-scale models, such as the direct numerical approach, the volume of fluid approach, and method, have been presented in literature based on the idea of pore network model. Since the transport phenomena are considered directly at the pore-scale, the better fundamental understanding on the interaction between solid structure and fluid flow can be provided. As discussed in Section 1, the mentioned interaction can be upscaled to determine the macroscopic transport coefficients of the continuous models such as the relative permeabilities of the liquid–gas mixture, the liquid diffusivity.

In this section, we focus on the non – isothermal drying pore network models and their application in estimating the macroscopic transport coefficients.

2.1 Pore network model for superheated steam drying

Based on the invasion percolation algorithm, several non-isothermal pore network models were developed to simulate the hot air drying process [2, 7, 8, 9]. These models accounted capillary, gravity effects, viscous effects, and the transport of vapor by diffusion in the gas phase under the thermal effect. Recently, the pore network model has been applied to describe the non-isothermal mixture vapor - liquid water transport during the superheated steam drying process [10]. The features of this pore network model is the fully treatment of condensed liquid by introducing the newly formulated liquid invasion rules in a two-dimensional domain. In this chapter, this two-dimensional pore network model is extended to three-dimensional model in this work to simulate the transport processes inside a capillary porous medium undergoing superheated steam drying with two different heating modes: (i) the convective heating at the top surface (convective- heating mode) and (ii) the simultaneous convective and conductive heating at the top and bottom surfaces (contact-heating mode) (c.f. Figure 2).

Figure 2.

Two different heating modes used in the pore network simulations.

Due to the absence of diffusion in the gas phase, the water vapor is transported by convective mechanism only and the Hagen-Poiseuille law is applied to calculate the vapor mass flow rate between two adjacent pores [10].

$${\dot{M}}_{v,\mathit{ij}}=\frac{\pi r_{\mathit{ij}}^4}{8{\upsilon }_{v}L}\left(p_i-p_j\right)$$E1

where the vapor mass flow rate in the cylindrical throat connecting pores i and j is denoted by ${\dot{M}}_{v,\mathit{ij}}$, p_i and p_j present the pressure of the vapor phase at pores I and j, respectively. The length and radius of throat are respectively indicated by L and r. The kinematic viscosity of the vapor phase is denoted by ${\upsilon }_v$. Under the quasi-steady assumption, the mass conservation of water vapor in a empty/partially pore is written as

$$\sum _j{\dot{M}}_{v,\mathit{ij}}-\sum _j\pi {r_{t,\mathit{ij}}}^2{\dot{m}}_{\mathrm{evp}/\mathrm{con},\mathit{ji}}-\pi {r_{p,i}}^2{\dot{m}}_{\mathit{evp}/\mathit{con},i}=0$$E2

Due to the non-uniform pore size distribution, the liquid within the network is transported under capillary action. Thus, the large pores/throats are preferentially emptied or filled to maintain the lowest liquid pressure according to Eq. 4. As a result, the liquid phase disintegrates into several liquid clusters. Since the liquid mass balance must be satisfied in the entire network as well as for each liquid cluster, the morphological information of liquid cluster is required. The Hoshen–Kopelman algorithm is applied for labeling the liquid clusters.

$$p_l=p_v-\frac{2\sigma \left(T\right)cos\theta }{r_{p/t}}$$E3

To fully model the phase transition including both evaporation and condensation, in addition the emptying and refilling events, the liquid invasion represented in Figure 3 is introduced in the present model. The time steps for the emptying and refilling events of meniscus pores or throats are computed by

Figure 3.

The 2D sketch of liquid invasion rules: Liquid invasion from single liquid throat to its neighboring pore (a), liquid invasion from fully filled liquid throat in a liquid cluster to its neighboring empty pore (b), liquid redistribution from fully filled pore in a liquid cluster to its neighboring empty throat (c).

$$\Delta t_{\mathit{\text{empty}},\mathit{cn}/\mathit{st},p/t}=\frac{V_{p/t}{\rho }_l\left(T\right)S_{p/t}}{{\sum }_{\mathit{cn}/\mathit{st}}{\dot{M}}_{\mathit{evp}/\mathit{con}}}\text{and}\Delta t_{\text{refill},\mathit{cn}/\mathit{st},p/t}=-\frac{V_{p/t}{\rho }_l\left(T\right)\left(1-S_{p/t}\right)}{{\sum }_{\mathit{cn}/\mathit{st}}{\dot{M}}_{\mathit{evp}/\mathit{con}}}$$E4

The time step for the first two liquid invasion events is calculated by

$$\Delta t_{\mathit{\text{in}},\mathit{cn}/\mathit{st}}=-\frac{V_{p/t,}{\rho }_l\left(T\right)}{{\sum }_{\mathit{cn}/\mathit{st}}{\dot{M}}_{\mathit{evp}/\mathit{con}}}$$E5

The pore/throat saturations after the capillary re-equilibration event are computed by

$$S_{t,\mathit{ih}}=1\text{and}{S}_{p,h}=\frac{V_{p,h}-V_{t,\mathit{ih}}}{V_{p,h}}$$E6

The thermal energy supplied to a control volume (c.f. Figure 4) lead to a change of the enthalpy in the control volume: the control volume temperature increases, and the phase transition occurs. Based on a fully implicit scheme, the energy balance equation written for the control volume around pore i in time step Δt is discretized as

Figure 4

Heat balance at a vapor pore. The arrows indicate the direction of the heat flux towards each face of the control volume.

$$\begin{array}{l}\sum _j-A_{\mathit{cv},\mathit{ij}}{\lambda }_{\mathit{eff},\mathit{ij}}\frac{T_{i,t+\mathrm{\Delta }t}-T_{j,t+\mathrm{\Delta }t}}{L}=\frac{{\rho }_{\mathit{eff},i}{c}_{\mathit{eff},i}{V}_i}{\mathrm{\Delta }t}\left(T_{i,t+\mathrm{\Delta }t}-T_i\right)\\ +\mathrm{\Delta }{H}_{\mathit{evp}}\sum _j\pi {r_{t,\mathit{ij}}}^2{\dot{m}}_{\mathit{evp}/\mathit{con},\mathit{ji}}+\mathrm{\Delta }{H}_{\mathit{evp}}\pi {r_{p,i}}^2{\dot{m}}_{\mathit{evp}/\mathit{con},i}\end{array}$$E7

The convective heat transfer boundary condition presented in Eq. 8 is applied for both heating modes. At the bottom of the network, perfect thermal insulation (Eq. 9) is imposed for the convective heating mode, whereas uniform constant temperature (Eq. 10) is set for the contact-heating mode.

$$\left[{\left.{\lambda }_{\mathit{eff}}\frac{\partial T}{\partial z}\right|}_{y=L_y}+{\left.\mathrm{\Delta }{H}_{\mathit{evp}}{\rho }_l{v}_l\right|}_{y=L_y}\right]•n=\alpha \left(T_{\mathit{\text{steam}}}-{\left.T\right|}_{y=L_y}\right)$$E8

$${\left.\frac{\partial T}{\partial y}\right|}_{y=0}=0$$E9

$${\left.T\right|}_{y=L_y}=T_{\mathit{\text{bottom}}}$$E10

The simulations are carried out for a 10 × 10 × 10 pore network with throat radius distribution 100 ± 10 μm, pore radius distribution 250 ± 10 μm and uniform distance between two adjacent pores L = 1000 μm. Because of the high computational demand required for a full simulation of drying, only one realization of network was considered here. The solid skeleton of the network is made of typical glass with thermal conductivity λ_s = 1 W/mK. The steam velocity at the network surface is 5 m/s. For the contact-heating mode, a uniform temperature boundary condition Tb_ottom = 105°C is imposed at the network bottom (c.f. Eq. 10), whereas the impermeable heat transfer boundary condition is set for the convective heating mode (c.f. Eq. 9). To depict the influence of heating mode on drying characteristics, the evolution of drying rate and of the local network saturation over network saturation obtained for these two heating modes are shown in Figures 5 and 6.

Figure 5.

Impact of heating mode on the local saturation evolution over network saturation.

Figure 6.

Impact of heating mode on the drying kinetic curve.

Since a large amount of thermal energy was supplied to the network with the contact-heating mode, the drying rate remains high and drying time is thus short compared to the network with the convective-heating mode. Also, since the surface temperature of the network with the contact-heating mode reaches boiling temperature fast, the surface condensation period is short and thus the amount of condensed liquid at the network surface is less, compared to the convective-heating mode. Moreover, it can be seen that the semi-constant drying rate period obtained from the network with the contact-heating mode is long compared to the convective heating model. This tendency can be explained by the temperature dependency of the surface tension; a lower surface tension obtained at high temperature leads to a lower capillary pressure and a higher liquid pressure at the liquid–vapor-solid interface. The pores/throats near the bottom of the network with the contact-heating mode are preferentially emptied. As a result, the vapor fingers appear markedly, and the bottom of the network is dried completely when the middle and top zones are still partly wet. The impact of vapor fingering on the phase distribution is shown in Figure 5. More number of liquid clusters is obtained with the contact-heating mode compared to the convective-heating mode. On the other hand, the network with the convective-heating mode is dried evenly from the top to the bottom zone. Since the top and middle zones remain saturated for long time, an extended quasi-constant drying rate period is obtained for the network with the contact-heating mode (c.f. Figure 6).

2.2 Effective parameters of continuous model assessed from the pore network modeling

As mentioned in the previous section, the heat and mass transfer inside the porous domain can be simulated by using pore network models. However, since the pore network model considers the non-isothermal fluid flow in the pore individually, the number of the discrete equations of the model increases with the size of the domain. Thus, the pore network model cannot be used to simulate the macroscale drying process of porous media which is consists of billions of pores. To describe the drying process of macroscopic porous material, the continuous models developed based on the macroscopic effective transport parameters should be used. These effective parameters are often correlated from a serial of experimental data. To provide a fundamental understanding of the link between the transport coefficient and the fluid transport mechanisms, the transport coefficients are revisited by using the pore network simulation results. For example, the influence of network saturation on the non-local equilibrium effect, the effective moisture transport is investigated in the works of Kharaghani et al. [11, 12, 13] for the isothermal drying process. The incorporation of these effective parameters in the highly non-isothermal drying process is still an open issue. In the future, the impact of non-isothermal pore scale fluid transport on the effective parameter should be taken into account.

3. Continuum-scale models

In continuum-scale models, the porous media is assumed isotropic and homogeneous. Both solid structure and fluid properties are averaged in a representative elementary volume (REV) based on the volume averaging technique. The size of REV should not large enough to avoid the fluctuations of the morphological properties due to micro heterogeneity [14, 15]. The energy and mass conservation equations are derived based on the first physical principles using effective parameters.

3.1 Diffusion model

In diffusion model, the mixture of liquid water and vapor water flow is considered as single phase named moisture and the gradient of moisture content is solely driving force of moisture transport. Since the water diffusion flow leads to the evolution of water concentration in the control volume, the mass conservation of water is derived as

$$\frac{\partial {\rho }_{0}X}{\partial t}+\frac{\partial }{\partial x}\cdot \left[-D_{\mathit{eff}}\frac{\partial }{\partial x}\left({\rho }_{0}X\right)\right]=0.$$E11

In Eq. 11, the apparent density of the dry porous medium where the void volume is taken into account is computed as ${\rho }_0=\frac{M_s}{V}$ (kg dry solid/m³), $X=\frac{M_l}{M_s}$ (kg water/kg dry solid) is the dry-based moisture content of the wet solid. D_eff (m²/s) is the effective diffusivity of moisture in the porous medium. Similarly, the energy conservation equation is written as the change of energy density caused by the enthalpy flow which consists of contributions from diffusive heat flow and heat conduction, namely,

$$\frac{\partial }{\partial t}\left[\left({\rho }_0{c}_{p,s}+{\rho }_0{c}_{p,l}X\right)T\right]-\nabla \cdot \left[D_{\mathit{eff}}{c}_{p,l}T\nabla \left({\rho }_{0}X\right)\right]-\nabla \cdot \left[{\lambda }_{\mathit{eff}}\nabla \left(T\right)\right]=0.$$E12

In Eq. 12, $c_{p,s}$ and $c_{p,l}$ (J/kg.K) denote the specific heat capacity of solid and liquid water, respectively; ${\lambda }_{\mathit{eff}}$ (W/m.K) denotes the effective thermal conductivity of the porous medium.

To perform the diffusional model simulation, the effective thermal conductivity and effective moisture diffusivity needed to be known. Several methods such as fitch method, laser flash, hot-disk method and hot-wire method, can be used to measure the effective thermal conductivity directly whereas measurement of effective moisture diffusivity is rather more complex. Generally, the effective diffusivity of wetted porous materials often decreases during the water dehydration process since the liquid water is strongly captured in the small pores under a higher capillary action compared to large pores [16, 17]. Recently, Khan et al. [18] reported that the moisture diffusivity for isotropic porous materials, therein food products, can be described as a function of moisture content as

$$D_{\mathit{eff}}=D_{\mathit{ref}}{\left(\frac{1+X_0}{1+X}\frac{{\rho }_l\left(1+X\right)+{\rho }_{0}X}{{\rho }_l\left(1+X_0\right)+{\rho }_0{X}_0}\right)}^2,$$E13

where D_ref (m²/s) is the reference diffusivity, X₀ is the initial moisture content. The question raised here is how to determine the reference effective diffusivity D_ref. Instead of measurement directly, it can be computed numerically by an optimization routine using invert method. The objective optimizing function $g\left(D_{\mathit{eff}}\right)$ is the sum of square of difference between numerical moisture content $X_{i,\mathit{sim}}$ and experimental $X_{i,exp}$ as

$$g\left(D_{\mathit{eff}}\right)=\sum _{i=1:\mathit{end}}{\left(X_{\exp }\left(t_i\right)-X_{\mathit{sim}}\left(t_i\right)\right)}^2,$$E14

where t_i denotes the time interval i of the data sampling. By minimizing the value of objective, the numerical moisture evolution is fitted to the experimental observation and the value of effective diffusivity is obtained. For an example, this method is applied for single wood particle drying process in superheated steam environment. The detail description of the geometrical information, initial and boundary conditions is presented in Le et al. [19]. The effective diffusivity of 6.633 × 10⁻⁹ m²/s determined using the experimental data observed at a drying temperature of 140°C can be used reliable for a large range of drying temperature as shown in Figure 7. It implies that the moisture diffusivity model presented in Eq. 13 can be used resonable for wooden porous media.

Figure 7.

The evolution of drying history during optimization process for D_ref determination (a) and experimental and numerical mean moisture content evolutions over time (b).

The validated diffusion model of wood drying is implemented in a CFD model is used to investigate the drying behavior of a pilot dryer. The sketch of dryer (L = 1000 mm, D = 3000 mm, and H = 500 mm) where 30 plates of wood (L = 900 mm, D = 200 mm, and H = 25 mm) is presented in Figure 8. The air with a temperature of 140°C and moisture content of 7 g ram vapor/kg dry air flows into the dryer with a velocity of 0.1 m/s. The initial moisture and temperature of wood are 0.61 kg water/kg dried solid and 20°C. The 2D velocity profile of drying agent is presented in Figure 8 and the evolution of moisture content of plate number (1), (2), (3), and (4) are plotted together with the mean moisture content of all plates in Figure 9. The results indicate a noticeable moisture content maldistribution of wood plates which should be remedied by wood plate disturbing during the drying process.

Figure 8.

Pilot dryer model and velocity profile inside the dryer.

Figure 9.

Moisture content evolution over time of wood.

3.2 Whitaker’s continuum model

In the 1970s, the simultaneous heat, mass, and momentum transfer for non-isothermal drying process is proposed by Whitaker based on the volume averaging technique. In this model, all pore-level mechanisms for heat and mass transfer are considered. The liquid flow due to capillary forces, vapor and gas flow due to convection and diffusion are taken into account in this model. Using the volume averaging technique, the properties of fluid and solid phases such as velocity, density, pressure are averaged in REV. Afterward, the macroscopic differential equations were defined in terms of average field quantities. The detail descriptions of Whitaker’s model can be found in Vu et al. [20, 21]. In this section, the mass and energy conservation equations are briefly recalled.

Mass conservation equation of liquid water

$$\frac{\partial }{\partial t}\left({\rho }_l{\epsilon }_l\right)+\nabla \cdot \left({\rho }_l{v}_l\right)+{\dot{m}}_v=0$$E15

Mass conservation equation of water vapor

$$\frac{\partial }{\partial t}\left({\rho }_v{\epsilon }_g\right)+\nabla \cdot \left({\rho }_v{v}_g\right)-\nabla \cdot \left({\rho }_g{D}_{\mathit{eff}}\nabla \left(\frac{{\rho }_v}{{\rho }_g}\right)\right)-{\dot{m}}_v=0$$E16

Mass conservation of water in both vapor and liquid phases (sum of Eqs. 15 and 16)

$$\frac{\partial }{\partial t}\left({\rho }_l{\epsilon }_l+{\rho }_v{\epsilon }_g\right)+\nabla \cdot \left({\rho }_l{v}_l+{\rho }_v{v}_g\right)-\nabla \cdot \left({\rho }_g{D}_{\mathit{eff}}\nabla \left(\frac{{\rho }_v}{{\rho }_g}\right)\right)=0$$E17

Mass conservation of air in the gas phase

$$\frac{\partial }{\partial t}\left({\rho }_a{\epsilon }_g\right)+\nabla \cdot \left({\rho }_a{v}_g\right)-\nabla \cdot \left({\rho }_g{D}_{\mathit{eff}}\nabla \left(\frac{{\rho }_a}{{\rho }_g}\right)\right)=0$$E18

Energy conservation equation

$$\begin{array}{l}\frac{\partial }{\partial t}\left({\rho }_s{\epsilon }_s{h}_s+{\rho }_l{\epsilon }_l{h}_l+{\rho }_v{\epsilon }_g{h}_v+{\rho }_a{\epsilon }_g{h}_a\right)+\nabla \cdot \left[{\rho }_l{h}_l{v}_l+\left({\rho }_v{h}_v+{\rho }_a{h}_a\right)v_g\right]\\ -\nabla \cdot \left[h_a{\rho }_g{D}_{\mathit{eff}}\nabla \left(\frac{{\rho }_a}{{\rho }_g}\right)\right]-\nabla \cdot \left[h_v{\rho }_g{D}_{\mathit{eff}}\nabla \left(\frac{{\rho }_v}{{\rho }_g}\right)\right]-\nabla \cdot \left[{\lambda }_{\mathit{eff}}\nabla \left(T\right)\right]=0\end{array}$$E19

In these equations, ρ_l, ρ_v and ρ_s (kg/m3) are the mass density of liquid, vapor and air, respectively. ${\epsilon }_l=\frac{V_l}{V}$, ${\epsilon }_g=\frac{V_g}{V}$ are the volume fraction of the liquid and gas phases. The porosity of the medium is computed as $\psi ={\epsilon }_l+{\epsilon }_g$ and ${\epsilon }_s=1-\psi$ denotes the volume fraction of the solid phase. The superficial velocity of the liquid and gas phased $v_l$ and $v_g$ (m/s) are calculated as

$$v_l=-\frac{{\mathit{KK}}_{r,l}}{{\mu }_l}\nabla p_l\text{and}{v}_g=-\frac{{\mathit{KK}}_{r,g}}{{\mu }_g}\nabla p_g$$E20

where K (m²) is the absolute permeability of the porous medium. $K_{r,l}$ and $K_{r,v}$ are the relative permeabilities of the liquid and gas phases, which are respectively calculated as $K_{r,l}={\left(\frac{{\epsilon }_l}{\psi }\right)}^2$ and $K_{r,g}={\left(1-\frac{{\epsilon }_l}{\psi }\right)}^2$ [21, 22]. In Eq. 19, $h_s$, $h_l$, $h_v$ and $h_a$ (J/kg) are the specific enthalpy of solid, liquid, vapor and air, respectively. Using the assumption of constant specific heat capacity, the specific enthalpy of these components can be calculated respectively by $h_s=c_{p,s}\left(T-T_{\mathit{ref}}\right)$, $h_l=c_{p,l}\left(T-T_{\mathit{ref}}\right)$, $h_a=c_{p,a}\left(T-T_{\mathit{ref}}\right)$ and $h_v=c_{p,v}\left(T-T_{\mathit{ref}}\right)+\Delta h_{\mathit{evp}}\left(T_{\mathit{ref}}\right)$ at reference temperature T_ref = 0°C.

The Whitaker model’s is suitable for the porous media made by impermeable solid phase. This model is applied to describe the drying process of wood material [23]. A good agreement between experimental and numerical data, shown in Figure 10, indicates the validity of the model.

Figure 10.

The distribution of moisture content (a), temperature (b), vapor pressure (c) over drying time for wood particle drying at 120 °C and the comparison between experimental and numerical mean moisture content obtained with different drying temperature (d).

For porous media structured by permeable solid, the liquid water is not only transported in the void space, but it also diffuses in the solid phase. In these hygroscopic porous materials, the liquid water is available in the form of both bound and free water:

$$M_l=M_{l,f}+M_{f,b}$$E21

where $M_l$, $M_{l,f}$ and $M_{f,b}$ denote the mass of total, free and bound liquid water, respectively. As a result, the total moisture content is the sum of the free and bound moisture content $X=X_{l,f}+X_{f,b}$. It should be noted that the bound water is mainly accumulated inside the cells of cellular porous products. A small amount of liquid water is located in small pores inside the solid matrix. In Le et al. [24], the bound water is assumed to be accumulated in the small pores of the macro-solid matrix. Thus, the bound water is removed by both local evaporation and liquid diffusion. The mass conservation equation of liquid can be rewrite as

$$\frac{\partial }{\partial t}\left({\epsilon }_{l,f}{\rho }_l+{\overline{\rho }}_b\right)+\nabla \cdot \left({\rho }_l{v}_{\mathrm{l},\mathrm{f}}+\overline{{\rho }_b{v}_b}\right)+{\dot{m}}_v=0$$E22

For the superheated steam drying process of rice seed, the bound water diffusivity can be calculated as

$${\overline{D}}_b={\overline{D}}_{b,\mathit{ref}}{\left(\frac{{\rho }_w+\frac{X}{1+X}{\rho }_s}{{\rho }_w+\frac{X_0}{1+X_0}{\rho }_s}\right)}^2\text{with}{\overline{D}}_{b,\mathit{ref}}=0.71\times 10^{-9}{\mathrm{m}}^2/\mathrm{s}$$E23

The internal structure of rice seed obtained by μ-CT and the comparison between the experimental observations and numerical results are presented in Figure 11. It indicates the validity of the proposed model.

Figure 11.

The morphology of rice seed (a) and the numerical temperature and moisture content evolution over time obtained from diffusion simulations (b).

Beside the porous material made by permeable solid phase, we consider the drying process of cellular plant product (i.e. fruits and vegetables) which is comprised of several types of cells such as parenchyma, collenchyma, and solute-conducting cells shown in Figure 12 [25, 26]. In these materials, the water is mainly accumulated in the cell space. For several cellular product; i.e. eggplant, cucumber; the amount of intracellular water make up more than 95% of total water [27]. During the drying process, both intracellular and extracellular water is removed. The extracellular water is transported to the medium surface due to both capillary flow and internal evaporation. The changing of extracellular water content leads to the different of water potential between the cell space and intercellular void space which is the driving force of the intracellular water transport across the cell membrane as shown in Figure 12. A advance heat and mass transfer model for superheated steam drying process of cellular porous material was developed in Le et al. [28]. The conservation equations of extracellular and intracellular water are recalled as

Figure 12.

Pictorial representation of extracelluar liquid water (in blue) and intracellular liquid water (in green) removal (reprinted from [28] with the permission).

$${\rho }_l\frac{\partial {\epsilon }_{l,\mathit{ex}}}{\partial t}+{\rho }_l\nabla .\left(v_{l,\mathit{ex}}\right)=-{\dot{m}}_v+j_{w,1}{A}_v\frac{{\epsilon }_{l,\mathit{ex}}}{{\epsilon }_{l,\mathit{ex}}+{\epsilon }_{v,\mathit{ex}}}$$E24

$$\frac{\partial \left({\rho }_v{\epsilon }_{v,\mathit{ex}}\right)}{\partial t}+\nabla .\left({\rho }_v{v}_{v,\mathit{ex}}\right)={\dot{m}}_v+j_{w,2}{A}_v\frac{{\epsilon }_{v,\mathit{ex}}}{{\epsilon }_{l,\mathit{ex}}+{\epsilon }_{v,\mathit{ex}}}$$E25

$$\begin{array}{l}{\rho }_l\frac{\partial {\epsilon }_{l,\mathit{ex}}}{\partial t}+{\rho }_l\nabla .\left(v_{l,\mathit{ex}}\right)+\frac{\partial \left({\rho }_v{\epsilon }_{v,\mathit{ex}}\right)}{\partial t}+\nabla .\left({\rho }_v{v}_{v,\mathit{ex}}\right)\\ -j_{w,1}{A}_v\frac{{\epsilon }_{l,\mathit{ex}}}{{\epsilon }_{l,\mathit{ex}}+{\epsilon }_{v,\mathit{ex}}}-j_{w,2}{A}_v\frac{{\epsilon }_{v,\mathit{ex}}}{{\epsilon }_{l,\mathit{ex}}+{\epsilon }_{v,\mathit{ex}}}=0\end{array}$$E26

In Eqs. 24–26, ${\epsilon }_{l,\mathit{ex}}=\frac{V_{l,\mathit{ex}}}{V}$, ${\epsilon }_{v,\mathit{ex}}=\frac{V_{v,\mathit{ex}}}{V}$ are the volume fractions of liquid water and water vapor in the intercellular void space, respectively. The volumetric evaporation flux is denoted by ${\dot{m}}_v$(kg/m³s). A_v (m²/m³) is the volumetric specific area of the medium, which is the exchange area between cells and intercellular void space in a unit volume. The term $A_v\frac{{\epsilon }_{l,\mathit{ex}}}{{\epsilon }_{l,\mathit{ex}}+{\epsilon }_{v,\mathit{ex}}}$refers to the cell surface area wetted by the liquid, where $\psi ={\epsilon }_{l,\mathit{ex}}+{\epsilon }_{v,\mathit{ex}}$ denotes the intercellular porosity of the medium. The diffusive water flux through the cell walls $j_{w,2}$ (kg/m²s) serve as vapor sources. The mass conservation equation for the extracellular water vapor can thus be written as

$$\frac{\partial \left({\rho }_v{\epsilon }_{v,\mathit{ex}}\right)}{\partial t}+\nabla .\left({\rho }_v{v}_{v,\mathit{ex}}\right)={\dot{m}}_v+j_{w,2}{A}_v\frac{{\epsilon }_{v,\mathit{ex}}}{{\epsilon }_{l,\mathit{ex}}+{\epsilon }_{v,\mathit{ex}}}$$E27

After rearranging of Eqs. 24 and 25, the overall mass conservation equation for the water in the intercellular void space can be written as

$$\begin{array}{l}{\rho }_l\frac{\partial {\epsilon }_{l,\mathit{ex}}}{\partial t}+{\rho }_l\nabla .\left(v_{l,\mathit{ex}}\right)+\frac{\partial \left({\rho }_v{\epsilon }_{v,\mathit{ex}}\right)}{\partial t}+\nabla .\left({\rho }_v{v}_{v,\mathit{ex}}\right)\\ -j_{w,1}{A}_v\frac{{\epsilon }_{l,\mathit{ex}}}{{\epsilon }_{l,\mathit{ex}}+{\epsilon }_{v,\mathit{ex}}}-j_{w,2}{A}_v\frac{{\epsilon }_{v,\mathit{ex}}}{{\epsilon }_{l,\mathit{ex}}+{\epsilon }_{v,\mathit{ex}}}=0\end{array}$$E28

We assumed that the water diffusion across the cell membrane is the mechanism of cell water removal, the mass conservation equation for intracellular water inside the cell-matrix reads

$${\rho }_l\frac{\partial {\epsilon }_{l,\text{in}}}{\partial t}+j_{w,1}{A}_v\frac{{\epsilon }_{l,\mathit{ex}}}{{\epsilon }_{l,\mathit{ex}}+{\epsilon }_{v,\mathit{ex}}}+j_{w,2}{A}_v\frac{{\epsilon }_{v,\mathit{ex}}}{{\epsilon }_{l,\mathit{ex}}+{\epsilon }_{v,\mathit{ex}}}=0$$E29

where ${\epsilon }_{l,\mathit{\text{in}}}=\frac{V_{l,\mathit{\text{in}}}}{V}$ denotes the volume fraction of cell water.

Using the local thermal equilibrium assumption, the heat balance equation of the control volume reads [26, 28].

$$\begin{array}{l}⟨\rho c_p⟩\frac{\partial T}{\partial t}+\left(c_{p,l}{\rho }_l{v}_{l,\mathit{ex}}+c_{p,v}{\rho }_v{v}_{v,\mathit{ex}}\right)\nabla \cdot T\\ =\nabla \cdot \left({\lambda }_{\mathit{eff}}\nabla T\right)-\mathrm{\Delta }{h}_{\mathit{evp}}\left({\dot{m}}_v+j_{w,2}{A}_v\frac{{\epsilon }_{v,\mathit{ex}}}{{\epsilon }_{l,\mathit{ex}}+{\epsilon }_{v,\mathit{ex}}}\right)\end{array}$$E30

In conservation equation system, the diffusive water fluxes through the cell membrane, i.e. $j_{w,1}$ and $j_{w,2}$ (red and violet arrows in Figure 12), need to be known. These diffusive fluxes are determined by using the concept of the water potential $\varphi$ (Pa) and water conductivity of the cell membrane as

$$j_{w,1}=-{\rho }_l{k}_p\left({\varphi }_{w,\mathit{\text{in}}}-{\varphi }_{l,\mathit{ex}}\right)$$E31

$$j_{w,2}=-{\rho }_l{k}_p\left({\varphi }_{w,\mathit{\text{in}}}-{\varphi }_{v,\mathit{ex}}\right)$$E32

The water potential $\varphi$ is computed as the difference of the total energy of one kilogram water molecules compared to liquid water molecule energy at the free liquid – water interface at 1 bar [25]. For the intracellular water, the water potential ${\varphi }_{w,\mathit{\text{in}}}$, which was empirically determined in [29] as a function of the moisture content ${\varphi }_{w,\mathit{\text{in}}}=f\left(X_{\mathit{\text{in}}}\right)$. The potential of liquid water in the intercellular void space ${\varphi }_{l,\mathit{ex}}$ is calculated as the capillary pressure, i.e. ${\varphi }_{l,\mathit{ex}}=-p_c$. The water vapor potential ${\varphi }_{v,\mathit{ex}}$ is computed from the gas theory based on the Gibbs free energy [30].

$${\varphi }_{v,\mathit{ex}}=\frac{\tilde{R}T{\rho }_l}{{\tilde{M}}_w}ln\frac{p_{v,\mathit{ex}}}{p_{v,\mathit{sat}}}$$E33

This model has been used to describe the superheated steam drying process of potato. As can be seen in Figure 13, the experimental observations can be reflected fairly by numerical results. It implies the predictive potential of the proposed model. Furthermore, the extracellular moisture is removed very soon when the dry process commences. The drying process is mainly driven by the dehydration of intracellular water by the diffusion mechanism. Although the Whitaker’s continuum approach is effectively used to predict the drying behavior, the parameters of the continuum model is seemingly difficult to be measured. Thus, in the future, with the help of pore-scale model, these effective parameters should be theoretically extracted from the pore level simulations.

Figure 13.

Simulated and experimental evolution of the normalized total, intracellular and extracellular moisture content over time for potato drying in superheated steam environment at 180°C (replotted from [28] with the permission).

4. Conclusions

In this chapter, some recent advance models on heat and mass transfer during non-isothermal drying process of porous media have been reviewed. It indicates that the drying process at pore-scale and macro-scale has been thoroughly investigated. However, a system study on the bridges between the pore-, macro- and plant scales shall be carried out. For example, several effective macroscopic parameters such as the capillary pressure curve or the relative permeabilities may be obtained from a post-processing of the Monte-Carlo simulation results of pore-scale modeling. The moisture diffusivity under the thermal effect also should be revisted and compared with the insothermal drying process. This upscaling statergy can help on both the product and system design, operation, and optimization.

Acknowledgments

This work is funded by the VIETNAMESE MINISTRY OF EDUCATION AND TRAINING (MOET) under research project B2021-BKA-12.