936

zaMP

The Effect of a Finite Response Time Upon the Propagation of Sinusoidal Oscillations of Fluids in Porous Media By PE~t~lJs A. C. RAATS, Agricultural Research Service, United States Department of Agriculture, Madison, Wisconsin, USA 1. I n t r o d u c t i o n

In exploring the limitations of the classical models describing the flow of fluids in porous media, most investigators have concentrated on the possibility that the relationship between flux and driving force is nonlinear. In the last decade some have realized that in transient flow local imbalance of the pressure of a fluid m a y be a more important source of deviation from the classical models (e.g., E1-51). The central idea in the new model that has been studied by this group is the division of a fluid phase in two parts, (1) fluid in a system of large pores, such as pores between aggregates, and (2) fluid in a system of small pores, such as pores in aggregates. The new features in the model are the ratio of the capacities of the two pore systems to store fluid and the finite re@onse time associated with the exchange of fluid between the two pore systems. The model can easily be adapted to describe a variety of physical problems, such as diffusion of a component of a fluid phase in a porous medium (cI., [5]) and conduction of heat in a structured medium. Thus far the implications of the new model have been studied mainly for systems of different geometry subject to step changes in the fluid pressure or flux at the boundary (e.g., I2-51). This paper deals with the propagation of a sinusoidal pressure oscillation at a plane boundary into a semi-infinite porous medium. The analysis shows some striking differences between the classical and the new model. It also suggests experimental methods for determining the various physical parameters in the new model. Since the partial differential equations describing the flow according to the new model, i.e., equations (9) and (10) or (12) and (44) of this paper, are linear, particular solutions of simple boundary value problems m a y be superposed to obtain solutions of more complex problems. In this sense the particular solution discussed here in detail is an important unit in m a n y other problems. Among the numerous possible applications are the propagation of tidal and river level fluctuations into aquifers [6-8], the aeration of soils with entrapped or dissolved gases [9], the response of plants [101 and of soil/plant systems to fluctuations in the environment, and the transfer of heat in structured media. 2. L i n e a r M o d e l f o r F l o w of a F l u i d i n a S t r u c t u r e d P o r o u s M e d i u m

Throughout this paper the subscript 1 refers to the fluid in the system of large pores and the subscript 2 to the fluid in the system of small pores. The balances of

Vol. 20, 1969

SinusoidaI Oscillations of Fluids in P o r o u s Media

937

mass for the fluid in pore systems 1 and 2 are, respectively: 0t + V - 5 1 v l - s l = 0 ,

(1)

002 0--5- + V "~2 v2 - s2 = 0 ,

(2)

where t is the time, V is the vector differential operator, 51 and 5~ are the densities, expressed as mass per unit bulk volume, 5i vi and 52 v~ are the fluxes, expressed as mass per unit bulk area per unit time, and sa and s 2 are the mass supplies expressed as mass per unit bulk volume per unit time. The exchange of fluid between the pore systems 1 and 2 is subject to the principle of conservation of mass, so t h a t : sl + s2 = 0 .

(3)

The key assumption in the model is the assignment of a pressure fit to the fluid in pore system 1 and a pressure p~ to the fluid in pore system 2. Using these pressures, the following constitutive assumptions are adopted: d51 = q d p l ,

(4)

d52

(5)

~2 d p 2 ,

51 vl -

-hi

(6)

Vpi ,

(7)

o2 v2 = 0 , sl =

s2 -

o: (P2 -

Pl) .

(S)

The capacities c 1 and c2, tile permeability k I and the exchange coefficient ~ are assumed to be constants for a given medium. The constitutive assumption (7) implies k2 = 0: it introduces an a s y m m e t r y into the model. The constitutive assumption (8) satisfies the principle of conservation of mass expressed in equation (3). Substitution of (4), (5), (6), (7) and (8) into (1) and (2) gives: pl Cl oOt

_ k~ v 2 p l + ~ (P2 -

#1),

(9)

3t5 2 C2 O t - -

= Cs (#1 - - P2) 9

(10)

T h e model expressed in equations (9) and (10) has been considered in numerous previous studies (e.g., [1-5]). Usually the coefficients q , c~, k 1 and ~ are expressed in terms of the properties of the fluid and solid phases. In order to keep the form of the equations as simple as possible, such more involved expressions for the coefficients are not introduced in this paper. The main purpose of this paper is to present a particular solution of (9) and (10) which satisfies : Pl Ix = 0, y, z, tJ = P0 + A sine) t (11) where x, y, z denotes a rectangular Cartesian coordinate system. The solution to this problem was suggested b y Truesdell's elegant discussion of diffusion of forced sinusoidal

938

PETRUS A. C. RAATS

ZAMP

shearing at the b o u n d a r y of a semi-infinite mass of an incompressible second-order fluid i11, 12]. The connection between these two problems will be discussed later. Some of the results in this paper are analogous to those obtained previously b y Steggewentz, Bosch and Wesseling in their studies of the transmission of tidal waves in an artesian basin in contact with a relatively impermeable layer which m a y exchange water with the artesian basin !6, 71 . In comparing the present development with those of Steggewentz, Bosch and Wesseling the system of large pores corresponds to the artesian basin and the system of small pores to the impermeable layer. Elsewhere I will discuss this particular application in more detail [8]. 3. S o l u t i o n for Pl Eliminating the pressure Ps from equations (9) and (10) gives the differential equation for the pressure Pl : @1 _ at

kt VSpl + c~ k1 c) cT+ c~ ~ c1+ c2 at Vsp~

