- Email: [email protected]
Reflection of an impulsive torsional wave by a moving boundary
Wave Motion 51 (2014) 1209–1224
Contents lists available at ScienceDirect
Wave Motion journal homepage: www.elsevier.com/locate/wavemoti
Reflection of an impulsive torsional wave by a moving boundary Kazumi Watanabe Sakuragaoka_1-4, Hikoshima, Shimonoseki, 750-0083, Japan
highlights • • • • •
Reflection problem of an impulsive torsional wave by a moving boundary is discussed exactly. A torsional wave source which produces the simple spherical wave is discussed. Exact closed form solution for the reflected wave is obtained. A new wave is found. Development of wave front shape is analytically obtained and illustrated.
article
info
Article history: Received 30 March 2014 Received in revised form 14 June 2014 Accepted 9 July 2014 Available online 17 July 2014 Keywords: Moving boundary Torsional wave Reflection Closed form solution
abstract Reflection of an impulsive torsional wave by a uniformly moving surface/boundary is discussed. Applying the Cagniard–de Hoop technique, the solution for the reflected wave is obtained in the closed form and a new wave that runs along the moving boundary has been found. Its wave front singularity is stronger than those of the incident and the regular reflected waves. This new wave has the cylindrical surface, but the shape of the disturbed region is a torus with a triangular-like cross section. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Wave reflection at a moving object/boundary is the substantial sensing principle for the radar systems. Numerous works have been carried out for electromagnetic waves (for example, [1–3]). The similar Doppler effects for acoustic waves have also been discussed theoretically and experimentally [4–6]. Vibration and waves in ropes and rods with moving supports are also important engineering problems and thus the vibration problem for the high speed elevator and the power transmission have been discussed [7–9]. However, as for the 2D and 3D elastic waves, not so many works have been done for the reflection problem by the moving boundary. Reflection of the 2D time-harmonic wave at a uniformly moving boundary and that of an impulsive SH-wave have been discussed by the author [10–12]. When the impulsive wave is reflected by the uniformly moving boundary, a new additional wave is produced and it runs along the moving boundary [11,12]. This new wave has been found in the case of the 2D anti-plane deformation. In order to confirm the appearance of the new wave for the 3D dynamic deformation, the reflection of the simplest 3D wave, i.e. the torsional wave, is discussed. This is the motivation of the present work. In order to discuss the appearance of the new wave clearly, a wave source which produces a simple spherical torsional wave is considered and then the reflection problem by the moving surface is discussed. Applying the regular integral
E-mail address: [email protected]. http://dx.doi.org/10.1016/j.wavemoti.2014.07.004 0165-2125/© 2014 Elsevier B.V. All rights reserved.
1210
K. Watanabe / Wave Motion 51 (2014) 1209–1224
Fig. 2.1. A needle-like torque source in a semi-infinite solid with a moving boundary.
transform and the Cagniard–de Hoop inversion technique, we have obtained an exact closed form solution for the reflected wave, i.e. the transient scattering Doppler effects. Thus, the total wave field is expressed by the simple algebraic functions. Wave propagation behaviors are discussed based on the closed form solution and the appearance of the new wave is reconfirmed for the 3D wave problem. It is shown that the new wave is produced by an incident ray with the critical incident angle. Wave front singularity is also discussed and a numerical example for the surface response is presented. Before starting the analysis, we would like to mention about our substantial assumption for the moving boundary problem: whenever an external force/load is applied to a solid medium, dynamic deformation is induced and its boundary moves slightly. If we input an additional wave to the deformed or deforming medium, some interactions between the initial deformation and the wave field would take place. In this exact sense, the linear elasticity theory is not applicable to the moving boundary problem of elastic media. However, if the amplitude of the incident wave is small enough and the propagation velocity of the incident wave is much larger than that of the boundary motion, the interaction may be neglected. Upon this assumption, we consider the reflection problem of the torsional wave by the moving flat boundary. 2. Solution method Let us take the cylindrical coordinate system (r , θ , z ) fixed in an elastic half space. The surface/boundary of the half space is z = l at an initial time t = 0 and then moves outward with the uniform subsonic velocity V . Thus, the half space occupies the region −∞ < z ≤ l + Vt (see Fig. 2.1). The pure torsional deformation is governed by the single differential equation for the torsional displacement uθ , 1 ∂ uθ uθ ∂ 2 uθ 1 ∂ 2 uθ ρ ∂ 2 uθ + − + − = − Bθ , 2 2 2 2 2 ∂r r ∂r r ∂z cs ∂ t µ
(2.1)
where µ, ρ and cs = (µ/ρ)1/2 are rigidity, density and the shear wave velocity. The circumferential body force Bθ is employed as a wave source. The standard wave source [13] for the torsional wave is a ring source with radius a. It is defined by Bθ = B0
δ(r − a) r
δ(z )f (t ),
(2.2)
where δ(.) and f (t ) are Dirac’s delta function and an arbitrary time function, and B0 is the magnitude of the source. This ring source produces the incident wave of torus (donut) shape, not simple sphere and thus the reflected wave is not simple shape. In order to simplify the wave analysis, if we take the zero ring radius a → 0, the resulting source does not produce any wave due to the torsional nature of the deformation. On the other hand, if we employ a moment couple composed of two forces as was used in the seismology, no pure torsional deformation is produced. In order to produce a simple spherical incident wave for the torsional deformation, a needle-like torque source
Bθ = − lim
a→0
T
δ(r − a)
2πµa
r
H (−z )H (t ),
(2.3)
is proposed by the author [14], where T is the magnitude of torque source (torque per unit length) and H (.) is Heaviside’s unit step function. This source produces two waves: one is a cylindrical wave with semi-infinite length (−∞ < z < 0) and the other is a spherical wave centered at the source edge (z = 0). The cylindrical wave disturbs the negative z region and does not interact with the moving boundary. However, the spherical wave radiates to all directions and catch up to the moving boundary. Thus, only the spherical wave interacts with the moving boundary and produces the reflected wave which has the simple spherical shape. The needle-like source is preferable source for the simpler wave analysis. Therefore, we employ the needle-like torque as the source of the incident wave and then discuss wave reflection at the moving boundary. A small comment on the semi-infinite length of the source have to be added. Replacing the step function H (−z ) in Eq. (2.3) with the delta function δ(z ),
Bθ = − lim
a→0
T
δ(r − a)
2πµa
r
δ(z )H (t ),
(2.4)
K. Watanabe / Wave Motion 51 (2014) 1209–1224
1211
we have a point torque source at the coordinate origin. If we employ this point torque, we have to differentiate the final solution (displacement). Its differentiation is very complicated, but the wave analysis is not simpler than that of the needlelike source. Consequently, it is the simplest way to employ the needle-like torque source defined by Eq. (2.3). A suddenly applied semi-infinite needle-like torque source (body force) [14] is assumed as the wave source. The source lies on the negative z-axis and its edge coincides with the coordinate origin. The 3D pure torsion problem with the needlelike torque source is governed by the simple wave equation for the torsional displacement
1 ∂ uθ uθ ∂ 2 uθ 1 ∂ 2 uθ T δ(r − a) ∂ 2 uθ + − + − = − lim H (−z )H (t ). a→0 2π µa ∂r2 r ∂r r2 ∂ z2 cs2 ∂ t 2 r
(2.5)
Hooke’s law for the torsional stress components is defined by
σz θ = µ
∂ uθ , ∂z
σr θ = µ
∂ uθ uθ − ∂r r
.
(2.6)
The stress-free condition on the moving boundary/surface, t > 0,
[σz θ ]z =l+Vt = 0;
(2.7)
the quiescent condition at an initial time
∂ uθ = 0; ∂t
uθ =
t = 0,
(2.8)
and the convergence condition at infinity
∂ uθ ∂ uθ = = 0; ∂z ∂r
uθ =
r 2 + z 2 → ∞,
(2.9)
are also employed. Consequently, our elastodynamic problem is constituted with Eqs. (2.5)–(2.9). This problem is decomposed into two parts: one is the incident wave field which is produced by the torque source and the other is the reflected wave produced by the moving boundary. We shall discuss these two wave fields, separately. (2.1) Incident wave The incident wave is given as a particular solution of the governing equation and its nonhomogeneous differential equation is rewritten as
(i) (i) (i) 1 ∂ uθ uθ ∂ 2 u(θi) 1 ∂ 2 uθ T δ(r − a) ∂ 2 u(θi) + − + − = − lim H (−z )H (t ) a→0 2π µa ∂r2 r ∂r r2 ∂ z2 cs2 ∂ t 2 r
(2.10)
(i)
where uθ denotes the particular solution. We apply the triple integral transform. Each integral transform is defined by (a) Hankel transform with respect to the radial distance r, f˜ (ξ ) =
∞
rf (r )J1 (ξ r )dr ⇔ f (r ) = 0
∞
ξ f˜ (ξ )J1 (ξ r )dξ ,
(2.11)
0
where J1 (.) is the Bessel function of the first order. (b) Laplace transform with respect to the time t, ∞
f (s) = ∗
f (t ) exp(−st )dt ⇔ f (t ) =
0
1 2π i
Br.(s)
f ∗ (s) exp(st )ds.
(2.12)
f¯ (ζ ) exp(−iζ z )dζ .
(2.13)
(c) Fourier transform with respect to the axial distance z, f¯ (ζ ) =
∞
f (z ) exp(+iζ z )dz ⇔ f (z ) = −∞
1 2π
∞ −∞
The triple integral transform ∞
t =0
∞ z =−∞
(i) (i) (i) ∂ 2 uθ(i) 1 ∂ uθ uθ ∂ 2 u(θi) 1 ∂ 2 uθ T δ(r − a) r + − 2 + − 2 = − lim H (−z )H (t ) a→0 2π µa ∂r2 r ∂r r ∂ z2 cs ∂ t 2 r r =0 ∞
× J1 (ξ r )dr exp(+iζ z )dz exp(−st )dt ,
(2.14)
is carried out with aids of the integration formulas ∞
f (x)δ(x − a)dx = f (a), 0
∞
H (−z ) exp(+iζ z )dz = δ− (ζ ), −∞
(2.15)
1212
K. Watanabe / Wave Motion 51 (2014) 1209–1224
where δ− (.) is Heisenberg’s delta function [15,16] defined by
δ− (ζ ) = π δ(ζ ) −
i
.
ζ
(2.16)
The limit in the right hand side of Eq. (2.14) is evaluated in the transformed domain, as
lim
T 2π µa
a→0
J1 (ξ a) = −
T 4π µ
ξ,
(2.17)
where we have employed the simple formula,
lim
J1 (x)
1
=
x
x→0
2
.
(2.18)
The triple transform leads Eq. (2.10) to a simple algebraic equation for the transformed displacement and the displacement in the transformed domain is given by ∗(i) u¯˜ θ =
ξ i π δ(ζ ) − . 4π µ s ζ 2 + ξ 2 + (s/cs )2 ζ T
(2.19)
Now, let us apply the Fourier inversion defined by Eq. (2.13)
ξT 1 = 4π µs 2π
∗(i) u˜ θ
∞
−∞
1
ζ 2 + ξ 2 + (s/cs )2
π δ(ζ ) −
i
ζ
exp(−iζ z )dζ .
(2.20)
This inversion integral is easily evaluated with use of the formulas, ∞
−∞
1
∞
π
δ(ζ ) 1 exp(−iζ z )dζ = 2 , 2 2 ζ + ξ + (s/cs ) ξ + (s/cs )2
(2.21)
2
0
1
1
ζ + ξ + (s/cs ) ζ 2
2
2
sin(ζ z )dζ =
sgn(z )
1
2 ξ + (s/cs ) 2
2
1 − exp −|z | ξ 2 + (s/cs )2
,
(2.22)
and then Eq. (2.20) yields ∗(i) u˜ θ
ξ /s ξ /s 2 2 = 2H (−z ) 2 + sgn(z ) 2 exp −|z | ξ + (s/cs ) . 8π µ ξ + (s/cs )2 ξ + (s/cs )2
T
(2.23)
The stress component which acts on the moving boundary is derived as 1
µ
σ˜ z∗(θ i) =
∗(i)
du˜ θ
dz
=−
ξ /s
T 8π µ
ξ 2 + (s/cs )2
exp −|z | ξ 2 + (s/cs )2 ,
(2.24)
and its formal Laplace inversion is given by 1
µ
σ˜ z(θi) = −
ξ 8π µ 2π i T
Br.(s)
1 exp −|z | ξ 2 + (s/cs )2 + st ds. s ξ 2 + (s/cs )2
(2.25)
As the stress-free condition on the moving boundary is discussed in the subsequent section, the above equation is retained in this form of Laplace inversion integral. But, the exact inversion of the displacement of Eq. (2.23) has already been carried out and the closed expressions for the displacement and stress have been obtained by the author [14]. Their final results are cited in Appendix A. (2.2) Reflected wave The reflected wave appears after the time at which the incident wave catches up to the moving boundary, and is given (r ) as the solution of the homogeneous differential equation for the displacement uθ , (r ) (r ) (r ) ∂ 2 u(θr ) 1 ∂ uθ uθ ∂ 2 u(θr ) 1 ∂ 2 uθ + − + − = 0. ∂r2 r ∂r r2 ∂ z2 cs2 ∂ t 2
(2.26)
Since the boundary is moving uniformly, we introduce the moving coordinate (r , Z , t ) along the reverse axial direction as r = r,
Z = l + Vt − z ,
t = t.
(2.27)
The differentials in Eq. (2.26) are replaced with new ones,
∂ ∂ = , ∂r ∂r
∂ ∂ =− , ∂z ∂Z
∂ ∂ ∂ =V + ∂t ∂Z ∂t
(2.28)
K. Watanabe / Wave Motion 51 (2014) 1209–1224
1213
and the differential equation (2.26) is rewritten as (r ) (r ) (r ) (r ) 2 (r ) ∂ 2 u(θr ) 1 ∂ uθ uθ M ∂ 2 uθ 1 ∂ 2 uθ 2 ∂ uθ + − + m − 2 − = 0, ∂r2 r ∂r r2 ∂Z2 cs ∂ Z ∂ t cs2 ∂ t 2
(2.29)
where Mach number is introduced as M = V /cs ,
m=
1 − M 2.
(2.30)
It should be noticed that the motion of the boundary is subsonic, M < 1. In order to obtain the general solution, we apply the double integral transform to Eq. (2.29). The double integral transform is composed of Hankel and Laplace transforms. The Hankel transform with respect to the radial distance is the same as that for the incident wave of Eq. (2.10), since the radial distance r is unchanged in spite of the moving coordinate. On the other hand, the Laplace transform with respect to the time t which is defined by ∞
f (p) = ⊗
f (t ) exp(−pt )dt ⇔ f (t ) = 0
1 2π i
f ⊗ (p) exp(pt )dp,
Br.(p)
(2.31)
is different from that for the incident wave since the new axial distance Z includes the time. Therefore, we employ another mark ‘‘⊗’’ for the present Laplace transform. Then, Eq. (2.29) is transformed to the simple ordinary differential equation for ⊗(r ) the transformed displacement u˜ θ , ⊗(r )
d2 u˜ θ
⊗(r ) u˜ θ
− 2M (p/cs )
⊗(r )
du˜ θ
⊗(r )
− {ξ 2 + (p/cs )2 }˜uθ = 0, dZ 2 dZ and its solution which satisfies the convergence condition is given by m2
(2.32)
(mξ )2 + (p/cs )2 − M (p/cs ) = A(ξ , p) exp − Z , 2
(2.33)
m
where A(ξ , p) is the unknown coefficient to be determined by the stress-free condition on the moving boundary. The associated stress is also obtained as 1
µ
r) σ˜ z⊗( θ
= A(ξ , p)
(mξ )2 + (p/cs )2 − M (p/cs )
exp −
m2
(mξ )2 + (p/cs )2 − M (p/cs ) m2
Z
.
(2.34)
Further, the formal Laplace inversion with respect to the parameter p is given by 1
µ
σ˜ z(θr )
=
1 2π i
Br.(p)
A(ξ , p)
(mξ )2 + (p/cs )2 − M (p/cs ) m2
(mξ )2 + (p/cs )2 − M (p/cs ) exp − Z + pt dp. (2.35) 2 m
(2.3) Stress-free condition on the moving boundary The total stress in the Hankel transformed domain is the sum of two stresses, incident and reflected waves,
σ˜ z θ = σ˜ z(θi) + σ˜ z(θr ) .
(2.36)
It is the sum of two different Laplace inversion integrals as 1
µ
σ˜ z θ = − +
ξ 8π µ 2π i T
Br.(s)
1
2π i
1
Br.(p)
s ξ 2 + (s/cs )2
A(ξ , p)
exp −|z | ξ 2 + (s/cs )2 + st ds
(mξ )2 + (p/cs )2 − M (p/cs ) m2
exp −
(mξ )2 + (p/cs )2 − M (p/cs ) m2
Z + pt
dp. (2.37)
We apply the stress-free condition on the moving surface,
1
µ
σ˜ z θ
z =l+Vt Z =0
= 0,
(2.38)
to Eq. (2.37) and have 1 2π i
Br.(p)
A(ξ , p)
ξ = 8π µ 2π i T
(mξ )2 + (p/cs )2 − M (p/cs ) m2
Br.(s)
exp −l ξ 2 + (s/cs )2
s ξ 2 + (s/cs )2
exp(pt )dp
exp
s−V
ξ 2 + (s/cs )2 t ds.
(2.39)
The integral in the left hand side is just in the form of Laplace inversion integral with respect to the parameter p, but, that in the right hand side is not the form of Laplace inversion integral. The latter integral must be converted to the form of the
1214
K. Watanabe / Wave Motion 51 (2014) 1209–1224
Laplace inversion integral with respect to the parameter p. Fortunately, we have already established a conversion formula between two Laplace inversion integrals [11,12,16] as shown in Appendix B. That is
1
exp −l ξ 2 + (s/cs )2
2π i
=
exp
s ξ 2 + (s/cs )2
Br.(s)
1
2π i
Br.(p)
s−V
exp −l ξ 2 + (s/cs )2
ξ 2 + (s/cs )2 t ds
ds
s ξ 2 + (s/cs )2
dp
exp(pt )dp,
(2.40)
s=s(p)
where
ξ 2 + (s/cs )2 , (p/cs ) + M (p/cs )2 + (mξ )2
p=s−V
(2.41)
. (2.42) m2 Applying the above formulas, Eqs. (2.40)–(2.42), to the right hand side of Eq. (2.39), the Laplace inversion integral with respect to the parameter s is converted to that with respect to the parameter p, s(p) = cs
1
2π i
1
Br.(s)
s ξ 2 + (s/cs )2
exp −l ξ 2 + (s/cs )2 exp
exp −
M (p/cs )+
√
(p/cs )2 +(mξ )2
m2
s−V
ξ 2 + (s/cs )2 t ds
l
m2 1 exp(pt )dp. 2π i Br.(p) cs (p/cs ) + M (p/cs )2 + (mξ )2 (p/cs )2 + (mξ )2 1
=
(2.43)
Then, the stress-free condition of Eq. (2.39) is rewritten in the identical form of the Laplace inversion integral with respect to the parameter p, 1
2π i
(mξ )2 + (p/cs )2 − M (p/cs ) m2
Br.(p)
A(ξ , p) exp(pt )dp =
T
ξ
8π µ 2π i
√ M (p/cs )+ (mξ )2 +(p/cs )2 exp − l 2 m m2 1 exp(pt )dp. × 2 2 2 2 cs (p/cs ) + M (mξ ) + (p/cs ) (mξ ) + (p/cs ) Br.(p)
(2.44)
This is just the form of the Laplace inversion integral with respect to the parameter ‘‘p’’ and thus we can determine the unknown function A(ξ , p) as A(ξ , p) =
T
m4
1
1
8π µ cs (p/cs ) + M (mξ ) + (p/cs ) (mξ )2 + (p/cs )2 − M (p/cs ) ξ M (p/cs ) + (mξ )2 + (p/cs )2 exp − × l . m2 (mξ )2 + (p/cs )2
2
2
(2.45)
After the slight modification in the denominator, we have the simpler expression as A(ξ , p) =
ξ 2 M ξ + (p/cs ) (mξ ) + (p/cs )2 (mξ )2 + (p/cs )2 M (p/cs ) + (mξ )2 + (p/cs )2 × exp − l . 2 T
m2
8π µ cs
2
m
⊗(r )
Substituting the above A(ξ , p) into the transformed reflected wave u˜ θ the reflected wave in the transformed domain, ⊗(r )
u˜ θ
=
in Eq. (2.33), we have the explicit expression for
ξ 2 M ξ + (p/cs ) (mξ ) + (p/cs )2 (mξ )2 + (p/cs )2 (Z + l) (mξ )2 + (p/cs )2 − M (p/cs )(Z − l) . × exp − 2 T
(2.46)
m2
8π µ cs
2
m
(2.47)
K. Watanabe / Wave Motion 51 (2014) 1209–1224
1215
3. Cagniard’s technique This section considers the inversion of the reflected wave given by Eq. (2.47). Our inversion technique is the Cagniard– de Hoop method [17] with Gakenheimer’ procedure [18]. As the first step, we apply the Hankel inversion integral to Eq. (2.47), ⊗(r ) uθ
=
ξ2 8π µ cs 0 M ξ 2 + (p/cs ) (p/cs )2 + (mξ )2 (p/cs )2 + (mξ )2 (Z + l) (mξ )2 + (p/cs )2 − M (p/cs )(Z − l) × exp − J1 (ξ r )dξ , 2 m2
T
∞
(3.1)
m
and introduce the variable transform from ξ to η as
ξ = (p/cs )η,
dξ = (p/cs )dη.
(3.2)
Eq. (3.1) is somewhat simplified as ⊗(r )
uθ
=
m2
T
8π µ cs
η2
∞
0
(mη)2 + 1 (Z + l) (mη)2 + 1 − M (Z − l)
M η2 +
× exp −(p/cs )
(mη)2 + 1
J1 ((p/cs )ηr )dη.
m2
(3.3)
Following the Gakenheimer’s procedure [18], we replace the Bessel function with its integral representation [19, p. 178], J1 ((p/cs )ηr ) = −
i
π
π
exp{+i(p/cs )ηr cos φ} cos φ dφ,
(3.4)
0
which is one modification of Hansen’s formula [20, p. 20]. Eq. (3.3) is transformed to the double integral, ⊗(r ) uθ
=
π
(−iη cos φ) 1 2 2 2 − 2 η + m−2 0 0 (M /m)η + η + m (Z + l) η2 + m−2 − (M /m)(Z − l) − i(mr )η cos φ × exp −(p/cs ) ηdφ dη. 1
T
∞
8π 2 µ cs
(3.5)
m
Since the supporting region for the double integral is the upper half plane in the polar coordinate system (η, φ), we can introduce rectilinear coordinates (u, v) in the plane (see Fig. 3.1) as u = η cos φ,
v = η sin φ,
dudv = ηdφ dη,
(3.6)
and change the double integral to that with these new variables as ⊗(r )
uθ
=
(−iu) 1 √ √ 2 + v 2 + m −2 2 2 − 2 u (M /m)( + v ) + u + v + m √ (Z + l) u2 + v 2 + m−2 − (M /m)(Z − l) − i(mr )u × exp −(p/cs ) dudv. 1
T
8π µ cs 2
∞
v=0
∞
u=−∞
u2
2
(3.7)
m
The Laplace transform parameter p is included only in the argument of the exponential function. So, we consider to transform the inner integral with respect to the variable u to the Laplace transform integral with respect to the time, i.e. the Cagniard–de Hoop technique. The inner integral with respect to the variable u is extracted as I (v; p) =
∞
F (u, v) exp −(p/cs )
√ (Z + l) u2 + v 2 + m−2 − (M /m)(Z − l) − i(mr )u m
u=−∞
du,
(3.8)
where F (u, v) =
(M /m)(
u2
(−iu) 1 . √ √ 2 2 − 2 2 +v )+ u +v +m u + v 2 + m −2 2
(3.9)
Now, let us consider the Cagniard path. We introduce the variable transform from u to the time t as t =
√ (Z + l) u2 + v 2 + m−2 − (M /m)(Z − l) − i(mr )u mcs
.
(3.10)
1216
K. Watanabe / Wave Motion 51 (2014) 1209–1224
Fig. 3.1. The supporting half plane for the double integral.
Fig. 3.2. The Cagniard path and the closed loop for the complex integral Φ .
Solving Eq. (3.10) for the variable u, we have the Cagniard path
+i(mr ){mcs t + (M /m)(Z − l)} ± (Z + l) {mcs t + (M /m)(Z − l)}2 − {(mr )2 + (Z + l)2 }(v 2 + m−2 ) u(±) = . (3.11) (mr )2 + (Z + l)2 This Cagniard path is composed of two semi-hyperbolic curves u(±) in the u-plane. As the time t varies from TL (v), {(mr )2 + (Z + l)2 }(v 2 + m−2 ) − (M /m)(Z − l) , (3.12) TL (v) = mcs
to infinity, u(+) moves along the semi-hyperbolic curve in the right hand side from its saddle point usaddle =
+i(mr )
(mr ) + (Z + l) 2
2
v 2 + m−2 ,
(3.13)
and u(−) does along that in the left hand side. Two semi-hyperbolas are connected at the saddle point usaddle . These are shown in Fig. 3.2. The gradient of Cagniard path is given by du(±) dt
= ±mcs
(Z + l){mcs t + (M /m)(Z − l)} ± i(mr ) {mcs t + (M /m)(Z − l)}2 − {(mr )2 + (Z + l)2 }(v 2 + m−2 ) . (3.14) (mr )2 + (Z + l)2 {mcs t + (M /m)(Z − l)}2 − {(mr )2 + (Z + l)2 }(v 2 + m−2 )
The Cagniard path has just been confirmed. We next consider the complex integral whose integrand is the same as that of the inner integral,
Φ=
F (u, v) exp −(p/cs )
√ (Z + l) u2 + v 2 + m−2 − (M /m)(Z − l) − i(mr )u m
L
du.
(3.15)
Examining the singular points in the integrand, we find the branch points
uB = ±i v 2 + m−2 , which are the zeroth of the radical,
uP = ±i v 2 + 1,
(3.16)
√
u2
+v + 2
m −2
= 0, and the simple pole (3.17)
K. Watanabe / Wave Motion 51 (2014) 1209–1224
1217
√ 2 2 −2 which is the zeroth of the denominator, (M /m)(u2 + v 2 ) + √ u + v + m = 0 in F (u, v). Then, two branch cuts along the imaginary axis are introduced so that the conditional Re u2 + v 2 + m−2 > 0 can be satisfied. For the complex integral Φ , the closed loop L shown in Fig. 3.2 is employed. Examining the denominator of F (u, v), we learn that the pole is within the loop only when |usaddle | > |uP |. This inequality gives two conditionals. They are 2 Mr − 1 and Mr > Z + l. (3.18) 0
√
> 0 and Im(u) > 0, and then the line integral along the real axis, i.e. the integral I (v; p), is converted √ to the sum of the integral along the Cagniard path and the residue at the pole uP = +i v 2 + 1. That is ∞ du(+) du(−) F (u(+) , v) I (v; p) = − F (u(−) , v) exp(−pt )dt Re
u2 + v 2 + m−2
dt
TL (v)
+H
Mr
dt
−1 H
Z +l
2
Mr
− 1 − v 2π i Res(uP ).
Z +l
(3.19)
The residue is evaluated exactly as
2π i Res(uP ) = 2π i
u − i v 2 + 1 F (u, v)
× exp −(p/cs ) 2π m
=
2
√ (Z + l) u2 + v 2 + m−2 − (M /m)(Z − l) − i(mr )u
1 + M2
√
m
u=+i
v 2 +1
2lM 2 exp −(p/cs ) +r v +1 . 2
(3.20)
m
Consequently, we have for the inner integral I (v; p) =
∞
TL (v)
+H
F (u(+) , v)
du(+) dt
Mr Z +l
− F (u(−) , v)
−1 H
dt
exp(−pt )dt
2π m2 2lM 2 −1−v exp −(p/cs ) +r v +1 . 1 + M2 m2
2
Mr
du(−)
Z +l
(3.21)
Substituting the above equation into the double integral of Eq. (3.7), we have for the Laplace transformed displacement as ⊗(r ) uθ
=
+H
8π 2 µcs
∞
T
∞
TL (v)
0
F (u(+) , v)
du(+)
− F (u(−) , v)
dt
du(−)
dt
exp(−pt )dt
2 Mr 2π m2 2lM −1 H − 1 − v exp −(p/cs ) + r v 2 + 1 dv. Z +l Z +l 1 + M2 m2
Mr
(3.22)
The order of integration is exchanged for the first double integral whose supporting region is shown in Fig. 3.3. As for the second term, the variable transform from v to the time t
√
t =
2lM /m + (mr ) v 2 + 1 cs m
,
v = vP (t ) =
cs t − 2lM /m2
2
r
− 1,
(3.23)
is introduced. Then, Eq. (3.22) yields to
⊗(r )
uθ
=
T 8π 2 µcs
+H
∞
exp(−pt )dt
tL
Mr Z +l
vL (t )
F (u(+) , v)
0
−1
2π m2 1 + M2
dt
− F (u(−) , v)
cs cs t − 2lM /m2
tQ tP
du(+)
r
cs t − 2lM /m2
2
du(−)
dt
dv
exp(−pt )dt ,
− r2
(3.24)
1218
K. Watanabe / Wave Motion 51 (2014) 1209–1224
where tL = TL (0) =
1
(mr )2 + (Z + l)2 − M (Z − l) ,
m2 cs
(3.25)
{(mcs t ) + (M /m)(Z − l)}2 1 − 2, (mr )2 + (Z + l)2 m 1 2Ml Mr 1 2Ml +r , tQ = + r . tP = cs m2 cs m2 Z +l vL (t ) =
(3.26)
(3.27)
Eq. (3.24) is the Laplace transformed displacement, and it is just in the form of the Laplace transform integral with respect to the parameter p. Thus, we can apply the simple Laplace inversion formulas, L
−1
∞
f (t ) exp(−pt )dt
= H (t − a)f (t ),
a
L
−1
b
f (t ) exp(−pt )dt
(3.28)
= H (t − a)H (b − t )f (t ),
a
and can perform the Laplace inversion by inspection. It yields
T
(r )
uθ (r , Z , t ) =
H (t − tL )
8π 2 µcs
F (u(+) , v)
du(+)
0
Mr
+ H
vL (t )
Z +l
− 1 H (t − tP )H (tQ − t )
− F (u(−) , v)
dt 2π m
dt
dv
cs cs t − 2lM /m2
2
1 + M2
du(−)
r
cs t − 2lM /m2
2
.
(3.29)
− r2
In order to discuss the wave front behavior, the torsional displacement of the reflected wave has to be simplified. Using the relation/identity, 1
du(±)
u2(±) + v 2 + m−2
dt
= ±
mcs
{mcs t + (M /m)(Z − l)} − {(mr )2 + (Z + l)2 }(v 2 + m−2 ) 2
,
(3.30)
the integrand in Eq. (3.29) is rewritten, and then we introduce the variable transform from v to ϕ ,
v = vL (t ) sin ϕ,
(3.31)
for the integral. We have the simpler expression for the integral as vL (t )
F (u(+) , v)
du(+) dt
0
− F (u(−) , v)
du(−)
dt
dv =
(mr ) + (Z + l) 2
π /2
mcs 2
G(u(+) , v) + G(u(−) , v) dϕ,
(3.32)
0
where
G(u, v) = F (u, v) u2 + v 2 + m−2
= (M /m)
√
u2
+v + 2
m −2
(−iu) √ . − M /m u2 + v 2 + m−2 + 1/(Mm)
(3.33)
The variable of Cagniard path u(±) is also simplified as
u(±)
+i(mr ) vL2 (t ) + m−2 ± (Z + l)vL (t ) cos ϕ = . (mr )2 + (Z + l)2
(3.34)
The second term in Eq. (3.29) is also reduced to cs cs t − 2lM /m2
2π m2 1+
M2
r
cs t − 2lM /m2
2
= − r2
2π m2 1+
M2
cs (cs t − cs tP + r ) . √ r (cs t − cs tP + 2r ) (cs t − cs tP )
(3.35)
K. Watanabe / Wave Motion 51 (2014) 1209–1224
1219
Fig. 3.3. The supporting region for the double integral in Eq. (3.22).
Then, Eq. (3.29) yields to the simpler form (r ) uθ ( r , Z , t )
=
T
H (t − tL )
8π µ 2
+H
Mr Z +l
π /2
m
(mr ) + (Z + l) 2
2
− 1 H (t − tP )H (tQ − t )
G(u(+) , v) + G(u(−) , v) dϕ
2π m2
cs (t − tP ) + r
0
1 + M2 r
. √ {cs (t − tP ) + 2r } cs (t − tP )
(3.36)
Fortunately, we can evaluate the integral in the above equation. Its evaluation procedure is described in Appendix C, since it takes a little bit labor time. The final result for the reflected wave is thus given by 8π µ (r ) 1 − M2 m uθ ( r , Z , t ) = H ( t − t L ) T 1 + M 2 Rm
{vm (t ) cos ψ − (M /m)} vm (t ) sin ψ + |vm (t ) cos ψ − (M /m)| cot ψ {vm (t ) cos ψ − (M /m)}2 + vL2 (t ) sin2 ψ vm (t ) sin ψ + {vm (t ) cos ψ + 1/(Mm)} cot ψ − {vm (t ) cos ψ + 1/(Mm)}2 + vL2 (t ) sin2 ψ
sgn ×
+ H (Mr − (Z + l)) H (t − tP )H (tQ − t )
1 − M2 1+
M2
2 {cs (t − tP ) + r } , √ r {cs (t − tP ) + 2r } cs (t − tP )
(3.37)
where the new variables (Rm , ψ) and the function vm (t ) are defined by Eqs. (C.3) and (C.6). 4. Wave front singularity and discussions Now, we shall discuss wave propagation behaviors based on the final closed form solution. The incident wave is produced by the needle-like source and its exact expression is given by Eq. (A.1) in Appendix A. Inspecting Eq. (A.1), the displacement is composed of two waves. One is the cylindrical wave with the product of two step functions, H (−z )H (t − r /cs ). This wave disturbs the negative z-region only and does not interact with the moving boundary as far as we assume the outward motion of the boundary, M = V /cs√> 0. The other is the spherical wave with H (t − R/cs ), where R is the spherical distance from the coordinate origin, R = r 2 + z 2 . This spherical wave radiates from the source edge to all directions and catches up to the moving boundary at the time t =
l cs (1 − M )
.
(4.1)
We call this wave ‘‘incident wave’’ for our moving boundary problem. Clearly, this incident wave shows the finite jump at the front
(i)
uθ
t →R/cs
=
T
r
8π µ zR
.
(4.2)
The wave reflection starts from the time of Eq. (4.1) and the reflected wave is given by the first term in Eq. (3.37). Its disturbed region is determined by the operation of the step function H (t − tL ) and its front t = tL is the spherical shape defined by
r 2 + Z − M (cs t ) −
1 + M2 m2
2 l
2 2Ml = (cs t ) − 2 , m
(4.3)
1220
K. Watanabe / Wave Motion 51 (2014) 1209–1224
in the moving coordinate system (r , Z , t ), or
2
r +
z−
2
2l
= cs t −
m2
2Ml
2
m2
,
(4.4)
in the original fixed coordinate system (r , z , t ). It is interesting that the center of the reflected spherical wave is fixed at z = 2l/m2 in spite of the reflector motion. Examining the first term in Eq. (3.37) near the wave front, we get the finite jump at the front,
(r )
uθ ( r , Z , t )
=
t → tL
T
1 − M 2 m sin ψ
8π µ 1 + M 2
1 cos ψ − M
Rm
−
1 cos ψ + 1/M
.
(4.5)
This jump becomes infinite when cos ψ − M = 0. This condition is equivalent to Z + l = Mr and thus the infinite jump takes place at the apex A of the reflected wave in Fig. 4.2. The second term in Eq. (3.37) shows a new wave whose disturbed region is determined by the product of three step functions,
H
Mr Z +l
− 1 H (t − tP )H (tQ − t ).
(4.6)
Examining the argument of each step function, we learn that the step function H (t − tP ) shows the wave nature and its front t = tP is the cylindrical shape r+
2Ml m2
= cs t .
(4.7)
Other step functions do not show the wave nature, but restrict the disturbed region of the new wave. The region by the product of two step functions, H (t − tP )H (tQ − t ) is included in the region by H (Mr /(Z + l) − 1). Thus the conditional Mr /(Z + l) > 1 does not give any meaning for the wave field. The substantial conditional for the new wave is the product of the former two step functions. When we considered the 2D anti-plane deformation [11,12], the conditional H (tQ − t ) did not appear and only the product H (Mr /(Z + l) − 1) H (t − tP ) determined the disturbed region. The conditional H (tQ − t ) is the characteristic of the 3D wave. The similar non wave front conditional can be found in the case of the 3D elastodynamic problem [18,21]. The new wave appears and starts when the front of the incident wave crosses the apex of the reflected wave just on the moving boundary (see Fig. 4.1(d)). It is the time t =
1 + M2 1−
M2
l cs M
,
(4.8)
when the moving boundary plane (flat vertical line) passes through the center of the spherical reflected wave as shown in Fig. 4.1(d). At this instant, the incident ray (the dotted thin arrow in Fig. 4.1(d)) has the critical incident angle sin−1
1 − M2
1 + M2
,
(4.9)
and its reflected ray is parallel to the plane of the moving boundary. This critical reflection has been shown for the timeharmonic wave [10]. Thus, the new cylindrical wave is produced by the incident ray with the critical incident angle. After the starting time of the new wave, the cylindrical front moves radially outward direction along the moving boundary. Figs. 4.1(e) and 4.2 show this new wave and the disturbed region is a torus with triangular-like cross section which is formed by the product of two step functions, H (t − tP )H (tQ − t ), and the flat moving boundary as shown by the region ABC in Fig. 4.2. In the figure the radial distances, rP , rL and rQ , on the moving boundary are defined by rP = cs t −
2Ml m2
,
rL =
m2
(cs t ) − 2Ml(cs t ) − , 2
l2
rQ =
l M
cs t −
2lM m2
.
(4.10)
If the boundary is immovable, the product of the three step functions does not give any region and thus this cylindrical wave disappears. The wave front behavior is easily extracted from the second term of Eq. (3.37). That is the inverse square root singularity, (r ) uθ ( r , Z , t )
t → tP
=
T
1 − M2
8π µ 1 + M 2
2 r
√
1
cs (t − tP )
.
(4.11)
It will be worth to mention that the new wave has the stronger singularity at its front than that of the regular reflected wave which has the finite jump. The similar wave nature is also found in the case of the 2D SH-wave [11,12]. For visual understanding, the wave front development is shown in Fig. 4.1 and the details around the new cylindrical wave is in Fig. 4.2, schematically.
K. Watanabe / Wave Motion 51 (2014) 1209–1224
(a) 0 < cs t < 1−l M .
M2 (c) 1−l M < cs t < 11+ −M 2
1221
(b) cs t = 1−l M .
l M
+M 2 (d) cs t = 11− M2
.
2
M (e) cs t > 11+ −M 2
l M
l M
.
.
Fig. 4.1. Schematic time development of wave front shape.
Finally, a numerical calculation is carried out for the surface response. The total surface response is given by (i)
(r )
uθ = uθ (r , z , t )|z =Mcs t +l + uθ (r , Z , t )|Z =0 .
(4.12)
Fig. 4.3 shows the response of the torsional displacement at several observing points on the moving surface. When the observing point (radius) is smaller than r /l = 1/M, the first arrival is the edge of the incident and regular reflected waves and the new cylindrical wave is not observed. The edge is at r = rL in Fig. 4.2. On the other hand, when the observing point is larger than r /l = 1/M, the first arrival is the new cylindrical wave and the second is the edge of the incident and reflected waves. The arrival of the wave edge is denoted by the vertical arrow and that of the new cylindrical wave is marked by the
1222
K. Watanabe / Wave Motion 51 (2014) 1209–1224
Fig. 4.2. Details around the new cylindrical wave H (t − tP )H (tQ − t ).
Fig. 4.3. Response of the torsional displacement at the moving surface (arrow: arrival of the wave edge tL , ‘‘∞’’: arrival of tP ).
infinity ‘‘∞’’ in the figure. At the critical observing point, r /l = 1/M, the two arrivals are at one time. We conclude that the new cylindrical wave cannot be observed, when the surface observing point is smaller than the critical radius r /l = 1/M. 5. Conclusion The reflection of the torsional wave by a uniformly moving flat boundary is discussed. An exact closed form solution for the reflected wave is obtained. Discussing wave propagation phenomena, a new cylindrical wave is found in addition to the regular reflected wave. It is produced by the critical ray in the incident wave and runs along the moving boundary after reflection. The new wave has the cylindrical front, but the disturbed region is one of torus with triangular-like cross section. The wave front singularity of the new wave is stronger than those of the incident and the regular reflected waves. The appearance of the new cylindrical wave is one of characteristic wave phenomena for the elastodynamic moving boundary problem. However, the energy break into the regular reflected and the new waves is not discussed here. So, the problem of this energy decomposition is open. Acknowledgment The author would like to express his thanks to the reviewer for many invaluable comments and suggestions which are included in the paper.
K. Watanabe / Wave Motion 51 (2014) 1209–1224
1223
Appendix A. Exact expressions for the incident wave [15]
(i) uθ
T
=
H (−z )H (t − r /cs )
8π µr
σz(θi) = −
T
r
8π R3
σr(θi) = −
2(cs t )
(cs t )2 − r 2
+ H (t − R/cs ) sgn(z )
(cs t ) (cs t )2 − r 2
−
z R
,
H (t − R/cs ),
(A.2)
T
2(cs t )
H (−z )H (t − r /cs )
8π r 2
(A.1)
2−
r2
(cs t )2 − r 2 (cs t )2 − r 2 r2 (cs t ) z z2 2− + H (t − R/cs ) sgn(z ) − 2+ 2 , (cs t )2 − r 2 R R (cs t )2 − r 2
(A.3)
where R=
r 2 + z2.
(A.4)
The displacement of the incident wave shows a finite jump on the spherical wave front as
(i)
uθ
t →R/cs
=
r
T
8π µ zR
.
(A.5)
Appendix B. Conversion of the Laplace inversion integral The conversion formula is established by the author [11,12] and its mathematical procedure is explained concisely in [16]. That is 1 2π i
Br.(s)
f (s) exp
1 ds 2 2 s − V ξ + (s/cs ) t ds = f (s) exp(pt )dp, 2π i Br.(p) dp s=s(p)
(B.1)
where the integrand f (s) vanishes at infinity, i.e. f (s) = 0; |s| → ∞ and the variable transform is defined by
s(p) =
ξ 2 + (s/cs )2 , p + cs M (mξ )2 + (p/cs )2
p=s−V
m2
(B.2)
.
(B.3)
Two Bromwich paths are Br.(s) : from γs − i∞ to γs − i∞, Br.(p) : from γp − i∞ to γp − i∞, where γs and γp are sufficiently large constants and have the relation
γp = (1 − M )γs .
(B.4)
Appendix C. Evaluation of the integral in Eq. (3.36) We consider the integral, I (r , Z , t ) =
π/2
m
(mr ) + (Z + l) 2
2
G(u(+) , v) + G(u(−) , v) dϕ.
(C.1)
0
Using partial fractions, we decompose the integrand as G(u, v) =
1 − M2 1 + M2
√
(−iu)
u2 + v 2 + m−2 − M /m
−√
(−iu)
u2 + v 2 + m−2 + 1/(Mm)
,
(C.2)
and introduce the polar coordinate-like variables (Rm , ψ) as mr = Rm sin ψ,
Z + l = Rm cos ψ;
Rm =
(mr )2 + (Z + l)2 .
(C.3)
1224
K. Watanabe / Wave Motion 51 (2014) 1209–1224
The variable −iu(±) and radicals are simplified as
−iu(±) = vm (t ) sin ψ ∓ ivL (t ) cos ψ cos ϕ, u2(±) + v 2 + m−2 = vm (t ) cos ψ ± ivL (t ) sin ψ cos ϕ,
(C.4) (C.5)
where
vm (t ) =
vL2 (t ) + m−2 .
(C.6)
Since
vm (t ) sin ψ − ivL (t ) cos ψ cos ϕ 1 + M2 vm (t ) cos ψ − (M /m) + ivL (t ) sin ψ cos ϕ vm (t ) sin ψ − ivL (t ) cos ψ cos ϕ , − Re vm (t ) cos ψ + 1/(Mm) + ivL (t ) sin ψ cos ϕ
G(u(+) , v) + G(u(−) , v) = 2
1 − M2
Re
(C.7)
the integral I (r , Z , t ) is rewritten as π/2
{vm (t ) cos ψ − (M /m)} vm (t ) − vL2 (t ) cos ψ cos2 ϕ dϕ 1 + M 2 Rm {vm (t ) cos ψ − (M /m)}2 + vL2 (t ) sin2 ψ cos2 ϕ 0 π/2 {vm (t ) cos ψ + 1/(Mm)} vm (t ) − vL2 (t ) cos ψ cos2 ϕ − d ϕ . {vm (t ) cos ψ + 1/(Mm)}2 + vL2 (t ) sin2 ψ cos2 ϕ 0
I (r , Z , t ) = 2
1 − M 2 m sin ψ
(C.8)
We apply the simple integration formulas, π/2
dϕ
√
π
,
π/2
cos2 ϕ
dϕ =
π
a2 + b2 cos2 ϕ 2b2 2a a2 + b2 0 to the integrals in Eq. (C.8). It results in the closed expression for the integral, 0
a2 + b2 cos2 ϕ
=
a 1− √ a2 + b 2
;
a > 0, b > 0, (C.9)
− (M /m)}vm (t ) sin ψ + |vm (t ) cos ψ − (M /m)| cot ψ sgn{vm (t ) cos ψ {vm (t ) cos ψ − (M /m)}2 vL2 (t ) sin2 ψ vm (t ) sin ψ + {vm (t ) cos ψ + 1/(Mm)} cot ψ − , {vm (t ) cos ψ + 1/(Mm)}2 + vL2 (t ) sin2 ψ
I (r , Z , t ) = π
1 − M2 m
1 + M 2 Rm
(C.10)
where the sign function is defined by sgn(x) =
+1; −1;
x > 0, x < 0.
(C.11)
References [1] D. Censor, Scattering of a plane wave at a plane interface separating two moving media, Radio Sci. 4 (11) (1969) 1079–1088. [2] D. Zutter, Reflection from linearly vibrating objects: plane mirror at oblique incidence, IEEE Trans. Antennas and Propagation AP-30 (5) (1981) 898–903. [3] Yao-Xiong Huang, Reflection and transmission of electromagnetic waves by a dielectric medium moving in an arbitrary direction, J. Appl. Phys. 76 (5) (1994) 2575–2581. [4] D. Censor, Scattering by time varying obstacles, J. Sound Vib. 25 (1) (1972) 101–110. [5] O. Bou Matt, J.P. Remenieras, C. Brunel, A. Roncin, F. Patat, Ultrasonic sensing of vibration, Ultarsonics 36 (1998) 391–396. [6] R. Wunenburger, N. Mujica, S. Fauve, Experimental study of the Doppler shift generated by a vibrating scatterer, J. Acoust. Soc. Am. 115 (2) (2004) 507–514. [7] D. Censor, The generalized Doppler effect and applications, J. Franklin Inst. 295 (2) (1973) 103–116. [8] S.-Y. Lee, M. Lee, A new wave technique for free vibration of a string with time-varying length, Trans. ASME, Ser. E, J. Appl. Mech. 69 (2002) 83–87. [9] H. Kimura, H. Ito, T. Nakagawa, Vibration analysis of elevator rope: forced vibration of rope with time-varying length, Trans. JSME, Ser. C 71 (706) (2005) 1871–1876 (in Japanese). [10] Kazumi Watanabe, Elastodynamic Doppler effects by a moving stress-free edge, Internat. J. Engrg. Sci. 48 (2010) 1995–2014. [11] Kazumi Watanabe, Reflection of an impulsive SH-wave at a moving edge, Quart. J. Mech. Appl. Math. 63 (3) (2010) 335–348. [12] Kazumi Watanabe, Transient SH-waves in an elastic plate with a moving edge, Acta Mech. 223 (8) (2012) 1823–1836. [13] K. Watanabe, R. Payton, Green’s function for torsional waves in a cylindrically monoclinic material, Internat. J. Engrg. Sci. 43 (2005) 1283–1291. [14] K. Watanabe, Torsional waves produced by a semi-infinite needle-like torque source, ZAMM Z. Angew. Math. Mech. (2014) http://dx.doi.org/10.1002/ zamm.201400043. in press. [15] I.N. Sneddon, Fourier Transforms, McGraw-Hill, London, 1951. [16] K. Watanabe, Integral Transform Techniques for Green’s Function, Springer, Switzerland, 2014. [17] A.T. De Hoop, A modification of Cagniard’s method for solving seismic pulse problems, Appl. Sci. Res., Sect. B 8 (1960) 349–356. [18] D.C. Gakenheimer, Response of an elastic half space to expanding surface loads, Trans. ASME, Ser. E, J. Appl. Mech. 38 (1) (1971) 99–110. [19] S. Moriguchi, K. Udagawa, S. Ichimatsu, Mathematical Formulas, Vol. III, Iwanami, Tokyo, 1987 (in Japanese). [20] G.N. Watson, Theory of Bessel Functions, second ed., Cambridge, 1966, p. 20. [21] F.R. Norwood, Exact transient response of an elastic half space loaded over a rectangular region of its surface, Trans. ASME, Ser. E, J. Appl. Mech. 36 (3) (1969) 516–522.
Contents lists available at ScienceDirect
Wave Motion journal homepage: www.elsevier.com/locate/wavemoti
Reflection of an impulsive torsional wave by a moving boundary Kazumi Watanabe Sakuragaoka_1-4, Hikoshima, Shimonoseki, 750-0083, Japan
highlights • • • • •
Reflection problem of an impulsive torsional wave by a moving boundary is discussed exactly. A torsional wave source which produces the simple spherical wave is discussed. Exact closed form solution for the reflected wave is obtained. A new wave is found. Development of wave front shape is analytically obtained and illustrated.
article
info
Article history: Received 30 March 2014 Received in revised form 14 June 2014 Accepted 9 July 2014 Available online 17 July 2014 Keywords: Moving boundary Torsional wave Reflection Closed form solution
abstract Reflection of an impulsive torsional wave by a uniformly moving surface/boundary is discussed. Applying the Cagniard–de Hoop technique, the solution for the reflected wave is obtained in the closed form and a new wave that runs along the moving boundary has been found. Its wave front singularity is stronger than those of the incident and the regular reflected waves. This new wave has the cylindrical surface, but the shape of the disturbed region is a torus with a triangular-like cross section. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Wave reflection at a moving object/boundary is the substantial sensing principle for the radar systems. Numerous works have been carried out for electromagnetic waves (for example, [1–3]). The similar Doppler effects for acoustic waves have also been discussed theoretically and experimentally [4–6]. Vibration and waves in ropes and rods with moving supports are also important engineering problems and thus the vibration problem for the high speed elevator and the power transmission have been discussed [7–9]. However, as for the 2D and 3D elastic waves, not so many works have been done for the reflection problem by the moving boundary. Reflection of the 2D time-harmonic wave at a uniformly moving boundary and that of an impulsive SH-wave have been discussed by the author [10–12]. When the impulsive wave is reflected by the uniformly moving boundary, a new additional wave is produced and it runs along the moving boundary [11,12]. This new wave has been found in the case of the 2D anti-plane deformation. In order to confirm the appearance of the new wave for the 3D dynamic deformation, the reflection of the simplest 3D wave, i.e. the torsional wave, is discussed. This is the motivation of the present work. In order to discuss the appearance of the new wave clearly, a wave source which produces a simple spherical torsional wave is considered and then the reflection problem by the moving surface is discussed. Applying the regular integral
E-mail address: [email protected]. http://dx.doi.org/10.1016/j.wavemoti.2014.07.004 0165-2125/© 2014 Elsevier B.V. All rights reserved.
1210
K. Watanabe / Wave Motion 51 (2014) 1209–1224
Fig. 2.1. A needle-like torque source in a semi-infinite solid with a moving boundary.
transform and the Cagniard–de Hoop inversion technique, we have obtained an exact closed form solution for the reflected wave, i.e. the transient scattering Doppler effects. Thus, the total wave field is expressed by the simple algebraic functions. Wave propagation behaviors are discussed based on the closed form solution and the appearance of the new wave is reconfirmed for the 3D wave problem. It is shown that the new wave is produced by an incident ray with the critical incident angle. Wave front singularity is also discussed and a numerical example for the surface response is presented. Before starting the analysis, we would like to mention about our substantial assumption for the moving boundary problem: whenever an external force/load is applied to a solid medium, dynamic deformation is induced and its boundary moves slightly. If we input an additional wave to the deformed or deforming medium, some interactions between the initial deformation and the wave field would take place. In this exact sense, the linear elasticity theory is not applicable to the moving boundary problem of elastic media. However, if the amplitude of the incident wave is small enough and the propagation velocity of the incident wave is much larger than that of the boundary motion, the interaction may be neglected. Upon this assumption, we consider the reflection problem of the torsional wave by the moving flat boundary. 2. Solution method Let us take the cylindrical coordinate system (r , θ , z ) fixed in an elastic half space. The surface/boundary of the half space is z = l at an initial time t = 0 and then moves outward with the uniform subsonic velocity V . Thus, the half space occupies the region −∞ < z ≤ l + Vt (see Fig. 2.1). The pure torsional deformation is governed by the single differential equation for the torsional displacement uθ , 1 ∂ uθ uθ ∂ 2 uθ 1 ∂ 2 uθ ρ ∂ 2 uθ + − + − = − Bθ , 2 2 2 2 2 ∂r r ∂r r ∂z cs ∂ t µ
(2.1)
where µ, ρ and cs = (µ/ρ)1/2 are rigidity, density and the shear wave velocity. The circumferential body force Bθ is employed as a wave source. The standard wave source [13] for the torsional wave is a ring source with radius a. It is defined by Bθ = B0
δ(r − a) r
δ(z )f (t ),
(2.2)
where δ(.) and f (t ) are Dirac’s delta function and an arbitrary time function, and B0 is the magnitude of the source. This ring source produces the incident wave of torus (donut) shape, not simple sphere and thus the reflected wave is not simple shape. In order to simplify the wave analysis, if we take the zero ring radius a → 0, the resulting source does not produce any wave due to the torsional nature of the deformation. On the other hand, if we employ a moment couple composed of two forces as was used in the seismology, no pure torsional deformation is produced. In order to produce a simple spherical incident wave for the torsional deformation, a needle-like torque source
Bθ = − lim
a→0
T
δ(r − a)
2πµa
r
H (−z )H (t ),
(2.3)
is proposed by the author [14], where T is the magnitude of torque source (torque per unit length) and H (.) is Heaviside’s unit step function. This source produces two waves: one is a cylindrical wave with semi-infinite length (−∞ < z < 0) and the other is a spherical wave centered at the source edge (z = 0). The cylindrical wave disturbs the negative z region and does not interact with the moving boundary. However, the spherical wave radiates to all directions and catch up to the moving boundary. Thus, only the spherical wave interacts with the moving boundary and produces the reflected wave which has the simple spherical shape. The needle-like source is preferable source for the simpler wave analysis. Therefore, we employ the needle-like torque as the source of the incident wave and then discuss wave reflection at the moving boundary. A small comment on the semi-infinite length of the source have to be added. Replacing the step function H (−z ) in Eq. (2.3) with the delta function δ(z ),
Bθ = − lim
a→0
T
δ(r − a)
2πµa
r
δ(z )H (t ),
(2.4)
K. Watanabe / Wave Motion 51 (2014) 1209–1224
1211
we have a point torque source at the coordinate origin. If we employ this point torque, we have to differentiate the final solution (displacement). Its differentiation is very complicated, but the wave analysis is not simpler than that of the needlelike source. Consequently, it is the simplest way to employ the needle-like torque source defined by Eq. (2.3). A suddenly applied semi-infinite needle-like torque source (body force) [14] is assumed as the wave source. The source lies on the negative z-axis and its edge coincides with the coordinate origin. The 3D pure torsion problem with the needlelike torque source is governed by the simple wave equation for the torsional displacement
1 ∂ uθ uθ ∂ 2 uθ 1 ∂ 2 uθ T δ(r − a) ∂ 2 uθ + − + − = − lim H (−z )H (t ). a→0 2π µa ∂r2 r ∂r r2 ∂ z2 cs2 ∂ t 2 r
(2.5)
Hooke’s law for the torsional stress components is defined by
σz θ = µ
∂ uθ , ∂z
σr θ = µ
∂ uθ uθ − ∂r r
.
(2.6)
The stress-free condition on the moving boundary/surface, t > 0,
[σz θ ]z =l+Vt = 0;
(2.7)
the quiescent condition at an initial time
∂ uθ = 0; ∂t
uθ =
t = 0,
(2.8)
and the convergence condition at infinity
∂ uθ ∂ uθ = = 0; ∂z ∂r
uθ =
r 2 + z 2 → ∞,
(2.9)
are also employed. Consequently, our elastodynamic problem is constituted with Eqs. (2.5)–(2.9). This problem is decomposed into two parts: one is the incident wave field which is produced by the torque source and the other is the reflected wave produced by the moving boundary. We shall discuss these two wave fields, separately. (2.1) Incident wave The incident wave is given as a particular solution of the governing equation and its nonhomogeneous differential equation is rewritten as
(i) (i) (i) 1 ∂ uθ uθ ∂ 2 u(θi) 1 ∂ 2 uθ T δ(r − a) ∂ 2 u(θi) + − + − = − lim H (−z )H (t ) a→0 2π µa ∂r2 r ∂r r2 ∂ z2 cs2 ∂ t 2 r
(2.10)
(i)
where uθ denotes the particular solution. We apply the triple integral transform. Each integral transform is defined by (a) Hankel transform with respect to the radial distance r, f˜ (ξ ) =
∞
rf (r )J1 (ξ r )dr ⇔ f (r ) = 0
∞
ξ f˜ (ξ )J1 (ξ r )dξ ,
(2.11)
0
where J1 (.) is the Bessel function of the first order. (b) Laplace transform with respect to the time t, ∞
f (s) = ∗
f (t ) exp(−st )dt ⇔ f (t ) =
0
1 2π i
Br.(s)
f ∗ (s) exp(st )ds.
(2.12)
f¯ (ζ ) exp(−iζ z )dζ .
(2.13)
(c) Fourier transform with respect to the axial distance z, f¯ (ζ ) =
∞
f (z ) exp(+iζ z )dz ⇔ f (z ) = −∞
1 2π
∞ −∞
The triple integral transform ∞
t =0
∞ z =−∞
(i) (i) (i) ∂ 2 uθ(i) 1 ∂ uθ uθ ∂ 2 u(θi) 1 ∂ 2 uθ T δ(r − a) r + − 2 + − 2 = − lim H (−z )H (t ) a→0 2π µa ∂r2 r ∂r r ∂ z2 cs ∂ t 2 r r =0 ∞
× J1 (ξ r )dr exp(+iζ z )dz exp(−st )dt ,
(2.14)
is carried out with aids of the integration formulas ∞
f (x)δ(x − a)dx = f (a), 0
∞
H (−z ) exp(+iζ z )dz = δ− (ζ ), −∞
(2.15)
1212
K. Watanabe / Wave Motion 51 (2014) 1209–1224
where δ− (.) is Heisenberg’s delta function [15,16] defined by
δ− (ζ ) = π δ(ζ ) −
i
.
ζ
(2.16)
The limit in the right hand side of Eq. (2.14) is evaluated in the transformed domain, as
lim
T 2π µa
a→0
J1 (ξ a) = −
T 4π µ
ξ,
(2.17)
where we have employed the simple formula,
lim
J1 (x)
1
=
x
x→0
2
.
(2.18)
The triple transform leads Eq. (2.10) to a simple algebraic equation for the transformed displacement and the displacement in the transformed domain is given by ∗(i) u¯˜ θ =
ξ i π δ(ζ ) − . 4π µ s ζ 2 + ξ 2 + (s/cs )2 ζ T
(2.19)
Now, let us apply the Fourier inversion defined by Eq. (2.13)
ξT 1 = 4π µs 2π
∗(i) u˜ θ
∞
−∞
1
ζ 2 + ξ 2 + (s/cs )2
π δ(ζ ) −
i
ζ
exp(−iζ z )dζ .
(2.20)
This inversion integral is easily evaluated with use of the formulas, ∞
−∞
1
∞
π
δ(ζ ) 1 exp(−iζ z )dζ = 2 , 2 2 ζ + ξ + (s/cs ) ξ + (s/cs )2
(2.21)
2
0
1
1
ζ + ξ + (s/cs ) ζ 2
2
2
sin(ζ z )dζ =
sgn(z )
1
2 ξ + (s/cs ) 2
2
1 − exp −|z | ξ 2 + (s/cs )2
,
(2.22)
and then Eq. (2.20) yields ∗(i) u˜ θ
ξ /s ξ /s 2 2 = 2H (−z ) 2 + sgn(z ) 2 exp −|z | ξ + (s/cs ) . 8π µ ξ + (s/cs )2 ξ + (s/cs )2
T
(2.23)
The stress component which acts on the moving boundary is derived as 1
µ
σ˜ z∗(θ i) =
∗(i)
du˜ θ
dz
=−
ξ /s
T 8π µ
ξ 2 + (s/cs )2
exp −|z | ξ 2 + (s/cs )2 ,
(2.24)
and its formal Laplace inversion is given by 1
µ
σ˜ z(θi) = −
ξ 8π µ 2π i T
Br.(s)
1 exp −|z | ξ 2 + (s/cs )2 + st ds. s ξ 2 + (s/cs )2
(2.25)
As the stress-free condition on the moving boundary is discussed in the subsequent section, the above equation is retained in this form of Laplace inversion integral. But, the exact inversion of the displacement of Eq. (2.23) has already been carried out and the closed expressions for the displacement and stress have been obtained by the author [14]. Their final results are cited in Appendix A. (2.2) Reflected wave The reflected wave appears after the time at which the incident wave catches up to the moving boundary, and is given (r ) as the solution of the homogeneous differential equation for the displacement uθ , (r ) (r ) (r ) ∂ 2 u(θr ) 1 ∂ uθ uθ ∂ 2 u(θr ) 1 ∂ 2 uθ + − + − = 0. ∂r2 r ∂r r2 ∂ z2 cs2 ∂ t 2
(2.26)
Since the boundary is moving uniformly, we introduce the moving coordinate (r , Z , t ) along the reverse axial direction as r = r,
Z = l + Vt − z ,
t = t.
(2.27)
The differentials in Eq. (2.26) are replaced with new ones,
∂ ∂ = , ∂r ∂r
∂ ∂ =− , ∂z ∂Z
∂ ∂ ∂ =V + ∂t ∂Z ∂t
(2.28)
K. Watanabe / Wave Motion 51 (2014) 1209–1224
1213
and the differential equation (2.26) is rewritten as (r ) (r ) (r ) (r ) 2 (r ) ∂ 2 u(θr ) 1 ∂ uθ uθ M ∂ 2 uθ 1 ∂ 2 uθ 2 ∂ uθ + − + m − 2 − = 0, ∂r2 r ∂r r2 ∂Z2 cs ∂ Z ∂ t cs2 ∂ t 2
(2.29)
where Mach number is introduced as M = V /cs ,
m=
1 − M 2.
(2.30)
It should be noticed that the motion of the boundary is subsonic, M < 1. In order to obtain the general solution, we apply the double integral transform to Eq. (2.29). The double integral transform is composed of Hankel and Laplace transforms. The Hankel transform with respect to the radial distance is the same as that for the incident wave of Eq. (2.10), since the radial distance r is unchanged in spite of the moving coordinate. On the other hand, the Laplace transform with respect to the time t which is defined by ∞
f (p) = ⊗
f (t ) exp(−pt )dt ⇔ f (t ) = 0
1 2π i
f ⊗ (p) exp(pt )dp,
Br.(p)
(2.31)
is different from that for the incident wave since the new axial distance Z includes the time. Therefore, we employ another mark ‘‘⊗’’ for the present Laplace transform. Then, Eq. (2.29) is transformed to the simple ordinary differential equation for ⊗(r ) the transformed displacement u˜ θ , ⊗(r )
d2 u˜ θ
⊗(r ) u˜ θ
− 2M (p/cs )
⊗(r )
du˜ θ
⊗(r )
− {ξ 2 + (p/cs )2 }˜uθ = 0, dZ 2 dZ and its solution which satisfies the convergence condition is given by m2
(2.32)
(mξ )2 + (p/cs )2 − M (p/cs ) = A(ξ , p) exp − Z , 2
(2.33)
m
where A(ξ , p) is the unknown coefficient to be determined by the stress-free condition on the moving boundary. The associated stress is also obtained as 1
µ
r) σ˜ z⊗( θ
= A(ξ , p)
(mξ )2 + (p/cs )2 − M (p/cs )
exp −
m2
(mξ )2 + (p/cs )2 − M (p/cs ) m2
Z
.
(2.34)
Further, the formal Laplace inversion with respect to the parameter p is given by 1
µ
σ˜ z(θr )
=
1 2π i
Br.(p)
A(ξ , p)
(mξ )2 + (p/cs )2 − M (p/cs ) m2
(mξ )2 + (p/cs )2 − M (p/cs ) exp − Z + pt dp. (2.35) 2 m
(2.3) Stress-free condition on the moving boundary The total stress in the Hankel transformed domain is the sum of two stresses, incident and reflected waves,
σ˜ z θ = σ˜ z(θi) + σ˜ z(θr ) .
(2.36)
It is the sum of two different Laplace inversion integrals as 1
µ
σ˜ z θ = − +
ξ 8π µ 2π i T
Br.(s)
1
2π i
1
Br.(p)
s ξ 2 + (s/cs )2
A(ξ , p)
exp −|z | ξ 2 + (s/cs )2 + st ds
(mξ )2 + (p/cs )2 − M (p/cs ) m2
exp −
(mξ )2 + (p/cs )2 − M (p/cs ) m2
Z + pt
dp. (2.37)
We apply the stress-free condition on the moving surface,
1
µ
σ˜ z θ
z =l+Vt Z =0
= 0,
(2.38)
to Eq. (2.37) and have 1 2π i
Br.(p)
A(ξ , p)
ξ = 8π µ 2π i T
(mξ )2 + (p/cs )2 − M (p/cs ) m2
Br.(s)
exp −l ξ 2 + (s/cs )2
s ξ 2 + (s/cs )2
exp(pt )dp
exp
s−V
ξ 2 + (s/cs )2 t ds.
(2.39)
The integral in the left hand side is just in the form of Laplace inversion integral with respect to the parameter p, but, that in the right hand side is not the form of Laplace inversion integral. The latter integral must be converted to the form of the
1214
K. Watanabe / Wave Motion 51 (2014) 1209–1224
Laplace inversion integral with respect to the parameter p. Fortunately, we have already established a conversion formula between two Laplace inversion integrals [11,12,16] as shown in Appendix B. That is
1
exp −l ξ 2 + (s/cs )2
2π i
=
exp
s ξ 2 + (s/cs )2
Br.(s)
1
2π i
Br.(p)
s−V
exp −l ξ 2 + (s/cs )2
ξ 2 + (s/cs )2 t ds
ds
s ξ 2 + (s/cs )2
dp
exp(pt )dp,
(2.40)
s=s(p)
where
ξ 2 + (s/cs )2 , (p/cs ) + M (p/cs )2 + (mξ )2
p=s−V
(2.41)
. (2.42) m2 Applying the above formulas, Eqs. (2.40)–(2.42), to the right hand side of Eq. (2.39), the Laplace inversion integral with respect to the parameter s is converted to that with respect to the parameter p, s(p) = cs
1
2π i
1
Br.(s)
s ξ 2 + (s/cs )2
exp −l ξ 2 + (s/cs )2 exp
exp −
M (p/cs )+
√
(p/cs )2 +(mξ )2
m2
s−V
ξ 2 + (s/cs )2 t ds
l
m2 1 exp(pt )dp. 2π i Br.(p) cs (p/cs ) + M (p/cs )2 + (mξ )2 (p/cs )2 + (mξ )2 1
=
(2.43)
Then, the stress-free condition of Eq. (2.39) is rewritten in the identical form of the Laplace inversion integral with respect to the parameter p, 1
2π i
(mξ )2 + (p/cs )2 − M (p/cs ) m2
Br.(p)
A(ξ , p) exp(pt )dp =
T
ξ
8π µ 2π i
√ M (p/cs )+ (mξ )2 +(p/cs )2 exp − l 2 m m2 1 exp(pt )dp. × 2 2 2 2 cs (p/cs ) + M (mξ ) + (p/cs ) (mξ ) + (p/cs ) Br.(p)
(2.44)
This is just the form of the Laplace inversion integral with respect to the parameter ‘‘p’’ and thus we can determine the unknown function A(ξ , p) as A(ξ , p) =
T
m4
1
1
8π µ cs (p/cs ) + M (mξ ) + (p/cs ) (mξ )2 + (p/cs )2 − M (p/cs ) ξ M (p/cs ) + (mξ )2 + (p/cs )2 exp − × l . m2 (mξ )2 + (p/cs )2
2
2
(2.45)
After the slight modification in the denominator, we have the simpler expression as A(ξ , p) =
ξ 2 M ξ + (p/cs ) (mξ ) + (p/cs )2 (mξ )2 + (p/cs )2 M (p/cs ) + (mξ )2 + (p/cs )2 × exp − l . 2 T
m2
8π µ cs
2
m
⊗(r )
Substituting the above A(ξ , p) into the transformed reflected wave u˜ θ the reflected wave in the transformed domain, ⊗(r )
u˜ θ
=
in Eq. (2.33), we have the explicit expression for
ξ 2 M ξ + (p/cs ) (mξ ) + (p/cs )2 (mξ )2 + (p/cs )2 (Z + l) (mξ )2 + (p/cs )2 − M (p/cs )(Z − l) . × exp − 2 T
(2.46)
m2
8π µ cs
2
m
(2.47)
K. Watanabe / Wave Motion 51 (2014) 1209–1224
1215
3. Cagniard’s technique This section considers the inversion of the reflected wave given by Eq. (2.47). Our inversion technique is the Cagniard– de Hoop method [17] with Gakenheimer’ procedure [18]. As the first step, we apply the Hankel inversion integral to Eq. (2.47), ⊗(r ) uθ
=
ξ2 8π µ cs 0 M ξ 2 + (p/cs ) (p/cs )2 + (mξ )2 (p/cs )2 + (mξ )2 (Z + l) (mξ )2 + (p/cs )2 − M (p/cs )(Z − l) × exp − J1 (ξ r )dξ , 2 m2
T
∞
(3.1)
m
and introduce the variable transform from ξ to η as
ξ = (p/cs )η,
dξ = (p/cs )dη.
(3.2)
Eq. (3.1) is somewhat simplified as ⊗(r )
uθ
=
m2
T
8π µ cs
η2
∞
0
(mη)2 + 1 (Z + l) (mη)2 + 1 − M (Z − l)
M η2 +
× exp −(p/cs )
(mη)2 + 1
J1 ((p/cs )ηr )dη.
m2
(3.3)
Following the Gakenheimer’s procedure [18], we replace the Bessel function with its integral representation [19, p. 178], J1 ((p/cs )ηr ) = −
i
π
π
exp{+i(p/cs )ηr cos φ} cos φ dφ,
(3.4)
0
which is one modification of Hansen’s formula [20, p. 20]. Eq. (3.3) is transformed to the double integral, ⊗(r ) uθ
=
π
(−iη cos φ) 1 2 2 2 − 2 η + m−2 0 0 (M /m)η + η + m (Z + l) η2 + m−2 − (M /m)(Z − l) − i(mr )η cos φ × exp −(p/cs ) ηdφ dη. 1
T
∞
8π 2 µ cs
(3.5)
m
Since the supporting region for the double integral is the upper half plane in the polar coordinate system (η, φ), we can introduce rectilinear coordinates (u, v) in the plane (see Fig. 3.1) as u = η cos φ,
v = η sin φ,
dudv = ηdφ dη,
(3.6)
and change the double integral to that with these new variables as ⊗(r )
uθ
=
(−iu) 1 √ √ 2 + v 2 + m −2 2 2 − 2 u (M /m)( + v ) + u + v + m √ (Z + l) u2 + v 2 + m−2 − (M /m)(Z − l) − i(mr )u × exp −(p/cs ) dudv. 1
T
8π µ cs 2
∞
v=0
∞
u=−∞
u2
2
(3.7)
m
The Laplace transform parameter p is included only in the argument of the exponential function. So, we consider to transform the inner integral with respect to the variable u to the Laplace transform integral with respect to the time, i.e. the Cagniard–de Hoop technique. The inner integral with respect to the variable u is extracted as I (v; p) =
∞
F (u, v) exp −(p/cs )
√ (Z + l) u2 + v 2 + m−2 − (M /m)(Z − l) − i(mr )u m
u=−∞
du,
(3.8)
where F (u, v) =
(M /m)(
u2
(−iu) 1 . √ √ 2 2 − 2 2 +v )+ u +v +m u + v 2 + m −2 2
(3.9)
Now, let us consider the Cagniard path. We introduce the variable transform from u to the time t as t =
√ (Z + l) u2 + v 2 + m−2 − (M /m)(Z − l) − i(mr )u mcs
.
(3.10)
1216
K. Watanabe / Wave Motion 51 (2014) 1209–1224
Fig. 3.1. The supporting half plane for the double integral.
Fig. 3.2. The Cagniard path and the closed loop for the complex integral Φ .
Solving Eq. (3.10) for the variable u, we have the Cagniard path
+i(mr ){mcs t + (M /m)(Z − l)} ± (Z + l) {mcs t + (M /m)(Z − l)}2 − {(mr )2 + (Z + l)2 }(v 2 + m−2 ) u(±) = . (3.11) (mr )2 + (Z + l)2 This Cagniard path is composed of two semi-hyperbolic curves u(±) in the u-plane. As the time t varies from TL (v), {(mr )2 + (Z + l)2 }(v 2 + m−2 ) − (M /m)(Z − l) , (3.12) TL (v) = mcs
to infinity, u(+) moves along the semi-hyperbolic curve in the right hand side from its saddle point usaddle =
+i(mr )
(mr ) + (Z + l) 2
2
v 2 + m−2 ,
(3.13)
and u(−) does along that in the left hand side. Two semi-hyperbolas are connected at the saddle point usaddle . These are shown in Fig. 3.2. The gradient of Cagniard path is given by du(±) dt
= ±mcs
(Z + l){mcs t + (M /m)(Z − l)} ± i(mr ) {mcs t + (M /m)(Z − l)}2 − {(mr )2 + (Z + l)2 }(v 2 + m−2 ) . (3.14) (mr )2 + (Z + l)2 {mcs t + (M /m)(Z − l)}2 − {(mr )2 + (Z + l)2 }(v 2 + m−2 )
The Cagniard path has just been confirmed. We next consider the complex integral whose integrand is the same as that of the inner integral,
Φ=
F (u, v) exp −(p/cs )
√ (Z + l) u2 + v 2 + m−2 − (M /m)(Z − l) − i(mr )u m
L
du.
(3.15)
Examining the singular points in the integrand, we find the branch points
uB = ±i v 2 + m−2 , which are the zeroth of the radical,
uP = ±i v 2 + 1,
(3.16)
√
u2
+v + 2
m −2
= 0, and the simple pole (3.17)
K. Watanabe / Wave Motion 51 (2014) 1209–1224
1217
√ 2 2 −2 which is the zeroth of the denominator, (M /m)(u2 + v 2 ) + √ u + v + m = 0 in F (u, v). Then, two branch cuts along the imaginary axis are introduced so that the conditional Re u2 + v 2 + m−2 > 0 can be satisfied. For the complex integral Φ , the closed loop L shown in Fig. 3.2 is employed. Examining the denominator of F (u, v), we learn that the pole is within the loop only when |usaddle | > |uP |. This inequality gives two conditionals. They are 2 Mr − 1 and Mr > Z + l. (3.18) 0
√
> 0 and Im(u) > 0, and then the line integral along the real axis, i.e. the integral I (v; p), is converted √ to the sum of the integral along the Cagniard path and the residue at the pole uP = +i v 2 + 1. That is ∞ du(+) du(−) F (u(+) , v) I (v; p) = − F (u(−) , v) exp(−pt )dt Re
u2 + v 2 + m−2
dt
TL (v)
+H
Mr
dt
−1 H
Z +l
2
Mr
− 1 − v 2π i Res(uP ).
Z +l
(3.19)
The residue is evaluated exactly as
2π i Res(uP ) = 2π i
u − i v 2 + 1 F (u, v)
× exp −(p/cs ) 2π m
=
2
√ (Z + l) u2 + v 2 + m−2 − (M /m)(Z − l) − i(mr )u
1 + M2
√
m
u=+i
v 2 +1
2lM 2 exp −(p/cs ) +r v +1 . 2
(3.20)
m
Consequently, we have for the inner integral I (v; p) =
∞
TL (v)
+H
F (u(+) , v)
du(+) dt
Mr Z +l
− F (u(−) , v)
−1 H
dt
exp(−pt )dt
2π m2 2lM 2 −1−v exp −(p/cs ) +r v +1 . 1 + M2 m2
2
Mr
du(−)
Z +l
(3.21)
Substituting the above equation into the double integral of Eq. (3.7), we have for the Laplace transformed displacement as ⊗(r ) uθ
=
+H
8π 2 µcs
∞
T
∞
TL (v)
0
F (u(+) , v)
du(+)
− F (u(−) , v)
dt
du(−)
dt
exp(−pt )dt
2 Mr 2π m2 2lM −1 H − 1 − v exp −(p/cs ) + r v 2 + 1 dv. Z +l Z +l 1 + M2 m2
Mr
(3.22)
The order of integration is exchanged for the first double integral whose supporting region is shown in Fig. 3.3. As for the second term, the variable transform from v to the time t
√
t =
2lM /m + (mr ) v 2 + 1 cs m
,
v = vP (t ) =
cs t − 2lM /m2
2
r
− 1,
(3.23)
is introduced. Then, Eq. (3.22) yields to
⊗(r )
uθ
=
T 8π 2 µcs
+H
∞
exp(−pt )dt
tL
Mr Z +l
vL (t )
F (u(+) , v)
0
−1
2π m2 1 + M2
dt
− F (u(−) , v)
cs cs t − 2lM /m2
tQ tP
du(+)
r
cs t − 2lM /m2
2
du(−)
dt
dv
exp(−pt )dt ,
− r2
(3.24)
1218
K. Watanabe / Wave Motion 51 (2014) 1209–1224
where tL = TL (0) =
1
(mr )2 + (Z + l)2 − M (Z − l) ,
m2 cs
(3.25)
{(mcs t ) + (M /m)(Z − l)}2 1 − 2, (mr )2 + (Z + l)2 m 1 2Ml Mr 1 2Ml +r , tQ = + r . tP = cs m2 cs m2 Z +l vL (t ) =
(3.26)
(3.27)
Eq. (3.24) is the Laplace transformed displacement, and it is just in the form of the Laplace transform integral with respect to the parameter p. Thus, we can apply the simple Laplace inversion formulas, L
−1
∞
f (t ) exp(−pt )dt
= H (t − a)f (t ),
a
L
−1
b
f (t ) exp(−pt )dt
(3.28)
= H (t − a)H (b − t )f (t ),
a
and can perform the Laplace inversion by inspection. It yields
T
(r )
uθ (r , Z , t ) =
H (t − tL )
8π 2 µcs
F (u(+) , v)
du(+)
0
Mr
+ H
vL (t )
Z +l
− 1 H (t − tP )H (tQ − t )
− F (u(−) , v)
dt 2π m
dt
dv
cs cs t − 2lM /m2
2
1 + M2
du(−)
r
cs t − 2lM /m2
2
.
(3.29)
− r2
In order to discuss the wave front behavior, the torsional displacement of the reflected wave has to be simplified. Using the relation/identity, 1
du(±)
u2(±) + v 2 + m−2
dt
= ±
mcs
{mcs t + (M /m)(Z − l)} − {(mr )2 + (Z + l)2 }(v 2 + m−2 ) 2
,
(3.30)
the integrand in Eq. (3.29) is rewritten, and then we introduce the variable transform from v to ϕ ,
v = vL (t ) sin ϕ,
(3.31)
for the integral. We have the simpler expression for the integral as vL (t )
F (u(+) , v)
du(+) dt
0
− F (u(−) , v)
du(−)
dt
dv =
(mr ) + (Z + l) 2
π /2
mcs 2
G(u(+) , v) + G(u(−) , v) dϕ,
(3.32)
0
where
G(u, v) = F (u, v) u2 + v 2 + m−2
= (M /m)
√
u2
+v + 2
m −2
(−iu) √ . − M /m u2 + v 2 + m−2 + 1/(Mm)
(3.33)
The variable of Cagniard path u(±) is also simplified as
u(±)
+i(mr ) vL2 (t ) + m−2 ± (Z + l)vL (t ) cos ϕ = . (mr )2 + (Z + l)2
(3.34)
The second term in Eq. (3.29) is also reduced to cs cs t − 2lM /m2
2π m2 1+
M2
r
cs t − 2lM /m2
2
= − r2
2π m2 1+
M2
cs (cs t − cs tP + r ) . √ r (cs t − cs tP + 2r ) (cs t − cs tP )
(3.35)
K. Watanabe / Wave Motion 51 (2014) 1209–1224
1219
Fig. 3.3. The supporting region for the double integral in Eq. (3.22).
Then, Eq. (3.29) yields to the simpler form (r ) uθ ( r , Z , t )
=
T
H (t − tL )
8π µ 2
+H
Mr Z +l
π /2
m
(mr ) + (Z + l) 2
2
− 1 H (t − tP )H (tQ − t )
G(u(+) , v) + G(u(−) , v) dϕ
2π m2
cs (t − tP ) + r
0
1 + M2 r
. √ {cs (t − tP ) + 2r } cs (t − tP )
(3.36)
Fortunately, we can evaluate the integral in the above equation. Its evaluation procedure is described in Appendix C, since it takes a little bit labor time. The final result for the reflected wave is thus given by 8π µ (r ) 1 − M2 m uθ ( r , Z , t ) = H ( t − t L ) T 1 + M 2 Rm
{vm (t ) cos ψ − (M /m)} vm (t ) sin ψ + |vm (t ) cos ψ − (M /m)| cot ψ {vm (t ) cos ψ − (M /m)}2 + vL2 (t ) sin2 ψ vm (t ) sin ψ + {vm (t ) cos ψ + 1/(Mm)} cot ψ − {vm (t ) cos ψ + 1/(Mm)}2 + vL2 (t ) sin2 ψ
sgn ×
+ H (Mr − (Z + l)) H (t − tP )H (tQ − t )
1 − M2 1+
M2
2 {cs (t − tP ) + r } , √ r {cs (t − tP ) + 2r } cs (t − tP )
(3.37)
where the new variables (Rm , ψ) and the function vm (t ) are defined by Eqs. (C.3) and (C.6). 4. Wave front singularity and discussions Now, we shall discuss wave propagation behaviors based on the final closed form solution. The incident wave is produced by the needle-like source and its exact expression is given by Eq. (A.1) in Appendix A. Inspecting Eq. (A.1), the displacement is composed of two waves. One is the cylindrical wave with the product of two step functions, H (−z )H (t − r /cs ). This wave disturbs the negative z-region only and does not interact with the moving boundary as far as we assume the outward motion of the boundary, M = V /cs√> 0. The other is the spherical wave with H (t − R/cs ), where R is the spherical distance from the coordinate origin, R = r 2 + z 2 . This spherical wave radiates from the source edge to all directions and catches up to the moving boundary at the time t =
l cs (1 − M )
.
(4.1)
We call this wave ‘‘incident wave’’ for our moving boundary problem. Clearly, this incident wave shows the finite jump at the front
(i)
uθ
t →R/cs
=
T
r
8π µ zR
.
(4.2)
The wave reflection starts from the time of Eq. (4.1) and the reflected wave is given by the first term in Eq. (3.37). Its disturbed region is determined by the operation of the step function H (t − tL ) and its front t = tL is the spherical shape defined by
r 2 + Z − M (cs t ) −
1 + M2 m2
2 l
2 2Ml = (cs t ) − 2 , m
(4.3)
1220
K. Watanabe / Wave Motion 51 (2014) 1209–1224
in the moving coordinate system (r , Z , t ), or
2
r +
z−
2
2l
= cs t −
m2
2Ml
2
m2
,
(4.4)
in the original fixed coordinate system (r , z , t ). It is interesting that the center of the reflected spherical wave is fixed at z = 2l/m2 in spite of the reflector motion. Examining the first term in Eq. (3.37) near the wave front, we get the finite jump at the front,
(r )
uθ ( r , Z , t )
=
t → tL
T
1 − M 2 m sin ψ
8π µ 1 + M 2
1 cos ψ − M
Rm
−
1 cos ψ + 1/M
.
(4.5)
This jump becomes infinite when cos ψ − M = 0. This condition is equivalent to Z + l = Mr and thus the infinite jump takes place at the apex A of the reflected wave in Fig. 4.2. The second term in Eq. (3.37) shows a new wave whose disturbed region is determined by the product of three step functions,
H
Mr Z +l
− 1 H (t − tP )H (tQ − t ).
(4.6)
Examining the argument of each step function, we learn that the step function H (t − tP ) shows the wave nature and its front t = tP is the cylindrical shape r+
2Ml m2
= cs t .
(4.7)
Other step functions do not show the wave nature, but restrict the disturbed region of the new wave. The region by the product of two step functions, H (t − tP )H (tQ − t ) is included in the region by H (Mr /(Z + l) − 1). Thus the conditional Mr /(Z + l) > 1 does not give any meaning for the wave field. The substantial conditional for the new wave is the product of the former two step functions. When we considered the 2D anti-plane deformation [11,12], the conditional H (tQ − t ) did not appear and only the product H (Mr /(Z + l) − 1) H (t − tP ) determined the disturbed region. The conditional H (tQ − t ) is the characteristic of the 3D wave. The similar non wave front conditional can be found in the case of the 3D elastodynamic problem [18,21]. The new wave appears and starts when the front of the incident wave crosses the apex of the reflected wave just on the moving boundary (see Fig. 4.1(d)). It is the time t =
1 + M2 1−
M2
l cs M
,
(4.8)
when the moving boundary plane (flat vertical line) passes through the center of the spherical reflected wave as shown in Fig. 4.1(d). At this instant, the incident ray (the dotted thin arrow in Fig. 4.1(d)) has the critical incident angle sin−1
1 − M2
1 + M2
,
(4.9)
and its reflected ray is parallel to the plane of the moving boundary. This critical reflection has been shown for the timeharmonic wave [10]. Thus, the new cylindrical wave is produced by the incident ray with the critical incident angle. After the starting time of the new wave, the cylindrical front moves radially outward direction along the moving boundary. Figs. 4.1(e) and 4.2 show this new wave and the disturbed region is a torus with triangular-like cross section which is formed by the product of two step functions, H (t − tP )H (tQ − t ), and the flat moving boundary as shown by the region ABC in Fig. 4.2. In the figure the radial distances, rP , rL and rQ , on the moving boundary are defined by rP = cs t −
2Ml m2
,
rL =
m2
(cs t ) − 2Ml(cs t ) − , 2
l2
rQ =
l M
cs t −
2lM m2
.
(4.10)
If the boundary is immovable, the product of the three step functions does not give any region and thus this cylindrical wave disappears. The wave front behavior is easily extracted from the second term of Eq. (3.37). That is the inverse square root singularity, (r ) uθ ( r , Z , t )
t → tP
=
T
1 − M2
8π µ 1 + M 2
2 r
√
1
cs (t − tP )
.
(4.11)
It will be worth to mention that the new wave has the stronger singularity at its front than that of the regular reflected wave which has the finite jump. The similar wave nature is also found in the case of the 2D SH-wave [11,12]. For visual understanding, the wave front development is shown in Fig. 4.1 and the details around the new cylindrical wave is in Fig. 4.2, schematically.
K. Watanabe / Wave Motion 51 (2014) 1209–1224
(a) 0 < cs t < 1−l M .
M2 (c) 1−l M < cs t < 11+ −M 2
1221
(b) cs t = 1−l M .
l M
+M 2 (d) cs t = 11− M2
.
2
M (e) cs t > 11+ −M 2
l M
l M
.
.
Fig. 4.1. Schematic time development of wave front shape.
Finally, a numerical calculation is carried out for the surface response. The total surface response is given by (i)
(r )
uθ = uθ (r , z , t )|z =Mcs t +l + uθ (r , Z , t )|Z =0 .
(4.12)
Fig. 4.3 shows the response of the torsional displacement at several observing points on the moving surface. When the observing point (radius) is smaller than r /l = 1/M, the first arrival is the edge of the incident and regular reflected waves and the new cylindrical wave is not observed. The edge is at r = rL in Fig. 4.2. On the other hand, when the observing point is larger than r /l = 1/M, the first arrival is the new cylindrical wave and the second is the edge of the incident and reflected waves. The arrival of the wave edge is denoted by the vertical arrow and that of the new cylindrical wave is marked by the
1222
K. Watanabe / Wave Motion 51 (2014) 1209–1224
Fig. 4.2. Details around the new cylindrical wave H (t − tP )H (tQ − t ).
Fig. 4.3. Response of the torsional displacement at the moving surface (arrow: arrival of the wave edge tL , ‘‘∞’’: arrival of tP ).
infinity ‘‘∞’’ in the figure. At the critical observing point, r /l = 1/M, the two arrivals are at one time. We conclude that the new cylindrical wave cannot be observed, when the surface observing point is smaller than the critical radius r /l = 1/M. 5. Conclusion The reflection of the torsional wave by a uniformly moving flat boundary is discussed. An exact closed form solution for the reflected wave is obtained. Discussing wave propagation phenomena, a new cylindrical wave is found in addition to the regular reflected wave. It is produced by the critical ray in the incident wave and runs along the moving boundary after reflection. The new wave has the cylindrical front, but the disturbed region is one of torus with triangular-like cross section. The wave front singularity of the new wave is stronger than those of the incident and the regular reflected waves. The appearance of the new cylindrical wave is one of characteristic wave phenomena for the elastodynamic moving boundary problem. However, the energy break into the regular reflected and the new waves is not discussed here. So, the problem of this energy decomposition is open. Acknowledgment The author would like to express his thanks to the reviewer for many invaluable comments and suggestions which are included in the paper.
K. Watanabe / Wave Motion 51 (2014) 1209–1224
1223
Appendix A. Exact expressions for the incident wave [15]
(i) uθ
T
=
H (−z )H (t − r /cs )
8π µr
σz(θi) = −
T
r
8π R3
σr(θi) = −
2(cs t )
(cs t )2 − r 2
+ H (t − R/cs ) sgn(z )
(cs t ) (cs t )2 − r 2
−
z R
,
H (t − R/cs ),
(A.2)
T
2(cs t )
H (−z )H (t − r /cs )
8π r 2
(A.1)
2−
r2
(cs t )2 − r 2 (cs t )2 − r 2 r2 (cs t ) z z2 2− + H (t − R/cs ) sgn(z ) − 2+ 2 , (cs t )2 − r 2 R R (cs t )2 − r 2
(A.3)
where R=
r 2 + z2.
(A.4)
The displacement of the incident wave shows a finite jump on the spherical wave front as
(i)
uθ
t →R/cs
=
r
T
8π µ zR
.
(A.5)
Appendix B. Conversion of the Laplace inversion integral The conversion formula is established by the author [11,12] and its mathematical procedure is explained concisely in [16]. That is 1 2π i
Br.(s)
f (s) exp
1 ds 2 2 s − V ξ + (s/cs ) t ds = f (s) exp(pt )dp, 2π i Br.(p) dp s=s(p)
(B.1)
where the integrand f (s) vanishes at infinity, i.e. f (s) = 0; |s| → ∞ and the variable transform is defined by
s(p) =
ξ 2 + (s/cs )2 , p + cs M (mξ )2 + (p/cs )2
p=s−V
m2
(B.2)
.
(B.3)
Two Bromwich paths are Br.(s) : from γs − i∞ to γs − i∞, Br.(p) : from γp − i∞ to γp − i∞, where γs and γp are sufficiently large constants and have the relation
γp = (1 − M )γs .
(B.4)
Appendix C. Evaluation of the integral in Eq. (3.36) We consider the integral, I (r , Z , t ) =
π/2
m
(mr ) + (Z + l) 2
2
G(u(+) , v) + G(u(−) , v) dϕ.
(C.1)
0
Using partial fractions, we decompose the integrand as G(u, v) =
1 − M2 1 + M2
√
(−iu)
u2 + v 2 + m−2 − M /m
−√
(−iu)
u2 + v 2 + m−2 + 1/(Mm)
,
(C.2)
and introduce the polar coordinate-like variables (Rm , ψ) as mr = Rm sin ψ,
Z + l = Rm cos ψ;
Rm =
(mr )2 + (Z + l)2 .
(C.3)
1224
K. Watanabe / Wave Motion 51 (2014) 1209–1224
The variable −iu(±) and radicals are simplified as
−iu(±) = vm (t ) sin ψ ∓ ivL (t ) cos ψ cos ϕ, u2(±) + v 2 + m−2 = vm (t ) cos ψ ± ivL (t ) sin ψ cos ϕ,
(C.4) (C.5)
where
vm (t ) =
vL2 (t ) + m−2 .
(C.6)
Since
vm (t ) sin ψ − ivL (t ) cos ψ cos ϕ 1 + M2 vm (t ) cos ψ − (M /m) + ivL (t ) sin ψ cos ϕ vm (t ) sin ψ − ivL (t ) cos ψ cos ϕ , − Re vm (t ) cos ψ + 1/(Mm) + ivL (t ) sin ψ cos ϕ
G(u(+) , v) + G(u(−) , v) = 2
1 − M2
Re
(C.7)
the integral I (r , Z , t ) is rewritten as π/2
{vm (t ) cos ψ − (M /m)} vm (t ) − vL2 (t ) cos ψ cos2 ϕ dϕ 1 + M 2 Rm {vm (t ) cos ψ − (M /m)}2 + vL2 (t ) sin2 ψ cos2 ϕ 0 π/2 {vm (t ) cos ψ + 1/(Mm)} vm (t ) − vL2 (t ) cos ψ cos2 ϕ − d ϕ . {vm (t ) cos ψ + 1/(Mm)}2 + vL2 (t ) sin2 ψ cos2 ϕ 0
I (r , Z , t ) = 2
1 − M 2 m sin ψ
(C.8)
We apply the simple integration formulas, π/2
dϕ
√
π
,
π/2
cos2 ϕ
dϕ =
π
a2 + b2 cos2 ϕ 2b2 2a a2 + b2 0 to the integrals in Eq. (C.8). It results in the closed expression for the integral, 0
a2 + b2 cos2 ϕ
=
a 1− √ a2 + b 2
;
a > 0, b > 0, (C.9)
− (M /m)}vm (t ) sin ψ + |vm (t ) cos ψ − (M /m)| cot ψ sgn{vm (t ) cos ψ {vm (t ) cos ψ − (M /m)}2 vL2 (t ) sin2 ψ vm (t ) sin ψ + {vm (t ) cos ψ + 1/(Mm)} cot ψ − , {vm (t ) cos ψ + 1/(Mm)}2 + vL2 (t ) sin2 ψ
I (r , Z , t ) = π
1 − M2 m
1 + M 2 Rm
(C.10)
where the sign function is defined by sgn(x) =
+1; −1;
x > 0, x < 0.
(C.11)
References [1] D. Censor, Scattering of a plane wave at a plane interface separating two moving media, Radio Sci. 4 (11) (1969) 1079–1088. [2] D. Zutter, Reflection from linearly vibrating objects: plane mirror at oblique incidence, IEEE Trans. Antennas and Propagation AP-30 (5) (1981) 898–903. [3] Yao-Xiong Huang, Reflection and transmission of electromagnetic waves by a dielectric medium moving in an arbitrary direction, J. Appl. Phys. 76 (5) (1994) 2575–2581. [4] D. Censor, Scattering by time varying obstacles, J. Sound Vib. 25 (1) (1972) 101–110. [5] O. Bou Matt, J.P. Remenieras, C. Brunel, A. Roncin, F. Patat, Ultrasonic sensing of vibration, Ultarsonics 36 (1998) 391–396. [6] R. Wunenburger, N. Mujica, S. Fauve, Experimental study of the Doppler shift generated by a vibrating scatterer, J. Acoust. Soc. Am. 115 (2) (2004) 507–514. [7] D. Censor, The generalized Doppler effect and applications, J. Franklin Inst. 295 (2) (1973) 103–116. [8] S.-Y. Lee, M. Lee, A new wave technique for free vibration of a string with time-varying length, Trans. ASME, Ser. E, J. Appl. Mech. 69 (2002) 83–87. [9] H. Kimura, H. Ito, T. Nakagawa, Vibration analysis of elevator rope: forced vibration of rope with time-varying length, Trans. JSME, Ser. C 71 (706) (2005) 1871–1876 (in Japanese). [10] Kazumi Watanabe, Elastodynamic Doppler effects by a moving stress-free edge, Internat. J. Engrg. Sci. 48 (2010) 1995–2014. [11] Kazumi Watanabe, Reflection of an impulsive SH-wave at a moving edge, Quart. J. Mech. Appl. Math. 63 (3) (2010) 335–348. [12] Kazumi Watanabe, Transient SH-waves in an elastic plate with a moving edge, Acta Mech. 223 (8) (2012) 1823–1836. [13] K. Watanabe, R. Payton, Green’s function for torsional waves in a cylindrically monoclinic material, Internat. J. Engrg. Sci. 43 (2005) 1283–1291. [14] K. Watanabe, Torsional waves produced by a semi-infinite needle-like torque source, ZAMM Z. Angew. Math. Mech. (2014) http://dx.doi.org/10.1002/ zamm.201400043. in press. [15] I.N. Sneddon, Fourier Transforms, McGraw-Hill, London, 1951. [16] K. Watanabe, Integral Transform Techniques for Green’s Function, Springer, Switzerland, 2014. [17] A.T. De Hoop, A modification of Cagniard’s method for solving seismic pulse problems, Appl. Sci. Res., Sect. B 8 (1960) 349–356. [18] D.C. Gakenheimer, Response of an elastic half space to expanding surface loads, Trans. ASME, Ser. E, J. Appl. Mech. 38 (1) (1971) 99–110. [19] S. Moriguchi, K. Udagawa, S. Ichimatsu, Mathematical Formulas, Vol. III, Iwanami, Tokyo, 1987 (in Japanese). [20] G.N. Watson, Theory of Bessel Functions, second ed., Cambridge, 1966, p. 20. [21] F.R. Norwood, Exact transient response of an elastic half space loaded over a rectangular region of its surface, Trans. ASME, Ser. E, J. Appl. Mech. 36 (3) (1969) 516–522.