Dispersion free wave splittings for structural elements

Computers and Structures 84 (2006) 514–527 www.elsevier.com/locate/compstruc

Dispersion free wave splittings for structural elements M. Johansson *, P.D. Folkow, P. Olsson Department of Applied Mechanics, Chalmers University of Technology, SE-412 96 Go¨teborg, Sweden Received 12 December 2003; accepted 28 September 2005

Abstract Wave splittings are derived for three types of structural elements: membranes, Timoshenko beams, and Mindlin plates. The Timoshenko beam equation and the Mindlin plate equation are inherently dispersive, as is each Fourier component of the membrane equation in an angular decomposition of the field. The distinctive feature of the wave splittings derived in the present paper is that, in homogeneous regions, they transform the dispersive wave equations into simple one-way wave equations without dispersion. Such splittings have uses both for radial scattering problems in the 2D cases and for scattering problems in dispersive media. As an example of how the splittings may be applied, a direct scattering problem is solved for a membrane with radially varying density. The imbedding method is utilized, and agreement is obtained with an FE simulation. Ó 2005 Published by Elsevier Ltd. Keywords: Wave splitting; Time domain methods; Green’s operator; Imbedding; Membrane; Timoshenko beam; Mindlin plate

1. Introduction During the last decades, time domain methods such as invariant imbedding, Green’s function techniques and propagator methods, have successfully been applied to a number of scattering and wave propagation problems. A survey of these methods is contained in [1]. Both direct and inverse scattering problems have been considered, as well as problems concerning, i.e., wave guides and slender structural elements. Much of the work has been concerned with electromagnetics, but acoustics and structural mechanics have also benefited from the time domain approaches, see e.g. [2–6]. One basic building block in such time domain methods is the wave splitting; i.e., a transformation that achieves a decomposition of the wave fields into one-way waves. A number of such splittings have been derived, with different properties and applications. For some purposes, the decomposition need only diagonalize *

Corresponding author. Tel.: +46 31 772 1521; fax: +46 31 772 3827. E-mail addresses: [email protected] (M. Johansson), [email protected] (P.D. Folkow), [email protected] (P. Olsson). 0045-7949/$ - see front matter Ó 2005 Published by Elsevier Ltd. doi:10.1016/j.compstruc.2005.09.006

the principal part of the operator matrix of the relevant differential equation. Often it is advantageous to have a transformation that exactly splits the PDE into one-way wave equations in the homogeneous region. If the original wave equation is dispersive, the split wave equations in general are so too even in the homogeneous regions. In regions where the original wave equations model inhomogeneities, the one-way wave equations couple, and the coupling may often be utilized for the reconstruction of the inhomogeneities. We refer to [1,7] for how this may be done numerically. Much of the work on wave splitting has been concerned with 1D wave propagation, or situations which may naturally be reduced to 1D propagation. Weston and collaborators have, however, done much to extend the wave splittings to more general situations, see e.g. [8,9]. Moreover, Karlsson and Kristensson [10] have treated a 3D radially symmetric scattering problem, and Kreider [11,12] has solved an inverse 2D radially symmetric scattering problem. In the present paper we consider 1D and 2D wave propagation in some structural elements. In all 2D cases considered, a distinguished radial direction allows us to reduce

M. Johansson et al. / Computers and Structures 84 (2006) 514–527

the problem to 1D propagation. However, common to both the 1D and 2D wave propagation for the structural elements considered is that the wave equations exhibit dispersion, either directly geometrical or indirectly geometrical. By directly geometrical dispersion we here denote dispersion due to space dependent coefficients of the wave equations, as in the case of radial wave equations. By indirectly geometrical dispersion we denote such dispersion that arises from the modelling of the structural elements as lower dimensional objects, reducing the wave equation from 3D to 1D or 2D in space variables. Our goal is then to construct wave splittings that not only diagonalize the original wave equations, but also remove both the directly and the indirectly geometrical dispersion. Preliminary considerations on the radial membrane showed that equations exhibiting directly geometrical dispersion can not be split by means of the standard methods utilized for indirectly geometrically dispersive equations. The method which is presented in this paper, not only achieves an exact wave splitting in homogeneous regions, but also provides some freedom to choose the split wave equations. Here, we may choose dispersion free oneway wave equations in a straightforward way. This method can also be utilized to obtain wave splittings for equations that are not directly geometrically dispersive, such as the Timoshenko beam equation. The same result may also be achieved by extending the standard wave splitting by means of Green’s operators, which is demonstrated in Section 2. It should be stressed, though, that this latter approach is not an option for the radial cases. A key feature of the dispersion free wave splittings is that the Green’s operators are incorporated into the wave splitting transformation. The computational difficulties involved when wave fields are propagated are thus not eliminated entirely by the dispersion free wave splitting. All information about the dispersion of the physical fields is contained in the wave splitting transformation, which of course must be computed. In Section 3 we derive the dispersion free wave splitting on the radial membrane. The derivation can be generalized and incorporated into the analysis of the radial Mindlin plate in Section 4. The plate problem is in a sense a combination of the beam and the membrane problems, as it contains both the indirectly geometrical dispersion of the beam, and the directly geometrical dispersion of the radial membrane In Section 5 we apply the wave splitting to an inhomogeneous membrane. The imbedding equations are discussed and implemented numerically for some direct problems. The results are compared to FE solutions. Some concluding remarks are found in Section 6. 2. The Timoshenko beam This section first presents a brief review of the original wave splitting transformation together with new results on how to transform to dispersion free fields. Then follows an alternative method of obtaining dispersion free fields,

515

which has the advantage of being more easily extended to the cases of the membrane and the plate. The Timoshenko beam equation may be written [13] 2 ox c ¼ c
2 1 ot w;

ð1Þ

o2x w

ð2Þ


fc ¼

2 c
2 2 ot w;

c ¼ ox w þ w;

ð3Þ

where w(x, t), w(x, t) and c(x, t) are the mean transverse deflection, the mean rotation and the mean shear angle of the cross section, respectively. Note that in [13–15] w is defined in the opposite direction. From now on, ox = o/ox etc. The quotient of the shearing and bending stiffnesses f, the velocities c1 (effective shear velocity) and c2 (rod velocity) are defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi f ¼ k 0 GA=IE; c1 ¼ k 0 G=q; c2 ¼ E=q. Here, A and I are the area and the moment of inertia of the cross section, while q is the density of the beam. E is the modulus of elasticity, G is the shear modulus and k 0 is the shear coefficient. 2.1. Standard wave splitting The original wave splitting for the Timoshenko beam was presented by Olsson and Kristensson [14], with an amendment given by Folkow et al. [15]. It takes its starting point from Eqs. (1)–(3), written as a spatially first order system according to 0 1 0 1 0
1 1 0 w B w C B 0 0 0 1C B C B C ox w ¼ Dw; D ¼ B
2 2 C. C; w ¼ B @ c A @ c 1 ot 0 0 0A 0

2 c
2 2 ot f 0

ox w ð4Þ

An operator P and its formal inverse P and chosen so as to diagonalize D,


1

are introduced

u ¼ Pw; w ¼ P
1 u ; PDP
1 ¼ diagð
k1 ;
k2 ; k1 ; k2 Þ; ð5Þ 

þ

T ðuþ 1 ; u1 ; u2 ; u2 Þ

where u ¼ is the vector of the split fields and the ki are eigenoperators of D. The transformation þ

operators are normalised so that w ¼ uþ 1 þ u1 þ u 2 þ u2 . Hence, the Timoshenko beam Eq. (4) may be written as dispersive one-way wave equations for the split fields as  ox u i  ki ui ¼ 0;

ki ¼

1 ot þ F i ðÞ; ci

i ¼ 1; 2.

ð6Þ

In [14] explicit representations in the time domain for the ki, Fi and the elements of P and P
1 are given. The star (*) denotes time convolution, which throughout this paper is defined by Z t f ðt0 Þgðt
t0 Þ dt0 ; f ðÞ  gðÞ ¼ 0

516

M. Johansson et al. / Computers and Structures 84 (2006) 514–527

since all fields are assumed quiescent at time t < 0. As the convolutions may involve distributions such as Dirac’s delta function, the integration interval may be interpreted to range from 0
to t+. For homogeneous beams, the left and right moving fields at a position x are related to the left and right moving fields at a position x 0 respectively, through linear operators called wave propagators or Green’s operators. These are defined by the following mapping properties [15]  0 0  0 u i ðx; t  ðx
x Þ=ci Þ ¼ Gi ðx ; xÞui ðx ; tÞ.

ð7Þ

Note that there are no restrictions on the relative magnitudes of x, and x 0 in the Green’s operators. Due to the symþ 0 0 metry of the homogeneous beam G
i ðx ; xÞ ¼ Gi ðx; x Þ ¼ 0 Gi ðx
xÞ. The inverse operator may then be written 0 0 G
1 i ðx
xÞ ¼ Gi ðx
x Þ. So, (7) may be summarised as 0 0  0 u i ðx; t  ðx
x Þ=ci Þ ¼ Gi ððx
x ÞÞui ðx ; tÞ.

