## Abstract

In Parnell & Abrahams (2008 *Proc. R. Soc. A* **464**, 1461–1482. (doi:10.1098/rspa.2007.0254)), a homogenization scheme was developed that gave rise to explicit forms for the effective antiplane shear moduli of a periodic unidirectional fibre-reinforced medium where fibres have non-circular cross section. The explicit expressions are rational functions in the volume fraction. In that scheme, a (non-dilute) approximation was invoked to determine leading-order expressions. Agreement with existing methods was shown to be good except at very high volume fractions. Here, the theory is extended in order to determine higher-order terms in the expansion. Explicit expressions for effective properties can be derived for fibres with non-circular cross section, without recourse to numerical methods. Terms appearing in the expressions are identified as being associated with the lattice geometry of the periodic fibre distribution, fibre cross-sectional shape and host/fibre material properties. Results are derived in the context of antiplane elasticity but the analogy with the potential problem illustrates the broad applicability of the method to, e.g. thermal, electrostatic and magnetostatic problems. The efficacy of the scheme is illustrated by comparison with the well-established method of asymptotic homogenization where for fibres of general cross section, the associated cell problem must be solved by some computational scheme.

## 1. Introduction

A classical problem in the mechanics of inhomogeneous media is to attempt to replace the two-dimensional potential problem ∇⋅(*μ*(** x**)∇

*w*(

**))=0, where**

*x***=(**

*x**x*

_{1},

*x*

_{2}) and

*μ*(

**) and**

*x**w*(

**) are periodic (scalar) functions that vary rapidly with**

*x***, by an equivalent problem of the form**

*x**w*

_{*}is the leading-order displacement field and

*** is the second-order tensor of effective shear moduli [1]. To do this, a so-called**

*μ**separation of scales*between the micro and macro lengthscales must be assumed. In Cartesian coordinates, (1.1) can be written in the indicial form

*ij*th component of the tensor in Cartesian coordinates. For orthotropic materials, for example,

*homogenization*and

*** depends strongly on the geometrical and physical properties of the medium in question [2,3]. Noting that the equations arise from the equilibrium equation ∇⋅**

*μ***=0, where**

*σ***=**

*σ***⋅**

*μ***, and**

*e***=∇**

*e**w*, it is stressed that the problem posed has broad applicability, as summarized in table 1. To fix ideas here, the application to antiplane elasticity shall be described so that

*w*is the out-of-plane displacement.

Assume now that *μ*(** x**) is piecewise constant, and a cross section of the medium takes the form as depicted in figure 1 so that the material can be classified as a unidirectional fibre-reinforced composite (FRC), so that fibres of general cross sections shall be considered. Such media are used in a multitude of applications where rather specific material properties are required in order to perform a task effectively and where naturally available homogeneous media are not effective or efficient [2]. In particular, materials of this form can provide high tensile stiffness and/or high directional conductivity while remaining relatively light by using only a small volume fraction of the fibre phase.

The microstructural lengthscales of such inhomogeneous media are becoming ever smaller with an increasing ability to engineer microstructures for improved macroscopic performance [4]. It is frequently convenient to consider the problem in the context of wave propagation so that an inertia term is added to the governing equation, and the effective medium is then governed by, for example
*ρ*_{*} is the effective mass density and where time-harmonic motion *e*^{−iωt} has been assumed, identifying *ω* as the angular frequency. By considering this dynamic context, a natural macroscopic lengthscale is introduced: the wavelength of the propagating effective wave. Intrinsic to the subject of homogenization then, where quasi-static properties are determined, is that the wavelength of the propagating wave is much longer than the microstructural lengthscale.

A number of methods have proved extremely successful at predicting ** μ*** in both the static and quasi-static regimes, including the method of asymptotic homogenization (MAH) [1,3,5–7], the equivalent inclusion method [8,9], boundary element solutions of integral equations [10] and the use of Fourier transforms [11–13]. All schemes rely on separation of scales and periodicity—the fact that a periodic cell is representative of the entire medium. We stress, however, that general fibre cross sections shall be considered here, and the semi-analytical approach developed is able to determine explicit expressions for the effective properties with minimal computational resource. This is in contrast with existing homogenization schemes where numerical methods are required and results must be computed each time a parameter is varied, e.g. volume fraction. Having general expressions in the volume fraction is of great utility and practicality in materials design and optimization.

