Scattering of Rayleigh and Stoneley waves by two adhering elastic wedges

Scattering of Rayleigh and Stoneley waves by two adhering elastic wedges

Wave Motion 33 (2001) 321–337 Scattering of Rayleigh and Stoneley waves by two adhering elastic wedges Bair V. Budaev, David B. Bogy∗ Department of M...

288KB Sizes 0 Downloads 55 Views

Wave Motion 33 (2001) 321–337

Scattering of Rayleigh and Stoneley waves by two adhering elastic wedges Bair V. Budaev, David B. Bogy∗ Department of Mechanical Engineering, University of California, Berkeley, CA 94720, USA Received 22 July 1999; received in revised form 24 July 2000; accepted 2 August 2000

Abstract The method of Sommerfeld integrals is used to study propagation of Rayleigh and Stoneley waves in a system of two bonded solid wedges with a common vertex. Numerical results are obtained for configurations with a wide range of angles of steel and aluminum wedges. © 2001 Elsevier Science B.V. All rights reserved.

1. Introduction Problems of elastic wave propagation in solid wedge-shaped structures have been considered for several decades with steady but very slow progress. Lamb [13] in the early 1900s studied propagation of elastic waves in a half-space which can be considered as a 180◦ wedge. This was the first solved problem of elastodynamics which is characterized by the existence of two types of elastic waves simultaneously propagating with different wave speeds and interacting with each other through the boundary conditions and other irregularities of the media. The next recognizable step was made in the late 1940s when Friedman [10] analyzed diffraction of a plane elastic impulse by a semi-infinite plane screen with traction-free boundary conditions. Later Fillipov [9] simplified Friedman’s solution and Maue [16] solved a time-harmonic version of this problem. In the 1960s new experimental techniques made it possible to measure Rayleigh wave reflection and transmission in a single solid wedge [12] but theoretical solutions of the problem were not obtained at that time. Starting from the 1980s new interest in problems of elastic wave propagation in wedges gave rise to various techniques based on the reduction of problems of diffraction to integral equations of different kinds [11,14]. Finally, a combination of the Sommerfeld–Maliuzhinetz technique [15,18] widely used in acoustics and electromagnetics with the methods based on the reduction to integral equations made it possible to develop a uniform approach to problems of diffraction of elastic waves in wedge-shaped regions [2–5,7]. In [3,4], reflection and transmission coefficients of Rayleigh waves in a single elastic wedge were computed, and they were shown to be in agreement with experimental data from [12]. In [15], the method from [3,4] was extended to the problem of elastic wave propagation in a system of two bonded solid wedges. It was shown that all structural features of the elastic wave fields in bonded wedges could be retrieved, and the performed numerical experiments provided convincing evidence that the method could be used for numerical analysis of realistic structures. ∗

Corresponding author. E-mail address: [email protected] (D.B. Bogy). 0165-2125/01/$ – see front matter © 2001 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 2 1 2 5 ( 0 0 ) 0 0 0 7 4 - 3

322

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337

Here we continue the analysis of elastic wave propagation in a system of two bonded solid wedges started in [5]. Unlike in [5], whose main goal was to develop a general approach which was able to handle the different particular questions, such as the behavior of wave fields near the vertex and the behavior of wave fields far away from the vertex, here we focus on the computation of surface Rayleigh and Stoneley waves generated in bi-material wedges made from steel and aluminum. Restriction to this narrow set of questions makes it possible to simplify both the theoretical and numerical considerations. Although we follow closely the theoretical scheme developed in [5], we take here a freedom to modify some of the notation. Admitting that such changes may cause difficulties for readers familiar with [5], we nevertheless believe that the theory is not finalized yet and hope that the modified notation is more transparent. In Sections 2 and 3 the problem is formulated and is converted into a system of singular integral equations that are to be solved numerically. We also discuss how the characteristics of the wave fields may be recovered from the solution of the integral equations. In Section 4, the numerical scheme employed for the solution of the singular integral equations is described, and in Section 5 the obtained numerical results are discussed.

2. Formulation Let there be two isotropic and homogeneous elastic media with material properties characterized by Poisson ratios σ1,2 , shear moduli µ1,2 and by material densities %1,2 . Let the first medium occupy the wedge Γ1 = {r ≥ 0, 0 ≤ θ ≤ 2α1 , −∞ < z < ∞} given in the cylindrical coordinates {r, θ, z}, and let the second medium occupy the adjacent wedge Γ2 = {r ≥ 0, 2α1 ≤ θ ≤ 2α1 + 2α2 , −∞ < z < ∞} (Fig. 1). The two wedges Γ1 and Γ2 are assumed to be in bonded contact along the interface θ = 2α1 , while the external faces θ = 0, θ = 2α1 + 2α2 of the composite wedge Γ = Γ1 ∪ Γ2 are assumed to be free of tractions. We are concerned with the plain strain motion in the composite wedge that is excited by an incident surface wave. In the first case, we consider an incident Rayleigh wave traveling along the traction-free face θ = 0 of the medium in Γ1 toward the apex r = 0. This wave interacts with the boundaries of the wedge and generates scattered waves that consist of a number of body and surface waves including, in particular, Rayleigh waves traveling along the free faces θ = 0 and θ = 2α1 + 2α2 in the direction away from the wedges’ apex. Also, under some conditions on the material parameters, there may be generated surface Stoneley waves propagating along the interface between the media outward from the apex. We also consider an incident Stoneley wave traveling along the interface θ = 2α1 , toward the apex r = 0. The problem, therefore, is a two-dimensional one of plain strain, and we are particularly interested in the calculation of the above-mentioned scattered Rayleigh and Stoneley waves.

Fig. 1. Bonded elastic wedges.

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337

323

The time-harmonic motions of the media are described by the fields φ1,2 (r, θ ) eiΩt , ψ1,2 (r, θ ) eiΩt with the longitudinal and shear potentials φ1,2 , ψ1,2 obeying the Helmholtz equations: 2 φn = 0 ∈ Γn , ∇ 2 φn + knp

2 ∇ 2 ψn + kns ψn = 0 ∈ Γn ,

n = 1, 2

(2.1)

with the wave numbers depending only on the frequency Ω and material parameters of the media: s r %n 1 − 2σn , n = 1, 2. , knp = kns kns = Ω µn 2(1 − σn ) (n)

