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. A464, 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.

Introduction

A classical problem in the mechanics of inhomogeneous media is to attempt to replace the two-dimensional potential problem ∇⋅(μ()∇w())=0, where =(x1,x2) and μ() and w() are periodic (scalar) functions that vary rapidly with , by an equivalent problem of the form 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 refers to the ijth component of the tensor in Cartesian coordinates. For orthotropic materials, for example, . 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 =⋅, and =∇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.

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) refers to the ijth component of the inverse of the acoustic density tensor.

application	σi	ei	w	μij
antiplane elasticity	antiplane stress vector (σ13,σ23)	displacement gradient ∇w	displacement w	shear moduli μij
thermal conductivity	heat flux qi	temperature gradient −∇T	temperature T	thermal conductivity kij
electrical conductivity	electrical current Ji	electric field Ei=∇Φ	electric potential Φ	electrical conductivity σ¯_ij
dielectrics	displacement field Di	electric field Ei=∇Φ	electric potential Φ	electric permittivity ϵij
magnetism	magnetic induction Bi	magnetic field Hi=∇Ψ	magnetic potential Ψ	magnetic permeability μij
porous media	weighted velocity ηvi	pressure gradient ∇p	pressure p	permeability kij
diffusion	diffusion flux ji	concentration gradient ∇c	concentration c	diffusivity Dij
acoustics	acceleration ∂v/∂t	pressure gradient −∇p	pressure p	inverse density tensor (ρ−1)ij

Assume now that μ() 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.

Figure 1.

Cross section of an inhomogeneous medium with piecewise constant material property μ(). The unidirectional fibre-reinforced composite medium depicted here has general periodic microstructure; in particular, the fibres have non-circular cross sections.

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 where in the context of antiplane elasticity ρ* is the effective mass density and where time-harmonic motion e−i 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, 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 and 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.

Figure 2.

Figure depicting the dimensional (a) and non-dimensional (b) periodic cell (and associated lengthscales) in the case of rectangular periodicity and fibre with cross section such that the material is at most macroscopically orthotropic.

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 and 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 x3-coordinate running parallel to the fibre axis and the x1x2-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 1,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 (e.g. figure 2). Each cell, therefore, consists of a fibre of general cross section occupying the domain V embedded within the host phase which is denoted by V0. The shear modulus and mass density of the host (fibre) is denoted as μ0 (μ) and ρ0 (ρ), 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. with f(θ)≥1 and for some θ∈[0,2π], f(θ)=1. Hence, 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 x3-direction is denoted by w(), where it is recalled that =(x1,x2). The elastic displacement w is governed by where ∇=(∂/∂x1,∂/∂x2) and where k0=ω/c0 and k=ω/c are the wavenumbers of the host and fibre, respectively, and where and are the squares of the phase velocities of the two phases. The so-called indicator function χ() is defined by The associated free-space Green’s function for the host domain satisfies where denotes the Hankel function of the first kind and zeroth order.

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 where ∇=(∂/∂y1,∂/∂y2). Non-dimensionalizing using the scalings and , with being a typical displacement magnitude, and noting that Green’s function is already non-dimensional, the lattice vector in scaled coordinates becomes =s1+t2 and the cross-sectional area of the periodic cell is . At this point, it appears convenient to define the volume fraction per unit span in the x3-direction, , 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=ρ/ρ0 and m=μ/μ0 are the contrasts in mass density and shear moduli, ε=qk0 and G(−)=(1/4i)H0(ε|−|). Note that as , where and γ≃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 and the volume fraction of the fibre cross section within the periodic cell is then Therefore, 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, ak0≪1 and ε=qk0≪1 and it can be assumed here that ak0=O(ε). In [21,22], the regime where ak0≪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 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=1,2 to yield Now take ∈V, i.e. within the (a, b)th fibre, which has position vector =a1+b2. By taking the Taylor expansion of Green’s function about the point ==(p1, p2)=s1+t2 (i.e. about the centre of the (s, t)th fibre), (3.1) becomes where denotes the ith derivative with respect to y and The variables introduced as can be thought of as displacement-gradient moments of order i+j, recalling that . 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 since 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 as the act of multiplying each side of equation (*) by (x1−r1)(x2−r2)/(δ!ξ!) and integrating in the plane over the domain V, where =(r1,r2). Hence, consider [(3.2)] and subsequently Taylor expand Green’s function and its derivatives about = to obtain, after some rearrangement where the property ∂G/∂x=−∂G/∂y has been employed. Here, the outermost summation notation refers to summing over every fibre location =(p1,p2) except for =. Furthermore, the following terms have been defined: and where the local polar coordinate system has been defined and where (2.8) is used in order to define . 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 factors , the form of which shall be considered shortly.