Moving away from the separation of scales regime, a significant amount of work has been undertaken that incorporates propagation at finite frequencies and particularly when the wavelength of the propagating wave is of the order of the microstructure, when dynamic effects become important. In this context, the microstructure can be designed to manipulate the wave-carrying capabilities of the medium. In particular for periodic materials, the microstructure can be chosen in order that the medium acts as a wave filter, for waves in certain frequency ranges (the so-called stop bands) are unable to propagate. Methods devised to determine the band gap structure are the plane wave expansion technique [14], multipole methods [15], and high-frequency homogenization [16]. See the useful review [17] for further details. Low-frequency effective properties can thus be deduced numerically from these schemes by considering propagation near the origin of the dispersion curves in question. If from the outset, however, one is interested purely in the low-frequency limit where homogenization applies, then there is no need to determine this full dynamic behaviour.

Here, then, a strict separation of scales shall be considered; the homogenization regime is assumed to hold and the material responds as an effective medium with uniform properties. The key novelty of the proposed scheme is the form of solutions that are derived as shall be illustrated below. The method to be discussed extends the work in [18] (referred to as PA below), where a new homogenization scheme was devised based on the integral equation form of the governing equation, considering antiplane wave propagation in the low-frequency limit, so that an equation of the form (1.2) was derived, and the leading-order result was determined. The work of PA was itself inspired by the method introduced in [19], in which expressions for the effective elastic properties of three-dimensional random particulate media were determined, but restricted to the dilute-dispersion limit [20]. Returning to the two-dimensional periodic medium considered here, although the method itself may be applied to a material of any macroscopic elastic symmetry, attention shall be restricted to microstructure that gives rise to orthotropic effective properties for clarity of exposition. In this case, *circular cylindrical fibres* the following result is derived in this work:

As in PA, it shall be assumed that all fibres in the composite are identical, that all phases are isotropic and that the lattice geometry and fibre cross sections are restricted such that the effective medium appears to be, at most, orthotropic on the macroscale. In §2, the governing equations are summarized and the associated integral equations are derived. The integral equation methodology is then described in §3 and the manner in which effective properties are determined is presented in §4. Results are given in §5 for a variety of media before concluding in §6. Some of the more detailed analysis is presented in appendices in order to aid the flow of the paper.

It is reiterated that although the method is described here in the context of antiplane elasticity and expressions for the effective antiplane shear moduli

## 2. Governing equations

The microstructure of the medium in question is illustrated in figure 2. Unidirectional fibres, considered so long as for the problem to be assumed two dimensional, are positioned on a periodic lattice and their cross section is taken as general with the restriction that the macroscopic anisotropy is at most orthotropic. The problem shall be formulated in Cartesian coordinates, with the *x*_{3}-coordinate running parallel to the fibre axis and the *x*_{1}*x*_{2}-plane being the plane of periodicity. The location of the centre of the (*s*, *t*)th periodic cell is defined by the lattice vector

for some two-dimensional vectors *l*_{1},*l*_{2}. Attention is restricted to the case where each periodic cell contains a single fibre. The lengthscale *q* can be considered as the characteristic lengthscale of the periodic cell and therefore of the microstructure of the medium. Consider a distribution of fibres that induces macroscopic *orthotropy* so that
*V*_{st} embedded within the host phase which is denoted by *V*_{0}. The shear modulus and mass density of the host (fibre) is denoted as *μ*_{0} (*μ*_{I}) and *ρ*_{0} (*ρ*_{I}), respectively.

The lengthscale *a* associated with the fibre cross section is introduced by defining the boundary of the cross section in terms of circular cylindrical coordinates, i.e.
*f*(*θ*)≥1 and for some *θ*∈[0,2*π*], *f*(*θ*)=1. Hence, *shape* of the fibre cross section into the analysis.

Horizontally polarized shear (SH) waves, otherwise known as antiplane waves, are considered to propagate through the medium. The associated time-harmonic elastic displacement, polarized in the *x*_{3}-direction is denoted by *w*(** x**), where it is recalled that

**=(**

*x**x*

_{1},

*x*

_{2}). The elastic displacement

*w*is governed by

*x*

_{1},∂/∂

*x*

_{2}) and where

*k*

_{0}=

*ω*/

*c*

_{0}and

*k*

_{I}=

*ω*/

*c*

_{I}are the wavenumbers of the host and fibre, respectively, and where

*indicator function*

*χ*(

**) is defined by**

*x*Combining (2.3) and (2.4) appropriately, imposing boundary conditions of continuity of displacement and traction on fibre/host interfaces, the problem is restated in integral equation form as
_{y}=(∂/∂*y*_{1},∂/∂*y*_{2}). Non-dimensionalizing using the scalings ** p**=

*s*

*l*_{1}+

*t*

*l*_{2}and the cross-sectional area of the periodic cell is

*x*

_{3}-direction,

