## Abstract

The dynamic propagation of a Heaviside input signal in a semi-infinite, dissipative/dispersive medium is considered. The exact solution to this problem, which corresponds to Stokes' first problem of fluid mechanics, is obtained and analysed using integral transform methods. Special/limiting cases are noted and numerical methods are used to illustrate the analytical findings. Specifically, the following results are presented: (i) critical values of the dispersion coefficient are noted and examined; (ii) for large values of time, the solution exhibits Taylor shock-like behaviour; (iii) the half-peak point of the Taylor shock exhibits a phase shift that depends on the dispersion coefficient. Lastly, links to other fields are noted and some associated mathematical relations are presented.

## 1. Introduction

Love (1944), using a variational approach, showed that if lateral displacements are taken into account, then the longitudinal motion of a slender, elastic rod is described by the fourth-order partial differential equation (PDE)(1.1)which is known today as *Love's equation* (see also Graff 1991). Here, is the displacement vector; the *x*-axis coincides with the centreline of the rod; *t* denotes time; is the wave speed, where the positive constants *E* and , respectively, denote Young's modulus and the mass density; the constant *ν* is Poisson's ratio; and the constant denotes the radius of gyration of a cross-section about the centreline. Additionally, it is assumed that the material of which the rod is made behaves elastically and follows Hooke's law, namely(1.2)where and , respectively, denote the axial stress and strain components.

It should be noted that Davis (1948) and Thomson (1960) obtained the exact solutions to equation (1.1) for the case of a finite rod and semi-infinite rod, respectively, subjected to a Heaviside jump in the stress at one end. Later, van Wijngaarden (1968) showed that equation (1.1) governs acoustic propagation in inviscid bubbly liquids. Whitham (1974) noted that it also arises in both the study of plasma waves and in the linear theory of water waves under the Boussinesq approximation for long waves.

In this article, we consider, in a general mathematical context, a dissipative version of equation (1.1), specifically, the PDE(1.3)Here, the positive constants , *ϰ* and *L* denote the wave speed, the damping coefficient and the dispersion coefficient, respectively, and we observe that *L* carries units of length. It should be noted that the inclusion of dissipation breaks the intrinsic temporal symmetry exhibited by equation (1.1), as evidenced by the presence of the third-order mixed derivative term in equation (1.3). Thus, unlike Love's equation, equation (1.3) is *not* invariant under the transformation .

The best-known context in which equation (1.3) arises is in van Wijngaarden's (1972) theory of acoustic propagation in viscous, isothermal bubbly liquids. Eringen (1990) re-derived the multi-dimensional form of this PDE based on a microcontinuum theory. Rubin *et al*. (1995) found that equation (1.3) also describes acoustic waves in a thermoelastic compressible Newtonian viscous fluid.1 Hayes & Saccomandi (2000) showed that it also governs the propagation of transverse plane waves in a particular class of viscoelastic media. Subsequently, Jordan & Feuillade (2004) obtained the exact solution to equation (1.3), which they termed the van Wijngaarden–Eringen equation, in the context of the compressible version of Stokes' second problem (see Schlichting 1968). More recently, Jordan & Feuillade (2006) also obtained the exact solution of this PDE in the context of the compressible version of Stokes' first problem (Schlichting 1968), but only for cases in which the dimensionless bubble radius, corresponding to the dimensionless dispersion coefficient of the present work (see equation (2.3)), is less than 1/2.

The primary aim of the present communication is to extend the work of Jordan & Feuillade (2006) by presenting the complete solution to the problem, i.e. one that is valid for arbitrary positive values of the dispersion coefficient. In doing so, we also generalize the work of Thomson (1960), which involved equation (1.1), and obtain the exact general solution to the corresponding problem for equation (1.3). In addition, we point out special/limiting cases and derive (possibly) new mathematical expressions related to the *J*-function. To this end, the present article is arranged as follows. In §2, we obtain the exact solution to the above-mentioned initial-boundary value problem (IBVP) using integral transform methods. In §3, we present a variety of analytical results including special/limiting cases of the solution. In §4, numerical results are presented. Finally, in §5, conclusions are stated and in §6 some open problems are noted.

## 2. Mathematical formulation and solution

### (a) Stokes' first problem: generic formulation