c2 Cl d2Pl o, c~ + c2 ate 9

(12)

I n (12) appear the diffusivity k l / ( q + ca), the response time c~/o~ and the capacity fraction c~/(q + @ . Of course, the capacity fraction m a y be expressed in terms of the capacity ratio q/c~. Note t h a t (12) is a linear third order partial differential equation. The expression : Pl = P0 + A e - a ' sin ((9 t - b x) (13) satisfies the b o u n d a r y condition for Pl expressed in (11). Substituting (13) into (12) gives : {k 1 a 2 --

k~ bs - 2 k 1 a b Ss + Cl ~o ~s} sin ((9 t - b x)

(14)

+ {ki as ~s - k~ bs Ss + 2 k I a b - c1 (9 - c s (9} cos (co t - b x) = 0 where the dimensionless frequency Ss is defined as:

Cs -

~

(15)

o~

For equation (14) to be satisfied in general the coefficients of sin ( o ) t - b x) and cos (a) t -- b x) must vanish:

kt(a s - b 2 ) - 2 k ~ s a b + c l o $ ~ kl ~2 (a s -

o,

b2) + 2 k~ a b - - q ( 9 -- caco = 0 .

(16) (17)

F r o m equations (16) and (17) it follows that" ab=

2~ ~

+Cl ,

(18)

(19)

VoI. 20, 1969

939

Sinusoidal Oscillations of Fluids in Porous Media

Solving (18) and (19) for a and b gives

(20) (21) The solution of (12) subject to (11) is given by (13), (20) and (21). Expressions analogous to (13), (18) and (19) were given b y WEss~sm~ [7]. For the purpose of presenting (20) and (21) graphically it is convenient to write them in the form: ~+

~

'

(22) (2s)

Plots of 2 k 1 a2/o: and 2 k 1 b~/~ versus ~ for several values of cl/c ~ are given in figures 1 and 2, respectively. For any given values of cl/c 2 and kl/oc the absorption coefficient a is a monotonically increasing function of the dimensionless frequency ~ . For any given values of ~2 and kl/o~ both the absorption coefficient a and the phase shift factor b are monotonically increasing functions of the capacity ratio c~/c2. According to equation (19) the difference between a ~ and b2 is independent of cl/c2. Before deriving the expressions for the pressure P2, the flux ~1 v~, and the mass supplies s 1 and s 2, I will discuss in the next section several limiting cases of the solution just presented.

Figure 1 Absorption as a function of frequency for several values of capacity ratio.

940

PETRUS A. C. RAATS

0 0

M 0 3

2

l

ZAMP

Figure 2 Phase-shift as a function of frequency for several values of c a p a c i t y ratio.

4. A b s o r p t i o n and P h a s e Shift in V a r i o u s L i m i t i n g C a s e s Negligible storage in the system of large pores I f C1 ~

C2,

i.e., if q/c 2 -> O, the differential equation (12) for the pressure Pl

reduces to : c)t

ca

TOt

V2Pl"

(24)

The solution of (24) subject to (11) is given by (13) with the absorption coefficient a and the phase shift factor b satisfying: lira

2 t~1 a 2 _

~

;~

[1/2 -~

lira

2 k1 b2 _ ~

~

[1/2

c,l~-~0

~

t f~--~[~ /

7"2 ~2

(25)

~

(26)

1+ ~ "

Equations (22) and (23) reduce to (25) and (26) if cl/c 2 + O. The curves for q/c 2 = 0 in figures 1 and 2 are in effect plots of equations (25) and (26). The solution of (24) subject to (11) was first given by STEG~EWEXTZ~6, 8] in his study of the transmission of tidal waves in a rigid artesian basin in contact with a relatively impermeable layer which exchanges water with the artesian basin. The solution was also given by MAI~XOVlTZand COLEMAN[13] and equations (25) and (26) were discussed in detail by TRUESDELL [11, 12]. These authors studied the mathematically analogous diffusion of forced sinusoidal shearing at the boundary of a semi-infinite mass of an incompressible second-order fluid. Recently CFIEN and GURTIN [14] incorporated the results of ~arkovitz and Coleman and of Truesdell in a study of heat conduction involving two temperatures. The implications of (25) and (26) with regard to porous media hydrodynamics will be discussed elsewhere [8].

Vol. 20, 1969

941

Sinusoidal Oscillations of Fluids in Porous Media

Low frequency approximation For a given value of the response time c2/a the dimensionless frequency ~a will approach zero if co is sufficiently small. When ~2 + 0 equations (22) and (23) reduce to: lim -2- -k-l a 2 -- lim 2 kl b~

- { c~ + 1} ~2

(27)

while (20) and (21) reduce to: lim a =

lira b =

C2---*o

C~-->o

{ (c1 +c2)~~ }1/2.

(28)

2 k~

For the special case when q/c 2 - i the straight line asymptote represented by equation (27) has been drawn in figure 3. Equation (28) shows that the low frequency approximation is the classical solution of the linear diffusion equation with diffusion coefficient kl/(q + c~) (cf. [151, section 2.6). Therefore, the diffusivity kt/(q + c~) can be determined

from measurements of the absorption coefficient a or the phase-shift factor b at low frequency. If k~ is known such an experiment will yield (c~ + @, and vice versa. In the low frequency range the propagation of the oscillations depends upon the capacity e2 but becomes independent of the exchange coefficient a. A e, execf 8 b, O e

~x~ct end41owfrequencyappr.

~

6 D a, h/ghtrequencyepp~ f 13,h/ghfrequencyopgl:, F e endAveryhighrrequencyeW.

5 5 e, dngsrrom'sproblem

l

/

/ p / / ~/

/

2

.

j~.

~

_.-~.~" ~*~.S.> ~ j ~ . ~

_~-~j~

/

8

Figure 3 Comparison of exact and approximate values of absorption and phase-shift as a function of frequency for c i / c 2 = 1.

The ratio m/b represents the speed at which the pressure oscillation penetrates into the porous medium. Therefore it is interesting to compare the phase-shift factor b as given by equation (23) with the low frequency approximation given by equation (27). Since for ~2 > 0: { (cl/c~+ 1)2+1 + ~'~ (c1 C2/C2)2 }1/2 < (Clw + 1) (29)

942

PETRUS A. C. RAATS

ZA.~IP

and:

r - > 0

(30>

1+ ~

it is evident that for any ~2 > 0 the value of 2 kx b2/c~as calculated from (23) is smaller than the low frequency approximation (27). In other words, for any ~2 the present model involving equation (72) as the differential equation for p~ predicts a greater speed

than the classical theory. High frequency approximation As ~ becomes large ~/(1 + ~) = y + 1, so that (19) reduces to" lim(a2-b 2)-

~

~,--~ 1

(31)

~1 "

Thus the ratio ~/k~ may be determined from the absorption coefficient a and the phase shift factor b at high frequency. As equation (31) indicates and as may be seen by comparing figures 1 and 2, for y + 1 and any c~/c2 the difference between 2 kl a2/z and 2 kl b~/c~ approaches 2, while equations (22) and (23) reduce to: y-+l

O~

~

~-1

-L ( el : 2 1 1 1 / 2 \

,-~1

~

-

~+1

c~

+1=

/)

+ ( - - ic~/ \ , J

c~2

--1=

+1

)

2

+1

+~ +~

1/2

+1

'

--1

(32) (33)

where the dimensionless frequency ~1 is defined as: e1 ~o

(34)

Plots of 2 kl a2/o~ and 2 hi b2/e versus ~2 according to (32) and (33) with q/c2 = 1 are given in figure 3. Note that for small values of ~2 equations (32) and (33) differ very much from the exact equations (22) and (23). If cl/c 2 4= 0 and ~2 -> 00 equations (32) and (33) may be approximated by: lim

,0-~

2 kx a~ ~ lim ~

cJc2 ~ 0

= ~--,oo

2 k~ b~ ~ ( c~ ) ~

=

~

~~ = ~

(35)

c~/c2 4- 0

or:

lim a ~ lira b =~ ~ cl ~o ~1/2

~-~oo

c~/c2 * 0

= ~-~oo

c,/c~ :# 0

[ 2 k 1j

"

(36)

The straight line asymptote represented by equation (35) has been drawn in figure 3 for q/c 2 = 1. Equation (36) shows that if cl/c 2 is finite the very high frequency approximation is the classical solution of the linear diffusion equation with diffusion coefficient kl/q (el. [151, section 2.6). Note that the capacity ca does not appear in (35) and (36).

Comparison with Angstrbm's solution for a rod that loses heat by radiation to a surrounding medium In terms of the notation of this paper ANOSTI~6W [16] (see also reference [15], section 4.4) assumes that p2(x, t) = P0. The system of equations (9) and (10) then

Vol. 20, 1969

943

Sinusoidal Oscillations of Fluids in Porous Media

reduces to the single equation: Opl = k

q ~-

02Pl

~ 0x~ - c~p~.

(38)

A~GSTR6M showed that the solution of (38) subject to (11) is given by (13) with the absorption coefficient a and the phase-shift factor b satisfying: a b -

q o) _ 2 kt a 2 _ b2

c~ ;. 2 k~ ~1 ,

(39) (40)

-

kl '

c~ - =

1 + \ cz / j + l = { l + Ct ~2 2

-- 1 -- {1 + r

+1,

(41)

-- 1.

(42)

Equations (18), (19), (22) and (23) reduce to (39), (40), (41), and (42) if r is large, so that ~/(1 + $~) = y + 1, and: 1

c~

1 + $22 < c~ < 1.

(43)

A n g s t r S m ' s s o l u t i o n a p p l i e s over a more n a r r o w range o f c o n d i H o n s t h a n the p r e v i o u s l y a p p r o x i m a t i o n : equation (43) expresses the additional restriction. While equations (31) and (40) are the same, (32) and (33) differ from (41) and (42). The capacity % does not appear in Angstrfm's solution. As was noted before, equation (31), and thus (40), implies that the ratio e / k ~ , may be determined from the measurement of the absorption coefficient a and the phase shift factor b. In addition to this, equation (39) may be used to calculate the diffusivity k l / q from the measurement of a and b at a certain frequency co. d i s c u s s e d high f r e q u e n c y

5. Solution for P2 Elimination of the pressure Pl from equations (9) and (10) gives the following differential equation for the pressure P2: op~ Ot

kl V2p2 + c1 + c e

c~ kl o ~ c~ + c 2 Ot V 2 p 2 -

c~ q ~ - c1 + c 2

o~p~ O/2-"

(44)

Comparison of equations (12) and (44) shows that the differential equations for Pl and P2 are the same. However Pl and pz satisfy different boundary conditions. From (10) and (11) it follows that P2 (x = O , y , z, t) must satisfy: @2

~

~

~ A

sine) t .

(45)

Equation (45) is a linear, first-order, ordinary differential equation. Noting that exp (~ t/c2) is an integrating factor for this equation, its solution is readily found to be : P2 + C e -~t/c~ = P0 + A sine) t - $2 eosco t 1+ r

(46)

where C is a constant of integration, which must be determined from the initial

944

PETRUSA. C. RAATS

ZAMP

condition. For sufficiently large t the influence of the initial condition will be negligible, so that, after using some trigonometric identities, equation (46) reduces to: A P2 (x = O, y, z, t) = 750 + (1 + ~)1,2 sin(e) t -- arctg~2) .

(47):

The solution of (44) subject to (47) is: A P2 = 751o+ (1 + ~)~,2 e-~x sin (o) t -- b x - a r c t g ~ ) .

(48)

A result analogous to (48) was obtained b y WESSELING [71. Expressions for a and b have been given already. F r o m (11), (13), (47), and (48) it follows t h a t at any point, including x = 0, the ratio of the amplitudes of p~ and of p~ is given b y : amplitude of p~ = (1 amplitude of p~

+ ~.~)1[2

(49)

The amplitude ratio is independent of q/c 2 . E q u a t i o n (49) was given before b y CHEN and GURTIN [14] in their discussion of the special case when q/c 2 = 0. An i m p o r t a n t implication of (49) is that the re@ome time c2/o~ can be determined from a measurement of the ratio of the amplitudes of Pl and p~ caused by a sinusoidal oscillation of known frequency co. As ~2 becomes large the amplitude ratio approaches to ~z. Equations (47) and (48) also show that at any x the pressure P2 lags behind the pressure Pr b y an angle arctg ~ . The reduction in amplitude and the phase lag of P2 as compared with pl indicate t h a t the oscillatory motion tends to bypass the fluid in the system of small pores. 6 . F l u x in P o r e S y s t e m 1 and the R a t e of E x c h a n g e b e t w e e n Pore S y s t e m s 1 and 2

The flux ~h vl of phase 1 calculated from (6) and (13) is: ~1 vl - hi A e-a x {a sin (e) t

-

b x) + b cos (~ t -- b x ) ) .

(50)

Using some trigonometric identities (50) m a y be put in the form: ~ol vl = hi A e- a " (a S + b~)lj2 s i n @ t - b x The factor (a 2 + bZ)1/2 m a y be evaluated from

(22)

arctg b + 2 ) "

(51)

and (23)"

+ (a~ + b~)lj~ = { CO(cl C2) {.1 -~-{CI/(C11+@~c')}2 ~ }112}1/2. k

(52)

F o r given values of o~, c~ and ca a smaller value of 0~will result in a larger value of ~2Since q/(c 1 + c2) < 1, such an increase in ~2 will, according to (52), result in a decrease .of (a ~ + bY)1/2. According to (20) an increase in ~2 will also result in a decrease of 9e x p ( - - a x ) . Therefore, equation (51) implies t h a t for given values of eo, q and ca .a decrease of the exchange coefficient a will cause a decrease of the amplitude of the flux ~1 v]. Paraphrasing TRUESDELL [11], who arrived at a similar conclusion for the :special case with c1 = O, a finite ~ reduces the capacity of the fluid in the porous medium .to be oscillated.

VoI. 20, 1969

Sinusoidal Oscillations of Fluids in Porous Media

945

The rate of exchange of liquid between the pore systems 1 and 2 calculated from (8), (13) and (48)is : sl

--s2=-aAe

'2(1+~)-

{cos(cot-bx)+~2sin(cot-bx)} (53)

[l+,~/

sin

-- b x -- arctan*2

2)

As ~ becomes large, so that ~/(1 + ~) = 7 + 1 and arctan ~2 + ~/2, this expression for the mass supplies reduces to -- ~ Pl.

7. Summary This paper presents an analysis of the propagation of sinusoidal pressure oscillations at a plane boundary into a semi-infinite, structured porous medium. The central idea in the model is the division of a fluid phase in two parts, (1) fluid in a system of large pores, and (2) fluid in a system of small pores. The ratio of the capacities of the two pore systems to store fluid and the finite response time associated with the exchange of fluid between the two pore systems are important parameters in the model. If tile frequency of the oscillation is small, the heterogeneity of the medium becomes unnoticeable and the ratio of the permeability and the total capacity Call be determined from either the absorption coefficient or the phase-shift factor. At any frequency the finite response time causes the pressure oscillation to penetrate into tile porous medium at a greater speed than would be expected from classical diffusion theory. At high frequency the difference of the squares of the absorption coefficient and the phase shift factor is equal t~ the ratio of the parameter describing the exchange of fluid between the two pore systems and the permeability of the system of large pores. Expressions have been derived for the pressure and the flux in the system of large pores, the pressure in the system of small pores, and the rate of exchange between the two pore systems. The response time associated with the exchange between the two pore systems can be determined from a measurement of the ratio of the amplitudes of the pressures caused by a sinusoidal oscillation of known frequency. Among the numerous possible applications are the propagation of tidal and river level fluctuations into aquifers, the response of plants and of soil-plant systems to fluctuations in the environment, aeration of soils with dissolved and/or entrapped gases, and the transfer of heat in heterogeneous media.

Acknowledgment I wish to thank Prof. C. TRUESDELL for pointing out the connection between the flow of fluids in structured porous media and the flow of second-order fluids. REFERENCES [1] G. I. BARENBLATa: and Yu. P. ZHELTOV, Fundamental Equations o/ Filtration o/ Homogeneous Liquids in tTissured Rocks, Doklady Akademii Nauk SSSR 732, 545548 [English translation, 522-5251 (1960). [2] G. I. BARENBLATT, YU. P. ZHELTOVand I. N. I4~OCHINA, Basic Concepts in the Theory o[ Seepage o[ Homogeneous Liquids in Fissured Rocks [Strata], Priklad. Mat. Mekh. 2d, 852-864 [English s 1286-1303] (1960}. ZAMP 20/60

946

PETRUS A. C. RAATS

ZA~IP

[3] R. C. OOODKNIGHT,W. _4.. KLIKOFF and I. FATT, Non-steady-state Fluid Plow and Di//usion in Porous Media Containing Dead-end Pore Volume, J. Phys. Chem. (Ithaca) 64, 1162-1168 (1960). [4] J. E. WARREN and P. J. RooT, The Behavior o~ Naturally Fractured Reservoirs, Trans. Soc. Petrol. Engin. 228, 245-255 (1963). [5] J. R. PHILIP, Diffusion, Dead-end Pores and Linearized Absorption in Aggregated Media, Austral. J. Soil Res. 6, 21-30 (1968). [6] J. H. STEGG~WENTZ, De invloed van de gelijbeweging van zeegn en getijrivieren op de sti]ghoogte van her grondwater (The Influence o/ the Tidal Motion o/ Seas and Tidal Rivers on the Height o/Rise o/the Groundwater), Thesis (Delft, The Netherlands 1933). [7] J. WESSELING, The Transmission o/ Tidal Waves in Elastic Artesian Basins, Neth. J. Agr. Sci. 7, 22-32 (1959). [8] P.A.C. IRAATS,Horizontal Transmission o/Pressure Fluctuations in Structured Porous Media (in preparation). [9] J.A. CIJRRIX,Gaseous Diffusion in the Aeration o/ Aggregated Soils, Soil Sci. 92, 40-45 (1961). [10] ]3. KLEPPER, Diurnal Pattern o/ Water Potential in Woody Plants, Plant Phys. 43, 1931-1934 (1968). [11] C. TRUESDt~LL, The Natural Time o/ a Visco-elastic Fluid." its Signi/icanee and Measurement, Phys. Fluids 7, 1134-1142 (1964). [121 C. TRUESDELL and W. NOLL, The Non-linear Field Theories o/Mechanics, in: Handbuch der Physik (S. Fliigge, ed.) Ill/3 (Springer, ]3erlin), p. 1-602 (in particular section 123) (1965). [1.3] H. MARKOVlTZ and ]3. D. COLEMAN, Nonsteady Helical Flows o/Second-order Fluids, Phys. Fluids 7, 833-841 (1964). [14] P. J. CItEN and M. E. GURTIN, On a Theory o/ Heat Conduction Involving two Temperatures, Z. angew. 2Xlath. Phys. 79, 614-627 (1968). [15] H. S. CARSLAWand J. C. JAEGER, Conduction o/Heat in Solids (Oxford at the Clarendon Press), 510 p. (1959). [16] A. J. ANGSTR6M, Neue Methode, das Wiirmeleilungsvermdgen der Kdrper zu bestimmen, Ann. Phys. und Chem. (PoggendoHf) (4/) 24, 513-530 (1861) [English translation: New Method of Determining the Thermal Conductibility o/ Bodies, The London, Edinburgh and Dublin Phil. Mag. J. Sci. (d) 25, 130-142 (1863)], Zusammen/assung

In dieser Arbeit wird die Fortpflanzung sinusoidaler Druckschwankungen in einem por6sen K6rper mit complexer Struktur betrachtet. (Received: April 28, 1969.)