*D*| is the (non-dimensional) cross-sectional area of the fibre.

Upon dropping the hat notation, the non-dimensional integral equation takes the form
*d*=*ρ*_{I}/*ρ*_{0} and *m*=*μ*_{I}/*μ*_{0} are the contrasts in mass density and shear moduli, *ε*=*qk*_{0} and *G*_{ε}(** x**−

**)=(1/4**

*y**i*)

*H*

_{0}(

*ε*|

**−**

*x***|). Note that as**

*y**γ*

_{e}≃0.577216… is Euler’s constant. Having already dropped hats, referring to (2.2), in non-dimensional coordinates the boundary of the fibre is, therefore, described by

*r*=ℓ

*f*(

*θ*),

*θ*∈[0,2

*π*], where ℓ=

*a*/

*q*. The non-dimensional fibre cross section is easily shown to be

*volume fraction*of the fibre cross section within the periodic cell is then

*ϕ*below, in some instances it is more convenient to work with ℓ initially. However, in the end using (2.8), general expressions are determined in terms of

*ϕ*. Attention here is also restricted to the scenario where

*τ*remains

*O*(1), with respect to

*ϕ*. Therefore, from (2.8),

*ak*

_{0}≪1 and

*ε*=

*qk*

_{0}≪1 and it can be assumed here that

*ak*

_{0}=

*O*(

*ε*). In [21,22], the regime where

*ak*

_{0}≪1 but

*ε*is not restricted to being small was considered. Note that this requires

*ϕ*≪1. If one wished, the method introduced in the next section could be modified in order to consider this regime.

It is important to note that as they stand the infinite sums in (2.5) are not absolutely convergent. This is an issue that arises frequently in effective media problems where the material in question is of infinite extent. Sensible results can be derived by a number of different approaches. In PA, this issue was dealt with by allowing a small amount of ‘loss’ in the system. This means taking wavenumbers of the form *k*+*iη* where *η*≪1. Non-zero *η* ensures convergence (in this case of integrals that arise when carrying out the sum to integral step, as developed in PA) and the limit

## 3. The integral equation method

In PA, it was shown that setting *m*=1 leads to the result *ρ*_{*}=(1−*ϕ*)+*dϕ* for the non-dimensional effective density (scaled on *ρ*_{0}) in the quasi-static limit. This is an exact result in the separation of scales regime. Here, without loss of generality and in order to determine the effective shear modulus, set *d*=1 in (2.5). Differentiate both sides of the resulting integral equation with respect to *x*_{k},*k*=1,2 to yield
** x**∈

*V*

_{ab}, i.e. within the (

*a*,

*b*)th fibre, which has position vector

**=**

*r**a*

*l*_{1}+

*b*

*l*_{2}. By taking the Taylor expansion of Green’s function about the point

**=**

*y***=(**

*p**p*

_{1},

*p*

_{2})=

*s*

*l*_{1}+

*t*

*l*_{2}(i.e. about the centre of the (

*s*,

*t*)th fibre), (3.1) becomes

*i*th derivative with respect to

*y*

_{k}and

*displacement-gradient moments*of order

*i*+

*j*, recalling that

*w*is a piecewise smooth function,

*w*and all its derivatives will be

*O*(1).

Note that the term for (*s*,*t*)=(*a*,*b*) is not included in the summation, nor has the Taylor series of Green’s function been taken in this term. This is because Green’s function is singular in the domain *V*_{ab} since ** x** is contained in this region. The assumption that one can Taylor expand Green’s function puts restrictions upon the parameters

*ε*and

*ϕ*. Either (i)

*ε*≪1 in which case

*ϕ*is unrestricted

*or*(ii)

*ε*=

*O*(1) and then

*ϕ*≪1 is required. Here, only scenario (i) is considered; see also the discussion at the end of the last section.

Proceed now by defining the operation *x*_{1}−*r*_{1})^{δ}(*x*_{2}−*r*_{2})^{ξ}/(*δ*!*ξ*!) and integrating in the ** x** plane over the domain

*V*

_{ab}, where

**=(**

*r**r*

_{1},

*r*

_{2}). Hence, consider

**=**

*x***to obtain, after some rearrangement**

*r**G*

_{ε}/∂

*x*

_{k}=−∂

*G*

_{ε}/∂

*y*

_{k}has been employed. Here, the outermost summation notation refers to summing over every fibre location

**=(**

*p**p*

_{1},

*p*

_{2}) except for

**=**

*p***. Furthermore, the following terms have been defined:**

*r**solely*in terms of the constants

### (a) The shape factor

The term ** y**=

