Abstract

The effect of a longitudinal stretch and a pressure-induced inhomogeneous radial deformation on the scattering of antiplane elastic waves from a cylindrical cavity is determined. Three popular nonlinear strain energy functions are considered: the neo-Hookean, the Mooney–Rivlin and a two-term Arruda–Boyce model. A new method is developed to analyse and solve the governing wave equations. It exploits their properties to determine an asymptotic solution in the far-field, which is then used to derive a boundary condition to numerically evaluate the equations local to the cavity. This method could be applied to any linear ordinary differential equation whose inhomogeneous coefficients tend to a constant as its independent variable tends to infinity. The effect of the pre-stress is evaluated by considering the scattering cross section. A longitudinal stretch is found to decrease the scattered power emanating from the cavity, whereas a compression increases it. The effect of the pressure difference depends on the strain energy function employed. For a Mooney–Rivlin material, a cavity inflation increases the scattered power and a deflation decreases it; for a neo-Hookean material, the scattering cross section is unaffected by the radial deformation; and for a two-term Arruda–Boyce material, both inflation and deflation are found to decrease the scattered power.

1. Introduction

Canonical wave scattering problems have been of great interest for many decades, with applications in numerous areas, including water waves, electromagnetics, acoustics, seismology and non-destructive evaluation. The objective is to determine the scattered field given information regarding the form of the incident wave and scattering obstacle or inhomogeneity. In particular, the type of boundary conditions imposed on the boundary between the obstacle and the surrounding host medium is key [1,2]. The case of cavities, cracks and other defects in an elastic medium is one of tremendous importance [3–6] due to the fact that their presence frequently leads to a degradation in material performance. The primary aim of non-destructive evaluation is to predict the presence of such flaws in a medium and, although this is always framed as an inverse model, the forward scattering problem must be understood in order to deduce conclusions about the inverse problem [7]. In many cases there may well be a number of such inhomogeneities present in the medium, and so much effort has gone into the prediction of the so-called effective wavenumber of inhomogeneous media. This theory has been developed in the case of flawed media as well as fibre-reinforced or particulate composites, which themselves may possess imperfections associated with their adherence to the host phase (e.g. [8]).

The propagation of elastic waves in an inhomogeneous material depends strongly on the distribution and properties of the inhomogeneities inside the material. At low frequency, the medium generally responds as an effective homogeneous medium with properties defined by the leading order scattering characteristics of the inhomogeneities, i.e. the Rayleigh limit [9–11]. At higher frequencies, complex frequency-dependent behaviour can develop, for example, attenuation in random media [12] and complete reflection in periodic media in the so-called band gaps [13]. Inhomogeneous materials are implicitly required in order to construct metamaterials, media that are generally strongly anisotropic and frequency dependent [14,15]. Such media are generally not found in nature and enable rather surprising effects such as cloaking and negative refraction [16–18].

Recently, nonlinear elastic materials have been used in order to construct tunable band gaps that can be shifted in real time [19–22]. This work employs the theory of small-on-large [23,24] to derive the governing pressure-dependent incremental wave equations by linearizing about the nonlinear equilibrium state. Shearer et al. [25] used this theory to consider the propagation of torsional waves through an inhomogeneously pre-stressed annular cylinder. In the main, however, interest has centred on the influence of homogeneous stretch distributions (and hence the influence of induced anisotropy alone) on subsequent wave propagation (e.g. [26,27]). When the medium in question is inhomogeneous (for example, a fibre-reinforced, or particulate composite, material) and the host phase is nonlinear elastic, pre-stress will almost always lead to non-homogeneous as well as anisotropic stretch distributions, except in very special cases (e.g. [19]).

It is, therefore, of interest to investigate how an initial static pre-stress, which leads to inhomogeneous stress and strain fields in a nonlinear elastic body, influences the propagation and scattering of small-amplitude waves. In this article, we consider a canonical scattering problem, i.e. that of scattering of horizontally polarized shear waves from a cylindrical cavity in an incompressible, pre-stressed, nonlinear elastic medium. Note that a similar problem involving homogeneous initial deformations was studied by Leungvichcharoen & Wijeyewickrema [28]. We extend the work of Parnell & Abrahams [29], in which antiplane wave scattering from a cylindrical cavity in a neo-Hookean material was investigated, by considering more complex constitutive models, namely the Mooney–Rivlin and Arruda–Boyce strain energy functions. Parnell & Abrahams [29] showed that, in the neo-Hookean case, the scattering coefficients are completely unaffected by the application of a hydrostatic pressure in the radial direction along with a longitudinal stretch. Parnell and co-workers [30,31] exploited this result to create a theoretical cloak for antiplane elastic waves. The benefit of this over classical metamaterial cloaks is that inhomogeneous density is not required and the necessary inhomogeneity in the incremental shear modulus is induced naturally by the imposed pre-stress. Norris & Parnell [32] extended the work to the case of in-plane elastic waves, and Parnell & Shearer [15] considered the effect of an imperfect cloaking material characterized by a Mooney–Rivlin strain energy function. The governing equations could not be solved analytically in this latter case; they were solved numerically over the finite domain of the cloak. It is more challenging to solve the governing equations over an infinite domain, the case that shall be considered in this work. The developed methodology can be employed in general for any strain energy function. It also seems to be a plausible scheme for the general solution of ordinary differential equations with inhomogeneous coefficients.

Because scattering coefficients associated with single, canonical scattering problems are employed in multiple scattering models [8,11], it is envisaged that this work will be useful for the prediction of effective wave propagation through pre-stressed nonlinear elastic inhomogeneous media (e.g. soft porous materials), certainly for the case of dilute dispersions.

In §2, we describe the finite deformation that ensues within a nonlinear elastic material when a far-field hydrostatic pressure, along with a longitudinal stretch, is imposed. In §3, we then carry out the small-on-large analysis, assuming that the incremental field is a horizontally polarized (antiplane) shear wave and we derive the governing equation for the antiplane displacement field in the pre-stressed medium. In §4, we describe a novel hybrid analytical/numerical method used to solve the governing equations over an infinite domain and use it to determine, in §5, the scattering cross sections of the waves in neo-Hookean, Mooney–Rivlin and Arruda–Boyce materials subjected to various levels of pre-stress. Concluding remarks are offered in §6.

2. Initial finite static deformation (pre-stress)

