A scheme for determining the acoustic and physical properties of saturated porous media

Journal of Sound and Vibration (1991) M(3),

A SCHEME PHYSICAL

FOR

DETERMINING

PROPERTIES

OF

M. F. Department

397-407

YEO

THE

SATURATED AND

ACOUSTIC

AND

POROUS

MEDIA

L. J. SCHMID

of Civil Engineering, The University of Adelaide, Adelaide, South Australia 5001

(Received

G.P.O. Box 498,

11 July 1989, and in revised form 27 March 1990)

A new technique, with existing finite elements, but with different elements occupying the same volume in space, has been used to model the acoustic and physical properties of saturated porous media. This new overlay technique allows the interphase drag between the solid and fluid phases to be modelled accurately. 1. INTRODUCTION constants of any medium are required to enable the theoretical acoustic speeds to be determined. The acoustic speeds in single-phase homogeneous materials are easily obtained as the material constants are usually uncomplicated. However, in twophase materials, and in particular for saturated porous media, which is the subject of this paper, the situation is more complex. Experimental and theoretical analyses have been reported that determine the material constants in saturated porous rock under dynamic excitation. References [ 1- 141 are typical and include consideration of equivalent material properties, boundary conditions, levels of saturation, loss mechanisms and acoustic properties. Many other references exist, particularly in the geophysical journals, as this topic is of particular interest to oil and mining companies carrying out seismic surveys, to enable a sensible geological picture to be obtained. In this paper a new approach is presented for obtaining the material properties and wave speeds in saturated porous media by using finite element analysis. An overlay technique is utilized in which two elements effectively occupy the same volume of space, one representing the fluid and the other the porous solid. Each element can act independently of the other, thus allowing relative motion between the fluid and solid phases. This overlay technique is new and has never been reported before. The technique accurately models the physical action of the flow of the fluid through the solid but does not directly attempt to model any chemical effects. Spencer [13], indicated that his experimental results could not be completely explained due to physical interaction but could indicate a chemical effect being superimposed upon the physical effects. The technique reported in this paper can model Spencer’s results exactly if two critical constants in the analysis are suitably selected. These constants are (a) the interphase drag and (b) the effective bulk modulus of the fluid. Spencer’s experiments were modeiled theoretically by Dunn [ 141. This paper therefore initially compares the results obtained from the overlay technique with Dunn’s results. It is acknowledged by Dunn in his theoretical modelling of Spencer’s experiments that some of the physical constants were not measured but were the theoretical book values. It is possible that the values used by Dunn were different from those of the actual The physical

397 0022-460X/91/030397+ 11/$03.00/O

0 1991 Academic

Press Limited

398

M. F. YE0

AND

L. J. SCHMID

specimens used by Spencer. Therefore by suitably adjusting the drag coefficient and bulk modulus, Dunn’s theoretical work is reproduced and subsequently Spencer’s experimental results are reproduced. The analysis is then extended to calculate the wave speeds in the saturated composite by using transfer matrix theory. 2. THEORETICAL OVERVIEW In a finite element static analysis the forces and stresses generated by nodal displacements in the fluid and saturated solid are caused by two interacting effects. Firstly, if the solid matrix is strained without a pressure build-up in the fluid then the solid nodal forces are due solely to the elastic constants of the solid, with no fluid nodal forces but usually with fluid nodal displacements. Secondly, if there is a pressure generated in the fluid due to an inflow of fluid, then the hydrostatic stresses generated in both the fluid and solid are due to the composite bulk modulus of the solid and fluid combination. This composite bulk modulus is required because the change in fluid pressure will cause a hydrostatic compression in the particles of the solid as well as in the fluid. A similar situation can occur when there is a change in the volume of the solid without fluid motion. This would change the pore volume and hence generate a fluid pressure which would then induce an identical pressure in the solid. In the dynamic situation in a saturated porous solid it is possible for there to be fluid motion independent of the motion of the solid. If a volume of the composite is defined by, say, a rectangular element then each of the nodes must have four degrees of freedom: (a) solid displacements in the x and y directions; (b) fluid displacements in the x and y directions. It is easier to visualize the physical situation if the concept of two overlayed elements is introduced. By overlaying is meant the occupation of the same volume by a solid element and by a fluid element. For this paper the nodes will be coincident. This concept lends itself ideally to the use of two existing elements: the eight-noded isoparametric solid element, to represent the porous solid [ 151 and the eight-noded isoparametric fluid element to represent the fluid [16]. In addition to the existing element matrices it is necessary to generate another matrix to model the drag between the two phases. In general the fluid displacements are uncoupled from the solid displacements other than through drag coupling between the phases. In this context drag coupling is used to imply all various types of coupling between the phases such as viscous drag and drag due to the shape and size of the pores. When an impervious layer is introduced the fluid displacements perpendicular to the interface are identical to those of the porous solid at the same point of contact. For an inviscid fluid the displacements parallel to the interface are not locked to those of the solid, but are free to move according to the mechanics of the situation. Both these effects can be achieved by modifying the element definitions to link the appropriate displacements. 3. ELEMENT MATRICES The element matrices can be obtained by considering the nodal displacements solid and the fluid separately and superimposing the effects.

of the

3.1. SOLID NODES DISPLACING If the solid nodes displace while the fluid nodes are held in position, i.e., there is no fluid motion, then two sets of forces can be generated. Firstly due to the constitutive

PROPERTIES

relationship

OF

of the solid the nodal

SATlJRATED

POROUS

forces generated

MEDIA

are given by

(1) where K, = j BTDB d( uol), 6, are solid displacements, 6, are fluid displacements, F” are forces at solid nodes, and FJ are forces at fluid nodes. K, is obviously the standard elastic element stiffness [ 151, D relates stresses to strains and B relates strains to displacements. Secondly, a volume change may occur in the solid, and hence in the pore volume, due to the solid nodes displacing. The total volumetric change, AV, causes a reduction in the pore volume of +AV, where 4 is the porosity. This implies that either an equal volume of fluid must flow out or else a pressure is generated in the fluid and hence in the solid. As the fluid nodes have been prevented from moving there can be no fluid outflow. The volumetric change is then @AV/cjV = A V/ V and the resulting pressure is given by composite. pressure = K, A V/ V, where K, is the effective bulk modulus of the fluid/solid The fluid element stiffness K,- can then be used to determine the fluid pressure and the associated nodal forces. Because of the porosity, part of the hydrostatic force is transmitted through the solid matrix material and part through the fluid. The solid to fluid transfer ratio is obviously (1 - 4): 4. The forces caused by the volume change in the solid are therefore given by 2) Unfortunately 3.2.

FLUID

this will result NODES

in the final matrices

being

unsymmetric.

DISPLACING

If there is a volume change in the fluid element, AV, due to the fluid nodes displacing, then the actual volume change is 4A V. This will cause a pressure change of K,I#JA V/ +V = K,AV/ V. Again this pressure exists in both the fluid and solid, and the nodal forces transmitted are given by (3) Combining

equations

K,, can be complex 3.3.

MASS

(l), (2) and (3) gives the element

stiffness of the overlayed

~K+W&M,

= l=;I.

if the Young’s

“---i”l modulus

I;;/

is complex.

K, is real for inviscid

system:

fluids.

MATRICES

mass matrices for finite element analysis have been well documented [ 171. The basic equation is [MI = m 1 NTN d(vol), where /Ml is the element mass matrix, m is the mass per unit volume and INI is the element shape function. Since the composite can be considered as two overlayed elements with independent displacements the mass matrix for the overlayed pair is The

4

[Ml=

I

NTN d( uol) 0

0

ml

where m,v is the mass of the solid per unit volume unit volume of the solid.

NTN d( uol) and m, is the mass of the fluid in a

400 3.4.

M. F. YE0 FLUID/SOLID

DRAG

AND

L. .I. SCHMID

MATRIX

When the fluid moves relative to the solid an interphase force exists between them. For saturated porous flow in this paper it is assumed that the force is proportional to the product of the relative velocity and the drag coefficient. The drag coefficient has therefore not only to represent the viscous forces but any inertial effects due to the pore size and tortuosity. Changes in any of the parameters would cause a change in the drag coefficient. It is therefore necessary to derive a relationship between the drag coefficient and the other parameters. In the case of steady state harmonic excitation the velocity is given by 8 = iw6, where w is the radian frequency and i = v’?. The nodal forces are therefore obtained from IFI=iWc,

1 N ‘( relative

displacement)

where C, is the drag coefficient between the phases. If the shape functions for the solid and fluid elements equation (4) becomes:

d( vol),

(4)

are iV, and N, respectively

then

= iwcn

The matrix then becomes a component of the imaginary part of the stiffness matrix. The units of the drag coefficient are force per unit volume per unit relative velocity. The calculation of the drag coefficient is complex and depends upon the porosity, the shape and size of the voids, the viscosity of the fluid and associated boundary layer thickness. Because the programs used for the analysis are ignoring boundary layer effects an equivalent drag coefficient must be obtained for the porous solid/fluid considered. This could be obtained experimentally or possibly through D’Arcy’s equation. For the purpose of this paper the value will be obtained by selecting a value to produce, for one of Dunn’s analyses, the correct peak attenuation frequency. This value will then be used for the remaining analyses as other parameters are varied.

4. SOLUTION The solution The governing

technique equation,

is nothing

TECHNIQUE

more than a simple

~K--w*M+iwC~)6~=~F~,

finite element

dynamic

analysis.

(9

is solved in the usual manner, where K is the real part of the stiffness matrix, w is the radian frequency, M is the mass matrix, C is the imaginary part of the stiffness matrix, S is the displacements vector, and F is the nodal forces vector. Spencer’s saturated cylinder, 140 mm long and 38 mm in diameter, was modelled by using axisymmetric elements. The mesh is fixed axially at one end and given an oscillatory axial displacement at the other end. All nodes were free to move radially (see Figure 1 (a)). The forces resulting at the oscillating nodes due to the oscillations are summed and divided by the cross-sectional area to give the complex stress. Dividing this stress by the applied strain gives the complex Young’s modulus as defined by Spencer. The ratio of the imaginary to real parts of the Young’s modulus gives the attenuation constant Qs as defined by Dunn.

PROPERTIES

OF SATURATED

POROUS

401

MEDIA

Rubber membrane,

(b) Figure

Dimensions

1. Mesh used in mm.

to model

a cylinder

5. COMPARISON

of saturated

WITH

sandstone

EXISTING

subjected

to oscillatory

axial

strains.

WORK

The proposed technique was initially tested by using the theoretically produced graphs of Dunn. He produced a theoretical scheme for calculating the attenuation and Young’s modulus of various porous media including Navajo sandstone, and compared his results with the experimental work of Spencer. Basically, a 38 mm diameter cylinder of saturated sandstone, 140 mm long, was fixed at one end and subjected to axial vibrations by shaking the other end of the cylinder. Experimentally, the maximum frequency was 400 Hz, while theoretically vibrations up to 100 000 Hz were applied. For the purpose of this test the axisymmetric finite element mesh used was that shown in Figure l(b). Initial runs were carried out to determine the mesh refinement for convergence. The most economic mesh was three overlayed elements radially for the saturated sandstone and a single element axially. The sandstone was surrounded by a thin lens of water, one element 0.1 mm thick, and a thin rubber membrane, one element O-1 mm thick, as shown. This mesh converged to within approximately 1% of other meshes with more radial and axial elements. The material properties listed by Dunn were used initially and are given in Table 1. Two key coefficients were not measured but were book values. These coefficients were the bulk modulus of the solid matrix and the permeability. These coefficients affect the calculations of the effective bulk modulus of the fluid and the drag coefficient respectively. Also the Young’s modulus of the solid was not given but could be determined from the plotted value at zero frequency. Dunn produced theoretical graphs for water saturated cylinders with radii of 9, 19 and 29 mm and with permeabilities of 1 x 10-9, 10x 10e9, 30x 10e9 and 100x 10e9 mm*. To calibrate the model the 19 mm radius cylinder with a permeability of 30 x 10T9mm2 was used. This calibration was used to verify the calculated values of the effective bulk modulus and to obtain the drag coefficient. Dunn used two models, one was quasistatic without inertia terms (i.e., zero mass), and the other included inertia. The initial comparison was carried out without inertia to limit the variables involved. To obtain exact agreement with Dunn’s curves the bulk modulus had to be reduced to 3.65 GPa from the calculated value of 13.26 GPa. The bulk modulus affects the peak value of the curve without significantly moving the curve horizontally. However, as will

402

M. F. YE0

AND

L. J. SCHMID

TABLE

Material

properties

used by Dunn

Property

KY G

4

KY P/

PS

k” 4 E

1

[ 141

Value

5.2 6.7

GPa GPat 2.2 GPa 35-OGPat 1000kg/m3 2650 kg/m3 0.001 kg/ms 1, 10,30,100 x 1O-9mm*

0.11 13.65 GPa$

t Book values. $ Scaled value from figure.

be seen, to model Spencer’s experimental work accurately the calculated value is required. This reduction is due to the fact that for some unknown reason Dunn’s graphs, although they have the same peak attenuation frequency, are significantly lower than Spencer’s experimental curves. The drag coefficient does not basically change the shape of the graph but moves the graph parallel to the horizontal axis. To obtain the correct frequency for peak attenuation a drag coefficient of 1.1 x lo9 N/m3 of composite per meter per second relative flow was required. With the model now calibrated the analysis was rerun with different radius cylinders. The results are shown in Figure 2 and compared with those of Dunn. As can be seen, almost exact agreement was obtained. The porosity was then varied and the equivalent drag coefficient calculated according to drag coefficient = calibrated drag coefficient x 30 x lo-‘/porosity. The authors of reference [18] found that fluid velocities varied linearly with porosity for porosities between 2% and 30%. Hence with our definition of the drag coefficient an inverse relationship will exist between the drag and porosity. The 19 mm radius cylinder was then re-analyzed with permeabilities of 1, 10, 100 x 10m9mm2 and the results compared with these of Dunn. As can be seen in Figure 3, excellent agreement was again obtained. With the overlay technique now proven once the coefficients have been calibrated, the effects of inertia can be introduced. With finite mass density included in the model, as will be shown, Spencer’s results can be matched to an acceptable degree of accuracy. It should be noted that the determination of complex Young’s modulus simply as the rate of stress to strain measured experimentally is theoretically correct only as frequency tends to zero. It may be a reasonable approximation in the stiffness controlled regime: i.e., for frequencies small compared to the resonant frequency. With mass effects included the ratio of stress to strain becomes purely imaginary at the resonant frequency. Re-running the calibrating example produced the results shown in Figure 4, which are consistent with the above theory. Spencer’s experiments obviously included the inertia effects. However, even with the same water-saturated specimen, completely different results could be obtained. He admitted that the second of the tests was carried out after the specimen was used for other saturants but had been thoroughly cleaned. This indicates that slight modifications to the

PROPERTIES

OF SATURATED

POROUS

403

MEDIA

0.12 0=19mm

3

0.09

-

:.

0.06

-

2 2 :=

O-03

-

\

: f

_--0

10’

102

104

10s

1 I5

Frequency (Hz) Figure 2. Comparison equal to 9 mm (-- . -.

of attenuation

and Young’s modulus with Dunn [14], for specimens with radius and 29 mm (---). Results from present analysis are shown by x

), 19 m (--)

0.12 k = 100 X lo-’

k= 1 X lQg mm* k=

lOO-’

ry;,/

mm*

%! 0094

\

x\x

;’

op

3

/'

0.06

/

-

\

X

x

/

\

/ "\,

003-

,/"

z

/ _/r 0

./+

,J

'.J\ "x \

'x, '1X

“<

x%x,x

/ .X *&/_:I 10’

k

! )4

/A

g 9 5

\

‘“\ “\

-+-&-___
10s (Hz)

104

1 10s

’

Figure 3. Comparison of attenuation and Young’s modulus with Dunn [ 141 for specimens with permeability k equal to 1 x 10d9 mm* (-), 10 x 10d9 mm* (---) and 100 x 10m9 mm2 (- . -. -. ). Results from the present analysis arc shown by x .

a

404

M. F. YE0

s

L. J. SCHMlD

0.12-

uT \ B .g 0.09 lu .-E 5

AND

-

0.06

-

0.03

-

: f

IC P

Frequency

(Hz)

Figure 4. Comparison of attenuation and Young’s modulus with Spencer [ 131 and Dunn [ 141, when inertia is included. Spencer’s experimental results are shown by 0. Dunns results (with no inertia) are shown by -. The key for the results from the present analysis for values of drag coefficient, C, and effective bulk modulus, K,, is as follows: ---, C, = 2.2 x lo9 NmW4 s, K, = 13.26 GPa; -.-. , C, = 1.1 x lo9 N mm4 s, K, = 13.26 GPa; .. . ., C, = 1.1 x lo9 N mm4 s, K, = 3.65 GPa.

parameters of one or both phases makes a significant difference to the result. However, Spencer indicated that the experiments were repeatable. Using the calculated values for the effective bulk modulus, instead of the calibrated value used to produce agreement with Dunn’s results, produces a graph that is similar to Spencer’s but shifted horizontally. Adjusting the drag coefficient by a factor of two, to 2.2 x lo9 N me4 s produces reasonable agreement, as can be seen from Figure 4. The original drag coefficient gave the same frequency for peak attenuation as Dunn, but his graph is offset vertically according to Spencer, who only estimated the frequency for peak attenuation. If Spencer’s estimate had been incorrect and is in fact the frequency indicated by our analysis, with the increased drag coefficient, then Dunn, when selecting his coefficients to obtain agreement with Spencer, could just as easily have obtained this frequency by modifying his coefficients and shifting all his curves horizontally. In this situation it would not have been necessary for us to adjust the drag coefficient after the initial calibration. Very similar results were obtained for the ethanol and N-decane saturants as were obtained by Spencer and by Dunn. Once the coefficients are calibrated, then both the experimental and theoretical results can be reproduced for these saturants to the same order of accuracy as was obtained above for the water-saturated sandstones. It is obvious that the drag and bulk modulus need to be determined accurately for a particular solid/fluid situation. Unfortunately, the authors do not have facilities for acoustic testing even if accurate determination of the coefficients were possible. Hence

PROPERTIES

POROUS

OF SATURATED

no attempt has been made to model an experiment been accurately determined.

405

MEDIA

in which these two parameters had

6. WAVE SPEEDS The properties of propagating waves in saturated porous media can also be obtained for an infinite prismatic section from the analysis of the axisymmetric mesh used above, but with axial length of 5 mm instead of 140 mm. The assembled dynamic matrix, equation (5), may be readily rearranged into a transfer matrix form. The eigensolutions of this transfer matrix provide propagating wave mode shapes, the eigenvectors, and propagation constants, the eigenvalues. Details of this procedure have been given in reference [16]. The results of the transfer function analysis for the 19 mm radius infinite cylinder are shown in Figure 5 for frequencies from 0 to 20 000 Hz. As can be seen, both the wave speed and attenuation rise and then decrease with increasing frequency. However, the fall-off in the wave speed is not significant and is relatively constant above 2000 Hz. For different radii infinite cylinders the wave speed is relatively unchanged but the attenuation changes significantly. For example at 10 000 Hz the attenuation and wave speed are as shown in Table 2. As the radii increases, for a particular drag coefficient, the saturant is less able the “escape” to the outside through the open pore boundary

1

3500

Wave

0

2

4

speed

I

1

I

6

8

10

Frequency

Figure 5. Effect on attenuation increased.

I

12

I

16

and wave speed for an infinitc4y

speed

1

18

2000 20

CkHz)

TABLE

Wave

I

14

long saturated

2

and attenuation for prism diameters

various

injinite

Speed

dB/m

11

2636

4.618

15 19 23 27

2652 2658 2658 2654

3.305 2.3157 1.635 1.182

Diameter

cylinder

as the frequency

is

406

M. F. YE0 AND

L. J. SCHMID

conditions. This therefore minimizes the relative movement of the fluid relative to the porous solid, and hence reduces the losses, thus the attenuation is reduced. In all cases it was found that the fluid moves in phase with the porous solid, with the fluid axial displacements being slightly smaller than those of the porous solid. For the frequency range considered the case of the fluid and porous solid displacements being out of phase was not observed. As the radius of the infinite cylinder increases to infinity the attenuation will obviously converge to some effectively constant value. However, to obtain the value a significant amount of computing with greatly refined finite element meshes would be required. This parametric study has not been attempted at this time. 7. CONCLUSIONS The overlay finite element technique, with a single drag coefficient being used to represent all the interphase effects, can accurately mode1 the physical and acoustic response of saturated media once the material constants for the composite have been accurately determined. Porosity variations can be modelled via a simple relationship with the drag coefficient. The effects of changes in pore size and shape on the drag coefficient have not been ascertained. Future research will therefore be necessary to derive the required relationship. REFERENCES 1. M. A. BIOT 1956 Journal ofthe Acoustical Society ofAmerica 28(2), 168-191. Theory of elastic waves in fluid saturated porous solids. 2. G. H. F. GARDNER 1962 Journal of the Acoustical Society of America 34(l), 36-40. Extension waves in fluid saturated porous cylinders. 3. J. G. BERRYMAN 1983 Journal of the Acoustical Society of America 34(l), 36-40. Extensional waves in fluid saturated cylinders. 4. W. F. MURPHY 1982 Journal of the Acoustical Society of America 71(6), 1458-1468. Effects of partial water saturation on attenuation in Massilon sandstone and Vycor porous glass. 5. T. G. JONES 1986 Geophysics 51(10), 1939-1953. Pore fluids and frequency dependent wave propagation. 6. B. R. TITTMAN, J. R. BULAN and M. ABEL-GAWAD 1983 in Physics and Chemistry of Porous Media (D. L. Johnson and P. N. Sen, editors) New York: American Institute of Physics. Conference Proceedings no. 107, Chapter 3, pp. 131-143. The role of viscous fluids in the attenuation and velocity of elastic waves in porous rocks. 7. R. J. O’CONNELL and B. BUDIANSKY 1977 Journal of Geophysical Research 82(36), 5719-5735. Visco-elastic properties of fluid saturated cracked soils. 8. W. F. MURPHY, K. W. WINKLER and R. L. KLEINBERG 1983 in Physics and Chemistry of Porous Media (D. L. Johnson and P. N. Sen, editors) New York: American Institute of Physics. Conference Proceedings no. 107, Chapter 3, pp. 176-190. Contact microphysics and viscous relaxation in sandstones. 9. M. N. TOKSOZ, D. H. JOHNSON and A. TIMUR Geophysics 44(4), 681-711. Attenuation of seismic waves in dry and saturated rocks. 10. T. J. PLONA 1980 Applied Physics Letters 36(4), 259-261. Observations of a second bulk compressional wave in a porous medium at ultrasonic frequencies. 11. 0. M. LOVERA 1987 Geophysics, 52(2), 174-178. Reflections at elastic medium and saturated porous medium. 12. J. E. WHITE 1986 Geophysics 51(3), 742-745. Biot-Gardner theory of extensional waves in porous rods. 13. J. W. SPENCER 1981 Journal of Geophysical Research 86(3), 1803-1812. Stress relaxations at low frequencies in fluid saturated rocks: attenuation and modulus dispersion. 14. K. J. DUNN 1986 Journal of the Acoustical Sociefy of America 79(6), 1709-1721. Acoustic attenuation in fluid saturated rocks at low frequencies. 15. Y. K. CHEUNG and M. F. YEO 1979 A Practical Znrroduction to Finite Element Analysis. London: Pitman.

PROPERTIES

OF SATURATED

POROUS

MEDIA

407

SCHMID 1989 Journal of Sound and Vibration 130(3), 439-445. Wave propagation in solid and fluid structures using finite element transfer matrices. 17. 0. C. ZIENKIEWICZ 1989 The Finite Element Method in Engineering Science. New York: McGraw Hill. 18. DE-hut HAN, A. NUR and D. MORGAN 1986 Geophysics 51(11), 2093-2107. Effects of porosity and clay content on wave velocities in sandstones.

16. M. F. YEO and L. J.