Ultrasound in wet materials subjected to drying: A modeling study

International Journal of Heat and Mass Transfer 84 (2015) 998–1007

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Ultrasound in wet materials subjected to drying: A modeling study S.J. Kowalski ´ University of Technology, Institute of Technology and Chemical Engineering, Department of Process Engineering, ul. Berdychowo 4, 60-965 Poznan ´ , Poland Poznan

a r t i c l e

i n f o

Article history: Received 11 November 2014 Received in revised form 15 January 2015 Accepted 15 January 2015

Keywords: Ultrasound Drying materials Wave decomposition Modeling

a b s t r a c t The aim of this paper is to present wave phenomena induced by the external action of ultrasound on porous wet materials subjected to drying. The purpose of these studies is to analyze the distribution of ultrasonic waves in such a complex medium and to discover the mechanism of ultrasonic interaction with both the solid skeleton and moisture in pores. The results obtained here should allow to state how ultrasonic waves are distributed in such material depending on ultrasonic frequency and on resistance between relative motion of the constituents. This knowledge may help to explain enhancement of the drying mechanism by ultrasound, and in particular in biological products such as fruits and vegetables where a significant elevation of temperature is not required. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction The problem of ultrasound (US) wave propagation in porous media has been a subject of interest in many sectors of science (e.g. geophysics, environmental science, hydrology) and engineering (e.g. civil engineering, oil industry), and others [1,2]. Recently, an increasing amount of literature has reported the very positive influence of higher-power ultrasound at frequencies 20–100 kHz on the drying efficiency of biological materials such as fruits and vegetables [3–7]. These reports show the ability how ultrasound enhances drying rate and improve quality of food products. However, they did not explain the heart of the US mechanism contributing to process intensification and improvement of product quality by drying with the help of ultrasound, which is the interest of present considerations. Judging from the published literature [8–11], high-power ultrasound are capable of improving heat- and mass-transfer processes in drying materials, and in particular in the drying of heat-sensitive biological materials such as fruits and vegetables. However, one should state that drying processes enhanced with ultrasounds have still not emerged from the laboratory research phase. The current quest is for the most economical versions of hybrid technological processes combining the ultrasound method with other drying techniques (convective, microwave, infrared) in such a way as to complement one another and to increase the moisture removal rate. The limited applicability of ultrasound to drying methods is primarily due to ultrasound high-power demand. The existing E-mail address: [email protected] http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.01.086 0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.

currently ultrasound radiators operating in gas media (e.g. siren) mostly do not fulfil this requirement. An original construction of a laboratory hybrid dryer based on an ultrasound airborne transducer was constructed in the author’s drying laboratory. This dryer enables hybrid drying in combinations of convective, microwave and ultrasound techniques. Fig. 1 presents, as a forerunner, the kinetics of convective drying enhanced with airborne ultrasound, which was applied to drying a layer of 16 apple slices of dimensions 40  20  5 mm, realized on this new equipment [12]. The moisture ratio (MR) as a function of time (drying curve) presented in Fig. 1 reveals that the drying time amounted to 235 min for pure convective drying (curve 1), while for drying with ultrasonic assistance at ultrasound power of 100 W and frequency 26 kHz it took only 160 min (curve 2). Besides, the ultrasonic assistance elevated the sample temperature only by about 1 °C on average as shown in Fig. 1, which is of significant importance in the drying of temperature sensitive biological materials. The main aim of this article is, first of all, to recognize the interaction mechanism between the ultrasonic wave and biological material, which could entail the reason for the intensification of moisture removal from dried products. The research hypothesis is a supposition that the periodical waves which are characteristic of ultrasound cause periodical changes of porosity and pore pressure. In this way they may evoke the moisture streaming from inside the material towards the surface, where it evaporates. Biological materials such as fruits and vegetables need very sublimed drying methods as they are very sensitive to temperatures higher than 60–70 °C and also to long drying time. The common drying methods (e.g. the hot convective air drying), may cause degradation of their valuable features (color, vitamin, minerals).

S.J. Kowalski / International Journal of Heat and Mass Transfer 84 (2015) 998–1007

For this reason the author of this paper has developed hybrid methods which are a combination of convective, microwave and infrared drying [13–17]. This work proposes an extension of hybrid drying additionally on ultrasound assistance. This needs, however, first a more detailed analysis of ultrasound wave distribution in wet materials subjected to drying. The positive outcome of these studies may contribute to essential changes important in the drying technology of biological material. Modified drying technology carried out in dryers supported with ultrasound equipment could find application in the industry, increase competitiveness by increasing productivity and decreasing energy consumption, and contribute to balanced research and development (R&D). 2. Theory 2.1. Fundamentals The theory of ultrasound propagation in wet porous materials is constructed on the basis of the mechanics of continua. Natural tools for development of such a theory are the balances of mass, momentum, moment of momentum, energy, entropy and the principles of irreversible thermodynamics. A version of such a drying theory was presented in [18, Chapter 3 and 4]. In this paper the dynamic terms such as inertia forces are introduced to the former theory which allow to extend it to ultrasonic enhancement of drying processes [19–21]. The present model of drying based on the mechanics of continua incorporates the following assumptions:  The considered material is a saturated porous body consisting of a solid skeleton (a = s) and moisture (a = m), which is a mixture of liquid (a = l) and gas (a = g) in pores.  Individual components are represented by mass concentrations and volume fractions, which are continuous functions of space and time due to an averaging procedure applied to the representative volume element (RVE)⁄.  The theory includes dynamic terms, such as accelerations, inertia forces, kinetic energy, etc. in the momentum and energy balance equations.