**due to the presence of derivatives of Green’s function in the integrand. However, as was shown in PA, this apparent singular contribution is found to be zero by splitting the domain**

*x**V*

_{ab}up into a non-singular part

*V*

_{ab}∖

*C*

_{ψ}and apparently singular part

*C*

_{ψ}, where

*C*

_{ψ}is a disc of radius

*ψ*≪1 with origin

**=**

*y***. Therefore, all that remains to consider are integrals of the type**

*x**G*

_{ε}/∂

*x*

_{k}=−∂

*G*

_{ε}/∂

*y*

_{k}, the

*x*

_{k}derivative may be taken inside the

**integral, and as the range of integration does not include the region where**

*y**G*(

**−**

*y***) is singular, exchanging the order of integration is permissable. Retaining the leading-order Green’s function as**

*x**G*

_{0}was defined in (2.6).

It is convenient here to represent (3.10) as a series expansion, i.e.
*f* involved. For elliptical fibres

Substituting (3.11) into (3.9) and adjusting the indices one obtains
*f*(*θ*)=1 and it follows that *δ*+*ξ*+*α*+*β* is odd. As

Equations (3.12) and (3.13) coupled with (3.4) give rise to an infinite homogeneous system of linear equations for the displacement gradient moments

### (b) Wave-like solutions

Noting (2.9), plane wave type solutions of (3.4) are sought of the form
*q*) in the direction of the angle subtended from the *x*_{1}-axis, *Θ*. Furthermore, define

Therefore, using (2.8), (3.7), (3.12)–(3.15) in (3.4) (recalling *k*=1,2) and recalling (2.9), at leading order (with respect to *ε*)

### (c) Lattice sums

First pick out the singular, non-integrable behaviour of the derivatives of Green’s function in the lattice sum (3.18), writing
*S*[*m*,*n*] as the singular part of the derivative. Then for a given *m*,*n*, (3.18) can be written as
*m*=*i*+*α*+*p* and *n*=*j*+*β*+*q*, such that (*p*,*q*)∈{(2,0),(1,1),(0,2)}, one is left with an integral that is *O*(1) with respect to *ε*, multiplied by a factor *ε*^{i+j+α+β}. Therefore, the only term from this step that remains in the limit *i*=*j*=*α*=*β*=0. Thus,
*p*,*q*)∈{(2,0),(1,1),(0,2)}.

Recalling the definition of ** Γ** in and below (3.14), taking

*Θ*=0 and defining

*γ*(0)=

*γ*

_{1}, so that the wavevector is

**=**

*Γ**ε*(

*γ*

_{1},0) with

*γ*

_{1}being the effective wavenumber in the

*x*

_{1}-direction, it was shown in the appendices to PA that (using the present notation)

If one wishes to determine the effective wavenumber in the *x*_{2}-direction (seeking *π*/2. Performing this action leaves the above integrals unchanged, the wave is considered to propagate in the (new) *x*_{1}-direction, merely using notation *γ*_{2} instead of *γ*_{1} to indicate the wavenumber associated with the new material direction of propagation.

One can determine a number of results for specific *L*_{0}[*m*,*n*] straightforwardly. Firstly, since the singular part of *G*_{0} satisfies Laplace’s equation (except at the singular point, which is not important in the lattice summations as they exclude this point), the lattice sums will satisfy
*both* *m* *and* *n* are even will give a non-zero *L*_{0}[*m*,*n*]. For example, consider the case of *L*_{0}[1,2]=−*L*_{0}[3,0] due to (3.26). The associated *S*_{0}[*m*,*n*] takes the form
*u*_{1} in the numerator ensure that the summation is zero. This same reasoning gives *L*_{0}[*m*,*n*]=0 whenever *m*+*n* is odd. It is shown in appendix B how *L*_{0}[*m*,*n*] can be straightforwardly determined when it is non-zero, i.e. when *m*+*n* is even.

### (d) The asymptotic system in *ϕ*

Using (3.23) in (3.16) and (3.17) the leading-order system, with respect to *ε*, is
*γ*_{1} and associated eigenvectors comprising the moments *γ*_{1} appears only in the term *g*(*γ*_{1}). To make progress take expansions in the volume fraction, posing
*g* is motivated by (3.25) and the fact that

Note the scaling in *ϕ* (or ℓ) of the wave-type solutions (3.14), so that moments whose indices total an odd number will have a fractional power at leading order in *ϕ*. For instance, *ϕ*^{1/2} should be included in (3.29)–(3.31). However, restricting attention to shapes featuring a rotational symmetry of *π*, henceforth referred to as *centrally symmetric* shapes, and by observing the properties of the integrals in appendix A it can be shown that such fractional powers are not required in the expansions. Proof of this is given in appendix C, and the remainder of the methodology shall be discussed in the context of the restriction of fibre cross sections to central symmetry.

Given the scalings involved, using (3.25) and *O*(*ϕ*) term. The concern of the next section is the determination of the coefficients *h*_{j} from the linear system.

## 4. Determining explicit forms for the effective properties

Equation (3.32) provides an explicit form for the effective shear modulus, as a rational function in *ϕ*, providing the coefficients *h*_{j} can be determined. Note that the approximation provided in PA was exactly this with both numerator and denominator truncated at *O*(*ϕ*). Here, the objective is to determine higher-order coefficients for a variety of fibre cross sections.

There is a straightforward algorithmic mechanism for determining the coefficients *h*_{j}. Hereon in the term ‘(*δξ*) equations’ shall refer to (3.27) and (3.28) for fixed *δ* and *ξ*. Start by noting that expressions for the *h*_{j} coefficients arise by considering the first (00) equation, (3.27), using expansions (3.29)–(3.31) and equating each order in *ϕ* of the resulting equation. To determine all coefficients up to *h*_{N}, say, the first (00) equation must be considered up to order *ϕ*^{N+1}. Each order provides an expression for a coefficient *h*_{j} in terms of the eigenvector components *δξ*) needs to be made and hence which additional (*δξ*) equations need to be considered in order to solve the linear system for the required terms *h*_{j}. The following examples illustrate the implementation of this scheme; orderings now refer to *ϕ*.

