An integral equation method for the homogenization of unidirectional fibre-reinforced media; antiplane elasticity and other potential problems

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.

.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 leadingorder 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. 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.

Introduction
A classical problem in the mechanics of inhomogeneous media is to attempt to replace the two-dimensional potential problem ∇ · (μ(x)∇w(x)) = 0, where x = (x 1 , x 2 ) and μ(x) and w(x) are periodic (scalar) functions that vary rapidly with x, by an equivalent problem of the form ∇ · (μ * · ∇w * (x)) = 0, (1.1) where 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 where μ * ij , (i, j = 1, 2) refers to the ijth component of the tensor in Cartesian coordinates. For orthotropic materials, for example, μ * ij = μ * 1 δ i1 δ j1 + μ * 2 δ i2 δ j2 . The path to (1.1) is the process of 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 σ = μ · e, and 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 ∇ · (μ * · ∇w * (x)) + ρ * ω 2 w * (x) = 0, (1.2) where in the context of antiplane elasticity ρ * is the effective mass density and where timeharmonic 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][6][7], the equivalent inclusion method [8,9], boundary element solutions of integral equations [10] and the use of Fourier transforms [11][12][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.  Table 1. The numerous application areas associated with the potential problem, together with corresponding variables. Here, attention is restricted to two-dimensional problems. It is noted that the acoustic scenario is applicable only in the dynamic setting of course and furthermore (ρ −1 ) ij refers to the ijth component of the inverse of the acoustic density tensor.   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 wavecarrying 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 highfrequency 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, μ * ij = δ 1i δ 1j μ * 1 + δ 2i δ 2j μ * 2 , (i, j = 1, 2) when written with respect to the principal axes of anisotropy. This means that the fibre cross section is restricted to having two planes of reflectional symmetry as indicated in figure 2. In PA, the effective moduli were shown to take the following rational function form at leading order  The coefficients C 1j and D 1j depend upon the shape of the fibre cross section, the periodic lattice geometry and the ratio of fibre to host shear moduli. Upon extending the method to higher orders, for circular cylindrical fibres the following result is derived in this work: and additional general expressions are determined for fibres of more complex cross section. In principle, the scheme presented can be extended to three dimensions and more complex microstructural geometries. 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 μ * 1 and μ * 2 are derived, the method is equally applicable to any of the applications summarized in table 1.

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 R(s, t) = q(sl 1 + tl 2 ), s, t ∈ Z, (2.1) 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 and the periodic cell is a parallelogram with area q 2 R = q 2 A 1 B 2 (e.g. figure 2). Each cell, therefore, consists of a fibre of general cross section occupying the domain 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.
) on the boundary of the fibre. This choice is convenient for calculations that will be carried out in appendix A associated with incorporating the 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 where ∇ = (∂/∂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 c 2 0 = μ 0 /ρ 0 and c 2 I = μ I /ρ I are the squares of the phase velocities of the two phases. The so-called indicator function χ (x) is defined by The associated free-space Green's function for the host domain satisfies where H where ∇ y = (∂/∂y 1 , ∂/∂y 2 ). Non-dimensionalizing using the scalingsx = qx,ŷ = qy andŵ = w/Ŵ, withŴ being a typical displacement magnitude, and noting that Green's function is already non-dimensional, the lattice vector in scaled coordinates becomes p = sl 1 + tl 2 and the crosssectional area of the periodic cell is R = A 1 B 2 . At this point, it appears convenient to define the volume fraction per unit span in the x 3 -direction, φ = |D|/R, where |D| is the (non-dimensional) cross-sectional area of the fibre. Upon dropping the hat notation, the non-dimensional integral equation takes the form where d = ρ I /ρ 0 and m = μ I /μ 0 are the contrasts in mass density and shear moduli, ε = qk 0 and G ε (x − y) = (1/4i)H 0 (ε|x − y|). Note that as ε → 0, where and the volume fraction of the fibre cross section within the periodic cell is then Although expansions are sought with respect to φ 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), In the homogenization regime, 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 η → 0 is then taken. The reader is referred to PA for further details.

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 Now take x ∈ V ab , i.e. within the (a, b)th fibre, which has position vector r = al 1 + bl 2 . By taking the Taylor expansion of Green's function about the point y = p = (p 1 , p 2 ) = sl 1 + tl 2 (i.e. about the centre of the (s, t)th fibre), (3.1) becomes where ∂ i y k denotes the ith derivative with respect to y k and . ij (p) can be thought of as displacement-gradient moments of order i + j, recalling that φ = |D|/R. The order of these moments has been deduced using the fact that, since 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 L δξ [( * )] as the act of multiplying each side of equation ( * ) by (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 L δξ [(3.2)] and subsequently Taylor expand Green's function and its derivatives about x = r to obtain, after some rearrangement where the property ∂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 = r. Furthermore, the following terms have been defined: and where the local polar coordinate system x 1 − r 1 = R cos Θ, x 2 − r 2 = R sin Θ has been defined and where (2.8) is used in order to defineĈ δξ αβ = O(1). As was shown in PA, the influence of the cross-sectional shape of the fibre is embedded solely in terms of the constantsĈ δξ αβ and the shape δξ , the form of which shall be considered shortly.
(a) The shape factor The term A (k) δξ incorporating the fibre cross section in (3.5) appears to possess a singularity at y = x 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 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 = x. Therefore, all that remains to consider are integrals of the type Once again employing the property ∂G ε /∂x k = −∂G ε /∂y k , the x k derivative may be taken inside the y integral, and as the range of integration does not include the region where G(y − x) is singular, exchanging the order of integration is permissable. Retaining the leading-order Green's function as ε → 0, recalling that G 0 was defined in (2.6). It is convenient here to represent (3.10) as a series expansion, i.e.
so that displacement gradient moments naturally arise in (3.9). A procedure for obtaining the coefficients D δξ pq for a given fibre cross section is outlined in appendix A. The order of truncation, P + 2 is governed by the shape function f involved. For elliptical fibres P = δ + ξ is an exact result for example. For arbitrarily shaped cross sections P → ∞ generally. It is noted that second derivatives of (3.10) make up the components of the generalized Hill tensor [23,24]. Its representation via a polynomial expansion as in (3.11) will be discussed in a forthcoming article by the authors.
Substituting (3.11) into (3.9) and adjusting the indices one obtains Note that D δξ 00 , D δξ 10 and D δξ 01 do not arise in the shape factor as they do not survive secondorder differentiation. In the case of fibres with circular cross sections, f (θ) = 1 and it follows that  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 ε) is a lattice sum, which is discussed in the next section.