Consider an isotropic, incompressible, nonlinear elastic material of infinite extent with a circular cylindrical cavity, of initial radius A, whose axis is parallel to the Z axis of a Cartesian coordinate system (X,Y,Z) with its centre in the (X,Y) plane located at the origin. The constitutive behaviour of the hyperelastic material is described by the strain energy function W=W(I1,I2) [23], where Ij, j=1,2,3 are the principal strain invariants, and we note that W is independent of I3 due to the constraint of incompressibility (i.e. I3=1 for all deformations).

The medium is deformed due to an imposed far-field hydrostatic pressure, an internal hydrostatic pressure and a uniform axial stretch, as depicted in figure 1. The symmetry of geometry and forcing therefore implies that the deformation is described by
R=R(r),Θ=θandZ=zζ,2.1
where (R,Θ,Z) and (r,θ,z) are cylindrical polar coordinates (X=Rcos⁡Θ, Y=Rsin⁡Θ, x=rcos⁡θ, y=rsin⁡θ) in the undeformed and deformed configurations, respectively, and R(r) is a function that is determined from the radial equation of equilibrium and incompressibility condition. Note the convention introduced in (2.1), i.e. that upper case variables correspond to the undeformed configuration, whereas lower case variables correspond to the deformed configuration. We are interested in incremental perturbations from the initial statically deformed state and, therefore, it will be convenient to derive equations in terms of coordinates in the deformed configuration. Position vectors in the undeformed and deformed configurations are, respectively,
X=RER(Θ)+ZEZ=R(r)ER(Θ)+zζEZ,x=rer(θ)+zez,2.2
where ER and EZ are the radial and longitudinal basis vectors in the undeformed configuration, and er and ez are the corresponding basis vectors in the deformed configuration.

(a) The undeformed cavity, with radius A. (b) The deformed cavity, with radius a. The deformation is due to an imposed far-field hydrostatic pressure p∞ and an internal pressure pa, along with an axial stretch ζ.

Using (2.1), it can be shown that the principal stretches for this deformation in the radial, azimuthal and longitudinal directions, respectively, are
λr=drdR=1R′(r),λθ=rR(r)andλz=ζ.2.3
The deformation gradient tensor F is given by
F=Gradx,2.4
where Grad represents the gradient operator in the undeformed configuration. In our case, we have
F=FiJei⊗EJ,FiJ=(λr(r)000λθ(r)000λz)=((R′(r))−1000rR(r)000ζ).2.5
For an incompressible material J=I3=detF=1 and so
λrλθλz=ζrR(r)R′(r)=1.2.6
This is easily solved to yield
R(r)=ζ(r2+M),2.7
where M is a constant defined by
M=A2ζ−a2.2.8

From Ogden [23,24], the Cauchy and nominal stress tensors for an incompressible material are, respectively, given, in terms of the deformation gradient and strain energy function, by
T=F∂W∂F+QIandS=∂W∂F+QF−1,2.9
where I is the identity tensor, Q is a Lagrange multiplier associated with the incompressibility constraint, often referred to as an arbitrary hydrostatic pressure, and we note that the differentiation of W with respect to F is defined component-wise as follows:
(∂W∂F)Ij=∂W∂FjI.2.10
The static equations of equilibrium are then given by
divT=0,2.11
where div signifies the divergence operator with respect to x. These reduce to
∂Trr∂r+1r(Trr−Tθθ)=0,∂Tθθ∂θ=0and∂Tzz∂z=0,2.12
where T=Tijei⊗ej. The second and third of these imply Q=Q(r). Integrating the first of (2.12), using (2.9) and applying Trr|r=a=−pa, we find that the radial stress is given by
Trr(r)+pa=∫ar1s(λθ(s)∂W∂λθ−λr(s)∂W∂λr)ds,2.13
then by applying Trr→−p∞ as r→∞, we obtain an equation relating the pressure difference to the resulting deformation,
pa−p∞=∫a∞1s(λθ(s)∂W∂λθ−λr(s)∂W∂λr)ds.2.14
To proceed, we must select a strain energy function W to substitute into equation (2.14). We shall use three different strain energy functions, the neo-Hookean:
WNH=μ2(I1−3)=μ2(λr2+λθ2+λz2−3),2.15
where μ is the ground-state shear modulus of the material under consideration, the Mooney–Rivlin:
WMR=μ2(C1(I1−3)+C2(I2−3))=μ2(C1(λr2+λθ2+λz2−3)+C2(λr2λθ2+λr2λz2+λθ2λz2−3)),2.16
where C1 and C2 are material constants that sum to unity, and the Arruda–Boyce model. The Arruda–Boyce strain energy function can be expressed as an infinite series in terms of powers of I1 [33]. For simplicity, we shall use the first two terms only:
WAB=5Nμ3+5N(12(I1−3)+120N(I12−9))=5Nμ3+5N(12(λr2+λθ2+λz2−3)+120N(λr4+λθ4+λz4+2(λr2λθ2+λr2λz2+λθ2λz2)−9)),2.17
where N is the assumed number of links in the polymer chains that form the molecular basis of rubber-like materials and we note that the coefficient in front of the shear modulus is necessary for this two-term model to be consistent with linear elasticity.

The resulting pressure–deformation relationship in the neo-Hookean case is
p∞−paμ=12ζ(log(A2ζa2)+A2ζa2−1),2.18
in the Mooney–Rivlin case is
p∞−paμ=C1+ζ2C22ζ(log(A2ζa2)+A2ζa2−1),2.19
and in the Arruda–Boyce case is
p∞−paμ=14(3+5N)ζ4(1−ζa2A2)(A4a4+2ζ2+A2a2ζ(3+10Nζ+2ζ3))+1+5Nζ+ζ32(3+5N)log(A2ζa2).2.20
In figure 2, we plot a/A as a function of (p∞−pa)/μ for each of the strain energy functions listed above, and for three separate values of ζ. The black curves represent the neo-Hookean strain energy function, the blue curves represent the Mooney–Rivlin strain energy function and the red curves represent the Arruda–Boyce strain energy function. The solid lines correspond to ζ=1, the dashed to ζ=12 and the dotted to ζ=2. The values chosen for the constants in the Mooney–Rivlin strain energy function were C1=0.724, C2=0.276, which were the values used by Mooney [34] to fit Gerke's [35] experimental data on the elongation of vulcanized rubber. The value chosen for N in the Arruda–Boyce strain energy function was N=26.5. This value of N was used by Arruda & Boyce [33] to match Treloar's [36] data on the uniaxial extension, biaxial extension and shear of vulcanized rubber. Note that, for ζ=1, the neo-Hookean and Mooney–Rivlin curves coincide and the Arruda–Boyce curve is so close to them that it also appears to lie on top of them in the figure.