### (a) Elliptical cylindrical fibres

Consider elliptical cylindrical fibres with semi-axes *a* and *b* in the *x*_{1} and *x*_{2} directions, respectively. At *O*(1), noting that *τ* is defined for a given shape in (2.8) and writing without loss of generality that *O*(*ϕ*) terms give (upon using the properties of the lattice sum)
*O*(*ϕ*) gives *L*_{0}[2,0]=−*L*_{0}[0,2]=−(*S*_{j}−1)/2, *C*_{4j} and *C*_{6j} depend upon the fourth- and sixth-order lattice sums, respectively, and the *p*_{ij} coefficients are rational functions in *m* and also depend on aspect ratio. Their forms are too lengthy to be given here but can straightforwardly be derived by implementing the algorithm above. As one should expect through checking for consistency with other leading-order approximations

### (b) Circular cylindrical fibres

To illustrate specific details of the derivation of higher-order terms, take the circular cross-sectional case, *a*=*b*=1 for which the details of the last section are clearly valid. Consider the first (00) equation to *O*(*ϕ*^{2}). After evaluating the *h*_{−1} and *h*_{0}, the equation reduces to
*ϕ* in order to obtain *h*_{1}. Using the naming convention that the collective set of (*ij*) equations where *i*+*j*=*n* shall be referred to as the order *n* equations, here the interest is therefore in the order 2 equations.

It transpires that these order 2 equations partition into two non-trivial decoupled sub-systems: the first (20) and (02) equations and the second (11) equation form a system in *h*_{1}. Hence, making use of (4.1) for *a*=*b*=1 and all other solutions from prior orders of *ϕ*, the leading-order equations of the former sub-system take the form
*h*_{1}=0.

Proceeding to general orders, with the aid of a symbolic package such as Mathematica, for general parallelogram lattices, results for circular cylindrical fibres take the form
*C*_{6j}=0 which is why order *ϕ*^{6} and *ϕ*^{7} terms appear in the general case and not the square lattice case.

In the special case when circular cylindrical fibres on a square lattice, expressions for the effective shear moduli take the form

### (c) More general fibre cross sections

Although rather lengthy, analytical forms for the effective shear moduli associated with general cross sections *can* be obtained. The procedure above can be followed and certainly in a symbolic package such as Mathematica, forms can be derived. However, such expressions are too lengthy and cumbersome in general to be provided here. Importantly, for a given cross-sectional shape, using the above procedure, a rational function approximation for the effective properties as a function of *ϕ* can be derived and subsequently used to great utility. In the following section, we discuss results obtained for shapes more general than elliptical and validate the scheme with the classical method of asymptotic homogenization (MAH). The argument for using the present scheme over the MAH (or other methods) is that this integral equation methodology can yield explicit forms, particularly when some aspects of the medium are fixed (e.g. square lattice, fixed *m*), of the effective moduli, retaining dependence on *ϕ*. Such forms can then be used with great rapidity in models without the need for recourse to finite-element simulations as soon as one changes, for example, the volume fraction, as is required in the MAH for general shapes, for example.

## 5. Results