ð8Þ

The functional form of the Green’s operators is obtained by applying Laplace transform techniques to (6) and solving the ordinary differential equations, which gives ~

0

0 ki ðsÞðx
x Þ ~ u u i ðx; sÞ ¼ ~ i ðx ; sÞe
1 þF ~

0 ðsci ¼~ u i ðx ; sÞe

i ðsÞÞðx
x

0Þ

ð9Þ

.

By defining the kernel ~ i ðx
x0 ; sÞ ¼ e
F~ i ðsÞðx
x0 Þ
1; G

ð10Þ

the Green’s operator may be represented by Gi ðx
x0 Þ ¼ 1 þ Gi ðx
x0 ; Þ. The kernels Gi are Green’s functions (in the sense of Krueger and Ochs [16]) and can be obtained from (10) in terms of Volterra equations of the second kind, see [15], Gi ðx
x0 ; tÞ ¼


x
x0 t

Z

t

t0 F i ðt0 ÞGi ðx
x0 ; t
t0 Þdt0
ðx
x0 ÞF i ðtÞ.

0

ð11Þ Dispersion free fields v i are obtained by using (9) 0

~


1 ðx
x0 Þ

F i ðsÞðx
x Þ 0 sci ~v u ¼~ u i ðx; sÞ ¼ ~ i ðx; sÞe i ðx ; sÞe

;

ð12Þ

which gives This means that the fields v i satisfy the dispersion free oneway equations 1  ot v ¼ 0; ci i

ð13Þ


so that the fields vþ i and vi propagate in opposite directions  at speeds ci without dispersion. The fields u i and vi are seen to be related through

v i ðx; tÞ u i ðx; tÞ

¼ ¼

Gi ððx
x ÞÞu i ðx; tÞ; 0 Gi ððx
x ÞÞv i ðx; tÞ.

v ¼ GPw;

w ¼ P
1 G
1 v;

ð15Þ

from (5) and (14). This transformation may now be used even in nonhomogeneous parts of the beam, where coupling occurs between the fields. All the usual approaches of imbedding, Green’s functions and propagators may then be utilized. However, in the present context we leave that aside and move on to another way of approaching the dispersion free wave splitting transformation. 2.2. Reduction to second order equations As in the previous section we will derive the dispersion free wave splitting in two steps. Here, this is accomplished by deriving a system of second order equations, resulting in an alternative factorization of the transform (15). Consider again Eqs. (1)–(3). Eliminating w and c, one obtains the other familiar form of the Timoshenko beam equation
2 2 2
2
2 2
2
2 4 ðo4x
ðc
2 ð16Þ 1 þ c2 Þox ot þ r0 c2 ot þ c1 c2 ot Þw ¼ 0; pffiffiffiffiffiffiffiffi where r0 ¼ I=E is the radius of gyration. Eq. (16) may be factorized and written in terms of the eigenoperators ki defined in (6),

ðo2x
k21 Þðo2x
k22 Þw ¼ 0.

ð17Þ

The differential operators in (17) could of course be factorized further, resulting in the one-way wave equations given in (6). However, we refrain from doing so since we later want to treat other cases in an analogous way, where the factorization is less obvious. The operators k2i can be represented by [15]   ð
1Þi 1 2 1 þ I 2 ð=sÞ ; ð18Þ k2i ¼ 2 o2t V; V ¼ 2r0 c2 s  ci
2 where s ¼ c2 r0 ðc
2 1
c2 Þ=2 is a characteristic time and I2 is the modified Bessel function of the first kind and order two. Introduce the fields u1(x, t) and u2(x, t) such that

ðo2x
k2i Þui ¼ 0;

i ¼ 1; 2;

w ¼ u1 þ u 2 .

ð19Þ

In analogy with Mindlin [17] and Vemula and Norris [18] we express the rotation as

 0 0 v i ðx; tÞ ¼ ui ðx ; t  ðx
x Þ=ci Þ.

ox v  i 


1 þ

T and G ¼ diagðG
1 Introducing v ¼ ðvþ 1 ; v2 ; v1 ; v2 Þ 1 ; G2 ; 0 G1 ; G2 Þ, where Gi ¼ Gi ðx
x Þ, the dispersion free wave splitting transformation becomes

0

ð14Þ

wi ¼ Ai ox ui ;

i ¼ 1; 2;

w ¼ w1 þ w2 ;

ð20Þ

where the (time domain) operators Ai are independent of the spatial coordinate. Since ui together with wi are solutions to the original Timoshenko beam equation (1)–(3) the operators become A1 ¼ Vk
2 1 ;


2
2 2
2 A2 ¼ ðc
2 1
c2 Þot k2
Vk2 ;

ð21Þ

using (1), (3), (19) and (20). Their explicit time domain expressions are given in Appendix A. Hence, the physical fields w = (w, w, c, oxw)T may be expressed in terms of u = (u1, oxu1, u2, oxu2)T

M. Johansson et al. / Computers and Structures 84 (2006) 514–527

w ¼ u1 þ u2 ;

c ¼ ð1 þ A1 Þox u1 þ ð1 þ A2 Þox u2 ;

w ¼ A 1 ox u1 þ A 2 ox u2 ;

ox w ¼ A1 k21 u1 þ A2 k22 u2 .

ð22Þ

Written on matrix form, the transformations become w ¼ S
1 u and u ¼ Sw where 1 0 1 0 1 0 C B B 0 0 A2 C A1 C B
1 S ¼B C; C B 0 1 þ A 0 1 þ A 1 2A @ A1 k21 0

1
N1

B B B S¼B B @

A2 k22

0 0

0

0


1
N2

1
N3

N1

0

0

0

N2

N3


N4

~ki ðsÞeF~ i ðsÞðx
x0 Þ 1
1 0 ~
~ki ðsÞeF i ðsÞðx
x Þ A. ~k
1 ðsÞe
F~ i ðsÞðx
x0 Þ i


~ki ðsÞe
F i ðsÞðx
x Þ 0

~

0

!

0

~

eF i ðsÞðx
x Þ

0

~

eF i ðsÞðx
x Þ ~ i ¼ 1@ B 2 e
F~ i ðsÞðx
x0 Þ

In the time domain, we have   Gþ G
1 i i B
1 ¼ ; Bi ¼ i
2
ki Gþ k G i i i

G
i Gþ i

Gþ i

0

We now derive the second transformation, which will split the fields ui into dispersion free left and right moving components. The two wave equations in (19) are treated simultaneously by considering the system     0 1 ui ox ui ¼ D i ui ; D i ¼ . ð23Þ ; ui ¼ k2i 0 ox ui Introduce the operator Bi , its inverse
T free fields vi ¼ ðvþ i ; vi Þ , such that

ð26Þ

!

k
1 i Gi ; þ k
1 i Gi

0

G
i

0

where ¼ Gi ðx
x Þ and ¼ Gi ðx
xÞ. This is in line
with the results in (14) using ui ¼ uþ and ox ui ¼ i þ ui
1 þ
ki ð
ui þ ui Þ. The operators ki are found in Appendix A. By defining v = (v1, v2)T and B ¼ diagðB1 ; B2 Þ we obtain

C 0 C C C. N4 C A

B
1 i

ð25Þ

;

ð27Þ

1

2.3. Transformation to dispersion free fields

ui ¼ B
1 i vi ;

0

~

e
F i ðsÞðx
x Þ

~
1 ¼ B i

0

The operators Ni are given in Appendix A. The operators S and S
1 may of course be derived using P and P
1 in
þ
(5), since ui ¼ uþ i þ ui and ox ui ¼ ki ð
ui þ ui Þ. The Timoshenko equation is thus reduced to two wave equations with different wave front speeds corresponding to the shear and bending modes respectively. At this point we have not yet performed the wave splitting operation but only decoupled the components of the wave field whose wave fronts propagate with different speeds.

v i ¼ B i ui ;

517

v ¼ BSw;

w ¼ S
1 B
1 v.

ð28Þ

The connection between (15) and (28) is readily seen if B
1 i and Bi are factorized as    Gi 0 I I
1 Bi ¼ ; 0 G
1
ki ki i ! !
1 I
k G
1 0 1 i i . Bi ¼ 0 Gi 2 I k
1 i Hence, the Green’s operator matrices carry the spatial dependences that deal with the dispersion, while the eigenoperator matrices are involved in the actual wave splitting. This is the reason why the original wave splitting transformation for the Timoshenko beam is found by means of a diagonalization, as the beam is left–right isotropic. A membrane or a plate with rotational symmetry lack the corresponding inward–outward symmetry. 3. The membrane

and dispersion