(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  Therefore, as |u| → 0. The first term of (3.20) can be turned into an integral in the same manner as in Sec. 3(a) of PA. It transpires that when this step is applied with 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 ε → 0 is when i = j = α = β = 0. Thus, where 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) One can show that the integral in (3.24) can be evaluated via stationary phase. Alternatively, as was discussed at the end of §2 one can add loss to the system. Further details of this approach can be found in the appendix of PA.
If one wishes to determine the effective wavenumber in the x 2 -direction (seeking μ * 2 instead of μ * 1 ), it is straightforward to rotate the material by π/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 (d) The asymptotic system in φ Using (3.23) in (3.16) and (3.17) the leading-order system, with respect to ε, iŝ This is an eigenvalue problem with eigenvalues γ 1 and associated eigenvectors comprising the momentsŴ (k) δξ , noting that γ 1 appears only in the term g(γ 1 ). To make progress take expansions in the volume fraction, posingŴ (1) where the form for g is motivated by (3.25) and the fact that γ 1 → 1 as φ → 0. 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, W (k) 10 (r) = O(φ 3/2 ). This would suggest that, in general, powers of φ 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 γ 2 1 = 1/μ * 1 , (3.31) gives Note the form of this expression, i.e. every coefficient at each order in the polynomials in the numerator and denominator is the same except for the O(φ) term. The concern of the next section is the determination of the coefficients h j from the linear system.

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)  in (2.8) and writing without loss of generality that u 0 00 = 1, it is determined that .
where 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 and For the specific case of fibres arranged on a square lattice, one finds that (4.5) Details of how higher-order terms are derived will be given explicitly in the next section for the case of circular cross sections.

(b) Circular cylindrical fibres
To illustrate specific details of the derivation of higher-order terms, take the circular crosssectional 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Ĉ 00αβ coefficients in the quadruple sum, using the relations This illustrates how higher-order moment terms arise, in this case u 0 20 , u 0 02 and v 0 11 , and therefore motivates which equations need to be studied subsequently. In this case, it means that the (11), (20) and (02) equations must be considered at leading order with respect 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)   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 Solving these equations gives 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 Note that the sixth-order lattice sum is zero for square lattices, hence 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 This is consistent with the form of higher-order terms outlined in the conclusion of PA.