The implementation of the above methodology is now described for a range of geometries in the case of a shear contrast of *m*=18.75 (corresponding to the case of graphite fibres in epoxy, for example). To fix ideas and since general cross sections are of principal interest here, attention is restricted to the case of square lattices. Extension to other lattices is straightforward. Results derived using the methodology shall be compared with those obtained using the MAH [1]. It transpires that the effective antiplane shear moduli (when scaled on that of the host medium) as determined by the MAH can be written in the form (see eqns (3.25) and (3.26) of [1])
**H** may be determined as
** n** is the outer unit normal to the fibre boundary ∂

*V*

_{ab},

**is the short lengthscale of the problem, and**

*ξ***=(**

*N**N*

_{1},

*N*

_{2}) is the solution to the associated

*cell problem*. In the case of orthotropy,

*H*

_{12}=

*H*

_{21}=0, and if the fibre cross section has a rotational symmetry of

*π*/2,

*H*

_{11}=

*H*

_{22}and thus

*H*

_{11}and

*H*

_{22}of the

**H**-tensor. To do this for the integral equation method (IEM), the shear moduli are determined using that scheme and then the components of the

**H**-tensor are computed using expressions (5.1). It will be shown that these components provide a very sensitive measure of the accuracy of the IEM and generally speaking, the method provides excellent accuracy even at extremely high volume fractions, for relatively low-order approximations of the rational function form. Note that when results are presented, we will use the term

*IEM order*to refer to the rational function form (3.32) when order

*n**n*polynomials in

*ϕ*are retained in the numerator and denominator. In particular,

*IEM order 1*refers to the scheme employed in PA, i.e. (1.3).

### (a) Circular cylindrical fibres

Consider circular fibres of radius *r*. Figure 3 compares results for *H*_{11}(*r*)=*H*_{22}(*r*) as obtained from the IEM at various orders, to those obtained using the MAH. Results are determined as a function of the scaled radius *r*=ℓ, related to the volume fraction *ϕ* via equation (2.8) for *f*(*θ*)=1 and *b* focuses on values of the radius close to this packing limit. The highest-order IEM employed here (order 12) begins to deviate from the MAH estimate around *r*=0.48, although this is the deviation when studying the components of **H**. Generally, the impact of this deviation on *effective properties* is weaker as will be seen shortly.

### (b) Elliptical cylindrical fibres

Figure 4 illustrates results obtained in the case of elliptical fibres with aspect ratio (major axis divided by minor axis) of *ϵ*=2 and here *r* is the major axis. As figure 4*a* illustrates, agreement is excellent in general and it is only near the packing limit where the IEM and MAH results show any deviation, but only for *H*_{11} since this is the component that is affected by the near-neighbour interaction. Figure 4*b* magnifies this large radius region to clearly illustrate this effect. The order 8 IEM expression begins to show the upturn required in order to match the MAH at high radius for *H*_{11}. The **H**-tensor components are a very sensitive measure of the accuracy of the IEM scheme however. It will be shown in the next section that when the effective properties are computed such high-order accuracy in these components is not required in order to provide good results.

### (c) Results for other fibre cross sections

The scheme is intended to work for fibres of general cross section. Having established promising results for fibres of elliptical cross section, it is now of interest to examine the efficacy of the method for more general cross sections and the polygonal case shall be considered first. In figure 5*a*, the results for *H*_{11}(*r*) and *H*_{22}(*r*) are examined in the case of rectangular-shaped fibres with an aspect ratio of 2 and where *r* is the long axis of the rectangle. Figure 5*b* illustrates results for the specific effective shear modulus **H**-tensor components exhibit, i.e. in general even for quite small-order approximations the IEM appears to provide exceptional predictions of the effective shear moduli, even up to extremely high-volume fractions.

As the rectangular fibres approach the edge of the cell clearly higher-order approximations are required. In this instance, as is shown in figure 5*a* the IEM order 6 expression is able to accurately replicate the upturn as *ϕ* of layered material–here this is clearly *ϕ*=1/2 and these results for *ϕ*=1/2 are indicated by the black squares in figure 6*a*,*b*.

Moving onto the consideration of a more complex fibre shape, figure 7 illustrates plots of the **H**-tensor in the case of a hexagonal cross section, with a magnified region of high radius plotted in figure 7*b*. The effective shear moduli are plotted in figure 8. Once again, there is extremely good agreement between the two methods and even the slight disagreement in calculations of the **H**-tensor makes little difference to the predictions of the effective moduli.

It should be noted that for a general polygonal shape (including the hexagonal case), the notion of ‘radius’, against which some results above are plotted, needs some interpretation. In fact for the case of square cells considered above, referring to (2.7) the scaled area is *r* here is chosen as *r*∈[0,0.5] for all plots.

