Tortuosity Perturbations Induced by Defects in Porous Media

Abstract

In this chapter, we describe the effects of defects in a homogeneous saturated porous medium. Defects are modelized by inclusions which disturb the motion of the viscous fluid flowing in the pore space of the medium. The seepage rate of the fluid in the host medium and in the inclusion is given by the Darcy’s law. Disturbances thus produced modify the shape of the stream lines from which we establish the tortuosity induced by the defects and its implications on the acoustic waves propagation in saturated porous media.

1. Introduction

Among the essential physical parameters to describe the microstructure of porous media, tortuosity is one of the most important parameters. For a review, we can refer to the paper of Ghanbarian et al. [1].

Tortuosity was introduced as a correction to the permeability of Kozeny’s model [2] of porous media defined by the Darcy’s law relating the fluidic characteristics and pore space of the medium [3]:

$$\mathbf{v}=-\frac{k}{\eta }\mathbf{\text{∇}}p,$$ E1

where $\mathbf{v}$ is the seepage rate of the fluid, $\eta$ the viscosity coefficient of the fluid, $\mathbf{\text{∇}}p$ is the pressure gradient applied to the medium, and $k$ is its permeability. The Kozeny’s model was developed in the framework of straight and parallel streamlines in porous media. Carman has generalized it to neither straight nor parallel streamlines by introducing the hydraulic tortuosity $\tau$ defined by:

$$\tau =\frac{<\lambda >}{L}.$$ E2

When a fluid flows through a porous medium from point $A$ to point $B$ distant from $L$ (Euclidean distance) (Figure 1), it follows different paths whose mean length is $<\lambda >$ , where $\lambda$ is the length of the different paths connecting these two points. In isotropic media, the tortuosity is a scalar number greater than unit ( $<\lambda >\ge L$ ), whereas for the low porous media, its values may be greater than 2; they range from 1 to 2 for high porosity media such as fibrous materials and some plastic foams.

Figure 1.

Some tortuous paths through a porous medium.

The lengthening of the field paths in porous media due to tortuosity does not only occur in the flow of fluids in porous media, but is a more general result. So we meet this concept in processes such as transport phenomena, particles diffusion, electric conductivity, or wave propagation in fluid saturated porous media. Researchers have thus developed many theoretical models adapted to their concerns to introduce the tortuosity, leading to unrelated definitions of this concept. For instance, Saomoto and Katgiri [4] presented numerical simulations to compare hydraulic and electrical tortuosities. Thus, using numerical models of fluid flow and electric conduction in same media, i.e., with the same local solid phase arrangements, the authors show that while electrical tortuosity remains close to the unit whatever the porosity and the shape of the grains, the stream lines of hydraulic flow are much more concentrated in some parts of the medium, leading to a much greater tortuosity.

This example shows that although the physical meaning of this parameter is obvious, in practice, it is not consistent and its treatment is often misleading. The conclusion that emerges from these observations is that tortuosity should not be viewed as an intrinsic parameter of the environment in which the transport process develops, but rather as a property of this process. This partly explains why there are different definitions of tortuosity, each with its own interpretation.

In acoustics of porous media, tortuosity has been introduced to take into account the frequency dependence of viscous and thermal interactions of fluid motion with the walls of pores. In [5], Johnson uses it to renormalize the fluid density ${\rho }_f$ . When the viscous skin depth is much larger than the characteristic dimensions of the pore, Lafarge et al. [6] have shown that the density of the fluid is equal to ${\rho }_f{\tau }_0$ , where ${\tau }_0$ is the static tortuosity for a constant flow ( $\omega =0$ ) defined by:

$${\tau }_0=\frac{<{\mathbf{v}}^2>}{<\mathbf{v}{>}^2},$$ E3

where $<.>$ denotes averaging over the pore fluid volume $V_f$ . Thereafter, in this chapter, we adopt this definition of tortuosity.

Through the definition (3), we see that the tortuosity is given as soon as the permeability of the porous medium is known in each of its points. As it is well known, many factors can affect the fluid flow in porous media, including pore shape, distribution of their radii, and Reynolds number to name a few. It follows that the presence of defects in an initially homogeneous medium (for instance, a local change of an intrinsic parameter) can be an important disturbance of the fluid motion, the result being a modification of the shape of the streamlines.

Taking into account the presence of defects that change the permeability of the porous medium leads to the notion of effective permeability ( $k_{\mathit{eff}}$ ). In general, the $k_{\mathit{eff}}$ value is not unique but depends on the chosen model for the homogenization of the porous medium. The homogenization process only makes sense for lower scales than the spatial variations of incident excitation, which therefore justifies that mobility is calculated for a low-frequency filtration rate (quasi-static regime). These considerations lead us to be interested only in the instantaneous individual response of defects to external solicitations. Since in our case only media with low levels of defect are considered, it is legitimate to ignore their mutual interactions.

The present chapter is organized as follows. Section 2 describes the mathematical model of the defects and gives the solution of the fluid flow in the presence of homogeneous and layered spherical and ellipsoidal defects. Then, the results are generalized to anisotropic defects. Finally, the hydraulic polarizability is introduced. Section 3 is relative to tortuosity. The expression of effective mobility is given for some particular defects. The induced tortuosity is deduced from the previous results and its effects on the wave propagation are given.

2. Defect model

In this chapter, what is called defect is a local change of permeability $k$ . Such a change is due, for instance, to variations in porosity in the microstructure of the medium. In this chapter, a defect is modelized as a porous inclusion $\mathrm{\Omega }$ characterized by its shape and own parameters: intrinsic permeability $k^{\left(i\right)}$ and porosity ${\varphi }^{\left(i\right)}$ . Intrinsic permeability is expressed in darcy: $1D=0.97\times {10}^{-12}{\mathrm{m}}^2$ . The porous media we are interested in have permeabilities of the order of $10D$ . Moreover, it is supposed that the fluid saturating the inclusion $\mathrm{\Omega }$ is the same (with viscosity coefficient $\eta$ ) as that flowing in the porous medium. Thereafter the mobility of the fluid defined by $\kappa =k/\eta$ is used. This notion combines one property of the porous medium (permeability) with one property of the fluid (viscosity). The inclusion is embedded in a porous medium with porosity ${\varphi }^{\left(o\right)}$ and permeability $k^{\left(o\right)}$ . The saturating fluid is subject to action of a uniform pressure gradient $\mathbf{\text{∇}}{p}^{\left(0\right)}$ . In the sequel, we use indifferently the words defect or inclusion.

2.1 Mathematical formulation

When the fluid flows through the porous medium, its motion is perturbed by the defects in the microstructure of the medium. Within the porous medium, the velocity $\mathbf{v}$ and the pressure gradient $\mathbf{\text{∇}}p$ are related by the Darcy’s law:

$${\mathbf{v}}^{\left(m\right)}\left(x\right)=\frac{k^{\left(m\right)}\left(x\right)}{\eta }\mathbf{\text{∇}}{p}^{\left(m\right)}\left(x\right),$$ E4

where $m=i$ if $x\in \mathrm{\Omega }$ and $m=o$ if $x$ is in the host medium ( $x\notin \mathrm{\Omega }$ ). These equations are subject to the following boundary conditions on $\mathrm{\partial \Omega }$ :

• continuity of fluid flow

$${\varphi }^{\left(o\right)}{\mathbf{v}}^{\left(o\right)}.{\mathbf{n}}^{\left(o\right)}={\varphi }^{\left(i\right)}{\mathbf{v}}^{\left(i\right)}.{\mathbf{n}}^{\left(i\right)},$$ E5

• continuity of normal stress component

$${\tau }_{\mathit{ij}}^{\left(o\right)}{n}_j^{\left(o\right)}={\tau }_{\mathit{ij}}^{\left(i\right)}{n}_j^{\left(i\right)}$$ E6

where ${\mathbf{n}}^{\left(o\right)}$ and ${\mathbf{n}}^{\left(i\right)}$ are unit vectors perpendicular to the interface. The hypothesis that Darcy’s law governs the dynamics of the flow of fluid in a porous excludes inclusions filled with fluid. Indeed, within such inclusions, the movement of the fluid is governed by the Navier-Stokes equations, which are impossible to reconcile with the law of Darcy in the porous medium with the available boundary conditions [7]. Figure 2 represents an oriented inclusion in a fluid in motion.

Figure 2.

Oriented inclusion in a porous medium.

A porous medium of infinite extension is considered in which a viscous fluid flows at a constant uniform velocity $U_{\infty }$ under the action of the pressure gradient along the $\mathit{Ox}$ axis. We want to determine the local changes of the fluid velocity when defects are present in the medium.

In the following, we give the solutions of Eqs. (4)–(6) for some particular defects in such situation. Analytical solutions are possible for homogeneous spherical and ellipsoidal inclusions. We show that the most important characteristic of these inclusions is their hydraulic dipole moment. For layered defects, i.e., when their permeability is a piecewise constant function, we give a matrix-based method to get their dipolar moment.

2.2 Isotropic homogeneous defect

In this section, we assume that the background and the defect (embedded inclusion) are homogeneous each with its own parameters: porosity ${\varphi }^{\left(o\right)}$ and ${\varphi }^{\left(i\right)}$ and permeability $k^{\left(o\right)}$ and $k^{\left(i\right)}$ which have constant values. A static incident pressure with a constant gradient along $\mathit{Ox}$ axis is applied to this system. What we seek is the pressure perturbation produced by the defect acting as a scatter and the expressions of the resulting seepage rate of the fluid inside and outside the porous inclusion.

2.2.1 Spherical defect

The simplest type of inclusion is the homogeneous spherical one, and we consider a porous sphere of radius $r=a$ , centered at the origin of axes with a constant permeability $k^{\left(o\right)}$ . Using the spherical coordinates ( $r$ , $\theta$ , $\phi$ ), Eqs. (4)–(6) become

