## Abstract

The complex variable method and conformal mapping are used to derive a closed-form expression for the compressibility of an isolated pore in an infinite two-dimensional, isotropic elastic body. The pore is assumed to have an *n*-fold axis of symmetry, and be represented by at most four terms in the mapping function that conformally maps the exterior of the pore into the interior of the unit circle. The results are validated against some special cases available in the literature, and against boundary-element calculations. By extrapolation of the results for pores obtained from three and four terms of the Schwarz–Christoffel mapping function for regular polygons, the compressibilities of a triangle, square, pentagon and hexagon are found (to at least three digits). Specific results for some other pore shapes, more general than the quasi-polygons obtained from the Schwarz–Christoffel mapping, are also presented. An approximate scaling law for the compressibility, in terms of the ratio of perimeter-squared to area, is also tested. This expression gives a reasonable approximation to the pore compressibility, but may overestimate it by as much as 20%.

## 1. Introduction

Modelling the effective elastic properties of porous media is an exercise that has important implications for rocks (Zimmerman 1991), ceramics (Rice 1998), and bone (Cowin 2001), as well as other materials, both natural and man-made. This problem can be considered to consist of two aspect: (i) estimating the intrinsic, shape-dependent effect of a single, isolated pore, and (ii) the development of methods, either analytical or numerical, to estimate the effect of finite concentrations of pores. Most treatments of this problem assume that the pores can be represented as ellipsoids or spheroids, or some degenerate special cases thereof, such as spheres, cylinders, or ‘penny-shaped’ cracks. For the entire family of ellipsoids, the solution to the first problem is essentially contained in the works of Eshelby (1957) and Wu (1966). Their results lead to exact expressions for the effective moduli, to first-order in porosity, that are consequently accurate for small values of porosity. Numerous schemes for extending these results to higher pore concentrations have been developed, although no single scheme is as yet universally accepted as being preferable for all cases (Hashin 1983; Christensen 1991; Nemat-Nasser & Hori 1999).

Although ellipsoidal pore models offer the advantage of mathematical tractability, examination of images of real materials reveals that pores are in many cases more irregular. Pores in ceramics often seem, at least visually in two-dimensional sections, to be better modelled as being polygonal. Pores in sedimentary rocks are much more irregular, and are often not convex. Analysis of the influence that such irregular three-dimensional pores would have on the elastic moduli is thus far intractable analytically, and only slight progress has been made numerically (Burnley & Davis 2004). But if the pores are modelled as two-dimensional, pores of essentially any shape can be analysed, using the complex variable methods developed by Muskhelishvili (1963) and others. These complex variable methods are indeed quite powerful, and have been used not only for linear elastic problems involving cavities, but can also be used, for example, for elastic inclusions (Chang & Conway 1968; Jasiuk 1995), problems involving finite elasticity (Ru *et al*. 2004) and thermopiezoelectric problems involving cavities (Qin *et al*. 1999).

The complex variable method involves mapping the interior of the unit circle into the exterior of the hole in the physical domain. Using this approach, Savin (1961) analysed several holes of polygonal and quasi-polygonal shape, but with emphasis on calculating the stress concentrations, rather than on the displacements or on the effect that the pores have on the macroscopic elastic moduli. With regards to the effect of pores on the macroscopic elastic properties, some results for quasi-polygons represented by three and four terms of the Schwarz–Christoffel mapping function have been obtained by Zimmerman (1986), Jasiuk *et al*. (1994) and Kachanov *et al*. (1994). In the present work, we develop explicit expressions for holes that have an *n*-fold axis of symmetry, and which can be represented by up to four terms in the mapping function. Our results are more general than those of Jasiuk *et al*. (1994) and Kachanov *et al*. (1994), in that we give explicit results for arbitrary (allowable) values of the mapping coefficients, rather than only for the case of the quasi-polygons represented by the Schwarz–Christoffel mapping. However, unlike Jasiuk *et al*. (1994) and Kachanov *et al*. (1994), we consider only the case of hydrostatic loading, and do not consider shear loading.

Specifically, we calculate the ‘pore compressibility’ (of an isolated pore), *C*_{pp}, which is the fractional change in pore area due to a hydrostatic pressure of unit magnitude acting along the pore walls. This parameter appears directly in equations such as the fluid diffusion equation for porous media, where it is added to the fluid compressibility in the denominator of the hydraulic diffusivity term (Bear 1967). But it is also directly related to the effective macroscopic compressibility, by , where *C*_{o} is the two-dimensional compressibility of the non-porous host material, and *ϕ* is the porosity (Zimmerman 1991). Hence, our calculations are related, more directly than might first appear, to the vast literature on effective elastic moduli.