(n)

The radial ur and angular uθ components of displacements are assumed to be bounded everywhere in the corresponding regions Γ1 or Γ2 which, in particular, guarantees boundedness of the energy density at the vertex (n) (n) r = 0 of the wedges. The stress components trθ , tθ θ are assumed to vanish along the faces θ = 0, θ = 2α1 + 2α2 , i.e. to satisfy the traction-free boundary conditions: (1)

(1)

(2)

trθ (r, α1∗ ) = tθθ (r, α1∗ ) = 0,

(2)

trθ (r, α2∗ ) = tθ θ (r, α2∗ ) = 0,

α1∗ = 0, α2∗ = 2α1 + 2α2 .

(2.2)

Both the displacements and tractions are assumed to obey the bonded contact conditions (1)

(2)

(1)

(2)

(2) u(1) r (r, 2α1 ) = ur (r, 2α1 ),

trθ (r, 2α1 ) = trθ (r, 2α1 ),

uθ (r, 2α1 ) = uθ (r, 2α1 ),

tθθ (r, 2α1 ) = tθ θ (r, 2α1 )

(1)

(2)

(2.3)

along the axis θ = 2α1 separating the regions Γ1 and Γ2 . The traction-free boundaries θ = α1∗ = 0 and θ = α2∗ = 2(α1 + α2 ) support propagation of surface Rayleigh waves that can be viewed as plane waves with complex angles of propagation. Rayleigh waves propagating along the boundaries θ = αn∗ toward the apex r = 0 are represented by pairs of plane waves ∗



Rp

φnR< (r, θ ) = eiknR cos(θ−αn −iθn ) ,

ψnR< (r, θ ) = CnR eiknR cos(θ−αn −iθn ) , Rs

(2.4)

Rp

where the coefficients CnR , the wave numbers knR and complex angles iθn , iθnRs are specific material constants depending only on the corresponding Poisson ratio σn . Similarly, Rayleigh waves propagating along the faces θ = αn∗ outward from the apex r = 0 are represented by the pair of plane waves ∗

Rp

φnR> (r, θ ) = eiknR cos(θ−π−αn −iθn ) ,



ψnR> (r, θ ) = CnR eiknR cos(θ−π −αn −iθn

Rs )

(2.5)

with the same parameters as above. If the material parameters satisfy certain conditions [20] then there exist Stoneley waves propagating along the interface θ = α ∗ ≡ 2α1 between the media. These waves have a structure ∗

Sp

S ikS cos(θ−αn −iθn ) e , φnS< (r, θ ) = Cnp ∗

Sp

S ikS cos(θ−π−αn −iθn ) e , φnS> (r, θ ) = Cnp



S ikS cos(θ−αn −iθn ) ψnS< (r, θ ) = Cns e , Ss



S ikS cos(θ−π −αn −iθn ψnS> (r, θ ) = Cns e

Ss )

(2.6)

similar to the structure (2.4) and (2.5) of the Rayleigh waves with such material constants as the Stoneley wave Sp S , C S , normalized, for definiteness, by C S = 1. number kS , ‘angles’ θn , θnSs and the coefficients Cnp np 1p Let the motion of the composite medium be excited by an incoming Rayleigh wave (φ1in , ψ1in ) arriving from infinity along the traction-free face θ = 0 of the wedge Γ1 . Then, the total fields generated in Γ = Γ1 ∪ Γ2 may be considered as superpositions φ1 = φ1in + CR φ1R> + CS φ1S> + φ1d ,

φ2 = CT φ2R> + CS φ2S> + φ2d ,

ψ1 = ψ1in + CR ψ1R> + CS ψ1S> + ψ1d ,

ψ2 = CT ψ2R> + CS ψ2S> + ψ2d

(2.7)

324

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337

involving incoming incident Rayleigh waves (φ1in , ψ1in ), diffracted body waves (φnd , ψnd ), outgoing Rayleigh waves CR (φ1R> , ψ1R> ), CT (φ2R> , ψ2R> ), propagating along the faces of Γ1 and Γ2 , respectively, and, under some conditions, outgoing Stoneley waves CS (φ1S> , ψ1S> , φ2S> , ψ2S> ), propagating along the interface between Γ1 and Γ2 . For obvious physical reasons the coefficient CR is usually referred to as the Rayleigh wave reflection coefficient, CT is referred to as the Rayleigh wave transmission coefficient, and CS is the Stoneley wave transmission coefficient. The goal of this paper is to compute the Rayleigh wave coefficients CR , CT and the Stoneley wave coefficient CS for various sets of wedges made from common materials.

3. Representation of the solution Following the general scheme from [5,6], we determine the bisectrices θ = θn : θ1 = α1 ,

θ2 = 2α1 + α2

of the wedges Γ1 , Γ2 , and represent potentials φn (r, θ ), ψn (r, θ ) in these wedges as combinations φn (r, θ ) = φn+ (r, θ − θn ) + φn− (r, θ − θn ),

ψn (r, θ ) = ψn+ (r, θ − θn ) + ψn− (r, θ − θn )

(3.1)

of the components φn+ (r, θ ) = 21 [φn (r, θn + θ) + φn (r, θn − θ )],

ψn+ (r, θ ) = 21 [ψn (r, θn + θ ) − ψn (r, θn − θ )],

φn− (r, θ ) = 21 [φn (r, θn + θ) − φn (r, θn − θ )],

ψn− (r, θ ) = 21 [ψn (r, θn + θ ) + ψn (r, θn − θ )],

