## Abstract

The volume oscillation of a cylindrical bubble in a microfluidic channel with planar elastic walls is studied. Analytical solutions are found for the bulk scattered wave propagating in the fluid gap and the surface waves of Lamb-type propagating at the fluid–solid interfaces. This type of surface wave has not yet been described theoretically. A dispersion equation for the Lamb-type waves is derived, which allows one to evaluate the wave speed for different values of the channel height *h*. It is shown that for *h*<λ_{t}, where λ_{t} is the wavelength of the transverse wave in the walls, the speed of the Lamb-type waves decreases with decreasing *h*, while for *h* on the order of or greater than λ_{t}, their speed tends to the Scholte wave speed. The solutions for the wave fields in the elastic walls and in the fluid are derived using the Hankel transforms. Numerical simulations are carried out to study the effect of the surface waves on the dynamics of a bubble confined between two elastic walls. It is shown that its resonance frequency can be up to 50% higher than the resonance frequency of a similar bubble confined between two rigid walls.

## 1. Introduction

This study is inspired by the use of bubbles in acoustically activated microfluidic systems [1]. An example of such a system is described by Rabaud *et al*. [2]: micrometre-size cylindrical bubbles are positioned in a fluid channel, 25 μm in depth, confined by elastic walls which are set into oscillation by a vibrating glass rod moulded in the upper wall. Under such conditions, the bubble dynamics exhibit a number of interesting effects, such as the generation of surface waves propagating along the solid–fluid interfaces and the self-organization of bubbles in crystal-like structures where equilibrium distances between bubbles are much smaller than distances that are predicted by the theory of secondary Bjerknes forces [3]. Another phenomenon of interest is acoustic microstreaming generated by bubbles in microfluidic devices [4]. This phenomenon is used for micromixing of fluids [5] and sorting of microparticles and cells [6,7]. A variety of experimental data on the dynamics of bubbles in narrow planar channels can be found in a series of recent papers [8–13], where both single bubbles and pairs of interacting bubbles are investigated.

Available theoretical models on the dynamics of cylindrical bubbles in planar channels consider only bulk scattered waves generated by the bubble oscillation in the fluid gap [14–16]. These models are appropriate if the walls of the channel are rigid. There are no models that allow for surface waves, which are generated at the solid–fluid interfaces if the channel walls are elastic, whereas experimental observations show that surface waves can play a decisive role in the motion of the fluid layer [2].

In acoustics, several types of surface waves are encountered [17]. The surface waves relevant to our study are Rayleigh, Scholte and Lamb waves, so it is pertinent to remind their definitions. Rayleigh waves are surface waves that travel at an interface between an elastic solid and a vacuum. The existence of these waves was theoretically predicted by Lord Rayleigh in 1885 [18]. The speed of Rayleigh waves is slightly less than that of shear waves propagating in the body of an elastic solid. Scholte waves are surface waves that propagate at an interface between an elastic solid and a fluid. They are named after Scholte, who discovered them in 1947 [19]. Scholte waves are slower than Rayleigh waves. Lamb waves propagate in an elastic plate placed in a vacuum or a fluid [20]. The mathematical analysis of these waves was first developed by Lamb in 1917 [21]. Lamb waves are divided into symmetric and antisymmetric types depending on whether the motion of the elastic medium is symmetric or antisymmetric about the middle plane of the elastic plate.

Although the type of surface waves considered in our study is intimately related to the types described above, it is different from them. In a certain sense, the case under study is opposite to the case considered by Lamb because we deal with surface waves propagating in a fluid channel embedded in an elastic solid. To the best of our knowledge, this case has not been considered as yet. Our study provides a dispersion equation for this type of surface waves, which allows one to evaluate the speed of the surface waves for different values of the channel height. The ultimate purpose of our study is to develop a theoretical model for the volume oscillation of a cylindrical bubble placed between two planar elastic walls, allowing for both bulk and surface waves generated in this system. It should be emphasized that we deal with a special case of bubble confinement. Indeed, the channel height is so small with respect to the bubble size that the bubble always remains in contact with the walls, which leads to a shape close to a cylinder, as depicted in figure 1. Therefore, there is no possibility for the bubble to take the spherical shape or be at a distance from the walls. This fact is a result of the design of microfluidic devices of our interest [2], which makes our case different from those implied by Ilinskii *et al*. [15] and Leighton *et al*. [22–24].