$${\mathbf{v}}^{\left(m\right)}=\frac{k^{\left(m\right)}}{\eta }\mathbf{\text{∇}}{p}^{\left(m\right)}\text{where}\left\{\begin{array}{l}m=i\text{if}r<a\\ m=o\text{if}r>a\end{array}\right.$$ E7

with the boundary conditions at $r=a$ :

$${\varphi }^{\left(i\right)}{v_r}^{\left(i\right)}\left(r=a\right)={\varphi }^{\left(i\right)}{v_r}^{\left(i\right)}\left(r=a\right),{\tau }_{\mathit{rr}}^{\left(o\right)}{\left|{}_{r=a}={\tau }_{\mathit{rr}}^{\left(i\right)}\right|}_{r=a}\forall \theta ,$$ E8

where

$${\tau }_{\mathit{rr}}=\varphi \left(-p+2\eta \frac{\partial v_r}{\partial r}\right),$$ E9

$v_r$ being the radial component of the velocity. For an incompressible fluid, Eq. (4) becomes

$$\mathrm{\Delta }{p}^{\left(m\right)}=0m=i,o.$$ E10

In spherical coordinates, this equation is written as:

$$\mathrm{\Delta }p=\frac{\partial }{\partial r}\left(r^2\frac{\partial p}{\partial r}\right)+\frac{1}{\text{sin}\theta }\frac{\partial p}{\partial \theta }\left(\text{sin}\theta \frac{\partial p}{\partial \theta }\right)+\frac{1}{{\text{sin}}^2\theta }\frac{{\partial }^{2}p}{{\partial }^2\phi }.$$ E11

The spherical symmetry of the problem (Figure 3) implies that the solution does not depend on $\phi$ . It follows that its solution is:

$$p^{\left(m\right)}\left(r,\theta \right)=\sum _l\left(A_l^{\left(m\right)}{r}^l+B_l^{\left(m\right)}{r}^{-\left(l+1\right)}\right)P_l\left(\text{cos}\theta \right),$$ E12

$P_l\left(x\right)$ being the Legendre polynomial of degree $l$ . The coefficients $A_l^{\left(m\right)}$ and $B_l^{\left(m\right)}$ are related by the boundary conditions (8), those at $r=0$ and when $r\to \infty$ . Inside the inclusion, the pressure must be finite at $r=0$ . This condition leads to $B_l^{\left(i\right)}=0$ for all $l$ . So $p^{\left(i\right)}$ becomes:

Figure 3.

Spherical inclusion in a fluid flow.

$$p^{\left(i\right)}\left(r,\theta \right)=\sum _l{A}_l^{\left(i\right)}{r}^l{P}_l\left(\text{cos}\theta \right).$$ E13

Outside, far from the inclusion, the pressure is:

$$p^{\left(o\right)}\sim -\frac{1}{{\kappa }^{\left(o\right)}}{U}_{\infty }\mathit{\text{rcosθ}}.$$ E14

It follows that the condition

$$p^{\left(o\right)}\left(r=a,\theta \right)=p^{\left(i\right)}\left(r=a,\theta \right)\forall \theta ,$$ E15

which implies that only the term with $l=1$ remains in the sum. Thus, we get:

$$A_1^{\left(o\right)}=-\frac{1}{{\kappa }^{\left(o\right)}}{U}_{\infty }.$$ E16

The expressions of the pressure are then:

$$p^{\left(i\right)}=A_1^{\left(i\right)}\mathit{\text{rcosθ}}\text{if}r<a,$$ E17

$$p^{\left(o\right)}=-\frac{U_{\infty }}{{\kappa }^{\left(o\right)}}\mathit{\text{rcosθ}}+B_1^{\left(o\right)}{r}^{-2}\mathit{\text{cosθ}}\text{if}r>a,$$ E18

$A_1^{\left(i\right)}$ and $B_1^{\left(o\right)}$ been given by the conditions (8). Finally, these expressions are:

$$p^{\left(i\right)}\left(r,\theta \right)=-\frac{{\varphi }^{\left(o\right)}}{{\varphi }^{\left(i\right)}}\frac{U_{\infty }}{{\kappa }^{\left(o\right)}}\left(\frac{3+12\frac{k^{\left(o\right)}}{a^2}}{2+\frac{k^{\left(i\right)}}{k^{\left(o\right)}}+12\frac{k^{\left(i\right)}}{a^2}}\right)r\text{cos}\theta ,r<a$$ E19

$$p^{\left(o\right)}\left(r,\theta \right)=-\frac{U_{\infty }}{{\kappa }^{\left(o\right)}}r\text{cos}\theta +\frac{U_{\infty }}{{\kappa }^{\left(o\right)}}\frac{a^3}{r^2}\left(\frac{\frac{k^{\left(i\right)}}{k^{\left(o\right)}}-1}{2+\frac{k^{\left(i\right)}}{k^{\left(o\right)}}+12\frac{k^{\left(i\right)}}{a^2}}\right)\text{cos}\theta ,r>a$$ E20

For defects with radius $r\sim {10}^{-2},{10}^{-3}\mathrm{m}$ , the quantities $k^{\left(i\right)}/a^2$ and $k^{\left(o\right)}/a^2$ are very small compared to unit ( $\simeq {10}^{-5},{10}^{-6}$ ) and can be neglected.

The velocity is deduced from the pressure due to Darcy’s law.

• The pressure inside the inclusion (19) describes a constant velocity field. When the porosities of host medium and inclusion are substantially equal, then the fluid velocity is uniform (constant and aligned with the applied pressure gradient) with value

$$v^{\left(i\right)}=U_{\infty }\frac{3}{1+2\frac{{\kappa }^{\left(o\right)}}{{\kappa }^{\left(i\right)}}}.$$ E21

If ${\kappa }^{\left(i\right)}>{\kappa }^{\left(o\right)}$ , then $v^{\left(i\right)}>U_{\infty }$ . In this case, the fluid passes preferentially through the inclusion and in the vicinity of the inclusion, the streamlines in the host medium are curved toward the inclusion. We have the inverse conclusion if ${\kappa }^{\left(i\right)}<{\kappa }^{\left(o\right)}$ . When ${\kappa }^{\left(o\right)}\to 0$ , i.e., for low permeability $k^{\left(o\right)}$ , then $v^{\left(i\right)}\sim 3U_{\infty }$ . If ${\kappa }^{\left(i\right)}\to 0$ , then $v^{\left(i\right)}\to 0$ .

Eq. (21) can also be written in the form:

$$v^{\left(i\right)}=U_{\infty }\frac{{\kappa }^{\left(i\right)}}{{\kappa }^{\left(o\right)}+\frac{1}{3}\left({\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}\right)}$$ E22

which will be generalized for the ellipsoidal inclusion.

• Outside the inclusion, the pressure is the sum of the applied pressure plus a dipolar contribution due to a induced dipole centered at the origin, the dipolar moment ${\mathcal{P}}_{\mathit{sph}}$ of which is

$${\mathcal{P}}_{\mathit{sph}}=4\pi a^3{U}_{\infty }{\kappa }^{\left(i\right)}\left(\frac{\frac{k^{\left(i\right)}}{k^{\left(o\right)}}-1}{2+\frac{k^{\left(i\right)}}{k^{\left(o\right)}}}\right).$$ E23

The corresponding density of induced “hydraulic” surface charges is

$${\sigma }_{\mathit{sph}}\left(\theta \right)=4\pi a^3{U}_{\infty }{\kappa }^{\left(i\right)}\left(\frac{\frac{k^{\left(i\right)}}{k^{\left(o\right)}}-1}{2+\frac{k^{\left(i\right)}}{k^{\left(o\right)}}}\right)\text{cos}\left(\theta \right).$$ E24

In the host medium ( $r>a$ ), the components of the seepage rate are:

$$v_r^{\left(o\right)}=-U_{\infty }\text{cos}\theta \left(1+2\frac{a^3}{r^3}\left(\frac{\frac{{\kappa }^{\left(i\right)}}{{\kappa }^{\left(o\right)}}-1}{2+\frac{{\kappa }^{\left(i\right)}}{{\kappa }^{\left(o\right)}}}\right)\right),$$ E25

$$v_{\theta }^{\left(o\right)}=U_{\infty }\text{sin}\theta \left(1-\frac{a^3}{r^3}\left(\frac{\frac{{\kappa }^{\left(i\right)}}{{\kappa }^{\left(o\right)}}-1}{2+\frac{{\kappa }^{\left(i\right)}}{{\kappa }^{\left(o\right)}}}\right)\right).$$ E26

Figure 4 represents the seepage rate in the porous medium for ${\kappa }^{\left(i\right)}>{\kappa }^{\left(o\right)}$ . Figure 4a shows the levels of the amplitude of the velocity, while Figure 4b shows its stream lines.

Figure 4.

Fluid flow through an inclusion in a porous medium: (a) velocity levels and (b) streamlines ( ${\kappa }^{\left(i\right)}/{\kappa }^{\left(o\right)}>1$ ).

So the response developed by a defect when submitted to a pressure gradient is an induced dipole $\mathcal{P}$ . When dealing with linear phenomena (low filtration speed), this response is proportional to its cause, the proportionality factor depending only on the shape of the defect, its volume, and on the ratio of its mobility to that of the host medium. In our case, the dipole moment is

$${\mathcal{P}}_{\mathit{sph}}=\alpha \left(-U_{\infty }\right),$$ E27

where $\alpha$ is a susceptibility which assesses the polarizability, i.e., the capacity of the porous inclusion to induce a dipole ${\mathcal{P}}_{\mathit{sph}}$ under the action of excitation $-U_{\infty }/{\kappa }^{\left(o\right)}$ . For a spherical inclusion of volume V, the susceptibility $\alpha$ is:

$$\alpha =3V{\kappa }^{\left(i\right)}\left(\frac{\frac{k^{\left(i\right)}}{k^{\left(o\right)}}-1}{2+\frac{k^{\left(i\right)}}{k^{\left(o\right)}}}\right).$$ E28

Often, the quantity $\chi =\alpha /{\kappa }^{\left(i\right)}V$ which does not depend on the volume of the inclusion is more relevant. Its variations as function of the ratio $\kappa ={\kappa }^{\left(i\right)}/{\kappa }^{\left(o\right)}$ are shown in Figure 5. They range from −3/2 when ${\kappa }^{\left(i\right)}<<{\kappa }^{\left(o\right)}$ to 3 when ${\kappa }^{\left(i\right)}>>{\kappa }^{\left(o\right)}$ . In the following, it is convenient to put the external pressure in the form:

Figure 5.

Hydraulic polarizability $\chi =\alpha /{\kappa }^{\left(i\right)}V$ of a spherical inclusion vs log (k).

$$p^{\left(o\right)}\left(r,\theta \right)=-\frac{U_{\infty }}{{\kappa }^{\left(o\right)}}\mathit{\text{rcosθ}}+\frac{P_{\mathit{sph}}}{4\pi }\frac{\mathit{\text{cosθ}}}{r^2}$$ E29

2.2.2 Ellipsoidal defect

In addition to the interest that ellipsoidal inclusion has an exact analytical solution, its study (its study) allows us to understand the effects of the shape of the defects on the fluid motion. Indeed, ellipsoidal surface can be seen as a generic element of a set of volumes comprising the disc, the sphere, and the oblong shape (needle) and so the nonsphericity can be appreciated through the values of its polarizability. The general ellipsoidal inclusion having semiaxes a, b, and c aligned with the axes of the Cartesian coordinates system and centered at the origin is described by the following equation:

$$\frac{x^2}{u+a^2}+\frac{y^2}{u+b^2}+\frac{z^2}{u+c^2}=1,$$ E30

where $x$ , $y$ , and $z$ are the position coordinates of any point on the surface of the ellipsoid. Eq. (30) has three roots $\xi$ , $\eta$ , and $\zeta$ which define the surfaces coordinates: surfaces with constant $\xi$ are ellipsoids, while surfaces with constant $\eta$ or $\zeta$ are hyperboloids. Surfaces of confocal ellipsoids are described by adjusting the scalar u. So, two ellipsoids defined by (30) with $u=u_1$ and $u=u_2$ are called confocal if their semi axes obey to the conditions

$$a_1^2-a_2^2=b_1^2-b_2^2=c_1^2-c_2^2$$ E31

where $a_i$ , $b_i$ , and $c_i$ are their respective semiaxes. Some results related to the ellipsoid are given in Appendix A.

Let $p\left(\mathbf{r}\right)$ be the pressure with a constant gradient directed along the $\mathit{Ox}$ axis applied to the porous medium (Figure 6). In absence of defect, its expression is:

Figure 6.

Ellipsoidal inclusion in a fluid flow.

$$p\left(\mathbf{r}\right)=p_0\left(x\right)=\mathit{Ex}=E\sqrt{\frac{\left(a^2+\xi \right)\left(a^2+\eta \right)\left(a^2+\zeta \right)}{\left(b^2-a^2\right)\left(c^2-a^2\right)}}$$ E32

where $E=-U_{\infty }/{\kappa }^{\left(o\right)}$ . In this relation, $\xi$ , $\eta$ , and $\zeta$ are the ellipsoidal coordinates given in Appendix A. Here, the field $p_0\left(x\right)$ can be viewed as an “incident pressure.”

When the inclusion is embedded in the porous medium, it produces perturbations in the fluid motion giving rise to the “scattered” pressure. When the filtration rate is low, the scattering by the inclusion is a linear phenomenon, leading to the following expression of the scattered pressure:

$$p_{\mathit{sc}}\left(\mathbf{r}\right)=p_0\left(x\right)F\left(\xi \right).$$ E33

where $F\left(\xi \right)$ is a proportional coefficient. For $p_{\mathit{sc}}\left(\mathbf{r}\right)$ to be a solution of the Laplace equation, it must verify the differential equation:

$$\frac{d^{2}F}{d{\xi }^2}+\frac{\mathit{dF}}{\mathit{d\xi }}\frac{d}{\mathit{d\xi }}\left(\text{ln}\left[R\left(\xi \right)\left(\xi +a^2\right)\right]\right)=0$$ E34

where $R\left(\sigma \right)=\left(a^2+\sigma \right)\left(b^2+\sigma \right)\left(c^2+\sigma \right)$ . $F\left(\xi \right)$ is then the sum of the two functions $F_1\left(\xi \right)+F_2\left(\xi \right)$ , where $F_1\left(x\right)=A$ is a constant and $F_2\left(x\right)$ is

$$F_2\left(\xi \right)={\int }_{\xi }^{\infty }\frac{\mathit{d\sigma }}{\left(\sigma +a^2\right)R\left(\sigma \right)}.$$ E35

Thus, the pressure outside the inclusion is

$$p\left(\mathbf{r}\right)=p_0\left(x\right)\left[A-\frac{B}{2}{\int }_{\xi }^{\infty }\frac{\mathit{d\sigma }}{\left(\sigma +a^2\right)R\left(\sigma \right)}\right].$$ E36

The two constants $A$ and $B$ are determined by the boundary conditions at the surface $\xi =0$ of the inclusion:

1. continuity of fluid flow at $\xi =0$ :

$${\varphi }^{\left(e\right)}{\kappa }^{\left(e\right)}\frac{1}{h_{\xi }}\frac{\partial p^{\left(e\right)}}{\partial \xi }={\varphi }^{\left(i\right)}{\kappa }^{\left(i\right)}\frac{1}{h_{\xi }}\frac{\partial p^{\left(i\right)}}{\partial \xi },$$ E37

2. continuity of the stress component ${\tau }_{\mathit{\xi \xi }}$ in fluid at $\xi =0$ :

$${\varphi }^{\left(e\right)}\left(-p^{\left(e\right)}+2\eta \frac{1}{h_{\xi }}\frac{\partial v_{\xi }^{\left(e\right)}}{\partial \xi }\right)={\varphi }^{\left(i\right)}\left(-p^{\left(i\right)}+2\eta \frac{1}{h_{\xi }}\frac{\partial v_{\xi }^{\left(i\right)}}{\partial \xi }\right)$$ E38

where $v_{\xi }$ is the normal component of the seepage flow to the surface $\xi =0$ given by the Darcy’s law:

$$v_{{\lambda }_i}=-\kappa \frac{1}{h_{{\lambda }_i}}{\partial }_{{\lambda }_i}p,{\lambda }_i=\xi ,\eta ,\zeta .$$ E39

In these relations, the coefficients $h_{{\lambda }_i}$ are the scale factors of ellipsoidal coordinates given in Appendix A. We obtain the values of the pressure such that:

1. within the inclusion:

$$p^{\left(i\right)}\left(x\right)=\frac{{\kappa }^{\left(o\right)}}{{\kappa }^{\left(o\right)}+N_x\left({\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}\right)}{p}_0\left(x\right),$$ E40

1. outside the inclusion:

$$p^{\left(o\right)}\left(x\right)=p_0\left(x\right)+p_{\mathit{sc}}\left(x\right)=p_0\left(x\right)-\frac{\mathit{abc}\left({\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}\right)}{{\kappa }^{\left(o\right)}+N_x\left({\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}\right)}{p}_0\left(x\right)F_2\left(\xi \right).$$ E41

Far from the center of the inclusion, $\xi \approx r^2$ , the scattered pressure can be approximate by:

$$\begin{array}{c}{p}_{\mathit{sc}}\left(x\right)\approx \frac{\mathit{abc}\left({\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}\right)}{{\kappa }^{\left(o\right)}+N_x\left({\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}\right)}{p}_0\left(x\right){\int }_{r^2}\frac{\mathit{dr}}{r^{5/2}}\\ \approx \frac{\mathit{abc}\left({\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}\right)}{{\kappa }^{\left(o\right)}+N_x\left({\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}\right)}\frac{x}{3r^3}E.\end{array}$$ E42

The right-hand side of Eq. (42) is the expression of pressure produced by a dipole aligned with the axis Ox. From the expression of the speed within the inclusion:

$$v^{\left(i\right)}=U_{\infty }\frac{{\kappa }^{\left(o\right)}}{{\kappa }^{\left(o\right)}+N_x\left({\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}\right)},$$ E43

the dipole moment is then

$${\mathcal{P}}_{\mathit{\text{elli}},a}=-{\mathit{VU}}_{\infty }{\kappa }^{\left(i\right)}\frac{\left({\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}\right)}{{\kappa }^{\left(o\right)}+N_x\left({\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}\right)}.$$ E44

Here, V is the volume of the ellipsoidal inclusion. The factor N_x,

$$N_x=\frac{\mathit{abc}}{3}{\int }_0^{\infty }\frac{\mathit{d\sigma }}{\left(\sigma +a^2\right)\sqrt{\sigma +a^2)\left(\sigma +b^2\right)\left(\sigma +c^2\right)}},$$ E45

