## Abstract

A mathematical model for mechanical and electrical dynamics in an anisotropic piezo-poroelastic (hereafter referred to as APP) medium is solved for three-dimensional propagation of harmonic plane waves. A system of modified Christoffel equations is derived to explain the existence and propagation of four waves in the medium of arbitrary anisotropy. This system is solved to calculate the phase velocities of four waves in an unbounded APP medium. Directional derivatives of phase velocity are derived analytically and are used to calculate the components of the ray velocity vector. A hypothetical numerical model is considered to compute the phase velocity for given (arbitrary) phase direction and then the ray velocity vector. Surfaces are plotted for the phase velocity and ray velocity of each wave in a saturated poroelastic medium in the absence/presence of piezoelectricity. The contributions of the piezoelectric activeness of the solid frame and pore-fluid to the phase and ray velocities are identified and analysed for each of the four waves in the medium.

## 1. Introduction

Piezoelectric materials show unique electromechanical coupling characteristics which produce electrical fields under the application of mechanical loads and elastic deformation on applying electrical loads. These materials find their utility as sensors and actuators in many applications involving signal transmission. Such materials are generally brittle in nature and efforts are often made (Taunaumang *et al.* 1994; Chan *et al.* 1999) to develop composite piezoelectric materials to mitigate their brittleness and to strengthen the coupling between their electrical and mechanical responses. For example, dispersion of a fibrous or a particulate phase in a polymer matrix results in a flexible as well as an effective piezoelectric structure. Due to the brittle nature of piezoelectric ceramics and possible defects of impurity, cavities, micro-voids, layer separations, inclusions and microcracks, failures of devices take place under mechanical or electrical loading. This limitation is overcome by reducing the material density through the addition of controlled porosity. The resulting piezoelectric porous materials are widely used for applications such as low-frequency, miniature hydrophones, accelerometers, vibratory sensors and contact microphones (Arai *et al.* 1991; Hayashi *et al.* 1991). Due to lower acoustic impedance, these materials can be used to improve mismatch of acoustic impedance at the interfaces of medical ultrasonic imaging devices or underwater sonar detectors (Mizumura *et al.* 1991).

In the last three decades, analytical and numerical models have been developed to predict the electromechanical behaviour of piezoelectric composites. Based on the connectivity of their constituent phases, Newnham *et al.* (1978) proposed a classification of piezoelectric composites. Hikita *et al.* (1983) confirmed that piezoelectric materials with large porosity levels display a high voltage output *g*-constant. Banno (1987) studied the effects of size, geometry and distribution of pores in piezocomposites. Bisegna & Luciano (1996, 1997) generated variational principles for electro-elastic characteristics and determined the bounds for the overall properties of piezoelectric materials. Gomez & Espinosa (1996) derived a model to predict the dielectric properties of porous piezoelectric materials and the dielectric coupling generated by the interstitial phase filling of the pores in a porous matrix. Gomez & Freijo (1996) and Gomez & Espinosa (1997) identified solutions for the effective electromechanical constants of composite piezoelectric materials. Li & Wang (2005) studied the wave localization in disorderly layered piezoelectric composite structures. In a recent study, Kar-Gupta & Venkatesh (2006) developed an analytical model to capture the complete electromechanical response of a piezoelectric composite system in which both matrix and fibre phases are considered to be elastically anisotropic and piezoelectrically active.