999

 Stress in moisture is represented by pressure only, as the stress deviator in moisture is considered insignificant with respect to that in the solid skeleton.  The skeleton is a deformable body and the deformations resulting from the ultrasounds are assumed to be small (strain fluctuations). The following definitions are used in the further considerations: Volume fraction /a – is a fraction of the RVE occupied by the constituent a (s-solid, m-moisture, being a mixture of l-liquid and g-gas). Porosity / – is a fraction of the RVE occupied by pores. One assumes that the pore space is filled by the moisture (liquid/ gas mixture), i.e. /  /m = /l + /g = 1  /s. Body saturation u – is a fraction of pore space occupied by the liquid phase, u = /l//. Constituent mass concentration qa = qar /a [kg/m3] – is the mass of constituent a in the RVE per volume of this element. True (intrinsic) constituent mass density qar [kg/m3] – is the mass of constituent a per volume of this constituent. Production of constituent mass ⁄qa [kg/m3 s] – is the rate of mass change between liquid and gas due to phase transitions within the RVE. Constituent velocity va [m/s] – is defined as the volumetric flux of constituent a per unit area of the body (filter velocity). 2.2. Balance equations The balance equations are developed by applying Euler’s description using the spatial coordinates x{x, y, z}. A control volume V(t) separated within the drying body is attributed to the deformable skeleton and enveloped in a smooth control surface A(t) oriented spatially with the outward-directed unit normal vector n, Fig. 2. The individual constituents (phases) located primarily in volume V at time t are displaced during drying with different velocities va and take different configurations V0 and V00 after time increment Dt. Consider an arbitrary physical quantity Wa(t) (mass,

Fig. 1. Drying curves and temperature of apple samples: 1-pure convective drying, 2-drying with ultrasound assistance.

1000

S.J. Kowalski / International Journal of Heat and Mass Transfer 84 (2015) 998–1007

a

qa

d ua ¼ raji v ai;j  qai;i þ qa r a þ ^ea dt

ð5Þ

qa

 a a q d sa qa ra ¼  ia þ a þ ^sa dt T ;i T

ð6Þ

In these equations da/dt denotes the substantial derivative with convection velocity of constituents a. Masses, momentums and energies of individual constituents do not have to be conserved, however, total mass, momentum and energy of the multicomponent medium as a whole must be conserved. The laws of conservation express the following constraints:

X

qa ¼ 0  for moisture constituents; qs B0  for the skeleton

a

ð7aÞ Fig. 2. Deformable control volume of saturated porous body.

X

momentum, energy, etc.) of constituent a that is present in the control volume V(t)

a

Z

Wa ðtÞ ¼

qa wa ðx; tÞdV

ð1Þ

 X a ^ ai þ qa v ai ¼ mi ¼ 0  for momentum m

  X X 1 ^ea þ m ^ ai v ai þ qa ua þ v ai v ai ¼ ea ¼ 0  for energy 2 a a ð7cÞ

VðtÞ

where qawa(x, t) is the density of Wa(t). The total balance of quantity Wa(t) reads

Z

@ @t

qa wa dV þ

V

I

qa wa v a  ndA ¼ A

I

Z waA dA þ qa waV dV A V Z a þ w dV

ð2Þ

V

The first integral on the left-hand side expresses the local variety of Wa(t) in time, the second integral describes the convection of this quantity; the first integral on the right-hand side expresses the supply of Wa(t) through the surface, the second integral expresses a volumetric supply of this quantity, and the third integral is the source of Wa(t). Individual balance equations of mass, momentum, energy, and entropy can be obtained by substituting suitable expressions as given in Table 1. Table 1 presents the definition of wa for specific balances, where: ua and sa denote the specific internal energy and entropy of constituent a, raji ; qai , ra, gi, and Ta denotes the partial stress tensor, heat flux vector, volumetric energy supply (radiation), gravity acceleration, and temperature; and q⁄a, mi a , e⁄a, and s⁄a denote the production terms (sources) of mass, momentum, energy and entropy, respectively. The balance of moment of momentum is not considered here. The only implication following from this balance is the symmetry of total stress tensor, which is true for non-magnetic media. The assumption of continuum allows us to write the following local forms of Eq. (2) for the balance of mass, momentum, energy, and entropy a

d qa þ qa v ai;i ¼ qa dt

qa

a a

d

vi

dt

ð3Þ

^ ai ¼ raji;j þ qa g i þ m

ð4Þ

Table 1 Component quantities of the individual balance equations. Balance

wa

waA

waV

w⁄ a

Mass Momentum

1

0 gi

q⁄a

Energy

ea ¼ ua þ 12 v ai v ai

0 tai ¼ raji nj  raji v ai  qaj nj  a a  qi =T ni

Entropy

ð7bÞ

a

v ai s

a

M i a

g i v ai þ r a a a

a

q r /T

e⁄a s⁄ a

The entropy of the multicomponent medium is conserved only in one special case, namely, for reversible processes. Drying processes, however, should be considered irreversible, so that the total entropy increases during drying and obeys the second law of thermodynamics which states that total entropy production is a positive quantity

X X ð^sa þ qa sa Þ ¼ sa P 0 a

ð7dÞ

a

By replacing the term qara in the entropy balance (6) with that separated from the balance of energy (5) one obtains



    a a a a ^ea 1 ad f a ad T a a a 1 þ q q þ q s  r v þ a ¼ ^sa a a ij i;j i dt dt T T ;i T

ð8Þ

where fa = ua  saTa denotes the free energy of a-constituent. 3. Simplified theory Bearing in mind examples of drying materials such as fruits and vegetables, one can now consider them as coherent bodies consisting of the solid skeleton and moisture (Fig. 3). Its structure is characterized mostly by small-sized pores (micropores) filled with liquid and containing gas bubbles which are initially almost invisible and which may grow in the course of drying. Differentiation of skeleton and moisture temperatures in such a porous material seems to be quite exaggerated. Therefore, it is assumed in further considerations that the constituent temperatures are equal in a given point of the material (Ta = T). Such an assumption significantly simplifies the theory (see Fig. 4). Furthermore, when considering the dynamic effects accompanying the propagation of ultrasounds in wet porous materials, it seems reasonable to consider two components, i.e. skeleton and moisture, instead of three components: skeleton, liquid, and gas. This is because the solid skeleton carries the greatest deal of dynamic forces applied to the medium, and the moisture, being a mixture of liquid (incompressible phase) and gas bubbles (compressible phase), has dynamic properties when is considered as one component. The moisture is described physically by density qm, velocity vm, and the spherical stress tensor rmij ¼ Pm dij ¼ ppor /dij , all of which are defined as follows:

qm ¼ ql þ qg ¼ qlr /l þ qgr /g ¼ qpor /

ð9aÞ

1001

S.J. Kowalski / International Journal of Heat and Mass Transfer 84 (2015) 998–1007

Fig. 3. Visualization of apple as a saturated porous material.