The shape factor

The term incorporating the fibre cross section in (3.5) appears to possess a singularity at = 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 up into a non-singular part V∖C and apparently singular part C, where C is a disc of radius ψ≪1 with origin =. Therefore, all that remains to consider are integrals of the type Once again employing the property ∂G/∂x=−∂G/∂y, the x derivative may be taken inside the integral, and as the range of integration does not include the region where G(−) is singular, exchanging the order of integration is permissable. Retaining the leading-order Green’s function as , where recalling that G0 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 for a given fibre cross section is outlined in appendix A. The order of truncation, is governed by the shape function f involved. For elliptical fibres is an exact result for example. For arbitrarily shaped cross sections 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 and Note that , and do not arise in the shape factor as they do not survive second-order differentiation. In the case of fibres with circular cross sections, f(θ)=1 and it follows that if δ+ξ+α+β is odd. As , equations (3.12)–(3.13) lead to the conclusion that .

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 . A wave-like ansatz for these moments is now posed.

Wave-like solutions

Noting (2.9), plane wave type solutions of (3.4) are sought of the form with being the non-dimensional effective wavenumber (scaled on q) in the direction of the angle subtended from the x1-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 ε) and where is a lattice sum, which is discussed in the next section.

Lattice sums

First pick out the singular, non-integrable behaviour of the derivatives of Green’s function in the lattice sum (3.18), writing identifying S[m,n] as the singular part of the derivative. Then for a given m,n, (3.18) can be written as with where Therefore, with as . 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 ε. Therefore, the only term from this step that remains in the limit is when i=j=α=β=0. Thus, where recalling that (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 x1-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 x2-direction (seeking instead of ), 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) x1-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 L0[m,n] straightforwardly. Firstly, since the singular part of G0 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 Furthermore, only cases where both m and n are even will give a non-zero L0[m,n]. For example, consider the case of L0[1,2]=−L0[3,0] due to (3.26). The associated S0[m,n] takes the form When this is employed in the lattice sum (3.22), the odd powers of u1 in the numerator ensure that the summation is zero. This same reasoning gives L0[m,n]=0 whenever m+n is odd. It is shown in appendix B how L0[m,n] can be straightforwardly determined when it is non-zero, i.e. when m+n is even.

The asymptotic system in ϕ