a/A as a function of (p∞−pa)/μ. The black curves represent the neo-Hookean strain energy function, the blue curves represent the Mooney–Rivlin strain energy function and the red curves represent the Arruda–Boyce strain energy function. The solid lines correspond to ζ=1, the dashed to ζ=12 and the dotted to ζ=2. The values chosen for the constants in the Mooney–Rivlin strain energy function were C1=0.724, C2=0.276, and the value chosen for N in the Arruda–Boyce strain energy function was N=26.5.

3. Incremental deformations

We now examine incremental deformations from the deformed state b. To this end, consider a finite deformation of the original body B to a new deformed state b¯ which is close to the configuration b. The position vector in the new deformed state b¯ is given by x¯ and we define
u=x¯−x,|u|≪13.1
as the difference between position vectors in b¯ and b. Since b¯ is close to b, u is called an incremental displacement which we assume to be time dependent. For the purposes of this article we assume that the incremental deformation is of antiplane type and time-harmonic, i.e.
u=ℜ[w(r,θ)exp⁡(−iωt)]ez,3.2
where the notation ℜ indicates that the real part of the expression inside the square brackets is to be taken. We next use the theory of small-on-large (see appendix A), to determine the governing equation of motion,
1r∂∂r(rμr(r)∂w∂r)+μθ(r)r2∂2w∂θ2+ρω2w=0,3.3
where
μr(r)=(λr(∂W/∂λr)−λz(∂W/∂λz)λr2−λz2)λr2andμθ(r)=(λθ(∂W/∂λθ)−λz(∂W/∂λz)λθ2−λz2)λθ2.}3.4
Given an incident field that takes plane wave form in the far-field, i.e.
wi∼eikx=∑n=0∞εninJn(kr)cos⁡(nθ)3.5
as r→∞, where Jn is the Bessel function of the first kind of order n and
εn={1n=0,2n≠0,3.6
we seek solutions to (3.3) in the form
w(r,θ)=∑n=0∞εninfn(r)cos⁡(nθ).3.7
As such for the nth mode, we have
1r∂∂r(rμr(r)fn′(r))+(ρω2−μθ(r)r2n2)fn(r)=0,3.8
in which prime notation denotes differentiation with respect to r. Clearly different incident fields could be considered as desired. By substituting (2.15)–(2.17) into (3.4), we can obtain explicit expressions for the anisotropic shear moduli μr(r) and μθ(r), which can then be used to determine the antiplane wave governing equation for each strain energy function. In the neo-Hookean case, we have
(1+Mr2)fn′′(r)+1r(1−Mr2)fn′(r)+(kNH2−n2r2+M)fn(r)=0,3.9
where
kNH2=ζK2andK2=ρω2μ,3.10
in the Mooney–Rivlin case we have
(1+mr2)fn′′(r)+1r(1−mr2)fn′(r)+(kMR2−n2r2(1−mr2+M))fn(r)=0,3.11
where
kMR2=ζ2K2T,m=MζC1T,T=1+(ζ−1)C1,3.12
and in the Arruda–Boyce case we have
(1+Mr2+χM2r4)fn′′(r)+1r(1−Mr2−3χM2r4)fn′(r)+(kAB2−n2(1r2+M+χM2r2(r2+M)2))fn(r)=0,3.13
where
kAB2=(3+5N)ζ2χK2andχ=12+5Nζ+ζ3,3.14
and we note that, when ζ=1, kNH=kMR=kAB=K.

By substituting (3.7) into (A 20), we see that the boundary condition on r=a reduces to
fn′(a)=0,3.15
for all n. This boundary condition remains the same regardless of the strain energy function used to characterize the material under consideration.

We note that the effect of the radial pre-stress on the governing wave equations (3.9)–(3.13) decreases as r→∞, and in the limit as M/r2→0, (3.9), (3.11) and (3.13) become
fn′′(r)+1rfn′(r)+(k2−n2r2)fn(r)=0,3.16
where k2=k2NH in the neo-Hookean case, k2=k2MR in the Mooney–Rivlin case, and k2=k2AB in the Arruda–Boyce case. Equation (3.16) is the standard governing equation for antiplane waves in a stress-free elastic medium (albeit with a modified wavenumber) and is readily solved:
fn(r)=fni(r)+fns(r),fni(r)=Jn(kr)andfns(r)=anHn(1)(kr),3.17
where we identify fni(r) as the nth component of the incoming plane wave and thus equate it to Jn via (3.5). Furthermore, fns(r) is the nth component of the scattered wave, represented by a Hankel function of the first kind, Hn(1)(kr), and an is the scattering coefficient associated with this nth outgoing wave term.

Owing to the complexity of the above incremental governing equations, they cannot be solved via standard methods, except for a neo-Hookean material (see §3a). In this case, exact solutions can be found [29]; however, in the Mooney–Rivlin (except when n=0, see §3b) and Arruda–Boyce cases, we must instead find an approximate solution. The method used to do this is described in §4.

(a) The special case of neo-Hookean media

Parnell & Abrahams [29] showed that (3.9) can be solved analytically to give
fn(r)=Jn(Kζ(r2+M))+anHn(1)(Kζ(r2+M))3.18
for the case of an incoming plane wave, where the scattering coefficients an are given by
an=−Jn′(Kζ(a2+M))Hn(1)′(Kζ(a2+M))=−Jn′(KA)Hn(1)′(KA).3.19
This is the form the scattering coefficients take in a stress-free medium containing a cavity of radius A, i.e. the radius of the undeformed cavity in the problem being considered in this paper. Therefore, in the neo-Hookean case, rather surprisingly, all the scattering coefficients are completely unaffected by the pre-stress.

(b) Analytical solution in the Mooney–Rivlin case when n=0

In the Mooney–Rivlin case, to the authors’ knowledge, a general solution for all mode numbers n cannot be found; however, we note that, when n=0, (3.11) reduces to
(1+mr2)f0′′(r)+1r(1−mr2)f0′(r)+k2f0(r)=0,3.20
which can be solved analytically (as for the neo-Hookean case) to give
f0(r)=c1J0(kr2+m)+c2H0(1)(kr2+m),3.21
where c1 and c2 are arbitrary constants. Upon matching this solution to the far-field form (3.17), we obtain c1=1 and c2=a0, and applying the boundary condition (3.15), we obtain
a0=−J0′(ka2+m)H0(1)′(ka2+m).3.22
This is a similar form to that which the neo-Hookean scattering coefficients take when n=0; however, since ka2+m≠KA, the above expression is dependent on the pre-stress. To determine the effect of the pre-stress on the other modes (n≠0), however, we must use the method described in the following section.

4. Solution method

As mentioned above, the Mooney–Rivlin (3.11) and Arruda–Boyce (3.13) governing equations are not readily solved via standard methods. Hence, to proceed, we recall that, in the limit as M/r2→0, an analytical solution (3.17) to the governing equations may be found. Therefore, for large r, we should be able to derive an approximate solution. For a given pre-stress parameter M, we shall choose a large radius b such that M/b2=ϵ2≪1 and derive an approximate solution to the governing equation for r≥b. Then, in the region r∈[a,b], we shall solve the governing equations numerically, with the boundary conditions on r=b being dictated by enforcing continuity of displacement and traction with the approximate solution. In what follows we describe a ‘general’ recipe for numerically evaluating the solution in the far-field to O(ϵ2), which will work for any strain energy function; for ease of exposition this is presented for a Mooney–Rivlin material, but the Arruda–Boyce case follows in a very similar fashion. Note, however, in the Mooney–Rivlin case we can in fact evaluate the integral expressions below analytically, as will be discussed at the end of this section. We begin by introducing new independent and dependent variables s=r/b and Fn(s)=fn(r), respectively, so that (3.11) can be rewritten as
(1+mMϵ2s2)d2Fnds2+1s(1−mMϵ2s2)dFnds+(κ2−n2s2(1−(m/M)ϵ2s2+ϵ2))Fn=0,4.1
where κ2=(kb)2. Given (3.17), we then assume the following (regular) expansion for Fn(s):
Fn(s)=Jn(κs)+anHn(1)(κs)+ϵ2Gn(s)+O(ϵ4),4.2
so that at O(1) we obtain
d2ds2(Jn(κs))+1sdds(Jn(κs))+(κ2−n2s2)Jn(κs)+an(d2ds2(Hn(1)(κs))+1sdds(Hn(1)(κs))+(κ2−n2s2)Hn(1)(κs))=0,4.3
which is satisfied trivially, and at O(ϵ2) we obtain
d2Gnds2+1sdGnds+(κ2−n2s2)Gn=mM(2s3dds(Jn(κs)+anHn(1)(κs))+(κ2s2−2n2s4)(Jn(κs)+anHn(1)(κs))).4.4
Now, if we let
Gn(s)=mMHn(1)(κs)gn(s),4.5
then (4.4) simplifies to
dds(s(Hn(1)(κs))2dgnds)=2s2Hn(1)(κs)dds(Jn(κs))+(κ2s−2n2s3)Hn(1)(κs)Jn(κs)+an(2s2Hn(1)(κs)dds(Hn(1)(κs))+(κ2s−2n2s3)(Hn(1)(κs))2).4.6
In the limit as s→∞, the governing equation (4.1) reduces to Bessel's equation (since the terms multiplying m/M disappear), whose solution is Fn(s)=Jn(κs)+anHn(1)(κs); hence Gn and gn must tend to 0 as s→∞. The next step of the method, therefore, is to integrate (4.6) from s to ∞, and divide by s(Hn(1)(κs))2, to obtain
dgnds=−1s(Hn(1)(κs))2∫s∞(2x2Hn(1)(κx)ddx(Jn(κx))+(κ2x−2n2x3)Hn(1)(κx)Jn(κx))dx+ans3−ans(Hn(1)(κs))2∫s∞(κ2x+2(1−n2)x3)(Hn(1)(κx))2dx,4.7
and, upon integrating again, we find
gn(s)=∫s∞1y(Hn(1)(κy))2∫y∞(2x2Hn(1)(κx)ddx(Jn(κx))+(κ2x−2n2x3)Hn(1)(κx)Jn(κx))dxdy−an2s2+an∫s∞1y(Hn(1)(κy))2∫y∞(κ2x+2(1−n2)x3)(Hn(1)(κx))2dxdy.4.8
So,
Fn(s)=Jn(κs)+anHn(1)(κs)+mMϵ2Hn(1)(κs)gn(s)+O(ϵ4),4.9
and, therefore,
fn(r)=Jn(kr)+anHn(1)(kr)+mMϵ2Hn(1)(kr)gn(rb)+O(ϵ4),4.10
and, hence,
fn(b)=Jn(kb)+anHn(1)(kb)+mb2Hn(1)(kb)gn(1)+O(ϵ4)4.11
and
dfndr|r=b=ddr(Jn(kr))|r=b+anddr(Hn(1)(kr))|r=b+mb3Hn(1)(kb)dgnds|s=1+mb2ddr(Hn(1)(kr))|r=bgn(1)+O(ϵ4).4.12
Now, in order to determine the appropriate boundary condition on r=b, we require a suitable choice for an. We note that we can write
fn(b)=αan+β,4.13
where
α=Hn(1)(kb)(1−m2b2+mb2∫1∞1y(Hn(1)(κy))2∫y∞(κ2x+2(1−n2)x3)(Hn(1)(κx))2dxdy),4.14β=Jn(kb)+mb2Hn(1)(kb)∫1∞1y(Hn(1)(κy))2×∫y∞(2x2Hn(1)(κx)ddx(Jn(κx))+(κ2x−2n2x3)Hn(1)(κx)Jn(κx))dxdy4.15andfn′(b)=dfndr|r=b=γan+δ,4.16
where
γ=ddr(Hn(1)(kr))|r=bαHn(1)(kb)+mb3Hn(1)(kb)×(1−1(Hn(1)(kb))2∫1∞(κ2x+2(1−n2)x3)(Hn(1)(κx))2dx)4.17
and
δ=ddr(Jn(kr))|r=b−ddr(Hn(1)(kr))|r=bJn(kb)Hn(1)(kb)+ddr(Hn(1)(kr))|r=bβHn(1)(kb)−mb3Hn(1)(kb)∫1∞(2x2Hn(1)(κx)ddx(Jn(κx))+(κ2x−2n2x3)Hn(1)(κx)Jn(κx))dx.4.18
Solving (4.13) and (4.14) for an, we obtain
an=1Δ(δfn(b)−βfn′(b)),4.19
where
Δ=αδ−βγ,4.20
and eliminating an from (4.13) and (4.16), we obtain our final boundary condition,
γfn(b)−αfn′(b)=−Δ.4.21
The method is now clear. We choose b such that M/b2≪1, and numerically solve equation (3.11) subject to (3.15) and (4.21). Once solved, we use the values of fn(b) and fn′(b) determined by our numerical solver in equation (4.19) to deduce the value of an. The predicted value of an will depend on the chosen value of b, therefore we must increase b until the solution has converged to a desired level of accuracy.

We note that, owing to their highly oscillatory nature, the integrals in equations (4.14)–(4.18) are difficult to evaluate numerically; however, in appendix B, we show how they can be rearranged into forms that are more readily evaluated. A similar procedure can be followed to obtain such integrals for other strain energy functions. In fact, for Mooney-Rivlin materials it transpires that the integrals can be evaluated exactly, leading to an explicit expression for the O(ϵ2) term:
Gn(s)=mM12sdds(Jn(κs)+anHn(1)(κs)).4.22
This has allowed us to verify that the numerical procedure is effective.

To demonstrate how the method described above converges with increasing b, in figure 3 we plot the fundamental scattering coefficient a0 in a pre-stressed neo-Hookean material with M=1, ζ=1, K=1 and a=1. The black line gives the values derived using the method above (in which we determined the far-field solution up to O(ϵ2)), the red dashed line shows the analytical solution (which is available since we are considering a neo-Hookean material), and the blue line represents the solution that is obtained when one neglects the O(ϵ2) correction to fn(r) (see (4.2)) in the far-field. We observe that the inclusion of the O(ϵ2) terms greatly improves the convergence. We discuss the convergence of the method further in appendix C.

The fundamental scattering coefficient a0 in a pre-stressed neo-Hookean material with M=1, ζ=1, K=1 and a=1, as a function of b, derived using the method described in §4 (black), the analytic expression (red, dashed) and a simpler version of the method in which only the O(1) terms in the far-field expansion are included (blue).

5. Scattering cross sections

In order to determine the effect of the pre-stress on the scattered waves, we plot the non-dimensionalized scattering cross section (see appendix D),
q=2kA∑n=0∞εn|an|2.5.1
This quantity gives the ratio of the time-averaged power emitted from a scatterer to the intensity of the incoming wave, scaled on the diameter of the scatterer. In the absence of a scatterer, this must obviously be equal to zero.

(a) The effect of the applied pressure difference

In figure 4, we plot the scattering cross section of a neo-Hookean material (black), a Mooney–Rivlin material with C1=0.724, C2=0.276 (blue) and an Arruda–Boyce material with N=26.5 (red), for three values of the non-dimensionalized pressure difference: (p∞−pa)/μ=0 (solid), (p∞−pa)/μ=1 (dashed) and (p∞−pa)/μ=−1 (dotted), as a function of KA. The values of the deformed radii a that result from these pre-stress values are shown in table 1. The value of the longitudinal stretch was chosen to be ζ=1 and the initial cavity radius was taken to be A=1. The numerical solver used to evaluate equations (3.11) and (3.13) was NDSolve in Mathematica 9.0 (Wolfram Research, Inc.) and the value selected for the outer radius of the numerical domain was b=80 in all cases.

The scattering cross section of a neo-Hookean (black), a Mooney–Rivlin with C1=0.724, C2=0.276 (blue) and an Arruda–Boyce material with N=26.5 (red) for three values of the non-dimensionalized pressure difference: (p∞−pa)/μ=0 (solid), (p∞−pa)/μ=1 (dashed) and (p∞−pa)/μ=−1 (dotted) as a function of KA. The value of the longitudinal stretch was chosen to be ζ=1 and the initial cavity radius was taken to be A=1.

The effect of the applied pressure differences on the deformed radius a for the three strain energy functions, given a longitudinal stretch ζ=1 and undeformed radius A=1.

We note that, as discussed earlier, all the neo-Hookean scattering cross sections and the stress-free scattering cross sections of the other materials coincide. In the Mooney–Rivlin material, a greater pressure at infinity than on r=a leads to a decrease in scattered energy, whereas a negative pressure difference increases the scattering. This is intuitive since a positive pressure difference decreases the size of the scatterer and a negative pressure difference increases it. Interestingly, however, in the Arruda–Boyce material any non-zero pressure difference (positive or negative) causes a small decrease in scattering (counterintuitively, the effect is more pronounced for a negative pressure difference). This can be seen in figure 5, which shows an enlarged version of the region 1.5≤KA≤2, 1.2≤q≤1.36 in figure 4. As the Arruda–Boyce strain energy function is essentially a higher order correction to a neo-Hookean strain energy function, involving the invariant I1 only, it is perhaps not surprising that the deviation of the behaviour of this material from neo-Hookean is small; however, the reversal of the effect of a negative pressure difference as compared with a Mooney–Rivlin material is somewhat surprising.

An enlarged version of the region 1.5≤KA≤2, 1.2≤q≤1.36 displayed in figure 4.

(b) The effect of the applied longitudinal stretch

In figure 6, we plot the scattering cross section of a neo-Hookean material (black), a Mooney–Rivlin material with C1=0.724, C2=0.276 (blue) and an Arruda–Boyce material with N=26.5 (red), for three values of the longitudinal stretch: ζ=1 (solid), ζ=12 (dashed) and ζ=2 (dotted), as a function of KA. The values of the deformed radii a that result from these pre-stress values are tabulated in table 2. The value of the non-dimensionalized pressure difference was chosen to be (p∞−pa)/μ=1 and the initial cavity radius was chosen to be A=1. Once again, the numerical solver used to solve equations (3.11) and (3.13) was NDSolve in Mathematica 9.0 (Wolfram Research, Inc.) and the value selected for the outer radius of the numerical domain was b=80.

The scattering cross section of a neo-Hookean (black), a Mooney–Rivlin with C1=0.724, C2=0.276 (blue) and an Arruda–Boyce material with N=26.5 (red) for three values of the longitudinal stretch: ζ=1 (solid), ζ=12 (dashed) and ζ=2 (dotted), as a function of KA. The value of the non-dimensionalized pressure difference was chosen to be (p∞−pa)/μ=1 and the initial cavity radius was chosen to be A=1.

The effect of the applied longitudinal stretches on the deformed radius a for each strain energy function used, given a non-dimensionalized pressure difference (p∞−pa)/μ=1 and undeformed radius A=1.

In all cases, a longitudinal stretch leads to a decrease in the scattered power, whereas a longitudinal compression increases it. This effect could have been predicted by observing that the scattering cross section is scaled on the deformed wavenumber k, which increases with increasing stretch and decreases with increasing compression (see equations (3.10), (3.12) and (3.14)). Once again, we observe that the Arruda–Boyce scattering cross sections take values very close to those found for the neo-Hookean strain energy function.

6. Discussion

In this paper, we have investigated the effect of pre-stress on the scattering of antiplane elastic waves from a cylindrical cavity. The governing equations have inhomogeneous coefficients, and, as a result, could not be solved analytically except for special cases. In order to analyse their behaviour, therefore, a new hybrid analytical–numerical method was developed which relies upon the fact that the inhomogeneities in the coefficients approach zero as r→∞, thus allowing an asymptotic solution to be derived for large r. The asymptotic solution was used to determine a boundary condition for a numerical solver in the region around the cavity. This hybrid method could be applied to any linear ordinary differential equation whose inhomogeneous coefficients tend to a constant as its independent variable tends to infinity.

In order to analyse the effect of the pre-stress, the scattering cross section was plotted for several values of the longitudinal stretch and pressure difference. It was shown that a positive longitudinal stretch led to a decrease in scattered power, whereas a compression caused an increase. The effect of the pressure difference, however, was more interesting. It was shown that, for a Mooney–Rivlin material, a positive pressure difference led to a decrease in scattered power, whereas a negative pressure difference led to an increase. This is intuitive as one would expect that inflating a cavity would increase its ability to scatter waves, whereas a deflation should decrease that ability. In the neo-Hookean case, however, the scattered power was completely unaffected by the pre-stress, whereas, in the Arruda–Boyce case, any non-zero pressure difference (positive or negative) led to a decrease in scattering. This result was not expected; however, it can be explained by considering the form of the strain energy function that was used. The two-term Arruda–Boyce model introduces an O(M2) modification to the neo-Hookean incremental equation (this can be seen by comparing equation (3.13) with equation (3.9)). Since M2 is positive for any real M, these new terms are insensitive to whether M is positive (corresponding to a cavity compression) or negative (corresponding to a cavity inflation) and therefore always reduce the scattering. It is possible that, by expanding the Arruda–Boyce model to third order in I1, the introduction of O(M3) terms in the governing incremental equation would reverse this behaviour; however, it is likely that these terms would be so small (due to the O(1/N2) coefficients that multiply the O(M3) terms—see equation (2.17) for reference) that the O(M2) terms would still dominate. In either case, the difference between the Arruda–Boyce and the neo-Hookean material response is expected to be small.

We note that the neo-Hookean and Arruda–Boyce strain energy functions are independent of the second strain invariant I2, and can therefore be classified as ‘I1 models’. It is well known that I1 models display unphysical behaviour in many circumstances. For example, they do not exhibit the Poynting effect [37]. It appears that the problem raised here for the Arruda–Boyce model is another example of the deficiencies of such strain energy functions.

Data accessibility

This work does not have any experimental data. Calculations were performed using Mathematica v. 9.0.

Author' contributions

T.S. applied the method described in §4 to obtain the results, analysed the results, participated in the design of the study and drafted the main body of the manuscript; W.J.P. conceived of and participated in the design of the study and drafted the Introduction section; I.D.A. conceived of and participated in the design of the study and developed the method described in §4.

Competing interests

We have no competing interests.

Funding

T.S. and W.J.P. would like to thank EPSRC for supporting this work via their Fellowship grant nos. EP/L017997/1 and EP/L018039/1, respectively. I.D.A. undertook part of this work while in receipt of a Royal Society Wolfson Research Merit Award.

Appendix A. Summary of small-on-large theory

For the small-on-large details, we follow the approach described in Ogden [23] and Destrade & Saccomandi [38], with slight modifications to notation. We consider incremental deformations from the deformed body b. To this end, consider a finite deformation of the original body B to a new deformed state b¯ which is close to the configuration b. We shall call this new configuration the perturbed configuration. Position vectors in the deformed states b and b¯ are defined by x and x¯, respectively, and we define
u=x¯−x,|u|≪1A 1
as the difference between position vectors in b¯ and b. Since b¯ is close to b, u is called an incremental displacement which, in our case, is time dependent. We assume that the time-harmonic incremental deformation is an antiplane wave, i.e.
u=w(x,y)e−iωtez,A 2
where ez is a unit vector in the z-direction in the body b.

Given the deformation gradient (2.4), we can define the additional deformation gradients f=gradx¯ and F¯=Gradx¯=fF, where grad represents the gradient operator with respect to the statically deformed configuration. Therefore,
Γ=Gradu=Grad(x¯−x)=(f−I)FA 3
and similarly γ=grad u=f−I, so that f=I+γ and Γ=γF.

Furthermore, given the stresses (2.9) in the configuration b, we can define the corresponding stresses T¯=T+τ and S¯=S+s in b¯, where τ and s are the incremental Cauchy and nominal stresses, respectively. By taking an expansion about the deformation F, it is straightforward to show that for an incompressible material
s=L:γF+qF−1−QF−1γ,A 4
where q is the increment of Q and L=∂2W/∂F2 defined componentwise via LIjKl=∂2W/∂FjI∂FlK. The colon notation indicates the double contraction A:b=Aijklblk. Defining the push-forward of the incremental nominal stress as ζ=Fs, from (A 4) this therefore has the form
ζ=M:γ+qI−Qγ,A 5
where the components of the instantaneous modulus tensor M are defined by
Mijkl=∂2W∂FjM∂FlNFiMFkN.A 6
Convenient forms for this tensor have been given by Ogden [23] as
2Miijj=λiλjWij=Mjjii,no sum oni,j,A 7Mijij=(λiWi−λjWjλi2−λj2)λi2,i≠j,λi≠λj,A 8Mijij=12(Miiii−Miijj+λiWi),i≠j,λi=λjA 9andMijji=Mjiij=Mijij−λiWi,i≠j,A 10
where Wi=∂W/∂λi and Wij=∂2W/∂λi∂λj. Analysis of the Cauchy stress gives rise to the relationship
τ=ζ+γT,A 11
for an incompressible material.

The equations of motion in body b¯ in terms of the Cauchy stress are
div¯T¯=ρ¯∂2U¯∂t2=ρ∂2u∂t2,A 12
where div¯ represents the divergence operator in the perturbed configuration. The second of the above equalities holds as the body is incompressible (i.e. ρ¯=ρ) and because U¯=U+u where U is not time dependent. Using the fact that div T=0, it is relatively straightforward to show that
div¯T¯=div(τ−γT)=divζ,A 13
where we have used the relationship between ζ and τ above. The incremental equations are therefore
divζ=ρ∂2u∂t2.A 14
If we substitute (A 2) into the above, two of the equations imply that q is independent of r and θ, and the single remaining governing equation of motion is
1r∂∂r(rζrz)+1r∂ζθz∂θ+ρω2w+∂q∂z=0,A 15
where
ζrz=Mrzrzγzrandζθz=Mθzθzγzθ.A 16
In this case τrz=ζrz and τθz=ζθz and so the instantaneous moduli Mrzrz and Mθzθz can be interpreted as pressure-dependent anisotropic shear moduli, given by
Mrzrz(r)=μr(r)=(λrWr−λzWzλr2−λz2)λr2andMθzθz(r)=μθ(r)=(λθWθ−λzWzλθ2−λz2)λθ2.A 17
Using (A 16) with (A 17) in (A 15) gives the governing equation of motion in terms of the displacement w, which is
1r∂∂r(rμr(r)∂w∂r)+μθ(r)r2∂2w∂θ2+ρω2w+∂q∂z=0.A 18
We now apply the condition that the pressure on the surface of the cavity (i.e. where r=a) remains unaltered by the perturbation, so that
T¯n¯=−pn¯,A 19
where n¯ is the outer unit normal in the perturbed configuration. This condition reduces to
∂w∂r|r=a=0,q=0,A 20
and the second of these equations allows us to reduce (A 18) to
1r∂∂r(rμr(r)∂w∂r)+μθ(r)r2∂2w∂θ2+ρω2w=0.A 21

Appendix B. Integrals

Here we describe how to numerically evaluate α, β, γ and δ defined in equations (4.14)–(4.18). We begin by noting that γ has two integrals of the form
Ip=∫1∞(Hn(1)(κx))2xpdx,p=1,3.B 1
The integrand here is highly oscillatory and so is not easy to evaluate numerically. We note, however, that the singularities of Hn(1)(κx) lie entirely in the lower half of the complex x-plane. Let us write u=κx, so that
Ip=κp−1∫κ∞(Hn(1)(u))2updu.B 2
Then, we note that, in the limit as |u|→∞,
(Hn(1)(u))2∼−2iπu(−1)ne2iu,B 3
therefore we can rotate the contour in Ip into the upper half-plane:
Ip=κp−1∫κi∞(Hn(1)(u))2updu.B 4
This form can now be easily evaluated numerically. Similarly, we note that α has integrals of the form
Lp=∫1∞1y(Hn(1)(κy))2∫y∞(Hn(1)(κx))2xpdxdy,p=1,3.B 5
Upon making the substitutions u=κx, v=κy and again rotating into the upper half-plane, we obtain
Lp=κp−1∫κ∞1v(Hn(1)(v))2∫v∞(Hn(1)(u))2updudv=κp−1∫κ∞1v(Hn(1)(v))2∫vi∞(Hn(1)(u))2updudv,B 6
which, again, is readily computable. Furthermore, δ contains integrals that can be written as
Jp=∫1∞Hn(κx)Jn(κx)xpdx=κp−1∫κ∞Hn(1)(u)Jn(u)updu=κp−12∫κ∞(Hn(1)(u))2+Hn(1)(u)Hn(2)(u)updu=Ip2+κp−12∫κ∞Hn(1)(u)Hn(2)(u)updu,p=1,3,B 7
where Hn(2)(u) is the Hankel function of the second kind. The last of these integrals is not oscillatory, and so can be numerically evaluated as is. The other integral in δ can be written as
J=∫1∞1x2Hn(1)(κx)ddx(Jn(κx))dx=κ2∫1∞1x2Hn(1)(κx)(Jn−1(κx)−Jn+1(κx))dx=κ24∫κi∞1u2Hn(1)(u)(Hn−1(1)(u)−Hn+1(1)(u))du+κ24∫κ∞1u2Hn(1)(u)(Hn−1(2)(u)−Hn+1(2)(u))du.B 8
Finally, the double integrals in β can be written as
Mp=∫1∞1y(Hn(1)(κy))2∫y∞Hn(1)(κx)Jn(κx)xpdxdy=κp−1∫κ∞1v(Hn(1)(v))2∫v∞Hn(1)(u)Jn(u)updudv=κp−12∫κ∞1v(Hn(1)(v))2∫vi∞(Hn(1)(u))2updudv+κp−12∫κ∞1v(Hn(1)(v))2∫v∞Hn(1)(u)Hn(2)(u)updudv=Lp2+κp−12∫κ∞1v(Hn(1)(v))2∫v∞Hn(1)(u)Hn(2)(u)updudv,p=1,3,B 9
and
P=∫1∞1y(Hn(1)(κy))2∫y∞1x2Hn(1)(κx)ddx(Jn(κx))dxdy=κ2∫κ∞1v(Hn(1)(v))2∫v∞1u2Hn(1)(u)ddu(Jn(u))dudv=κ24∫κ∞1v(Hn(1)(v))2(∫vi∞Hn(1)(u)(Hn−1(1)(u)−Hn+1(1)(u))u2du+∫v∞Hn(1)(u)(Hn−1(2)(u)−Hn+1(2)(u))u2du)dv,B 10
which completes the rewriting of all the integrals in α–δ to make them easy to evaluate.

Appendix C. Convergence of the method described in §4

The asymptotic expansion (4.2) includes terms up to O(ϵ2), and so one might expect the error in determining each scattering coefficient to be of the order of the discarded terms, i.e. O(ϵ4). However, the asymptotic expansion is not the only source of error in the method. With reference to equation (4.19), we note that an can be represented schematically as
an=(An(0)+ϵ2An(2))fn(b)+(Bn(0)+ϵ2Bn(2))fn′(b)+O(ϵ4),C 1
where An(0), An(2), Bn(0) and Bn(2) can be derived by series expanding δ/Δ and −β/Δ.

Because fn(b) is derived computationally, its value will include some level of numerical error en(1)(b), say, so that we can write
fn(b)=f^n(b)+en(1)(b),C 2
where f^n(b) is the exact value of fn(b). The function fn(b) is derived using a shooting method, so the error en(1)(b) will remain small until, at some very large value of b, it will start to grow. Upon substituting (C 2) into (C 1), we obtain
an=(An(0)+ϵ2An(2))(f^n(b)+en(1)(b))+(Bn(0)+ϵ2Bn(2))(f^n′(b)+en(2)(b))+O(ϵ4),C 3
where en(2)(b) is the numerical error in determining the value of fn′(b), which we assume to be of the same size as en(1)(b).

Now let us series expand the true value of the nth-order scattering coefficient a^n:
a^n=(A^n(0)+ϵ2A^n(2))f^n(b)+(B^n(0)+ϵ2B^n(2))f^n′(b)+O(ϵ4),C 4
where A^n(0), A^n(2), B^n(0) and B^n(2) are the true values of An(0), An(2), Bn(0) and Bn(2), respectively. The terms An(0) and Bn(0) in (C 3) should be exact as their evaluation does not require the numerical solution of any integrals. The terms An(2) and Bn(2) will, however, involve some level of numerical error (dn(1) and dn(2), say, which we assume are of the same order of magnitude as each other), and since the evaluation of the relevant integrals does not depend on a shooting method, we do not expect these to change with increasing b. Hence, for our purposes we have
An(0)=A^n(0),An(2)=A^n(2)+dn(1),Bn(0)=B^n(0)andBn(2)=B^n(2)+dn(2).C 5
Therefore, if we subtract (C 3) from equation (C 4), we obtain the total error that arises from the method described in §4:
En=|anexact−an|=|An(0)en(1)(b)+ϵ2en(1)(b)An(2)+ϵ2dn(1)f^n(b)+ϵ2en(1)(b)dn(1)+Bn(0)e^n(2)(b)+ϵ2e^n(2)(b)Bn(2)+ϵ2dn(2)f^n′(b)+ϵ2en(2)(b)dn(2)+O(ϵ4)|.C 6
For small values of b, the O(en(j)(b)) and O(ϵ2dn(j)), j=1,2, terms in equation (C 6) will be negligible and so the O(ϵ4) error terms will dominate. Thus, the method will converge like 1/b4; however, as b increases the O(ϵ4) terms will rapidly decrease in size until they are smaller than the O(ϵ2dn(j)) terms. Therefore, the solution will then converge at a rate proportional to 1/b2. For larger values of b, these terms will also decay to the point where they are smaller than O(en(j)(b)) quantities. The latter will offer a very small constant error to the solution until, at a very large value of b, the shooting method will fail and the error will start to grow.

The key to selecting b is to choose a value that is large enough so that the O(ϵ4) and O(ϵ2dn(j)) errors are insignificant, but small enough that the O(en(j)(b)) error is also insignificant. The magnitude of the O(en(j)(b)) terms can be decreased (and therefore the error of the overall scheme decreased) by increasing the accuracy of the numerical solver; however, this will, of course, increase computation times.

Appendix D. Scattering cross section

The scattered power flow through a surface S surrounding the cavity may be written [39] as
Pout=∫SPout⋅dA,D 1
where dA=n dA is an incremental area element, n is an outer unit normal to the surface S and
Pout=−∂uout∂t⋅τoutD 2
is a vector that is analogous to the Poynting vector for electromagnetic waves, where uout represents the scattered part of u, and τout is the part of the incremental Cauchy stress (defined in (A 11)) associated with the scattered waves. We shall choose S to be a cylinder of unit axial length that is concentric with the cavity and with radius r*≫a, so that
n=eranddA=r∗dθdz.D 3
Substituting (D 2) and (D 3) into (D 1), we obtain
Pout=∫01∫02π(−∂wout∂tτzroutr∗)dθdz=∫02π(−∂wout∂tτzroutr∗)dθ,D 4
where wout is defined via uout=woutez, and τzrout is the (z,r)th component of τout. Because, in general, we do not have explicit expressions for u and τ, we shall use their far-field forms, so that
wout∼Re(∑n=0∞ϵninanHn(1)(kr∗)cos⁡(nθ)e−iωt)andτout=μr∗∂wout∂r,D 5
where
μr∗=limr→∞μr(r).D 6
The time-averaged scattered power is defined as follows:
P¯out=ω2π∫02π/ωPoutdtD 7=ω2π∫02π/ω∫02π(−∂wout∂tτzroutr∗)dθdtD 8=−ωμr∗r∗2π∫02π/ω∫02π(∂wout∂t∂wout∂r)dθdt.D 9
Upon substituting (D 5) into (D 6), and letting r∗→∞, it can be shown that, after substantial algebraic simplification,
P¯out=2μr∗ω∑n=0∞ϵn|an|2.D 10
The intensity of the incoming wave is defined as the time-averaged power passing through a unit square, which using the same methods as above can be obtained as follows:
I=P¯in=μr∗kω2.D 11
The scattering cross section, then, is defined as follows:
Q=P¯outI=4k∑n=0∞ϵn|an|2.D 12

As in Lewis et al. [40], however, we prefer the non-dimensionalized form of the scattering cross section:
q=Q2A=2kA∑n=0∞ϵn|an|2.D 13
Note that we have non-dimensionalized on the undeformed cavity radius A.

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.

. 1926VI. On the concentration of stress in the neighbourhood of a small spherical flaw; and on the propagation of fatigue fractures in ‘statistically isotropic’ materials. Lond. Edin. Dublin Phil. Mag. J. Science1, 71–97. (doi:10.1080/14786442608633614)