where mint can be interpreted as the skeleton–moisture interaction i body force, while the force of the opposite sign is exerted by moisture upon the skeleton. The partial stress tensor for the moisture is formally expressed now as follows

X

Pm ¼

 Pa  qa v ai v ai ¼ ppor / where v ai ¼ v ai  v m i for a ¼ l; g

a¼l;g

ð13Þ For further reasoning it is convenient to multiply the momentum balance equations (12a) by v si and (12b) by v m i , and integrate them by parts over control volume V. The sum of these integrals gives

Lint þ K ¼ Lext

ð14Þ

The mechanical energy balance (14) states that the external work increment Lextdt provided during time dt by the surface and body forces acting on the matter inside the control volume is transformed into kinetic energy Kdt and internal mechanical energy Lintdt, where

Lint ¼

Z h

Fig. 4. Ultrasonic waves of different velocity propagating in a saturated porous domain.

V

K¼

qm v m ¼ ql v l þ qg v g

ð9bÞ

 Pm ¼ Pl þ Pg ¼  pl /l þ pg /g ¼ ppor /

Lext ¼ ð9cÞ

The gas bubbles may change their dimensions periodically as a result of compression and decompression of moisture caused by periodical changes of pore pressure due to ultrasound waves. The pore pressure is a mean pressure of liquid and gas depending on body saturation u = /l// por

p

/l /g ¼p þ pg ¼ pg þ ðpl  pg Þu ¼ pg þ pcap u / / l

ð10Þ

where pcap = pl  pg denotes capillary pressure as the difference between liquid and gas pressure. The mass balance equations for skeleton and moisture are expressed as follows: s

d qs þ qs v si;i ¼ 0 dt

ð11aÞ

m

d qm þ qm v m i;i ¼ 0 dt

ð11bÞ

The momentum balance equations for the solid skeleton and moisture read s

d v si q ¼ rsji;j þ qs g i  mint i dt s

m

m

q

d

v mi

dt

¼

Pm ;i

m

þ q gi þ

mint i

ð12aÞ

ð12bÞ

Z 

1 2

qs

V

Z



A

i



rsji dsji þ Pm dmii þ mint v si  v mi dV ¼ i

Z

lint dV

 s m  d  s s d  v i v i þ qm v mi v mi dV dt dt

 m tsi v si þ t m i v i dA þ

Z



V

ð15Þ

V

ð16Þ 

qs v si g i þ qm v mi g i dV

ð17Þ

a

and dij denotes the symmetric strain rate tensor, t ai ¼ raji nj is the stress vector of constituent a = (s, m) acting on boundary surface A enveloping control volume V. The mechanical power density lint expresses the power associated with deformation of the matter   s  rsji dsji þ Pm dmii and interaction mint v i  v mi of internal force mint i i when considering the two-component body as a whole. The global balance equation of energy (2) applied for a twocomponent body consisting of the skeleton and moisture as a whole can be written as follows

U þ K ¼ Lext þ Q

ð18Þ

Formula (18) formally presents the first law of thermodynamics which expresses energy conservation when considering the twocomponent medium as a whole, i.e. the constraint (7c) is held. In this equation U denotes the rate of internal energy of the medium as a whole, and Q is the non-mechanical power (heat) supplied to the control volume. The new terms in (18) are expressed as follows

U¼

Z 

qs

V

Q ¼

Z A

 s m d us d um dV þ qm dt dt

qi ni dA þ

Z V

qrdV ¼

ð19Þ

Z



 qi;i þ qr dV

ð20Þ

V

where us and um denote the density of internal energy per unit mass associated with the skeleton and moisture, qi ¼ qsi þ qm i is the overall heat flux; while qr = qsrs + qmrm denotes the rate of

1002

S.J. Kowalski / International Journal of Heat and Mass Transfer 84 (2015) 998–1007

volumetric energy supply (e.g. microwave or ultrasonic radiation) per unit mass and q = qs + qm is the mass density of the medium as a whole. By substituting expressions (16), (17), (19) and (20) into (18) and omitting the volumetric integral one obtains the first law of thermodynamics in the local form s

qs

m

d us d um s m s m þ qm ¼ rsji dji þ Pm dii  qi;i þ qr þ mint i ðv i  v i Þ dt dt

ð21Þ

4. Rate equations for heat and mass transfer Heat and mass transfer phenomena in drying processes are considered to be irreversible, so that expressions (27) and (28) have to be positive independently of each other, i.e. Ts⁄q P 0 and Ts⁄w P 0. By calculating the gradient of chemical potential lm(Pm, T) (25) and including the local form of the first thermodynamic law for moisture one obtains

1

The second law of thermodynamics states that entropy production s⁄ by irreversible processes must always be positive. According to Eqs. (6) and constraint (7d), by assumption that Ts = Tm = T, the balance of entropy for the skeleton and moisture as a whole reads

lm;i ¼ 

d ss d sm qi qr  qs þ qm þ  ¼s P0 T dt dt T ;i

Tsw ¼ 

s

m

 s s m m s df d f dT q s m  qs þ qm þ rsji dji þ Pm dii  i T;i  qs s dt dt dt T   int m þ wm sm T ;i þ im i ¼ Ts P 0

qm ð23Þ

where fa = ua  saT is the free energy of constituent a = (s, m), s = ss + Xmsm is the total entropy of the two-component medium per unit mass of the solid, Xm = qm/qs denotes the mass of moisture m m s content per unit mass of the solid, and wm i ¼ q ðv i  v i Þ is the moisture flux with respect to the skeleton reference frame. The second law in the form (23) is satisfied implicitly when specifying the expressions of the free energies and the constitutive relations linking the partial stresses to the strains and temperatures. For further usage the inequality can be rewritten in a more convenient form to develop the constitutive relations s

 qs f_  qs sT_ þ rij dji þ qs lm X_ m  

¼ Ts P 0 s

 qi lm;i þ sm T ;i þ ami  g i T ;i  wm i T ð24Þ

m m

where f = f + X f denotes the total free energy of the body per unit mass of the skeleton, rij ¼ rsij þ P m dij is the total stress tensor of the saturated porous body, lm is the chemical potential of moisture, and am i is the acceleration of moisture motion. The dot over the symbols denotes the substantial time derivative with convection velocity of the skeleton and the comma between the bottom indexes denotes the partial derivative with respect to space coordinates. Note that the chemical potential lm has appeared in Eq. (24) in a natural way by taking the form