Despite numerous experimental studies on the dynamic behaviour of piezoelectric composites, theoretical studies of wave propagation in these materials are few. It was in the early 1990s that researchers in the field of applied mechanics realized the importance of wave propagation in analysing the electromechanical response of elastic materials. Nayfeh & Chien (1992) studied the interaction of ultrasonic waves with piezoelectric substrates and derived expressions for reflection and transmission coefficients. Zinchuk & Podlipenets (2000, 2001) discussed the dispersion relations for the propagation of shear waves and Rayleigh-type surface waves in piezoelectric layered structures. Recently Wang & Yuan (2007) investigated, theoretically and experimentally, the propagation characteristics of Lamb waves in composites with emphasis on group velocity and characteristic wave curves. Green’s functions for piezoelectric composites were obtained by Pan & Han (2004, 2005). The evaluation of Green’s function for an elastic medium is a fundamental problem in the field wave propagation. A medium’s response to an arbitrary point force is obtained through its convolution with Green’s tensor for the medium. Hence, the availability of Green’s tensor enables one to calculate the numerical solutions for a wide selection of media and sources (Ben-Menahem & Sena 1990). In an anisotropic medium, the ray velocities of waves in the medium are required to calculate Green’s tensor. Moreover, the geometrical properties of slowness (or phase velocity) surfaces are used in computing the displacement fields and radiation patterns from point sources in anisotropic media (Gajewski 1993). More recently, Vashishth & Gupta (2009*a*) studied the wave propagation in transversely isotropic porous piezoelectric materials and analysed the piezoelectric effect on slowness surfaces. Biot’s (1955, 1956) theory of poroelasticity was used to derive an analytical model that governs the electrical and mechanical dynamics in a porous piezoelectric medium. By considering transverse isotropy, it may be possible to obtain simplified expressions but the likely candidates of anisotropic materials may have the lower symmetries. Then the arbitrary anisotropy provides the liberty to choose any desired combinations of elastic as well as dielectric properties.

The work presented herein studies the slowness surfaces as well as wave surfaces for the propagation of harmonic plane waves in anisotropic piezo-poroelastic (APP) materials. Biot’s theory of poroelasticity is used to derive the analytical model that governs the electrical and mechanical dynamics in a piezo-poroelastic medium of arbitrary (general) anisotropy. This mathematical model is solved for the propagation of plane harmonic waves. Phase velocities of the four waves are determined as the roots of a quartic polynomial. The phase velocity and the corresponding phase direction of a wave are used to define its ray velocity (wave) surface. Ray velocity surfaces and phase velocity surfaces of the four waves are computed for a particular numerical model of an APP material (barium titanate). The piezoelectric contributions of the solid and fluid phases of piezo-poroelastic medium in changing the phase and ray velocities are exhibited as functions of phase direction.

## 2. Anisotropic piezo-poroelasticity

The influence of nonlinearity, dissipation and dispersion on wave propagation in poroelastic material can be estimated from the porosity, the character pore size and the frequency (or wavelength) of the wave. A nonlinear mathematical model is required to represent the dynamics in a poroelastic medium when the length of the wave is assumed to be much greater than the distance between the pores, and the distance between the pores, in turn, is much greater than the radius of the pores. Poroelastic dissipation is induced by the motion of pore-fluid relative to a solid frame. It is governed by Darcy’s law with dependence on hydraulic permeability and the porosity. The field equations of Darcy’s law become nonlinear when a variable permeability is taken into account. Pore-size is used to describe the dependence of viscous resistance to fluid-flow on the size and shape of the pores. A study of Chapman & Liu (2007) suggests that the nonlinear wave propagation results from the existence of a second fluid in the porous solid. The nonlinear effects are assumed to be due to interphase interactions and, in the case of a single saturating fluid, linear wave propagation is considered. The partial gas saturation is expected to have a strong impact on the attenuation and dispersion of waves. Chapman *et al.* (2002) considered the pore space as consisting of a collection of thin ellipsoidal cracks and spherical pores which are connected by some effective permeability. As a seismic wave passes, pressure gradients are created between the spherical pores and the thinner cracks. Fluid flow takes place to relieve these pressure gradients, giving rise to dispersion and attenuation.