Using (3.23) in (3.16) and (3.17) the leading-order system, with respect to ε, is and This is an eigenvalue problem with eigenvalues γ1 and associated eigenvectors comprising the moments , noting that γ1 appears only in the term g(γ1). To make progress take expansions in the volume fraction, posing where the form for g is motivated by (3.25) and the fact that as .

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, . 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 , (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 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 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. Hereon in the term ‘(δξ) equations’ shall refer to (3.27) and (3.28) for fixed δ and ξ. Start by noting that expressions for the h 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, say, the first (00) equation must be considered up to order ϕ. Each order provides an expression for a coefficient h in terms of the eigenvector components , lattice sums, coefficients of the shape factors and possibly components from other moments and . By tracking which terms from the displacement gradient moments appear in this highest-order equation, one can then see which truncation with respect to the choices of indices (δξ) needs to be made and hence which additional (δξ) equations need to be considered in order to solve the linear system for the required terms and and hence subsequently h. The following examples illustrate the implementation of this scheme; orderings now refer to ϕ.

Elliptical cylindrical fibres

Consider elliptical cylindrical fibres with semi-axes a and b in the x1 and x2 directions, respectively. At O(1), noting that , where we recall that τ is defined for a given shape in (2.8) and writing without loss of generality that , it is determined that It transpires that the second (00) equation yields . The O(ϕ) terms give (upon using the properties of the lattice sum) and after using (4.1) this gives while the second (00) equation to O(ϕ) gives . Continuing in this vein, (3.32) takes the form where L0[2,0]=−L0[0,2]=−(S−1)/2, C4 and C6 depend upon the fourth- and sixth-order lattice sums, respectively, and the p 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 For the specific case of fibres arranged on a square lattice, one finds that Details of how higher-order terms are derived will be given explicitly in the next section for the case of circular cross sections.

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 coefficients in the quadruple sum, using the relations (3.26) and the rule involving odd indices within the lattice sum, in addition to the result from the previous section that and the forms (4.1) and (4.2) for h−1 and h0, the equation reduces to This illustrates how higher-order moment terms arise, in this case , and , 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 h1. 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 , and , while the second (20) and (02) equations and the first (11) equation form a system in , and . The latter subsystem is homogeneous with respect to the leading-order coefficients of the moments featured, while (4.6) illustrates how the former contributes to the process of obtaining h1. 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 and hence (4.6) gives h1=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 where . Note that the sixth-order lattice sum is zero for square lattices, hence C6=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.

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=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]) where the tensor components of H may be determined as Here, is the outer unit normal to the fibre boundary ∂V, is the short lengthscale of the problem, and =(N1,N2) is the solution to the associated cell problem. In the case of orthotropy, H12=H21=0, and if the fibre cross section has a rotational symmetry of π/2, H11=H22 and thus . 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 H11 and H22 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 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).

Consider circular fibres of radius r. Figure 3 compares results for H11(r)=H22(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 . 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 3.

Plot of H11(r) for m=18.75 when circular fibres are arranged on a square lattice, with associated solutions of the cell problem from the MAH inset (a) and a magnified plot for large radii (b). (Online version in colour.)

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 H11 since this is the component that is affected by the near-neighbour 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 H11. 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.

Figure 4.

Plot of H11(r) and H22(r), where r is the major axis of the ellipse, for m=18.75 when elliptical fibres of aspect ratio 2 are arranged on a square lattice, with associated solutions of the cell problem from the MAH inset (a) and a magnified plot for large radii (b). (Online version in colour.)

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 H11(r) and H22(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 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 and 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.

Figure 5.

Plot of H11(r) and H22(r) for the situation when for m=18.75 and rectangular fibres of aspect ratio 2 are arranged on a square lattice, with associated solutions of the cell problem from the MAH inset (a) and a plot of the effective property versus volume fraction for circular, rectangular and elliptical fibre cross sections (b). (Online version in colour.)

Figure 6.

Plot of the effective antiplane shear moduli (a) and (b) versus fibre ‘radius’ for the situation when circular, elliptical and rectangular fibres are arranged on a square lattice, m=18.75. The black squares at r=1/2 indicate the layered media results as given in (5.3). (Online version in colour.)

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 . Furthermore, in the limit as 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 and 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.

Figure 7.

Plot of H11(r) and H22(r) for m=18.75 when hexagonal fibres are arranged on a square lattice, with associated solutions of the cell problem from the MAH inset (a) and a close-up for high radii (b). (Online version in colour.)

Figure 8.

Plot of effective antiplane shear moduli and for m=18.75 when hexagonal-shaped fibres are arranged on a square lattice. (Online version in colour.)

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 , and referring to (2.2), the radius r here is chosen as . This ensures that r∈[0,0.5] for all plots.

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, n1=100, m=6 and n2=n3=62. Alternatively for the rectangle, a=2, b=1, m=4 and n1=n2=n3=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 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.

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.

α	β	i	j	reason lattice sums are trivial
even	even	odd	even	i+α, i+α+2 and j+β+1 are all odd
odd	odd	even	odd	i+α, i+α+2 and j+β+1 are all odd
even	even	even	odd	i+α+1, j+β and j+β+2 are all odd
odd	odd	odd	even	i+α+1, j+β and j+β+2 are all odd