## 2. Problem formulation and solution

Let us consider a cylindrical gas bubble located in a fluid channel of height *h* between two planar elastic walls (figure 1). The bubble undergoes volume oscillation in response to an imposed acoustic pressure field. The bubble oscillation is considered to be linear and dependent on time as exp(−i*ωt*). For simplicity, the time factor will be omitted from equations. Quantities related to the upper and lower walls will be denoted by subscripts 1 and 2, respectively. The action of the bubble on the fluid–solid interfaces is modelled as a normal harmonic point load. Because the problem under consideration is axially symmetric, cylindrical coordinates are used, whose origin is located in the middle of the bubble and the *z*-axis is perpendicular to the surface of the walls, as shown in figure 1. It should be emphasized that the wave field in the fluid channel consists of two components. One component is the bulk scattered wave caused by the bubble radial oscillation. The second component is caused by the surface waves that propagate in the elastic walls and which, in the process of this propagation, induce perturbations in the fluid. This component will be called surface waves in the fluid channel.

The calculation consists of the following stages. The equations for the surface waves in the elastic walls and in the fluid channel are derived in §2a and §2b, respectively. The derivation is based on Hankel transforms. The boundary conditions on the solid–fluid interfaces are obtained in §2c. The dispersion equation for the surface wave speed is considered in §2d. The resulting expressions for the surface waves in the space domain are calculated by the inverse Hankel transforms in §2e. The equations for the bulk scattered wave in the fluid channel are given in §2f. The boundary conditions on the side bubble surface are obtained in §2g. The equation of bubble oscillation is derived in §2h.

### (a) Waves in the elastic walls

As shown in figure 1, the walls occupy the space with *j*th wall (*j*=1,2) can be written as
*φ*_{j} and *ψ*_{j} being called the scalar and vector potentials, respectively. Equation (2.1) is used as the most general mathematical representation for an arbitrary vector field (e.g. [17] or [25]). In the view of axial symmetry, the dependence on the angular coordinate *θ* is absent and hence *φ*_{j} and *ψ*_{j} can be written as *φ*_{j}=*φ*_{j}(*r*,*z*) and *ψ*_{j}=*ψ*_{j}(*r*,*z*)*e*_{θ}, respectively, where *e*_{θ} is the unit azimuth vector. The function *φ*_{j} and *ψ*_{j} obey the following equations [25]:
*k*_{l}=*ω*/*c*_{l} and *k*_{t}=*ω*/*c*_{t} are the wavenumbers of the longitudinal and transverse waves in the solid, respectively. The speeds of these waves, *c*_{l} and *c*_{t}, are given by [17]
*μ* are the Lamé coefficients of the walls (*μ* is also called shear modulus) and *ρ*_{s} is the wall density. The components of the displacement vector *u*_{j} are calculated by [25]

Solutions to equations (2.2) and (2.3) can be found by applying Hankel transforms with respect to the radial variable *r*. The Hankel transform of order *n* of a function *f*(*r*) is defined as [26]
*k* is the parameter of the Hankel transform and *J*_{n} is the Bessel function of the first kind of order *n*. The inverse transform is given by
*z*| are
*Φ*_{j}(*k*) and *Ψ*_{j}(*k*) are functions to be determined from the boundary conditions on the solid–fluid interfaces at *z* = ±*h*/2.

Application of the Hankel transform to equations (2.5) yields

### (b) Surface waves in the fluid channel

The fluid area corresponds to −*h*/2≤*z*≤*h*/2. The displacement vector in the fluid can be written as
*φ*_{f} obeys the equation [27]
*k*_{f} = *ω*/*c*_{f} and *c*_{f} is the speed of sound in the fluid. The fluid velocity corresponding to *u*_{f} is *v*_{f}=−i*ω**u*_{f}. The components of *u*_{f} are calculated by
*Φ*_{f1}(*k*) and *Φ*_{f2}(*k*) are functions to be determined from the boundary conditions at *z* = ±*h*/2.

Application of the Hankel transform to equations (2.22) yields

### (c) Boundary conditions on the solid–fluid interfaces