## 2. Formulation of the problem

The stresses and displacements under two-dimensional plane strain or plane stress conditions can be represented in terms of two complex potential functions, *ϕ*(*z*) and *ψ*(*z*) as follows (Godfrey 1959; Muskhelishvili 1963; England 1971):(2.1)where *G* is the shear modulus of the host material, and *κ*=3−4*ν* for plane strain and (3−*ν*)/(1+*ν*) for plane stress, where *ν* is the Poisson ratio. It is sometimes convenient to refer to the ‘complex displacement’, *U*≡*u*_{x}+i*u*_{y}. Consider an infinite elastic body containing a single, simple closed contour *Γ*, with no stresses acting at infinity, and with a uniform hydrostatic pressure of magnitude *p* acting along *Γ*. For problems in which the tractions are specified along the contour *Γ*, the boundary condition for these potentials along *Γ* can be written as(2.2)where *F* is equal to *i* times the integral of the complex traction vector, *t*_{x}+i*t*_{y}, along the boundary contour, starting from some arbitrary point *z*_{0} on *Γ*.

The solution of this problem using the complex potential method proceeds as follows. First, the region outside of *Γ* in the *z*-plane is mapped into the interior of the unit circle *γ* in the *ζ*-plane by a conformal mapping of the form (figure 1)(2.3)If only two terms in the mapping are taken, i.e. , the hole is a hypotrochoid, which is a quasi-polygon having *n*+1 equal ‘sides’. In order for the mapping to be single-valued, and for *Γ* not to contain any self-intersections, *a*_{n} must satisfy the restriction 0≤*a*_{n}≤1/*n*. The choice *a*_{n}=0 gives a circle, whereas *a*_{n}=1/*n* gives a pore with *n*+1 pointed cusps. For the particular choice , the mapping coincides with the first two terms of the Schwarz–Christoffel mapping for an (*n*+1)-sided equilateral polygon, and resembles a polygon with slightly rounded corners.

If the hole contour possesses an (*n*+1)-fold axis of symmetry, only powers that differ by (*n*+1) will appear in the mapping function, i.e.(2.4)The problem of the hypotrochoidal hole (i.e. the two-term mapping function) under uniform normal traction has been solved previously (Zimmerman 1986). Some special cases, such as holes represented by the first three or four terms in the Schwarz–Christoffel mapping for an (*n*+1)-sided equilateral polygon, have been discussed by Jasiuk *et al*. (1994) and Kachanov *et al*. (1994). In the present work, we will consider the somewhat more general case of holes having (*n*+1)-fold symmetry that can be represented by the first three or four terms in this mapping function. The solution will be derived in detail for the three-term case, whereas for space considerations the solution will be presented without detailed derivation for the four-term case, which poses no fundamental additional difficulties aside from increased algebraic complexity. To simplify the notation slightly, we write our mapping function as(2.5)

In order for the mapping to be conformal, and for the contour not to have any self-intersections, it must never be the case that *ω*′(*ζ*)=0 along the contour. This poses restrictions on the allowable range of values for the *m*_{i} coefficients. If *ω*′(*ζ*)=0 for some value of *ζ* on the unit circle, as the *m*_{i} values increase, this will first occur at (*n*+1) equally spaced points that include the point corresponding to *ζ*=1. Hence, the restrictions for the *m*_{i} can be found by setting *ω*′(1)=0. For the two-term mapping, this leads to the restriction *m*_{1}≤1/*n*. For the three-term mapping, the condition is , and for the four-term mapping, the condition is .

## 3. Determination of the complex potentials

If all the terms in the boundary condition (2.2) are expressed in terms of *ζ*, this equation becomes(3.1)where we use *σ* to denote values of *ζ* along the unit circle *γ* in the *ζ*-plane. To avoid a cumbersome notation, hereafter we write *ϕ*[*z*(*σ*)]=*ϕ*(*σ*), and *ψ*[*z*(*σ*)]=*ψ*(*σ*), in which case (3.1) becomes(3.2)

If the traction acting on the boundary is a hydrostatic pressure of magnitude *p*, it can easily be shown (Sokolnikoff 1956) that *F*=−*pz*=−*pω*(*σ*). Taking *p* to be of unit magnitude, equation (3.2) takes the form(3.3)

For a problem with no body forces and no stresses acting at infinity, both potentials will be represented (Sokolnikoff 1956, p. 279) by power series that converge for *ζ*≤1, i.e.(3.4)(3.5)We now expand the second term in (3.3) in a power series in *σ*. Considering first the case of a three-term mapping function, we will need the following results:(3.6)(3.7)(3.8)in which case(3.9)

After inserting (3.9) into (3.3), the coefficients of each power on *σ* can be equated on both sides of the equation. This procedure yields(3.10)In principal, this procedure could yield all the coefficients *b*_{k}, even those for *k*>*n*, but it would require generating all terms in the series expansion (3.9), which is not practical. So, to find the additional non-zero terms in *ϕ*, if any, we divide both sides of equation (3.3) by 2*π*i(*σ*−*ζ*), and then integrate around the unit circle *γ*:(3.11)Evaluation of the integrals, ignoring any resulting constant terms (which correspond to rigid-body displacements), yields(3.12)

If four terms are taken in the mapping function, the same procedure as used previously eventually leads to(3.13)where(3.14)

(3.15)

To calculate the second potential function, *ψ*, we start by taking the conjugate of equation (3.3), making use of the fact that *ϕ* is now known, and expand out the product, to obtain(3.16)We now divide both sides of this equation by 2*π*i(*σ*−*ζ*), and integrate around the unit circle in the *ζ*-plane. As the mapping is single-valued, *ω*′(*σ*) can never vanish on or within *γ*. Most of the integrals therefore contain only simple poles. After ignoring any resulting constants, the evaluation of these integrals yields(3.17)The remaining integrand has poles at *σ*=*ζ* and 0. To find the residue at *σ*=0, we expand the denominator as a power series in the numerator:(3.18)The residue at *σ*=0 is the coefficient of 1/*σ* in the above series, namely, 1/*ζ*^{n−1}. Including the terms that multiply this integral in (3.18), we find a contribution of . Similarly, the residue associated with the simple pole at *ζ*=*σ* leads to a contribution . Therefore, the second potential is(3.19)

The full expression for the second potential associated with the *four*-term mapping function is similarly found to be(3.20)where the *B*_{i} terms are defined in (3.14) and (3.15). Although these expressions seem to have singularities at *ζ*=0, algebraic manipulation shows that these apparent singularities are removable, and *ψ*(ζ) is indeed analytic inside the unit circle, as originally claimed.

## 4. Determination of the pore compressibility

Two pore compressibility parameters can be defined for a porous material: *C*_{pc}, the fractional decrease in the area of the pore due to an external hydrostatic confining pressure of unit magnitude, and *C*_{pp}, the fractional increase in the area of the pore due to an internal hydrostatic pore pressure of unit magnitude (Zimmerman 1991; Detournay & Cheng 1992). Simple superposition arguments show that these two parameters are related by , where (*κ*−1)/2*G* is the two-dimensional areal bulk compressibility (Zimmerman 1986). For plane strain, this coefficient is equal to (1−2*ν*)/*G*.

By definition, . In the context of linear elasticity, and bearing in mind our choice of *p*=1, we see that *C*_{pp}=Δ*A*/*A*, where *A* is the initial area of the pore. The initial area *A* can be expressed in terms of the mapping function, as follows (Zimmerman 1986), using Green's theorem:(4.1)where the minus sign is introduced to ensure that the integration is performed in the anti-clockwise sense with respect to the hole in the *z*-plane. For the three-term mapping function, the functions needed for this last integral are(4.2)(4.3)Therefore, the initial area is(4.4)Similarly, it can be shown that the original area of the pore represented by the four-term mapping function is given by(4.5)