describes how the dipole moment of the inclusion changes with its shape and its orientation in relation with the incident pressure field. The geometric parameters N_x, N_y, and N_z appear for the first time in hydrodynamics [8] to describe the disturbance brought by a solid immersed in an infinite fluid in uniform motion. Their values were computed by Stoner [9] and Osborn [10]. The name “depolarization factors” comes from electromagnetism (see, for example, Landau and Lifchitz [11]).

From Eq. (163), it is possible to find the values of the depolarization factors of some particular inclusions such as:

• spherical inclusion $N_x=N_y=N_z=1/3$ ,

• inclusion in disc form (axis Oz) $N_x=N_y=0,N_z=1$ , and

• oblong inclusion $N_x=N_y=1/2,N_z=0$ .

Figure 7 is a plot of depolarization factors of a ellipsoidal inclusion having two equal semiaxes according to their ratio.

Figure 7.

Depolarization factors for a ellipsoidal inclusion with semiaxes (a, a, b) as functions of a/b.

The expression of $p_{\mathit{sc}}\left(x\right)$ given by (42) is not exact since it is a result of the approximation $\xi \sim r^2$ , i.e., far from the inclusion, where only the dipolar effects are relevant.

In the case of an ellipsoidal inclusion, the polarizability is no longer a scalar but is a tensor. Its eigenvalues are polarizabilities along the axes of the ellipsoid. So, we can write the dipole moment (44) as ${\mathcal{P}}_{\mathit{\text{elli}},a}={\alpha }_a{U}_{\infty }$ , where ${\alpha }_a$ is the eigenvalue of the tensor polarizability along its principal direction Ox which defines the polarizability along this axis. In Figure 8, we depict the variations of the susceptibility $\chi$ as function of ${\kappa }^{\left(i\right)}/{\kappa }^{\left(o\right)}$ for different sets of the depolarization factors when incident pressure is along Ox.

Figure 8.

Hydraulic polarizability $\chi ={\alpha }_x/V$ of an ellipsoidal inclusion vs. $\text{log}\left({\kappa }^{\left(i\right)}/{\kappa }^{\left(o\right)}\right)$ for different values of the depolarization factors ( $N_x$ , $N_y$ , $N_z$ ) ( $-•-$ : $N_x=0.1,N_y=0.1,N_z=0.8$ ); ( $-\diamond -$ : $N_x=0.2,N_y=0.2,N_z=0.6$ ); ( $-◀-$ : $N_x=0.4,N_y=0.4,N_z=0.2$ ); ( $-■-$ : $N_x=0.45,N_y=0.45,N_z=0.1$ ); (– – –: $N_x=1/3,N_y=1/3,N_z=1/3$ ).

The pressure outside the inclusion is then:

$$p^{\left(o\right)}\left(r,\theta \right)=-\frac{U_{\infty }}{{\kappa }^{\left(o\right)}}r\text{cos}\theta +\frac{{\mathcal{P}}_{\mathit{\text{elli}},a}}{4\pi }\frac{\text{cos}\theta }{r^2}.$$ E46

This relation is similar to (29), differing from it only by the expression of the dipole moment. The fundamental difference between the spherical and ellipsoidal inclusions is that the pressure scattered by the sphere contains only a dipolar field, whereas in strictness, the ellipsoid also scatters high-order multipolar fields. We can then deduce from this remark that the more the shape of the inclusion is distant from that of the sphere, the more the scattered pressure contains high-order multipolar terms.

The result obtained in (42) does not show these terms since the calculation of the integral F₂ is an approximate computation when $\xi /r\gg 1$ . When we move away from the ellipsoidal inclusion, the multipolar terms of order greater than 2 decrease very quickly, leaving only the contribution of the incident pressure and the dipole term of the scattered pressure.

2.3 Inhomogeneous defect

In fact, defects rarely exhibit homogeneous structure. The parameter that characterizes the defect (in our case the permeability) is generally a variable varying according to a law which depends on the way in which the defect develops.

For the spherical defect, the simplest situation is the radial variation of the permeability. The fundamental difference between homogeneous and inhomogeneous spherical inclusions is that in the latter case, the velocity field loses its uniformity. The determination of the dipole moment requires a different approach from that previously developed. Two cases are considered: (i) the permeability is a piecewise constant function and (ii) the permeability is a continuously varying function.

2.3.1 Layered spherical defect

Consider an inhomogeneous sphere of porous medium embedded in a homogeneous host medium. We assume that the permeability of the sphere depends only on the radius and is a piecewise constant function, i.e., the sphere is a set of nested spherical layers. The permeability of the background medium is $k^{\left(o\right)}$ , that of the outermost layer is $k^{\left(1\right)}$ , and so on to the central sphere whose permeability is $k^{\left(i\right)}$ .

To calculate the perturbation of the incident pressure due to the sphere seen as a scatter, we proceed the following: the pressure field is calculated in each layer of the sphere. The pressure field in the layer number n is related to those in the layer number $n+1$ and $n-1$ by the boundary conditions. It is assumed that the defect is a set of N concentric spherical layers in which the value of permeability is constant $k^{\left(n\right)}$ . Let ${\kappa }^{\left(n\right)}$ denotes the ratio $k^{\left(n\right)}/\eta$ in the layer number n delimited by the spheres of radii $r_n$ and $r_{n+1}$ such that $r_n>r_{n+1}$ . The core of the sphere has index $i=N$ (Figure 9). The determination of pressure and velocity in this type of inclusion consists in solving the Laplace equation $\mathrm{\Delta }{p}^{\left(n\right)}=0$ in each layer and in connecting the solutions using the boundary conditions: continuity of fluid flow and that of the radial component of the stress.

Figure 9.

Layered sphere.

In the layer number n, the pressure is noted:

$$p^{\left(n\right)}=A^{\left(n\right)}r\text{cos}\theta +B^{\left(n\right)}{r}^{-2}\text{cos}\theta .$$ E47

The coefficients A and B of two consecutive layers are connected by the following conditions in $r=r_{n+1}$ :

$${\varphi }^{\left(n\right)}{v}_r^{\left(n\right)}\left(r=r_{n+1}\right)={\varphi }^{\left(n+1\right)}{v}_r^{\left(n+1\right)}\left(r=r_{n+1}\right),$$ E48

$${\tau }_{\mathit{rr}}^{\left(n\right)}\left(r=r_{n+1}\right)={\tau }_{\mathit{rr}}^{\left(n+1\right)}\left(r=r_{n+1}\right).$$ E49

Transfer matrix: The linearity of the problem makes it possible to write that the pairs of coefficients ( $A_n$ , $B_n$ ) and ( $A_{n+1}$ , $B_{n+1}$ ) are linked by a matrix equation such that:

$$\left(\begin{array}{l}{A}^{\left(n\right)}\\ B^{\left(n\right)}\end{array}\right)=T_{n,n+1}\left(\begin{array}{l}{A}^{\left(n+1\right)}\\ B^{\left(n+1\right)}\end{array}\right),$$ E50

where $T_{n,n+1}$ is the transfer matrix between the two consecutive layers $n$ and $n+1$ , the entries of which are:

$$T_{11}=\frac{{\varphi }^{\left(n+1\right)}}{{\varphi }^n}\left(\frac{k^{\left(n+1\right)}}{k^n}+2\frac{1-\frac{k^{\left(n+1\right)}}{k^n}}{3+12k^n{r}_n^{-2}}\right),$$ E51

$$T_{12}=2\frac{{\varphi }^{\left(n+1\right)}}{{\varphi }^n}{r}_n^{-3}\left(-\frac{k^{n+1}}{k^n}+\frac{1+2\frac{k^{\left(n+1\right)}}{k^n}+12k^{\left(n+1\right)}{r}_n^{-2}}{3+12k^n{r}_n^{-2}}\right),$$ E52

$$T_{21}=\frac{{\varphi }^{\left(n+1\right)}}{{\varphi }^n}{r}_n^3\left(\frac{1-\frac{k^{\left(n+1\right)}}{k^n}}{3+12k^n{r}_n^{-2}}\right),$$ E53

$$T_{22}=\frac{{\varphi }^{\left(n+1\right)}}{{\varphi }^n}\left(\frac{1+2\frac{k^{\left(n+1\right)}}{k^n}+12k^{\left(n+1\right)}{r}_n^{-2}}{3+12k^n{r}_n^{-2}}\right).$$ E54

From $A^0=-U_{\infty }/{\kappa }^{\left(o\right)}$ and $B^{\left(i\right)}=0$ , it is possible to get the pressure in each layer of the sphere.

The effects of two consecutive layers are obtained by the product of the transfer matrices of each of these layers. So from the matrix equations:

$$\left(\begin{array}{l}{A}^{\left(n\right)}\\ B^{\left(n\right)}\end{array}\right)=T_{n,n+1}\left(\begin{array}{l}{A}^{\left(n+1\right)}\\ B^{\left(n+1\right)}\end{array}\right),$$ E55

$$\left(\begin{array}{l}{A}^{\left(n+1\right)}\\ B^{\left(n+1\right)}\end{array}\right)=T_{n+1,n+2}\left(\begin{array}{l}{A}^{\left(n+2\right)}\\ B^{\left(n+2\right)}\end{array}\right),$$ E56

we get:

$$\left(\begin{array}{l}{A}^{\left(n\right)}\\ B^{\left(n\right)}\end{array}\right)=T_{n,n+2}\left(\begin{array}{l}{A}^{\left(n+2\right)}\\ B^{\left(n+2\right)}\end{array}\right),$$ E57

where $T_{n,n+2}=T_{n,n+1}{T}_{n+1,n+2}$ .

