Antiplane wave scattering from a cylindrical cavity in pre-stressed nonlinear elastic media

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.

TS, 0000-0001-7536-5547 The effect of a longitudinal stretch and a pressureinduced 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 prestress 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.

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][4][5][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][10][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][17][18].
Recently, nonlinear elastic materials have been used in order to construct tunable band gaps that can be shifted in real time [19][20][21][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.

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(I 1 , I 2 ) [23], where I j , j = 1, 2, 3 are the principal strain invariants, and we note that W is independent of I 3 due to the constraint of incompressibility (i.e. I 3 = 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 where (R, Θ, Z) and (r, θ, z) are cylindrical polar coordinates (X = R cos Θ, Y = R sin Θ, x = r cos θ, y = r sin θ ) 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, where E R and E Z are the radial and longitudinal basis vectors in the undeformed configuration, and e r and e z are the corresponding basis vectors in the deformed configuration. Using (2.1), it can be shown that the principal stretches for this deformation in the radial, azimuthal and longitudinal directions, respectively, are , λ θ = r R(r) and λ z = ζ . 3) The deformation gradient tensor F is given by  The deformed cavity, with radius a. The deformation is due to an imposed far-field hydrostatic pressure p ∞ and an internal pressure p a , along with an axial stretch ζ .
where Grad represents the gradient operator in the undeformed configuration. In our case, we have For an incompressible material J = √ I 3 = det F = 1 and so This is easily solved to yield where M is a constant defined by 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 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: The static equations of equilibrium are then given by div T = 0, (2.11) where div signifies the divergence operator with respect to x. These reduce to where T = T ij e i ⊗ e j . The second and third of these imply Q = Q(r). Integrating the first of (2.12), using (2.9) and applying T rr | r=a = −p a , we find that the radial stress is given by then by applying T rr → −p ∞ as r → ∞, we obtain an equation relating the pressure difference to the resulting deformation, (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: where μ is the ground-state shear modulus of the material under consideration, the Mooney-Rivlin: where C 1 and C 2 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 I 1 [33]. For simplicity, we shall use the first two terms only: (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 in the Mooney-Rivlin case is and in the Arruda-Boyce case is  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.

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 stateb which is close to the configuration b. The position vector in the new deformed stateb is given byx and we define as the difference between position vectors inb and b. Sinceb 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.
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, 1 r Given an incident field that takes plane wave form in the far-field, i.e. we seek solutions to (3.3) in the form ε n i n f n (r) cos(nθ). (3.7) As such for the nth mode, we have 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 in the Mooney-Rivlin case we have and in the Arruda-Boyce case we have where and we note that, when ζ = 1, k NH = k MR = k AB = K. By substituting (3.7) into (A 20), we see that the boundary condition on r = a reduces to f n (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/r 2 → 0, (3.9), (3.11) and (3.13) become where k 2 = k 2 NH in the neo-Hookean case, k 2 = k 2 MR in the Mooney-Rivlin case, and k 2 = k 2 AB in the Arruda-Boyce case. Equation stress-free elastic medium (albeit with a modified wavenumber) and is readily solved: , f i n (r) = J n (kr) and f s n (r) = a n H (1) n (kr), (3.17) where we identify f i n (r) as the nth component of the incoming plane wave and thus equate it to J n via (3.5). Furthermore, f s n (r) is the nth component of the scattered wave, represented by a Hankel function of the first kind, H (1) n (kr), and a n 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 f n (r) = J n K ζ (r 2 + M) + a n H (1) n K ζ (r 2 + M) (3.18) for the case of an incoming plane wave, where the scattering coefficients a n are given by . (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 which can be solved analytically (as for the neo-Hookean case) to give where c 1 and c 2 are arbitrary constants. Upon matching this solution to the far-field form (3.17), we obtain c 1 = 1 and c 2 = a 0 , and applying the boundary condition (3.15), we obtain This is a similar form to that which the neo-Hookean scattering coefficients take when n = 0; however, since k a 2 + 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.

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/r 2 → 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/b 2 = 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 F n (s) = f n (r), respectively, so that (3.11) can be rewritten as where κ 2 = (kb) 2 . Given (3.17), we then assume the following (regular) expansion for F n (s): F n (s) = J n (κs) + a n H (1) 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 F n (s) = J n (κs) + a n H (1) n (κs); hence G n and g n 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(H n (κx)J n (κx) dx dy − a n 2s 2 + a n and, therefore, f n (r) = J n (kr) + a n H (1) and, hence, Now, in order to determine the appropriate boundary condition on r = b, we require a suitable choice for a n . We note that we can write where n (κx)) 2 dx dy , (4.14) n (κx)J n (κx) dx dy (4.15) and f n (b) = df n dr r=b = γ a n + δ,  Solving (4.13) and (4.16) for a n , we obtain and eliminating a n from (4.13) and (4.16), we obtain our final boundary condition, The method is now clear. We choose b such that M/b 2 1, and numerically solve equation (3.11) subject to (3.15) and (4.21). Once solved, we use the values of f n (b) and f n (b) determined by our numerical solver in equation (4.19) to deduce the value of a n . The predicted value of a n 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: d ds (J n (κs) + a n H (1) n (κ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 a 0 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 f n (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.

Scattering cross sections
In order to determine the effect of the pre-stress on the scattered waves, we plot the nondimensionalized scattering cross section (see appendix D), 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. We note that, as discussed earlier, all the neo-Hookean scattering cross sections and the stressfree 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 I 1 only, it is perhaps not surprising that the deviation of the behaviour of this material from neo-Hookean   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.

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 twoterm Arruda-Boyce model introduces an O(M 2 ) modification to the neo-Hookean incremental equation (this can be seen by comparing equation (3.13) with equation (3.9)). Since M 2 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 I 1 , the introduction of O(M 3 ) 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/N 2 ) coefficients that multiply the O(M 3 ) terms-see equation (2.17) for reference) that the O(M 2 ) 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 I 2 , and can therefore be classified as 'I 1 models'. It is well known that I 1 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. Furthermore, given the stresses (2.9) in the configuration b, we can define the corresponding stressesT = T + τ andS = S + s inb, 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 = ∂ 2 W/∂F 2 defined componentwise via L IjKl = ∂ 2 W/∂F jI ∂F lK . The colon notation indicates the double contraction A : b = A ijkl b lk . Defining the push-forward of the incremental nominal stress as ζ = Fs, from (A 4) this therefore has the form where the components of the instantaneous modulus tensor M are defined by n =Â (2) n + d (1) n , B n =B (0) n and B (2) n =B (2) n + d (2) n .
Therefore, if we subtract (C 3) from equation (C 4), we obtain the total error that arises from the method described in §4: E n = |a exact n − a n | = |A n ), 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/b 4 ; however, as b increases the O( 4 ) terms will rapidly decrease in size until they are smaller than the O( 2 d (j) n ) terms. Therefore, the solution will then converge at a rate proportional to 1/b 2 . For larger values of b, these terms will also decay to the point where they are smaller than O(e (j) n (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 n (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.