In the present study, the poroelastic medium is assumed to consist of a porous solid matrix with its interconnected pores filled with a single compressible fluid. The deformations are assumed small enough, which guarantee the linearity of the mechanical processes. The fluid is assumed to be without viscosity, which also represents the limiting behaviour of wave propagation at very high frequencies. It also represents a case of no relative motion of fluid in the pores to the solid frame such that dissipation due to fluid friction disappears. The wave propagation is governed by Biot’s (1956) theory of poroelasticity. According to this theory, the fluid-phase is supposed to be continuous through connected pores. Thermomechanical phenomena are neglected. The macroscopic system is described for a volume element, which is assumed small relative to the wavelength of the elastic waves and the size of pores is assumed small when compared with the size of the element. Following Biot (1956), the governing equations for a fluid-saturated porous media in the absence of body forces and dissipation are
2.1
In these equations, *u*_{i} and *U*_{i} are the components of displacements of the solid and fluid particles, respectively. The dot notation is used to represent differentiation with respect to time. Indices can take the values 1, 2 and 3. The summation convention is valid for repeated indices. The comma (,) before an index represents partial space differentiation. *ρ*_{11},*ρ*_{12} and *ρ*_{22} are the dynamical constants related to the porosity of solid (*f*), tortuosity of pores (), densities of solid particles and interstitial fluid (*ρ*_{s},*ρ*_{f}), through the relations: . The absence of dissipation in equation (2.1) implies that the skeletal frame is elastic and the pore-fluid is non-viscous or there is no fluid motion in the pores relative to the frame.

For an APP material, the constitutive equations for stresses in solid frame (i.e. *σ*_{ij}) and pore-fluid (i.e. *σ*) are (Vashishth & Gupta 2009*b*)
2.2

The gradients of electric potentials *ψ* and *χ* define the electric fields in the solid and fluid phases, respectively. The coefficients *c*_{ijkl},*m*_{ij} and *R* are the material constants of a linear porous material. The piezoelectric activeness of the solid phase and fluid phase of poroelastic medium affect the stresses in the solid frame through the third rank tensors of coupling coefficients *η*_{ijk} and *ζ*_{ijk} (symmetric in indices *j* and *k*), respectively. The coefficients *a*_{j} and *b*_{j} control the equivalent effect on isotropic stresses in pore-fluid.

The electric field displacements in the solid frame (*D*_{i}) and interstitial fluid () are given by (Vashishth & Gupta 2009*b*)
2.3
where second rank tensors *α*_{ik},*β*_{ik} and *γ*_{ik} define the dielectric coefficients for the two phases of the poroelastic medium. The electric field displacements in these phases of the poroelastic medium are governed by the Maxwell equations, i.e.
2.4