which may be considered as wave potentials symmetrized with respect to the bisectrix θ = θn . The potentials φn+ , ψn+ (n) (n) are referred to as symmetric fields, because they determine displacement fields ur , uθ that are symmetric with respect to the corresponding bisectrix θ = θn . Similarly, the potentials φn− , ψn− are referred to as antisymmetric fields because they determine displacements that are antisymmetric with respect to θ = θn . We seek symmetrized potentials φ± , ψ± in the form of the Sommerfeld integrals Z 1 [Φ ± (ω + θ) ± Φn± (ω − θ )] eiknp r cos ω dω, φn± (r, θ ) = 2πi C n Z 1 [Ψ ± (ω + θ) ∓ Ψn± (ω − θ )] eikns r cos ω dω, (3.2) ψn± (r, θ ) = 2πi C n taken over the standard ∪-shaped contour C in the complex ω-plane running from −π/2 + i∞ to 3π/2 + i∞ as shown in Fig. 2. From the properties of the Sommerfeld integrals [2,5,15], one may formulate certain conditions on the analytic amplitudes Φn± (ω), Ψn± (ω) that guarantee the corresponding potentials φn± (r, θ ), ψn± (r, θ ) solve the original problem of elastodynamics. Thus, functions Φ± (ω), Ψ± (ω) must behave at infinity as Φn± (ω), Ψn± (ω) ≈ O(e−p|Im ω| ),

Im ω → ∞, Re p > −1

(3.3)

have the symmetry Φn± (−ω) = ∓Φn± (ω),

Ψn± (−ω) = ±Ψn± (ω),

(3.4)

and be regular in the strips |Re ω| ≤ αn , except for a few specified singularities related with the incident waves propagating in the wedge Γ1 . This means that functions Φn± (ω), Ψn± (ω) must have primitive poles with specified

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337

325

Fig. 2. Sommerfeld contours and singularities of the amplitude functions Ψ1± (ω).

residues p

Φn± (ω) : ω = ±ωn

Ψn± (ω) : ω = ±ωns

p

with Res[Φn± (ωn )] = A± n, with Res[Ψn± (ωns )] = Bn± .

(3.5)

p

In the case of the incident Rayleigh waves, poles ωn , ωns are located at the points p

Rp

ωn ≡ αn − iθn , where

θnRs

ωns ≡ αn − iθnRs ,

(3.6)

are the positive roots of the equations

det(TRn (iξ ))

= 0,

TRn (ω)

=

n (ω) t11

n (ω) t12

n (ω) t21

n (ω) t22

!

with coefficients (n)

t11 (ω) =

