- Email: [email protected]
Conformal mapping methods for variable parameter elastodynamics
Wave Motion 36 (2002) 185–202
Conformal mapping methods for variable parameter elastodynamics George D. Manolis a,∗ , Tsviatko V. Rangelov b , Richard Paul Shaw c b
a Department of Civil Engineering, Aristotle University, Thessaloniki GR-54006, Greece Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, Sofia 1113, Bulgaria c Department of Civil Engineering, State University of New York, Buffalo, NY 14260, USA
Received 4 December 2001; accepted 15 January 2002
Abstract This work is a generalization of earlier efforts to develop free-space fundamental solutions (Green’s functions) for a certain class of non-homogenous materials (i.e., materials with position-dependent elastic parameters and density) under time-harmonic conditions by means of conformal mapping methods. These methods were shown to be very effective in conjunction with the scalar wave equation, but their extension to the 2D and 3D vector wave equations of elastodynamics is hampered by the following two factors: (a) the governing equations of motion are vectorial, which implies that any transformation of coordinates affects the basis vectors used in the original (i.e., Cartesian) coordinate system; and (b) an extension of conformal mapping to a 3D coordinate system is not readily obvious. In the present work, both issues raised above are addressed. Furthermore, the necessary background is established whereby the construction of Green’s functions for 2D and 3D elastodynamics can be systematically accomplished in a mapped space, followed by an inversion to the original Cartesian coordinate system. Of course, conformal mapping as used here is inherently an inverse method, in the sense that the material parameter profiles for a given problem are recovered as constraints during the solution procedure and therefore depend on the type of mapping used. © 2002 Elsevier Science B.V. All rights reserved.
1. Introduction The development of fundamental solutions for elastic wave propagation in solids that are not homogeneous and isotropic, but instead exhibit heterogeneity, anisotropy, layering, randomness, etc., are important in their own right [1]. At the same time, fundamental solutions for the free-space (known as Green’s functions) comprise an essential part of any boundary integral equation method (BIEM) formulation [2] and, of course, of its corresponding numerical solution technique, the boundary element method (BEM). The latter method has been very successful in solving boundary-value problems of engineering importance in transient elastodynamics [3,4]. Since the preferred formulation is invariably a displacement-traction approach based on Somigliana’s identity, the corresponding Green’s functions must come from a solution of the equations of elastodynamics for a point force and under radiation-type boundary conditions. This can be done by a variety of ways, e.g., use of potentials (Helmholtz potentials, Stokes potentials), displacement vector decomposition into dilatational and rotational components, use of the dynamic equivalent to Galerkin’s vector, integral transforms or other specialized approaches. ∗
Corresponding author. Tel.: +30310-995-663; fax: +30310-995-769. E-mail address: [email protected] (G.D. Manolis). 0165-2125/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 2 1 2 5 ( 0 2 ) 0 0 0 1 1 - 2
186
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
Any extension of the aforementioned methodologies to problems involving other categories of media, which go beyond the homogeneous, isotropic and linear elastic continuum, has proven difficult. This, however, would be very useful in various engineering fields, such as seismology, earthquake engineering, non-destructive testing evaluation, acoustics, electromagnetic signal transmission, ocean engineering, etc., where the materials through which various types of signals propagate cannot be satisfactorily modeled as homogeneous [5–8]. In this respect, we briefly mention some recent activity in the area of time-harmonic acoustic waves [9] and elastic waves [10] in heterogeneous media using a coordinate transformation based on conformal mapping. In the latter case of elastic waves, which are governed by vector differential equations, the class of solvable problems was limited to variable (i.e., position-dependent) density and constant elastic parameters. Even then, an additional assumption regarding decomposition of the displacement vector into pseudo-dilatational and pseudo-rotational parts was necessary in order to achieve a closed-form solution through reduction to the anti-plane strain case. While this approach at first sight appears to be limited to 1D (or, at best, to 2D) problems and is inverse in form, i.e., a given mapping will lead to a particular type of heterogeneity, it represents a step forward in the development of Green’s functions of the type required by the BIEM and BEM. In this work, we examine in detail the coordinate transformations implied through conformal mapping for the vector wave equations of elastodynamics under both 3D and plane strain conditions, whereby the latter case is proved to be a special case of the former. Furthermore, we seek to establish the conditions under which it is possible to obtain a solution in the transformed, orthogonal curvilinear coordinate system into which the original (Cartesian) coordinates are mapped. The key advantage of using what are essentially coordinate transformations is that the inverse mapping, which is required to give the final Green’s function in Cartesian coordinates, is straightforward to perform. Although the variable density case represents the simplest material heterogeneity, it is possible to look at the general case as well, where all material parameters of the medium in question vary with position.
2. The governing equations of elastodynamics in 3D Navier’s equations [11] of dynamic equilibrium for time-harmonic elastodynamics written in the reference Cartesian coordinate system {X} are L(U) + ρω2 U ≡ (λ + 2µ)∇∇ · U − µ∇ × ∇ × U + ρω2 U = F ,
(1)
where λ, µ are the Lame constants, ρ is the mass density, and ω the frequency of vibration. Furthermore, U = (U1 , U2 , U3 ) is the displacement vector, F the body force per unit mass, and the two vector operators used in the above formulation are ∇∇ · U = grad div U and ∇ × ∇ × U = curl curl U. Let I map {X} into an orthogonal curvilinear coordinate system {x}, where xj = xj (X), with index j = 1, 2, 3 as shown in Fig. 1. The Jacobian matrix of the transformation is defined as J = (∂x/∂X), provided that det J = |J| = 0. Also, g = JJT is the metric tensor in the new coordinate system {x} given as h1 0 0 g = 0 h2 0 0 0 h3 with hj > 0. If I is a conformal mapping, then h1 = h2 = h3 = h2 , where h = |J|1/3 [12]. Furthermore, we need to define as {Ei } the orthogonal basis in the reference system {X}, which maps onto the curvilinear orthogonal basis {ei } of {x} through I. Obviously, what complicates transformations between coordinate systems is that the base vectors are dependent on the spatial coordinates, which implies use of the chain rule for differentiation of vector functions [13]. The transformed Navier’s equations in the new coordinates {x} are obtained using the appropriate representations for the differential operators acting on vector functions. Specifically, let the displacement vector be expressed either
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
187
Fig. 1. Conformal mapping between two orthogonal curvilinear coordinate systems in: (a) 3D space; (b) 2D space.
in terms of the original coordinate system {X} as U = U1 E1 + U2 E2 + U3 E3 or in terms of the transformed coordinate system {x} as u = u1 e1 + u2 e2 + u3 e3 . If a solution can be obtained in the transformed coordinates, then the solution of Eq. (1) is simply recovered by applying the inverse transformation I−1 so that U(X) = u(x(X)). There is an inherent duality in these transformations, because although the {X} space is originally viewed as Cartesian and the {x} as orthogonal curvilinear, any fundamental solution obtained in the latter space can be treated as the Cartesian coordinate system solution. Then, the inverse mapping can be used to return this solution back to the {X} coordinates, which now comprise an orthogonal curvilinear system as seen from {x}. In what follows, we list the transformation of all key operators, noting that the summation convention is implied for repeated indices. Case 1: The gradient operator for a scalar displacement function is ∇U = Ej ·
∂U ∂Xj
in {X},
1 ∂u ej · h ∂xj
in {x}.
(2)
1 ∂(uj h2 ) h3 ∂xj
in {x}.
(3)
∇u =
Case 2: The divergence operator is ∇ ·U =
∂Uj ∂Xj
in {X},
∇ ·u=
Case 3: The curl operator is E1 E2 E3 ∂ ∂ ∂ ∇ × U = ∂X ∂X ∂X 1 2 3 U U U 1 2 3
in {X},
1 ∇ × u = 3 h
he1 ∂ ∂x1 hu1
he2 ∂ ∂x2 hu2
he3 ∂ ∂x3 hu3
in {x}.
(4)
Case 4: The grad div operator is computed by applying the second equation (2) to the scalar function defined by the second equation (3). Specifically, we first recover the intermediate result ∇ · u = (1/ h3 )∂k (uk h2 ) = (1/ h)(∂k uk + 2∂k H · uk ), where the notation ∂/∂xj = ∂j and ln|h| = H has been introduced. Finally, 1 1 ∇∇ · u = ej ∂j (∂k uk + 2∂k H · uk ) h h 1 = 2 ej [∂jk2 uk − ∂j H · ∂k uk + 2∂k H · ∂j uk + (−2∂j H · ∂k H + 2∂jk2 H )uk ]. (5) h
188
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
Case 5: The curl curl operator is computed by applying the second equation (4) to the vector function curl u. Since the cross product yields ∇ × u = (1/ h2 )ej wj , where the three components are w1 = (1/ h2 )[∂2 (hu3 ) − ∂3 (hu2 )], w2 = (1/ h2 )[∂3 (hu1 ) − ∂1 (hu3 )] and w3 = (1/ h2 )[∂1 (hu2 ) − ∂2 (hu1 )], then he1 he2 he3 1 ∂ ∂ ∂ ∇ × ∇ × u = 3 h ∂x1 ∂x2 ∂x3 hw hw hw 1 2 3 1 = 2 {e1 [∂2 (hw3 ) − ∂3 (hw2 )] + e2 [∂3 (hw1 ) − ∂1 (hw3 )] + e3 [∂1 (hw2 ) − ∂2 (hw1 )]}. (6) h The coefficient multiplying vector component e1 in the above equation is given by the following expression: 1 1 1 1 1 [∂2 (hw3 ) − ∂3 (hw2 )] = 2 ∂2 h · w3 + ∂2 w3 − 2 ∂3 h · w2 − ∂3 w2 = W11 + W12 + W13 + W14 . (7) h h h2 h h Carrying through with the differentiations yields 1 1 1 1 ∂2 h · 2 [∂1 (hu2 ) − ∂2 (hu1 )] 2 ∂2 H · [∂1 h · u2 + h∂1 u2 − ∂2 h · u1 − h∂2 u1 ] h h2 h h 1 = 2 [−∂2 H · ∂2 u1 + ∂2 H · ∂1 u2 − (∂2 H )2 u1 + ∂1 H · ∂2 H · u2 ], h
1 1 W12 = 2 h∂2 2 [∂1 (hu2 ) − ∂2 (hu1 )] h h
1 2 1 = − 3 ∂2 h[∂1 (hu2 ) − ∂2 (hu1 )] + 2 ∂2 [∂1 h · u2 + h∂1 u2 − ∂2 h · u1 − h∂2 u1 ] h h h 1 = 2 [−2∂2 H (∂1 H · u2 + ∂1 u2 − ∂2 H · u1 − ∂2 u1 )] h
2 h 2 h ∂22 1 ∂21 2 2 + 2 · u2 + ∂1 H ∂2 u2 + ∂2 H ∂1 u2 + ∂12 u2 − · u1 − ∂2 H ∂2 u1 − ∂2 H ∂2 u1 − ∂22 u1 h h h W11 =
=
1 2 2 2 (−∂22 u1 + ∂12 u2 − ∂2 H · ∂1 u2 + ∂1 H · ∂2 u2 − [∂22 H − (∂2 H )2 ] · u1 h2 2 +[∂21 H − (∂1 H.∂2 H )] · u2 ),
(8)
where the equality ∂ij2 h/ h = ∂ij2 H + ∂i H · ∂j H has been used. The remaining two parts are obtained in a similar fashion, i.e., 1 1 1 ∂3 h · 2 [∂3 (hu1 ) − ∂1 (hu3 )] = 2 [−∂3 H · ∂3 u1 + ∂3 H · ∂1 u3 − (∂3 H )2 u1 + ∂1 H · ∂3 H · u3 ], 2 h
h h 1 1 W14 = − 2 h∂3 2 [∂3 (hu1 ) − ∂1 (hu3 )] h h 1 2 2 2 u1 + ∂13 u3 + ∂1 H · ∂3 u3 − ∂3 H · ∂1 u3 − [∂33 H − (∂3 H )2 ] · u1 = 2 (−∂33 h 2 +[∂13 H − (∂1 H ) · (∂3 H )] · u3 ). (9) W13 = −
Since we have evaluated the first component of the curl curl operator, namely (∇×∇×u)1 = W11 +W12 +W13 +W14 , the remaining two components are derived by cyclic interchange of the indices as (∇ × ∇ × u)j =
1 2 {−∇ 2 uj + ∂kj2 uk − ∂k H · ∂k uj + ∂j H · ∂k uk + ∂jk2 H · uk − ∂kk H · uj }. h2
(10)
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
189
Case 6: The Laplacian ∇ 2 in the transformed coordinates {x} can be computed from the grad div and the curl curl operators as follows: ∇ 2 u = ∇∇ · u − ∇ × ∇ × u.
(11)
Thus, the jth component (j = 1, 2, 3) is equal to (∇ 2 u)j =
1 2 {∇ 2 uj + ∂k H · ∂k uj + 2∂k H · ∂j uk − 2∂j H · ∂k uk + ∂kk H · uj + [∂jk2 H − 2(∂k H )(∂j H )]uk }. h2 (12)
2.1. The transformed Navier’s equations By combining all the previous results, the components of the differential operator L that appeared in Eq. (1) are now (Lu)j =
1 {µ∇ 2 uj + (λ + µ)∂jk2 uk + µ∂k H · ∂k uj + 2(λ + 2µ)∂k H · ∂j uk − (λ + 3µ)∂j H · ∂k uk h2 2 +µ∂kk Huj + [(2λ + 3µ)∂jk2 H − 2(λ + 2µ)∂j H · ∂k H ]uk }.
(13)
The first two terms in the above equation comprise Lame’s operator, but we see the emergence of additional terms that can be viewed as the influence of inhomogeneities in the original material. The key objective is to try and match these extra terms to a particular density and/or elastic parameter profile. For instance, in previous work [14], it was possible to recover the standard form of Navier’s equation in the {x} coordinates by prescribing a density profile equal to ρ(x, y) = ρ0 |J(x, y)|. Navier’s system of equations in the transformed coordinate system {x} is now (Lu)j + ρω2 uj = fj ,
(14)
where (Lu)1 ≡
1 2 2 2 2 2 {(λ + 2µ)∂11 u1 + µ(∂22 u1 + ∂33 u1 ) + (λ + µ)(∂12 u2 + ∂13 u3 ) + (λ + 2µ)∂1 H · ∂1 u1 h2 +µ(∂2 H · ∂2 u1 + ∂3 H · ∂3 u1 ) + 2(λ + 2µ)(∂2 H · ∂1 u2 + ∂3 H · ∂1 u3 )
−(λ + 3µ)(∂2 H · ∂2 u2 + ∂3 H · ∂3 u3 ) + (K10 + K11 )u1 + K12 u2 + K13 u3 }, 1 2 2 2 2 2 (Lu)2 ≡ 2 {(λ + 2µ)∂22 u2 + µ(∂11 u2 + ∂33 u2 ) + (λ + µ)(∂21 u1 + ∂23 u3 ) + (λ + 2µ)∂2 H · ∂2 u2 h +µ(∂1 H · ∂1 u2 + ∂3 H · ∂3 u2 ) + 2(λ + 2µ)(∂1 H · ∂2 u1 + ∂3 H · ∂2 u3 ) −(λ + 3µ)(∂1 H · ∂1 u1 + ∂3 H · ∂3 u3 ) + K21 u1 + (K20 + K22 )u2 + K23 u3 }, 1 2 2 2 2 2 (Lu)3 ≡ 2 {(λ + 2µ)∂33 u3 + µ(∂11 u3 + ∂22 u3 ) + (λ + µ)(∂31 u1 + ∂32 u2 ) + (λ + 2µ)∂3 H · ∂3 u3 h +µ(∂1 H · ∂1 u3 + ∂2 H · ∂2 u3 ) + 2(λ + 2µ)(∂1 H · ∂3 u1 + ∂2 H · ∂3 u2 ) −(λ + 3µ)(∂1 H · ∂1 u1 + ∂2 H · ∂2 u2 ) + K31 u1 + K32 u2 + (K30 + K33 )u3 }.
(15)
2 H + ∂ 2 H + ∂ 2 H ] and K = (2λ + 3µ)∂ 2 H − 2(λ + In the above equations, we have that Kj 0 = µ[∂11 jk 22 33 jk 2µ)(∂j H )(∂k H ), where j, k = 1, 2, 3. Note at this point that index j in the first expression is superfluous, since K10 = K20 = K30 , but is retained for notational convenience. From a mechanical point of view, the Kj k terms in the transformed equations correspond to spatial curvature, while the first order derivatives of the dependent variable u indicate dispersion phenomena. The effect of the curvature terms is localized and, from a practical viewpoint, it is possible to approximate their influence by converting the original wave numbers into complex quantities. The dispersive terms, on the other hand, produce damping that is responsible for signal de-amplification and dispersion.
190
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
This effect can noticeably alter the propagating signal, especially as the distance from source increases and at high frequencies of vibration. 2.2. Particular cases In this section, we will examine some special cases concerning the Navier’s equations in the transformed coordinate system {x}. Case 1: We assume here that each displacement component depends on one and only one spatial variable. Ordinarily, this would lead to three separate 1D equations of motion, but as will be shown here, there is coupling due to the spatial curvature and the first order terms mentioned above. Specifically, for u = (u1 , u2 , u3 ), we have that u1 = u1 (x1 ), u2 = u2 (x2 ), u3 = u3 (x3 ) and any solutions in this form must satisfy the following system: 1 2 {(λ + 2µ)∂11 u1 + (λ + 2µ)∂1 H · ∂1 u1 − (λ + 3µ)(∂2 H · ∂2 u2 + ∂3 H · ∂3 u3 ) h2 +(K10 + K11 )u1 + K12 u2 + K13 u3 } + ρω2 u1 = f1 , 1 2 (Lu)2 + ρω2 u2 ≡ 2 {(λ + 2µ)∂22 u2 + (λ + 2µ)∂2 H · ∂2 u2 − (λ + 3µ)(∂1 H · ∂1 u1 + ∂3 H · ∂3 u3 ) h +K21 u1 + (K20 + K22 )u2 + K23 u3 } + ρω2 u2 = f2 , 1 2 (Lu)3 + ρω2 u3 ≡ 2 {(λ + 2µ)∂33 u3 + (λ + 2µ)∂3 H · ∂3 u3 − (λ + 3µ)(∂1 H · ∂1 u1 + ∂2 H · ∂2 u2 ) h +K31 u1 + K32 u2 + (K30 + K33 )u3 } + ρω2 u3 = f3 . (Lu)1 + ρω2 u1 ≡
(16)
Case 2: This parallels the standard plane strain conditions of classical elastodynamics. The relevant assumptions are well known, i.e., we have that u1 = u1 (x1 , x2 ), u2 = u2 (x1 , x2 ), u3 = u3 (x3 ). Solutions in this form must satisfy the following system of coupled equations: (Lu)1 + ρω2 u1 ≡
1 2 2 2 {(λ + 2µ)∂11 u1 + µ∂22 u1 + (λ + µ)∂12 u2 + (λ + 2µ)∂1 H · ∂1 u1 + µ∂2 H · ∂2 u1 h2 +2(λ + 2µ)∂2 H · ∂1 u2 − (λ + 3µ)(∂2 H · ∂2 u2 + ∂3 H · ∂3 u3 ) + (K10 + K11 )u1
+K12 u2 + K13 u3 } + ρω2 u1 = f1 , 1 2 2 2 (Lu)2 + ρω2 u2 ≡ 2 {(λ + 2µ)∂22 u2 + µ∂11 u2 + (λ + µ)∂21 u1 + (λ + 2µ)∂2 H · ∂2 u + µ∂1 H · ∂1 u2 h +2(λ + 2µ)∂1 H · ∂2 u1 − (λ + 3µ)(∂1 H · ∂1 u1 + ∂3 H · ∂3 u3 ) + K21 u1 +(K20 + K22 )u2 + K23 u3 } + ρω2 u2 = f2 , 1 2 u3 + (λ + 2µ)∂3 H · ∂3 u3 − (λ + 3µ)(∂1 H · ∂1 u1 + ∂2 H · ∂2 u2 ) (Lu)3 + ρω2 u3 ≡ 2 {(λ + 2µ)∂33 h +K31 u1 + K32 u2 + (K30 + K33 )u3 } + ρω2 u3 = f3 .
(17)
Case 3: This last case parallels the anti-plane strain condition of elastodynamics, whereby the displacement vector components are u1 = 0, u2 = 0, u3 = u3 (x1 , x2 ). Then, any solution must satisfy the system 1 {(λ + µ)∂3 H · ∂1 u3 + K13 u3 } = f1 , h2 1 (Lu)2 + ρω2 u2 ≡ 2 {(λ + µ)∂3 H · ∂2 u3 + K23 u3 } = f2 , h 1 2 2 (Lu)3 + ρω2 u3 ≡ 2 {µ(∂11 u3 + ∂22 u3 ) + µ(∂1 H · ∂1 u3 + ∂2 H · ∂2 u3 ) + (K30 + K33 )u3 } + ρω2 u3 = f3 . h (18) (Lu)1 + ρω2 u1 ≡
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
191
3. Conformal mappings in 3D In the previous section, conformal mapping was used to transform Navier’s equations from orthogonal to curvilinear coordinates. In order to quantify the terms h, ∂j H, ∂jk2 H, Kjk that appear in Eq. (15), it is necessary to classify conformal mappings that are possible in 3D [15]. At first, the mapping from {X} to {x} is conformal if the metric G in {X} is proportional to the metric g in {x}, i.e., there exists a function f (x) = 0 such that gαβ = f (x)Gαβ . In our case, 1, α = β, Gαβ = 0, α = β, since it is a Cartesian metric. Thus, f (x), α = β, gαβ = 0, α = β, and f (x) = h2 (x) = |J|1/3 , where J is the Jacobian matrix of the transformation I at x. Next, we list the following three basic cases of conformal mappings in 3D: (a) Rotation: Let x = AX, where A is an orthogonal matrix, i.e., AT = A−1 . Then J = A, g = JJT = E, where E is the identity matrix, h(x) ≡ 1, and finally gαβ = Gαβ . (b) Dilatation: Let x = h0 X, where h0 is a constant; then J = h0 E, h(x) ≡ h0 , and gαβ = h20 Gαβ . Note that both rotation and dilatation are linear transformations. (c) Inversion: Let x = (X − X0 )/|X − X0 |2 , where X0 is a fixed point. Such a transformation is defined in R3 \{X 0 } and h(X) = −|X − X0 |−2 , |J| = −|X − X0 |−6 , and gαβ = |X − X0 |−4 Gαβ . A theorem by Liouville [15] classifies all smooth conformal mappings in 3D as a superposition of rotation, dilatation and inversion. It therefore follows that there can be no conformal transformations in 3D other than those composed from cases (a)–(c), plus rigid body motion. This is in sharp contrast to 2D, where every analytic function in R2 defines a conformal mapping. Next, we present some details regarding the particular form assumed by operator L in Eq. (14) for the aforementioned three basic cases. Specifically, for case (a), transformation I is defined as x = AX, the inverse transformation I−1 is X = A−1 x, and we have that h = 1, H = ln h = 0, ∂j H = 0, Kjk = 0 and L is the same as the operator that appears in Eq. (1). Next, for case (b), transformation I is defined as x = h0 X, the inverse transformation −4 2 I−1 is X = h−1 0 x, h = h0 , H = ln h = 0, ∂j H = 0, Kjk = 0 and L is that of Eq. (1) scaled by h0 . Finally, transformation I for case (c) is x = X/|X|2 , where we let X0 = 0 for simplicity. The inverse transformation I−1 is X = x/|x|2 and h(x) = −1/|x|2 , H = −2 ln |x|. There is no simple relation between the transformed operator L in Eq. (13) and the original one; details for this case will be given in the next section.
4. Solutions for the reduced anti-plane strain case We will pursue some solutions of the system of equations that pertain to the anti-plane strain case as reduced from 3D. Specifically, if the force components in the transformed coordinates f1 , f2 , f3 are continuous, Eq. (18) have a solution u3 locally for the following forms of function H (=ln|h|): (a) H is constant in a neighborhood Ω of field point x0 . (b) The derivative of H with respect to the third spatial coordinate ∂3 H = 0 at x0 and consequently in the neighborhood of that point.
192
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
In the first case, ∂j H = 0 for j = 1, 2, 3 in Ω and the left-hand sides of the first two equations (18) are identically zero, which implies that components f1 , f2 must also be zero in order to have a solution. Then, Eq. (18) are reduced to just the last equation 1 2 2 µ(∂11 u3 + ∂22 u3 ) + ρω2 u3 = f3 . h2
(19)
Let u˜ 3 be a fundamental solution of Eq. (19) for a unit point force (Dirac’s delta function), i.e., 2 2 ∂11 u˜ 3 + ∂22 u˜ 3 + k02 u˜ 3 = δ,
(20) (1)
where k02 = e2H ρω2 /µ is a wave number and e2H = h2 . The fundamental solution is u˜ 3 = −(i/4)H0 (k0 |x|), while that of Eq. (19) is u3 = u˜ 3 ∗ {f3 e2H /µ}, where ∗ is the convolution operator [16]. The actual density profile for which this solution is valid is non-constant and equal to ρ(x) = µ(x){k02 /ω2 e2H }. More specifically, we recover that ρ(x)/µ(x) = ρ0 /µ0 , where ρ 0 , µ0 are reference (or background) values identified with a key point, such as the origin of the coordinate system and that the true wave number is k = k0 / h. A similar profile, whereby the density is proportional to the shear modulus, was derived by Manolis and Shaw [17], who worked directly with the 2D Helmholtz equation. In the second case, the first and second equations (18) yield ∂1 u3 = fˆ1 − Kˆ 13 u3 , ∂2 u3 = fˆ2 − Kˆ 23 u3 ,
h 2 f1 K13 , Kˆ 13 = , (λ + µ)∂3 H (λ + µ)∂3 H h 2 f2 K23 fˆ2 = , Kˆ 23 = , (λ + µ)∂3 H (λ + µ)∂3 H
fˆ1 =
(21)
while the third one is 2 2 ∂11 u3 + ∂22 u3 + kˆ02 u3 = fˆ3 ,
(22)
where the modified forcing function and the modified wave number are equal to h fˆ3 = f3 − ∂1 H · fˆ1 − ∂2 H.fˆ2 , µ 2
ρω2 h2 1 kˆ02 = − ∂1 H · Kˆ 13 − ∂2 H · Kˆ 23 + (K30 + K33 ). µ µ
(23)
Eq. (22) is of the Helmholtz type, and given that u˜ 3 is a fundamental solution for a unit point force, then the complete solution is u3 = u˜ 3 ∗ fˆ3 . We note here that since h(x) is known, density ρ(x) and shear modulus µ(x) in Eq. (23) must combine in such a way so as to give a constant value for kˆ02 . Remark 1. Following the smooth conformal mapping I from {X} to {x} in the anti-plane strain case, Liouville’s theorem implies that at an arbitrary field point x0 , one of the following possibilities must hold true: (i) ∂1 H = ∂2 H = ∂3 H = 0, in which case Eq. (18) is solvable as in case (a). Essentially, the problem uncouples if the forcing function components corresponding to zero displacements are also zero. (ii) ∂j H = 0 for j = 1, 2, 3. Then, there exists a neighborhood Ω of x0 such that ∂j H = 0 and the system of equations (18) are solvable for uj as in case (b). Specifically, if ∂1 H = 0, we have a solution for u1 (x2 , x3 ), u2 = u3 = 0, and if ∂2 H = 0, the solution is for u2 (x1 , x3 ), u1 = u3 = 0. The last case, where ∂3 H = 0, is the same as case (b). In what follows, we will further examine the anti-plane strain case in conjunction with two specific types of conformal mappings, namely the linear (dilatation plus rotation) and the inverse that were presented in Section 3.
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
193
4.1. The anti-plane strain case for a conformal mapping of the linear type Let I be a linear conformal mapping, defined by a constant 3 × 3 Jacobian matrix J, det J = |J| = 0, such that (1/|J|)J is an orthogonal matrix. Therefore, x = JX, |J| = h30 and the metric tensor is h20 0 0 2 g = JJT = 0 h0 0 , 0 0 h20 which implies that h = h0 and H = ln|h| = 0 (h0 is a dimensionless factor). More precisely, the Jacobian matrix is J = [ajk ] and the orthogonality condition can be summarized as 3 h20 , i = j, aik ajk = 0, i = j. k=1 For this type of mapping, Eq. (18) read as (Lu)1 + ρω2 u1 ≡ 0 = f1 , (Lu)2 + ρω2 u2 ≡ 0 = f2 , µ 2 2 (∂11 u3 + ∂22 u3 ) + ρω2 u3 = f3 . (Lu)3 + ρω2 u3 ≡ h20
(24)
We distinguish the following two possibilities: (i) If (f1 )2 + (f2 )2 = 0, then Eq. (24) has no solution. (1) (ii) If√f1 = f2 = 0, the Helmholtz solution previously given as u3 = u˜ 3 ∗ {f3 h20 /µ}, where u˜ 3 = −(i/4)H0 { ρ0 /µ0 h0 ω|x|}, is valid. Therefore, the corresponding solution to Navier’s equations (1) in the original Cartesian coordinate system {X} is U = u(x(X)) with U1 = 0, U2 = 0 and 3 3 U3 = u3 (x1 (X), x2 (X)) = u3 ((JX)1 , (JX)2 ) = u3 a1j Xj , (25) a2j Xj . j =1
j =1
4.2. The anti-plane strain case for a conformal mapping of the inverse type Let conformal mapping I be an inversion defined from R3 \{0} to R3 \{0} as x = X/|X|2 , while I−1 is defined as X = x/|x|2 . If the radial distance in {X} is denoted by R = |X|, then R,j = X,j /R and the elements of the Jacobian matrix J are Jjk = (1/R 2 )(δjk − 2R,j R,k ). Furthermore, det J = |J| = −1/R 6 , h(x) = −r −2 , H = −2 ln|x| and the metric tensor is h2 (x), α = β, gαβ = 0, α = β. For this mapping, the following terms appearing in the transformed Navier’s equations (18) are computed: 2 2 4 2µ ∂j H = − r,j , ∂jk2 H = − 2 (δjk − 2r,j r,k ), (∂j H ) · (∂k H ) = 2 r,j r,k , Kj 0 = − 2 , r r r r 2 (26) Kjk = − 2 [(2λ + 3µ)δjk + 2µr,j r,k ]. r In the above, the radial distance in {x} is defined as r = |x| with r,j = xj /r. As defined, the {x} coordinate system has dimensions of inverse length. We conclude that Eq. (18) are solvable in every domain B ⊂ R3 \{x3 = 0}, since
194
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
∂3 H = −(2/r)r,3 = 0. The solution is u3 = u˜ 3 ∗ fˆ3 , where u˜ 3 is a fundamental solution of Eq. (22) and fˆ3 , kˆ02 are given in Eq. (23). The corresponding solution to Navier’s equations (1) in Cartesian coordinates {X} is X1 X2 U3 = u3 U1 = U2 = 0, . (27) , |X|2 |X|2 5. The governing equations of elastodynamics in 2D In 2D elastodynamics, the key problem that arises as one seeks to define a suitable transformation between the reference Cartesian coordinates {X} and new curvilinear coordinates {x} is that the standard definition for the cross product introduces a third dimension (see Fig. 1). To circumvent this problem, Navier’s equations for the plane strain, time-harmonic case are presented in a form that no longer involves the curl curl operator. Specifically, for a simply connected domain B ⊂ R2 , we have that LU + ρω2 U ≡ µ∇ 2 U + (λ + µ)∇∇ · U + ρω2 U = F ,
(28)
where λ, µ, ρ, ω have all been previously defined, while U = (U1 , U2 ) and F = (F1 , F2 ), respectively, are the displacement and body force per unit mass vectors. The standard definitions for the vector operators are used, i.e., ∂U1 ∂U2 ∂U ∂U , div U = ∇ · U = , + , grad U = ∇U = ∂X1 ∂X2 ∂X1 ∂X2 ∂ 2 U1 ∂ 2 U2 ∂ 2 U1 ∂ 2 U2 grad div U = ∇∇ · U = + , , + ∂X1 ∂X2 ∂X1 ∂X2 ∂X12 ∂X22 ∂ 2 U1 ∂ 2 U 1 ∂ 2 U2 ∂ 2 U2 2 . (29) + , + div grad U = ∇ U = ∂X12 ∂X22 ∂X12 ∂X22 Let I map domain B into a 2D manifold M embedded in R3 . The coordinate system used in B is {X}, while that used in M is {x}, where xj = xj (X). The Jacobian matrix of the transformation is J = (∂x/∂X), where det J = 0. If we define gαβ = JJT as the contravariant metric tensor, then gαβ is the corresponding covariant tensor β and g αβ · gαβ = δα (i.e., the Kronecker delta). If I is an orthogonal mapping, then at every point in M the following conditions are fulfilled: hα , α = β, (hα )−1 , α = β, gαβ = and g αβ = 0, α = β, 0, α = β. Furthermore, if I is conformal, then we have the additional condition h1 = h2 = h2 = |J| at every point in M. As before, we denote by {Ei } the orthogonal basis in {X}, which is mapped into the curvilinear orthogonal basis {ei } in {x}. The objective is to solve Navier’s system of equations in the transformed coordinates and, by applying the inverse mapping I−1 , obtain the solution of Eq. (28) as U(X) = u(x(X)). It is now necessary to re-define the vector operations in Eq. (29) for the transformed coordinate system. In this respect, we follow Ref. [11] and also use tensor notation [12], where the summation rules for rising and lowering indices is implied along with summation for repeated symbols. The relevant expressions are given below for a scalar function ϕ and a vector function w on M (note that the formulas employing tensor notation use Greek symbols, while standard indices are used for vector notation; also, ∂s = ∂/∂xs ). Case 1: Gradient operator for scalar functions: (grad ϕ)α = (∇ϕ)α = g αβ ∂β ϕ,
(grad ϕ)j = (∇ϕ)j =
1 ej ∂j ϕ. h
(30)
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
195
Case 2: Gradient operator for vector functions: (grad w)αβ = (∇w)αβ = g βγ ∂γ wα ,
grad w = ∇w =
1 ej ∂j w. h
(31)
Case 3: Divergence operator: 1 β div w = ∇ · w = √ (∂α w γ + Γνβ w ν ), g
div w = ∇ · w =
1 ∂j (wj h). h2
(32)
k is In the above, wα is a covariant vector, wα is a contravariant vector and w α = g αβ wβ , wδ = gδγ w γ . Also, Γαβ 1 k ks the Christoffel tensor and Γαβ = Γsαβ g , where Γsαβ = 2 (∂β gsα + ∂α gβs − ∂s gβα ). The next step is to recover the grad div and Laplace operators. First, define ln|h| = H and note that ∂ij2 h/ h = ∂ij2 H + ∂i H · ∂j H . Then, by applying Eq. (30) to the vector function defined by Eq. (32), the grad div operator is
(∇∇ · w)j =
1 2 [∂ wk + ∂k H · ∂j wk − ∂j H · ∂k wk + (∂jk2 H − ∂j H · ∂k H )wk ]. h2 jk
(33)
For the Laplacian, we use Eqs. (31) and (32) and get (∇ 2 w)j =
1 2 [∇wj − 2∂j H · ∂k wk + ∂k H · ∂k wj + ∂k H · ∂j wk + ∂kk Hwj − (∂j H ) · (∂k H ) · wk ]. h2
(34)
By substituting the vector operators of Eqs. (33) and (34) in Eq. (28), Naveir’s system of equation in {x} has the form: (Lu)j + ρω2 uj = fj , where (Lu)j =
1 {µ∇ 2 uj + (λ + µ)∂jk2 uk + µ∂k H · ∂k uj − (λ + 3µ)∂j H · ∂k uk h2 +2(λ + 2µ)∂k H · ∂j uk + K0 uj + Kjk · uk }.
2 H µ[∂11
2 H] ∂22
+ and Kjk = (2λ + Let K0 = operator L in the above equation are
3µ)∂jk2 H
(35)
− 2(λ + 2µ)(∂j H ) · (∂k H ). Then, the two components of
1 2 2 2 2 [(λ + 2µ)(∂11 u1 + ∂22 u1 ) − (λ + µ)(∂22 u1 − ∂12 u2 ) + (λ + 2µ)∂1 H · ∂1 u1 + µ∂2 H · ∂2 u1 h2 +2(λ + 2µ)∂2 H · ∂1 u2 − (λ + 3µ)∂1 H · ∂2 u2 + (K0 + K11 )u1 + K12 u2 ] + ρω2 u1 = f1 , 1 2 2 2 2 u2 + ∂22 u2 ) − (λ + µ)(∂11 u2 − ∂21 u1 ) + (λ + 2µ)∂2 H · ∂2 u2 + µ∂1 H · ∂1 u2 (Lu)2 ≡ 2 [(λ + 2µ)(∂11 h +2(λ + 2µ)∂1 H · ∂2 u1 − (λ + 3µ)∂2 H · ∂1 u1 + K21 u1 + (K0 + K22 )u2 ] + ρω2 u2 = f2 . (36) (Lu)1 ≡
As with the 3D case, the Kj k terms in the transformed operator correspond to spatial curvature, while the first order derivatives of the dependent variable indicate dispersion phenomena. Using matrix notation, the Navier’s equations are re-written as 1 1 2 2 2 2 [(λ + 2µ)(∂11 + ∂22 ) − (λ + µ)∂22 [(λ + µ)∂ + 2(λ + 2µ)∂ H · ∂ 2 1 12 2 l h h2 +(λ + 2µ)∂ H · ∂ + µ∂ H · ∂ + K + K ] + ρω2 −(λ + 3µ)∂1 H · ∂2 + K12 ] 1 1 2 2 0 11 1 2 2 2 1 [(λ + 2µ)(∂ + ∂ ) − (λ + µ)∂ 11 22 11 2 h2 2 [(λ + µ)∂21 + 2(λ + 2µ)∂1 H · ∂2 h +(λ + 2µ)∂2 H · ∂2 + µ∂1 H · ∂1 −(λ + 3µ)∂2 H · ∂1 + K21 ] +K0 + K22 ] + ρω2 u1 f1 × = . (37) u2 f2
196
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
Remark 2. It is possible to represent the transform Navier system in unified form for both 2D and 3D cases. Let In (n = 2, 3) be a conformal mapping in Rn ; then, the transform Navier operator Ln acting on vector-valued functions u(n ) has the general form (Ln u(n) )j =
1 (n) (n) (n) (n) {µ∇n2 uj + (λ + µ)∂jk2 uk + µ∂k Hn · ∂k uj + (n − 1)(λ + 2µ)∂k Hn · ∂j uk h2n (n)
(n)
2 −(λ + 3µ)∂j Hn · ∂k uk + µ∂kk Hn uj + [((n − 1)λ + (2n − 3)µ)∂jk2 Hn (n)
−(n − 1)(λ + 2µ)∂j Hn · ∂k Hn ]uk }.
(38)
In the above, Hn = ln|hn |, hn = |J n |1/n and J n is the Jacobian corresponding to conformal mapping In . The differential operators that appear in the Navier’s equations are as follows: (div gradn u(n) )j =
1 2 (n) (n) (n) {∂ u + (n − 1)∂k Hn · ∂j uk − ∂j Hn · ∂k uk + (n − 1)[∂jk2 Hn h2n jk k (n)
−(n − 1)∂j Hn · ∂k Hn ]uk }, (curl2n u(n) )j =
1 (n) (n) (n) (n) (n) (n) 2 {−∇n2 uj + ∂jk2 uk − ∂k Hn · ∂k uj + ∂j Hn · ∂k uk + ∂jk2 Hn uk − ∂kk Hn uj }, h2n
(∇n2 u(n) )j = (div gradn u(n) )j − (curl2n u(n) )j =
1 (n) (n) (n) (n) {∇ 2 u − 2∂j Hn · ∂k uk + (n − 1)∂k Hn · ∂j uk + ∂k Hn · ∂k uj h2n n j (n)
(n)
2 +[(n − 2)∂jk2 Hn − (n − 1)∂j Hn · ∂k Hn ]uk + ∂kk Hn · uj }.
(39)
5.1. Examples for 2D conformal mappings Let I be a conformal mapping in 2D defined by an analytical, complex function f (z) : B → C. This implies that f is one-to-one and regular in B, i.e., f (z) = 0 for z ∈ B. If we define f = x1 (X) + ix2 (X), then x = x(X), J = (∂x/∂X) and det J 2 (X0 ) = |f (z0 )|, where z0 = X10 + iX20 . Conformal mappings preserve the angle between any two curves that pass through point z0 , while the dilatation (expansion) is equal to |f (z0 )| in every direction [18]. The latter property means that as z → z0 locally, we have |f (z) − f (z0 )| = |z − z0 | · |f (z0 )|, which implies (x1 − x10 )2 + (x2 − x20 )2 = cases of conformal mappings:
(X1 − X10 )2 + (X2 − X20 )2 |J(X 0 )|2 . In what follows, we will present two basic
(a) Exponential mapping: Let mapping I be defined by function f (z) = ez = eX1 (cos X2 + i sin X2 ), which is conformal in B = {(X1 , X2 ) : X1 ∈ R1 , X2 ∈ (−π/2, π/2)} and I(B) = R2 \{(x1 , 0) : x1 > 0}, as shown X1 X1 in Fig. 2(a). The direct transformation is x1 = e cos X2 , x2 = e sin X2 , and the inverse transformation is X1 = ln x12 + x22 , X2 = arctan x2 /x1 . Also, the Jacobian matrix is eX1 cos X2 −eX1 sin X2 x1 −x2 J= = , x2 x1 eX1 sin X2 eX1 cos X2 its determinant is |J| = x12 + x22 , the metric tensor is 2 x1 + x22 0 T g = JJ = , 0 x12 + x22
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
197
Fig. 2. Domain of definition for conformal mapping in the complex plane; mappings defined by: (a) f (z) = exp(z); (b) f (z) = z2 .
and finally h = x12 + x22 and H = 0.5 ln(x12 + x22 ). If the radial distance in the transformed coordinate system is r = x12 + x22 and its derivatives with respect to the spatial coordinates are r,j = xj /r, then the coefficients multiplying the dependent variable uk and its derivatives ∂j uk in Eq. (37) are given below as 1 1 ∂jk2 H = 2 (δjk − 2r,j · r,k ), r,j , r r 1 Kjk = 2 [(2λ + 3µ)δjk − (6λ + 10µ)r,j · r,k ]. r
∂j H =
K0 = 0, (40)
(b) Quadratic mapping: In this case, let mapping I be defined by f (z) = z2 = X12 − X22 + i2X1 X2 . As shown in Fig. 2(b), I is conformal in B = {(X1 , X2 ) : X1 ∈ R1 , X2 > 0} and I(B) = R2 \{(x1 , 0) : x1 > 0}. The transformation is x1 = X12 − X22 , x2 = 2X1 X2 , the inverse transformation is −x + x 2 + x 2 1 x2 1 2 , X2 = X1 = 2 2(−x1 + x12 + x22 ) and the remaining properties concerning I are 2X1 −2X2 , |J| = 4 x12 + x22 , J= 2X1 2X2
h = 2(x12 + x22 )1/4 ,
and
H = 0.25 ln(x12 + x22 ) + ln 2. As before, the coefficients multiplying uk and ∂j uk in the transform system of equation (37) are given as 1 1 r,j , ∂jk2 H = 2 (δjk − 2r,j · r,k ), 2r 2r 1 Kjk = 2 [(2λ + 3µ)δjk − (6λ + 10µ)r,j · r,k ]. 2r ∂j H =
K0 = 0, (41)
198
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
6. Numerical examples In this section, we present numerical results for the anti-plane strain case under two conformal mappings, namely the linear and the inverse, which were outlined in Sections 4.1 and 4.2, respectively. The geological medium through which the wave signals propagate has the following background (or reference) material properties: λ0 = 5.0 × 109 N/m2 , cp0 = 2.0 km/s,
µ0 = 2.5 × 109 N/m2 ,
ρ0 = 2500 kg/m3 ,
cs0 = 1.0 km/s,
(42)
where cp0 and cs 0 are the pressure (P) and shear (S) wave velocities. We will further assume that the aforementioned waves travel at an intermediate frequency of ω = 2.0 Hz, which yields the following values for the P and S wave number and wavelength: kp0 =
ω = 6.28 km−1 , cp0
λp0 =
2π = 1.0 km, kp0
ks0 =
ω = 12.56 km−1 , cs0
λs0 =
2π = 0.5 km. ks0 (43)
Fig. 3. Wave amplitude |U3 | versus X1 , X2 for the following linear mapping cases: (a) homogeneous material: h0 = 1, ϕ = 0◦ ; (b) inhomogeneous material: h0 = 2, ϕ = 0◦ ; (c) inhomogeneous material: h0 = 1; ϕ = 45◦ .
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
199
Example 1. The linear conformal mapping, which is a superposition of rotation (with angle of rotation ϕ) and dilatation (with stretching factor h0 ), yields the following transformation from {X} to {x}: x1 = h0 (X1 cos ϕ − X2 sin ϕ),
x2 = h0 (X1 sin ϕ + X2 cos ϕ),
x3 = h0 X3 .
(44)
The solution for U3 (see Eq. (25)) in terms of the modified Bessel function of the second kind [19] for f3 = δ(x) at the origin of the coordinate system is
1 2 2 2 U3 = − K0 −iks0 h0 x1 + x2 , (45) 2π where ks0 h20 is the effective wave number. Fig. 3 plots |U3 | for an {X1 , X2 } grid of dimensions comparable to the S wave length (i.e., 0.5 × 0.5 km2 ) and at the x3 = X3 = 0 level. All plots here have been obtained by using the Mathematica package [20]. Specifically, Fig. 3(a) is for the homogeneous case, which results when h0 = 1, ϕ = 0◦ , while Fig. 3(b) and (c), respectively, plot the same information, but for the variable parameter cases that result when h0 = 2, ϕ = 0◦ and h0 = 1, ϕ = 45◦ . We observe a change in magnitude in |U3 | when h0 = 2 due to the ‘stretching’ of coordinates, while pure rotation (ϕ = 45◦ ) produces no noticeable effect.
Fig. 4. Wave amplitude |U3 | for the inverse mapping case versus X1 , X2 at the following X3 levels: (a) 0.01 km; (b) 0.25 km; (c) 0.50 km.
200
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
Example 2. For this case, the inverse mapping x = X/|X|2 yields the following values for the curvature terms Kˆ 13 and Kˆ 23 in Eq. (21) and for the wave number kˆ0 in Eq. (23): 2µ 1 2µ 1 r,1 , Kˆ 23 = r,2 , λ+µr λ+µr
2 ˆk0 = ρω · 1 1 − 4µ [(r,1 )2 + (r,2 )2 ] − 2 1 + 2λ + 3µ + 2(r,3 )2 . µ r4 λ+µ µ Kˆ 13 =
(46)
We note that the effective wave number in this coordinate system has units of length. If kˆ0 is to remain constant (equal to ks 0 ), then the density profile ρ(x) is equal to
µ 2 4µ 2λ + 3µ ρ(x) = k02 r 2 − r . (47) [(r,1 )2 + (r,2 )2 ] + 2 1 + + 2(r,3 )2 λ+µ µ ω2 We may assume for the purpose of this example that the shear modulus µ remains constant, although it is obvious that it is possible to satisfy Eq. (46) by varying both ρ and µ. The solution |U3 | for a point load at the origin has the same form as that given by Eq. (45), but with the modified wave number kˆ0 replacing the effective wave number and with the inverse transformation now defined as xj = Xj /|X|2 . Fig. 4(a)–(c) plot the amplitude of U3 for the same X3 , X2 spatial grid as before, and at three levels of X3 , namely X3 = 0.01, 0.25 and 0.50 km, respectively. We see a very unusual, almost flat profile in Fig. 4(a), while in Fig. 4(c) the profile approaches that of a homogeneous material. The second profile can be viewed as an intermediate one. Concurrently plotted in Fig. 5(a)–(c) are the corresponding density profiles, where it is observed that the effect of inhomogeneity is less pronounced the farther
Fig. 5. Variables density profiles for the inverse mapping case at the following X3 levels: (a) 0.01 km; (b) 0.25 km; (c) 0.50 km.
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
201
away one moves form the (X1 , X2 ) reference plane. Note that these densities are normalized with respect to their peak values. Specifically, in Fig. 5(a) the material becomes denser around the origin, while in Fig. 5(c) the density is almost uniform. As before, Fig. 5(b) plots an intermediate situation.
7. Conclusions In this work, we have shown the applicability of conformal mapping as a method for obtaining fundamental solutions to both 2D and 3D problems in elastodynamics that involve non-homogeneous elastic media. Conformal mapping is chosen here because the underlying Cauchy–Riemann conditions act as constraints, limiting on one hand the number of possible mappings but, on the other hand, yielding transformed vector operators in the new coordinate system which are simpler in form than they would otherwise be. In earlier work, we showed that for 1D cases the transformed operator obtained by conformal mapping was essentially a scaled version of the original one, so that it was possible to obtain a Green’s function by inspection. The inverse transformation, which completed the solution, was an algebraic one and thus simple to perform since it involved a change of coordinates. Also, the density profile that was recovered depended on the conformal mapping used and was physically plausible. As demonstrated here, all the above advantages of conformal mapping are retained for the vector wave (Navier) equations of elastodynamics, despite the fact that the ensuing coordinate transformations involved a change of the basis vectors. Furthermore, it is possible to extend the conformal mapping method the 3D case, despite the underlying complex variable formalism. Finally, some specific cases were examined in the context of anti-plane as well as of plane strain elastodynamics and the necessary conditions for obtaining closed-form solutions were identified.
Acknowledgements The authors wish to acknowledge the financial support provided by the NATO Science and Environment Affairs Division through grant No. EST.CLG 977774/2001. References [1] W.M. Ewing, W.S. Jardetsky, F. Press, Elastic Waves in Layered Media, McGraw-Hill, New York, 1957. [2] T.A. Cruse, F.A. Rizzo, A direct formulation and numerical solution of the general transient elastodynamic problem. I, J. Math. Anal. Appl. 22 (1968) 244–259. [3] G.D. Manolis, D.E. Beskos, Boundary Element Methods in Elastodynamics, Allen & Unwin, London, 1988. [4] J. Dominguez, Boundary Elements in Dynamics, Computational Mechanics Publications, Southampton, 1993. [5] W.C. Chew, Waves and Fields in Inhomogeneous Media, Van Nostrand Reinhold, New York, 1990. [6] K. Helbig (Ed.), Modeling the Earth for Oil Exploration, Pergamon/Elsevier, Oxford/Amsterdam, 1993. [7] A.H. Nayfeh, Wave Propagation in Layered Anisotropic Media, North-Holland, Amsterdam, 1995. [8] H. Sato, M.C. Fehler, Wave Propagation and Scattering in the Heterogeneous Earth, Springer, New York, 1997. [9] R.P. Shaw, G.D. Manolis, Conformal mapping solutions for 2D heterogeneous Helmholtz equations, in: M. Marchetti, C.A. Brebbia, H.H. Aliabadi (Eds.), Proceedings of the 19th BEM International Conference, Computational Mechanics Publications, Southampton, 1997, pp. 411–418. [10] G.D. Manolis, R.P. Shaw, Green’s function for the vector wave equation in a mildly heterogeneous medium, Wave Motion 24 (1996) 59–83. [11] J.D. Achenbach, Wave Propagation in Elastic Solids, North-Holland, Amsterdam, 1973. [12] J. Gerretsen, Lectures on Tensor Calculus and Differential Geometry, Noordhoff, Groningen, 1962. [13] P.M. Morse, H. Feshbach, Methods of Theoretical Physics, Vol. 1, McGraw-Hill, New York, 1953. [14] G.D. Manolis, R.P. Shaw, Fundamental solutions for variable density two-dimensional elastodynamic problems, Eng. Anal. Bound. Elements 24 (10) (2000) 739–750. [15] B. Dubrovin, S. Novikov, A. Fomenko, Modern Geometry: Methods and Applications, Nauka, Moscow, 1979 (in Russian). [16] V.S. Vladimirov, Equations of Mathematical Physics, Mir, Moscow, 1984.
202
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
[17] G.D. Manolis, R.P. Shaw, Wave motions in stochastic heterogeneous media: a Green’s function approach, Eng. Anal. Bound. Elements 15 (1995) 225–234. [18] Y. Sidorov, M. Fedoriuk, M. Shabunin, Lectures on the Theory of Complex Functions, Nauka, Moscow, 1982 (in Russian). [19] I.S. Gradsteyn, I.M. Ryzhik, Table of Integrals, Series and Products, Academic Press, New York, 1980. [20] S. Wolfram, A System for Doing Mathematics by Computer, Version 3.0, Cambridge University Press, Cambridge, 1996.
Conformal mapping methods for variable parameter elastodynamics George D. Manolis a,∗ , Tsviatko V. Rangelov b , Richard Paul Shaw c b
a Department of Civil Engineering, Aristotle University, Thessaloniki GR-54006, Greece Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, Sofia 1113, Bulgaria c Department of Civil Engineering, State University of New York, Buffalo, NY 14260, USA
Received 4 December 2001; accepted 15 January 2002
Abstract This work is a generalization of earlier efforts to develop free-space fundamental solutions (Green’s functions) for a certain class of non-homogenous materials (i.e., materials with position-dependent elastic parameters and density) under time-harmonic conditions by means of conformal mapping methods. These methods were shown to be very effective in conjunction with the scalar wave equation, but their extension to the 2D and 3D vector wave equations of elastodynamics is hampered by the following two factors: (a) the governing equations of motion are vectorial, which implies that any transformation of coordinates affects the basis vectors used in the original (i.e., Cartesian) coordinate system; and (b) an extension of conformal mapping to a 3D coordinate system is not readily obvious. In the present work, both issues raised above are addressed. Furthermore, the necessary background is established whereby the construction of Green’s functions for 2D and 3D elastodynamics can be systematically accomplished in a mapped space, followed by an inversion to the original Cartesian coordinate system. Of course, conformal mapping as used here is inherently an inverse method, in the sense that the material parameter profiles for a given problem are recovered as constraints during the solution procedure and therefore depend on the type of mapping used. © 2002 Elsevier Science B.V. All rights reserved.
1. Introduction The development of fundamental solutions for elastic wave propagation in solids that are not homogeneous and isotropic, but instead exhibit heterogeneity, anisotropy, layering, randomness, etc., are important in their own right [1]. At the same time, fundamental solutions for the free-space (known as Green’s functions) comprise an essential part of any boundary integral equation method (BIEM) formulation [2] and, of course, of its corresponding numerical solution technique, the boundary element method (BEM). The latter method has been very successful in solving boundary-value problems of engineering importance in transient elastodynamics [3,4]. Since the preferred formulation is invariably a displacement-traction approach based on Somigliana’s identity, the corresponding Green’s functions must come from a solution of the equations of elastodynamics for a point force and under radiation-type boundary conditions. This can be done by a variety of ways, e.g., use of potentials (Helmholtz potentials, Stokes potentials), displacement vector decomposition into dilatational and rotational components, use of the dynamic equivalent to Galerkin’s vector, integral transforms or other specialized approaches. ∗
Corresponding author. Tel.: +30310-995-663; fax: +30310-995-769. E-mail address: [email protected] (G.D. Manolis). 0165-2125/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 2 1 2 5 ( 0 2 ) 0 0 0 1 1 - 2
186
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
Any extension of the aforementioned methodologies to problems involving other categories of media, which go beyond the homogeneous, isotropic and linear elastic continuum, has proven difficult. This, however, would be very useful in various engineering fields, such as seismology, earthquake engineering, non-destructive testing evaluation, acoustics, electromagnetic signal transmission, ocean engineering, etc., where the materials through which various types of signals propagate cannot be satisfactorily modeled as homogeneous [5–8]. In this respect, we briefly mention some recent activity in the area of time-harmonic acoustic waves [9] and elastic waves [10] in heterogeneous media using a coordinate transformation based on conformal mapping. In the latter case of elastic waves, which are governed by vector differential equations, the class of solvable problems was limited to variable (i.e., position-dependent) density and constant elastic parameters. Even then, an additional assumption regarding decomposition of the displacement vector into pseudo-dilatational and pseudo-rotational parts was necessary in order to achieve a closed-form solution through reduction to the anti-plane strain case. While this approach at first sight appears to be limited to 1D (or, at best, to 2D) problems and is inverse in form, i.e., a given mapping will lead to a particular type of heterogeneity, it represents a step forward in the development of Green’s functions of the type required by the BIEM and BEM. In this work, we examine in detail the coordinate transformations implied through conformal mapping for the vector wave equations of elastodynamics under both 3D and plane strain conditions, whereby the latter case is proved to be a special case of the former. Furthermore, we seek to establish the conditions under which it is possible to obtain a solution in the transformed, orthogonal curvilinear coordinate system into which the original (Cartesian) coordinates are mapped. The key advantage of using what are essentially coordinate transformations is that the inverse mapping, which is required to give the final Green’s function in Cartesian coordinates, is straightforward to perform. Although the variable density case represents the simplest material heterogeneity, it is possible to look at the general case as well, where all material parameters of the medium in question vary with position.
2. The governing equations of elastodynamics in 3D Navier’s equations [11] of dynamic equilibrium for time-harmonic elastodynamics written in the reference Cartesian coordinate system {X} are L(U) + ρω2 U ≡ (λ + 2µ)∇∇ · U − µ∇ × ∇ × U + ρω2 U = F ,
(1)
where λ, µ are the Lame constants, ρ is the mass density, and ω the frequency of vibration. Furthermore, U = (U1 , U2 , U3 ) is the displacement vector, F the body force per unit mass, and the two vector operators used in the above formulation are ∇∇ · U = grad div U and ∇ × ∇ × U = curl curl U. Let I map {X} into an orthogonal curvilinear coordinate system {x}, where xj = xj (X), with index j = 1, 2, 3 as shown in Fig. 1. The Jacobian matrix of the transformation is defined as J = (∂x/∂X), provided that det J = |J| = 0. Also, g = JJT is the metric tensor in the new coordinate system {x} given as h1 0 0 g = 0 h2 0 0 0 h3 with hj > 0. If I is a conformal mapping, then h1 = h2 = h3 = h2 , where h = |J|1/3 [12]. Furthermore, we need to define as {Ei } the orthogonal basis in the reference system {X}, which maps onto the curvilinear orthogonal basis {ei } of {x} through I. Obviously, what complicates transformations between coordinate systems is that the base vectors are dependent on the spatial coordinates, which implies use of the chain rule for differentiation of vector functions [13]. The transformed Navier’s equations in the new coordinates {x} are obtained using the appropriate representations for the differential operators acting on vector functions. Specifically, let the displacement vector be expressed either
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
187
Fig. 1. Conformal mapping between two orthogonal curvilinear coordinate systems in: (a) 3D space; (b) 2D space.
in terms of the original coordinate system {X} as U = U1 E1 + U2 E2 + U3 E3 or in terms of the transformed coordinate system {x} as u = u1 e1 + u2 e2 + u3 e3 . If a solution can be obtained in the transformed coordinates, then the solution of Eq. (1) is simply recovered by applying the inverse transformation I−1 so that U(X) = u(x(X)). There is an inherent duality in these transformations, because although the {X} space is originally viewed as Cartesian and the {x} as orthogonal curvilinear, any fundamental solution obtained in the latter space can be treated as the Cartesian coordinate system solution. Then, the inverse mapping can be used to return this solution back to the {X} coordinates, which now comprise an orthogonal curvilinear system as seen from {x}. In what follows, we list the transformation of all key operators, noting that the summation convention is implied for repeated indices. Case 1: The gradient operator for a scalar displacement function is ∇U = Ej ·
∂U ∂Xj
in {X},
1 ∂u ej · h ∂xj
in {x}.
(2)
1 ∂(uj h2 ) h3 ∂xj
in {x}.
(3)
∇u =
Case 2: The divergence operator is ∇ ·U =
∂Uj ∂Xj
in {X},
∇ ·u=
Case 3: The curl operator is E1 E2 E3 ∂ ∂ ∂ ∇ × U = ∂X ∂X ∂X 1 2 3 U U U 1 2 3
in {X},
1 ∇ × u = 3 h
he1 ∂ ∂x1 hu1
he2 ∂ ∂x2 hu2
he3 ∂ ∂x3 hu3
in {x}.
(4)
Case 4: The grad div operator is computed by applying the second equation (2) to the scalar function defined by the second equation (3). Specifically, we first recover the intermediate result ∇ · u = (1/ h3 )∂k (uk h2 ) = (1/ h)(∂k uk + 2∂k H · uk ), where the notation ∂/∂xj = ∂j and ln|h| = H has been introduced. Finally, 1 1 ∇∇ · u = ej ∂j (∂k uk + 2∂k H · uk ) h h 1 = 2 ej [∂jk2 uk − ∂j H · ∂k uk + 2∂k H · ∂j uk + (−2∂j H · ∂k H + 2∂jk2 H )uk ]. (5) h
188
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
Case 5: The curl curl operator is computed by applying the second equation (4) to the vector function curl u. Since the cross product yields ∇ × u = (1/ h2 )ej wj , where the three components are w1 = (1/ h2 )[∂2 (hu3 ) − ∂3 (hu2 )], w2 = (1/ h2 )[∂3 (hu1 ) − ∂1 (hu3 )] and w3 = (1/ h2 )[∂1 (hu2 ) − ∂2 (hu1 )], then he1 he2 he3 1 ∂ ∂ ∂ ∇ × ∇ × u = 3 h ∂x1 ∂x2 ∂x3 hw hw hw 1 2 3 1 = 2 {e1 [∂2 (hw3 ) − ∂3 (hw2 )] + e2 [∂3 (hw1 ) − ∂1 (hw3 )] + e3 [∂1 (hw2 ) − ∂2 (hw1 )]}. (6) h The coefficient multiplying vector component e1 in the above equation is given by the following expression: 1 1 1 1 1 [∂2 (hw3 ) − ∂3 (hw2 )] = 2 ∂2 h · w3 + ∂2 w3 − 2 ∂3 h · w2 − ∂3 w2 = W11 + W12 + W13 + W14 . (7) h h h2 h h Carrying through with the differentiations yields 1 1 1 1 ∂2 h · 2 [∂1 (hu2 ) − ∂2 (hu1 )] 2 ∂2 H · [∂1 h · u2 + h∂1 u2 − ∂2 h · u1 − h∂2 u1 ] h h2 h h 1 = 2 [−∂2 H · ∂2 u1 + ∂2 H · ∂1 u2 − (∂2 H )2 u1 + ∂1 H · ∂2 H · u2 ], h
1 1 W12 = 2 h∂2 2 [∂1 (hu2 ) − ∂2 (hu1 )] h h
1 2 1 = − 3 ∂2 h[∂1 (hu2 ) − ∂2 (hu1 )] + 2 ∂2 [∂1 h · u2 + h∂1 u2 − ∂2 h · u1 − h∂2 u1 ] h h h 1 = 2 [−2∂2 H (∂1 H · u2 + ∂1 u2 − ∂2 H · u1 − ∂2 u1 )] h
2 h 2 h ∂22 1 ∂21 2 2 + 2 · u2 + ∂1 H ∂2 u2 + ∂2 H ∂1 u2 + ∂12 u2 − · u1 − ∂2 H ∂2 u1 − ∂2 H ∂2 u1 − ∂22 u1 h h h W11 =
=
1 2 2 2 (−∂22 u1 + ∂12 u2 − ∂2 H · ∂1 u2 + ∂1 H · ∂2 u2 − [∂22 H − (∂2 H )2 ] · u1 h2 2 +[∂21 H − (∂1 H.∂2 H )] · u2 ),
(8)
where the equality ∂ij2 h/ h = ∂ij2 H + ∂i H · ∂j H has been used. The remaining two parts are obtained in a similar fashion, i.e., 1 1 1 ∂3 h · 2 [∂3 (hu1 ) − ∂1 (hu3 )] = 2 [−∂3 H · ∂3 u1 + ∂3 H · ∂1 u3 − (∂3 H )2 u1 + ∂1 H · ∂3 H · u3 ], 2 h
h h 1 1 W14 = − 2 h∂3 2 [∂3 (hu1 ) − ∂1 (hu3 )] h h 1 2 2 2 u1 + ∂13 u3 + ∂1 H · ∂3 u3 − ∂3 H · ∂1 u3 − [∂33 H − (∂3 H )2 ] · u1 = 2 (−∂33 h 2 +[∂13 H − (∂1 H ) · (∂3 H )] · u3 ). (9) W13 = −
Since we have evaluated the first component of the curl curl operator, namely (∇×∇×u)1 = W11 +W12 +W13 +W14 , the remaining two components are derived by cyclic interchange of the indices as (∇ × ∇ × u)j =
1 2 {−∇ 2 uj + ∂kj2 uk − ∂k H · ∂k uj + ∂j H · ∂k uk + ∂jk2 H · uk − ∂kk H · uj }. h2
(10)
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
189
Case 6: The Laplacian ∇ 2 in the transformed coordinates {x} can be computed from the grad div and the curl curl operators as follows: ∇ 2 u = ∇∇ · u − ∇ × ∇ × u.
(11)
Thus, the jth component (j = 1, 2, 3) is equal to (∇ 2 u)j =
1 2 {∇ 2 uj + ∂k H · ∂k uj + 2∂k H · ∂j uk − 2∂j H · ∂k uk + ∂kk H · uj + [∂jk2 H − 2(∂k H )(∂j H )]uk }. h2 (12)
2.1. The transformed Navier’s equations By combining all the previous results, the components of the differential operator L that appeared in Eq. (1) are now (Lu)j =
1 {µ∇ 2 uj + (λ + µ)∂jk2 uk + µ∂k H · ∂k uj + 2(λ + 2µ)∂k H · ∂j uk − (λ + 3µ)∂j H · ∂k uk h2 2 +µ∂kk Huj + [(2λ + 3µ)∂jk2 H − 2(λ + 2µ)∂j H · ∂k H ]uk }.
(13)
The first two terms in the above equation comprise Lame’s operator, but we see the emergence of additional terms that can be viewed as the influence of inhomogeneities in the original material. The key objective is to try and match these extra terms to a particular density and/or elastic parameter profile. For instance, in previous work [14], it was possible to recover the standard form of Navier’s equation in the {x} coordinates by prescribing a density profile equal to ρ(x, y) = ρ0 |J(x, y)|. Navier’s system of equations in the transformed coordinate system {x} is now (Lu)j + ρω2 uj = fj ,
(14)
where (Lu)1 ≡
1 2 2 2 2 2 {(λ + 2µ)∂11 u1 + µ(∂22 u1 + ∂33 u1 ) + (λ + µ)(∂12 u2 + ∂13 u3 ) + (λ + 2µ)∂1 H · ∂1 u1 h2 +µ(∂2 H · ∂2 u1 + ∂3 H · ∂3 u1 ) + 2(λ + 2µ)(∂2 H · ∂1 u2 + ∂3 H · ∂1 u3 )
−(λ + 3µ)(∂2 H · ∂2 u2 + ∂3 H · ∂3 u3 ) + (K10 + K11 )u1 + K12 u2 + K13 u3 }, 1 2 2 2 2 2 (Lu)2 ≡ 2 {(λ + 2µ)∂22 u2 + µ(∂11 u2 + ∂33 u2 ) + (λ + µ)(∂21 u1 + ∂23 u3 ) + (λ + 2µ)∂2 H · ∂2 u2 h +µ(∂1 H · ∂1 u2 + ∂3 H · ∂3 u2 ) + 2(λ + 2µ)(∂1 H · ∂2 u1 + ∂3 H · ∂2 u3 ) −(λ + 3µ)(∂1 H · ∂1 u1 + ∂3 H · ∂3 u3 ) + K21 u1 + (K20 + K22 )u2 + K23 u3 }, 1 2 2 2 2 2 (Lu)3 ≡ 2 {(λ + 2µ)∂33 u3 + µ(∂11 u3 + ∂22 u3 ) + (λ + µ)(∂31 u1 + ∂32 u2 ) + (λ + 2µ)∂3 H · ∂3 u3 h +µ(∂1 H · ∂1 u3 + ∂2 H · ∂2 u3 ) + 2(λ + 2µ)(∂1 H · ∂3 u1 + ∂2 H · ∂3 u2 ) −(λ + 3µ)(∂1 H · ∂1 u1 + ∂2 H · ∂2 u2 ) + K31 u1 + K32 u2 + (K30 + K33 )u3 }.
(15)
2 H + ∂ 2 H + ∂ 2 H ] and K = (2λ + 3µ)∂ 2 H − 2(λ + In the above equations, we have that Kj 0 = µ[∂11 jk 22 33 jk 2µ)(∂j H )(∂k H ), where j, k = 1, 2, 3. Note at this point that index j in the first expression is superfluous, since K10 = K20 = K30 , but is retained for notational convenience. From a mechanical point of view, the Kj k terms in the transformed equations correspond to spatial curvature, while the first order derivatives of the dependent variable u indicate dispersion phenomena. The effect of the curvature terms is localized and, from a practical viewpoint, it is possible to approximate their influence by converting the original wave numbers into complex quantities. The dispersive terms, on the other hand, produce damping that is responsible for signal de-amplification and dispersion.
190
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
This effect can noticeably alter the propagating signal, especially as the distance from source increases and at high frequencies of vibration. 2.2. Particular cases In this section, we will examine some special cases concerning the Navier’s equations in the transformed coordinate system {x}. Case 1: We assume here that each displacement component depends on one and only one spatial variable. Ordinarily, this would lead to three separate 1D equations of motion, but as will be shown here, there is coupling due to the spatial curvature and the first order terms mentioned above. Specifically, for u = (u1 , u2 , u3 ), we have that u1 = u1 (x1 ), u2 = u2 (x2 ), u3 = u3 (x3 ) and any solutions in this form must satisfy the following system: 1 2 {(λ + 2µ)∂11 u1 + (λ + 2µ)∂1 H · ∂1 u1 − (λ + 3µ)(∂2 H · ∂2 u2 + ∂3 H · ∂3 u3 ) h2 +(K10 + K11 )u1 + K12 u2 + K13 u3 } + ρω2 u1 = f1 , 1 2 (Lu)2 + ρω2 u2 ≡ 2 {(λ + 2µ)∂22 u2 + (λ + 2µ)∂2 H · ∂2 u2 − (λ + 3µ)(∂1 H · ∂1 u1 + ∂3 H · ∂3 u3 ) h +K21 u1 + (K20 + K22 )u2 + K23 u3 } + ρω2 u2 = f2 , 1 2 (Lu)3 + ρω2 u3 ≡ 2 {(λ + 2µ)∂33 u3 + (λ + 2µ)∂3 H · ∂3 u3 − (λ + 3µ)(∂1 H · ∂1 u1 + ∂2 H · ∂2 u2 ) h +K31 u1 + K32 u2 + (K30 + K33 )u3 } + ρω2 u3 = f3 . (Lu)1 + ρω2 u1 ≡
(16)
Case 2: This parallels the standard plane strain conditions of classical elastodynamics. The relevant assumptions are well known, i.e., we have that u1 = u1 (x1 , x2 ), u2 = u2 (x1 , x2 ), u3 = u3 (x3 ). Solutions in this form must satisfy the following system of coupled equations: (Lu)1 + ρω2 u1 ≡
1 2 2 2 {(λ + 2µ)∂11 u1 + µ∂22 u1 + (λ + µ)∂12 u2 + (λ + 2µ)∂1 H · ∂1 u1 + µ∂2 H · ∂2 u1 h2 +2(λ + 2µ)∂2 H · ∂1 u2 − (λ + 3µ)(∂2 H · ∂2 u2 + ∂3 H · ∂3 u3 ) + (K10 + K11 )u1
+K12 u2 + K13 u3 } + ρω2 u1 = f1 , 1 2 2 2 (Lu)2 + ρω2 u2 ≡ 2 {(λ + 2µ)∂22 u2 + µ∂11 u2 + (λ + µ)∂21 u1 + (λ + 2µ)∂2 H · ∂2 u + µ∂1 H · ∂1 u2 h +2(λ + 2µ)∂1 H · ∂2 u1 − (λ + 3µ)(∂1 H · ∂1 u1 + ∂3 H · ∂3 u3 ) + K21 u1 +(K20 + K22 )u2 + K23 u3 } + ρω2 u2 = f2 , 1 2 u3 + (λ + 2µ)∂3 H · ∂3 u3 − (λ + 3µ)(∂1 H · ∂1 u1 + ∂2 H · ∂2 u2 ) (Lu)3 + ρω2 u3 ≡ 2 {(λ + 2µ)∂33 h +K31 u1 + K32 u2 + (K30 + K33 )u3 } + ρω2 u3 = f3 .
(17)
Case 3: This last case parallels the anti-plane strain condition of elastodynamics, whereby the displacement vector components are u1 = 0, u2 = 0, u3 = u3 (x1 , x2 ). Then, any solution must satisfy the system 1 {(λ + µ)∂3 H · ∂1 u3 + K13 u3 } = f1 , h2 1 (Lu)2 + ρω2 u2 ≡ 2 {(λ + µ)∂3 H · ∂2 u3 + K23 u3 } = f2 , h 1 2 2 (Lu)3 + ρω2 u3 ≡ 2 {µ(∂11 u3 + ∂22 u3 ) + µ(∂1 H · ∂1 u3 + ∂2 H · ∂2 u3 ) + (K30 + K33 )u3 } + ρω2 u3 = f3 . h (18) (Lu)1 + ρω2 u1 ≡
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
191
3. Conformal mappings in 3D In the previous section, conformal mapping was used to transform Navier’s equations from orthogonal to curvilinear coordinates. In order to quantify the terms h, ∂j H, ∂jk2 H, Kjk that appear in Eq. (15), it is necessary to classify conformal mappings that are possible in 3D [15]. At first, the mapping from {X} to {x} is conformal if the metric G in {X} is proportional to the metric g in {x}, i.e., there exists a function f (x) = 0 such that gαβ = f (x)Gαβ . In our case, 1, α = β, Gαβ = 0, α = β, since it is a Cartesian metric. Thus, f (x), α = β, gαβ = 0, α = β, and f (x) = h2 (x) = |J|1/3 , where J is the Jacobian matrix of the transformation I at x. Next, we list the following three basic cases of conformal mappings in 3D: (a) Rotation: Let x = AX, where A is an orthogonal matrix, i.e., AT = A−1 . Then J = A, g = JJT = E, where E is the identity matrix, h(x) ≡ 1, and finally gαβ = Gαβ . (b) Dilatation: Let x = h0 X, where h0 is a constant; then J = h0 E, h(x) ≡ h0 , and gαβ = h20 Gαβ . Note that both rotation and dilatation are linear transformations. (c) Inversion: Let x = (X − X0 )/|X − X0 |2 , where X0 is a fixed point. Such a transformation is defined in R3 \{X 0 } and h(X) = −|X − X0 |−2 , |J| = −|X − X0 |−6 , and gαβ = |X − X0 |−4 Gαβ . A theorem by Liouville [15] classifies all smooth conformal mappings in 3D as a superposition of rotation, dilatation and inversion. It therefore follows that there can be no conformal transformations in 3D other than those composed from cases (a)–(c), plus rigid body motion. This is in sharp contrast to 2D, where every analytic function in R2 defines a conformal mapping. Next, we present some details regarding the particular form assumed by operator L in Eq. (14) for the aforementioned three basic cases. Specifically, for case (a), transformation I is defined as x = AX, the inverse transformation I−1 is X = A−1 x, and we have that h = 1, H = ln h = 0, ∂j H = 0, Kjk = 0 and L is the same as the operator that appears in Eq. (1). Next, for case (b), transformation I is defined as x = h0 X, the inverse transformation −4 2 I−1 is X = h−1 0 x, h = h0 , H = ln h = 0, ∂j H = 0, Kjk = 0 and L is that of Eq. (1) scaled by h0 . Finally, transformation I for case (c) is x = X/|X|2 , where we let X0 = 0 for simplicity. The inverse transformation I−1 is X = x/|x|2 and h(x) = −1/|x|2 , H = −2 ln |x|. There is no simple relation between the transformed operator L in Eq. (13) and the original one; details for this case will be given in the next section.
4. Solutions for the reduced anti-plane strain case We will pursue some solutions of the system of equations that pertain to the anti-plane strain case as reduced from 3D. Specifically, if the force components in the transformed coordinates f1 , f2 , f3 are continuous, Eq. (18) have a solution u3 locally for the following forms of function H (=ln|h|): (a) H is constant in a neighborhood Ω of field point x0 . (b) The derivative of H with respect to the third spatial coordinate ∂3 H = 0 at x0 and consequently in the neighborhood of that point.
192
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
In the first case, ∂j H = 0 for j = 1, 2, 3 in Ω and the left-hand sides of the first two equations (18) are identically zero, which implies that components f1 , f2 must also be zero in order to have a solution. Then, Eq. (18) are reduced to just the last equation 1 2 2 µ(∂11 u3 + ∂22 u3 ) + ρω2 u3 = f3 . h2
(19)
Let u˜ 3 be a fundamental solution of Eq. (19) for a unit point force (Dirac’s delta function), i.e., 2 2 ∂11 u˜ 3 + ∂22 u˜ 3 + k02 u˜ 3 = δ,
(20) (1)
where k02 = e2H ρω2 /µ is a wave number and e2H = h2 . The fundamental solution is u˜ 3 = −(i/4)H0 (k0 |x|), while that of Eq. (19) is u3 = u˜ 3 ∗ {f3 e2H /µ}, where ∗ is the convolution operator [16]. The actual density profile for which this solution is valid is non-constant and equal to ρ(x) = µ(x){k02 /ω2 e2H }. More specifically, we recover that ρ(x)/µ(x) = ρ0 /µ0 , where ρ 0 , µ0 are reference (or background) values identified with a key point, such as the origin of the coordinate system and that the true wave number is k = k0 / h. A similar profile, whereby the density is proportional to the shear modulus, was derived by Manolis and Shaw [17], who worked directly with the 2D Helmholtz equation. In the second case, the first and second equations (18) yield ∂1 u3 = fˆ1 − Kˆ 13 u3 , ∂2 u3 = fˆ2 − Kˆ 23 u3 ,
h 2 f1 K13 , Kˆ 13 = , (λ + µ)∂3 H (λ + µ)∂3 H h 2 f2 K23 fˆ2 = , Kˆ 23 = , (λ + µ)∂3 H (λ + µ)∂3 H
fˆ1 =
(21)
while the third one is 2 2 ∂11 u3 + ∂22 u3 + kˆ02 u3 = fˆ3 ,
(22)
where the modified forcing function and the modified wave number are equal to h fˆ3 = f3 − ∂1 H · fˆ1 − ∂2 H.fˆ2 , µ 2
ρω2 h2 1 kˆ02 = − ∂1 H · Kˆ 13 − ∂2 H · Kˆ 23 + (K30 + K33 ). µ µ
(23)
Eq. (22) is of the Helmholtz type, and given that u˜ 3 is a fundamental solution for a unit point force, then the complete solution is u3 = u˜ 3 ∗ fˆ3 . We note here that since h(x) is known, density ρ(x) and shear modulus µ(x) in Eq. (23) must combine in such a way so as to give a constant value for kˆ02 . Remark 1. Following the smooth conformal mapping I from {X} to {x} in the anti-plane strain case, Liouville’s theorem implies that at an arbitrary field point x0 , one of the following possibilities must hold true: (i) ∂1 H = ∂2 H = ∂3 H = 0, in which case Eq. (18) is solvable as in case (a). Essentially, the problem uncouples if the forcing function components corresponding to zero displacements are also zero. (ii) ∂j H = 0 for j = 1, 2, 3. Then, there exists a neighborhood Ω of x0 such that ∂j H = 0 and the system of equations (18) are solvable for uj as in case (b). Specifically, if ∂1 H = 0, we have a solution for u1 (x2 , x3 ), u2 = u3 = 0, and if ∂2 H = 0, the solution is for u2 (x1 , x3 ), u1 = u3 = 0. The last case, where ∂3 H = 0, is the same as case (b). In what follows, we will further examine the anti-plane strain case in conjunction with two specific types of conformal mappings, namely the linear (dilatation plus rotation) and the inverse that were presented in Section 3.
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
193
4.1. The anti-plane strain case for a conformal mapping of the linear type Let I be a linear conformal mapping, defined by a constant 3 × 3 Jacobian matrix J, det J = |J| = 0, such that (1/|J|)J is an orthogonal matrix. Therefore, x = JX, |J| = h30 and the metric tensor is h20 0 0 2 g = JJT = 0 h0 0 , 0 0 h20 which implies that h = h0 and H = ln|h| = 0 (h0 is a dimensionless factor). More precisely, the Jacobian matrix is J = [ajk ] and the orthogonality condition can be summarized as 3 h20 , i = j, aik ajk = 0, i = j. k=1 For this type of mapping, Eq. (18) read as (Lu)1 + ρω2 u1 ≡ 0 = f1 , (Lu)2 + ρω2 u2 ≡ 0 = f2 , µ 2 2 (∂11 u3 + ∂22 u3 ) + ρω2 u3 = f3 . (Lu)3 + ρω2 u3 ≡ h20
(24)
We distinguish the following two possibilities: (i) If (f1 )2 + (f2 )2 = 0, then Eq. (24) has no solution. (1) (ii) If√f1 = f2 = 0, the Helmholtz solution previously given as u3 = u˜ 3 ∗ {f3 h20 /µ}, where u˜ 3 = −(i/4)H0 { ρ0 /µ0 h0 ω|x|}, is valid. Therefore, the corresponding solution to Navier’s equations (1) in the original Cartesian coordinate system {X} is U = u(x(X)) with U1 = 0, U2 = 0 and 3 3 U3 = u3 (x1 (X), x2 (X)) = u3 ((JX)1 , (JX)2 ) = u3 a1j Xj , (25) a2j Xj . j =1
j =1
4.2. The anti-plane strain case for a conformal mapping of the inverse type Let conformal mapping I be an inversion defined from R3 \{0} to R3 \{0} as x = X/|X|2 , while I−1 is defined as X = x/|x|2 . If the radial distance in {X} is denoted by R = |X|, then R,j = X,j /R and the elements of the Jacobian matrix J are Jjk = (1/R 2 )(δjk − 2R,j R,k ). Furthermore, det J = |J| = −1/R 6 , h(x) = −r −2 , H = −2 ln|x| and the metric tensor is h2 (x), α = β, gαβ = 0, α = β. For this mapping, the following terms appearing in the transformed Navier’s equations (18) are computed: 2 2 4 2µ ∂j H = − r,j , ∂jk2 H = − 2 (δjk − 2r,j r,k ), (∂j H ) · (∂k H ) = 2 r,j r,k , Kj 0 = − 2 , r r r r 2 (26) Kjk = − 2 [(2λ + 3µ)δjk + 2µr,j r,k ]. r In the above, the radial distance in {x} is defined as r = |x| with r,j = xj /r. As defined, the {x} coordinate system has dimensions of inverse length. We conclude that Eq. (18) are solvable in every domain B ⊂ R3 \{x3 = 0}, since
194
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
∂3 H = −(2/r)r,3 = 0. The solution is u3 = u˜ 3 ∗ fˆ3 , where u˜ 3 is a fundamental solution of Eq. (22) and fˆ3 , kˆ02 are given in Eq. (23). The corresponding solution to Navier’s equations (1) in Cartesian coordinates {X} is X1 X2 U3 = u3 U1 = U2 = 0, . (27) , |X|2 |X|2 5. The governing equations of elastodynamics in 2D In 2D elastodynamics, the key problem that arises as one seeks to define a suitable transformation between the reference Cartesian coordinates {X} and new curvilinear coordinates {x} is that the standard definition for the cross product introduces a third dimension (see Fig. 1). To circumvent this problem, Navier’s equations for the plane strain, time-harmonic case are presented in a form that no longer involves the curl curl operator. Specifically, for a simply connected domain B ⊂ R2 , we have that LU + ρω2 U ≡ µ∇ 2 U + (λ + µ)∇∇ · U + ρω2 U = F ,
(28)
where λ, µ, ρ, ω have all been previously defined, while U = (U1 , U2 ) and F = (F1 , F2 ), respectively, are the displacement and body force per unit mass vectors. The standard definitions for the vector operators are used, i.e., ∂U1 ∂U2 ∂U ∂U , div U = ∇ · U = , + , grad U = ∇U = ∂X1 ∂X2 ∂X1 ∂X2 ∂ 2 U1 ∂ 2 U2 ∂ 2 U1 ∂ 2 U2 grad div U = ∇∇ · U = + , , + ∂X1 ∂X2 ∂X1 ∂X2 ∂X12 ∂X22 ∂ 2 U1 ∂ 2 U 1 ∂ 2 U2 ∂ 2 U2 2 . (29) + , + div grad U = ∇ U = ∂X12 ∂X22 ∂X12 ∂X22 Let I map domain B into a 2D manifold M embedded in R3 . The coordinate system used in B is {X}, while that used in M is {x}, where xj = xj (X). The Jacobian matrix of the transformation is J = (∂x/∂X), where det J = 0. If we define gαβ = JJT as the contravariant metric tensor, then gαβ is the corresponding covariant tensor β and g αβ · gαβ = δα (i.e., the Kronecker delta). If I is an orthogonal mapping, then at every point in M the following conditions are fulfilled: hα , α = β, (hα )−1 , α = β, gαβ = and g αβ = 0, α = β, 0, α = β. Furthermore, if I is conformal, then we have the additional condition h1 = h2 = h2 = |J| at every point in M. As before, we denote by {Ei } the orthogonal basis in {X}, which is mapped into the curvilinear orthogonal basis {ei } in {x}. The objective is to solve Navier’s system of equations in the transformed coordinates and, by applying the inverse mapping I−1 , obtain the solution of Eq. (28) as U(X) = u(x(X)). It is now necessary to re-define the vector operations in Eq. (29) for the transformed coordinate system. In this respect, we follow Ref. [11] and also use tensor notation [12], where the summation rules for rising and lowering indices is implied along with summation for repeated symbols. The relevant expressions are given below for a scalar function ϕ and a vector function w on M (note that the formulas employing tensor notation use Greek symbols, while standard indices are used for vector notation; also, ∂s = ∂/∂xs ). Case 1: Gradient operator for scalar functions: (grad ϕ)α = (∇ϕ)α = g αβ ∂β ϕ,
(grad ϕ)j = (∇ϕ)j =
1 ej ∂j ϕ. h
(30)
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
195
Case 2: Gradient operator for vector functions: (grad w)αβ = (∇w)αβ = g βγ ∂γ wα ,
grad w = ∇w =
1 ej ∂j w. h
(31)
Case 3: Divergence operator: 1 β div w = ∇ · w = √ (∂α w γ + Γνβ w ν ), g
div w = ∇ · w =
1 ∂j (wj h). h2
(32)
k is In the above, wα is a covariant vector, wα is a contravariant vector and w α = g αβ wβ , wδ = gδγ w γ . Also, Γαβ 1 k ks the Christoffel tensor and Γαβ = Γsαβ g , where Γsαβ = 2 (∂β gsα + ∂α gβs − ∂s gβα ). The next step is to recover the grad div and Laplace operators. First, define ln|h| = H and note that ∂ij2 h/ h = ∂ij2 H + ∂i H · ∂j H . Then, by applying Eq. (30) to the vector function defined by Eq. (32), the grad div operator is
(∇∇ · w)j =
1 2 [∂ wk + ∂k H · ∂j wk − ∂j H · ∂k wk + (∂jk2 H − ∂j H · ∂k H )wk ]. h2 jk
(33)
For the Laplacian, we use Eqs. (31) and (32) and get (∇ 2 w)j =
1 2 [∇wj − 2∂j H · ∂k wk + ∂k H · ∂k wj + ∂k H · ∂j wk + ∂kk Hwj − (∂j H ) · (∂k H ) · wk ]. h2
(34)
By substituting the vector operators of Eqs. (33) and (34) in Eq. (28), Naveir’s system of equation in {x} has the form: (Lu)j + ρω2 uj = fj , where (Lu)j =
1 {µ∇ 2 uj + (λ + µ)∂jk2 uk + µ∂k H · ∂k uj − (λ + 3µ)∂j H · ∂k uk h2 +2(λ + 2µ)∂k H · ∂j uk + K0 uj + Kjk · uk }.
2 H µ[∂11
2 H] ∂22
+ and Kjk = (2λ + Let K0 = operator L in the above equation are
3µ)∂jk2 H
(35)
− 2(λ + 2µ)(∂j H ) · (∂k H ). Then, the two components of
1 2 2 2 2 [(λ + 2µ)(∂11 u1 + ∂22 u1 ) − (λ + µ)(∂22 u1 − ∂12 u2 ) + (λ + 2µ)∂1 H · ∂1 u1 + µ∂2 H · ∂2 u1 h2 +2(λ + 2µ)∂2 H · ∂1 u2 − (λ + 3µ)∂1 H · ∂2 u2 + (K0 + K11 )u1 + K12 u2 ] + ρω2 u1 = f1 , 1 2 2 2 2 u2 + ∂22 u2 ) − (λ + µ)(∂11 u2 − ∂21 u1 ) + (λ + 2µ)∂2 H · ∂2 u2 + µ∂1 H · ∂1 u2 (Lu)2 ≡ 2 [(λ + 2µ)(∂11 h +2(λ + 2µ)∂1 H · ∂2 u1 − (λ + 3µ)∂2 H · ∂1 u1 + K21 u1 + (K0 + K22 )u2 ] + ρω2 u2 = f2 . (36) (Lu)1 ≡
As with the 3D case, the Kj k terms in the transformed operator correspond to spatial curvature, while the first order derivatives of the dependent variable indicate dispersion phenomena. Using matrix notation, the Navier’s equations are re-written as 1 1 2 2 2 2 [(λ + 2µ)(∂11 + ∂22 ) − (λ + µ)∂22 [(λ + µ)∂ + 2(λ + 2µ)∂ H · ∂ 2 1 12 2 l h h2 +(λ + 2µ)∂ H · ∂ + µ∂ H · ∂ + K + K ] + ρω2 −(λ + 3µ)∂1 H · ∂2 + K12 ] 1 1 2 2 0 11 1 2 2 2 1 [(λ + 2µ)(∂ + ∂ ) − (λ + µ)∂ 11 22 11 2 h2 2 [(λ + µ)∂21 + 2(λ + 2µ)∂1 H · ∂2 h +(λ + 2µ)∂2 H · ∂2 + µ∂1 H · ∂1 −(λ + 3µ)∂2 H · ∂1 + K21 ] +K0 + K22 ] + ρω2 u1 f1 × = . (37) u2 f2
196
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
Remark 2. It is possible to represent the transform Navier system in unified form for both 2D and 3D cases. Let In (n = 2, 3) be a conformal mapping in Rn ; then, the transform Navier operator Ln acting on vector-valued functions u(n ) has the general form (Ln u(n) )j =
1 (n) (n) (n) (n) {µ∇n2 uj + (λ + µ)∂jk2 uk + µ∂k Hn · ∂k uj + (n − 1)(λ + 2µ)∂k Hn · ∂j uk h2n (n)
(n)
2 −(λ + 3µ)∂j Hn · ∂k uk + µ∂kk Hn uj + [((n − 1)λ + (2n − 3)µ)∂jk2 Hn (n)
−(n − 1)(λ + 2µ)∂j Hn · ∂k Hn ]uk }.
(38)
In the above, Hn = ln|hn |, hn = |J n |1/n and J n is the Jacobian corresponding to conformal mapping In . The differential operators that appear in the Navier’s equations are as follows: (div gradn u(n) )j =
1 2 (n) (n) (n) {∂ u + (n − 1)∂k Hn · ∂j uk − ∂j Hn · ∂k uk + (n − 1)[∂jk2 Hn h2n jk k (n)
−(n − 1)∂j Hn · ∂k Hn ]uk }, (curl2n u(n) )j =
1 (n) (n) (n) (n) (n) (n) 2 {−∇n2 uj + ∂jk2 uk − ∂k Hn · ∂k uj + ∂j Hn · ∂k uk + ∂jk2 Hn uk − ∂kk Hn uj }, h2n
(∇n2 u(n) )j = (div gradn u(n) )j − (curl2n u(n) )j =
1 (n) (n) (n) (n) {∇ 2 u − 2∂j Hn · ∂k uk + (n − 1)∂k Hn · ∂j uk + ∂k Hn · ∂k uj h2n n j (n)
(n)
2 +[(n − 2)∂jk2 Hn − (n − 1)∂j Hn · ∂k Hn ]uk + ∂kk Hn · uj }.
(39)
5.1. Examples for 2D conformal mappings Let I be a conformal mapping in 2D defined by an analytical, complex function f (z) : B → C. This implies that f is one-to-one and regular in B, i.e., f (z) = 0 for z ∈ B. If we define f = x1 (X) + ix2 (X), then x = x(X), J = (∂x/∂X) and det J 2 (X0 ) = |f (z0 )|, where z0 = X10 + iX20 . Conformal mappings preserve the angle between any two curves that pass through point z0 , while the dilatation (expansion) is equal to |f (z0 )| in every direction [18]. The latter property means that as z → z0 locally, we have |f (z) − f (z0 )| = |z − z0 | · |f (z0 )|, which implies (x1 − x10 )2 + (x2 − x20 )2 = cases of conformal mappings:
(X1 − X10 )2 + (X2 − X20 )2 |J(X 0 )|2 . In what follows, we will present two basic
(a) Exponential mapping: Let mapping I be defined by function f (z) = ez = eX1 (cos X2 + i sin X2 ), which is conformal in B = {(X1 , X2 ) : X1 ∈ R1 , X2 ∈ (−π/2, π/2)} and I(B) = R2 \{(x1 , 0) : x1 > 0}, as shown X1 X1 in Fig. 2(a). The direct transformation is x1 = e cos X2 , x2 = e sin X2 , and the inverse transformation is X1 = ln x12 + x22 , X2 = arctan x2 /x1 . Also, the Jacobian matrix is eX1 cos X2 −eX1 sin X2 x1 −x2 J= = , x2 x1 eX1 sin X2 eX1 cos X2 its determinant is |J| = x12 + x22 , the metric tensor is 2 x1 + x22 0 T g = JJ = , 0 x12 + x22
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
197
Fig. 2. Domain of definition for conformal mapping in the complex plane; mappings defined by: (a) f (z) = exp(z); (b) f (z) = z2 .
and finally h = x12 + x22 and H = 0.5 ln(x12 + x22 ). If the radial distance in the transformed coordinate system is r = x12 + x22 and its derivatives with respect to the spatial coordinates are r,j = xj /r, then the coefficients multiplying the dependent variable uk and its derivatives ∂j uk in Eq. (37) are given below as 1 1 ∂jk2 H = 2 (δjk − 2r,j · r,k ), r,j , r r 1 Kjk = 2 [(2λ + 3µ)δjk − (6λ + 10µ)r,j · r,k ]. r
∂j H =
K0 = 0, (40)
(b) Quadratic mapping: In this case, let mapping I be defined by f (z) = z2 = X12 − X22 + i2X1 X2 . As shown in Fig. 2(b), I is conformal in B = {(X1 , X2 ) : X1 ∈ R1 , X2 > 0} and I(B) = R2 \{(x1 , 0) : x1 > 0}. The transformation is x1 = X12 − X22 , x2 = 2X1 X2 , the inverse transformation is −x + x 2 + x 2 1 x2 1 2 , X2 = X1 = 2 2(−x1 + x12 + x22 ) and the remaining properties concerning I are 2X1 −2X2 , |J| = 4 x12 + x22 , J= 2X1 2X2
h = 2(x12 + x22 )1/4 ,
and
H = 0.25 ln(x12 + x22 ) + ln 2. As before, the coefficients multiplying uk and ∂j uk in the transform system of equation (37) are given as 1 1 r,j , ∂jk2 H = 2 (δjk − 2r,j · r,k ), 2r 2r 1 Kjk = 2 [(2λ + 3µ)δjk − (6λ + 10µ)r,j · r,k ]. 2r ∂j H =
K0 = 0, (41)
198
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
6. Numerical examples In this section, we present numerical results for the anti-plane strain case under two conformal mappings, namely the linear and the inverse, which were outlined in Sections 4.1 and 4.2, respectively. The geological medium through which the wave signals propagate has the following background (or reference) material properties: λ0 = 5.0 × 109 N/m2 , cp0 = 2.0 km/s,
µ0 = 2.5 × 109 N/m2 ,
ρ0 = 2500 kg/m3 ,
cs0 = 1.0 km/s,
(42)
where cp0 and cs 0 are the pressure (P) and shear (S) wave velocities. We will further assume that the aforementioned waves travel at an intermediate frequency of ω = 2.0 Hz, which yields the following values for the P and S wave number and wavelength: kp0 =
ω = 6.28 km−1 , cp0
λp0 =
2π = 1.0 km, kp0
ks0 =
ω = 12.56 km−1 , cs0
λs0 =
2π = 0.5 km. ks0 (43)
Fig. 3. Wave amplitude |U3 | versus X1 , X2 for the following linear mapping cases: (a) homogeneous material: h0 = 1, ϕ = 0◦ ; (b) inhomogeneous material: h0 = 2, ϕ = 0◦ ; (c) inhomogeneous material: h0 = 1; ϕ = 45◦ .
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
199
Example 1. The linear conformal mapping, which is a superposition of rotation (with angle of rotation ϕ) and dilatation (with stretching factor h0 ), yields the following transformation from {X} to {x}: x1 = h0 (X1 cos ϕ − X2 sin ϕ),
x2 = h0 (X1 sin ϕ + X2 cos ϕ),
x3 = h0 X3 .
(44)
The solution for U3 (see Eq. (25)) in terms of the modified Bessel function of the second kind [19] for f3 = δ(x) at the origin of the coordinate system is
1 2 2 2 U3 = − K0 −iks0 h0 x1 + x2 , (45) 2π where ks0 h20 is the effective wave number. Fig. 3 plots |U3 | for an {X1 , X2 } grid of dimensions comparable to the S wave length (i.e., 0.5 × 0.5 km2 ) and at the x3 = X3 = 0 level. All plots here have been obtained by using the Mathematica package [20]. Specifically, Fig. 3(a) is for the homogeneous case, which results when h0 = 1, ϕ = 0◦ , while Fig. 3(b) and (c), respectively, plot the same information, but for the variable parameter cases that result when h0 = 2, ϕ = 0◦ and h0 = 1, ϕ = 45◦ . We observe a change in magnitude in |U3 | when h0 = 2 due to the ‘stretching’ of coordinates, while pure rotation (ϕ = 45◦ ) produces no noticeable effect.
Fig. 4. Wave amplitude |U3 | for the inverse mapping case versus X1 , X2 at the following X3 levels: (a) 0.01 km; (b) 0.25 km; (c) 0.50 km.
200
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
Example 2. For this case, the inverse mapping x = X/|X|2 yields the following values for the curvature terms Kˆ 13 and Kˆ 23 in Eq. (21) and for the wave number kˆ0 in Eq. (23): 2µ 1 2µ 1 r,1 , Kˆ 23 = r,2 , λ+µr λ+µr
2 ˆk0 = ρω · 1 1 − 4µ [(r,1 )2 + (r,2 )2 ] − 2 1 + 2λ + 3µ + 2(r,3 )2 . µ r4 λ+µ µ Kˆ 13 =
(46)
We note that the effective wave number in this coordinate system has units of length. If kˆ0 is to remain constant (equal to ks 0 ), then the density profile ρ(x) is equal to
µ 2 4µ 2λ + 3µ ρ(x) = k02 r 2 − r . (47) [(r,1 )2 + (r,2 )2 ] + 2 1 + + 2(r,3 )2 λ+µ µ ω2 We may assume for the purpose of this example that the shear modulus µ remains constant, although it is obvious that it is possible to satisfy Eq. (46) by varying both ρ and µ. The solution |U3 | for a point load at the origin has the same form as that given by Eq. (45), but with the modified wave number kˆ0 replacing the effective wave number and with the inverse transformation now defined as xj = Xj /|X|2 . Fig. 4(a)–(c) plot the amplitude of U3 for the same X3 , X2 spatial grid as before, and at three levels of X3 , namely X3 = 0.01, 0.25 and 0.50 km, respectively. We see a very unusual, almost flat profile in Fig. 4(a), while in Fig. 4(c) the profile approaches that of a homogeneous material. The second profile can be viewed as an intermediate one. Concurrently plotted in Fig. 5(a)–(c) are the corresponding density profiles, where it is observed that the effect of inhomogeneity is less pronounced the farther
Fig. 5. Variables density profiles for the inverse mapping case at the following X3 levels: (a) 0.01 km; (b) 0.25 km; (c) 0.50 km.
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
201
away one moves form the (X1 , X2 ) reference plane. Note that these densities are normalized with respect to their peak values. Specifically, in Fig. 5(a) the material becomes denser around the origin, while in Fig. 5(c) the density is almost uniform. As before, Fig. 5(b) plots an intermediate situation.
7. Conclusions In this work, we have shown the applicability of conformal mapping as a method for obtaining fundamental solutions to both 2D and 3D problems in elastodynamics that involve non-homogeneous elastic media. Conformal mapping is chosen here because the underlying Cauchy–Riemann conditions act as constraints, limiting on one hand the number of possible mappings but, on the other hand, yielding transformed vector operators in the new coordinate system which are simpler in form than they would otherwise be. In earlier work, we showed that for 1D cases the transformed operator obtained by conformal mapping was essentially a scaled version of the original one, so that it was possible to obtain a Green’s function by inspection. The inverse transformation, which completed the solution, was an algebraic one and thus simple to perform since it involved a change of coordinates. Also, the density profile that was recovered depended on the conformal mapping used and was physically plausible. As demonstrated here, all the above advantages of conformal mapping are retained for the vector wave (Navier) equations of elastodynamics, despite the fact that the ensuing coordinate transformations involved a change of the basis vectors. Furthermore, it is possible to extend the conformal mapping method the 3D case, despite the underlying complex variable formalism. Finally, some specific cases were examined in the context of anti-plane as well as of plane strain elastodynamics and the necessary conditions for obtaining closed-form solutions were identified.
Acknowledgements The authors wish to acknowledge the financial support provided by the NATO Science and Environment Affairs Division through grant No. EST.CLG 977774/2001. References [1] W.M. Ewing, W.S. Jardetsky, F. Press, Elastic Waves in Layered Media, McGraw-Hill, New York, 1957. [2] T.A. Cruse, F.A. Rizzo, A direct formulation and numerical solution of the general transient elastodynamic problem. I, J. Math. Anal. Appl. 22 (1968) 244–259. [3] G.D. Manolis, D.E. Beskos, Boundary Element Methods in Elastodynamics, Allen & Unwin, London, 1988. [4] J. Dominguez, Boundary Elements in Dynamics, Computational Mechanics Publications, Southampton, 1993. [5] W.C. Chew, Waves and Fields in Inhomogeneous Media, Van Nostrand Reinhold, New York, 1990. [6] K. Helbig (Ed.), Modeling the Earth for Oil Exploration, Pergamon/Elsevier, Oxford/Amsterdam, 1993. [7] A.H. Nayfeh, Wave Propagation in Layered Anisotropic Media, North-Holland, Amsterdam, 1995. [8] H. Sato, M.C. Fehler, Wave Propagation and Scattering in the Heterogeneous Earth, Springer, New York, 1997. [9] R.P. Shaw, G.D. Manolis, Conformal mapping solutions for 2D heterogeneous Helmholtz equations, in: M. Marchetti, C.A. Brebbia, H.H. Aliabadi (Eds.), Proceedings of the 19th BEM International Conference, Computational Mechanics Publications, Southampton, 1997, pp. 411–418. [10] G.D. Manolis, R.P. Shaw, Green’s function for the vector wave equation in a mildly heterogeneous medium, Wave Motion 24 (1996) 59–83. [11] J.D. Achenbach, Wave Propagation in Elastic Solids, North-Holland, Amsterdam, 1973. [12] J. Gerretsen, Lectures on Tensor Calculus and Differential Geometry, Noordhoff, Groningen, 1962. [13] P.M. Morse, H. Feshbach, Methods of Theoretical Physics, Vol. 1, McGraw-Hill, New York, 1953. [14] G.D. Manolis, R.P. Shaw, Fundamental solutions for variable density two-dimensional elastodynamic problems, Eng. Anal. Bound. Elements 24 (10) (2000) 739–750. [15] B. Dubrovin, S. Novikov, A. Fomenko, Modern Geometry: Methods and Applications, Nauka, Moscow, 1979 (in Russian). [16] V.S. Vladimirov, Equations of Mathematical Physics, Mir, Moscow, 1984.
202
G.D. Manolis et al. / Wave Motion 36 (2002) 185–202
[17] G.D. Manolis, R.P. Shaw, Wave motions in stochastic heterogeneous media: a Green’s function approach, Eng. Anal. Bound. Elements 15 (1995) 225–234. [18] Y. Sidorov, M. Fedoriuk, M. Shabunin, Lectures on the Theory of Complex Functions, Nauka, Moscow, 1982 (in Russian). [19] I.S. Gradsteyn, I.M. Ryzhik, Table of Integrals, Series and Products, Academic Press, New York, 1980. [20] S. Wolfram, A System for Doing Mathematics by Computer, Version 3.0, Cambridge University Press, Cambridge, 1996.