To seek a solution for the propagation of harmonic plane waves in an APP medium, we write
2.5
where *ω* is the angular frequency and (*p*_{1},*p*_{2},*p*_{3}) the slowness vector. For *p*_{k}=*n*_{k}/*v*, the slowness vector defines *v* as the phase velocity of a wave propagating along the phase direction (*n*_{1},*n*_{2},*n*_{3}), a unit vector. Substituting equation (2.5) in equations (2.1) and (2.4), through equations (2.2) and (2.3), we get a system of eight homogeneous equations in *S*_{1},*S*_{2},*S*_{3},*F*_{1},*F*_{2},*F*_{3},*G*,*H*. Algebraic manipulation of this system of equations yields the relations
2.6
and a system of six linear homogeneous equations given by
2.7
The matrices in equation (2.7) are expressed as
2.8
where **I** is an identity matrix of order three and row matrix **N**=(*n*_{1},*n*_{2},*n*_{3}). Various notations used in the above equations are defined as follows.
With the replacement of matrices **A**^{e}, **B**^{e}, **C**^{e}, **D**^{e} in equation (2.8) with a null matrix, the system (2.7) defines the propagation of plane harmonic waves in an anisotropic poroelastic medium (Sharma 2004). This implies that the piezoelectric effect on the dynamics of solid and fluid particles in an APP medium is represented through these four matrices only. For the special case of a poroelastic medium saturated with a piezoelectric fluid, these matrices reduce to
2.9
When the interstitial pore-fluid in a piezoelectric porous frame does not support piezoelectric effects, we have
2.10
The system of equations (2.7) is solved further to get a relation between displacements of solid and fluid particles, i.e.
2.11
and a system that defines the modified Christoffel equations for the wave propagation in an APP medium. The inverse of matrix **D** is expressed as
2.12
where **X**=**I**+**QN**^{T}**N**, **Y**=*r*_{22}**X**^{−1}, ** Φ**=

**adj**(

**Y**), ,

*d*

_{1}=−

**N**

*Φ***N**

^{T}. Non-symmetric matrix

**can be defined by the cyclic permutation of indices from two of its elements given by The modified Christoffel equations for the propagation of plane harmonic waves in an APP medium are given by, 2.13 where**

*Ψ*

*Δ*_{j}, (

*j*=1,2,3,4), the square matrices of order 3, are defined as follows:

### (a) Dissipation in porous medium

The expressions derived above are only valid for non-dissipative poroelastic solids. This implies that energy loss is not considered, i.e. the skeletal frame is elastic and the pore-fluid is non-viscous or there is no fluid motion relative to the solid frame. In the presence of viscosity in the pore-fluid, the scalar dynamical constants (*ρ*_{11},*ρ*_{12},*ρ*_{22}) change into second-order tensors as follows.
2.14
where tensor {*α*} represents the tortuosity of pores for viscous pore-fluid flow. In terms of tortuosity () for a non-viscous fluid, . Here **d** is a dissipation matrix, given by
2.15
where *ω*_{c} is the characteristic frequency (Biot 1956), *μ* the dynamic viscosity of pore-fluid, *Λ* the viscous characteristic length (Albert 1993) and tensor {*χ*_{a}} represents hydraulic (permeability) anisotropy. Then, the coefficients of the quartic equation (3.1) will be complex and yield complex velocities *v*_{j},(*j*=1,2,3,4), from the roots of this equation. These results apply to wave propagation as long as the wavelengths of the elastic waves are much larger than the linear dimensions of the pores and the grains in the medium. It is also assumed that the properties of the pore-fluid are unaffected by its proximity to the walls of the solid.

Another representation of viscous loss in terms of permeability, the pore size and the fluid–solid inertial coupling is also used (Sharma 2005). Wave propagation is discussed for two different frequency regimes. The low-frequency regime remains valid only when the flow in the pores is of Poiseuille type. In the high-frequency regime, only asymptotic cases are considered for very small or very large values of pore-fluid viscosity to pore-size ratio.

## 3. Phase velocity

For the non-trivial solution of the system of equations (2.13), the determinant of matrix *W* must vanish. This gives a quartic equation in *h* (=*ρv*^{2}/*R*) given by
3.1
In the absence of dissipation in the APP medium, the real coefficients *c*_{j}, (*j*=0,1,2,3,4) are derived (Sharma 2005) from matrices *Δ*_{j} in equation (2.13) and scalars *d*_{0} and *d*_{1} in equation (2.12). Only positive values of these coefficients ensure that all the four roots of this quartic equation are positive. Let *h*_{j}, (*j*=1,2,3,4), denote the four (positive) roots of this equation in a descending order. Then the phase velocities of the four waves, given by , will vary with the phase direction (*n*_{1},*n*_{2},*n*_{3}). Keeping in mind the propagation of compressional and shear waves in an isotropic poroelastic medium, the waves represented by *j*=1,2,3 and 4, are called the *qP*1, *qS*1, *qS*2 and *qP*2 waves, respectively. The prefix ‘*q*’ implies that these waves are either quasi-longitudinal (*qP*) or quasi-transverse (*qS*) because the polarizations of these four waves may not be along or normal to the phase direction.

## 4. Ray velocity vector

The phase velocities calculated in the previous section and the corresponding phase directions define the phase velocity surfaces of the four waves in the APP media. It is the plane of constant phase that propagates with phase velocity along the corresponding phase direction. In an anisotropic medium, there are different planes of constant phase for propagation along different phase directions. Let the unit vector in the Cartesian co-ordinate system represent the general direction in a local spherical co-ordinate system (*r*,*θ*,*ϕ*). So, corresponding to each wave in the medium, we have a two-parameter (i.e. *θ*,*ϕ*) family of constant-phase planes and the envelope of this family defines a ray velocity (or wave) surface. The energy (or envelope) of a wave in this medium travels with ray velocity but along a ray making an angle to the propagation direction. Let *v*(*θ*,*ϕ*) be the phase velocity of a wave along the general direction. Then, the ray direction at a point on the wave surface is determined from the components of the ray velocity vector at the point. Following Ben-Menahem & Sena (1990), these components *w*_{j}, (*j*=1,2,3) are expressed as
4.1
The magnitude of the ray velocity is defined as
4.2
and the ray direction (*θ*_{g}*ϕ*_{g}) is given by
4.3
The partial derivatives of *v*, for each of the four waves, derived from the quartic equation (3.1), are given by
4.4
The derivatives of the coefficients (i.e. , *k*=0,1,2,3,4) are derived analytically from the relations defined in §3.

## 5. Numerical results and discussion

The analytical expressions derived in the previous sections represent a mathematical model for propagation of elastic waves in an APP medium. These expressions constitute a procedure to calculate phase velocities for the propagation of four waves along a given (arbitrary) phase direction in the medium. The phase velocity of each wave along with phase direction is used to calculate its ray velocity and ray direction. The variations of phase velocities and wave velocities with phase or ray direction can only be analysed through numerical examples. However, a numerical example with real-time data may yield results of both qualitative and quantitative significance. But, due to the non-availability of real data, a hypothetical numerical model is chosen to explain the variations of velocities. Numerical results from such a particular model may not carry any quantitative significance and hence cannot be generalized. However, a qualitative analysis of these results may certainly be useful in understanding the wave-induced piezoelectric character of composite materials.

In this study, barium titanate is chosen as the anisotropic porous material with 20 per cent pore-space of tortuosity () 6 occupied by a fluid of density 50 kg m^{−3}. This symbolic fluid-density is chosen to lie between 1.2 kg m^{−3}, the density of air, and 103 kg m^{−3}, the partial density of CO_{2} (Garg & Nayfeh 1986). Numerical values of various constants and parameters involved in the derivations of previous sections are chosen as follows. The elastic matrix (in GPa) for the poroelastic medium is given by
The symmetric matrix *c*_{IJ} of order six represents the elastic tensor *c*_{ijkl} in two-suffix notation. The *ρ*_{11}, *ρ*_{12} and *ρ*_{22} are calculated for *ρ*_{s}=5700 kg m^{−3}, *ρ*_{f}=50 kg m^{−3} and *f*=0.2. The values (in GPa) assumed for parameters of elastic coupling are {*m*_{11},*m*_{22},*m*_{33},*m*_{12},*m*_{13},*m*_{23}}={2.5, 0.06, 0.05; 0.06, 2, 0.07; 0.05, 0.07, 1.5}. Positive definite symmetric matrices *α*_{ij}=0.1*γ*_{ij}, *β*_{ij}=0.5*γ*_{ij}, and *γ*_{ij}={1.3,0.03,0.04;0.03,0.9,0.02;0.04,0.02,1.5}, define dielectric permittivities in the unit of 10^{−8} *C*^{2} *N*^{−1}m^{−2}. The arrays *η*_{1jk}={0.5, 1.3, 11.3; 1.3, −1.5, 1; 11.3, 1, 0.5}; *η*_{2jk}={0.6,0.4,1.5;0.4,0.5,1.4;1.5,1.4,−1.5}; *η*_{3jk}={−4.1, 0.5, 1.5; 0.5, −1, −1; 1.5, −1, 17.4}, define third-order tensors of piezoelectric moduli *η*_{ijk} in the unit of *C* *m*^{−2} and then *ζ*_{ijk}=0.1*η*_{ijk}. Values (in *C* *m*^{−2}) for first order piezoelectric tensors are, arbitrarily, chosen as *a*_{i}=(0.002,−0.002,0.2) and *b*_{i}=(−0.001,0.001,−0.1).

Using the above numerical values, the phase velocity and ray velocity components are computed for each of the four (*qP*1, *qS*1, *qS*2, *qP*2) waves along the given (arbitrary) phase directions. The variations of the velocities with direction in three-dimensional space are displayed in the form of velocity surfaces. The locus of point { in the three-dimensional Cartesian coordinate system defines the phase velocity surface of a wave propagating with phase velocity *v* along the (*θ*, *ϕ*) direction. Similarly, the points (*w*_{1},*w*_{2},*w*_{3}) defined in equation (4.1) constitute the ray velocity surface to this wave. The (phase and ray) velocity surfaces of the four waves in poroelastic medium are exhibited in figure 1. The corresponding surfaces in the APP medium are plotted in figure 2. The comparisons of corresponding plots in these figures identify the effect of piezoelectric activeness on the velocities of four waves in the medium. For the chosen numerical model, the anisotropy in the velocities of the qP1 and qP2 waves is not very significant. The effect of anisotropy is observed much on the ray velocity, particularly, for the propagation of the qS1 and qS2 waves.

There may not be a purpose in discussing the velocity surfaces, in detail, for the considered hypothetical model. Moreover, the effect of piezoelectric response of the APP medium to these surfaces is not quite clear as seen from the comparisons of figures 1 and 2. The APP medium being a two-phase medium, the separation of piezoelectric effect due to each of its constituents may be interesting. To explore the piezoelectric response of the solid frame and pore-fluid, the variations of the phase velocities and ray velocities with phase direction are presented in figures 3–6, for the qP1, qS1, qS2 and qP2 waves, respectively. In these figures, the plots in (*a*) exhibit the variations of phase velocity and ray velocity of the qP1, qS1, qS2 and qP2 waves in a poroelastic medium without any piezoelectric presence. The effect of the piezoelectric character of the solid frame only on velocities is presented through the plots in (*b*) of these figures whereas the effect of piezoelectricity from the pore-fluid only is shown in the plots in (*c*). The plots in (*d*) show the effect of piezoelectric activeness of both constituents of the poroelastic medium on the phase and ray velocities.

It is noted that for this chosen model of an APP medium the major piezoelectric change to the velocities of the three faster waves comes from the piezoelectric activeness of the solid frame. Piezoelectric fluid in the pores has a negligible effect on the phase as well as the ray velocity of any of the four waves. However, the piezoelectric response of the solid frame on the slowest (qP2) wave is smaller than even that of the pore-fluid. It was checked that the piezoelectric effect of pore-fluid does not improve even when the pores are filled with a high-density fluid (say, liquid) with the same piezoelectric properties. However, the presence of high-density fluid in the pores may decrease the velocity of qS2 wave to half. In general, the effect of the piezoelectric activeness of the solid frame or pore-fluid or both may be a little more on ray velocities than on phase velocities.

## 6. Concluding remarks

This study presents the mathematical description of piezoelectric wave propagation in saturated poroelastic materials. This enables one to understand the effects of piezoelectric properties on wave propagation in realistic media. The arbitrary anisotropy provides the liberty to choose any combination of elastic or piezoelectric constants to represent the anisotropic symmetries existing in crystals. The analytical derivations in this work constitute a procedure to draw the phase velocity and ray velocity (or, wave) surfaces for any of the four waves in a poroelastic medium with the option of choosing the piezoelectric activeness of the pore-fluid or solid frame or both. Ray velocity vectors are tangent to the rays constructed from a point source. The loci of the end points of these vectors from the source form the wave surface. Hence, a wave surface shows the actual form of the wave front of the wave generated by a point source.

The analytical expressions derived in this work provide the essentials to derive the characteristic equation for surface waves and to study the reflection/refraction phenomena in APP media. The derivations for ray velocity surfaces provide the essential input to study the scattering from a small inhomogeneity in anisotropic media. Solving equations (2.11)–(2.13) for polarizations of solid and fluid particles in porous medium, the relation (2.6) can be used to calculate the wave-induced piezoelectric fields in the medium. Piezoelectric poroelastic materials certainly have some advantages over dense monolithic piezoceramics. The formulation based on Biot’s theory may not be as convenient as the micromechanics theory of composite materials but with the presence of interaction parameters it represents a more general mathematical model in the present study. Another difficulty is to compute the specific values for these parameters. Then, it is felt that while studying the effect of a parameter, sometimes a bracketed range may provide more useful interpretation when compared with a specific value (keeping in mind the errors associated with its measurement). However, with the availability of real data, this work can be used for simulation studies aimed at improving the designs of devices used in signal processing and communication. This study may also be helpful to structural engineers working on exploration of the internal structure of composite materials through travel times or velocity anisotropy. Researchers in the field of non-destructive evaluation may use the expressions derived in this work to improve their models when interpreting complex data from real-time experiments.

## Footnotes

- Received October 9, 2009.
- Accepted January 8, 2010.

- © 2010 The Royal Society