Wave Motion 17 (1993) 319-327 Elsevier
319
Coupled extensional and bending motion in elastic waveguides Douglas S. Drumheller Geothermal Research Department, Sandia National Laboratories, Albuquerque, NM 87185, USA Received 20 March 1992, Revised 22 December 1992
A variety of different linear elastic waves propagate in long slender rods. These waves do not interact in straight rods; however, extensional and bending waves do interact in slightly curved rods. If curved rods are used as waveguides for extensional waves, wave energy is exchanged between extensional and bending waves. The result has the superficial appearance of viscous attenuation. The petroleum and geothermal industries have an interest in using drill pipe as a waveguide to transmit data by acoustic waves. Attempts have been made to transmit data from transducers in the vicinity of a drill bit to the drill rig at the surface. These efforts have not been successful. Poor signal strength is a major problem. Historically, the observed reduction in signal strength is attributed to viscous dissipation. This paper proposes an alternative model - mode conversion from extensional waves to bending waves. This model is analyzed both as an eigenvalue problem and a transient finite-difference problem. The results predict measured attenuation levels, and account for many anomalous observations.
I. Introduction
Long and slender rods act as waveguides for acoustic waves with predominantly axial motion. For decades it was known that these extensional waves propagated with relatively low attenuation for thousands of meters in steel pipe. These phenomena were observed for pipe surrounded by both air and liquid. The oilfield drilling community is interested in exploiting this phenomenon. Drill strings used to drill deep wells are assembled from 10-meter long sections of hollow steel pipe. During drilling operations, navigation data is often collected at the bottom of the well and transmitted to the surface by means of modulated pressure pulses in the flow of the drilling mud. Very low frequencies are used (approximately 1 Hz), and the data transmission rates are extremely low. Modulation of extensional waves in the drill string itself offers the promise of data transmission at much higher frequencies (approximately 1 kHz), however, attempts to develop such a system have been unsuccessful. These efforts are discussed in the technical and patent literature as well as government reports; for
example, see [ 1-3]. The biggest impediment to the success of these projects is the lack of understanding of the mechanisms of propagation of extensional waves in long slender pipe. Drill strings do not have a uniform cross-sectional area, individual pipes are joined together at 10-meter intervals with heavy, threaded couplings called tool joints. These periodic and abrupt changes in cross-sectional area form a one-dimensional elastic lattice which exhibits the classical acoustical patterns of Brillouin scattering. Detailed analysis and observations of these phenomena are given by DrumheUer [ 1 ]. Random lengths of drill pipe alter these scattering patterns and produce attenuation of extensional waves, but only a portion of observed attenuation is explained by this [ 1-4]. The typical test uses a steel drill string about 125 mm in diameter and 400 m in length. Extensional stress waves are generated at one end of the string either by hammer blows or with continuous-wave transducers, and the elastic waves are measured at the opposite end. Reported estimates of attenuation range from as low as 7 d B / k m to as high as 60 dB/km, but there is no consistent correlation between attenuation
0165-2125/93/$06.00 © 1993 - Elsevier Science Publishers B.V. All rights reserved
320
D.S. Drumheller / E&stic waveguides
and frequency. Neither is there any correlation of the attenuation levels to the properties of the fluid surrounding the pipe, but it is reported that signal attenuation in new pipe is usually lower than that in worn pipe. The analysis of attenuation is limited and incomplete, however, one work deserves comment. Squire and Whitehouse [4] attempt to quantify two potential mechanisms, viscous drag and sound radiation. The two mechanisms are quite different. The axial motion of extensional waves produces a boundary layer effect in the fluid surrounding the pipe. Energy is dissipated through shearing motion of the fluid contained in this boundary layer. Squire and Whitehouse estimate 6 dB/ km attenuation at 1 kHz for a typical drill string in a typical drilling mud and negligible attenuation for a drill string surrounded by air. They then estimate the attenuation due to sound radiation into the surrounding fluid. This radiation results from the radial motion in the drill string due to the Poisson effect. Their attenuation estimate is more than 300 d B / k m at 1 kHz for a drill string surrounded by fluid, but again the attenuation is negligible for a drill string surrounded by air. This is also contrary to observation which shows that attenuation is neither strongly dependent upon the presence of fluid nor as large as this analysis suggests. The value of 300 d B / k m is actually one-quarter of the estimate in [4]. This adjustment is made because eq. (29) in [4] is missing a factor of 1/2, and due to the omission of the kinetic energy from the total energy expression, eq. (33) is missing a factor of 2. Furthermore, the use of this radiation term can be questioned on two points. First, the term is derived from uniform radial motion of an infinitely long cylinder, whereas the radial motion of the drill string is not uniform. As a periodic wave propagates along the drill string, half the drill string contracts while the other half expands, and the average radial motion of the drill string is zero. Far-field radiation is strongly coupled to this average motion; for example, see Morse [ 5 ] page 305. Second, the analysis applies to radiation into an infinite fluid, while the drill string actually radiates energy into fluid con:.ned to the small volume of the wellbore. One obvious source of attenuation has been overlooked - mode conversion from extensional motion to
bending and torsional motions. It is the objective of this work to examine this mechanism by limiting the analysis to long slender rods with curvature in only one plane. Coupling to torsional motion is omitted, and only large curvatures are considered. Furthermore, the analysis is limited to rods of constant density, elastic modulus, and cross-sectional geometry. Only the radius of curvature can change along the length of the rod. Thus Brillouin scattering phenomena is ignored in order to simplify the physics and isolate the essential features of this attenuation mechanism.
2. G o v e r n i n g equations
The flow of time is denoted by the variable t. In a two-dimensional Euclidean space, with spatial coordinates (x, y) the function f ( x ) denotes the reference line of the rod. The condition )9(x) = 0 defines a straight rod. For slightly curved rods
d2y dx 2
-
1 r(x) '
(1)
where r(x) is the radius of curvature of the rod. Equation ( 1 ) is only valid for large magnitudes of r(x). At any value of x, the tangent to the reference line is a normal to a plane that intersects the rod. The region of the plane contained within the rod is called the cross section of the rod. The parameter a denotes the area of the cross section. Only rods where a is independent of x are considered. The reference line is positioned so that the first moment of a about 39(x) is zero. The parameter I denotes the second moment of a about
y(x). Deformation of the rod is described by displacement of the reference line. The functions u(x, t) and ~/(x, t) are the components of displacement in the x and y directions. Deformation produces a normal-strain component acting in the direction of the reference line of the rod. To simplify the strain description, the usual kinematical assumption of beam theory is employed; that is, a plane cross section of the undeformed rod maps to a plane cross section of the deformed rod. The resulting axial strain e(x, y, t) is given by the familiar relationship
321
D.S. Drumheller / Elastic waveguides
0U
77
0277
E= ~x + -r + ( ) S - y ) Ox2 .
(2)
From Hook' s law, the axial stress is equal to Young' s Modulus E times the axial strain. The total force F ( x , t) is the integral of this stress over the cross section, F=Ea
+
(3)
,
and 7/. To simplify the resulting expressions, the sound speed c and the radius of gyration R are defined as
02r/ (4)
Ox----7 .
Figure 1 illustrates an element of the rod contained between two cross sections. The axial force F, shear force S, and bending moment M are depicted for each cross section. The balance of forces in the axial direction is OF 0x
S 02u r = p a Or2,
(5)
where p denotes the mass density of the material; the balance of forces normal to this direction is 0S 0x
F r
02r/. =pa --
Ol 2 '
(6)
and the balance of moments about points A is S=-
~M --~-x"
(7)
These relations can be reduced to two equations in u
(8)
R e = l/a.
(9)
and
Then eqs. (3), (4) and (7) are substituted into eqs. (5) and (6), to obtain
and the bending moment M ( x , t) is the first moment of the stress about )~(x), M=-El
c2=E/p,
02U ~0x
0 ['O~
R2 03r/
+ ~ trJ
r ~0x
1 02u -
c 2 0t 2 '
(10)
and _ R 2 04r/ OX4
l(0u ~) 1 02r/ r ~X + = c 2 0t 2 "
(11)
3. Harmonic analysis
For an infinitely long rod with uniform curvature, the harmonic solutions to eqs. (10) and ( 11 ) are (u, r / ) = (Uo, r/o) exp i(oJt+kx)
(12)
where i 2 = - 1, ¢o is the circular frequency, and k is the wave number. The real part of the wave number R e a l ( k ) is related to the wavelength A by Real (k) = 2~r/h. Substitution of solutions (12) into eqs. (10) and ( 11 ) yields the following eigenvalue problem: +_i(l+k2R2)k/r ][Uo] W2/C 2 - 1/r2--RZk4JLr/o_] = u "
0)2/c2-k Tik/r
(13) For large values of r, the eigenvalues of the determinant of the coefficients of eq. (13) are
I
s+(os/ox)~x
7"
I
+
Fig. 1. Elemental volume of the rod.
f o)/c k='~ ~
.
(14)
Li The second equation of (13) gives r/__£o= 5: ikr UO _kaR2r2 + to2r2/c2_ l
(15)
D.S. Drumheller / Elastic waveguides
322
When the first eigenvalue in eq. 14) is substituted into this expression, rh~ = Uo
+_ikr -k4RZrZ+kZr
2-
(16)
1
resents a stationary bending wave. For both eqs. (19) and (22), the ratios of the components of the motion are equal to kr= 2rrr/A. A rod with uniform curvature forms a ring. The parameter kr represents the ratio of the circumference of this ring over the wavelength.
For wavelengths where r >> A >> R, (kr)2>> 1 ,
(17)
4. Transient analysis
and (kR)2<< 1 .
(18)
Equation (16) becomes u~ =T-ikr ~o
fork 2=
.
(19)
This expression coupled with eq. (17) show that the first eigenvalue of eq. (14) represents a wave with predominantly extensional motion. The phase velocity c~ of this wave is given by the ratio oJ/k, or ca = c .
(20)
The group velocity Ca of this wave is given by dw/dk, or
Ca = c .
(21)
Since both the phase and group velocities are equal and constant, waves of this type with multiple frequency components will not disperse. The first root represents a steady propagating extensional wave. When the last two eigenvalues of eq. (14) are substituted into eq. (15), rl~° = -T-ikr Uo
for k 4= R c
"
(22)
For these waves, the motion is predominantly bending. Since the second root is real, the phase velocity Cb and group velocity Cb are given by Cb= WV~C,
(23)
C b = 2 wgrwRcc.
(24)
and
The field equations and geometry of this problem are relatively simple; however, one aspect of the analysis of these equations is an especially difficult issue. This is the accurate calculation of relatively small levels of wave attenuation. Many analysis methods produce numerical dissipation and dispersion which exceeds the physical attenuation levels of this linear problem. This is especially true for conventional finite-difference and finite-element methods. The algorithm outlined here is an adaptation of the work by Drumheller [6] which is designed to overcome these difficulties. Both the axial position x and time t are discretized into uniform intervals by defining the dimensionless parameters n a n d j so that x = x ( n ) and t = t ( j ) . Any function F(x, t) can then be recast as F(7, J)- Two difference operators A and 8 are defined as A F ( n , j ) = F ( n + 7,J) ~ ' - F(n-~_,~ j) ,
(25)
aF(n,j)=I~(n,j+½)--I~(n,j--½)
(26)
and
Repeated application of these operators is denoted as A k+ I F = A(AkF), and a k+ I F = a ( a * r ) . The variables x and t are required to be linearly dependent upon n and j: x = n A x + x ( O ) , and t = j & + t(0). Because of this linear relationship, the quantities ~ and at are constants. They are often called the mesh intervals of the finite-difference algorithm. In general, each can be arbitrarily specified; however, the key to the present method is the application of the following constraint (see [ 6 ] ): gx=c&.
Because of the frequency dependence of these velocities, the second root represents a dispersive propagating bending wave. Since the third root is imaginary, it rep-
.
(27)
For convenience, the radius of curvature and the radius of gyration are normalized to the mesh interval Aac.
323
D.S. Drumheller / Elastic waveguides
R=R/Ax.
(28)
F=r/Ax.
(29)
The finiS-difference approximations of the partial derivatives with respect to n andj are OkF -= AkF,
(30)
and g2-{2-4/~2[l-cos(k)]2}g+ 1=0. For all real values of k,
OkF
(31)
--
ojk = S k F .
By using eq. (27), the chain-rule expansion, and these finite-differenceapproximations, the equations of motion (10) and ( 11 ) become A2U --[-A
-- 7
A3~ ~ ~2u'
(39)
Ial = 1 , and Igl~
and
(38)
for R~< 1/2.
(40)
Equation (39) shows that the solution for the extensional motion is unconditionally stable. This result is directly attributable to the constraint (27) placed on the mesh intervals. In contrast, eq. (40) indicates that the bending solution is stable or decays only when the mesh interval Ax is greater than twice the radius of gyration R.
(32) 5. Nondimensional quantities
and (33) These relationships are used to evaluate u and r/at integer values of n andj. This is accomplished by direct expansion of the difference operators coupled with the averaging approximation
The mesh intervals Ax and & provide convenient scales for defining nondimensional quantities. The nondimensional radius of curvature g and the nondimensional radius of gyration/~, have already been defined using Ax. In a similar fashion, nondimensional quantities are defined for axial displacement ti, radial displacement 9, frequency t3, and wavelength ,~.
(34)
ff=u/Ax .
(41)
The result is an explicit difference algorithm which is applied in a cascading fashion to compute u and r/at j + 1 from known values atj and j - 1. This finite-difference algorithm is conditionally stable. To examine the conditions that influence stability, the behavior of a specific set of solutions is studied. These solutions are
r~= r//Ax.
(42)
o3= to&.
(43)
X=A/Ax.
(44)
2F(n,j)=F(n+
1,~ j : ) + F ( n - ½ , j ) .
u(n, j) = G j exp(ikn) ,
(35)
7/( r/, j) = gJ exp(ikn).
(36)
and
This leads to the nondimensional wave number/~, /~= 2~r/X = kAx.
Nondimensional expressions for the phase and group velocities of the propagating extensional and bending waves are:
h the algorithm is simplified by ignoring the terms contaiaing the curvature 1 /g, substitution of these solutions int~ eqs. (32) and (33) gives G Z - 2 cos(k)G+ 1 = 0 ,
(37)
(45)
and
~a = o3//~= 1,
(46)
t~a = do3/d/~= 1,
(47)
Cb -~" ~ ' ~
(48)
,
324
D.S. Drumheller / Elastic" waveguides
Cb = 2 ~ - R .
(49) "" I
/fi!
,1;: 5 F
I
I '1
J i
'
I
'
6. Results
The accuracy of the finite-difference algorithm for a straight rod with pure extensional motion is documented in Drumheller [6]. However, an illustration of the algorithm for pure bending of a straight rod is not contained in that work. This section begins with such an example. Then solutions for both extensional and bending waves in a uniformly curved rod are compared to results from the harmonic analysis. This is followed by a treatment of waves in a rod with continually changing curvature. 6.1. Pure bending in a straight rod
A rod occupies the region 0 ~
,jt
0 ~;:
'li
'L
/
/"i
j
i
•u7
0 ~' k
' I,'1'
')
!,~1 ;'~O,,
'] ;, r:c:i,zeJ ;r-,c 17: Fig, 2. Displacement history at n
=
50.
main wave packet. Within the central portion of the packet, the period of the wave is j = 72, and equal to the period of the driving signal. In uniform materials without dispersion, the displacement history measured at a given position is quite similar to the displacement profile measured at a given time. However, because the group and phase velocities differ by a factor of 2, the characteristics of these two plots for the present example are quite different. Figure 3 illustrates the displacement profile at the timej = 780. The result is plotted against the normalized position n. The two vertical dashed lines shown in this plot again represent the positions of the signals created at the beginning and the end of the 5-cycle disturbance. Approximately 10 waves appear in the central portion of the wave packet. The wavelength of this signal matches the predicted value of n = 12, and the ratio of 10 cycles in Fig. 3 to 5 cycles in Fig. 2 is directly connected to the factor of 2 difference between phase and group velocity. 6.2. Rods with uniform curvature
A uniformly curved rod is now subjected to a displacement of the left boundary. A bending motion r/is applied first. Except for the curvature of the rod, this problem is identical to the previous problem. However, because of the curvature, the bending and extensional motions are now coupled. Each of the three types of waves described in Section 3 are stimulated by this
325
D.S. Drumheller/ Elastic waveguides ] .0
"
~
,--T
-
i i
oo
q
--~ \/",,/1 i
tl lil
3
k
10 [
i
.....
0
20C
400
Naturalized
Position
600
(~.)
Fig. 3. Displacementprofileatj = 780.
I
120xu T/
+d o @ I(3 o_ u~
-1
dE] 2
--3E 0
I
,
200 Normalized
i
J
I
400
600
Position
800
ets are evident. The propagation velocities of these packets closely match the predicted group velocities. Both packets also have the same dominant frequency. The slower bending wave contains the majority of the energy produced by the boundary condition. The ratio of the bending to extensional motion closely matches the predicted value of 120. These motions are also 90 ° out of phase. The faster traveling extensional wave has only a small fraction of the total energy. The magnification of the extensional motion does not allow a visual comparison to the coupled bending motion of this wave; however, the ratio of the magnitudes and phases for this wave also match the predictions of the harmonic analysis. The solution is now recomputed with a new boundary condition. It is changed to one involving pure extensional motion, u. Again, because of the curvature, the extensional driving motion is coupled to bending motion. Figure 5 contains the results. In this case the dominant motion is extensional, and the bending motion has been magnified by the factor kr = 20. The group velocities, amplitude ratios, and phase relationships between the motions closely match the predicted values of the harmonic analysis. 6.3. Rods with changing curvature
(~)
Fig. 4. Displacementprofile atj = 750 for applied bendingmotion. boundary condition. The radius of curvature is F= 270/~r, which is 720 times greater than the radius of gyration. The boundary displacement r/is a 5-cycle sine wave with unit amplitude and frequency 03= ~r/ 36. For this frequency, the phase and group velocities for the two propagating waves are ?a = Ca = 1, and 2 ? b = C b = l / 3 . From eq. (19) for the extensional wave, the ratio of the extensional motion and the bending motion is kr=20. From eq. (22) for the bending wave, the ratio of the bending motion to the extensional motion is kr = 120. For both types of waves the motions are 90 degrees out of phase. Figure 4 contains calculated displacement profiles at the timej = 750. The extensional motion has been magnified by the factor k r = 120 as means of comparing these results to the harmonic analysis. Two wave pack-
Both the harmonic analysis and the calculations for uniformly curved rods describe extensional waves which propagate without dispersion or energy loss. While coupling to bending motion is present, exten2
r
r
20 x'r/
•
c ~
0
:
::!
o cl_
:
?5
7
:. ::
....
.
i ii: -2
i
0
i
200
i
400 Normalized
Position
600
800
(n)
Fig.5. Displacementprofileatj = 750forappliedextensionalmotion.
326
D.S. Drumheller / Elastic waveguides
sional waves propagate as steady waves. In order to establish curvature and coupling to bending motion as a significant source of distortion and attenuation of extensional waves, an additional feature is added to the problem. That feature is irregular curvature. The harmonic analysis illustrates that the eigenvectors of both the extensional and bending waves depend upon the curvature. Transient calculations also show that an abrupt change in curvature causes a partial reflection of both bending and extensional waves. Therefore continually changing curvature will produce continually changing waveforms. To illustrate this effect, the following curvature is selected: 1
1
#(n)
#w
sin(2"rrn / W) ,
(50)
where the constants ?w = rw/zix and W are arbitrary. Integration of eq. ( 1 ) yields
Y(n)(W) rw = ~
2 sin(2"rrn/W).
(51)
This curvature expression is applied to the rod described in the previous section. The selected curvature parameters are W = 10, and ~w= 720/~r. From eq. (51 ), the maximum value of y ( n ) for these parameters is 0.035R. Thus, the deviation of the reference line from a straight line is only a small fraction of the rod diameter. An extensional motion is applied to the left boundary. It is chosen to be a half cycle of a sine wave with a unit amplitude and a period of 18&. The solution is computed to j = 800 and the displacement histories at n = 400 are recorded. The digital Fourier transforms H(u) and H(r/) of the resulting histories are shown in Fig. 6. These calculations are labeled as W = 10, and the frequency is normalized to the frequency of an extensional wave with the same period as the curvature, tOw= 27r/10. An additional calculation is also shown. For this calculation, the curvature of the rod is the sum of two sine components. The first component has the original parameters, W = 10 and ~w=720/Tr, and the second component has the parameters, W = 7 and ew = 720/'rr. The general trend of H ( u ) to decrease with increasing
5 0 [ - - ~
....
,
,
020 Frequency
(~/~.)
,
20'
000
0 10
050
Fig. 6. Fourier transforms of histories at n
0 40
=
0 50
400.
frequency is a property of the original half-cycle waveform of the boundary condition. In fact, at all but the three distinct dips in these profiles, the calculated amplitudes are nearly equal to the amplitude of the boundary condition. The dips in the curves for extensional motion correspond to peaks in the response of the bending motion. These are frequencies where the curvature results in strong coupling between the extensional and bending motions. For the simple curvature with one sine component, two peaks in the bending response are evident. The coupling to bending is much stronger at the higher frequency. For the more complex curvature with two sine components, the original peaks are still present, but two additional peaks also exist - a peak at 0.94 is not shown. Evidently, increasing the complexity of the curvature of the rod broadens the frequency range over which the strong coupling between extensional and bending motion occurs. An additional calculation using the curvature with a single-sine component will help to illustrate this phenomenon further. For this calculation, the boundary condition is changed to a 5-cycle sine wave. The amplitude is again unity, but the period is changed to 29.5&. This concentrates the energy of the applied extensional motion into a narrow frequency band centered at the normalized frequency w~tow = 0.34. The time-domain results for j = 200 and j = 750 are shown in Fig. 7. At j = 200, the extensional wave still retains the 5-cycle shape of the boundary condition, and only a small amount of bending motion is present. However, at
D.S. Drumheller / Elastic waveguides
1.0
•
.
i
'
'
q
'
'
F
0.5
....
0.0 09 0 0_ CO
05
-1.0 0
200 400 600 Normalized Position (n)
800
Fig. 7. Displacement profiles forj = 200 andj = 750.
j = 750, the extensional wave is distorted into an unusual form, and the bending motion has greatly increased. Because the bending waves propagate at 1/3 of the speed of the extensional waves, this cloud of bending wave energy increases in length in the region behind the extensional wave. The bending motion represents a loss of energy from the extensional wave. In this example the extensional wave has an energy loss of about 1.5 dB over a distance of 750At.
7. Conclusions The purpose of this work is to demonstrate the relationship between extensional and bending motion in long slender rods with slight curvature. The motivation behind the work is to demonstrate that the observed attenuation of extensional waves in oilfield drill pipe can be explained by this relationship. This point is demonstrated, albeit in the somewhat cryptic fashion of nondimensional solutions. A specific dimensional example is in order. Given the worldwide oilfield use of English units, this example is presented in inches and feet. Five-inch diameter steel drill pipe is in common use for off-shore drilling. It has a cross-sectional area of 6.63 in 2 and a radius of gyration of 1.8 in. The results of Section 6.3 can be applied to this pipe by selecting A t = 1.8~r =5.67 in. The single-component curvature then has a periodic length of 10At=56.7 in, and a maximum eccentricity of 0.035R -- 0.063 in. The extensional wave speed of steel is 202,000 in/s. Therefore
327
& = A t / c = 28 I~s, and the dominant frequency of the applied 5-cycle boundary function is f = 1/29.5&= 1.2 kHz. The propagation distance used in this calculation is 750At = 354 ft. Therefore the calculated attenuation is 1.5 dB/354 ft = 4.24 dB/1000 ft = 13.9 dB/km, which is well within an acceptable range. In this simple model, wave scattering from pipe joints is ignored, and curvature is restricted to a single spatial plane. Even so, a variety of phenomena are explained: 1. Calculated attenuation levels roughly correspond to observation. 2. Attenuation is present even when drilling fluids are absent. 3. As illustrated by Fig. 6, there is no simple correlation between attenuation and frequency. 4. The model predicts a dramatic shift of wave profiles from simple wave packets at certain frequencies to the waves illustrated in Fig. 7. This makes the characterization of attenuation for a particular drill string extremely difficult and confuses efforts to correlate observations from different drill strings.
Acknowledgment This work was supported by the U.S. Department of Energy at Sandia National Laboratories under Contract DE-AC04-76DP00789.
References [1] D.S. Drumheller, "Acoustical properties of drill strings", J. Acoustical Society of America 85, 1048-1064 (1989). [2] W.H. Cox and P.E. Chaney, Telemetry system, U.S. Patent No. 4,293,936, 1981. [3] W.H. Kent, P.G. Mitchell and R.V. Row, Geothermal downwell instrumentation (during drilling), Sperry Research Center, Sudbury Massachusetts 01776, June 1979, SCRC-CR-79-41. [4] W.D. Squire and H.J. Whitehouse, "A new approach to drillstring acoustic telemetry", Society of Petroleum Engineers of AIME (1979) SPE 8340. [5] P.M. Morse, Vibration and Sound, McGraw-Hill, New York (1948). [6] D.S. Dramheller, "Time-domain computations of onedimensional elastic waves in liquids, solids, and ferroelectric ceramics", Wave Motion 17, 63-88 (1993).