We consider the following dimensionless IBVP involving equation (1.3) in the half-space *x*>0:(2.1)(2.2)where denotes the Heaviside unit step function, , and the dimensionless dispersion coefficient is given by(2.3)In addition, we used the following non-dimensional transformations: , and , where the constant is the amplitude coefficient of in the dimensional formulation, and where all primes have been omitted but are understood.

### (b) Exact general solution using integral transform methods

To solve the IBVP, we employ a dual transform approach (e.g. Duffy 2004). Specifically, we first apply the spatial sine transform, thus reducing equation (2.1) to an ordinary differential equation (ODE) with independent variable *t*, and then use the temporal Laplace transform to solve this ODE. This yields, after simplifying, the dual transform domain solution(2.4)Here, *ξ* and *s* are the sine and Laplace transform parameters, respectively, a hat (resp. bar) superposed over a quantity denotes the image of that quantity in the sine (resp. Laplace) transform domain, and(2.5)From equation (2.5), it is apparent that , where , is a critical value of the dimensionless dispersion coefficient. This means that the roots are either complex for all *ξ*>0, when ; or complex on the finite interval , and real on , when , where . Hence, as we will soon see, the temporal part of our solution has a different analytical character for than it does for .

Obtaining first the Laplace inverse of equation (2.4) using a table of inverses (e.g. Carslaw & Jaeger 1963; Duffy 2004), multiplying the result by , and then integrating with respect to *ξ* from zero to infinity, we find the exact *xt*-domain solution to be(2.6)where(2.7)(2.8) and are given by(2.9)and we observe that *ξ*^{*} is the integration breakpoint.

## 3. Analytical results

### (a) Special case solution:

Returning to equation (2.1), and applying *only* the temporal Laplace transform both to it and the boundary conditions given in equation (2.2), we find, on employing the initial conditions and solving the resulting subsidiary equation, that the (Laplace) transform domain solution is given by(3.1)where(3.2)Now, if , then and equation (3.1) reduces to(3.3)Inverting equation (3.3) using a table of inverses (Luke 1962) gives us(3.4)where denotes the *J*-function (e.g. Goldstein 1953; Luke 1962) and denotes the modified Bessel function of the first kind of order *q*. If, instead, we invert equation (3.3) using the Laplace inversion formula, making use of the residue theorem (e.g. Duffy 2004) to evaluate the resulting contour integral, we obtain the double-series solution(3.5)which through the use of equation (A 2) can be expressed more simply as(3.6)where is the incomplete Gamma function (see appendix A).

Let us now, assuming *t* to be very large, use equation (A 4) in equation (3.6), where *κ* and *ζ* become *m* and 2*t*, respectively. Then, using the fact that(3.7)and observing that each term of this series is positive, it is evident that(3.8)for every positive integer *M*. Clearly, this difference can be made less than any given positive number by taking *M* sufficiently large. Consequently, it is readily established that for large *t*, equation (3.6) behaves like(3.9)Next, we observe that equation (3.6) can be re-expressed as(3.10)If is now taken to be very small, then, using equation (A 4) once again, but this time in equation (3.10), where *κ* and *ζ* now become *m* and 2*x*, respectively, one can, using an approach similar to that described above, obtain the small-*t* approximation (P. Amore 2006, personal communication)(3.11)which we observe does *not* satisfy the boundary condition (BC) at *x*=0.

### (b) Temporal limits and small-time behaviour

As shown by Jordan & Feuillade (2006), , for all . This follows from the fact that and . They also showed that at start-up, *p* suffers a jump discontinuity of amplitude(3.12)where we observe that the initial condition is *not* satisfied from the right.