ox vi ¼ Ki vi ;

ð24Þ

þ
where Ki ¼ c
1 i diagf
ot ; ot g. Now, as vi and vi satisfy first order one-way wave equations they may be expressed as  0 0 v i ðx; tÞ ¼ vi ðx ; t  ðx
x Þ=ci Þ; þ where x 0 is some reference position. Using ui ¼ B
1 i;11 vi þ
1
Bi;12 vi in (23) and adopting Laplace transform techniques, the ordinary differential equations are solved as 0
1 ~ 0 ~ ~
1 ¼ a1 ðsÞeðsc
1 i þki ðsÞÞðx
x Þ þ b ðsÞeðsci
ki ðsÞÞðx
x Þ ; B 1 i;11 0
1 ~ 0 ~ ~
1 ¼ a2 ðsÞe
ðsc
1 i þki ðsÞÞðx
x Þ þ b ðsÞe
ðsci
ki ðsÞÞðx
x Þ . B 2 i;12


As vþ i propagates in the positive x-direction while vi propagate in the negative x-direction, we should have a1(s) = a2(s) = 0. If we set b1(s) = b2(s) = 1, the same normalization is used as in Section 2.1. It is straightforward to calcu~
1 and B ~
1 using (6). Then, late B i;21 i;22

In this section the dispersion free transformation is derived for radial waves in a homogeneous membrane. The objective is to obtain radially in- and outgoing waves that are uncoupled and propagate without dispersion. Since we have cylindrical geometry we will encounter geometrical dispersion. So, the transformation is performed in line with the methods presented in Section 2.3. The wave equation for the membrane is r2 wðr; h; tÞ
c
2 o2t wðr; h; tÞ ¼ 0;

ð29Þ 2

where w(r, h, t) is the deflection, pffiffiffiffiffiffiffiffiffi $ is the two-dimensional Laplace operator, c ¼ T =q is the wave speed and q, T are mass per unit area and tension per unit length respectively. In order to make comparisons to the plate in Section 4 more apparent, the membrane equation is written ðr2
k2 Þw ¼ 0;

k2 ¼ c
2 o2t .

ð30Þ

As the aim is to find a transformation for radial waves, w is expanded according to

518

M. Johansson et al. / Computers and Structures 84 (2006) 514–527

wðr;h;tÞ ¼ we0 ðr;tÞ þ

1 X ðwem ðr;tÞcosðmhÞ þ wom ðr;tÞ sinðmhÞÞ. m¼1

ð31Þ For each Fourier component of the physical variables we may write (30) as    e;o  0 1 wm 2 ¼ ; w . or wm ¼ Dm wm ; Dm ¼ m k2 þ mr2
1r or we;o m ð32Þ The parity index e,o is dropped hereafter, since the even and odd variables satisfy the same equations. Introduce the transformation operator Bm , its inverse B
1 m and dis
T persion free fields vm ¼ ðvþ ; v Þ , such that m m vm ¼ Bm wm ; wm ¼ B
1 v or vm ¼ Kvm ; m m;
where K = c
1diag{
ot, ot}. Since vþ m and vm satisfy first order one-way wave equations they may be expressed as  0 0 v m ðr; tÞ ¼ vm ðr ; t  ðr
r Þ=cÞ;

ð33Þ þ B
1 m;11 vm þ

0

where r is some reference radius. Inserting wm ¼
B
1 m;12 vm in (32) and using Laplace transform techniques, the solutions of the ordinary differential equations are obtained as ~ ðsc
1 ðr
r0 ÞÞ þ b1 ðsÞK m ðkrÞe ~ ðsc
1 ðr
r0 ÞÞ ; ~
1 ¼ a1 ðsÞI m ðkrÞe B m;11


1 0
1 0 ~
1 ¼ a2 ðsÞI m ð~ B krÞe
ðsc ðr
r ÞÞ þ b2 ðsÞK m ð~ krÞe
ðsc ðr
r ÞÞ ; m;12

ð34Þ where Im and Km are the modified Bessel functions of the first and second kind respectively. Reasonable physical properties of w are that the incoming fields, due to v
m, are regular at the origin, while the outgoing fields, due to vþ m , are regular at infinity. This is accomplished by setting a1(s) = b2(s) = 0. If we choose a2(s) = b1(s) = 1 and set r 0 = 0 for all the fields, the Laplace transforms of the operators B
1 m and Bm are ! sr
src ~ ~ c ð krÞe I ð krÞe K
1 m m ~ ¼ B ; sr sr m ~ kK 0m ð~ kI 0m ð~ krÞe c ~ krÞe
c ! 0 ~
src ~ ~
sr kI ð krÞe
I m ðkrÞe c m ~m ¼r . ð35Þ B sr sr
~ kK 0m ð~ krÞe c K m ð~ krÞe c
1 In the time domain B
1 m ¼ Bm ðr; Þ and Bm ¼ Bm ðr; Þ. The kernel matrices may be written as [19,20]  þ   þ
 L Um
U
Um U
m m B
1 ¼ ¼ r ; B ; m m
L
Uþ Uþ L
Uþ Lþ U
m m m m

ð36Þ ±


1

where L = or ± c ot. Using t 0 = ct/r, [19,20] 8 T m ðt0 þ 1Þ
U m

are given by

where Tm are the Chebyshev polynomials of the first kind. 0 Note that the kernels U m are singular at t = 0 and that the
0 Um are singular at t = 2 and have compact support. There are reasons to comment on the fields v m . Suppose a line source is located at a radius r = r1, say. The transformations from the physical fields are such that this source þ will solely cause fields v
m in r < r1 and vm in r > r1. Reflec
tions from the origin are part of vm for r < r1. Hence, the fields v m are not in- and outgoing fields in the usual sense. Moreover, vþ 0 ð0; tÞ may be interpreted as the magnitude of a point force at the origin, as K 0 ðsrcÞ is seen to be the Green’s function for the modified Helmholtz equation. The present normalization of the transformation is such that v
m ð0; tÞ ¼ wm ð0; tÞ, in the case of no sources at the origin. There are of course other normalization alternatives that may be a better choice when considering numerical implementation. For the calculations in Section 5.4 we make the pffiffi choice a2 ðsÞ ¼ b1 ðsÞ ¼ s. 4. The Mindlin plate As for the Timoshenko beam, the derivation of the dispersion free transformation will take place in steps. The first step is to decouple the equations obtaining three second order partial differential equations, which is in line with Section 2.2. Up to this point, we do not specify any particular coordinate system in the plane of the plate. However, for the wave splitting operation, the cylindrical geometry is chosen in order to obtain dispersion free radial waves. Hence, the dependent variables are expanded in the angular direction. From here, the derivation of the dispersion free fields is performed in a similar manner as for the membrane, Section 3. The Mindlin plate equation for a plate of thickness h and with density q may be written [17] D qh3 2 fð1
mÞr2 w þ ð1 þ mÞrr  wg
j2 Ghc ¼ o w; ð38Þ 2 12 t 2 ð39Þ j2 Ghr  c ¼ qhot w; c ¼ rw þ w;

ð40Þ

where D = Eh3/(12(1
m2)), E is Young’s modulus, G is the shear modulus, m is Poisson’s ratio and j2 is the shear coefficient. The displacement w is transverse to the central plane of the plate, that is, in the z-direction. w is the twodimensional vector of rotations and c is the two-dimensional vector of shearing angles. 4.1. Reduction to second order equations In analogy with [17,18], we make the general ansatz

ð37Þ

w ¼ A rw
ez  ru3 ; ð41Þ where A is an operator which is independent of the spatial coordinates. To separate w and u3, we first take the curl of (38), using (38)–(41), which reduces to ðr2
k23 Þu3 ¼ 0;

2 2 2 k23 ¼ c
2 3 ot þ j =r0 ;

ð42Þ

M. Johansson et al. / Computers and Structures 84 (2006) 514–527

pffiffiffiffiffi where r0 ¼ h= 12 is the radius of gyration of the plate pffiffiffiffiffiffiffiffiffi cross section and c3 ¼ G=q is the wave front speed. The equation for w is obtained by taking the divergence of Eq. (38), using (38)–(41), which turns into the often used equation for the plate
2 2 2
2
2 2
2
2 4 ðr2 r2
ðc
2 1 þ c2 Þot r þ r0 c2 ot þ c1 c2 ot Þw ¼ 0.

ð43Þ Hence, (43) is the two-dimensional counterpart of the beam pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi equation (16), where c1 = jc3 and c2 ¼ E=ðqð1
m2 ÞÞ are the velocities of the shearing and bending modes respectively. In analogy with (17), Eq. (43) can be factorized as ðr2
k21 Þðr2
k22 Þw ¼ 0;

where the operators are obtained from (18), but with c1, c2 and r0 for the plate. Introduce new fields u1 and u2 which, in accordance with (19), gives the set of equations ðr2
k2i Þui ¼ 0;