The area change can be determined by integrating the normal component of the displacement around the pore boundary. The outward unit normal vector to the contour is given by , the arc-length of the contour is d*s*=|*z*′(*α*)|d*α*, the displacement vector is ** u**=(

*u*,

*v*), and the normal component of the displacement is

*u*

_{n}=

*.*

**u***, so(4.6)In vector notation, the displacement is represented by*

**n***=(*

**u***u*,

*v*), whereas in the complex plane this same displacement is written as

*U*=

*u*+i

*v*; both notations have been used in (4.6), where appropriate.

Inserting the complex displacement (2.1) into (4.6) shows that the area change can be broken up into three parts, which we evaluate separately. In each case the integral can be evaluated by elementary means:(4.7)(4.8)and similarly, but omitting the details,(4.9)Hence, the total area change of the ‘three-term’ hole is given by(4.10)The corresponding terms in the area change for the four-term mapping function are(4.11)(4.12)(4.13)where *B*_{1} and *B*_{2} are as defined in (3.14) and (3.15).

Finally, the pore compressibility with respect to pore pressure is found from *C*_{pp}=Δ*A*/*A*, where *A* is given by (4.4) or (4.5), and Δ*A* is given by (4.10) or (4.11)–(4.13). The pore compressibility with respect to the far-field stress, *C*_{pc}, is related to *C*_{pp} by *C*_{pc}=*C*_{pp}+*C*_{o}, which is to say, . The pore compressibility *C*_{pc} is related to the effective macroscopic compressibility by *C*_{eff}=*C*_{o}+*ϕC*_{pc}. In order to relate our results to previous work on the effective moduli problem, the following discussion will be phrased in terms of *C*_{pc}.

## 5. Results and discussion

To validate our results, we compare them against several special cases that have been obtained previously, and also compare them to the values obtained by boundary-element calculations. The boundary-element calculations were performed using a code written by Martel (Martel & Muller 2000), which is a simplified version of the more general two-dimensional BEM code from Crouch & Starfield (1983) that is based on the displacement discontinuity method. Martel's code is in a sense optimized for the problem of a single void or crack in an infinite elastic body, with many of the options included in the original code removed, thus rendering it easier to use for our problem.

In our calculations, all far-field stresses are set to zero, as are the body forces. A uniform normal traction of unit magnitude is prescribed over the surface of the hole. The cavity boundary is discretized into a number of equal-length elements. The number of elements is always taken to be a multiple of *n*+1, the number of ‘sides’ of the pore, to ensure that two boundary elements meet precisely at each corner or cusp, so that the corners are not chopped off. We generally found that roughly 300 boundary elements are sufficient to achieve convergence of the computed compressibilities.

As mentioned previously, the results for pores represented by a two-term mapping function (hypotrochoid) have been obtained by Zimmerman (1986), Jasiuk *et al*. (1994) and Kachanov *et al*. (1994): . The present results indeed reduce to these values identically when *m*_{2}=*m*_{3}=0. Moreover, it can be shown by tedious algebraic manipulation of our solution that *C*_{pc} *always* depends on the elastic moduli of the host material only through the multiplicative factor (1−*ν*)/*G*. Hence, the purely geometric effect of pore shape can be discussed in terms of the dimensionless pore compressibility, *GC*_{pc}/(1−*ν*). For a circular hole, this dimensionless compressibility is exactly 2.

Next, consider the quasi-polygons represented by the first three or four terms of the Schwarz–Christoffel mapping:(5.1)Table 1 shows our results for the three-term and four-term holes, compared with the values obtained by Jasiuk *et al*. (1994) and Kachanov *et al*. (1994). With the exception of the result obtained by Kachanov *et al*. (1994) for the three-term quasi-triangle, all the values agree to three decimal places. A more sensitive test of the results is shown in table 2, where the percentage changes in *C*_{pc} obtained by going from two to three terms, and then from three to four terms, are shown. The agreement between the present analytical results, the boundary-element results and the results reported by Jasiuk *et al*. (1994) are excellent.

It would be of interest to have ‘exact’ values for the compressibilities of the equilateral polygons that are represented by an infinite number of terms of the Schwarz–Christoffel mappings. If the calculated compressibilities are plotted against 1/*N*^{2}, where *N* is the number of terms taken in the mapping function, the results seem to be converging linearly with 1/*N*^{2} (figure 2). We therefore calculate the values for *N*→∞ by extrapolating from the values for *N*=3 and 4, using 1/*N*^{2} as the independent variable, and letting 1/*N*^{2}→0. These extrapolated values are given in table 3, along with the values that are computed using the boundary-element method. The extrapolated analytical values generally agree with the boundary-element calculations to three digits, thus giving confidence that these values are accurate to at least three digits. It should also be noted that the extrapolated values are in each case much closer to the BEM values than are the values calculated using only four terms in the mapping function, which provides some justification for the extrapolation procedure.

These results show that the normalized compressibility of an equilateral triangle is 3.25, or about 63% larger than that of a circle. The compressibility of these equilateral polygons decreases rapidly as the number of sides increases: a square is 21% more compressible than a circle, a pentagon 10% more compressible and a hexagon only 5% more compressible. Regular polygons with more than six sides will therefore have compressibilities that are essentially the same as that of a circle.

Although both the conformal mapping approach and the BEM method can in principle be used to compute the compressibility of a pore of any shape, calculation of the mapping coefficients for complex pore shapes is extremely tedious, and computationally non-trivial (Sisavath *et al*. 2001; Tsukrov & Novak 2002). So, it would also be useful if the pore compressibility could be calculated from some simple geometric attributes of the pore shape, without requiring elaborate analytical or numerical calculations. Such a capability would be useful in attempts to estimate elastic moduli from images of heterogeneous media (e.g. Tsukrov *et al*. 2005).

Zimmerman (1986) suggested that the pore compressibility *C*_{pc} scales (approximately) with *P*^{2}/*A*, where *P* is the perimeter and *A* is the area. Forcing this scaling law to be exact for a circular hole leads to the approximation(5.2)Zimmerman (1991) showed that this approximation has an error of less than 8% for all hypotrochoids (i.e. two terms in the mapping function), and an error of about 23% for thin, crack-like pores, which he suggested might be the ‘worst-case’ shape. Tsukrov & Novak (2002) verified that (5.2) had an error of only 8% for a single, arbitrarily drawn irregular pore. In order to assess the usefulness of this approximation for other pore shapes, we have tested it against some of the shapes that can be generated using our two-term, three-term and four-term mappings.

Tables 4 and 5 show the pore compressibility values for some two-term and three-term pores having threefold symmetry (i.e. *n*=2). The compressibility increases as the hole becomes ‘less circular’, in some macroscopic sense. The presence of pointed cusps also causes the compressibility to increase. The scaling law (5.2) accounts in a rough sense for the effect of pore shape on compressibility, but in general may have an error of as much as 20%.

## 6. Conclusions

We have derived closed-form expressions for the compressibilities of isolated pores in an infinite, isotropic two-dimensional elastic body that have *n*-fold axes of symmetry, and which can be represented by at most four terms in the mapping function that conformally maps the exterior of the pore into the interior of the unit circle. Our results were validated against some special cases that have previously been derived by Zimmerman (1986), Jasiuk *et al*. (1994) and Kachanov *et al*. (1994), and also against boundary-element calculations. By extrapolation of the results for pores obtained from three and four terms of the Schwarz–Christoffel mapping function for regular polygons, we explicitly found (to at least three digits) the compressibilities of a triangle, square, pentagon and hexagon. Specific results for some other pore shapes, more general than the quasi-polygons obtained from the Schwarz–Christoffel mapping, were also presented.

For the family of pores analysed in this work, the compressibility of the pore with respect to a far-field hydrostatic stress, *C*_{pc}, was explicitly found to depend on the elastic moduli of the host material only through the multiplicative factor (1−*ν*)/*G*. In fact, it can be shown, using a result of Jasiuk (1995), that the compressibility of a pore of *any* shape is proportional to (1−*ν*)/*G*; however, the proof is not straightforward, and will not be given here. It follows that the ratio of the compressibility of a given pore to that of a circular pore, which is 2(1−*ν*)/*G*, depends only on the pore shape, and not on the elastic moduli.

An approximate scaling law proposed by Zimmerman (1986) for the compressibility, in terms of the ratio of perimeter-squared to area, was also tested. This expression was generally found to overestimate the compressibility, by an amount that varied from about 0 to 21%. Although the error could perhaps be reduced, on average, by using a different numerical constant in the correlation, it is clearly not entirely true that the dimensionless compressibility depends only on the dimensionless parameter *P*^{2}/*A*.

The analysis presented previously, insofar as it relates to the effective moduli problem, is applied only to hydrostatic far-field stresses. In the hydrostatic case, the orientation of the pores with respect to any arbitrary coordinate system is immaterial. Hence, for each specific pore shape, our results can be used to yield the effective bulk modulus of a material containing a collection of such pores whose orientations are randomly distributed, through the relation . In order to find the effective *shear* modulus of such a material, we would need to solve the problem of an infinite elastic body, containing a single pore, with a stress state of pure shear at infinity. If we were only interested in the effective shear modulus, rather than the details of the stress distribution around the pore, we would only need to calculate the excess strain energy due to the presence of the pore. As shown by Kachanov *et al*. (1994), for all symmetrical pores except those having fourfold symmetry, this energy term is independent of the orientation of the pore with respect to the direction of shear. Hence, for all symmetric pores except those having fourfold symmetry, we would only need to consider a single specific orientation of the pore with respect to the far-field shear stress. The solution to this problem could be obtained by following roughly the same steps as were used previously for the case of hydrostatic loading, although the algebraic details would probably be somewhat more complicated. The main difference in the problem formulation would be a different expression for the function *F* appearing on the right-hand side of (3.2); the procedure for calculating this function *F* is given by Muskhelishvili (1963) and England (1971).

## Footnotes

- Received October 9, 2005.
- Accepted January 6, 2006.

- © 2006 The Royal Society