2) sin ω(2 cos2 ω − γns q , 2 − cos2 ω γnp

(n)

t21 (ω) = sin 2ω,

(n)

t22 (ω) = −

(n)

t12 (ω) = sin 2ω, 2) sin ω(2 cos2 ω − γns p , 2 − cos2 ω γns

(3.7)

depending on the ratios knp kns , γns = k k of the wave numbers kmp , kms from (2.1) and of a reference wave number k that may be defined, e.g. as k = max{k1s , k2s }. Rp The value of θn in (3.6) is defined as   kns Rp cosh(θnRs ) , θn = arccosh knp γnp =

± and the residues A± n , Bn must be proportional to each other as nontrivial solutions of the degenerate homogeneous system of algebraic equations ! A± n n Rs = 0, TR (iθn ) Bn±

326

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337

so that the values of Bn± may be pre-selected. Thus, the incident Rayleigh wave arriving along the open boundary θ = 0 of the wedge Γ1 corresponds to the following choice of the residues Bn± : B1± = 1,

B2± = 0.

In the case of an incident Stoneley wave, the poles ωns from (3.5), are determined as p

Sp

ωn ≡ αn − iθn ,

ωns ≡ αn − iθnSs ,

(3.8)

where θ1Ss is the positive root of the equation det(TS (iξ )) = 0 with the matrix 

(1)

µ1 t11 (ω)

(1)

µ1 t12 (ω)

  (1) (1)  µ1 t21 (ω) µ1 t22 (ω)  TS (ω) =   (1) (1) t12 (ω)  s11 (ω)  (1)

t21 (ω)

(1)

(2)

µ2 t11 (ω) (2)

−µ2 t21 (ω) (2)

s11 (ω) (2)

s22 (ω)

−t21 (ω)

(2)

µ2 t12 (ω)



  (2) −µ2 t22 (ω)   ,  (2) t12 (ω)   (2)

−s22 (ω)

whose coefficients are defined by (3.7) and by formulas (n) (n) (n) (ω) = tmm (ω) + mm (ω), smm

γ 2 sin ω (n) , 1 = q ns 2 − cos2 ω 2 γnp

−γ 2 sin ω (n) 2 = p ns , 2 − cos2 ω 2 γns

(n = 1, 2).

(3.9)

Sp

Values of θ2Ss , θ2 in (3.8) are defined as  θ2Ss = arccosh

 k1s cosh(θ1Ss ) , k2s

 Sp

θn = arccosh

 k1s cosh(θ1Ss ) , knp

± and the residues A± 1,2 , B1,2 are defined as a nontrivial solution of the degenerate homogeneous system of linear algebraic equations



A± 1



   B±    1  TS (iθ1Ss )   ±=0  A2    B2± formulated for each sign ± and normalized by conditions B1± = 1. The next and the most important step of the Maliuzhinetz approach employed here is to express the boundary conditions through some functional equations with respect to analytic functions Φ± (ω), Ψ± (ω). The derivation of such equations is discussed in detail in [2,5,15] and, for brevity, we will not reproduce it here. According to [5], the functional equations representing the boundary conditions may be conveniently formulated in terms of the following vector columns:

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337







Φ1+ (g11 (ω) ± α1 )   +  Ψ (g12 (ω) ± α1 )    1   −  Φ (g11 (ω) ± α1 )    1   −  Ψ (g12 (ω) ± α1 )    1  Z± (ω) =   Φ + (g21 (ω) ± α2 )  ,   2   +  Ψ (g22 (ω) ± α2 )    2    Φ − (g21 (ω) ± α2 )    2   − Ψ2 (g22 (ω) ± α2 )



 +   Ψ (ω)   1   −   Φ (ω)   1   −   Ψ (ω)   1  Z(ω) =  + ,  Φ (ω)   2   +   Ψ (ω)   2     Φ − (ω)   2 



 γ11  γ12     γ11     γ12   , γ˜ =    γ21   γ22     γ21  γ22



 α1  α1     α1     α1   α˜ =   α2  , (3.10)    α2     α2  α2

Ψ2− (ω)

and auxiliary functions



gmn (ω) ≡ g(ω; γmn ),

Φ1+ (ω)

327

g(ω; γ ) = arccos

 1 cos ω . γ

(3.11)

In the introduced notation the Maliuzhinetz functional equations equivalent to the boundary conditions of the original problem of diffraction have the following vector form: [T+ (ω) + E+ (ω)]Z+ (ω) + [T− (ω) + E− (ω)]Z− (ω) = Q(ω)

(3.12)

with matrices T± (ω) having the block-structure     T1 (n) (n) t t   T1 11 12  t  , , T± (ω) = S± Tn (ω) =    (n) (n) T2 t21 t22 T2   E1 ±E1 −E2 ±E2     (n)  E1 ±E1 −E2 ±E2   0 1  e  , En (ω) =  E± (ω) = S±  , (n)  E1 ±E1 −E2 ±E2  0    2 E1

±E1

−E2

±E2

t and S e are constant diagonal matrices: where S± t = diag[1, ±1, ±1, 1, 1, ±1, ±1, 1], S± µn , n = 1, 2. µˆ n = (µ2 − µ1 )

e S± = diag[µˆ 2 , µˆ 2 , µˆ 2 , µˆ 2 , µˆ 1 , µˆ 1 , −µˆ 1 , −µˆ 1 ],

As for the right-hand side Q(ω) of (3.12) this is a vector-column whose elements Qn (ω) are arbitrary odd second-order trigonometric polynomials Qn (ω) = Cn sin ω + Dn sin 2ω.

(3.13)

Functional equations (3.12) provide a clear method of analytical continuation of the functions Φ± (ω), Ψ± (ω) from some fundamental strip onto the entire complex plane. Multiply Eq. (3.12) from the right by the inverse of T+ (ω) + E+ (ω) and get the formulas Z+ (ω) = K(ω)Z− (ω) + Kc (ω)Q(ω), Kc (ω) = [T+ (ω) + E+ (ω)]−1 , K(ω) = −Kc [T− (ω) + E− (ω)],

(3.14)

328

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337

which explicitly express the values Φn± (gn1 (ω) + αn ), Ψn± (gn2 (ω) + αn ), n = 1, 2, in terms of Φn± (gn1 (ω) − αn ), Ψn± (gn2 (ω) − αn ). As for the matrix K, its elements have a physical meaning as various reflection and refraction coefficients of a plane wave incident on the boundaries θ = 0, θ = 2α1 and θ = 2α1 + 2α2 . c (ω), as well as to denote Let us agree to denote the elements of the matrices K(ω) and Kc (ω) by Kmn (ω) and Kmn the elements of the vectors from (3.10) by Zm , γ˜m , and α˜ m . Then, Eq. (3.14) may be represented in the scalar form Zm (ω) =

8 X c [K˜ mn (ω)Zn (g˜ mn (ω) − α˜ n ) + K˜ mn (ω)Qn (g˜ mn (ω))],

1 ≤ m ≤ 8,

(3.15)

n=1

where K˜ mn (ω) = Kmn (g(ω − α˜ m ; 1/γ˜m )),

c c K˜ mn (ω) = Kmn (g(ω − α˜ m ; 1/γ˜m )),

g˜ mn (ω) = g(ω − α˜ m ; γ˜n /γ˜m ). Formulas (3.15) may be regarded [2,5,6] as a formalization of the well-known Snell’s rule of geometric reflection– transmission of plane waves, and, if applied recursively, (3.15) provide an explicit way to continue the functions Zm (ω) from any strips like |Re ω − a| ≤ α˜ m onto the entire complex plane. In particular, from relations (3.15) one may calculate the singularities of Xn (ω) a priori, before computing these functions. Indeed, let Zn (ω) have a ‘zero-level’ pole ω = ωn0 with the unit residue, then from (3.15) it follows that every function Zm (ω) has a ‘first-level’ pole 1 = g(ωn0 + α˜ n ; γ˜n /γ˜m ) + α˜ m ω = ωmn 1 − α 1 ) they are also explicitly at the root of the equation g(ωmn ˜ m ; γ˜n /γ˜m ) = ωn0 . As for the residues at g(ωmn determined by (3.15). Similarly, every ‘first-level’ pole of each function Zn (ω) generates ‘second-level’ poles for all functions Zn (ω), n = 1, . . . , 8, etc. Such a recursive process has a natural start up point. Note that formulas (3.5) and the definition of the auxiliary vector function Z(ω) require that functions Z5 (ω)–Z8 (ω) are regular in the strip |Re ω| ≤ α2 , and functions Z1 (ω)–Z4 (ω) must have the following poles |Re ω| ≤ α1 related to the incident waves:

Z1 (ω) : ω10 = ±(α1 − iθp1 ), Res Z1 (ω10 ) = 1, Z2 (ω) : ω20 = ±(α1 − iθs1 ),

Res Z2 (ω20 ) = ±C1R ,

Z3 (ω) : ω30 = ±(α1 − iθp1 ), Res Z3 (ω30 ) = ±1, Z4 (ω) : ω40 = ±(α1 − iθs1 ),

Res Z2 (ω20 ) = C1R

with previously defined parameters θp1 , θs1 , C1R . If we start up the recursive process from the listed poles that are located in the semi-plane Re ω > 0, then we can calculate in a finite number of steps all the singularities of all the functions Zn (ω) located in the region 0 ≤ Re ω < π + α˜ n . Furthermore, the singularities of these functions in the mirror domain −π − α˜ n < Re ω ≤ 0 may be located by use of the symmetry and (3.4), which requires the functions Z1 (ω), Z4 (ω), Z5 (ω), Z8 (ω) to be odd, and the functions Z2 (ω), Z3 (ω), Z6 (ω), Z7 (ω) to be even. After computing all singularities of Zn (ω) located in the strip |Re ω| ≤ π/2 + αˆ n we define the function X 1 Res Zn (ω∗ )σ (ω − ω∗ ), σ (ω) = , (3.16) Zn∗ (ω) = 4 sin(ω/4) ω∗ ∈Ωn (a,b)

and represent the unknown functions Zn (ω) as sums Zn (ω) = Zn∗ (ω) + Z˜ n (ω)

(3.17)

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337

329

of explicitly known pre-defined singular terms Zn∗ (ω) and of terms Z˜ n (ω) that are regular in the strip |Re ω| ≤ π/2 + αˆ n . This is the decomposition that makes it possible to convert the system of functional equations (3.12) to a system of singular integral equations available for efficient numerical analysis. Note that σ (ω) does not need to have the given form, it can be any function which is bounded at infinity and has no singularities in a sufficiently wide strip except for the primitive pole ω = 0 with a unit residue. Substitute (3.17) into (3.10), and (3.10) into the functional equations (3.12). Then after translating all pre-defined terms to the right-hand sides we obtain new inhomogeneous functional equations with respect to the analytic functions Z˜ n (ω), which are assumed to be regular in |Re ω| ≤ π/2 + α˜ n . The advantage of such inhomogeneous equations is that they may be conveniently converted to singular integral equations that determine the auxiliary functions Xn (ω) = Z˜ n (g(ω; γ˜n ) + α˜ n ) + Z˜ n (g(ω; γ˜n ) − α˜ n ),

n = 1, . . . , 8,

(3.18)

whose values along the line Re ω = π/2 are coupled with the values of the similar combinations Hn Xn (ω) = Z˜ n (g(ω; γ˜n ) + α˜ n ) − Z˜ n (g(ω; γ˜n ) − α˜ n ),

n = 1, . . . , 8,

through the operators Hn ≡ H(γ˜n ; α˜ n ) defined by Cauchy-type singular integrals Z X(ξ )g 0 (ξ ; γ ) dξ 1 π/2+i∞ , H(γ ; α)X(ω) = 2α π/2−i∞ sin(π/2α)(g(ξ ; γ ) − g(ω; γ ))

(3.19)

discussed in detail in [2,3,5]. ∗ (ω), m = 1, 2, . . . , 8, be the value of the left-hand side of the mth functional equation of the system (3.12) Let Rm with the functions Zn (ω) replaced by the known functions Zn∗ (ω), n = 1, 2, . . . , 8. Then, with use of the operators Hn we readily obtain the following system of eight singular integral equations along the line Re ω = π/2 (1)

(1)

t11 X1 + t12 X2 + µˆ 2 E1 (X) = Q1 − R1∗ , (1)

(1)

(1)

(1)

t11 H3 X3 + t12 H4 X4 + µˆ 2 E1 (X) = Q3 − R3∗ , (2)

(2)

t11 X5 + t12 X6 + µˆ 1 E1 (X) = Q5 − R5∗ , (2)

(1)

t21 H1 X1 + t22 H2 X2 + µˆ 2 E2 (X) = Q2 − R2∗ , (1)

t21 X3 + t22 X4 + µˆ 2 E2 (X) = Q4 − R4∗ , (2)

(2)

t21 H5 X5 + t22 H6 X6 − µˆ 1 E2 (X) = Q6 − R6∗ ,

(2)

t11 H7 X7 + t12 H8 X8 − µˆ 1 E1 (X) = Q7 − R7∗ ,

(2)

(2)

t21 X7 + t22 X8 + µˆ 1 E2 (X) = Q8 − R8∗

(3.20) (3.21) (3.22) (3.23)

with asymptotically small interconnecting singular integral operators (1)

(2)

E1 (X) = 11 (X1 + H3 X3 ) − 11 (X5 − H7 X7 ),

(1)

(2)

E2 (X) = 22 (X4 + H2 X2 ) − 22 (X8 − H6 X6 ).

Integral equations (3.20)–(3.23) are to be solved numerically, by a numerical scheme discussed in Section 4, for example. These equations determine values of Xn (w), n = 1, . . . , 8, along the line Re ω = π/2, as well as second-order trigonometric polynomials Qn (ω) of the type (3.13). Therefore, all we need to do to complete the solution of the diffraction problem is to recover from these data analytical functions Zn (ω), n = 1, . . . , 8, in the domain π − αˆ n < Re ω < π + αˆ n where the Sommerfeld contour of integration in (3.2) is located. We recover first the functions Z˜ n (ω) from (3.18) which are supposed to be regular in the strips |Re ω −π/2| ≤ α˜ n . As shown in [2,3,5] these functions may be recovered by the integrals Z π/2+i∞ π 1 Xn (ξ )g 0 (ξ, γ˜n ) dξ , Re ω − < α˜ n , (3.24) Z˜ n (ω) = 4α˜ n i π/2−i∞ cos(π/2α˜ n )(g(ξ, γ˜n ) − ω)) 2 which together with (3.16) and (3.17) constructively determine functions Zn (ω) in the corresponding strips |Re ω − π/2| ≤ αˆ n . After that a recursive application of functional relations (3.15) analytically continues functions Zn (ω) to the entire complex plane.