i ¼ 1; 2; 3.

ð45Þ

By using (40)–(42) and (45) the new fields are related to the physical fields through relations similar to (22) w ¼ u1 þ u 2 ; w ¼ A1 ru1 þ A2 ru2
ez  ru3 ; c ¼ ð1 þ A1 Þru1 þ ð1 þ A2 Þru2
ez  ru3 ; r  w ¼ A1 k21 u1 þ A1 k22 u2 ; r  w ¼
k23 u3 ez .

0

1

direction. In the present case, the fields are expanded in the angular direction. It is clear from (46) that the fields may be collected according to g = {w, u1, u2, wr, cr} and t = {wh, ch, u3}. For each term, the fields may be expanded through gðr; h; tÞ ¼ ge0 ðr; tÞ þ

ðgem ðr; tÞ cosðmhÞ þ gom ðr; tÞ sinðmhÞÞ;

m¼1

ð47Þ tðr; h; tÞ ¼ to0 ðr; tÞ þ

1 X

ðtom ðr; tÞ cosðmhÞ
tem ðr; tÞ sinðmhÞÞ.

m¼1

0

1

ð46Þ

ð48Þ The even and odd components decouple, so the parity index e, o will be suppressed in the following. Expansions (47) and (48) enable us to eliminate, from the nine first order variables, the ones which are differentiated with respect to the angle h. Thus, all in all we have six dependent variables for each m and parity. Define the vectors of the physical variables wm = (wm, wr,m, cr,m, orwr,m, wh,m, orwh,m)T and the new variables um = (u1,m, oru1,m, u2,m, oru2,m, u3,m, oru3,m)T. By (46), these fields are related through um ¼ S m w m ;


rm2 M2

0

0

wm ¼ S
1 m um ;

ð49Þ

where

0

0

B 0 A1 0 A2
mr B B B 0 1 þ A1 0 1 þ A2
mr B S
1 ¼ B 2 2 m m B A1 ðk21 þ mr2 Þ
1r A1 A2 ðk22 þ mr2 Þ
1r A2 r2 B m m B A1 0 A2 0 0 @ r r 2 m m
rm2 A1 A1
rm2 A2 A2
ðk23 þ mr2 Þ r r 0 m
1r N4 0
N4 N4 1
N1 r B m2 m 0
1
N2 þ r2 M4 1
N3 0
r2 M4 B B 1 B N1 N4 0 N4
mr N4 B r Sm ¼ B 2 m B 0 N2
mr2 M4 N3 0 M4 r2 B B m 1 0 M3 0 0
r M3 @ r m M1 r

1 X

ð44Þ

k2i

519


mr M2

4.2. Reduction of dependent variables The aim is to reduce the equation system of nine dependent variables to a one-dimensional problem. Folkow [21] expressed the state of a plate in a Cartesian coordinate system where the number of dependent variables was reduced to six after performing a Fourier transform in one spatial

2

1 þ mr2 M2

1

0 C C C 0 C C ; m C C r C 1 C A 1
r 0

1

C
mr M4 C C 0 C C C. m M 4 C r C
M3 C A 0

The operators Ni and Ai are obtained from the beam case in Section 2.2 but with c1, c2 and r0 for the plate. The Mi are given in Appendix A. Note that u3,0 decouples from u1,0 and u2,0, in accordance with [21]. The same kind of decoupling also occurs when r ! 1, in which case u1,m and u2,m transform as in the Timoshenko case presented in Section 2.2.

520

M. Johansson et al. / Computers and Structures 84 (2006) 514–527

So far the Fourier components of the physical fields wm have been transformed to a new set of variables um, which obey three decoupled second order wave equations with different wave front speeds. We now perform the wave splitting on each of these wave equations, giving six oneway wave equations, with the additional condition that the one-way wave equations shall be free of dispersion. 4.3. Transformation to dispersion free fields The transformation to dispersion free fields is performed in the same manner as for the membrane in Section 3. For each Fourier component, Eq. (45) may be written ! 0 1 or ui;m ¼ Di;m ui;m ; Di;m ¼ ; 2 k2i þ mr2
1r ð50Þ   ui;m . ui;m ¼ or ui;m Introduce the transformation operator Bi;m , its inverse þ
T B
1 i;m and dispersion free fields vi;m ¼ ðvi;m ; vi;m Þ such that vi;m ¼ Bi;m ui;m ;

ui;m ¼

B
1 i;m vi;m ;

or vi;m ¼ Ki vi;m ;

ð51Þ

c
1 i diagf
ot ; ot g.

The following steps are idenwhere Ki ¼ tical to the membrane case and Eqs. (33)–(36) holds for the fields ui,m, i = 1, 2, 3, by adding the index i to c; k; a1 ; a2 ;   b1 ; b2 ; Bm ; B
1 m ; L ; Um . The coefficients ai,1 . . . bi,2 are chosen on the same basis as for the membrane noting that the leading terms of the ki are of the same order as in the membrane case. Hence, ai,1 = bi,2 = 0 and ai,2 = bi,1 = 1. The inverse Laplace transforms of these elements are not as easy as in the case of the membrane. There is, however, at least one case where an explicit expression for the elements of the matrices exist, namely when i = 3 and m = 0 [20]. Then, 8  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > t0 ð2 þ t0 Þ < c3 cospaffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 < t0 ; þ 0 ð2 þ t 0 Þ U3;0 ¼ r t > : 0 t0 < 0; 0  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c3 cosh a t0 ð2
t0 Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t0 2 ð0; 2Þ; B U
ð52Þ rp t0 ð2
t0 Þ 3;0 ¼ @ 0

0 t 62 ð0; 2Þ; where a = jr/r0 and t 0 = c3t/r. In conformity with the 0 membrane case, the kernels U 3;0 are singular at t = 0 and
0 U3;0 is singular at t = 2 and has compact support. Note that U
3;0 grows exponentially with increasing radius, in contradistinction to the case of the u3-mode in Cartesian coordinates [21]. For the other cases, the behavior in the time domain for small times can be obtained approximately by studying the asymptotics of the elements of ~ i;m and B ~
1 for large values of the Laplace variable, see B i;m Appendix B. From the asymptotic expansions it is indicated that all the kernels U
i;m have compact support. þ þ


T With vm ¼ ðvþ ; v ; v ; 1;m 1;m 2;m v2;m ; v3;m ; v3;m Þ and B ¼ diag  ðB1;m ; B2;m ; B3;m Þ, the combined transform is vm ¼ Bm Sm wm ;


1 wm ¼ S
1 m B m vm .

ð53Þ

5. Application to an inhomogeneous membrane 5.1. The dynamic equations We now consider the case when the mass per unit area of the membrane varies as a function of the radius. Hence, the wave speed will be a function of r. We use the same transformation matrices Bm and B
1 m as in Section 3, however with a varying velocity c(r). Insertion of wm ¼ B
1 m vm into Eq. (32) gives
1 or vm ¼ ðBDm B
1 m
Bor Bm Þvm .

In order to use the results from the homogeneous membrane in Section 3, we introduce the following decomposition
1
1 0 or B
1 m ¼ or jcðrÞ¼const. Bm þ c ðrÞocðrÞ Bm .

ð54Þ

Hence, the dynamics of the split fields becomes or vm ¼ ðK þ Lm Þvm ;

K ¼ cðrÞ
1 diagð
ot ; ot Þ;

Lm ¼
c0 ðrÞBm ocðrÞ B
1 m .

ð55Þ

Thus, at each point we perform the wave splitting with respect to the local wave speed. In the Laplace domain, the operator elements of Lm are 0

~ m;11 ¼ c f2 ðI m ðfÞ
I 0 ðfÞÞðK m ðfÞ þ K 0 ðfÞÞ
m2 I m ðfÞK m ðfÞ; L m m c ð56Þ 0 ~ m;12 ¼ c f2 ðI m
1 ðfÞI mþ1 ðfÞ
I 2 ðfÞÞe
2f ; ð57Þ L m c 0 ~ m;21 ¼ c f2 ð
K m
1 ðfÞK mþ1 ðfÞ þ K 2 ðfÞÞe2f ; L ð58Þ m c ~ m;22 ¼
L ~ m;11 ; L ð59Þ where c = c(r), c 0 = c 0 (r) and f = sr/c. The elements Lm,ij(r, t) can be obtained in the time domain by taking the inverse Laplace transform of each modified Bessel function separately, applying the convolution theorem and finally differentiate twice with respect to time. The dynamics for the split fields (55) become
1 þ þ
or v þ m ¼
c ot vm þ Lm;11  vm þ Lm;12  vm ;
1
þ
or v
m ¼ c ot vm þ Lm;21  vm
Lm;11  vm ;

ð60Þ