lm ¼ f m 

Pm

q

m

¼ um 

Pm

qm

 sm T

ð25Þ

On the basis of expression (24) one can state that the production of entropy consists of three elements: s⁄ = s⁄f + s⁄q + s⁄w, where s Tsf ¼ qs f_  qs sT_ þ rij dji þ qs lm X_ m

qi T ;i T  ¼ wm lm;i þ sm T ;i þ ami  g i i

ð26Þ

Tsq ¼ 

ð27Þ

Tsw

ð28Þ

The entropy production terms express energy dissipation associated with s⁄f – porous body deformation, s⁄q – heat conduction, and s⁄w – moisture transport.

ð29Þ

h wm i

q

m

 m i m Pm ai  g i P 0 ;i þ q

ð30Þ

Entropy production (30) associated with moisture transport presents the product of moisture flux wm i and the driving forces consisting of the pore pressure gradient and the difference between gravity and inertia forces. Inequality (30) will express a positively defined quadratic form if

wm

q

m Pm ;i  s T ;i

By insertion (29) into (28) one can write

ð22Þ

Eliminating the a priori defined radiation term qr in the inequality (22) by using the balance of energy (21), one obtains the thermodynamic inequality of the form



qm

¼ ðv m  v s Þ ¼ K ½gradP m þ qm ðam  gÞ

ð31Þ

where K is a permeability tensor for the anisotropic body, being of scalar form for an isotropic body. It is dependent on the physical properties of both moisture and the porous body structure. For the isotropic body K is a constant and its physical significance is stated as K = k/g, where k is the specific permeability of the porous medium and g is the viscosity of the moisture [22]. Formula (31) performs the rate equation for the moisture flow, revealed as Darcy’s law if am = 0 and porosity / = const. In order to identify the unknown interaction force mint present in the momentum balance equations (12a,b), the equation of momentum for moisture (12b) is used to replace the gradient of pore pressure gradPm with that existing in the moisture flux (31). After some operations are conducted, one finds

mint ¼

g k

ðv s  v m Þ

ð32Þ

The skeleton–moisture interaction force mint [N/m3] can now be interpreted as a drag force being proportional to the difference between the skeleton and moisture velocities, with the coefficient being the ratio of moisture viscosity g to pore structure k (permeability). Entropy production associated with heat conduction (27) is a positive quantity

Tsq ¼ 

qi T ;i P 0 T

ð33Þ

The expression (33) appears as the product of heat flux and the temperature gradient, so it will be positively defined if

q ¼ k  gradT

ð34Þ

where k > 0 is the thermal conductivity being a tensor for anisotropic or a scalar for isotropic bodies. Formula (34) performs the rate equation known as Fourier’s law. 5. Constitutive equations for a porous body The porous body deformation is irreversible when the skeleton of the porous body reveals viscoelastic properties, or elastic–plastic when the stresses overcome the yield criterion.In this section linear deformations of the material are assumed, which can be both reversible and irreversible. The geometrically linearized theory s assumes that the strain rate tensor dji of the skeleton is equal to s the time derivative of the strain tensor, i.e. dji ¼ e_ sji . The strain

1003

S.J. Kowalski / International Journal of Heat and Mass Transfer 84 (2015) 998–1007 ðrÞ

tensor esji consists then of two parts, the reversible eij and the irreversible s ji

e

ls ¼ f s 

eðirÞ strains of the skeleton ij ð35Þ

where u denotes the displacement vector of the skeleton. The expression concerning dissipation of the mechanical energy (26) can now be written as follows

ðrÞ

D

F ¼ F D ðeij Þ and FS = FS(T, qs, qm) and after setting the time derivative of these functions to (41) and (42), the following equations of state are obtained

@F D

sij ¼

 ðrÞ ðirÞ ðrÞ qs f_  qs sT_ þ sij e_ ji þ e_ ji þ re_ ii þ qs lm X_ m ¼ TsðrÞ þ TsðirÞ P 0 where sij 

ð42Þ

By applying the chain-rule expansion to free energy functions

  1 @ui @uj ðrÞ ðirÞ ¼ eij þ eij ¼ þ 2 @xj @xi

ssij

rs rs ¼ us  s  ss T qs q

!

ðrÞ

@eij

@F S ; S¼ @T

!

l

; qs ;qm

s

@F S ¼ @ qs

! ;

m

l

T;qm

@F S ¼ @ qm

¼ rij  Pdij is the stress deviator of the porous body, ðirÞ

ðirÞ

ðirÞ

ekk dij =3 is the deviator of the irreversible strain. The volumetric

T;qs

ð43Þ

ð36Þ

equal to the stress deviator in the skeleton, and eij ¼ eij

!

By developing the individual state equations with respect to the parameters of state and by using the differential for chemical potentials (25) and (42), one obtains the following total differentials

ðirÞ

irreversible strain is assumed to be zero (ekk ¼ 0), as the spherical part of the stress tensor does not cause irreversible deformations [23]. Thus, the volumetric reversible strain represents the total volðrÞ kk

s kk ).