330

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337

4. Numerical scheme The numerical evaluation of the solution developed in the preceding sections consists of two major steps. First, the system (3.20)–(3.23) of singular integral equations has to be solved. Then, the solution of the integral equations has to be transformed to characteristics of wave fields. Although both of these steps require numerical work there is a significant difference in the nature of them. Thus, the second step involves recursive applications of explicit algebraic formulas and it also requires some numerical integration of pre-defined functions. Such operations are rather simple and their evaluation is neither time nor memory consuming, although some numerical errors may be caused by numerical integration. The first step involves numerical solution of an 8 × 8 system of singular integral equations which may be both time and memory consuming and, potentially, may cause unpredictable numerical effects such as large errors or instability. Therefore, the limitations of the approach are mostly determined by the algorithms used to solve the singular integral equations (3.20)–(3.23). Numerical analysis of singular integral equations (3.20)–(3.23) may be based on one of several methods. One possible way employs the observation that the changes of variables π π  π π  ξ− , ω → x = −i tan ω− , ξ → t = −i tan 2α 2 2α 2 transform the singular integral operator Z F (ξ ) dξ π 1 π/2+i∞ , Re ω = , F (ω) → 2α π/2−i∞ sin(π/2α)(ξ − ω) 2 on the infinite line to the standard Cauchy singular integral operator on the standard unit interval [−1, 1]: s Z 1 1 − x 2 f (t) dt. f (x) → 2 −1 1 − t t − z Then, applying this transform to the system (3.20)–(3.23), we get new equations which fit the well-studied class of singular integral equations [17] for which a broad range of numerical recipes are available. Another approach to numerical solution of (3.20)–(3.23) employed here is based on a possibility to approximate singular integrals of the type (3.19) by the sums Z π/2+i∞ X f (ξ ) dξ f (ξk )∆ξk ∆ξk →0 , → sin(π/2α)(ξk − ω) π/2−i∞ sin(π/2α)(ξ − ω) where {ξk } and {ωk } are different equally spaced grids shifted with respect to each other as follows: ξk = ξ0 + ∆k,