Scattering matrix: Another way of linking the coefficients A and B of two consecutive layers is the use of the scattering matrix $S_{n,n+1}$ . The scattering matrix is sometimes more efficient than the transfer matrix for calculating the amplitudes of the waves reflected and transmitted by an object subjected to incident waves. It is such that:

$$\left(\begin{array}{l}{B}^{\left(n\right)}\\ A^{\left(n+1\right)}\end{array}\right)=S_{n,n+1}\left(\begin{array}{l}{A}^{\left(n\right)}\\ B^{\left(n+1\right)}\end{array}\right).$$ E58

Its entries are:

$$S_{n,n+1}=\frac{1}{{\kappa }^{\left(n+1\right)}+2{\kappa }^{\left(n\right)}}\left(\begin{array}{ll}\left({\kappa }^{\left(n+1\right)}-{\kappa }^{\left(n\right)}\right)r_{n+1}^3& 3{\kappa }^{\left(n+1\right)}\\ 3{\kappa }^{\left(n\right)}& 2\left({\kappa }^{\left(n+1\right)}-{\kappa }^{\left(n\right)}\right)r_{n+1}^{-3}\end{array}\right).$$ E59

The scattering matrix of two consecutive layers is given by their Redheffer product. From

$$\left(\begin{array}{l}{B}^{\left(n\right)}\\ A^{\left(n+1\right)}\end{array}\right)=S_{n,n+1}\left(\begin{array}{l}{A}^{\left(n\right)}\\ B^{\left(n+1\right)}\end{array}\right),$$ E60

$$\left(\begin{array}{l}{B}^{\left(n+1\right)}\\ A^{\left(n+2\right)}\end{array}\right)=S_{n+1,n+2}\left(\begin{array}{l}{A}^{\left(n+1\right)}\\ B^{\left(n+2\right)}\end{array}\right),$$ E61

we get:

$$\left(\begin{array}{l}{B}^{\left(n\right)}\\ A^{\left(n+2\right)}\end{array}\right)=S_{n,n+2}\left(\begin{array}{l}{A}^{\left(n\right)}\\ B^{\left(n+2\right)}\end{array}\right),$$ E62

where $S_{n,n+2}={S_{n,n+1}}^{\ast }{S}_{n+1,n+2}$ . In this relation, the right-hand side is the Redheffer star product. For more details about the Redheffer star product, we can refer to [12].

2.3.2 Spherical inclusion with continuously variable permeability

When $\kappa$ is a continuous function of the variable r, the matrix Eq. (50) becomes the system of differential equations:

$$\frac{d}{\mathit{dr}}\left(\frac{r^3{\kappa }^2\mathit{dA}/\mathit{dr}}{\mathit{d\kappa }/\mathit{dr}}\right)+r^2\mathit{\kappa A}=0$$ E63

$$\frac{d}{\mathit{dr}}\left(\frac{{\kappa }^2\mathit{dB}/\mathit{dr}}{r^3\mathit{d\kappa }/\mathit{dr}}\right)-\frac{2\kappa }{r^4}B=0.$$ E64

2.3.3 Layered ellipsoidal inclusion

The generalization of the radial variation of the permeability of the spherical inclusion to the ellipsoidal requires that the permeability only depends on $\xi$ . This is true in orthogonal directions at its surface. This condition entails that inside the ellipsoid, the strata are limited by confocal ellipsoidal surfaces $\xi ={\xi }_k$ , i.e., having the same foci as the surface $\xi$ ; hence, their semiaxes are related by the following relations:

$$a^2-a_k^2=b^2-b_k^2=c^2-c_k^2.$$ E65

Consider a porous inhomogeneous ellipsoidal inclusion having the permeability $k^{\left(i\right)}$ embedded in a background medium of homogeneous mobility ${\kappa }^{\left(o\right)}$ . We assume that the mobility of the inclusion is stratified, i.e., it is a constant piecewise function and each layer has its own mobility ${\kappa }^{\left(n\right)}$ . The mobility in the outermost layer is ${\kappa }^{\left(1\right)}$ , that of the next layer is ${\kappa }^{\left(2\right)}$ , etc. to the central layer whose the mobilty is ${\kappa }^{\left(i\right)}$ . The layers are limited by the confocal surfaces $\xi ={\xi }_k$ whose semiaxes obey (65) and are numbered from 1 to $i=N$ from the outside to the inside, such that the strata $n$ and $n+1$ have the common boundary $\xi ={\xi }_{n+1}$ (Figure 10). In each of these strata, the pressure is the solution of the Laplace equation given by:

Figure 10.

Layered ellipsoidal inclusion.