(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.

Results
The implementation of the above methodology is now described for a range of geometries in the case of a shear contrast of m  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]) where the tensor components of H may be determined as Here, 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 μ * 1 = μ * 2 . Generally for fibres of non-circular cross section, the finite-element method (FEM) is employed to solve the cell problem. Indeed here, COMSOL multiphysics is used to solve the cell problem for antiplane shear. The methods are compared by plotting the components 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 n to refer to the rational function form (3.32) when order 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 R = 1. Results for the MAH when using the multipole method [1] and the FEM are also plotted to illustrate that the FEM remains consistent with the multipole method. Agreement between the MAH and the IEM estimates is, in general, excellent until the radius begins to approach 0.5 when fibres become very close. Figure 3b 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. 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 4a 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 nearneighbour interaction. Figure 4b 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 5a, 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 5b illustrates results for the specific effective shear modulus μ * 1 for all three fibre shapes considered thus far (circular, elliptical and rectangular), when plotted against volume fraction, as opposed to fibre 'radius'. Figure 6 plots both effective shear moduli μ * 1 and μ * 2 against fibre 'radius'. As well as returning the results to the context of the actual physical-effective properties, these plots illustrate the extra sensitivity that the 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 5a the IEM order 6 expression is able to accurately replicate the upturn as r → 1/2. Furthermore, in the limit as r → 1/2 this medium becomes a layered structure. In this limit, there are exact expressions for the two shear moduli, which are [25,26] for a volume fraction φ of layered material-here this is clearly φ = 1/2 and these results for μ * 1 and μ * 2 at φ = 1/2 are indicated by the black squares in figure 6a,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 7b. 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 = 1, and referring to (2.2), the radius r here is chosen as max θ∈ [0,2π]   Of great utility when describing general-shaped fibres is the so-called superellipse function defined as [27] In the hexagonal case above, 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.

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 cellproblem 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 In general, it is clearly not possible to obtain an analytical form for 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 where m ≥ 1 in the latter. One then finds that where integration by parts was used for the ζ integration and a different notation has been associated with this integral form so that it may be referred to later. Also, where Equating these with the relations linking a δξ n , b δξ n (evaluated at ρ = 1) to the coefficients D δξ pq , such as (A 3), one obtains a linear system in D δξ pq . However, it transpires that more conditions are needed in order to close the system. These are obtained by applying the Laplacian operator to the forms of J δξ in (3.10) and (3.11), equating and exploiting the fact that ∇ 2 G 0 = δ(y − x), to find Comparing coefficients of (y 1 − r 1 ) i (y 2 − r 2 ) j provides the further equations relating the coefficients D δξ pq . Some thought leads to the conclusion that for p, q ≥ 2 This additional set of equations closes the linear system. (a) Linear system The linear system consists of (A 8) and equations, such as e.g. It may be written in the form and the vector 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: P = 0, δ = ξ = 0. The above expressions clarify better why P = 0 is referred to as the second-order case: the sum of the lower indices of the D constants never exceed 2. In this case, the only non-trivial coefficients that have an influence on the shape factor term are and D 00 02 = a (a + b) , so, using (3.12) and (3.13), the shape factor for this system has the form (ii) Circular fibre fourth-order case example: P = 2, δ = 2, ξ = 0 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 aŝ 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. 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 ∂ 2 u 1 G 0 (u) and is Employing this in the lattice sum and considering a rectangular lattice, this gives By isolating the parts of the sum where u = 0 and v = 0 in turn, the above may be written as where sums from −∞ to ∞ excluding 0 are replaced by twice the sum from 1 to ∞ since the powers of u and v are all even. By first performing the summation with respect to v, Note that in the special case of square lattices, 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 Thus via the identity (3.26), L 0 [m, n] will be zero for any even m and n whose sum is an odd number multiplied by 2. [cosech 2 (π Bu)(2 coth 4 (π Bu) + 11 coth 2 (π Bu)cosech 2 (π Bu) + 2cosech 4 (π Bu))] .
When 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Ĉ δξ αβ as defined in (3.7), dividing the range of integration into two parts-one forθ ∈ [0, π ) and the other forθ ∈ [π , 2π ), and using the substitutionθ = ψ + π in the integral from π to 2π . By making use of trigonometric sum identities and the fact that f (ψ + π ) = f (ψ) due to central symmetry,Ĉ δξ αβ can be transformed into an integral multiplied by the factor 1 + (−1) δ+ξ +α+β . Therefore, with δ + ξ being even, (−1) δ+ξ +α+β = −1 in the case when α + β is odd and soĈ δξ αβ = 0. Then, for centrally symmetric fibre cross sections with δ + ξ 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 α + β  Table 2. List of possible scenarios where both α + β is even and i + j is odd and how each scenario causes the lattice sum terms in (3.27) and (3.28) to be trivial. is even, at least one of the inputs of the lattice sum multiplying W (k) ij is odd. All such scenarios are outlined in table 2. As established earlier, the lattice sum evaluates to zero if either of its inputs is odd. Therefore, no moments of odd order will appear in the quadruples sum terms in even-ordered equations (i.e. when δ + ξ 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 D δξ ij coefficients are zero when 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 A As was the case when examiningĈ δξ αβ , splitting the range of integration into two parts-one forθ ∈ [0, π ) and the other forθ ∈ [π , 2π ) then substitutingθ = ψ + π into the latter and making further use of the central symmetry of f and using trigonometric sum identities, the integral forms forÂ δξ m andB δξ m can be rewritten as single integrals multiplied by the term (1 + (−1) m+δ+ξ ). Since m is odd and δ + ξ is even, 1 + (−1) m+δ+ξ = 0 and therefore A δξ m and B δξ m are zero. Therefore, any subsystem of equations in coefficients D δξ ij where 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 crosssectional 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 φ.