ωk = ξk + 21 ∆

k ∈ Z.

Direct application of the double-grid quadrature formulas to Eqs. (3.20)–(3.23) is made possible by the observation that no one line of this system contains any of the functions Xn (ω), n = 1, . . . , 8, together with their transforms Hm Xn (ω). Let us fix the two parameters ∆, N , and introduce two grids ξk = ∆k,

ωk = ξk + 21 ∆,

k = −N, . . . , N

containing 2N + 1 nodes each. Then we introduce 8(2N + 1) unknown numbers: ( Xn (ξk ) if n ∈ {1, 2, 5, 6}, Xnk = Xn (ωk ) if n ∈ {3, 4, 7, 8}, which means that the grid {ξk } is used for interpolation of the functions X1 , X2 , X5 , X6 , while the other functions X3 , X4 , X7 , X8 are interpolated on the grid {ωk }.

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337

331

To use the double-grid approximation of singular integrals involved in (3.20)–(3.23), we compute the values of the first, fourth, fifth, and eighth lines of this system on the grid {ωk } and compute the other lines on the grid {ξk }. Then we get a system of linear algebraic equations 8,N X

ml 1 2 Tnk Xnk − Cm Rml − Dm Rml = −Rml ,

1 ≤ m ≤ 8, −N ≤ l ≤ N

(4.1)

n=−1, k=−N

with the vectors Rml = Rm (ζlm ),

1 Rml = sin(ζlm ),

2 Rml = sin(2ζlm ),

depending on the used grid according to the rule ( ζlm

=

ωl

if m ∈ {1, 4, 5, 8},

ξl

if m ∈ {2, 3, 6, 7}.

(4.2)

It was discussed with regard to (3.13) that the constants C1 , . . . , C8 , D1 , . . . , D8 do not contribute to the wave fields and, therefore, they can be selected for convenience. All these constants, e.g., can be set to zeros. A better option, however, is not to fix the indefinite constants C1 , . . . , C8 , D1 , . . . , D8 but to consider (4.1) as an under-determinate system of 16N + 8 equations with 16N + 32 unknowns. The presence of extra unknowns relaxes the system and makes its inversion more stable in exchange for the necessity to consider thereafter the full version of functional equations (3.15) where the trigonometric polynomials Qm might disappear if the constants C1 , . . . , C8 , D1 , . . . , D8 , were set to zeros. In general, the singular integral equations (3.20)–(3.23) may have a weak singularity at infinity [5] which may require special handling to provide acceptable accuracy of the numerical procedures. If the wedges Γ1 , Γ2 are filled by such common materials as steel, aluminum or brass, and the angles of these wedges are in the ranges from 140 to 240◦ , 60 to 150◦ , respectively, then the weak singularities turn out to be negligible and it is sufficient to discretize the integral equations (3.20)–(3.23) only on the interval −7 ≤ Im ω ≤ 7 that leads to the linear algebraic equations (4.1) with conditioning numbers of the order 10−4 , providing stable inversion with accuracy sufficient for practical needs. If it is necessary to consider certain combinations of materials and wedge angles [1,5,8] for which the weak singularities of the singular integral equations (3.20)–(3.23) are not so weak as to be neglected then special techniques should be employed [19] in order to provide sufficient accuracy of the solutions. Here, however, we do not discuss such techniques because our numerical experiments are restricted to configurations which do not require special care of the weak singularities. When the solution X˜ n (ω), n = 1, . . . , 8, of the system (3.20)–(3.23) is computed then the integrals (3.24) can be approximated by quadrature formulas Z˜ n (ω) =

N Xn (ζkn )g 0 (ζkn , γ˜n )∆ 1 X , 4α˜ n i cos(π/2α˜ n )(g(ζkn , γ˜n ) − ω) k=−N

π Re ω − < α˜ n , 2

based on the same grids (4.2) where the functions Xn (ω) were defined from the numerical solution of the integral equations. Then the formulas (3.16) and (3.17) explicitly determine the approximation of functions Zn (ω) in the corresponding strips |Re ω − π/2| ≤ αˆ n , and recursive application of the functional relations (3.15) determines the approximate values of Zn (ω) in the entire complex plane, which effectively solves the discussed problem of diffraction.

332

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337

Table 1 Parameter (kg/m3 )

Mass density Young modulus (GN/m2 ) Shear modulus (GN/m2 ) Lame’s parameter (GN/m2 ) Modulus of rigidity Poisson ratio Wave speeds ratio

Symbol

Steel

Aluminum

% E µ λ G σ γ

7850 207 80 115 0.386 0.294 0.540

2700 70 27 39 0.386 0.296 0.543

5. Numerical results We simulated scattering of surface waves in bi-material wedges made from steel and aluminum characterized by the material parameters given in Table 1. The listed media parameters are not independent, but are related with each other by the formulas r 1 σE 2(1 − σ ) µ σ = − 1, λ= , γ = , G= , E 2G (1 + σ )(1 − 2σ ) 1 − 2σ so that only three constants, say σ, E, %, may be considered as fundamental material parameters. Steel and aluminum have very different mass densities and Young moduli but they have almost identical Poisson ratios σ , moduli of rigidity G and wave speed ratios γ , that somewhat simplify the computations because equal values of γ means that instead of four different special functions g(ω; γ˜n ) with n = 1, 2, 5, 6 we need to consider only two such functions with n = 1, 2. This proportionality of parameters also insure the existence of Stoneley waves propagating along the steel–aluminum interface. We run four series of numerical experiments with these metals forming a bi-material wedge structure with the base wedge angle 2α1 ranging from 140 to 240◦ and with the attached wedge having the angle 2α2 ranging from 65 to 120–150◦ depending on the base wedge angle. In the first series, we considered scattering of a Rayleigh incident wave arriving from infinity along the traction-free face of the steel base. In the second series, we switched materials and considered scattering of a Rayleigh wave coming along the traction-free surface of the aluminum base. In the third series, we considered scattering of an incident Stoneley wave between the steel base and the aluminum attachment. In the last series, we considered scattering of a Rayleigh incident wave in a mono-material steel or aluminum wedge of the angle from 120 to 280◦ which is treated as a degenerate case of the bi-material structure. It is well-known that the Rayleigh wave scattering in a mono-material medium depends on only one material parameter which is the Poisson ratio. Therefore, since steel and aluminum have similar Poisson ratios we did not need to consider Rayleigh wave scattering in these two material separately. The results of the first three series of numerical experiments are presented in Figs. 3–5, each of which consists of six subplots with the horizontal axes representing angles of the attached wedge and with the vertical axes representing the computed quantities. Each subplot contains several individual graphs corresponding to different base wedge angles. Each series has one solid line and several dashed and dotted lines. The solid lines correspond to the base wedge angle of 180◦ , e.g. to the half-space base. The dashed and dotted lines correspond to base wedges with angles greater and less than 180◦ , respectively. All graphs in Figs. 3–6 represent amplitudes or phases of the reflection or transmission coefficients whose definitions are somewhat ambiguous and, therefore, should be specified for clarity. We compute the reflection and transmission coefficients as the ratios of the longitudinal potentials of the corresponding surface waves to the longitudinal potentials of the incident waves. These are the ratios of the surface normal displacements caused by secondary waves to the normal displacements of the incident waves with orientations of the normals fixed with respect to the direction of wave propagation. Thus, the positive normal displacements of the incident and reflected Rayleigh waves have opposite directions because these waves propagate in opposite directions.

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337

Fig. 3. Scattering of the Rayleigh wave incident from the steel.

Fig. 4. Scattering of the Rayleigh wave incident from the aluminum.

333

334

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337

Fig. 5. Scattering of the Stoneley wave.

It is important to notice that the energy of a surface wave is determined by surface displacements and by parameters of materials forming this surface. Therefore, the Rayleigh waves propagating along the surfaces θ = 0, θ = 2α1 + 2α2 of different media and the Stoneley wave propagating along the interface θ = 2α1 may carry different energy even if they have equal amplitudes of displacements, the result of which is that a small incident wave arriving from a stronger material (steel) may cause significant displacements in a weaker material (aluminum). Consider first the solid lines in Fig. 3 demonstrating scattering of the incident Rayleigh waves by an aluminum wedge bonded to the top of the steel half-space. It is evident that in these configurations the reflected Rayleigh waves are very weak and practically independent of the angle of the attached aluminum wedge. Unlike the reflected

Fig. 6. Scattering of the Raleigh wave in the steel wedge (σ = 0.294).

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337

335

waves, the transmitted Rayleigh and Stoneley waves are not negligible. The amplitudes of the Stoneley waves grow while the angle of the aluminum wedge increases up to 90◦ and then they stabilize at a constant level. Independent of the angle of the attached aluminum wedge, the Stoneley waves have practically zero phase shifts which makes it possible to consider them as direct continuations of the incident Rayleigh waves. The amplitudes of the transmitted Rayleigh waves have rather complicated dependence on the attached wedge’s angle. The amplitude has a definite minimum near the angle of 90◦ , then it grows to a maximum near 100◦ and after that it steadily decays, stabilizing at an observable level of about 50% of the maximum value. Despite the complicated behavior of the amplitudes the phases of the transmitted Rayleigh waves have very simple behavior. They have a constant value for attached wedge angles below 90◦ then they rapidly show a 180◦ phase shift and remain at the new constant level for all angles above 90◦ . Consider now the dashed and dotted lines in Fig. 3, corresponding to the cases with the steel base is a wedge of angle greater and less than a half-space, respectively. The amplitudes of the reflected Rayleigh waves steadily grow as the base wedge angle changes from 180◦ in either direction, but these waves still have a relatively low amplitude, which is practically independent of the angle of the aluminum attachment. The phases of the reflected waves also demonstrate independence of the attached wedge angle but have positive or negative shifts if the base wedge angle increases or decreases. The Stoneley waves react similarly to changes of the base wedge angle. Their amplitudes steadily decay as the base wedge angle deviates from 180◦ in either direction, and their phases have shifts proportional to the deviations of the base’s angle from 180◦ . The transmitted waves are more sensitive to changes of the base wedge angle especially if the attached wedge angles are less than 90◦ . In this case the amplitudes of the transmitted waves grow considerably as the base wedge angle decreases, and decay considerably as the base wedge angle increases. But when the attached wedge angle becomes larger than approximately 110◦ an increase of the base wedge angle does not change the amplitude of the transmitted Rayleigh wave. Similarly, if the attached wedge angle is greater than 110◦ then decrease of the base wedge angle does not change the phase of the transmitted Rayleigh wave but generates some decrease in the amplitudes. It should also be noticed that if the base wedge angle is smaller than a half-space (dashed lines) then the amplitudes of the transmitted Rayleigh waves may become very small at some values of the attached wedge angles so that the corresponding phases become indefinite. This explains the abnormal behavior of three dashed lines of the phase graphs. The second performed series of numerical experiments was similar to the first series but with the steel and aluminum interchanged, e.g. we simulated scattering of Rayleigh waves arriving along the traction-free surface of the aluminum wedge of 140–240◦ with an attached steel wedge whose angle changes from 65 to 150◦ . The results are presented in Fig. 4, which is organized in a way similar to Fig. 3. Comparing these figures one may conclude that interchanging the materials does not change the qualitative picture of the Rayleigh wave scattering. The only visible qualitative difference is that the phase lines of the transmitted waves get 180◦ shifts in opposite directions as the attached wedge angle increases beyond the level of 90◦ . There are also some quantitative differences between the two cases represented in Figs. 3 and 4. For instance, in the cases with steel bases and aluminum attachments the reflected wave propagating in steel is weak, but the transmitted wave propagating in aluminum and the Stoneley wave propagating between the materials have considerable amplitudes. On the contrary, in the cases with aluminum bases and steel attachments the reflected wave propagating in aluminum is stronger, while the Stoneley wave and the Rayleigh wave propagating in aluminum are weaker. These observation may be interpreted by a physically obvious conjecture that in the presence of two materials the weaker material (aluminum) is more responsive and has greater displacements than the stronger material (steel). It should also be noticed that the irregular behavior in Fig. 4 of the phases of the Stoneley waves near the left-hand side of the plot are not representative because the amplitudes of the corresponding Stoneley waves become very small and comparable with the numerical errors of the computations. The third series of numerical experiments simulated scattering of an incident Stoneley wave by configurations consisting of a steel base with an attached aluminum wedge. The results of these computations are presented in Fig. 5, which is organized similar to Figs. 3 and 4.

336

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337

It is interesting to notice that the Rayleigh waves generated in the stronger material (steel) of the base wedge Γ1 do not depend on the angle of the attached wedge made from the weaker material (aluminum). However, the Rayleigh wave radiated into the weaker material and the reflected Stoneley wave propagating along the interface between the two materials depend on the angle of the aluminum attachment, and they depend on the angle of the steel base. Results of the last series of computations are plotted in Fig. 6. Here we consider reflection and transmission of an incident Rayleigh wave in a wedge made from steel or any other material with Poisson’s ratio σ = 0.294. This mono-material problem has already been studied in the literature, but here it is treated as a particular case of the bi-material configuration which creates an opportunity to test the employed algorithms and codes. Rayleigh wave reflection and transmission coefficients are computed for the wedges from the interval 120–280◦ which is chosen to avoid specific complications arising in cases of smaller or larger angles. For example, if the wedge angle is narrower than 110–120◦ then some precautions must be made for evaluation of the integrals (3.24) which may have poles close to the path of integration. The analytical continuation provides a simple way to overcome such difficulties but implementation of the corresponding algorithms was beyond our current goals. Another kind of difficulty arises when the wedge angle is greater than 270–280◦ . The problem is caused by the already mentioned fact that if the wedge’s angle exceeds 180◦ then the introduced integral equations have a weak singularity at infinity. If the angle is close to 180◦ then the singularity has a negligible order, but the singularity becomes strong at the limit 360◦ . Numerical experiments show that while the wedge angle is below 270–280◦ then the integral equations can be solved without any special treatment of the weak singularities at infinity. However, for larger angles direct computations become unstable which suggests that a specific treatment of the singularities at infinity should be employed. In this regard, it should be noticed that computations reported in our previous paper [4] were inaccurate for wedges larger than 180◦ . 6. Conclusion The problem of surface wave scattering by a system of two adhering elastic waves is considered by the extension of the Sommerfeld-Maliuzhinetz technique which may be considered as a kind of boundary integral equations method. This technique appears to be an adequate tool for analysis of surface waves in wedges which are of an essentially asymptotic nature because these waves are formed only at a significant distance from the wedge’s tip measured in wave lengths. Such structure of surface waves in wedges restricts use of direct numerical methods, but the employed technique providing an analytical description of surface waves generated in wedges is shown to be effective, although it also relies on the numerical solution of some integral equations. It is shown that a simple double-grid method of inversion of singular integral equations provides sufficient accuracy for moderate configuration, excluding cases with small wedges where multiply reflected waves are generated, and excluding also cases with a strong stress concentration near the vertex. Acknowledgements This research was supported by NSF Grant 9730779. DBB holds William S. Floyd Jr. distinguished Professorship in Engineering, which also supported the work. References [1] D.B. Bogy, Two edge-bonded elastic wedges of different materials and wedge angles under surface tractions, J. Appl. Mech. 38 (1971) 377–385. [2] B.V. Budaev, Diffraction by wedges, Pitman Research Notes in Mathematics Series, Vol. 322, Longman, Essex, 1995.

B.V. Budaev, D.B. Bogy / Wave Motion 33 (2001) 321–337 [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20]

337

B.V. Budaev, D.B. Bogy, Rayleigh wave scattering by a wedge I, Wave Motion 22 (1995) 239–257. B.V. Budaev, D.B. Bogy, Rayleigh wave scattering by a wedge II, Wave Motion 24 (1996) 307–314. B.V. Budaev, D.B. Bogy, Rayleigh wave scattering by two adhering wedges, Proc. Roy. Soc. London. Ser. A 454 (1979) (1998) 2949–2996. B.V. Budaev, D.B. Bogy, Rigorous solutions of acoustic wave diffraction by penetrable wedges, J. Acoust. Soc. Am. 105 (1999) 74–83. B.V. Budaev, D.B. Bogy, Rayleigh wave scattering in a wedge with mixed boundary conditions, Proc. Roy. Soc. London. Ser. A., in press. J. Dundurs, Edge-bonded dissimilar orthogonal elastic wedges under normal and shear loading, J. Appl. Mech. 36 (1969) 650–652. A.S. Fillipov, Some problems of diffraction of elastic waves, Appl. Math. Mech. 20 (1956) 688–703. M.M. Friedman, Diffraction of a plane elastic wave by a straight strain-free screen, Dokl. Acad. Nauk. USSR 66 (1949) 21–24. A.K. Gautesen, Scattering of Rayleigh wave by an elastic quarter space, J. Appl. Mech. 52 (1985) 664–668. L. Knopoff, A.F. Gangi, Transmission and reflection of Rayleigh waves by a wedge, J. Geophys. 25 (1960) 1203–1214. H. Lamb, On propagation of tremors over the surface of an elastic solid, Phil. Trans. Roy. Soc. London A 203 (1904) 1–42. Z.L. Li, J.D. Achenbach, I. Komsky, Y.C. Lee, Reflection and transmission of obliquely incident surface waves by an edge of a quarter space: theory and experiment, J. Appl. Mech. 59 (1992) 349–355. G.D. Maliuzhinetz, Radiation of sound by oscillating faces of an arbitrary wedge, Acoust. J. 1 (1955) 144–164. A.W. Maue, Die bengung elastcher wellen an der halbebene, Z. Angew. Math. Mech. 33 (1953) 1–10. N.I. Muskhelishvili, Singular Integral Equations: Boundary Problems of Function Theory and their Application to Mathematical Physics, 2nd Edition, Dover, New York, 1992. A.V. Osipov, A.N. Norris, The Malyuzhinets theory for scattering from wedge boundaries: a review, Wave Motion 29 (1999) 313–340. S. Prossdorf, Numerical Analysis for Integral and Related Operator Equations, Birkhauser, Boston, MA, 1991. S. Timoshenko, Theory of Elasticity, 3rd Edition, McGraw-Hill, New York, 1970.