$$p^{\left(j\right)}\left(\mathbf{r}\right)=-\mathit{Ex}\left[A^{\left(j\right)}-\frac{B^{\left(j\right)}}{2}{\int }_{\xi }^{\infty }\frac{\mathit{d\sigma }}{(\left(\sigma +a^2\right)R^{\left(j\right)}\left(\sigma \right)}\right]$$ E66

where

$$R^{\left(j\right)}\left(\sigma \right)=\left(a_j^2+\sigma \right)\left(b_j^2+\sigma \right)\left(c_j^2+\sigma \right).$$ E67

In the layers $n$ and $n+1$ , the solutions of the Laplace equations are:

$$p^{\left(n\right)}\left(\mathbf{r}\right)=-\mathit{Ex}\left[A^{\left(n\right)}-\frac{B^{\left(n\right)}}{2}{\int }_{\xi }^{\infty }\frac{\mathit{d\sigma }}{(\left(\sigma +a^2\right)R^{\left(n\right)}\left(\sigma \right)}\right],$$ E68

$$p^{\left(n+1\right)}\left(\mathbf{r}\right)=-\mathit{Ex}\left[A^{\left(n+1\right)}-\frac{B^{\left(n+1\right)}}{2}{\int }_{\xi }^{\infty }\frac{\mathit{d\sigma }}{(\left(\sigma +a^2\right)R^{\left(n+1\right)}\left(\sigma \right)}\right]$$ E69

with the boundary condition at $\xi ={\xi }_{n+1}$ :

$${\varphi }^{\left(n\right)}\frac{1}{h_{\xi \left(n\right)}}\frac{\partial p^{\left(n\right)}}{\partial \xi }={\varphi }^{\left(n+1\right)}\frac{1}{h_{\xi \left(n+1\right)}}\frac{\partial p^{\left(n+1\right)}}{\partial \xi },$$ E70

and

$${\varphi }^{\left(n\right)}\left(-p^{\left(n\right)}+2\eta \frac{\partial v_{\xi }^{\left(n\right)}}{h_{\xi \left(n\right)}\partial \xi }\right)={\varphi }^{\left(n+1\right)}\left(-p^{\left(n+1\right)}+2\eta \frac{\partial v_{\xi }^{\left(n+1\right)}}{h_{\xi \left(n+1\right)}\partial \xi }\right).$$ E71

By proceeding in the same way as for the spherical cavity, one finds the transfer matrix and the scattering matrix.

Transfer matrix: The transfer matrix of the ellipsoidal inclusion is:

$$T_{n,n+1}=\frac{1}{{\kappa }^{\left(n\right)}}\left(\begin{array}{ll}{\kappa }^{\left(n\right)}+N_{x,n+1}\left({\kappa }^{\left(n+1\right)}-{\kappa }^{\left(n\right)}\right)& \frac{N_{x,n+1}\left(1-N_{x,n+1}\right)}{a_{n+1}{b}_{n+1}{b}_{n+1}}\left({\kappa }^{\left(n+1\right)}-{\kappa }^{\left(n\right)}\right)\\ a_{n+1}{b}_{n+1}{c}_{n+1}\left({\kappa }^{\left(n+1\right)}-{\kappa }^{\left(n\right)}\right)& {\kappa }^{\left(n+1\right)}+N_{x,n+1}\left({\kappa }^{\left(n\right)}-{\kappa }^{\left(n\right)}\right)\end{array}\right).$$ E72

Scattering matrix: The scattering matrix of the ellipsoidal inclusion is:

$$\begin{array}{l}{S}_{n,n+1}=\frac{1}{{\kappa }^{\left(n\right)}+N_{x,n+1}\left({\kappa }^{\left(n+1\right)}-{\kappa }^{\left(n\right)}\right)}\\ \left(\begin{array}{ll}\left({\kappa }^{\left(n+1\right)}-{\kappa }^{\left(n\right)}\right)V_{n+1}& {\kappa }^{\left(n+1\right)}\\ {\kappa }^{\left(n\right)}& N_{x,n+1}\left(1-N_{x,n+1}\right)\left({\kappa }^{\left(n+1\right)}-{\kappa }^{\left(n\right)}\right)V_{n+1}^{-1}\end{array}\right)\end{array}$$ E73

where $V_n=a_n{b}_n{c}_n$ .

Dipole moment: When we are only interested in the scattered far field, the inclusion can be replaced by an equivalent dipole. When r is large in front of the lengths of the axes of the ellipsoid ( $r\gg a,b,c$ ), $\xi$ is of the order of r²:

$$\xi \simeq r^2$$ E74

and the dipolar term is

$$\frac{B^{\left(o\right)}}{2}{\int }_{\xi }^{\infty }\frac{\mathit{d\sigma }}{\left(\sigma +a^2\right)R\left(\sigma \right)}\simeq \frac{B^{\left(o\right)}}{2}{\int }_{r^2}^{\infty }\frac{\mathit{d\sigma }}{{\sigma }^{5/2}}=\frac{B^{\left(o\right)}}{3r^3}$$ E75

The component ${\mathcal{P}}_{\mathit{\text{elli}},a}$ of dipole moment along the direction of the fluid flow is obtained by identification of the terms in $r^{-2}$ in relations (46) and (75), namely:

$${\mathcal{P}}_{\mathit{\text{elli}},a}=\frac{4\pi }{3}{U}_{\infty }{B}^{\left(o\right)}.$$ E76

2.3.4 Inclusion with continuously variable permeability

For ellipsoidal inclusion, we assume that mobility depends only on the variable $\xi$ . From (33), the problem comes down to finding of the differential equation of the function $F\left(\xi \right)$ . Eq. 33 is then:

$$\frac{d}{\mathit{d\xi }}\left(\kappa \left(\xi \right)\left(\xi +a^2\right)R\left(\xi \right)\frac{\mathit{dF}}{\mathit{d\xi }}\right)+\frac{\mathit{d\kappa }\left(\xi \right)}{\mathit{d\xi }}\frac{R\left(\xi \right)}{2}F\left(\xi \right)-\frac{\mathit{d\kappa }\left(\xi \right)}{\mathit{d\xi }}\frac{R\left(\xi \right)}{2}E=0.$$ E77

It is easy to verify that when $\kappa$ is constant, we find the case of the homogeneous inclusion, and that if we put $a=b=c$ , then we find the result of spherical inclusion.

2.4 Anisotropic defects

Often the defects occurring in porous media are anisotropic, i.e., some of their physical parameters like permeability are no longer scalar quantities but are tensors. For an anisotropic porous medium, assuming the Einstein convention, the Darcy’s law is

$$v_i=-\frac{k_{\mathit{ij}}}{\eta }{\partial }_{j}p\text{where}{\partial }_{j}p=\frac{\partial p}{\partial x_j}.$$ E78

The permeability is then defined by nine components k_ij, i.e., it has different values in different directions of the space. Liakopoulos [13] had shown that the permeability is a symmetric tensor of second rank. This leads to great simplifications for the study of such porous media. If in isotropic media the fluid velocity is aligned with the hydraulic gradient, in anisotropic media, this is true only along the principal directions of the tensor. It is therefore not surprising that the flow movement of the fluid is seriously disturbed by this type of defects.

In a 3D space, the permeability tensor has three principal directions perpendicular to each other and for which the permeability corresponds to the tensor eigenvalues. In the coordinates system defined by these directions, the permeability tensor is diagonal:

$$\mathbf{k}=\left(\begin{array}{lll}{k}_1& 0& 0\\ 0& k_2& 0\\ 0& 0& k_3\end{array}\right).$$ E79

By its definition, mobility inherits properties of symmetry of permeability and is therefore a symmetric tensor such that ${\kappa }_{\mathit{ij}}=k_{\mathit{ij}}/\eta$ .

The principal directions of the permeability tensors of the host medium and of the inclusion define coordinate systems which generally do not coincide. The system linked to the host environment is called the primary system, while that of the inclusion is called the secondary system.

In this part, the study of the effects of an anisotropic spherical inclusion in a porous medium explores three different configurations:

1. the host medium is isotropic and the defect is anisotropic,

2. the host medium is anisotropic and the defect is isotropic, and

3. the host medium and the defect are anisotropic.

2.4.1 Anisotropic defect in isotropic medium

In this configuration, the primary system is such that the incident pressure gradient is along the Oz axis and the secondary coordinate system is defined by the principal directions of the permeability tensor. Then, let ( ${\theta }_0$ , ${\phi }_0$ ) and ( $\theta$ , $\phi$ ) be the angular directions of the incident pressure gradient $\mathbf{\text{∇}}{p}^{\left(o\right)}$ and of the observation vector $\mathbf{OM}=\mathbf{r}$ in the primary system. We note $\beta$ as the angle between these two directions.

The external pressure verifies the classical Laplace equation, while the internal one is solution of the following equation:

$${\kappa }_1^{\left(i\right)}\frac{\partial p^{\left(i\right)}}{\partial x^2}+{\kappa }_2^{\left(i\right)}\frac{\partial p^{\left(i\right)}}{\partial y^2}+{\kappa }_3^{\left(i\right)}\frac{\partial p^{\left(i\right)}}{\partial z^2}=0.$$ E80

This one is transformed into a Laplace equation as follows: at first, we dimensionalize the mobilities ${\kappa }_j^{\left(i\right)}$ by introducing the scalar quantity ${{\kappa }^{\prime }}^{\left(i\right)}$ . Eq. (80) then becomes

$${{\kappa }^{\prime }}^{\left(i\right)}\left({\kappa }_{r,1}^{\left(i\right)}\frac{\partial p^{\left(i\right)}}{\partial x^2}+{\kappa }_{r,2}^{\left(i\right)}\frac{\partial p^{\left(i\right)}}{\partial y^2}+{\kappa }_{r,3}^{\left(i\right)}\frac{\partial p^{\left(i\right)}}{\partial z^2}\right)=0,$$ E81

where ${\kappa }_{r,j}^{\left(i\right)}={\kappa }_j^{\left(i\right)}/{{\kappa }^{\prime }}^{\left(i\right)}$ . Using the linear transformation:

$$\frac{\partial }{\partial x_i^{\prime }}=\sqrt{{\kappa }_{r,i}^{\left(i\right)}}\frac{\partial }{\partial x_i},$$ E82

Eq. (80) takes the form:

$${{\kappa }^{\prime }}^{\left(i\right)}\left(\frac{{\partial }^2{p}^{\left(i\right)}}{{\partial }^2{x}^{\prime }}+\frac{{\partial }^2{p}^{\left(i\right)}}{{\partial }^2{y}^{\prime }}+\frac{{\partial }^2{p}^{\left(i\right)}}{{\partial }^2{z}^{\prime }}\right)=0.$$ E83

The solutions $p^{\left(i\right)}$ and $p^{\left(o\right)}$ , respectively, are:

$$p^{\left(i\right)}\left(r,\theta ,\phi \right)=\sum _{m,n}{A}_{m,n}^{\left(i\right)}{r}^n{P}_n^m\left(\cos \theta \right)\mathit{cos}\left(\mathit{m\phi }\right)+\sum _{m,n}{B}_{m,n}^{\left(i\right)}{r}^n{P}_n^m\left(\text{cos}\theta \right)\mathit{sin}\left(\mathit{m\phi }\right),$$ E84

$$p^{\left(o\right)}\left(r,\theta ,\phi \right)=\sum _{m,n}\left(A_{m,n}^{\left(o\right)}{r}^n+\frac{B_{m,n}^{\left(o\right)}}{r^{n+1}}\right)P_n^m\left(\text{cos}\theta \right)\mathit{cos}\left(\mathit{m\phi }\right)+\sum _{m,n}\left(C_{m,n}^{\left(o\right)}{r}^n+\frac{D_{m,n}^{\left(o\right)}}{r^{n+1}}\right)P_n^m\left(\text{cos}\theta \right)\mathit{sin}\left(\mathit{m\phi }\right).$$ E85

In (84), only the finite terms at $r=0$ appear.

The amplitudes $A_{m,n}^{\left(o\right)}$ , $B_{m,n}^{\left(o\right)}$ , $C_{m,n}^{\left(o\right)}$ , and $D_{m,n}^{\left(o\right)}$ , are determined by the boundary conditions at $r=a$ and when $r\to \infty$ .

When $r\to \infty$ , $p^{\left(o\right)}$ is:

$$p^{\left(o\right)}\to E_{\infty }r\text{cos}\beta ,\text{with}{E}_{\infty }=-\frac{U_{\infty }}{{\kappa }^{\left(o\right)}},$$ E86

where $\beta$ is the angle between vector ${\mathbf{E}}_{\infty }$ and the direction of the observer $\mathbf{OM}=\mathbf{r}$ . If ( $\theta$ , $\phi$ ) resp. ( ${\theta }_0$ , ${\phi }_0$ ) are the angular coordinates of the observer (resp. of the incident pressure gradient), then

$$\text{cos}\beta =\text{sin}\theta \text{sin}{\theta }_0\text{cos}\left(\phi -{\phi }_0\right)+\text{cos}\theta \text{cos}{\theta }_0.$$ E87

The expressions (86) and (87) show that, in the expansion (85), only nonzero terms are those for which $n=1$ . Taking into account the relations $P_1^0\left(\text{cos}\theta \right)=P_1\left(\text{cos}\theta \right)=\text{cos}\theta$ and $P_1^1\left(\text{cos}\theta \right)=\text{sin}\theta$ , we obtain:

$$p^{\left(o\right)}\left(r,\theta ,\phi \right)=A_{0,1}^{\left(o\right)}r\text{cos}\theta +A_{1,1}^{\left(o\right)}r\text{sin}\theta \text{cos}\phi +C_{1,1}^{\left(o\right)}r\text{sin}\theta \text{sin}\phi +\sum _{m,n}\frac{B_{m,n}^{\left(o\right)}}{r^{n+1}}{P}_n^m\left(\text{cos}\theta \right)\mathit{cos}\left(\mathit{m\phi }\right)+\sum _{m,n}\frac{D_{m,n}^{\left(o\right)}}{r^{n+1}}{P}_n^m\left(\text{cos}\theta \right)\mathit{sin}\left(\mathit{m\phi }\right).$$ E88

By identification with (86) with help of (87), we find:

$$A_{0,1}^{\left(o\right)}=E_{\infty }\text{cos}{\theta }_0,$$ E89

$$A_{1,1}^{\left(o\right)}=E_{\infty }\text{sin}{\theta }_0\text{cos}{\phi }_0,$$ E90

$$C_{1,1}^{\left(o\right)}=E_{\infty }\text{sin}{\theta }_0\text{sin}{\phi }_0.$$ E91

The boundary conditions at $r=a$ (continuity of the stress component ${\tau }_{\mathit{rr}}$ and conservation of the fluid flow through the inclusion surface) are:

$$p^{\left(o\right)}{\left|{}_{r=a}=p^{\left(i\right)}\right|}_{r=a}$$ E92

$${\kappa }^{\left(o\right)}\frac{\partial p^{\left(o\right)}}{\partial r}{\left|{}_{r=a}=\left({\kappa }_{11}^{\left(i\right)}\frac{\partial p^{\left(i\right)}}{\partial r}+{\kappa }_{12}^{\left(i\right)}\frac{\partial p^{\left(i\right)}}{r\partial \theta }+{\kappa }_{13}^{\left(i\right)}\frac{\partial p^{\left(i\right)}}{r\text{sin}\theta \partial \phi }\right)\right|}_{r=a}.$$ E93

In these equations, ${\kappa }_{\mathit{nm}}^{\left(i\right)}$ are the components of the tensor ${\kappa }^{\left(i\right)}$ in the spherical coordinates given in Appendix D. These relations lead to the following expressions of the pressure:

$$p^{\left(i\right)}\left(r,\theta ,\phi \right)=\frac{3A{\kappa }_{\left(o\right)}}{2{\kappa }^{\left(o\right)}+{\kappa }_3^{\left(i\right)}}r\text{cos}\theta +\frac{3B{\kappa }_{\left(o\right)}}{2{\kappa }^{\left(o\right)}+{\kappa }_1^{\left(i\right)}}r\text{sin}\theta \text{cos}\phi +\frac{3D{\kappa }_{\left(o\right)}}{2{\kappa }^{\left(o\right)}+{\kappa }_2^{\left(i\right)}}r\text{sin}\theta \text{sin}\phi ,$$ E94

$$p^{\left(o\right)}\left(r,\theta ,\phi \right)=\mathit{Ar}\text{cos}\theta +\mathit{Br}\text{sin}\theta \text{cos}\phi +\mathit{Dr}\text{sin}\theta \text{sin}\phi +\frac{{\mathit{Aa}}^3}{r^2}\frac{{\kappa }^{\left(o\right)}-{\kappa }_3^{\left(i\right)}}{2{\kappa }^{\left(o\right)}+{\kappa }_3^{\left(i\right)}}\text{cos}\theta +\frac{{\mathit{Ba}}^3}{r^2}\frac{{\kappa }^{\left(o\right)}-{\kappa }_1^{\left(i\right)}}{2{\kappa }^{\left(o\right)}+{\kappa }_1^{\left(i\right)}}\text{sin}\theta \text{cos}\phi +\frac{{\mathit{Da}}^3}{r^2}\frac{{\kappa }^{\left(o\right)}-{\kappa }_2^{\left(i\right)}}{2{\kappa }^{\left(o\right)}+{\kappa }_2^{\left(i\right)}}\text{sin}\theta \text{sin}\phi ,$$ E95

where

$$A=E_{\infty }\text{cos}{\theta }_0,B=E_{\infty }\text{sin}{\theta }_0\text{cos}{\phi }_0,D=E_{\infty }\text{sin}{\theta }_0\text{sin}{\phi }_0.$$ E96

The first three terms of the right-hand side of (95) are due to the pressure gradient applied to the porous medium. The last three terms are the pressure induced by the hydraulic dipoles directed along the principal directions of the anisotropic sphere.

When ${\mathbf{E}}_{\infty }$ is along the $\mathit{Oz}$ axis and for ${\phi }_0=0,{\theta }_0=0$ and ${\kappa }_1^{\left(i\right)}={\kappa }_2^{\left(i\right)}={\kappa }_3^{\left(i\right)}={\kappa }^{\left(i\right)}$ , we find the internal and external pressures of isotropic spherical inclusions (19) and (20).

Moreover, from the relations (94) and (95), it is possible to obtain the directions of the pressure gradient and of the velocity field inside the defect.

The inside fluid velocity results from (94); its components are given by $v_j^{\left(i\right)}=-{\kappa }_j^{\left(i\right)}{\partial }_j{p}^{\left(i\right)}$ , from which we obtain:

$${\mathbf{v}}^{\left(i\right)}=-\frac{3B{\kappa }_1^{\left(i\right)}}{2{\kappa }^{\left(o\right)}+{\kappa }_1}\mathbf{i}-\frac{3D{\kappa }_2^{\left(i\right)}}{2{\kappa }^{\left(o\right)}+{\kappa }_2}\mathbf{j}-\frac{3A{\kappa }_3^{\left(i\right)}}{2{\kappa }^{\left(o\right)}+{\kappa }_3}\mathbf{k}.$$ E97

This is the generalization to the 3D case of the result obtained for the spherical inclusion when the pressure gradient is along the axis Ox (19).

The inner product of ${\mathbf{v}}^{\left(i\right)}$ and of the incident field ${\mathbf{U}}_{\infty }$ , gives the angle $\gamma$ whose the internal fluid velocity is deflected by the anisotropy of the inclusion:

$$\text{cos}\gamma =\frac{{\mathbf{U}}_{\infty }\cdot {\mathbf{v}}^{\left(i\right)}}{\parallel {\mathbf{U}}_{\infty }\parallel \parallel {\mathbf{v}}^{\left(i\right)}\parallel },$$ E98

$$=\frac{\frac{{\text{sin}}^2{\theta }_0{\text{sin}}^2{\phi }_0{\kappa }_1^{\left(i\right)}}{2{\kappa }^{\left(e\right)}+{\kappa }_2^{\left(i\right)}}+\frac{{\text{sin}}^2{\theta }_0{\text{cos}}^2{\phi }_0{\kappa }_2^{\left(i\right)}}{2{\kappa }^{\left(e\right)}+{\kappa }_1^{\left(i\right)}}+\frac{{\text{cos}}^2{\theta }_0{\kappa }_3^{\left(i\right)}}{2{\kappa }^{\left(e\right)}+{\kappa }_2^{\left(i\right)}}}{\sqrt{\frac{{\text{sin}}^2{\theta }_0{\text{sin}}^2{\phi }_0{\kappa }_1^{\left(i\right)}}{{\left(2{\kappa }^{\left(e\right)}+{\kappa }_2^{\left(i\right)}\right)}^2}+\frac{{\text{sin}}^2{\theta }_0{\text{cos}}^2{\phi }_0{\kappa }_2^{\left(i\right)}}{{\left(2{\kappa }^{\left(e\right)}+{\kappa }_1^{\left(i\right)}\right)}^2}+\frac{{\text{cos}}^2{\theta }_0{\kappa }_3^{\left(i\right)}}{{\left(2{\kappa }^{\left(e\right)}+{\kappa }_2^{\left(i\right)}\right)}^2}}}.$$ E99

2.4.2 Isotropic defect in anisotropic porous medium

Consider an isotropic sphere of radius r whose mobility is ${\kappa }^{\left(i\right)}$ which is included in an anisotropic host medium with its own mobility ${\mathit{\kappa }}^{\left(o\right)}$ . The incompressibility of the saturating fluid imposes that the outside pressure is the solution of the equation:

$${\partial }_i\left({\kappa }_{\mathit{ij}}^{\left(o\right)}{\partial }_{j}p\right)=0.$$ E100

In the system of Cartesian coordinate defined by the principal directions of the tensor ${\mathit{\kappa }}^{\left(o\right)}$ , this equation is written as:

$${{\kappa }^{\prime }}^{\left(o\right)}\left(\frac{{\kappa }_x^{\left(o\right)}}{{{\kappa }^{\prime }}^{\left(o\right)}}\frac{{\partial }^2{p}^{\left(o\right)}}{\partial x^2}+\frac{{\kappa }_y^{\left(o\right)}}{{{\kappa }^{\prime }}^{\left(o\right)}}\frac{{\partial }^2{p}^{\left(o\right)}}{\partial y^2}+\frac{{\kappa }_z^{\left(o\right)}}{{{\kappa }^{\prime }}^{\left(o\right)}}\frac{{\partial }^2{p}^{\left(o\right)}}{\partial z^2}\right)=0,$$ E101

where ${\kappa }_j^{\left(o\right)}$ are the eigenvalues of the mobility tensor and ${{\kappa }^{\prime }}^{\left(o\right)}$ is an arbitrary scalar such that the ratio ${\kappa }_{r,j}^{\left(o\right)}={\kappa }_j^{\left(o\right)}/{{\kappa }^{\prime }}^{\left(o\right)}$ is a dimensionless quantity. Using the linear transformation of coordinates:

$$x^{\prime }=\frac{x}{\sqrt{{\kappa }_{r,x}^{\left(o\right)}}},y^{\prime }=\frac{y}{\sqrt{{\kappa }_{r,y}^{\left(o\right)}}},z^{\prime }=\frac{z}{\sqrt{{\kappa }_{r,z}^{\left(o\right)}}},$$ E102

Eq. (101) becomes a Laplace equation. Correspondingly, the sphere is transformed into an ellipsoid with the semiaxes $a_x=r/\sqrt{{\kappa }_{r,x}^{\left(o\right)}}$ , $a_y=r/\sqrt{{\kappa }_{r,y}^{\left(o\right)}}$ , and $a_z=r/\sqrt{{\kappa }_{r,z}^{\left(o\right)}}$ . Since the principal directions of the inside permeability ${\kappa }^{\left(i\right)}$ coincide with the axes of the ellipsoid, for each direction j, we find, for each of the components of the pressure gradient, the result of the ellipsoidal inclusion (40). The internal pressure gradient is then:

$${\partial }_j{p}^{\left(i\right)}=\frac{{{\kappa }^{\prime }}^{\left(o\right)}}{{{\kappa }^{\prime }}^{\left(o\right)}+N_j\left({\kappa }^{\left(i\right)}/{\kappa }_{\mathit{rj}}^{\left(o\right)}-{{\kappa }^{\prime }}^{\left(o\right)}\right)}{\partial }_j{p}^{\left(o\right)},$$ E103

or

$${\partial }_j{p}^{\left(i\right)}=\frac{{\kappa }_j^{\left(o\right)}}{{\kappa }_j^{\left(o\right)}+N_j\left({\kappa }^{\left(i\right)}-{\kappa }_j^{\left(o\right)}\right)}{\partial }_j{p}^{\left(o\right)}.$$ E104

In this equation, the depolarization factor N_j is

$$N_j=\frac{a_x{a}_y{a}_z}{2}{\int }_O^{\infty }\frac{\mathit{ds}}{\left(s+a_j^2\right)\sqrt{\left(s+a_x^2\right)\left(s+a_y^2\right)\left(s+a_z^2\right)}}\text{for}j=x,y,z.$$ E105

Thus, the anisotropy induced in the sphere by the change of variables appears through the depolarization factor N_j.

2.4.3 Anisotropic defect in anisotropic porous medium

We assume now that the host medium and the defect have their own anisotropic microstructure with the mobilities tensors ${\kappa }_{\mathit{ij}}^{\left(o\right)}$ and ${\kappa }_{\mathit{ij}}^{\left(e\right)}$ . The velocity of the fluid flowing in each part of the porous medium is given by equations:

$$v_i^{\left(o\right)}=-{\kappa }_{\mathit{ij}}^{\left(o\right)}{\nabla }_j{p}^{\left(o\right)},v_i^{\left(i\right)}=-{\kappa }_{\mathit{ij}}^{\left(i\right)}{\nabla }_j{p}^{\left(i\right)}.$$ E106

Without restricting the generality of the problem, the first relation of (106) can be written as:

$$v_i^{\left(o\right)}=-{\kappa }_i^{\left(o\right)}\frac{\partial p^{\left(o\right)}}{\partial x_i},$$ E107

where ${\kappa }_j^{\left(o\right)}$ , $j=\mathrm{1,2,3}$ , are the eigenvalues of the tensor ${\mathit{\kappa }}^{\left(\mathit{o}\right)}$ and $v_i^{\left(o\right)}$ and $\frac{\partial p^{\left(o\right)}}{\partial x_i}$ are the components of the velocity and of the pressure gradient along the principal directions of this tensor.

The incompressibility of the fluid implies the condition:

$${\nabla }_i{v}_i^{\left(o\right)}=0,$$ E108

or

$${\kappa }_1^{\left(o\right)}\frac{{\partial }^2{p}^{\left(o\right)}}{\partial x_1^2}+{\kappa }_2^{\left(o\right)}\frac{{\partial }^2{p}^{\left(o\right)}}{\partial x_2^2}+{\kappa }_3^{\left(o\right)}\frac{{\partial }^2{p}^{\left(o\right)}}{\partial x_3^2}=0.$$ E109

To transform this equation into a Laplace equation, we proceed as before by using the change of variables

$$\frac{\partial }{\partial x_i^{\prime }}=\sqrt{{\kappa }_{r,i}^{\left(o\right)}}\frac{\partial }{\partial x_i}.$$ E110

Then, the external environment becomes an isotropic medium and the outside pressure is the solution of the Laplace equation:

$$\frac{\partial p^{\left(o\right)}}{\partial {x^{\prime }}_1^2}+\frac{\partial p^{\left(o\right)}}{\partial {x^{\prime }}_2^2}+\frac{{\partial }^2{p}^{\left(o\right)}}{\partial {x^{\prime }}_3^2}=0.$$ E111

The new $x_i^{\prime }$ variables constitute a new coordinate system. The host medium is transformed into an isotropic medium, while the inclusion medium becomes anisotropic. In the new coordinate system, the pressure gradient is transformed according to:

$${\mathbf{\text{∇}}}^{\prime }p=\sqrt{{\kappa }_{r,j}^{\left(o\right)}}\mathbf{\text{∇}}p,$$ E112

while the components of the position vector become:

$$r_i^{\prime }={\left({\kappa }_{r,i}^{\left(o\right)}\right)}^{-1/2}{r}_i.$$ E113

Using the Darcy’s law and (112), the incompressibility of the fluid inside the inclusion

$${\nabla }_i\left({\kappa }_{\mathit{ij}}^{\left(i\right)}{\nabla }_j{p}^{\left(i\right)}\right)=0$$ E114

implies the new equation:

$${\nabla }_i^{\prime }\left[{\left({\kappa }_i^{\left(o\right)}\right)}^{-1/2}{\kappa }_{\mathit{ij}}^{\left(i\right)}{\left({\kappa }_j^{\left(o\right)}\right)}^{-1/2}{\nabla }_j^{\prime }{p}^{\left(i\right)}\right]=0.$$ E115

In the new coordinates system, the mobility in the inclusion defined by the equation:

$${\nabla }_i^{\prime }\left({\kappa }_{\mathit{ij}}^{{\left(i\right)}^{\prime }}{\nabla }_i^{\prime }{p}_j^{\left(i\right)}\right)=0$$ E116

is such that:

$${\kappa }_{\mathit{ij}}^{\left(i\right)\text{'}}={\left({\kappa }_{r,i}^{\left(o\right)}\right)}^{-1/2}{\kappa }_{\mathit{ij}}^{\left(i\right)}{\left({\kappa }_{r,j}^{\left(o\right)}\right)}^{-1/2}.$$ E117

In the coordinates $x_i^{\prime }$ , the semiaxes of the inclusion can be calculated from the equation of the ellipsoidal inclusion surface written in matrix form as $R^t\mathit{AR}$ , where R is the position vector of a point of this surface ( $R^t=\left(x,y,z\right)$ ) and A is the diagonal matrix whose entries are the lengths of the half-axes:

$$A=\left(\begin{array}{lll}{a}_1& 0& 0\\ 0& a_2& 0\\ 0& 0& a_3\end{array}\right).$$ E118

Then, the linear transformation (113) changes A into $A^{\prime }$ :

$$A^{\prime }=\left(\begin{array}{lll}{a}_1^{\prime }& 0& 0\\ 0& a_2^{\prime }& 0\\ 0& 0& a_3^{\prime }\end{array}\right),$$ E119

with

$${a^{\prime }}_i^2=={\left({\kappa }_{r,i}^{\left(o\right)}\right)}^{-1/2}{a}_i^2{\left({\kappa }_{r,i}^{\left(o\right)}\right)}^{-1/2}.$$ E120

Thus, the operation that transforms the anisotropic host medium into an isotropic one transforms the ellipsoidal inclusion with the semiaxes (a₁, a₂, a₃) into another one with the new semiaxes ( $a_1^{\prime }$ , $a_2^{\prime }$ , $a_3^{\prime }$ ) and the new mobility ${\kappa }_{\mathit{ij}}^{{\left(i\right)}^{\prime }}$ given respectively by (120) and.

We recover the previous case where the outer medium is isotropic and the inner medium is anisotropic. So, in accordance with (104):

$${\partial }_j{p}^{\left(i\right)}=\frac{{{\kappa }^{\prime }}^{\left(o\right)}}{{{\kappa }^{\prime }}^{\left(o\right)}+N_j^{\prime }\left({\kappa }^{\left(i\right)}/{\kappa }_{\mathit{rj}}^{\left(o\right)}-{{\kappa }^{\prime }}^{\left(o\right)}\right)}{\partial }_j{p}^{\left(o\right)},$$ E121

where the depolarization factors of the new inclusion are given by:

$$N_i^{\prime }=\frac{\mathit{det}{A}^{\prime }}{2}{\int }_0^{\infty }\frac{\mathit{d\sigma }}{\left(a_i^{\prime }+\sigma \right)\sqrt{\mathit{det}\left({A^{\prime }}^2+\mathit{\sigma I}\right)}}.$$ E122

2.5 Hydraulic polarisability

As mentioned above, the reaction of a saturated porous inclusion subject to a pressure gradient is to induce a hydraulic dipole whose dipole moment is $\mathcal{P}$ . This dipole results from the appearance of pressure discontinuities at the inclusion-host interface. They have different signs depending on whether the flow is incoming or outgoing, but have the same absolute value. They are the hydraulic analogues of electrostatic charges induced by an electric field in a dielectric medium. The resulting hydraulic polarization is only nonzero if the contrast between the mobility of the host environment and that of inclusion is itself nonzero.

$$\mathcal{P}={\int }_{\mathrm{\Omega }}\left({\kappa }^{\left(i\right)}\left(\mathbf{r}\right)-{\kappa }^{\left(o\right)}\right){\kappa }^{\left(i\right)}\mathbf{\text{∇}}{p}^{\left(i\right)}\left(\mathbf{r}\right)\mathit{dV}.$$ E123

For spherical or ellipsoidal inclusions and for low filtration rates, we have seen that the internal pressure gradient is proportional to the incident one. For this type of inclusions, the dipole moment is written as:

$$\mathcal{P}=\alpha {\mathbf{v}}^{\left(o\right)},$$ E124

where the value of the susceptibility α measures the ability of the inclusion to induce a dipole under the action of a pressure gradient. α can be seen as the “hydraulic polarisability” of the defect. For a spherical defect of volume V, we have:

$$\mathcal{P}=\alpha \left(-{\mathbf{U}}_{\infty }\right)$$ E125

with

$$\alpha =3V{\kappa }^{\left(i\right)}\left(\frac{\frac{k^{\left(i\right)}}{k^{\left(o\right)}}-1}{2+\frac{k^{\left(i\right)}}{k^{\left(o\right)}}}\right).$$ E126

For an ellipsoidal inclusion, hydraulic polarisability is not a scalar since the response of the inclusion is a function of the direction of pressure incidence, but a second rank tensor whose eigenvalues are the susceptibilities along the three axes of the ellipsoid:

$${\alpha }_i=V{\kappa }^{\left(i\right)}\frac{{\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}}{{\kappa }^{\left(o\right)}+N_i\left({\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}\right)},i=x,y,z.$$ E127

Figure 11 represents the surface hydraulic charges induced by the hydraulic polarization on a spherical inclusion and on ellipsoidal inclusions with different orientations with respect to the direction of the incident flux. The red areas represent the surface “hydraulic charge” density ${\sigma }_{\mathit{pol}}$ . Its expression depends on the direction of the incident pressure gradient and is written as the sum of the contributions of the dipoles along the three axes of the ellipsoid.

Figure 11.

Sketch of polarization surface “charge” density ${\sigma }_{\mathit{pol}}$ for a spherical inclination and an ellipsoidal inclusion with different orientations with respect to the incident flow.

When the permeability of the inclusion is stratified, the dipole moment is given by the dipolar term ( $B^{\left(o\right)}$ ) of the external pressure field obtained by the transfer matrix method or by the scattering matrix method.

3. Tortuosity induced by defects

In this section, we determine the hydraulic effects of defects on the permeability of porous media. As mentioned above, the shape of the defects is one of the most important factors for the modification of the current lines of the seepage rate in the whole porous medium and thus contributes to its acoustic properties.

3.1 Homogenization: generalities

Experiments show that a nonhomogeneous medium subject to excitation behaves in the same way as its different components, but with different parameter values. The homogenization of an inhomogeneous porous medium consists in replacing it with an effective homogeneous medium with the permeability k_eff. This operation is only possible at a fairly large observation scale. Determining the value of the effective permeability from the mobility values of the structure components and their relative positions is not a simple averaging operation. The calculation of the global mobility of a mixture of porous inclusions immersed in a homogeneous medium is a topic widely addressed in many research fields such as hydrology, oil recovery, chemical industry, etc. As a consequence, a considerable number of works deal with this problem based on various methods: renormalization theory, variational methods, T-Matrix method, field theory methods, nonperturbative approach based on Feynman path integral. To quote some of authors, we can refer to the works of Prakash and Raja-Sekhar [14], King [15, 16], Drummond and Horgan [17], Dzhabrailov and Meilanov [18], Teodorovich [19, 20], Stepanyants and Teodorovich [21], and Hristopulos and Christakos [22].

In the case of media subject to a variable field action, homogenization requires defining a length below which it is no longer relevant. For example, for a periodic field, acting on a medium whose average distance between inhomogeneities is a, this length is the wavelength $\lambda$ if $\lambda /a\gg 1$ . In our case, effective mobility being essentially a low-frequency concept, this remark justifies that the effective mobility should then be calculated from a steady filtration velocity.

3.2 Effective mobility

Darcy’s law is often used as the definition of the mobility of a porous medium, and the easiest way to introduce the effective mobility ${\kappa }_{\mathit{eff}}$ is to use it as follows:

$$<\mathbf{v}>={\kappa }_{\mathit{eff}}<\mathbf{E}>,$$ E128

where $\mathbf{E}=-\mathbf{\text{∇}}p$ and $<\cdot >$ is the averaging operation. The mean values of the filtration rate and the pressure gradient are given by:

$$<\mathbf{v}>=f<{\kappa }^{\left(i\right)}{\mathbf{E}}^{\left(i\right)}>+\left(1-f\right)<{\kappa }^{\left(o\right)}{\mathbf{E}}^{\left(o\right)}>,$$ E129

$$<\mathbf{E}>=f<{\mathbf{E}}^{\left(i\right)}>+\left(1-f\right)<{\mathbf{E}}^{\left(o\right)}>,$$ E130

where f is the volumic fraction of the defect. Putting ${\mathbf{E}}^{\left(i\right)}=\mathcal{A}{\mathbf{E}}^{\left(o\right)}$ we show that:

$${\kappa }_{\mathit{eff}}=\frac{f\mathcal{A}{\kappa }^{\left(i\right)}+\left(1-f\right){\kappa }^{\left(o\right)}}{f\mathcal{A}+\left(1-f\right)}.$$ E131

For spherical defects, $\mathcal{A}=3{\kappa }^{\left(o\right)}/(2{\kappa }^{\left(o\right)}+{\kappa }^{\left(i\right)}$ , Eq. (131) leads to the result:

$${\kappa }_{\mathit{eff}}={\kappa }^{\left(o\right)}\frac{f\frac{3{\kappa }^{\left(i\right)}}{2{\kappa }^{\left(o\right)}+{\kappa }^{\left(i\right)}}+\left(1-f\right)}{f\frac{3{\kappa }^{\left(o\right)}}{2{\kappa }^{\left(o\right)}+{\kappa }^{\left(i\right)}}+\left(1-f\right)}.$$ E132

When $f\to 0$ , then ${\kappa }_{\mathit{eff}}\sim {\kappa }^{\left(o\right)}$ , and when $f\to 1$ , then ${\kappa }_{\mathit{eff}}\sim {\kappa }^{\left(i\right)}$ . Finally when $f<<1$ , then

$${\kappa }_{\mathit{eff}}\sim {\kappa }^{\left(o\right)}+3f{\kappa }^{\left(o\right)}\frac{{\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}}{{\kappa }^{\left(i\right)}+2{\kappa }^{\left(o\right)}}.$$ E133

For anisotropic inclusion, mobility is a second rank tensor defined by the relationship:

$$<v_i>={\kappa }_{\mathit{eff},\mathit{ij}}<E_j>.$$ E134

As a result, for each main direction, we have:

$$<v_j>={\kappa }_{\mathit{eff},j}<E_j>,j=x,y,z,$$ E135

where ${\kappa }_{\mathit{eff},j}$ are the eigenvalues of ${\kappa }_{\mathit{eff}}$ . Taking into account that $E_j^{\left(i\right)}={\mathcal{A}}_j{E}_j^{\left(o\right)}$ , Eq. (135) shows that:

$${\kappa }_{\mathit{eff},j}=\frac{f{\mathcal{A}}_j{\kappa }^{\left(i\right)}+\left(1-f\right){\kappa }^{\left(o\right)}}{f{\mathcal{A}}_j+\left(1-f\right)},$$ E136

with

$${\mathcal{A}}_j=\frac{3{\kappa }^{\left(o\right)}}{2{\kappa }^{\left(o\right)}+{\kappa }_j^{\left(i\right)}},j=\mathrm{1,2,3}.$$ E137

When the environment has several defects, the calculation of k_eff is more complicated because their mutual influence must be taken into account. The excitation pressure gradient ${\mathbf{E}}^e$ defined from the filtration rate is introduced by the equation:

$${\kappa }^{\left(o\right)}{\mathbf{E}}^e={\kappa }^{\left(o\right)}\mathbf{E}+L\cdot P.$$ E138

In this relationship, L is an operator that takes into account the shape of the defect and its orientation with respect to the fluid flow, and P is due to the induced polarization in the inclusion. P defined by (123) is related to the dipole moment induced by the interaction between the fluid moving in the porous medium and the defect. When the medium contains n identical defects per unit volume, $P=\mathit{np}$ , p being the dipole moment of each defect, and since p is proportional to the applied field ( $p=\alpha {\mathbf{E}}^e$ ), we have $P=\mathit{n\alpha }{\mathbf{E}}^e$ . For an ellipsoidal defect, L is reduced to depolarization factors, i.e., $L=N_k$ , $k=x,y,z$ , which takes into account the direction of fluid flow. In this case, the excitation field is:

$$E^e=\frac{E}{1-\mathit{n\alpha }{N}_k/{\kappa }^{\left(o\right)}},$$ E139

leading to the following ${\kappa }_{\mathit{eff}}$ expression:

$${\kappa }_{\mathit{eff}}={\kappa }^{\left(o\right)}+\frac{\mathit{n\alpha }}{1-\mathit{n\alpha }{N}_k/{\kappa }^{\left(o\right)}},$$ E140

leading, for spherical defects, to expression:

$${\kappa }_{\mathit{eff}}={\kappa }^{\left(o\right)}+\frac{\mathit{n\alpha }}{1-\mathit{n\alpha }/3{\kappa }^{\left(o\right)}}.$$ E141

It is then possible to calculate the effective mobility of a set of ellipsoidal inclusions in different geometries (Figure 12):

• the ellipsoids are aligned with the direction x of the fluid flow:

Figure 12.

Sets of ellipsoidal inclusions: (a) aligned and (b) randomly oriented.

$${\kappa }_{\mathit{eff}}={\kappa }^{\left(o\right)}+\frac{n{\alpha }_x}{1-n{\alpha }_x{N}_x/{\kappa }^{\left(o\right)}}.$$ E142

• the ellipsoids are randomly oriented:

$${\kappa }_{\mathit{eff}}={\kappa }^{\left(o\right)}+\frac{1/3\sum _{i=x,y,z}n{\alpha }_i}{1-\sum _{i=x,y,z}n{\alpha }_i{N}_i/{\kappa }^{\left(o\right)}}.$$ E143

3.3 Induced tortuosity

We restrict ourselves here to the calculation of the tortuosity induced by homogeneous spherical inclusions (Figure 13). Since the dipole moment is the essential element for this calculation element for this calculation, it is easy to generalize the results obtained with other types of inclusions: inhomogeneous spherical, ellipsoidal, etc.

Figure 13.

Domains used in numerical simulations for the evaluation of $<v^{\left(o\right)2}>$ .

The tortuosity induced by the presence of defects ${\tau }_d$ is defined by:

$${\tau }_d=\frac{<v^{\left(o\right)2}>}{<v^{\left(o\right)}{>}^2}$$ E144

where $v^{\left(o\right)}$ is the perturbation of filtration rate due to the defects. When the ratio $k/a^2$ where a is the characteristic size of the defects is small relative to the unity, it is legitimate to neglect the volume of the defects for the calculation of $<v^2>$ , whereas it is taken into account for that of $<v>$ .

With the pressure scattered field by the inclusions being limited to the dipolar terms, the expression of $<v^{\left(o\right)2}>$ is then:

$$<v^{\left(o\right)2}>=<v_r^{\left(o\right)2}>+<v_{\theta }^{\left(o\right)2}>$$ E145

where

$$v_r^{\left(o\right)}=-{\kappa }^{\left(o\right)}\frac{\partial p^{\left(o\right)}}{\partial r},\text{and}{v}_{\theta }^{\left(o\right)}=-{\kappa }^{\left(o\right)}\frac{\partial p^{\left(o\right)}}{r\partial \theta }.$$ E146

For ellipsoidal inclusion, the external pressure is:

$$p^{\left(o\right)}\left(r,\theta \right)=-\frac{U_{\infty }}{{\kappa }^{\left(o\right)}}r\text{cos}\theta +\frac{{\mathcal{P}}_d}{4\pi }\frac{\text{cos}\theta }{r^2}$$ E147

$$=-\frac{U_{\infty }}{{\kappa }^{\left(o\right)}}\left(r-\frac{\alpha }{r^2}\right)\text{cos}\theta$$ E148

where ${\mathcal{P}}_d$ and $\alpha$ are, respectively, the dipol moment and the polarisability of the inclusion. By keeping only the terms greater than or equal to $r^{-2}$ , one obtains:

$$v^{\left(o\right)2}\approx {\frac{U_{\infty }}{{\kappa }^{\left(o\right)}}}^2\left[1-\frac{2\alpha }{r^3}\left(2{\text{cos}}^2\theta -{\text{sin}}^2\theta \right)\right]$$ E149

The average value $<v^{\left(o\right)2}>$ is calculated by integration on the volume between two spheres of radius a (characteristic size of the defect) and R sufficiently large so that the dipolar effects are negligible. For a spherical inclusion, it results:

$$<v^{\left(o\right)2}>=U_{\infty }^2+U_{\infty }^2\left(1-\frac{a^3}{R^3}\right){\left(\frac{\frac{{\kappa }^{\left(i\right)}}{{\kappa }^{\left(o\right)}}-1}{2+\frac{{\kappa }^{\left(i\right)}}{{\kappa }^{\left(o\right)}}}\right)}^2.$$ E150

$<v^{\left(o\right)}{>}^2$ is calculated from the definition of effective mobility:

$$<v>={\kappa }_{\mathit{eff}}<E>$$ E151

where ${\kappa }_{\mathit{eff}}$ and $<E>$ are given by (131) and (130). When $f<<1$ , we have

$${\kappa }_{\mathit{eff}}\approx {\kappa }^{\left(o\right)}+f\mathcal{A}\left({\kappa }^{\left(i\right)}-{\kappa }^{\left(o\right)}\right)\text{and}<E>\approx \frac{U_{\infty }}{{\kappa }^{\left(o\right)}}\left(1+f\mathcal{A}\right).$$ E152

From these two relations, we obtain the expression of the induced tortuosity:

$${\tau }_d=\left[1+\left(1-\frac{a^3}{R^3}\right){\left(\frac{\frac{{\kappa }^{\left(i\right)}}{{\kappa }^{\left(o\right)}}-1}{2+\frac{{\kappa }^{\left(i\right)}}{{\kappa }^{\left(o\right)}}}\right)}^2\right]\left[1-2f\frac{3\frac{{\kappa }^{\left(i\right)}}{{\kappa }^{\left(o\right)}}}{2+\frac{{\kappa }^{\left(i\right)}}{{\kappa }^{\left(o\right)}}}\right].$$ E153

Results of numerical simulations: The results of a numerical simulation for ${\kappa }^{\left(i\right)}/{\kappa }^{\left(o\right)}=10$ and ${\kappa }^{\left(i\right)}/{\kappa }^{\left(o\right)}=100$ are shown in Figure 14. The tortuosity value is calculated on square domains around the inclusion (Figure 13). Inside the inclusion, ${\tau }_b=1$ . As x increases, the tortuosity increases to reach its maximum value at $x=1.7$ for ${\kappa }^{\left(i\right)}/{\kappa }^{\left(o\right)}=10$ and $x=1.6$ when ${\kappa }^{\left(i\right)}/{\kappa }^{\left(o\right)}=10$ . For larger values of x, it decreases toward 1 since, far from inclusion, the field lines again become parallel to the direction of the incident pressure gradient. This result confirms the behavior of the field lines of Figure 4b.

Figure 14.

Evaluation of the tortuosity induced as a function of the distance to the defect.

4. Conclusion

In this chapter, we studied the effect of defects on the circulation of the fluid saturating a porous medium. We have shown that the modification of the stream lines of the filtration velocities leads to a modification of the value of the tortuosity and thus on the local velocity of the waves susceptible to propagate in such media. The induced tortuosity was calculated from the pressure field scattered by the inclusions. The model used is based on the Darcy’s law. in addition to being general, its major interest is to lead to a very practical mathematical expression of tortuosity

Acknowledgments

C. Depollier is supported by Russian Science Foundation grant number 14-49-00079.

The ellipsoidal coordinates ( $\xi$ , $\eta$ , $\zeta$ ) are the solutions of the cubic equation:

$$\frac{x^2}{a^2+u}+\frac{y^2}{b^2+u}+\frac{z^2}{c^2+u}=1.$$ E154

They are connected to the Cartesian coordinates ( $x$ , $y$ , $z$ ) by the relations:

$$x^2=\frac{\left(a^2+\xi \right)\left(a^2+\eta \right)\left(a^2+\zeta \right)}{\left(b^2-a^2\right)\left(c^2-a^2\right)},$$ E155

$$y^2=\frac{\left(b^2+\xi \right)\left(b^2+\eta \right)\left(b^2+\zeta \right)}{\left(a^2-b^2\right)\left(c^2-b^2\right)},$$ E156

$$z^2=\frac{\left(c^2+\xi \right)\left(c^2+\eta \right)\left(c^2+\zeta \right)}{\left(a^2-c^2\right)\left(b^2-c^2\right)},$$ E157

subject to the conditions $-\xi <c^2<-\eta <b^2<-\zeta <a^2$ .

The scalar factors are the vector norms:

$$h_{q_i}=\parallel \frac{\partial \mathbf{r}}{\partial q_i}\parallel \mathrm{o}\grave{\mathrm{u}}{q}_i=\xi ,\eta ,\zeta .$$ E158

Their values are:

$$h_{\xi }=\frac{1}{2}\sqrt{\frac{\left(\eta -\xi \right)\left(\zeta -\xi \right)}{\left(a^2-\xi \right)\left(b^2-\xi \right)\left(c^2-\xi \right)}},$$ E159

$$h_{\eta }=\frac{1}{2}\sqrt{\frac{\left(\xi -\eta \right)\left(\zeta -\eta \right)}{\left(a^2-\eta \right)\left(b^2-\eta \right)\left(c^2-\eta \right)}},$$ E160

$$h_{\zeta }=\frac{1}{2}\sqrt{\frac{\left(\eta -\zeta \right)\left(\xi -\zeta \right)}{\left(a^2-\zeta \right)\left(b^2-\zeta \right)\left(c^2-\zeta \right)}}.$$ E161

The depolarization factors are important quantities for the expression of solutions of the Laplace equation. They take into account the form of the domain in which this solution is sought and its orientation in relation to the excitation field. Their expression is:

$$N_k=\frac{\mathit{abc}}{3}{\int }_0^{\infty }\frac{\mathit{d\sigma }}{\left(\sigma +q_k^2\right)\sqrt{\left(\sigma +a^2\right)\left(\sigma +b^2\right)\left(\sigma +c^2\right)}}$$ E162

where $k=x$ (resp. y, z) et $q_k=a$ , (resp. b, c) and satisfy the relation:

$$N_x+N_y+N_z=1.$$ E163

Consider the rectangular coordinate systems ( $x$ , $y$ , $z$ ) and ( $x^{\prime }$ , $y^{\prime }$ , $z^{\prime }$ ). We are looking for the relations between the spherical coordinates ( $r$ , $\theta$ , $\phi$ ) and ( $r^{\prime }$ , ${\theta }^{\prime }$ , ${\phi }^{\prime }$ ) associated with each of them. From

$$x^{\prime }=r^{\prime }\text{sin}{\theta }^{\prime }\text{cos}{\phi }^{\prime }x=r\text{sin}\theta \text{cos}\phi ,$$ E164

$$y^{\prime }=r^{\prime }\text{sin}{\theta }^{\prime }\text{sin}{\phi }^{\prime }y=r\text{sin}\theta \text{sin}\phi ,$$ E165

$$z^{\prime }=r^{\prime }\text{cos}{\theta }^{\prime }z=r\text{sin}\theta ,$$ E166

one deduces

$$\begin{array}{cc}{r^{\prime }}^2& ={x^{\prime }}^2+{y^{\prime }}^2+{z^{\prime }}^2\\ & =\frac{x^2}{{\kappa }_1}+\frac{y^2}{{\kappa }_2}+\frac{z^2}{{\kappa }_3}\\ & =r^2\left(\frac{{\text{sin}}^2\theta {\text{cos}}^2\phi }{{\kappa }_1}+\frac{{\text{sin}}^2\theta {\text{sin}}^2\phi }{{\kappa }_2}+\frac{{\text{cos}}^2\theta }{{\kappa }_3}\right)\end{array}$$

or

$$r^{\prime }=r\mathrm{\Delta }\mathrm{o}\grave{\mathrm{u}}\mathrm{\Delta }=\sqrt{\frac{{\text{sin}}^2\theta {\text{cos}}^2\phi }{{\kappa }_1}+\frac{{\text{sin}}^2\theta {\text{sin}}^2\phi }{{\kappa }_2}+\frac{{\text{cos}}^2\theta }{{\kappa }_3}}.$$ E167

From (164) and (165), one has:

$$\begin{array}{l}{x}^{\prime }=\frac{x}{{\kappa }_1}\Rightarrow r^{\prime }\text{sin}{\theta }^{\prime }\text{cos}{\phi }^{\prime }=\frac{r}{{\kappa }_1}\text{sin}\theta \text{cos}\phi \\ y^{\prime }=\frac{y}{{\kappa }_2}\Rightarrow r^{\prime }\text{sin}{\theta }^{\prime }\text{sin}{\phi }^{\prime }=\frac{r}{{\kappa }_2}\text{sin}\theta \text{sin}\phi .\end{array}$$

By eliminating $\phi$ , one finds:

$$\text{sin}{\theta }^{\prime }=\text{sin}\theta \frac{\delta }{\mathrm{\Delta }}\text{avec}\delta =\sqrt{\frac{{\text{cos}}^2\phi }{{\kappa }_1}+\frac{{\text{sin}}^2\phi }{{\kappa }_2}}.$$ E168

In the same way, from (166), one can establish

$$\text{cos}{\theta }^{\prime }=\frac{1}{\mathrm{\Delta }}\frac{\text{cos}\theta }{\sqrt{{\kappa }_3}}.$$ E169

Similar relationships between angles $\phi$ and ${\phi }^{\prime }$ are deduced from (164) and (165):

$$\text{sin}{\phi }^{\prime }=\frac{1}{\delta }\frac{\text{sin}\phi }{\sqrt{{\kappa }_2}},$$ E170

$$\text{cos}{\phi }^{\prime }=\frac{1}{\delta }\frac{\text{cos}\phi }{\sqrt{{\kappa }_1}}.$$ E171

D.

Let $\mathit{\kappa }$ be a tensor of rank 2. We denote by ${\mathit{\kappa }}_r$ , its expression in the system of rectangular coordinates defined by its principal directions, and ${\mathit{\kappa }}_s$ , its expression in the corresponding spherical coordinates system. So we have

$${\mathit{\kappa }}_r=\left(\begin{array}{lll}{\kappa }_1& 0& 0\\ 0& {\kappa }_2& 0\\ 0& 0& {\kappa }_3\end{array}\right){\mathit{\kappa }}_s=\left(\begin{array}{lll}{\kappa }_{11}& {\kappa }_{12}& {\kappa }_{13}\\ {\kappa }_{21}& {\kappa }_{22}& {\kappa }_{23}\\ {\kappa }_{31}& {\kappa }_{32}& {\kappa }_{33}\end{array}\right)$$ E172

with

$${\kappa }_{11}={\kappa }_1{\text{sin}}^2\theta {\text{cos}}^2\phi +{\kappa }_2{\text{sin}}^2\theta {\text{sin}}^2\phi +{\kappa }_3{\text{cos}}^2\theta ,$$ E173

$${\kappa }_{12}={\kappa }_1\text{cos}\theta \text{sin}\theta {\text{cos}}^2\phi +{\kappa }_2\text{cos}\theta \text{sin}\theta {\text{sin}}^2\phi -{\kappa }_3\text{cos}\theta \text{sin}\theta ,$$ E174

$${\kappa }_{13}=\left({\kappa }_2-{\kappa }_1\right)\text{sin}\theta \text{cos}\phi \text{sin}\phi ,$$ E175

$${\kappa }_{22}={\kappa }_1{\text{cos}}^2\theta {\text{cos}}^2\phi +{\kappa }_2{\text{cos}}^2\theta {\text{sin}}^2\phi +{\kappa }_3{\text{sin}}^2\theta ,$$ E176

$${\kappa }_{23}=\left({\kappa }_2-{\kappa }_1\right)\text{cos}\theta \text{cos}\phi \text{sin}\phi ,$$ E177

$${\kappa }_{33}={\kappa }_1{\text{sin}}^2\phi +{\kappa }_2{\text{cos}}^2\phi ,$$ E178

$${\kappa }_{21}={\kappa }_{12},$$ E179

$${\kappa }_{31}={\kappa }_{13},$$ E180

$${\kappa }_{32}={\kappa }_{23}.$$ E181

Or, alternatively in the matrix form:

$$\mathit{k}={\kappa }_1\mathbf{I}+\left({\kappa }_2-{\kappa }_1\right)\mathbf{A}+\left({\kappa }_3-{\kappa }_1\right)\mathbf{B},$$ E182

where $\mathbf{I}$ is the unit matrix $3\times 3$ and $\mathbf{A}$ and $\mathbf{B}$ are given by:

$$\mathbf{A}=\left(\begin{array}{lll}{\text{sin}}^2\theta {\mathit{sin}}^2\phi & \text{cos}\theta \text{sin}\theta {\text{sin}}^2\phi & \text{sin}\theta \text{cos}\phi \text{sin}\phi \\ \text{cos}\theta \text{sin}\theta {\text{sin}}^2\phi & {\text{cos}}^2\theta {\text{sin}}^2\phi & \text{cos}\theta \text{cos}\phi \text{sin}\phi \\ \text{sin}\theta \text{cos}\phi \text{sin}\phi & \text{cos}\theta \text{cos}\phi \text{sin}\phi & {\text{cos}}^2\theta \end{array}\right),$$ E183

$$\mathbf{B}=\left(\begin{array}{lll}{\text{cos}}^2\theta & -\text{cos}\theta \text{sin}\theta & 0\\ -\text{cos}\theta \text{sin}\theta & {\text{sin}}^2\theta & 0\\ 0& 0& 0\end{array}\right).$$ E184