þ where Lm;11  vþ m ¼ Lm;11 ðr; Þ  vm ðr; Þ etc. Hence, for an inhomogeneous region the split fields satisfy coupled partial integro-differential equations. Using these equations one can exploit various time domain methods to solve direct and inverse problems, such as the imbedding method, Green’s functions or the propagator technique. In the example considered here we utilize the imbedding method. In Section 5.3 explicit expressions for the Lm,ij-kernels are given in the case m = 0.

5.2. The imbedding equation We seek an equation for the operator that, for a certain r in the inhomogeneous region, relates the internal fields to

M. Johansson et al. / Computers and Structures 84 (2006) 514–527

each other. Due to linearity, causality and time translation invariance this relation has the form
vþ m ðr; tÞ ¼ Rm ðr; Þ  vm ðr; Þ.

ð61Þ

The kernel Rm(r, t) will be referred to as the reflection kernel as usual, although v m are not physical wave fields travelling inwards and outwards in the usual sense. To obtain an equation for Rm, vþ m is eliminated from (60) by insertion of (61). After use of both equations, partial integration and the causality of Rm and v
m , the equation for Rm(r, t) is obtained as   2 or þ ot Rm
2Lm;11  Rm þ Rm  Lm;21  Rm
Lm;12 ¼ 0. c ð62Þ In the following the analysis will be restricted to the case m = 0, so the index m is suppressed hereafter. It is also assumed that c(r) is continuous, while c 0 (r) is piecewise continuous. 5.3. Certain known components of the reflection kernel for the case m = 0 In the time domain the kernels elements of L0 become, after lengthy calculations, ( 0

L11 ¼ c b11

Eðd 211 Þ
ð1
d 211 ÞKðd 211 Þ;

ct < 2r;

d 11 Eðd
2 11 Þ;

ct > 2r; 8
1 dðctÞ þ b12 fð4r
ctÞEðd 12 Þ þ 2ð2r
ctÞKðd
1 > 12 Þg; > > 0< ct < 2r;
c L12 ¼ 2p > dð4r
ctÞ þ b12 H ð4r
ctÞfctEðd 12 Þ þ 2ð2r
ctÞKðd 12 Þg; > > : ct > 2r;
pc0
1 fdðctÞ þ b21 fð4r þ ctÞEðd
1 L21 ¼ 21 Þ
2ð2r þ ctÞKðd 21 Þgg; ð63Þ 2

where H is the Heavyside step function, d is the Dirac delta function while K and E are the complete elliptic integrals of the first and second kind, respectively. The bij and dij are given by b11 ¼ b12 ¼ b21 ¼ d 11

8r3 ; pc2 t2 ð4r2
c2 t2 Þ 8r2 pctð4r
ctÞð2r
ctÞ

2

;

8r2

; 2 pctð4r þ ctÞð2r þ ctÞ  2 ct 4r
ct ¼ ; d 12 ¼ ; 2r ct

d 21

 2 4r þ ct ¼ . ct

The kernels Lij of the coupling operator involve terms of E and K with algebraic coefficients and arguments. The coefficients are generally singular. In addition, K has a logarithmic singularity when the argument approaches unity. As the kernels also contain step functions and delta functions, the reflection kernel R(r, t) (obtained through (62)) may contain these various types of singularities as well. Billger

521

and Folkow [2] give a proof which, in effect, states that multiplicative and convolutional parts of (62) must be satisfied independently. For example, the coefficients for collected delta functions with a specific argument appearing in (62) must then be zero. It may also be shown that coefficients for power and logarithmic singularities must be zero. This allows for analytical determination of some of the singularities and distributions of the reflection kernel. For continuous wave speed c(r), the reflection kernel R(r, t) contains no delta functions. Hence, the reflection kernel can be decomposed according to X ½Ri ðrÞH ðt
bi ðrÞÞ Rðr; tÞ ¼ R1 ðr; tÞ þ X X
1 þ fi ðrÞðt
ci ðrÞÞ þ gi ðrÞ lnðjt
ci ðrÞjÞ X þ hi ðrÞðt
ci ðrÞÞ lnðjt
ci ðrÞjÞ. ð64Þ Here R1(r, t) is the bounded and smooth part of R(r, t), [R]i(r) is the jump discontinuities across t = bi(r), that is [R]i(r) = R(r, bi(r)+)
R(r, bi(r)
). The functions bi(r), ci(r), fi(r), gi(r) and hi(r) are yet to be determined. The jump discontinuities give rise to delta functions when inserted into (62), but so do also convolutions of power singularities. Consider for instance [22] Z t dt0 0 0 0 ðt
t 1 Þðt
t
t 0 Þ   ðt
t0 Þðt
t1 Þ 1 
p2 dðt
t0
t1 Þ; ð65Þ ¼ ln   t
t0
t1 t0 t1 where t0, t1 > 0. The initial jump discontinuity at t = 0 can be determined directly. As the kernel L12 appearing in (62) contains d(ct), R has a jump discontinuity [R]0 given by ½R0 ¼
c0 =4p;

b0 ðrÞ ¼ 0.

ð66Þ

Before studying other jump discontinuities, consider some of the power singularities. By studying Lij in the neighborhood of the singularities, i.e. by examining expansions about singular points, the coefficients of the various singular parts can be identified. These can be compared to the corresponding terms due to R, obtained by inserting R from (64) into the integro-differential equation (62). By collecting singular terms of the same order propagating along specific characteristics, some of the unknown functions appearing in (64) can be determined directly. For power singularities of order one the convolution integrals in (61) and (62) are interpreted as CPV integrals since they are valid for all t P 0. A power singularity of order two is found in L12. The only term that can cancel this term comes from the principal part of (62). Thus,
2


2


f0 ðt
c0 Þ ð2=c
c00 Þ
ð
2c0 rp
2 ðct
2rÞ Þ ¼ 0; which gives c0(r) = 2r/c(r) and f0(r) = f0 = p
2. Investigation of the other terms in (62) shows that there are no further singularities in R(r, t) along that characteristic, g0(r) = h0(r) = 0, but from inspection of the principal part it is seen that there may be power singularities propagating along t = ci>0(r), where

522

ci>0 ðrÞ ¼

M. Johansson et al. / Computers and Structures 84 (2006) 514–527

Z 0

r

2 dr0 þ C ci ; cðr0 Þ

ð67Þ

and C ci are to be determined. To determine these, consider a case where c is constant for r < r 0 and a function of r on the outside. Then R(r < r 0 , t) = 0. From continuity considerations for v± and (61) it follows that there must be a power singularity f1(r 0 ) =
f0, that propagates from the point (r, t) = (r 0 , c0(r 0 )) along the characteristic t = c1(r), where c1(r) is given by (67) with C c1 ¼ 0, see II in Fig. 1. For continuous c, there are no further power singularities. From the other terms in (62) we may conclude also that f1(r) =
f0 and that g1(r) satisfy (2p)2org1 = (c 0 )2/c
c 0 /r. Since (62) is non-linear in R and the kernels Lij and R are singular, there are also delta functions propagating along t = 2c0, 2c1, c0 + c1. It turns out that when each of these non-differentiated terms in (62) are collected, including d(4r
ct) from L12, the only one with a non-vanishing coefficient propagates along t = 2c1(r) and has the coefficient c 0 /(2pc). Hence, another jump discontinuity appears in R along t = b1(r) = 2c1(r), which is given by ½R1 ð2=c
2c01 Þ þ c0 =ð2pcÞ ¼ 0. Thus, [R]1 =
[R]0. There may be additional jump discontinuities [R]i>1 propagating along t = bi>1(r), where bi>1(r) is given by the right hand side of (67). These are initiated where c 0 has jump discontinuities, see I and III in Fig. 1. From considerations based on the continuity of v± and Eq. (62), initial values and differential equations for these additional jump discontinuities may be derived. The additional logarithmic singularities that these jump discontinuities lead to cancel out in (62). Hence, they do not imply further singularities in R. There are, however, additional singularities in individual terms in (62) and they run along characteristics such as t = c1 + bi>0, c2 + bi>0. Finally, we note that c1(r) is the true travel time from radius r to the origin and back, while c0(r) is the travel time with the constant wave speed c(r) everywhere, i.e. the travel time for the homogeneous case.

Fig. 1. Characteristics for an inhomogeneous ring with c 0 (r) discontinuous at the boundaries. Note that [R]0 runs along t = 0 for r 0 < r < r00 .

5.4. Numerical example In this section we present numerical results from simulated direct scattering problems for inhomogeneous membranes. All results concern the case with circular symmetry; w(r, h, t) = w0(r, t), see Section 3. We will first give an outline of the numerical implementation of the radial wave splitting and the computation of the reflection kernel. The direct problem is then considered for three wave speed profiles, where the result is compared to a FE solution. 5.4.1. Numerical implementation of the wave splitting and the reflection kernel As was mentioned at the end of Section 3, one may consider other normalizations of the transforms than the one presented pffiffi there. Here, we make the choice a2 ðsÞ ¼ b1 ðsÞ ¼ s. This makes the transform operators more suited for numerical implementation in the time domain, see Appendix A. The kernels are still singular but integrable, interpreted as CPV integrals. This re-normalization does not influence the coupling terms (63), and hence not the dynamics of the split fields, (60). As a consequence, the reflection operator is not altered. In Section 5.3 it was demonstrated how parts of the reflection kernel could be determined analytically. The remaining part of R may then be calculated numerically with simple methods, taking into consideration the known singularities. Here, only a subset of the analytical information in Section 5.3 has been utilized, namely the information about the power singularities travelling along t = c0 and t = c1 and the jump discontinuity [R]0. The rest of the jumps and singularities are being part of the numerical solution, and have thus been handled with less accuracy. The reason for this is that in many scattering situations, the main objective is to be able to determine the fields up to and past one travel time through the inhomogeneous medium and back. This is for example the case when studying inverse problems, where information about the varying media is to be recovered from reflection data. The reflection kernel is therefore divided into a regular part, Rreg, and a singular part, whereby we mean that the regular part does not contain power singularities and the jump discontinuity [R]0. The regular part is then calculated numerically by means of finite differences. Since all terms in the imbedding equation contain singularities it is inevitable to have to perform convolutions with singular kernels numerically. This problem has been solved by replacing all singularities with discrete functions that give an approximate result when being integrated, see Appendix C. Travel time coordinates are introduced to simplify the calculations. If the membrane is homogeneous outside a radius r = r1 say, the non-dimensional travel time coordinates are Z Z r1 1 r dr0 dr0 ; g ¼ t=2l; l ¼ . ð68Þ nðrÞ ¼ l 0 cðr0 Þ cðr0 Þ 0

M. Johansson et al. / Computers and Structures 84 (2006) 514–527

This makes the characteristic equation t = c1(r) become g = n in non-dimensional space. By introducing Lmod as 12 L12 without the power singularity of order two and the first delta function at t = 0, (62) may be written ðon þ og ÞRreg
2L11  R þ R  L21  R
Lmod 12 ¼ 0.

ð69Þ

Then, (69) does not contain the power singularity of order two and no delta functions for t < 2 min(c0, c1). The inhomogeneous region is then divided into N subintervals and Rreg(1, g) is solved numerically from knowledge of the coupling terms and R(n, 0), which is known from the jump condition (66). The integrals are all approximated with the trapezoidal rule, where the logarithmic singularities and the power singularities of order one are represented according to Appendix C. 5.4.2. Comparisons with FE solutions Solutions obtained by means of the commercial software FEMLAB are used for comparisons. The membranes that are considered are of radius r = r2, with a possible region of inhomogeneity from the origin to r = r1 < r2. In all the examples discussed below, the outer radii r2 are large enough so that no reflections from these boundaries have to be considered. The membrane is initially at rest, and waves are generated at the boundaries according to ( sin2 ð5tÞ; t 2 ð0; p=5Þ; wðr2 ; tÞ ¼ ð70Þ 0; t 62 ð0; p=5Þ. We start by noting that the properties of the wave splitting are such that no v+-waves are generated when boundary condition (70) is applied to a homogeneous membrane. The v
-waves may then be computed solely from knowledge of w. As v
propagates without dispersion, it is straightforward to obtain w(r, t) and orw(r, t) everywhere in the membrane using the transformations. When compared to the solutions using the FE solver, the fields w and wFE agree very well. However, the w-solutions were obtained much quicker for a certain accuracy. For the inhomogeneous membranes, the objective is to obtain the fields at the outer boundaries of the inhomogeneous regions, r = r1. Using the FE solution, vþ FE ðr1 ; tÞ and v
ðr ; tÞ are transformed from w (r , t) and its spatial 1 FE 1 FE derivative. v+(r1, t) may also be computed from v
ðr FE 1 ; tÞ by means of (61) and is denoted vþ ðr ; tÞ. Comparisons of vþ 1 R R þ and vFE are made, where the time is measured from when the incoming waves reach r = r1. Note that it is possible to calculate v
(r2, t) over a finite time interval from knowledge of only w(r2, t), using v+(r2, t) = 0 as in the homogeneous case. This solution, which is valid up until the reflected fields from the inhomogeneity reach r = r2, may then be used at r = r1 instead of
v
FE ðr 1 ; tÞ. So, the fields v (r1, t) are the same no matter the inhomogeneity, over a finite time interval. This means that v
for the homogeneous case can be used in the inhomogeneous cases over this finite time interval.

523

2 1.5 1 0.5 0 –0.5 –1 –1.5

0

0.5

1

1.5

2

2.5

3

3.5

Fig. 2. Different components of w at r = r1: (dash-dotted) wSING, (dashed) þ w
FE , (solid) wFE, (dotted) wFE .

5.4.2.1. Interior with a lower speed of propagation. Here the wave speed profile has a linear variation, where c = 0.5 for r < 0.5 and c = 1 for r > 0.7 = r1. We first take a look at the different components of the wave field in Fig. 2. By setting vþ FE ¼ 0 and using the wave splitting transformation we obtain the wave field for a homogeneous membrane, w
FE . Similarly, by setting þ v
FE ¼ 0 we obtain the wave field wFE , which is the difference between the physical fields for a homogeneous membrane and an inhomogeneous membrane. Of course, wFE ¼
wþ FE þ wFE , which is also indicated in Fig. 2. A rough estimation of the wave field in an inhomogeneous membrane is obtained by approximating the reflection kernel with the known power singularities. The compensation for the reflection at the origin is then taken care of while the reflections at the inhomogeneities, which are small in comparison to the total field, are neglected. In Fig. 2 this field is denoted by wSING, and it is seen to agree remarkably well. We now turn to the comparisons presented in Fig. 3, þ þ where the differences vþ FE
vR are compared to 1% of vFE for different number of subintervals N of the inhomogeneous region. As is seen in Fig. 3(a), the error is small up to approximately t = 3. For this particular wave speed profile we have that c0(0.7) = 1.4 and c1(0.7) 2.55. The influences from these singularities are seen clearly for N = 100. The increase in error around t = 3 is probably due to distributions propagating along t = 2c0(r) = 2.8. The deviation between t 0.6 and t 1.1 stem from the jump discontinuity, which propagates along t = b2. The convergence is not that clear in this region since it depends on both the location of the jump discontinuity as well as the approximation of it. It is rather the absolute value of the differences between successive approximations that converges to zero. The structure of the error for other times is due to convolutions with power singularities and the transformation from the physical fields to the split fields. 5.4.2.2. Interior with a higher speed of propagation. Here, the wave speed profile has a linear variation where c = 2

524

M. Johansson et al. / Computers and Structures 84 (2006) 514–527 0.01

0.01

0.005

0.005

0

0

–0.005

–0.005

–0.01

0

0.5

1

1.5

2

2.5

–0.01

3

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

(b)

(a)

þ þ Fig. 3. The difference vþ FE
vR for slow and fast interior at r = r1: (dash-dotted) N = 100, (dashed) N = 200, (solid) N = 400, (dotted) 1% of vFE . (a) Slow interior, (b) fast interior.

5

x 10–3

0.04 0.02

2.5 0 0

–0.02 –0.04

–2.5 –0.06 –5

–0.08 0

0.5

1

1.5

2

2.5

(a)

0

1

2

3

4

5

6

(b)

þ Fig. 4. The difference vþ FE
vR and the regular part of R for stratified interior at r = r1. In (a), (dash-dotted) N = 100, (dashed) N = 200, (solid) N = 400, (dotted) 1% of vþ . (a) Stratified interior, (b) the regular part of R. FE

for r < 0.5 and c = 1 for r > 0.7 = r1. We then have that c0(0.7) = 1.4 and c1(0.7) 0.77. As is seen in Fig. 3(b), the error is small up to approximately t = 1.4. The time of the first influence from distributions propagating along t = 2ci(r) and t = c0(r) + c1(r) is in this case t 1.27. 5.4.2.3. Stratified interior. In this case r1 = 1 and for r < r1 there are several layers where the wave speed vary linearly between the radii r = 0, 0.2, 0.5, 0.7, 1. The wave speeds at these radii are c = 1.25, 1.25, 1.2, 1.1, 1 respectively. At the boundary of the inhomogeneous region c0(1) = 2 and c1(1) 1.73. From Fig. 4(a) it is seen that the error is small for times up to approximately t = 2. The first influence from the ‘‘later’’ characteristics occur for t 2.05. The regular part of the reflection kernel, including the logarithmic singularity propagating along t = c1 1.73, is shown in Fig. 4(b) for the stratified case. Here, some of the numerical difficulties are apparent around time t = 2.5.

6. Conclusions This paper deals with wave splittings that transform to one-way wave equations. For the radial membrane and the radial Mindlin plate, which both exhibit directly geometric dispersion, the splitting was obtained using non-standard methods. This resulted in dispersion free split fields in a straightforward way. This feature could also be exploited for other structures that do not require the presented method to obtain a wave splitting, such as the Timoshenko beam. The main advantage of deriving dispersion free wave splittings is that when the split fields at some point are determined, translation through homogeneous regions is trivial. This property can be utilized when solving inverse problems since one then do not want to waste computational resources on homogeneous regions, especially when the fields are measured far from the inhomogeneous region. But, as is demonstrated in the case of the Timoshenko beam, this is done to the price of a wave splitting transformation that is more resource demanding than the original one.

M. Johansson et al. / Computers and Structures 84 (2006) 514–527

The method agrees with the extension of a previous wave splitting for the Timoshenko beam. For the radial membrane the transformation is a decomposition into the field from a source at the origin and the regular part of the physical field at the origin. (This decomposition is in analogy with, for instance, the T matrix method. See, e.g., Ref. [23] which incidentally gives the first time domain formulation of the T matrix method.) By combining these two cases, it is indicated how the wave splitting transformation for the Mindlin plate may be obtained. To validate the method for inhomogeneous media the inhomogeneous membrane has been investigated by means of the invariant imbedding technique for a direct problem. The results from a numerical implementation has been compared to solutions from FE computations, and agreement was obtained within time intervals relevant for solving inverse problems. The corresponding problem for the Mindlin plate can be expected to be more complex due to the existence of several wave front speeds and lack of analytical representations of certain operators. Note, however, that the major difficulties in the case of the inhomogeneous membrane were due to the varying wave speed resulting from the type of inhomogeneity. A thickness variation in the Mindlin plate, for instance, would not change the wave front speeds. Note also, that the difficulties of having several wave front speeds have successfully been handled for the case of Timoshenko beams [2–4,7]. Acknowledgements This work has been supported by the Swedish Research Council for Engineering Sciences (TFR) and the Swedish Research Council (VR). This is gratefully acknowledged. Appendix A. Operators The explicit representation of the operators Ai are   c2 c21 2 þ I 2 ð=sÞ sinðc1  =r0 Þ; A1 ¼ 2c1 s c22   2    c2 c2 c22 2 þ I 2 ð=sÞ sinðc1  =r0 Þ  . A2 ¼ 2
1
c1 2c1 s c21  The operators Ni can be given explicitly using the operators C1 ¼ C 1 ðÞ; C2 ¼ C1 ð1
C1 Þ
1 ¼ C 2 ðÞ; Z t=s 1 I 1 ðnÞ dn 1 ; C 2 ðtÞ ¼ I 1 ðt=sÞ. C 1 ðtÞ ¼ 2s 0 n 2s Then, with q ¼ ðc22 þ c21 Þ=ðc22
c21 Þ > 1, the Ni are q
1 þ qC2 ; N1 ¼
C2 ; N2 ¼ 2 q
1 2 C2 þ C1 C2 ; N3 ¼ qþ1 qþ1 q
1
2 ot ð1 þ 2C2 Þ. N4 ¼ c22 2

525

are The operators k
1 i k
1 1 ¼ c1 Q1  þF 2  Q2 ; Q1 ðtÞ ¼ J 0 ðc1 t=r0 Þ;

k
1 2 ¼ c2 Q1  þF 1  Q2 ;

Q2 ðtÞ ¼ c1 c2 o
1 t Q1 ðtÞ.

The Mi are given by Mi ¼ M i ðÞ, where c22 r20 M 1 ðtÞ; c21 Z t=s c2 r0 M 4 ðtÞ ¼ 1 2 I 0 ðnÞ dn. 2c2 j 0

M 1 ðtÞ ¼

c1 sinðc1 t=r0 Þ; r0

M 3 ðtÞ ¼

r20 M 1 ðtÞ; j2

M 2 ðtÞ ¼

The multiplication by s±1/2 to the wave splitting transfor~ 0;ij and B ~
1 given in Section 3 is equivamation kernels B 0;ij lent to a 1/2 order time integral [24]. After some simplifications, the new transform kernels Bsij are given for the case m = 0 with Bs0;ij ¼ Bsij , 8 pffiffiffiffi pffiffiffi  c0 c > > 4p0 dðtÞ þ 2p3=2 t ð1
c10 =tÞ Eðt=c0 Þ > > >  > > < þðEðt=c Þ
Kðt=c ÞÞ ; t < c0 ; 0 0 Bs11 ¼   > > > 1 > 3=21 pffi ð1
t=c Eðc =tÞ þ 2ðEðc =tÞ
Kðc =tÞÞ ; > 0 0 0 t 2p > 0Þ > : t > c0 ; ( cpffiffiffi c
p3=20 Kðt=c0 Þ; t < c0 ; Bs12 ¼ cc0 pffi Kðc =tÞ; t > c ;
p3=2 0 0 t rffiffiffiffiffiffiffi  pc0 1 Bs21 ¼ dðtÞ þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Eðt=ðt þ c0 ÞÞ 4 pðt þ c0 Þ  c0 þ ðEðt=ðt þ c0 ÞÞ
Kðt=ðt þ c0 ÞÞÞ ; 2t cc0 s B22 ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Eðt=ðt þ c0 ÞÞ; pðt þ c0 Þ rffiffiffiffi p 1 Bs
1 ¼ dðtÞ
1=2 pffiffiffiffiffiffiffiffiffiffiffi ðKðt=ðt þ c0 ÞÞ
Eðt=ðt þ c0 ÞÞÞ; 11 c0 p t t þ c0  8 qffiffiffiffiffi c0 1 1 > pffiffiffi dðtÞ
ðEðt=c0 Þ > 3=2 pc0 c0 ðt
c0 Þ t p > > <  Bs
1
Kðt=c0 ÞÞ þ Kðt=c0 Þ ; t < c0 ; 12 ¼ > > > > : p
1 ffi Eðc0 =tÞ; t > c0 . p3=2 t ðt
c Þ 0

The remaining two elements, Bs21
1 and Bs22
1 of the inverse transformation have not been explicitly written, but due to the re-normalization, they are essentially obtained by differentiating Bs21 and Bs11 with respect to t, respectively. Appendix B. Inverse Laplace transform by use of asymptotics The time domain behavior of the wave splitting transformation kernels, discussed at the end of Section 4, can be obtained approximately for small values of the time variable by studying the asymptotics for large values of the Laplace variable. The starting point is the uniform asymptotic expressions for large arguments for the modified Bessel functions Im [25]

526

M. Johansson et al. / Computers and Structures 84 (2006) 514–527 rð~ kk
s=ck Þ

e I m ð~ kk rÞe
sr=ck pffiffiffiffiffiffiffiffiffiffiffiffi 2p~ kk r

j 1 X ð
1Þ Aj;m ~j r j k j¼0


rð~ kk þs=ck Þ

e
ie
mpi pffiffiffiffiffiffiffiffiffiffiffiffi 2p~ kk r

k

1 X j¼0

Aj;m ; j ~ kk r j

ð71Þ

for
p 6 arg ~ kk 6 0 , k = 1, 2, 3. For 0 6 arg ~ kk 6 p, one kk rÞ. uses the connection formula I m ð~ kk rÞ ¼ empi I m ðe
pi ~ The ~ kk are replaced by their asymptotic expansions for large s. Hence, (71) may then be written as an asymptotic series in inverse powers of s. Each term is transformed back to the time domain analytically, giving a series representation in the time domain for small values of the time variable. To know how to close the contours, the analyticity properties of the ~ kk must be investigated. The analyticity of the ~ kl , l = 1, 2 for the Timoshenko beam, which are essentially the same as for the Mindlin plate, has been studied before by several authors, among them [26]. They have however presented sets of branch cuts that are not suitable for our purposes since we may expect to integrate also in the left half plane. In the Laplace domain the ~ kk are given by ðqs2 ð
1Þlþ1 sðs2
s
2 Þ1=2 Þ1=2 ~ kl ¼ ; 1=2 c2 ðq
1Þ ðs2 þ c21 =r20 Þ1=2 ~ k3 ¼ ; c3

ð72Þ

where q ¼ ðc22 þ c21 Þ=ðc22
c21 Þ > 1. Consider first ~ k3 . It has the branch points s = ± ic1/r0 and is analytic for jsj > c1/ r0. The branch cut may be chosen in such a way that c3 ~k3 =s ! 1 for s ! 1. Similarly, the inner square root of ~ kl is analytic for jsj > s
1 and the branch cut may be chosen so that (s2
s
2)1/2/s ! for sffi ! 1. Then, ~ kl are analytic p1ffiffiffiffiffiffiffiffiffiffiffiffi
1 2 for jsj > s maxð1; 1= q
1Þ and the cuts may be chosen 1=2 kl =s ! 1, when s ! 1. For ~ kk , the principal so that cl ~ branch is used. Noting that the leading terms of ~ kl are s/cl, the leading terms of the exponentials of the exponential functions in (71) for large s are, for the first factor, at most a constant, and for the second factor, at most a constant after extraction of a factor e
2sr=ck . These factors may then too be expanded in negative powers of s. Combining the formulas (71) for the two half planes and using the asymptotic expansions of the ~ kk , p affiffiffiffiffiffiffiffiffiffiffiffi power ffi series, valid for s 2 U ¼ fs : jsj > s
1 maxð1; 1= q2
1; sc1 =rÞg, is obtained as 1 X kk rÞe
sr=ck

h1kjm s
ðjþ1=2Þ þð
iÞe
mpi e
2sr=ck signðImðsÞÞ I m ð~ j¼0



1 X

h2kjm s
ðjþ1=2Þ ;

j¼0

where the h1kjm and the h1kjm are coefficients, not specified here. The exponential factor e
2sr=ck will just give a time shift by 2r/ck in the time domain so the first sum will give an expansion for values of t near 0 while the second sum

gives an expansion for T values of t near 2r/ck. The first sum is analytic in s 2 U fs : Res > 0g. The terms in the first sum may then formally be integrated along a path to the right of the imaginary axis. Similarly, the second sum T is analytic in s 2 U fs : Res < 0g. The powers of the second sum may then formally be integrated along a path to the left of the imaginary axis. By a change in variable according to s = xe
ip/2 the corresponding Fourier transforms are obtained and are identified as the distributions (x ± i0)
(j + 1/2), which inverse Fourier transforms are given by [27]. So, by means of the formulas j
1=2

L
1 ½s
ðjþ1=2Þ  ¼ tþ

=Cðj þ 1=2Þ; j

j
1=2 =Cðj þ 1=2Þ; L
1 ½signðImðsÞÞs
ðjþ1=2Þ  ¼ ið
1Þ t


ð73Þ

expansions for t near 0 and 2r/ck respectively, are found as 1 X j¼0

h1kjm j
1=2 t ; Cðj þ 1=2Þ þ

e
mpi

j 1 X ð
1Þ h2kjm j
1=2 ðt
2r=ck Þ
. Cðj þ 1=2Þ j¼0

It turns out, though we do not prove it, that h2kjm = (
1)jh1kjm. Then, L
1 ½I m ð~kk rÞe
sr=ck  is an even function about t = r/ck for even m and an odd function for odd m, which simplifies the computation of the functions in the time domain somewhat. The time domain representation of K m ð~kk rÞesr=ck for small t may be computed in an analogous way. The asymptotic expression for large arguments is 1 e
zk þsr=ck X Aj;m . K m ð~kk rÞesr=ck pffiffiffiffiffiffiffiffiffiffiffiffi j 2p~kk r j¼0 ~kk rj From the T analysis for the ~kk it is seen this series is analytic for s 2 U fs : Res > 0g. Expansion into powers of s and use of (73) gives 1 X h3ijm ðj
1=2Þ tþ ; Cðj þ 1=2Þ j¼0 where h3ijm are coefficients similar to h1ijm and h2ijm. Appendix C. Discretization of singularities in the numerical computation The convolution of power singularities with smooth functions are interpreted as CPV integrals. In (62), some kernels contain power singularities, which propagate along known characteristics. When the computational region is discretized, the characteristic lines will pass through discrete points, as well as between them. We thus want to be able to treat both cases within the frame of the trapezoidal rule. The numerical implementation of these integrals has been done by modifying the discretization of the singular function in the neighborhood of the singular point. Consider the singular function g(t) = (t
t 0 )
1. At a certain discrete point in space, t 0 = kh ± d, where h is the time step, k is an integer and 0 6 d 6 h/2. Then, g(t) has been approximated by

M. Johansson et al. / Computers and Structures 84 (2006) 514–527
1

gapp ðtÞ ¼ ðt
kh  dÞ

w1 ðt
khÞ
1 þ w2 ðt
kh  h=2Þ
1 ;

ð74Þ

where w1 = (h
2d)/h and w2 = 2d/h. The absolute error is then large for t kh ± d, but the integrated error is small. If gapp(t) is convolved with some smooth function f using the CPV definition, Taylor series expansion and the trapezoidal rule, it will result in I, say. The modified function gmod app is then defined as the discrete function that give the result I, when convolved with f, using only the trapezoidal rule and the same discretization. The discrete function gmod app is then given by
w1
w2 mod þ ; gapp ¼ gð0Þ; gðhÞ; . . . ; gððk
2ÞhÞ; 2h hð1  1=2Þ

2w2 w1 w2  ; þ ; gððk þ 2ÞhÞ; . . . . ð75Þ h 2h hð1  1=2Þ Logarithmic singularities have been represented in terms of the integral of a discretized and modified power singularity. The reason for this is to obtain a uniform representation of logarithmic singularities in the imbedding equation. References [1] He S, Stro¨m S, Weston VH. Time domain wave-splittings and inverse problems. Monographs in electrical and electronic engineering. Oxford: Oxford University Press; 1998. [2] Billger DVJ, Folkow PD. The imbedding equations for the Timoshenko beam. J Sound Vib 1998;209(4):609–34. [3] Billger DVJ, Folkow PD. Wave propagators for the Timoshenko beam. Wave Motion 2003;37(4):313–32. [4] Billger DVJ, Wall DJN. A time domain algorithm for the reflection of waves on a viscoelastically restrained Timoshenko beam. Q J Mech Appl Math 1999;52(2):211–36. [5] Folkow PD, Kreider K. Direct and inverse problems on nonlinear rods. Math Comput Simul 1999;50(5–6):577–95. [6] Romeo M. Wave splitting in linear viscoelasticity. Eur J Mech A— Solids 1999;18:539–55. [7] Folkow PD. Time domain inversion of a viscoelastically restrained Timoshenko beam. Inverse Problems 1999;15(2):551–62. [8] Weston VH. Invariant imbedding for the wave equation in three dimensions and the application to the direct and inverse problems. Inverse Problems 1990;6(6):1075–105.

527

[9] Weston VH. Invariant imbedding and wave splitting in R3 : II. The Green function approach to inverse scattering. Inverse Problems 1992;8(6):919–47. [10] Karlsson A, Kristensson G. Wave splitting in the time domain for a radially symmetric geometry. Wave Motion 1990;12:197–211. [11] Kreider KL. Time dependent direct and inverse electromagnetic scattering for the dispersive cylinder. Wave Motion 1989;11:427– 40. [12] Kreider KL. A wave splitting approach to time dependent inverse scattering for the stratified cylinder. SIAM J Appl Math 1989;49(3): 932–43. [13] Timoshenko SP. On the correction for shear of the differential equation for transverse vibrations of prismatic bars. Philos. Mag. 1921;XLI:744–6. reprinted in: The Collected Papers of Stephen P. Timoshenko, McGraw-Hill, London, 1953. [14] Olsson P, Kristensson G. Wave splitting of the Timoshenko beam equation in the time domain. Z Angew Math Phys 1994;45: 866–81. [15] Folkow PD, Olsson P, Kristensson G. Time domain Green functions for the homogeneous Timoshenko beam. Q J Mech Appl Math 1998; 51(1):125–41. [16] Krueger RJ, Ochs RL. A Green’s function approach to the determination of internal fields. Wave Motion 1989;11:525–43. [17] Mindlin RD. Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates. J Appl Mech 1951;18:31–8. [18] Vemula C, Norris AN. Flexural wave propagation and scattering on thin plates using Mindlin theory. Wave Motion 1997;26:1–12. [19] Erde´lyi A. Tables of integral transforms. The Bateman manuscript project, vol. 1. New York: McGraw-Hill; 1954. [20] Oberhettinger F, Badii L. Tables of Laplace transforms. Berlin: Springer-Verlag; 1973. [21] Folkow PD. Wave propagation in structural elements—direct and inverse problems in the time domain. PhD thesis. Chalmers University of Technology. 1998. [22] Jones DS. Infinite integrals and convolutions. Proc R Soc London A 1980;371:479–508. [23] Bostro¨m A. Time-dependent scattering by a bounded obstacle in three dimensions. J Math Phys 1982;23(8):1444–50. [24] Oldham KB, Spanier J. The fractional calculus. Mathematics in science and engineering, vol. 111. London: Academic Press; 1974. [25] Wilcox CH, Langer RE, editorsAsymptotic solutions of differential equations and their applications. New York: Publication of the Mathematics Research Center, United States Army, University of Wisconsin, Wiley; 1964. [26] Miklowitz J. The theory of elastic waves and waveguides. 2nd ed. North-Holland series in applied mathematics. Amsterdam: NorthHolland; 1978. [27] Gel’fand IM, Shilov GE. Generalized functions, vol. 1. London: Academic Press; 1964.