Let us now examine the evolution of *p* immediately after start-up; i.e. immediately after assuming the profile . This we do by expanding the argument of the exponential in equation (3.1) in powers of 1/*s* and then neglecting terms . After inverting term-by-term, we find that for sufficiently small, *p* is approximately given by(3.13)This expression, which is a refinement of the (one-term) small-time approximation given by Jordan & Feuillade (2006), in fact, satisfies the BC at *x*=0, unlike equation (3.11). In addition, we observe that for and , where , equation (3.13) reduces to and , respectively.

### (c) Limiting cases of *R*

Since a number of other studies devoted to the (i.e. ) limiting case of equation (2.1), which is known as Stokes' equation, appear in the literature (e.g. Jordan *et al*. 2000 and references therein), we will do no more here than mention that in the limit , equation (2.8) reduces to(3.14)where and . We should also mention that, since it was obtained using a dual transform approach, equation (3.14) is possibly a new form of the solution to the present problem involving Stokes' equation.

If, on the other hand, we allow *R* to become very large, as a result of , then it can be shown that(3.15)If we now neglect terms of , then and it follows that equation (3.1) becomes(3.16)which we note is the exact Laplace transform domain solution to the present IBVP involving equation (1.1) in a suitable dimensionless form2 (see also Thomson (1960)). Expanding equation (3.16) in powers of 1/*s*, and then inverting term-by-term, it is not difficult to obtain the small-*t* (relative to *R*^{2}) expression(3.17)which clearly exhibits an extremely weak *t*-dependence. If, instead, we invert equation (3.16) using the Laplace inversion formula, we find that(3.18)where and . This is the dimensionless form of the solution given by Thomson (1960) for the case of equation (1.1).

### (d) New expression for the *J*-function?

Finally, using the fact that the r.h.s. of equation (3.4) must be equal to that of equation (3.6), it is not difficult to establish the relation(3.19)To the best of our knowledge, this is a new expression for the *J*-function. Numerical evaluations of equation (3.19) for various values of *η*, (not shown) indicate that this series converges very rapidly, as might be inferred from the factor in the denominator, and is therefore an efficient method for computing the *J*-function. In the same way, we can also evaluate the integral that appears in equation (3.4):(3.20)Of course, using equations (3.9) and (3.11), one can derive parallel asymptotic expressions for the quantities in equations (3.19) and (3.20).

## 4. Numerical results

In figures 1–3, the curves for the cases were calculated using the special case solution given in equation (3.6), while the other curves shown were calculated using the general solution given in equation (2.6). In addition, the graph of , which is the exact solution to the corresponding IBVP(4.1)(4.2)supplemented with Sommerfeld's out-going radiation condition, is superimposed for reference.

We should also note that and that we verified our computational results by numerically inverting equation (3.1) using Tzou's Riemann inversion algorithm (Tzou 1997), i.e.(4.3)for , where and denotes the real part of a complex quantity.

### (a) Temporal evolution and the effects of varying *R*

In figure 1, we have plotted *p* versus *x* for four different values of *R*, and a relatively small value of . Note that *p*=1 when *x*=0, as expected, and that, for small values of time, *p* is an increasing function of *R*, as may also be inferred from equation (3.13).

Figure 2 was plotted for the same values of *R* used in figure 1, but for a larger value of . As *t* is increased, two phenomena are seen to emerge. First, the highly exponential character of the profiles observed in figure 1 changes, such that they begin to develop a shoulder, a behaviour previously reported by Jordan & Feuillade (2006) for the special case . Second, while the development of the shoulder occurs more rapidly for smaller values of *R*, as may be discerned in the curve corresponding to in figure 2, our present analysis shows that the same phenomenon occurs for other values of *R*. In figure 2, this is seen to lead to a crossover of the profiles, such that *p* becomes a *decreasing* function of *R* to the left of the crossover point.

### (b) Taylor shock-like behaviour and half-peak point

In figure 3, the behaviour of *p* versus *x* is examined for a much larger value of *t*=100. In this case, however, for reasons of clarity, curves are plotted for only three of the four values of *R* used in figures 1 and 2. For such a large-time value, the shoulder seen emerging in figure 2 is now fully developed, and each of the profiles has evolved to take a form similar to a *Taylor shock*3, a finding reported earlier by Jordan & Feuillade (2006) for the case . In particular, we observe that as *x* increases, *p* falls from very nearly 1 to essentially zero, over a relatively narrow transition region, for every value of *R* considered.

We also see, by comparing the three curves plotted in figure 3, that the speed (and shape) of the profiles, as they penetrate deeper into the half-space, is virtually independent of *R*. This is shown more explicitly in figure 4. Here, the *x*-coordinate of the half-peak point , which is defined as the solution of the equation(4.4)where *t* and *R* are assumed to be given, is plotted as a function of *t* for the same three values of *R* used in figure 3. While a minor departure from linear behaviour is observed for small values of *t*, due to the time delay as the Taylor shock develops, in each case, it is evident from the slopes that these waveforms propagate with a speed that tends to unity as *t* is increased; i.e. , as , for all values of *R* considered, where a superposed dot denotes d/d*t*. Note also that there is a slight, but noticeable, *R*-dependent *phase differential*, such that the half-peak points of the profiles corresponding to smaller *R* values are shifted further to the right.

### (c) Depth of influence

In figure 5, we have plotted versus *R*, where the ‘depth of influence’ satisfies the equation(4.5)for various values of *t*, together with the limiting case (see equation (3.12)). Clearly, is an increasing function of both *R* and *t*. Moreover, the large-*R* behaviour of appears to be both linear, i.e. as , and virtually independent of *t*.

## 5. Conclusions

As was previously shown to be true for the case of Stokes' second problem involving equation (1.3) (Jordan & Feuillade 2004), and , where and , are critical values of the dimensionless dispersion coefficient.

When , the exact solution can be expressed in terms of the

*J*-function (see equations (3.4)–(3.6)).For , the small-time expressions given in equations (3.11) and (3.13) provide complementary approximations of

*p*for large and small values of , respectively.As

*R*becomes very large, i.e. , the solution given by Thomson (1960) for the undamped rod is recovered (see equation (3.18)).While

*p*initially exhibits a decaying exponential character, numerical studies indicate that*p*begins to behave in a*wave-like*manner as*t*increases, at least for values of*R*in the range . Specifically, a shoulder forms and the profile assumes a shape similar to that of a Taylor shock as (see figures 1–3).For the values of

*R*considered, it was found that the propagation speed of the Taylor shock's half-peak point rapidly tends to unity as (see figures 3 and 4); in other words, , as , in terms of the dimensional variables. Additionally, the Taylor shock profile assumes an approximate ‘’ shape as*t*increases (see figure 3).During the Taylor shock stage, a slight, but noticeable,

*R*-dependent phase differential appears, such that the half-peak points of profiles corresponding to smaller*R*values are shifted further to the right (see figures 3 and 4).As suggested by figure 5, is, practically speaking, independent of

*t*for large*R*. Additionally, the versus*R*profile suggests the asymptotic relation(5.1)which follows from solving the equation (see §4

*c*).

## 6. Closing remarks

While the present work has successfully addressed a number of issues regarding the solution of IBVP (2.1) and (2.2) that were not covered by Jordan & Feuillade (2006), there are still several aspects of this problem, which could form the basis for future investigations, that remain either unresolved or unaddressed. For example:

Determine, for arbitrary

*R*>0, the exact inverse of equation (3.1) using the Laplace inversion formula.Derive, using only analytical methods, an asymptotic expression which accurately represents the Taylor shock stage of the profile. (Here, we should mention that a numerically motivated

*ad hoc*expression for the Taylor shock stage was presented by Jordan & Feuillade (2006).)Again using only analytical methods, derive an expression for the phase shift exhibited by the half-peak point of the Taylor shock profile.

Examine the effects of non-zero initial conditions.

## Acknowledgments

The authors thank Drs Stanley A. Chin-Bing and Edward R. Franchi for their encouragement and support. The authors also thank Prof. P. Amore for his kind assistance. This research was supported by ONR/NRL funding (PE 061153N). All numerical work was carried out using the software package Mathematica (v. 4.0).

## Footnotes

↵Here, we take issue with the description ‘thermoelastic compressible Newtonian viscous fluid’, since fluid viscosity+elasticity implies viscoelastic, i.e. non-Newtonian behaviour. What Rubin

*et al*. (1995) appear to be describing is a thermally conducting generalization of the*isothermal*bubbly liquid medium of van Wijngaarden's (1972) theory.↵Note also that the start-up jump (3.12) also occurs when damping is

*absent*.↵Strictly speaking, ‘Taylor shock’ refers to an inherently nonlinear phenomenon (e.g. Dodd

*et al*. 1982). However, because of its remarkable similarity to an actual Taylor shock, and for convenience, we have adopted this term to describe the large-time behaviour of the solution profile in the present work.- Received November 19, 2005.
- Accepted January 13, 2006.

- © 2006 The Royal Society