umetric strain (e ¼ e The mechanical power now splits into the deviatory (shear) power, which concerns both reversible and irreversible strains, and into the volumetric power which concerns the reversible strain only. Thus, it was assumed in relation (36) that the irreversible deformations concern shape deformations only. This is because the irreversible deformations (plastic, viscous) refer first of all to shape deformations. Obviously, dissipation of the mechanical work may concern only irreversible strains, i.e. ðirÞ

sij e_ ji ¼ TsðirÞ P 0

ð37Þ

The inequality (37) will be satisfied if the product of the stress tensor and the irreversible strain rate constitute a positively defined quadratic form. The sufficient condition to fulfill this requirement provides the following relation ðirÞ sij ¼ C ijkl e_ kl

ð38Þ

where Cijkl is a tensor of material coefficients constructed in such a way as to ensure a positive definition of inequality (37). The rate equation (38) refers to plastic or viscoelastic materials. The inelastic deformations are not considered here by ultrasound wave propagation. Therefore, the details concerning inelastic deformation can be found in the relevant literature [23,24]. For the thermo-elastic matrix the intrinsic dissipation is zero. s_

dsij ¼

ðrÞ ðrÞ q f  q sT_ þ sij e_ ji þ re_ ii þ qs lm X_ m ¼ TsðrÞ ¼ 0 s

ð39Þ

@2FD ðrÞ

ðrÞ

ðrÞ

@eij @ekl

dekl

ð44aÞ

dS ¼ ss dqs þ sm dqm 

@ 2 F SP @T 2

! ! @ 2 F SP @ 2 F SP s s m d dq m þ s q  þ s @T@ qs @T@ qm

dT 

ð44bÞ





!

@ 2 F SP @ 2 F SP @ 2 F SP s s þ s d q þ dqm dT þ @ qs @T @ qs @ qm ð@ qs Þ2

1

drs ¼ s

q

1

! @ 2 F SP @ 2 F SP @ 2 F SP m dT þ þ s dqm þ m s dqs 2 m @ q @T @q @q ð@ qm Þ

m

dP ¼ m

q

ðrÞ F_ D ¼ sij e_ ij

ð40Þ

F_ S ¼ ST_ þ ls q_ s þ lm q_ m s

D

S

ð41Þ s s

where F = q f = F + F is the free energy per unit volume, S = q s + qmsm is the total entropy of the body per unit volume, and qa is the mass concentration (partial density) of constituent a = (s, m). The chemical potential for moisture is expressed by formula (25), and that for the solid skeleton, which appeared in equation (42), is expressed as follows

ð44dÞ

The respective second order derivatives represent the material coefficients. They are denoted as follows:

@2FD ðrÞ

ðrÞ

@eij @ekl @ 2 F SP ð@ qs Þ2

qm

¼

¼ Mijkl ¼ 2M s dik djl ; Ks ðqs Þ2

;

@ 2 F SP þ sm @T@ qm

@ 2 F SP ð@ qm Þ2 !

¼



@ 2 F SP

Km ðqm Þ2

@T 2 ;

q

@ 2 F SP Q ¼ ; @ qs @ qm qs qm

¼

Cv ; T

s

@ 2 F SP þ ss @T@ qs

¼ cmðTÞ

! ¼ csðTÞ ;

ð45Þ

Thus, the following physical relations for the drying body in an incremental form are obtained: ðrÞ

For development of physical relations in elastic case it is convenient to separate the mechanical power produced by the stress deviator from that produced by the spherical stress. This is because the changes of temperature and moisture do not influence the shape deformations but volumetric strains only. Then, the free energy in expression (39) splits into the deviatory (D) and spherical (S) parts as follows

ð44cÞ

dsij ¼ 2Ms deij

qs dss þ qm dsm ¼

ð46aÞ Cv dqs dqm dT  csðTÞ s  cmðTÞ m T q q

drs ¼ csðTÞ dT  K s

m

dqs

q

dP ¼ cmðTÞ dT  K m

s

Q

dqm

qm

dqm

qm

Q

dqs

qs

ð46bÞ

ð46cÞ

ð46dÞ

where cs(T) = 3Ksas(T), cm(T) = 3Kmam(T) express the heat and mass transfer parameters, Ks and Km are the volumetric elastic modules for the porous matrix and moisture, Q is the coefficient of the coupling between the volume changes of the solid and that of the moisture, as(T) and am(T) are the coefficients of thermal expansion of the porous matrix and moisture, and Cv is the heat capacity for the medium as a whole, respectively.

1004

S.J. Kowalski / International Journal of Heat and Mass Transfer 84 (2015) 998–1007



Using the mass balance equations (11a,b), one can write

dqs

qs

dqm

¼ v ss;i dt des ;

qm

m ¼ v m i;i dt ¼ de

qm U€ i ¼ Quk;ki þ K m U k;ki  cmðTÞ T ;i þ qm g i þ B u_ i  U_ i ð47Þ

where des and dem denote the volume changes of the skeleton and moisture. It is assumed that the material coefficients (45) are constants in the first approximation and that the material parameters, such as volume strains of the skeleton and moisture and the temperature, undergo only infinitesimal variations with respect to their initial values (fluctuations caused by ultrasounds). Then, the physical relations corresponding to isotropic elastic material can be written as follows: ðrÞ

sij  ssij ¼ s0ij þ 2M s eij S ¼ S0 þ

ð48aÞ

Cv ðT  T 0 Þ þ csðTÞ es þ cmðTÞ em T0

ð48bÞ

m

P ¼

Pm 0

s

þ Qe þ K

m m

e c

mðTÞ

r ¼

ssij

s

þ r dij ¼ r

s 0ij

þ

2Ms sij



s s

ðirÞ

Ts ¼ sij e_ ji 

ð48dÞ

m

sðTÞ

e þ A e þ Qe  c

 ðT  T 0 Þ dij

2 qi g T ;i þ v si  v m i T k

ð52Þ

Substituting (52) into the balance of entropy (22) leads to the following heat balance

  s m d ss d sm ¼ qi;i þ qr þ D T qs þ qm dt dt

ð53Þ

where

ð48cÞ

The sum of physical relations (48a) for the stress deviator and (48c) for the spherical stress in the skeleton gives the full physical relation for the skeleton in the elastic range s ij

In order to develop the differential equation for determination of the temperature, it is necessary to use the thermodynamic expression (24) and eliminate from it the reversible part by using the identity (40) and the expressions (12b) and (32). Thus, one obtains the formula describing the dissipated energy in the form

ðirÞ

ðT  T 0 Þ

ð51bÞ

where B ¼ gk.

D ¼ sij e_ ji þ

rs ¼ rs0 þ K s es þ Q em  csðTÞ ðT  T 0 Þ



g_ k

ui  U_ i

2

ð54Þ

is the mechanical energy which is dissipated and changed into heat. By substituting the differentials in the bracket of Eq. (53) by differentials (46b) and using the rate equation for heat (34), one gets the differential equation describing the evolution of temperature in a wet porous body

 C v T_ þ T csðTÞ u_ k;k þ cmðTÞ U_ k;k ¼ kr2 T þ qr þ D

ð55Þ

where A = K  2M /3 is the bulk elastic module for the porous matrix.

where qr = audPu denotes the heat generated by the ultrasonic waves, au [–] is the absorption coefficient of ultrasonic energy, d [–] is the working efficiency of the ultrasonic transducer, and Pu [W/m3] is the ultrasound power per unit volume.

6. Equations of drying theory with ultrasonic assistance

7. Ultrasonic waves in the drying body

The balance and constitutive equations developed above can serve for an analysis of ultrasound interaction between the skeleton and moisture. The ultrasounds cause harmonic oscillations of material particles around stagnant points. These oscillations (vibrations) are transferred through the material in the form of waves. As the material particles oscillate around the stagnant points, and as such are not displaced in space, one can neglect the convective terms in the substantial derivatives and consider only local time derivatives in heat and mass transport equations. The differential equations describing the ultrasonic waves are based on mass balance equations (11a,b), momentum balance equations (12a,b), entropy balance equations (22), rate equation for heat conduction (34), and the constitutive equations (48b), (48d) and (49). The mass balance equations for the skeleton and moisture can now be rewritten as follows

Below the effect of ultrasound wave propagation through a porous medium whose pores are filled with moisture is considered. The influence of dynamic coupling between the fluid and solid and the US frequency on the phase velocity of ultrasonic waves is determined. Therefore, in this analysis the absence of an external heat supply is assumed. So, the thermodynamic state of the body is changed due to ultrasound action only. Because of the high US frequency f (f = 26 kHz), the process of heat generation due to absorption of ultrasound energy can be considered adiabatic. In such circumstances the temperature gradients as well as gravity forces in Eqs. (51a,b) can be neglected. The most interesting, from the point of view of drying theory, are the fluctuations of constituent densities and strains caused by periodically changeable acoustic pressure. The governed functions in this analysis are the displacements and densities of the solid skeleton and moisture, i.e.

ð49Þ s

s

s

q_ s þ qs u_ i;i ¼ 0

ð50aÞ

q_ m þ qm U_ i;i ¼ 0

ð50bÞ

g ¼ u; U; qs; qm

ð56Þ

Each of these functions is assumed to be of the form

where u and U denote the displacements vectors for the skeleton and moisture, respectively. The dot over the symbol denotes a partial time derivative. The momentum balance equations for the solid skeleton and moisture after substituting the physical relation (48d) and (49) and the interaction force (32) give the following differential equation for displacements ui and Ui

where g denotes the local and instant fluctuation of quantity g around the equilibrium value g ¼ const. In the present consider s and q0m q m. ations is assumed q0s q The system of equations (50a,b) and (51a,b) after application of (57) and the above-mentioned assumptions reads

qs u€ i ¼ Ms r2 ui þ ðMs þ As Þuk;ki þ QU k;ki  csðTÞ T ;i þ qs g i

q_ s0 þ q s u_ 0i;i ¼ 0

ð58aÞ

q_ m0 þ q m U_ 0i;i ¼ 0

ð58bÞ

  B u_ i  U_ i

ð51aÞ

gðx; y; z; tÞ ¼ g þ g 0 ðx; y; z; tÞ

ð57Þ

0

1005

S.J. Kowalski / International Journal of Heat and Mass Transfer 84 (2015) 998–1007

 q s u€ 0i ¼ Ms r2 u0i þ ðMs þ As Þu0k;ki þ QU 0k;ki  B u_ 0i  U_ 0i 

q m U€ 0i ¼ Qu0k;ki þ K m U 0k;ki þ B u_ 0i  U_ 0i



ð59aÞ ð59bÞ

Henceforth, the dilatational waves will be considered as the most relevant for drying purposes. They can involve periodical compression and rarefaction of the drying material and thus may be responsible for cavitation or sponge effect in drying material. By applying the divergence operation to Eqs. (58a,b) and (59a,b) one obtains

q_ 0s þ q s e_ 0s ¼ 0

ð60aÞ

q_ 0m þ q m e_ 0m ¼ 0

ð60bÞ

q s €e0s ¼ r2 ðNs e0s þ Q e0m Þ  Bðe_ 0s  e_ 0m Þ

ð61aÞ

q m €e0m ¼ r2 ðQ e0s þ K m e0m Þ þ Bðe_ 0s  e_ 0m Þ

ð61bÞ

s

s

s

where N = 2M + A is the dilatational elastic modulus in the skeleton, and e0 s = divu0 and e0 m = divU0 denote the dilatational fluctuations caused by ultrasonic waves, which are periodical in character. Following Biot’s suggestion [25,26], a reference state is introduced in which the coupling between the motion of the solid and fluid (moisture) is so strong that the relative motion between the solid and fluid is completely prevented in some way, i.e. u0 = U0 and divu0 = divU0 . In such a state the wave equations (61a,b) are reduced to one equation of the form



  fr  s s asm es0 þ am em i e0  em 0 z þ cm 0  cm em ¼ 0 f

ð67bÞ

where z = (c/v)2 denotes the square of the wave velocity ratio, v = x/ b is the complex wave velocity, and f = x/2p is the frequency of ultrasonic waves. Eliminating the constants es0 and em 0 from (67a,b) yields the algebraic equation for z in the form



 f as am  a2sm z2  ðas cm þ am cs Þz þ cm cs þ cm r iðz  1Þ ¼ 0 f

ð68Þ

In the case of a purely elastic wave, and when fr = 0, the reduced equation (68) has two roots, z1 and z2, expressed as follows

z1 ¼

z2 ¼

as cm þ am cs 

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðas cm  am cs Þ2 þ 4a2sm cm cs

2ðas am  a2sm Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi as cm þ am cs þ ðas cm  am cs Þ2 þ 4a2sm cm cs 2ðas am  a2sm Þ

ð69aÞ

ð69bÞ

Note that if the coupling coefficient asm reflecting the volume interactions between the skeleton and moisture is neglected, then one gets

z1 ¼

z2 ¼

c2

v 21 c2

v 22

¼

¼

cs as

cm am

and

and

v1 v2

sffiffiffiffiffiffi Ns ¼ q s

ð70aÞ

sffiffiffiffiffiffiffi Km ¼ q m

ð70bÞ

where c is the velocity of wave propagation in a composition of two constituents with the absence of relative motion between them, q ¼ q s þ q m is the mass density of such a body, and H = Ns + Km + 2Q characterizes this body’s elastic property. By applying this reference state one can rewrite the equations (61a,b) as follows

It follows from (70a,b) that if the mechanical coupling asm and the frictional coupling fr equal zero, then two waves propagate in the medium, one in the skeleton and the other in the moisture, and their velocities are expressed by (70a) and (70b), respectively. If the mechanical coupling asm is not zero, then from (69a,b) still follows that two coupled waves propagate through the medium as a whole with different velocities, i.e. one is faster than the other. Noting that if z1 < z2 then v1 > v2. Equation (68) after including the solutions (69a,b) reads

cs €e0s ¼ c2 r2 ðas e0s þ asm e0m Þ  2pcm f r ðe_ 0s  e_ 0m Þ

ð63aÞ

z2  ðz1 þ z2 Þz þ z1 z2 þ sðz  1Þi ¼ 0 with

cm €e0m ¼ c2 r2 ðasm e0s þ am e0m Þ þ 2pcm f r ðe_ 0s  e_ 0m Þ

ð63bÞ

€e0s ¼ c2 r2 e0s with c2 ¼

H q

ð62Þ

where the following non-dimensional parameters are introduced s

as ¼

N ; H

m

am ¼

K ; H

asm ¼

Q ; H

q s cs ¼  ; q

q m cm ¼  q

ð64Þ

and

fr ¼

  B 1 m s 2p q

ð65Þ

is a reference frequency. Furthermore, the one-dimensional problem of wave propagation is considered. Thus the dilatational fluctuations are assumed as propagating in x-directions in the form of harmonic attenuated waves, i.e.

e0a ðx; tÞ ¼ ea0 exp iðbx  xtÞ; a ¼ ðs; mÞ

ð66Þ

where x is the frequency of the ultrasonic waves and b = br + ibi is the complex number consisting of the wave number br and the attenuation coefficient bi. Inserting the proposed solution (66) into Eqs. (63a,b) will yield after rearranging the following system of algebraic equations



  fr  s s as es0 þ asm em i e0  em 0 z  cm 0  cs e0 ¼ 0 f

ð67aÞ

s¼

fr cm   f as am  a2sm ð71Þ

Equation (71) has two complex roots, zI and zII, and their form points to the fact that ultrasounds in a porous material filled with viscous moisture split into two waves of different velocity. The complex form of these roots suggests that the waves are attenuated. In order to determine the wave velocities one has to calculate the square roots of zI and zII, which read

pffiffiffiffi pffiffiffiffiffi zI ¼ RI þ II i and zII ¼ RII þ III i

ð72Þ

The relative phase velocities of the first kind and the second kind, V1 and V2, are given by

V1 ¼

v1 c

¼

1 j RI j

and V 2 ¼

v2 c

¼

1 jRII j

ð73Þ

The attenuation coefficients of the first- and second-kind waves are given by

  bi jII j ¼ br I jRI j

and

  bi jIII j ¼ br II jRII j

ð74Þ

The dilatational fluctuations (66) can now be expressed as h h x i x i e0a ðx;tÞ ¼ ea0I exp II x expi RI ðx  v 1 tÞ c c h h x i x i þ ea0II exp III x expi RII ðx  v 2 tÞ for a ¼ ðs;mÞ ð75Þ c c

1006

S.J. Kowalski / International Journal of Heat and Mass Transfer 84 (2015) 998–1007

The wave number br and the attenuation coefficient bi in formula (66) are substituted by

brI ¼ RI

x c

;

brII ¼ RII

x c

;

biI ¼ II

x c

;

biI ¼ III

x c

Note that the amplitudes of the dilatational fluctuations in Eq. (75) are attenuated with distance x. This means that the ultrasonic waves may disappear after some distance from the start point. 8. Results and conclusions By using the MatCad program for numerical calculus one can draw the plots of the velocities and the attenuation coefficients for ultrasonic waves propagating through a porous medium whose pores are filled with moisture. The values of material coefficients used for numerical calculations are presented in Table 1 [25].The velocities and the attenuated coefficients of the first- and second-kind dilatational waves are presented below as a function of the coupling parameter s (Eq. (71)). Note that this parameter is constructed of such quantities as moisture viscosity g (for water 103 Pa s), porous medium permeability k ( 1010 m2 on average), moisture density qm ( 103 kg/m3), and ultrasonic frequency f (>2  104 1/s). These quantities give the ratio fr/f 1, on average [27,28]. Therefore, in the analysis below the following range of the coupling parameter s is assumed

06s¼

fr cm 6 1:5 f ðas am  a2sm Þ

Fig. 6. Phase velocities V2 = v2/c of the second kind dilatational waves as a function of coupling parameter s.

Table 2 Mechanical properties of saturated porous medium. No.

as

am

asm

cs

cm

1 2 3 4

0.60 0.60 0.40 0.70

0.30 0.30 0.60 0.30

0.04 0.04 0.04 0.40

0.35 0.65 0.35 0.35

0.65 0.35 0.65 0.65

Note that s ? 0 when g ? 0 or f ? 1. As follows from Figs. 5 and 6, the velocities of first-kind waves are faster than those of second-kind waves for cases 1, 2 and 4 (Table 2). This is because the modulus of skeleton compressibility is greater than that of moisture. The opposite effect was reached in case 3, where the moisture revealed a greater modulus of compressibility than the skeleton. Besides, the smaller constituent mass density is favorable for greater velocity of the first-kind wave but not for the second-kind wave, as is seen in cases 1 and 2 (Table 2). The velocities of both kinds of waves slow down along with an increase of the coupling parameter s, and those of the first-kind wave (cases 1, 2 and 4) tends to velocity c. The opposite effect is reached in case 3. It is seen that for the growing parameter s the two-component medium tends to a single kinematics. This may take place when there is no relative motion between the constituents because of strong coupling due to high moisture viscosity. Fig. 7. Attenuation coefficients b1 = (bi/br)I of the first kind dilatational waves as a function of coupling parameter s.

Fig. 5. Phase velocities V1 = v1/c of the first kind dilatational waves as a function of coupling parameter s.

Figs. 7 and 8 show the attenuation coefficients of first- and second-kind waves with an increase of the coupling parameter s. It is seen from the plots on these figures that the waves are not attenuated for s = 0. The first-kind waves are attenuated insignificantly with an increase of s, except in case 3 (Fig. 7). The secondkind waves are stronger attenuated (Fig. 8), except that in case 3 characterized with greater modulus of moisture compressibility. It means that the waves of the second kind will disappear faster than those of the first kind, except that in case 3. They may disappear after short time from the start point, so that the distance of their propagation may be insignificant. The waves of the first kind, on the other hand, are weakly attenuated, so the distance of their propagation can be much longer. Based on the above presented continuous drying model the effect of ultrasound waves on moisture transport through drying media will be investigated theoretically and validated experimentally. The increase of moisture flow and the drying rate as well as

S.J. Kowalski / International Journal of Heat and Mass Transfer 84 (2015) 998–1007

Fig. 8. Attenuation coefficients b2 = (bi/br)II of the second kind dilatational waves as a function of coupling parameter s.

several other relevant effects will be studied in such materials as fruits and vegetables. One can mention, for instance, the reduction in adherence between pore wall and the liquid, viscosity decrease due to energy dissipation, acoustic streaming, acoustic cavitation, in-pore turbulence, shrinkage of drying bodies, drying induced stresses, cracks phenomena, etc. The analytical assessment of such effects as ‘‘heating effect’’, ‘‘vibration effect’’, and ‘‘synergistic effect’’ will be presented in separate article [12] and proved that the estimated in this way effects contribute to drying efficiency. Hybrid dryer with the air-born ultrasound transducer of high power up to 200 W and frequency of 26 kHz will be used to carry out the experimental drying test and to validate the agreement between the predicted and experimental data. Conflict of interest No any conflict to publish this work. Acknowledgments This work was carried out as part of a research Project, No. 2012/05/B/ST8/01773, sponsored by the National Science Centre in Poland. References [1] A. S´liwin´ski, Ultrasounds and their Application, WNT, Warszawa, 2001 (in Polish).

1007

[2] M.J.W. Povey, Ultrasonic Techniques for Fluids Characterization, Academic Press, San Diego, London, 1997. [3] A. Mulet, J.A. Carcel, N. Sanjuan, J. Bon, New food drying technologies – use of ultrasound, Food Sci. Technol. Int. 9 (2003) 215–221. [4] H.T. Sabarez, J.A. Gallego-Juárez, E. Riera, Ultrasonic-assisted convective drying of apple slices, Drying Technol. 30 (2012) 989–997. [5] K. Schössler, H. Jäger, D. Knorr, Effect of continuous and intermittent ultrasound on drying time and effective diffusivity during convective drying of apple and red bell pepper, J. Food Eng. 108 (2012) 103–110. [6] J.V. Garcia-Perez, J.A. Carcel, J. Benedito, A. Mulet, Power ultrasound mass transfer enhancement in food drying, Food Bioprod. Process. 85 (2007) 247– 254. [7] Ó. Rodríguez, J.V. Santacatalina, S. Simal, J.V. Garcia-Perez, A. Femenia, C. Rosselló, Influence of power ultrasound application on drying kinetics of apple and its antioxidant and microstructural properties, J. Food Eng. 129 (2014) 21– 29. [8] J.V. Santacatalina, J.V. García-Pérez, A. Mulet, J.A. Cárcel, Combined effect of temperature and air velocity on low-temperature drying of apple by power ultrasound, in: Proceedings of 19th International Drying Symposium (IDS 2014), Lyon, France, 24–27 August 2014. [9] J.A. Carcel, J.V. Garcia-Perez, J. Benedito, A. Mulet, Food process innovation through new technologies: use of ultrasound, J. Food Eng. 110 (2012) 200–207. [10] A.C. Soria, M. Villamiel, Effect of ultrasound of the technological properties of bioactivity of food: a review, Trends Food Sci. Technol. 21 (2010) 323–331. [11] C. Ozuna, A. Mulet, J.V. García-Pérez, E. Riera, J.A. Cárcel, Influence of vegetable structure on ultrasonically assisted convective drying, in: Proceedings of 19th International Drying Symposium (IDS 2014), Lyon, France, 24–27 August 2014. [12] S.J. Kowalski, A. Pawłowski, Intensification of apple drying due to ultrasound enhancement, J. Food Eng. (in press). [13] S.J. Kowalski, K. Rajewska, Effectiveness of hybrid drying, Chem. Eng. Process. Process Intensification 49 (2009) 1302–1309. [14] S.J. Kowalski, G. Musielak, J. Banaszak, Heat and mass transfer during microwave-convective drying, AIChE 56 (1) (2010) 24–35. [15] S.J. Kowalski, Chemical engineering science, control of mechanical processes in drying, Theory Exp. 65 (2010) 890–899. [16] S.J. Kowalski, K. Rajewska, Convective drying enhanced with microwave and infrared radiation, Drying Technol. 27 (7–8) (2009) 878–887. [17] S.J. Kowalski, D. Mierzwa, Hybrid drying of red bell pepper: energy and quality issues, Drying Technol. 29 (10) (2011) 1195–1203. [18] S.J. Kowalski, Thermomechanics of Drying Processes, Springer Verlag, Heilderberg-Berlin, 2003. [19] P. Piyasena, E. Mohareb, R.C. McKellar, Inactivation of microbes using ultrasound: a review, Int. J. Food Microbiol. 87 (2003) 207–216. [20] Ye. Yao, Using power ultrasound for regeneration of dehumidizers in desiccant air-conditioning systems: a review of prospective studies and unexplored issues, Renew. Sustainable Energy Rev. 14 (2010) 1860–1873. [21] S.J. Kowalski, Ultrasonic waves in diluted and densified suspensions, Ultrasonics 43 (2004) 101–111. [22] A.E. Scheidegger, The Physics of Flow through Porous Media, third ed., University of Toronto Press, 1974. [23] R. Hill, Mathematical Theory of Plasticity, University Press, Oxford, 1971. [24] Y.C. Fung, Fundamental of Solid Mechanics, PWN, Warszawa, 1969 (in Polish). [25] M.A. Biot, Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low frequency range, J. Acoust. Soc. Am. 28 (1956) 168–178. [26] M.A. Biot, Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range, J. Acoust. Soc. Am. 28 (1956) 179–191. [27] S.J. Kowalski, Engineering of Porous Materials, Poznan´ University of Technology Press, Poznan´, 2004 (in Polish). [28] G.A. Aksielrud, M.A. Altszuler, Mass Transport in Porous Materials, WNT, Warszawa, 1987 (in Polish).