The boundary conditions on the surfaces of the elastic walls are written as follows [25,28]:
*T*_{jrz} and *T*_{jzz} are the shear and normal stresses in the walls, *p*_{f} is the fluid pressure corresponding to *u*_{f}, *Q* is the magnitude of the point load produced by the bubble at the points with *r*=0 and *z*=±*h*/2 and *δ*(*r*) is the Dirac delta function. Equations (2.29) follow from the continuity of the normal displacement at the solid–fluid interfaces, equations (2.30) mean that the shear stress vanishes at the interfaces and equations (2.31) describe the normal stress balance at the interfaces.

The expressions for *T*_{jrz}, *T*_{jzz} and *p*_{f} are given by [25,27]
*ρ*_{f} is the fluid density. Applying the Hankel transforms to equations (2.29)–(2.34), one obtains
*Φ*_{j}(*k*), *Ψ*_{j}(*k*), *Φ*_{f1}(*k*) and *Φ*_{f2}(*k*):
*z*=0. Therefore, the following conditions should be met [17]:

### (d) Dispersion equation for surface waves

The theoretical analysis of Rayleigh waves produced by a point load, based on the use of Hankel integral transforms, is given in a book by Achenbach [25]. An analogous consideration of Scholte waves is presented by Zhu *et al*. [28]. These works show that at large distances from the loading point, where general solutions similar to equations (2.53) and (2.54) turn into far-field expressions giving the Rayleigh and Scholte waves, the dispersion equations for these waves are identical to the denominators of the general solutions similar to equations (2.53) and (2.54). This follows from the fact that a vanishing denominator gives rise to a dominant term in the spectrum of emitted waves. Analogously, in our case the dispersion equation for the surface waves of our interest is given by
*ρ*_{f}→0, equation (2.57) reduces to the equation for Rayleigh waves [17,25], while for infinite channel height, *h*. Examples of such evaluations are presented in §3a. It is noteworthy that equation (2.57) is derived for the first time so it can be considered as one of the main results of this study.

### (e) Expressions for *u*_{1}, *u*_{f} and *p*_{f} in the space domain

*u*

*u*

Applying the inverse Hankel transform to equations (2.16), (2.17), (2.27), (2.28) and (2.40) and using equations (2.53)–(2.55), one obtains
*c* stands for the speed that corresponds to the wavenumber *k*. The normalization by *k*_{t} provides the simplest form of the equations. The method of numerical calculation of integrals (2.63)–(2.67) is described in appendix A.

### (f) Bulk wave in the fluid channel

According to Ilinskii *et al*. [15], the bulk scattered wave generated by the bubble oscillation in the fluid channel can be represented by
*v*_{b} is the fluid velocity and the velocity potential *φ*_{b} is given by the equation
*A* is a constant to be determined from the boundary conditions on the bubble surface. Equations (2.71) and (2.72) satisfy the wave equation for the fluid motion in the case of cylindrical symmetry. The fluid pressure corresponding to *v*_{b} is calculated by
_{f}=2*π*/*k*_{f}, is large compared with the wavelength of the surface waves as well as to the dimensions of the microfluidic system under study. This means that the spatial scale of the variations of *p*_{b} is much greater than that of the perturbations produced by the surface waves. For this reason, *p*_{b} does not contribute to the above perturbations and hence does not appear in equation (2.31).

### (g) Boundary conditions on the bubble surface

Let us assume that the time-varying bubble radius can be represented as
*R*_{0} is the equilibrium bubble radius. Considering that the net fluid velocity is the sum of the velocities of the bulk and surface waves, the boundary condition for the normal velocity at the side bubble surface can be written as
*v*_{br} is the radial component of *v*_{b},*v*_{fr} is the radial component of the surface wave velocity, given by *v*_{fr}=−i*ωu*_{fr}, and the right-hand term is the amplitude of the velocity of the bubble oscillation. Assuming that *h* is small compared with the transverse wavelength λ_{t}=2*π*/*k*_{t}, the value of *u*_{fr} can be taken at *z* = 0; see equations (2.60) and (2.65). Substitution of equations (2.60), (2.71) and (2.72) into equation (2.75) yields
*p*_{g} is the gas pressure within the bubble, *P*_{0} is the hydrostatic pressure in the fluid, *p*_{ac} is the driving acoustic pressure and *p*_{st} is the surface tension pressure. Calculating the bubble volume as
*P*_{g0} is the equilibrium gas pressure and *γ* is the ratio of specific heats of the gas. The acoustic pressure *p*_{ac} is specified by
*P*_{a} denotes the acoustic pressure amplitude. For *p*_{st}, one has
*σ*_{f} is the surface tension coefficient for the fluid–gas interface. Substituting equations (2.62), (2.72), (2.73) and (2.79)–(2.81) into equation (2.77) and separating constant and time-varying terms, one obtains
*μ*, one finds the following expression for *a*:

where
*α*_{f}=*k*_{f}*R*_{0} and *α*_{t}=*k*_{t}*R*_{0}.

### (h) Equation for bubble oscillation

To derive an equation for the amplitude of the bubble oscillation *a*, we need an expression for the load magnitude *Q* in terms of *a*. Clearly, *Q* should be proportional to the variation of the gas pressure in the bubble, given by equation (2.79). It is admissible to assume that *Q* is also proportional to the area of the lateral section of the bubble, *Q* can be represented as
*ε* is a proportionality coefficient that allows for the fact that the effective area that describes the action of the bubble on the wall can be different from *ε* can be a complex quantity because a phase shift can exist between the variation of the gas pressure and the reaction of the walls to it. This occurs because the walls may behave as a viscoelastic material with complex elastic and shear moduli, which leads to a phase lag between stress and strain [30].

Substitution of equation (2.87) into equations (2.84) and (2.86) yields
*ε* for different values of the material and acoustic parameters. It should be mentioned that for *ε*=0, equation (2.88) reduces to the result obtained by Ilinskii *et al*. [15] for a cylindrical bubble between two rigid walls.

## 3. Numerical simulations and experimental verification

### (a) Speed of surface waves

The speed of the surface waves propagating at the solid–fluid interfaces is given by solutions of equation (2.57). Numerical evaluations were performed for the following values of the physical parameters: Young’s modulus *E*=1.6 MPa, Poisson’s ratio *σ*=0.499,*ρ*_{s}=970 kg m^{−3},*c*_{f}=1481 m s^{−1} and *ρ*_{f}=998 kg m^{−3}. These values are typical of the elastic walls of a microfluidic channel made of polydimethylsiloxane (PDMS) and filled with water. The transverse wave speed for these parameters is *c*_{t}=23.456 m s^{−3}. This quantity sets the upper limit for the speed of surface waves [17]. The evaluations show that equation (2.57) has one real root which gives a frequency-dependent speed *c*_{s}<*c*_{t}. Figure 2 exemplifies the frequency dependence of *c*_{s} at *h*=25 μm. The dependence of *c*_{s} on *h* is demonstrated by figure 3. Note that the values of *c*_{s} are normalized by *c*_{t}, and those of *h*, by λ_{t}=*c*_{t}/*f*, where *f* is the driving frequency. With this normalization, the plot shown in figure 3 remains the same for all frequencies. Figure 3 reveals that *c*_{s} increases with increasing *h* and tends to the speed of the Scholte wave, *c*_{Sch}, when *h* tends to λ_{t}. However, at small *h*, the presence of the second boundary makes *c*_{s} much smaller than *c*_{Sch}.

### (b) Solutions at small and large distances from the bubble

The theory of surface waves predicts that the behaviour of integrals (2.63)–(2.67) should coincide with the behaviour of the Hankel functions of zero and first orders at large distances from the loading point [17]. The calculation of the integrals by the method described in appendix A conforms to this expectation. As an example, figure 4 compares the behaviour of the integral *I*_{sz}, which describes the vertical displacement of wall 1, and the Hankel function *k*_{s}=*ω*/*c*_{s} is the wavenumber of the surface wave, given by equation (2.57). The calculations were made at *f*=30 kHz and *h*=25 μm, the material parameters being the same as in the preceding subsection. Figure 4*a* shows Re{*I*_{sz}(*r*,*h*/2)} (solid) and *r* which is normalized by λ_{s}=2*π*/*k*_{s}. The phases of these functions are compared in figure 4*b*. Note that the vertical lines in this figure do nothing but show that the sign of the phase is changed by jump. The comparison reveals that the functions give the same result for distances exceeding 0.4λ_{s}. However, for smaller distances, there is a considerable difference in their behaviour. Therefore, the approximation by the Hankel function is not suitable for use in equations (2.88) and (2.89) which describe the bubble oscillation.

### (c) Resonance curves

Comparison with experimental measurements (presented in the next subsection) shows that agreement with the theoretical model is achieved if the parameter *ε* is considered as a complex quantity whose phase is proportional to the bubble radius *R*_{0}. Therefore, examples of resonance curves were calculated assuming that
*ε*_{1} and *ε*_{2} serve as variation parameters.

Resonance curves shown in figure 5 were calculated by equation (2.88) for *h*=25 μm and *f*=30 kHz, varying *R*_{0} in the range from 16 to 83 μm. The material parameters for the solid and the fluid were the same as in the preceding subsections. The parameters related to the gas behaviour and the acoustic excitation were as follows: *σ*_{f}=0.072 N m^{−1}, *γ*=1.4, *P*_{0}=101.3 kPa and *P*_{a}=10 kPa. Figure 5*a* shows resonance curves for increasing values of *ε*_{1}, assuming *ε*_{2}=0. It is seen that at *ε*_{1}=0 (no surface waves), the resonance peak is at *R*_{0}*f*=1.1 m s^{−1}, which corresponds to the result obtained by Ilinskii *et al*. [15]. When *ε*_{1} is increased, the resonance peak shifts to larger values of *R*_{0} and increases in magnitude. At *ε*_{1}=1, the resonance peak is found to be at *R*_{0}*f*=1.5 m s^{−1} and its magnitude is maximal. Note that for *ε*_{1}=1, the effective surface area by which the bubble acts on the wall is equal to the area of the lateral section of the bubble. Further increase of *ε*_{1} decreases the magnitude of the resonance peak, moving its position to larger values of *R*_{0}*f*, which, however, do not exceed 1.7 m s^{−1}.

Figure 5*b* shows the evolution of resonance curves with increasing |*ε*_{2}| at *ε*_{1}=1. Increasing |*ε*_{2}| can be related to an increase of damping in the wall material, which leads to flattening of resonance curves and decreasing resonance bubble radius. This behaviour is observed if *ε*_{2} is negative.

### (d) Experimental verification

The experimental data shown by circles in figures 6–8 were adopted from work by Mekki-Berrada *et al*. [31], where they were acquired by a microfluidic set-up described by Rabaud *et al*. [2]. This set-up has been already mentioned in the Introduction. The experimental data show the normalized amplitude of bubble oscillation for cylindrical bubbles of different radii, placed in a fluid channel, 25 μm in depth, confined by elastic walls made of PDMS. The bubble oscillation is induced by a vibrating glass rod moulded in the upper channel wall about 150 μm above the solid–fluid interface. The experimental data were obtained at three values of the driving frequency: 30, 40 and 50 kHz. For each frequency, the measurements were made for a certain range of bubble radii, which was dictated by experimental conditions.

The solid lines in figures 6–8 show the fitting of the experimental points by equation (2.88) with the fitting parameter *ε* given by equation (3.1). The material parameters were taken as pointed out in the preceding subsections. The theoretical curves were calculated at the following values of *ε*_{1} and *ε*_{2}: for 30 kHz, *ε*_{1}=1.76,*ε*_{2}=−0.451; for 40 kHz, *ε*_{1}=2.15,*ε*_{2}=−0.465; for 50 kHz, *ε*_{1}=1.95,*ε*_{2}=−0.756. The parts (*a*) of the figures show fitting in the experimental measurement range and the parts (*b*) show the view of the full theoretical curve. The normalization of the data is performed by a maximal value in a depicted range. For this reason, the normalization is different in figure 8*a*,*b*. The experimental data indicate that there is strong damping in the PDMS walls since the resonance curves are very flat. As one can see, agreement between the theoretical and experimental data is achieved if the values of the fitting parameters are changed with frequency. This result suggests that both effective surface area of the bubble–wall contact and the damping characteristics of the PDMS walls depend on frequency by a more complicated way than that specified by equation (3.1). This problem calls for further consideration.

## 4. Conclusion

A theoretical model has been developed that describes the volume oscillation of an ultrasound-activated cylindrical bubble situated in a fluid channel between two planar elastic walls. The model includes both the bulk scattered wave, which propagates in the fluid gap, and the surface waves of Lamb-type, which propagate at the solid–fluid interfaces. The force exerted by the bubble oscillation on the channel walls, which is the source of the Lamb-type waves, was treated as a normal harmonic point load. A dispersion equation for the above surface waves was derived, which allows one to examine the dependence of the speed of these waves on the channel height *h*. Such an equation was calculated for the first time, which makes it one of the main results of the present study. It shows that for *h*<λ_{t}, where λ_{t} is the wavelength of the transverse wave in the elastic walls, the speed of the Lamb-type waves decreases with decreasing *h*, while for *h*→λ_{t}, their speed tends to the Scholte wave speed.

The solutions for the wave fields produced by the Lamb-type waves in the elastic walls and in the fluid channel were obtained using the Hankel transforms and expressed in terms of improper integrals. A method was proposed for numerically evaluating the integral solutions. Numerical simulations were carried out to examine the effect of the surface waves on the bubble dynamics. It was shown that the resonance frequency of a cylindrical bubble confined between two identical elastic walls could be up to 50% higher than the resonance frequency of a similar bubble confined between two rigid walls.

Experimental verification of the theoretical model has been carried out using measured values of the bubble oscillation amplitude, which were obtained for bubbles of different radii at three values of the driving frequency. The fitting of experimental data was performed by using two parameters. The first parameter describes the magnitude of the load exerted on the channel walls by the variation of the bubble gas pressure. The second parameter describes the phase shift between the variation of the bubble gas pressure and the reaction of the walls. It was found that agreement between the theoretical and experimental data is achieved on condition that the fitting parameters are considered frequency-dependent quantities.

## Authors' contributions

A.D. carried out the development of the theory and drafted the manuscript; F.M.B. carried out the experimental measurements and helped draft the manuscript; P.T. participated in the experimental measurements and the interpretation of the results and helped draft the manuscript; P.M. participated in the development of the theory and the experimental measurements, coordinated the study and helped draft the manuscript. All authors gave final approval for publication.

## Competing interests

We have no competing interests.

## Funding

This research has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007–2013)/ERC grant agreement no. 614655 ‘Bubbleboost’.

## Appendix A. Calculation of the integral expressions

The denominator functions of the improper integrals (2.63)–(2.67) have zeros, which makes direct numerical integration impossible. This difficulty can be overcome by using contour integration and the residue theorem as proposed by Vrettos [32]. The calculation method will be demonstrated here with integrals (2.65) and (2.67) as these integrals are used in the equation of bubble oscillation.

Expressing the Bessel function *J*_{0} as a sum of the Hankel functions of the first and second kind
*ξ* with the complex variable *s*=*ξ*+*iτ*, one obtains

Note that the variable *ξ* is also replaced with *s* in the radical functions *s*=*ξ*_{f}, *ξ*_{l},1, where *ξ*_{f}<*ξ*_{l}<1 because it is assumed that *c*_{f}>*c*_{l}>*c*_{t}; see equations (2.68). The pole of the integrand *F*(*s*,*z*) is given by the root of the equation *ξ*_{s}, corresponds to a Scholte-type wave and *ξ*_{s}>1 as the speed of the Scholte wave is smaller than *c*_{t}; see the first of equations (2.68). To make the integrand single valued and take into account the presence of the pole, the integration contours *Γ*_{1} and *Γ*_{2} are taken as shown in figure 9.

Integrating *I*_{1} along the contour *Γ*_{1} and using the fact that

the superscripts + and − denote that the function *F*(*s*,*z*) is calculated for the positive and negative values of the imaginary part of the radicals _{p}(*ξ*_{s}) is the residue at the pole *ξ*_{s} which is calculated by [33]

Integrating *I*_{2} along the contour *Γ*_{2} and using the fact that

Equations (A 7)–(A 9) can be simplified using the fact that *F*^{−}(*ξ*,*z*)=[*F*^{+}(*ξ*,*z*)]*, where the asterisk denotes the complex conjugate. Further, the sum of the integrals *I*_{11} and *I*_{2} can be transformed using the following identities:
*K*_{n}(*x*) is the modified Bessel function of the second kind of order *n* [29]. As a result of these operations, the final expression for the integral *I*_{p}(*r*,*z*) is found to be
*x*} means the imaginary part of *x*. These integrals are nonsingular and convergent and hence they can be calculated numerically without any difficulty.

Application of the same method to integral (2.65) yields
*ξ*_{s} is calculated by

- Received January 13, 2016.
- Accepted March 11, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.