Of great utility when describing general-shaped fibres is the so-called *superellipse* function defined as [27]
*a* and *b* are 1, *n*_{1}=100, *m*=6 and *n*_{2}=*n*_{3}=62. Alternatively for the rectangle, *a*=2, *b*=1, *m*=4 and *n*_{1}=*n*_{2}=*n*_{3}=200.

## 6. Conclusion

New, explicit expressions for the effective properties of inhomogeneous media have been derived in the case of the scalar wave problem associated with periodic two-dimensional FRC in which the fibres have non-circular cross section. Results can be interpreted in the sense of low-frequency antiplane elastic waves, acoustics, etc. and are also even more broad in applicability due to the link with potential problems (table 1). The fibres in question can be of general cross section and results obtained have been validated successfully with the classical MAH. In the latter method, the cell-problem for general shapes must always be calculated numerically and the cell problem solution must be recomputed whenever parameters are altered. The advantage of the present scheme is that the resulting expressions can be used for general volume fractions without recourse to computational methods each time the volume fraction is modified. The method is designed with general fibre cross-sectional shapes in mind, but in the special case of centrally symmetric shapes, results such as (4.3) and (4.7) may be obtained. Extensions to the case of full elastodynamics and three-dimensional scenarios for both the potential problem and elastodynamics are currently underway.

## Data accessibility

This work is associated with methodology. Therefore, all research data supporting this article are directly available within this publication.

## Competing interests

The authors claim no conflict of interest.

## Funding

D.J. is grateful to Thales UK and the Engineering and Physical Sciences Research Council (EPSRC) for funding via a CASE PhD studentship. W.J.P. acknowledges the EPSRC for funding his research fellowship (EP/L018039/1). I.D.A. undertook part of this work whilst in receipt of a Royal Society Wolfson Research Merit Award, and part was supported by the Isaac Newton Institute under EPSRC Grant Number EP/K032208/1. R.C.A. would like to acknowledge support from the EPSRC (EP/N013719/1).

## Authors' contributions

W.J.P. and I.D.A. proposed the methodology and developed the leading-order method; D.J. extended the work to the higher-order scenario with assistance from W.J.P. and I.D.A.; D.J. and W.J.P. drafted the manuscript which was subsequently edited by all authors. R.C.A. computed the results associated with the MAH and produced figures 3–8 associated with results, with assistance from D.J. All authors gave final approval for publication.

## Appendix A. Two-dimensional expansion for *J*_{δξ}

Employ the following plane polar coordinates *J*_{δξ}(**y**) and it must be evaluated numerically for a given **y**, hence the motivation for a Taylor expansion to approximate such functions. Therefore, pose expansion (3.11). In terms of the local coordinates, this may be rewritten as
*m*≥0 and *m*≥1 for integer *m* and then integrate over *θ*∈[0,2*π*] where
*J*_{δξ} in (6.1) and (6.2), equate and set *ρ*=1, to obtain
*m*≥1 in the latter. One then finds that
*ζ* integration and a different notation has been associated with this integral form so that it may be referred to later. Also,
*ρ*=1) to the coefficients *J*_{δξ} in (3.10) and (3.11), equating and exploiting the fact that ∇^{2}*G*_{0}=*δ*(** y**−

**), to find**

*x**y*

_{1}−

*r*

_{1})

^{i}(

*y*

_{2}−

*r*

_{2})

^{j}provides the further equations relating the coefficients

*p*,

*q*≥2

**(a) Linear system**

The linear system consists of (A 8) and equations, such as e.g.
**R** is essentially a sparse vector including some constants that arise due to the Laplacian simplification on the left-hand side. The following examples are given due to their relevance to the example systems in §4.

**(i) Elliptical fibre second-order case: **

For this case, only *J*_{00}(*ρ*,*θ*) is sought as a polynomial expansion of the form
*D* constants never exceed 2. In this case, the only non-trivial coefficients that have an influence on the shape factor term are

**(ii) Circular fibre fourth-order case example: **

Note that in order to fully consider equations (3.27) and (3.28) to fourth order the cases *δ*=*ξ*=1 and *δ*=0, *ξ*=2 must also be considered. The case *δ*=2, *ξ*=0 has been chosen merely as an example to illustrate how the system for the shape factor works at this order.

In this instance, the shape factor may be written as

The fact that this shape factor depends on displacement gradient moments for *δ*=*ξ*=1 and *δ*=0, *ξ*=2 illustrates why those choices for *δ* and *ξ* must also be considered in order to fully establish this fourth-order problem.

## Appendix B. Lattice sums

Recall the definition of *L*_{0}[*m*,*n*] from (3.22). While it has already been shown that this sums to zero when either *m*, *n* or both *m* and *n* are odd, and that many sums are equivalent to each other through the identity (3.26), it still remains to evaluate the non-trivial sums, when both *m* and *n* are even. Attention is restricted to parallelogram-shaped lattices, which will in turn restrict the forms of ** u** that appear in the summation. For ease of exposition, consider the further restriction of rectangular-shaped lattices, so that summations are over integer

*u*,

*v*such that

**=(**

*u**Bu*,

*v*), where

*B*≥1 is the aspect ratio.

Consider the case of the lattice sum *L*_{0}[2,0] (equivalent to −*L*_{0}[0,2] by the identity (3.26)), which will involve the singular part of *u*=0 and *v*=0 in turn, the above may be written as
*u* and *v* are all even. By first performing the summation with respect to *v*,
*B*=1, performing the summation with respect to *u* yields *L*_{0}[2,0]=1/2. It transpires that in the case of square lattices, even more of the lattice sums become trivial, since
*L*_{0}[*m*,*n*] will be zero for any even *m* and *n* whose sum is an odd number multiplied by 2.

As a further example, in the general rectangular lattice case, *L*_{0}[6,0] has the form
*B*=1 here, it transpires that the summation term equals 8/62 and thus *L*_{0}[6,0]=0 for a square lattice, which is consistent with the results derived above. Additionally, in the case of square lattices, for example

## Appendix C. Justifying the form of expansion for centrally symmetric shapes

Recall that central symmetry of the inclusion shape is defined such that *f*(*θ*+*π*)=*f*(*θ*). First, consider the summation terms on the left sides of (3.27) and (3.28). When *δ*+*ξ* is even, including the case *δ*+*ξ*=0 of course, displacement gradient moments of odd order could appear in these terms, but only when *α*+*β* is even. This can be observed by examining the coefficient *π* to 2*π*. By making use of trigonometric sum identities and the fact that *f*(*ψ*+*π*)=*f*(*ψ*) due to central symmetry, ^{δ+ξ+α+β}. Therefore, with *δ*+*ξ* being even, (−1)^{δ+ξ+α+β}=−1 in the case when *α*+*β* is odd and so *δ*+*ξ* being even, only terms within the quadruple sum where *α*+*β* is even will be potentially non-trivial.

With this established, notice next that each moment in the quadruple sum term is multiplied by a lattice sum. It transpires that for every combination of *α*, *β*, *i* and *j* where *i*+*j* is odd and *α*+*β* is even, at least one of the inputs of the lattice sum multiplying *δ*+*ξ* is even).

This leaves the displacement gradient moment expansion of the shape factor (3.12) and (3.13) as the only remaining terms that can potentially involve moments of odd order. Note from the form of its expansion that the shape factor term will not include any such moments if the *i*+*j* is odd (hereon in described as coefficients of odd order). It transpires from the process discussed in appendix A that the coefficients of odd order and those of even order separate into two subsystems. Therefore, if all coefficients of odd order are to be zero, the subsystem of such coefficients should be homogeneous—i.e. no terms independent of the coefficients with indices totalling an odd number should be present in the subsystem.

Firstly note that since *δ*+*ξ* is even, the only inhomogeneous Laplace type condition that one obtains from (A 8) is guaranteed not to involve any odd-ordered coefficients. Therefore, it just remains to be shown that the equations in the system from appendix A related to odd-ordered coefficients are homogeneous. This requires the *m* is odd. When *m* is odd these reduce to
*f* and using trigonometric sum identities, the integral forms for ^{m+δ+ξ}). Since *m* is odd and *δ*+*ξ* is even, 1+(−1)^{m+δ+ξ}=0 and therefore *i*+*j* is odd will be homogeneous in these coefficients and hence the solution for this subset of coefficients will be trivial. This means that, when *δ*+*ξ* is even, the expansion of *J*_{δξ}(** y**) will not feature any powers of the coordinate basis of odd order.

With regard to the shape factor therefore, the above results in no displacement gradient moments of odd order appearing in (3.12) and (3.13) when *δ*+*ξ* is even and the fibre cross-sectional shape is centrally symmetric. Therefore, in this case, no fractional powers of the volume fraction *ϕ* shall appear in the system, which means that asymptotically expanding the moments and *g*(*γ*_{1}) in powers of *ϕ*^{1/2} is unnecessary. Instead, one can expand in powers of *ϕ*.

- Received February 3, 2017.
- Accepted April 3, 2017.

- © 2017 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.