A Mathematical Analysis of the Strip-Element Method for the Computation of Dispersion Curves of Guided Waves in Anisotropic Layered Media

Abstract

Dispersion curves of elastic guided waves in plates can be efficiently computed by the Strip-Element Method. This method is based on a finite-element discretization in the thickness direction of the plate and leads to an eigenvalue problem relating frequencies to wavenumbers of the wave modes. In this paper we present a rigorous mathematical background of the Strip-Element Method for anisotropic media including a thorough analysis of the corresponding infinite-dimensional eigenvalue problem as well as a proof of the existence of eigenvalues.

1. Introduction

In recent years there has been considerable interest in the study of the behaviour of elastic guided waves in plates due to their potential use in Nondestructive Evaluation (NDE) and Structural Health Monitoring (SHM); see, for example, the comprehensive books of Giurgiutiu [1] or Rose [2]. Already in isotropic plates Lamb waves are dispersive and the dispersion relations expressed by the Rayleigh-Lamb equations must be computed numerically; see, the study by Achenbach in [3]. Elastic wave propagation in layered and anisotropic media is an even more complex problem, and efficient numerical methods are required to obtain dispersion curves. Among those methods are the Transfer Matrix Method and the Global Matrix Method and we refer to the study by Lowe in [4] for an overview. One of the most efficient and flexible methods is based on a finite-element discretization in the thickness direction of the plate and leads to a generalized eigenvalue problem relating frequencies to wavenumbers of the Lamb wave modes. This method is known, amongst others, as the strip-element method (SEM) or layer-element method or semianalytic finite-element-method; see, for example, the early works of Dong and Nelson [5] and Aalami [6], or Kausel [7], Galán and Abascal [8], the excellent book of Liu and Xi [9], and the many references therein. Also the works of Gavrić [10], Bartoli et al. [11], Marzani et al. [12], and Treyssède [13] could be of interest.

So far, there seems to be no strict mathematical analysis of the underlying infinite-dimensional eigenvalue problem in the anisotropic case. For the isotropic case we recommend reading the paper by Bouhennache in [14], who also uses a more abstract setting. Since the SEM has been successfully applied in practice for several years, we think that it is worthwhile to start such an analysis and give here a mathematical proof of the existence of eigenvalues.

In the next section we recall some basic facts about generalized eigenvalue problems in Hilbert spaces. The differential equation governing the wave propagation in laminated plates is formulated in Section 3. In the main Section 4 we analyse the weak form of those equations for the elementary Lamb wave modes, show that weak and strong solutions coincide and are layerwise 𝒞^∞, and prove existence of weak solutions and eigenvalues of the related eigenvalue problem. Finally we present some numerical results illustrating the increasingly direction-dependent behaviour of transversely isotropic material with an increasing degree of anisotropy.

2. Generalized Eigenvalue Problems in Hilbert Spaces

Although the results of this section are wellknown, they might not always be explicitly found in the given form and therefore we prove some of them for the convenience of the reader. For further reading we suggest the books of Lax [15] and Conway [16].

Let ℋ,ℋ₁,ℋ₂ be real or complex Hilbert spaces with scalar products ⟨⋅∣⋅⟩. By ℒ(ℋ₁,ℋ₂) we denote the space of continuous linear operators from ℋ₁ to ℋ₂ and set ℒ(ℋ)∶=ℒ(ℋ,ℋ). For 𝑇∈ℒ(ℋ₁,ℋ₂) we denote by 𝑇^∗∈ℒ(ℋ₂,ℋ₁) the adjoint operator defined by ⟨𝑇𝑣∣𝑢⟩_ℋ2=⟨𝑣∣𝑇^∗𝑢⟩_ℋ1,∀𝑣∈ℋ₁,𝑢∈ℋ₂. (2.1) Range and nullspace are denoted by ℛ(𝑇) and 𝒩(𝑇), respectively. We say that 𝑇∈ℒ(ℋ) is self-adjoint if 𝑇^∗=𝑇, and positive if ⟨𝑢∣𝑇𝑢⟩_ℋ≥0,∀𝑢∈ℋ. (2.2)

Lemma 2.1. Let 𝑇∈ℒ(ℋ) be self-adjoint. If 𝑇 is injective and has closed range, then it is bijective and hence continuously invertible.

Proof. For 𝒳⊂ℋ, let 𝒳^⟂={𝑥∈ℋ∣⟨𝑥∣𝑦⟩=0,forall𝑦∈ℋ} be the orthogonal complement of 𝒳. The assertion follows from ℛ(𝑇)=ℛ(𝑇)=𝒩(𝑇^∗)^⟂=𝒩(𝑇)^⟂={0}^⟂=ℋ, (2.3) which means that 𝑇 is also surjective.

Lemma 2.2. An operator 𝑇∈ℒ(ℋ₁,ℋ₂) is injective and has closed range if and only if there exists some constant 𝑐>0 such that ‖𝑇𝑢‖_ℋ2≥𝑐‖𝑢‖_ℋ1,∀𝑢∈ℋ₁. (2.4)

Proof. See, for example, the study by Schröder in [17].

We obtain the following.

Corollary 2.3. An operator 𝑇∈ℒ(ℋ) is injective and has closed range if there exists some constant 𝑐>0 such that ⟨𝑢∣𝑇𝑢⟩_ℋ≥𝑐‖𝑢‖²_ℋ,∀𝑢∈ℋ. (2.5)

Proof. Together with the Cauchy-Schwartz inequality, we obtain from (2.5) ‖𝑢‖_ℋ‖𝑇𝑢‖_ℋ≥⟨𝑢∣𝑇𝑢⟩_ℋ≥𝑐‖𝑢‖²_ℋ,∀𝑢∈ℋ. (2.6) With Lemma 2.2, we can conclude the assertion.

For 𝑀,𝑆∈ℒ(ℋ), consider the generalized eigenvalue problem 𝑀𝑢=𝜆𝑆𝑢. (2.7) That is, “for which 𝜆∈ℝ,ℂ do nontrivial solutions 𝑢∈ℋ to (2.7) exist?” In case 𝑆 is bijective, (2.7) is equivalent to the standard eigenvalue problem 𝑆⁻¹𝑀𝑢=𝜆𝑢, (2.8) and we readily obtain the following proposition which we state for the case dimℋ=∞. It applies accordingly to the finite-dimensional case with finite sequences 𝜆_𝑛, 𝑢_𝑛.

Proposition 2.4. Let 𝑀,𝑆∈ℒ(ℋ) be self-adjoint and positive. Assume that 𝑀 is compact and injective and that 𝑆 is bijective. Then there exists a decreasing sequence 0<𝜆_𝑛→0 (counted with multiplicity) and a sequence 𝑢_𝑛∈ℋ such that 𝑀𝑢_𝑛=𝜆_𝑛𝑆𝑢_𝑛,1𝜆_𝑛⟨𝑢_𝑛∣𝑀𝑢_𝑚⟩_ℋ=𝛿_𝑛,𝑚=⟨𝑢_𝑛∣𝑆𝑢_𝑚⟩_ℋ, (2.9)(2.10) and each 𝑢∈ℋ can be uniquely represented by an ℋ-convergent series 𝑢=^∞_𝑛=1𝛼_𝑛𝑢_𝑛, (2.11) with an ℓ₂-sequence 𝛼_𝑛=⟨𝑢∣𝑆𝑢_𝑛⟩_ℋ. Especially, ⟨𝑣∣𝑢⟩_𝑆∶=⟨𝑣∣𝑆𝑢⟩_ℋ defines a scalar product on ℋ which induces an equivalent norm on ℋ, and 𝑢_𝑛 is an orthonormal basis with respect to this scalar product.

Proof. Since 𝑆 is self-adjoint, bijective, and positive, ⟨𝑣∣𝑢⟩_𝑆∶=⟨𝑣∣𝑆𝑢⟩_ℋ indeed defines a scalar product on ℋ which induces an equivalent norm on ℋ. With respect to this scalar product, the injective and compact operator 𝑀∶=𝑆⁻¹𝑀 is self-adjoint and positive. Therefore (2.9) is equivalent to the standard eigenvalue problem 𝑀𝑢_𝑛=𝜆_𝑛𝑢_𝑛 for the injective, compact, self-adjoint, and positive operator 𝑀∈ℒ{ℋ,⟨⋅∣⋅⟩_𝑆}.

3. Differential Equations in Matrix Notation

We recall some relations of elasticity theory in matrix notation which is especially suited for the finite-element formulation [9]. For a profound study of the mathematical theory of elasticity in tensor notation, we refer to the book of Marsden and Hughes in [18]. (i)Differential operator matrix 𝐿 is given as, 𝐿=𝐿_𝑥𝜕_𝑥+𝐿_𝑦𝜕_𝑦+𝐿_𝑧𝜕_𝑧, (3.1) with the constant matrices 𝐿_𝑥=100000000000001010,𝐿_𝑦=000010000001000100,𝐿_𝑧=000000001010100000. (3.2)(ii)Displacement vector 𝑢 is given as 𝑢=𝑢(𝑥,𝑦,𝑧,𝑡)=(𝑢_𝑥,𝑢_𝑦,𝑢_𝑧)^𝑇. (3.3)(iii)Strain vector 𝜖 is given as 𝜖=(𝜖_𝑥,𝜖_𝑦,𝜖_𝑧,𝛾_𝑦𝑧,𝛾_𝑥𝑧,𝛾_𝑥𝑦)^𝑇. (3.4)(iv)Strain-displacement relation is given as 𝜖=𝐿𝑢. (3.5)(v)Stress vector 𝜎 is given as 𝜎=(𝜎_𝑥,𝜎_𝑦,𝜎_𝑧,𝜏_𝑦𝑧,𝜏_𝑥𝑧,𝜏_𝑥𝑦)^𝑇. (3.6)(vi)And the generalized Hooke Law is given as 𝜎=𝐶𝜖, (3.7) with 𝐶=(𝐶_𝑖𝑗)_{𝑖,𝑗=1,…,6} being the matrix of material constants. In the following, 𝐶 is supposed to be real, symmetric, and positive definite (spd). For an isotropic material we have, for example, 𝐶_isotropic=𝜆+2𝜇𝜆𝜆000𝜆𝜆+2𝜇𝜆000𝜆𝜆𝜆+2𝜇000000𝜇000000𝜇000000𝜇, (3.8) with Lamé's constants 𝜆, 𝜇.

Now consider a laminated plate of thickness 2𝐻 in direction of the 𝑧-axis (middle 𝑧=0, top 𝑧=+𝐻, bottom 𝑧=−𝐻) and infinite in the 𝑥-𝑦 plane; see Figure 1. The plate consists of 𝑁 layers. Layer 𝑙=1,…,𝑁 has thickness ℎ_𝑙=𝐻_𝑙−𝐻_𝑙−1 with −𝐻=𝐻₀<𝐻₁<⋯<𝐻_𝑁=+𝐻 and is supposed to consist of homogenous, anisotropic, elastic material with density 𝜌_𝑙 and material constants 𝐶_𝑙.

Figure 1 

Laminate and coordinate systems.

Let 𝑢=(𝑢_𝑥,𝑢_𝑦,𝑢_𝑧)^𝑇=𝑢(𝑥,𝑦,𝑧,𝑡) be the displacement vector of a wave travelling in the plate in the absence of external forces. In each layer the elastic wave equation in matrix notation is 𝜌_𝑙𝑢̈=𝐿^𝑇𝐶_𝑙𝐿𝑢,∀𝑥,𝑦,𝑡andall𝑧∈(𝐻_𝑙−1,𝐻_𝑙). (3.9) The traction-free boundary conditions on the top and bottom surfaces are 𝐿^𝑇_𝑧𝐶_𝑁𝐿𝑢|_𝑧=+𝐻=𝐿^𝑇_𝑧𝐶₁𝐿𝑢|_𝑧=−𝐻=0. (3.10) Besides continuity of the displacement vector, the following interface conditions concerning continuity of the stresses are supposed to hold: 𝐿^𝑇_𝑧𝐶_𝑙𝐿𝑢|_{𝑧=𝐻𝑙}=𝐿^𝑇_𝑧𝐶_𝑙+1𝐿𝑢|_{𝑧=𝐻𝑙},∀𝑙=1,…,𝑁−1. (3.11) As ansatz for a wave mode, we take a plane harmonic wave travelling in the 𝑥-𝑦 plane 𝑢(𝑥,𝑦,𝑧,𝑡)=𝑢(𝑧)𝑒^{𝑖(𝑘𝑥𝑥+𝑘𝑦𝑦±𝜔𝑡)}, (3.12) with 𝑧-dependent amplitude vector 𝑢(𝑧)=(𝑢_𝑥(𝑧),𝑢_𝑦(𝑧),𝑢_𝑧(𝑧))^𝑇, real circular frequency 𝜔, and real wave vector 𝑘=(𝑘_𝑥,𝑘_𝑦)^𝑇. For such a wave mode, the wave equation (3.9), boundary (3.10) and interface conditions (3.11) reduce to −𝜌_𝑙𝜔²𝑢=𝐿^𝑇_𝑘𝐶_𝑙𝐿_𝑘𝑢,∀𝑧∈(𝐻_𝑙−1,𝐻_𝑙𝐿),𝑙=1,…,𝑁,_𝑇𝑧𝐶_𝑁𝐿_𝑘𝑢|_𝑧=+𝐻=𝐿^𝑇_𝑧𝐶₁𝐿_𝑘𝑢|_𝑧=−𝐻𝐿=0,_𝑇𝑧𝐶_𝑙𝐿_𝑘𝑢|_{𝑧=𝐻𝑙}=𝐿^𝑇_𝑧𝐶_𝑙+1𝐿_𝑘𝑢|_{𝑧=𝐻𝑙},∀𝑙=1,…,𝑁−1, (3.13)(3.14)(3.15) with the differential operator matrix 𝐿_𝑘 given as 𝐿_𝑘=𝑖𝑘_𝑥𝐿_𝑥+𝑖𝑘_𝑦𝐿_𝑦+𝐿_𝑧𝜕_𝑧. (3.16) Obviously for complex conjugation we have 𝐿_𝑘=−𝑖𝑘_𝑥𝐿_𝑥−𝑖𝑘_𝑦𝐿_𝑦+𝐿_𝑧𝜕_𝑧. (3.17) The question is “for which combinations of circular frequencies 𝜔 and wave vectors 𝑘 do nontrivial solutions 𝑢 of (3.13), (3.14), and (3.15) exist?” The answer leads to the dispersion relations 𝜔(𝑘) and is given in the next section. For better readability we will write 𝜔 instead of 𝜔(𝑘) in the following.

4. Weak Form of the Reduced Wave Equation and Related Eigenvalue Problem

We define the piecewise constant functions of density 𝜌 and matrix of material constants 𝐶 as 𝜌(𝑧)=𝜌_𝑙,𝐶(𝑧)=𝐶_𝑙,∀𝑧∈(𝐻_𝑙−1,𝐻_𝑙),𝑙=1,…,𝑁. (4.1) For suitable virtual displacements 𝑣=𝑣(𝑧), which we define later, the weak form of (3.13), (3.14), and (3.15) can then be written as 𝜔^2+𝐻_−𝐻𝑣^𝑇𝜌𝑢𝑑𝑧=^+𝐻_−𝐻(𝐿_𝑘𝑣)^𝑇𝐶(𝐿_𝑘𝑢)𝑑𝑧. (4.2)

The next proposition shows that strong and weak solutions coincide for smooth enough functions.

Proposition 4.1. Suppose that 𝑢=𝑢(𝑧) is continuous and layerwise 𝒞²; that is, 𝑢∈𝒞([−𝐻,+𝐻]) and 𝑢|_{(𝐻𝑙−1,𝐻𝑙)}∈𝒞²([𝐻_𝑙−1,𝐻_𝑙]), 𝑙=1,…,𝑁. Then 𝑢 fulfills (4.2) for all 𝑣∈𝒞¹([−𝐻,+𝐻]) if and only if 𝑢 is a solution of (3.13), (3.14), and (3.15).

Proof. Note that we shortly write 𝑢|_{(𝐻𝑙−1,𝐻𝑙)}∈𝒞²([𝐻_𝑙−1,𝐻_𝑙]) meaning that there is a function 𝑢∈𝒞²([𝐻_𝑙−1,𝐻_𝑙]) such that 𝑢|_{(𝐻𝑙−1,𝐻𝑙)}=𝑢|_{(𝐻𝑙−1,𝐻𝑙)}. By the smoothness assumptions on 𝑢, relation (3.17), and layerwise partial integration on the right-hand side of (4.2), we see that for all 𝑣∈𝒞¹([−𝐻,+𝐻]) (4.2) is equivalent to 𝜔^2+𝐻_−𝐻𝑣^𝑇𝜌𝑢𝑑𝑧=−^+𝐻_−𝐻𝑣^𝑇𝐿^𝑇_𝑘𝐶𝐿_𝑘𝑢𝑑𝑧+^𝑁_𝑙=1[𝑣^𝑇𝐿^𝑇_𝑧𝐶_𝑙𝐿_𝑘𝑢]^𝐻𝑙_𝐻𝑙−1. (4.3) Let 𝑢 fulfill this equation for all 𝑣∈𝒞¹([−𝐻,+𝐻]). At first we fix a layer 𝑙=1,…,𝑁 and choose arbitrary functions 𝑣 that have compact support in the interior (𝐻_𝑙−1,𝐻_𝑙) of the layer. For such functions, (4.3) reduces to 𝜔^2𝐻𝑙_𝐻𝑙−1𝑣^𝑇𝜌_𝑙𝑢𝑑𝑧=−^𝐻𝑙_𝐻𝑙−1𝑣^𝑇𝐿^𝑇_𝑘𝐶_𝑙𝐿_𝑘𝑢𝑑𝑧 (4.4) from which we infer that 𝑢 fulfills (3.13). Repeating this for all layers, we conclude that the integrals on the left- and right-hand side of (4.3) coincide and hence (4.3) reduces to 0=^𝑁_𝑙=1[𝑣^𝑇𝐿^𝑇_𝑧𝐶_𝑙𝐿_𝑘𝑢]^𝐻𝑙_𝐻𝑙−1. (4.5) Now by successively choosing functions 𝑣 that equal 1 in a vicinity of one of the points −𝐻=𝐻₀,𝐻₁,…,𝐻_𝑁=+𝐻 and that are 0 everywhere else, we see that 𝑢 also fulfills (3.14) and (3.15).
The converse assertion is obvious.

To prove existence of nontrivial weak solutions, we will transform (4.2) into a generalized eigenvalue problem in a suitable Hilbert space. Let 𝐻¹=𝐻¹(−𝐻,+𝐻) be the Sobolev space of all complex-valued square integrable functions on (−𝐻,+𝐻) whose distributional first derivative can also be identified with a square integrable function; see, the study by Adams in [19]. For simplicity in the following we use the same symbol 𝐻¹ also for (𝐻¹)³=𝐻¹×𝐻¹×𝐻¹,(𝐻¹)⁶,…, and likewise for other spaces like the spaces of square integrable functions 𝐿²=𝐿²(−𝐻,+𝐻) and continuously differentiable functions 𝒞¹=𝒞¹([−𝐻,+𝐻]),…, since it will become clear from the context how many components the vector-valued functions have. Endowed with the scalar product ⟨𝑣∣𝑢⟩_𝐻¹=⟨𝑣∣𝑢⟩_𝐿²+⟨𝜕_𝑧𝑣∣𝜕_𝑧𝑢⟩_𝐿2=^+𝐻_−𝐻𝑣(𝑧)^𝑇𝑢(𝑧)𝑑𝑧+^+𝐻_−𝐻𝜕_𝑧𝑣(𝑧)^𝑇𝜕_𝑧𝑢(𝑧)𝑑𝑧, (4.6) the space 𝐻¹ is a complex Hilbert space. The space 𝒞¹ is continuously embedded and dense in 𝐻¹ and the inclusions 𝐻¹↪𝒞 and 𝐽∶𝐻¹↪𝐿² are compact. See, for example, the study by Maz'ja in [20]. Multiplications with the piecewise constant function 𝜌>0 and the piecewise constant spd-matrix 𝐶, respectively, define bijective, positive, self-adjoint, continuous linear operators on the respective 𝐿²-spaces. Furthermore, the differential operator matrix 𝐿_𝑘 (3.16) defines a continuous linear operator 𝐿_𝑘∶𝐻¹→𝐿². For functions 𝑢,𝑣∈𝐻¹, which we assume in the following, the left-hand side of (4.2) can then be written as 𝜔^2+𝐻_−𝐻𝑣^𝑇𝜌𝑢𝑑𝑧=𝜔²⟨𝐽𝑣∣𝜌𝐽𝑢⟩_𝐿2=𝜔²⟨𝑣∣𝐽^∗𝜌𝐽𝑢⟩_𝐻1=𝜔²⟨𝑣∣𝑀𝑢⟩_𝐻¹ (4.7) with the compact, injective, positive, and self-adjoint operator 𝑀∶=𝐽^∗𝜌𝐽, and the right-hand side can be written as ^+𝐻_−𝐻(𝐿_𝑘𝑣)^𝑇𝐶(𝐿_𝑘𝑢)𝑑𝑧=⟨𝐿_𝑘𝑣∣𝐶𝐿_𝑘𝑢⟩_𝐿2=⟨𝑣∣𝐿^∗_𝑘𝐶𝐿_𝑘𝑢⟩_𝐻1=⟨𝑣∣𝑆_𝑘𝑢⟩_𝐻1 (4.8) with the positive and self-adjoint operator 𝑆_𝑘∶=𝐿^∗_𝑘𝐶𝐿_𝑘. Thus, a function 𝑢∈𝐻¹ fulfills (4.2) for all 𝑣∈𝐻¹ if and only if 𝜔²⟨𝑣∣𝑀𝑢⟩_𝐻¹=⟨𝑣∣𝑆_𝑘𝑢⟩_𝐻1,∀𝑣∈𝐻¹, (4.9) which is equivalent to the generalized eigenvalue problem 𝜔²𝑀𝑢=𝑆_𝑘𝑢⇔(𝛾+𝜔²)𝑀𝑢=(𝛾𝑀+𝑆_𝑘)𝑢⇔𝜆𝑀𝑢=𝑆_𝛾𝑢 (4.10) with 𝜆∶=𝛾+𝜔² and the positive and self-adjoint operator 𝑆_𝛾∶=𝛾𝑀+𝑆_𝑘 for an arbitrary 𝛾>0. We will need this disturbance with 𝛾𝑀 to prove that the operator 𝑆_𝛾 is also bijective and hence can apply Proposition 2.4 to guarantee nontrivial solutions. By Lemma 2.1 and Corollary 2.3, bijectivity of 𝑆_𝛾 follows by showing that there exists some constant 𝑐>0 such that ⟨𝑢∣𝑆_𝛾𝑢⟩_𝐻1≥𝑐‖𝑢‖²_𝐻1,∀𝑢∈𝐻¹. (4.11) To show this, we write ⟨𝑢∣𝑆_𝛾𝑢⟩_𝐻1=𝛾⟨𝑢∣𝑀𝑢⟩_𝐻¹+⟨𝑢∣𝑆_𝑘𝑢⟩_𝐻1=^+𝐻_−𝐻𝛾𝑢^𝑇𝜌𝑢+(𝐿_𝑘𝑢)^𝑇𝐶(𝐿_𝑘𝑢)𝑑𝑧. (4.12) Since layerwise 𝐶 is an spd-matrix and 𝜌>0, there is a constant 𝑐₁>0 such that pointwise a.e. 𝛾𝑢^𝑇𝜌𝑢+(𝐿_𝑘𝑢)^𝑇𝐶(𝐿_𝑘𝑢)≥𝑐₁(𝛾𝑢^𝑇𝑢+(𝐿_𝑘𝑢)^𝑇(𝐿_𝑘𝑢)), (4.13) and thus we have ⟨𝑢∣𝑆_𝛾𝑢⟩_𝐻1≥𝑐₁^+𝐻_−𝐻𝛾𝑢^𝑇𝑢+(𝐿_𝑘𝑢)^𝑇(𝐿_𝑘𝑢)𝑑𝑧=𝑐₁^+𝐻_−𝐻𝑢^𝑇𝑄_𝑘𝑢𝑑𝑧, (4.14) with 𝑢=(𝑢_𝑥,𝑢_𝑦,𝑢_𝑧,𝜕_𝑧𝑢_𝑥,𝜕_𝑧𝑢_𝑦,𝜕_𝑧𝑢_𝑧)^𝑇 (4.15) and the Hermitian matrix 𝑄_𝑘∶=|𝑘|²+𝛾𝑘_𝑥𝑘_𝑦𝑘0000_𝑥𝑘_𝑦|𝑘|²+𝛾000000|𝑘|²+𝛾−𝑖𝑘_𝑥−𝑖𝑘_𝑦000𝑖𝑘_𝑥10000𝑖𝑘_𝑦010000001, (4.16) where |𝑘|²=𝑘²_𝑥+𝑘²_𝑦≥0. For all real 𝑘=(𝑘_𝑥,𝑘_𝑦), the matrix 𝑄_𝑘 is positive definite since for 𝛾>0 all eigenvalues 𝜆_1,2=1,𝜆_3,4=|𝑘|²+𝛾±𝑘_𝑥𝑘_𝑦,𝜆_5,6=12(1+|𝑘|²+𝛾±(1+|𝑘|²+𝛾)²−4𝛾) (4.17) are strictly positive. Hence there is a constant 𝑐₂>0 such that pointwise a.e. 𝑢^𝑇𝑄_𝑘𝑢≥𝑐₂𝑢^𝑇𝑢, (4.18) and by the definition of the scalar product (4.6) we arrive at ⟨𝑢∣𝑆_𝛾𝑢⟩_𝐻1≥𝑐₁𝑐₂^+𝐻_−𝐻𝑢^𝑇𝑢𝑑𝑧=𝑐₁𝑐₂‖𝑢‖²_𝐻1. (4.19) Now we can apply Proposition 2.4 to the pair 𝑀, 𝑆_𝛾 and get the following.

Proposition 4.2. There exists an increasing sequence 0≤𝜔²_𝑛→∞ and a sequence 𝑢_𝑛∈𝐻¹ such that 𝜔²_𝑛𝑀𝑢_𝑛=𝑆_𝑘𝑢_𝑛,(𝛾+𝜔²_𝑛)⟨𝑢_𝑛∣𝑀𝑢_𝑚⟩_𝐻1=𝛿_𝑛,𝑚=𝛾+𝜔²_𝑛𝜔²_𝑛⟨𝑢_𝑛∣𝑆_𝑘𝑢_𝑚⟩_𝐻1, (4.20)(4.21) and each 𝑢∈𝐻¹ can be uniquely represented by an 𝐻¹-convergent series 𝑢=^∞_𝑛=0𝛼_𝑛𝑢_𝑛, (4.22) with an ℓ₂-sequence 𝛼_𝑛=((𝛾+𝜔²_𝑛)/𝜔²_𝑛)⟨𝑢∣𝑆_𝑘𝑢_𝑛⟩_𝐻1.

Proof. Let the assertions of Proposition 2.4 appropriately hold for the pair 𝑀, 𝑆_𝛾. At first we observe that by (2.10) we have 𝜆₁=⟨𝑢₁∣𝑆_𝛾𝑢₁⟩_𝐻1⟨𝑢₁∣𝑀𝑢₁⟩_𝐻1⟨𝑢=𝛾+¹∣𝑆_𝑘𝑢₁⟩_𝐻1⟨𝑢₁∣𝑀𝑢₁⟩_𝐻1≥𝛾, (4.23) and hence 𝜔²_𝑛∶=𝜆_𝑛−𝛾≥0 and 𝜔²_𝑛→∞ increasingly. By (2.10), on one hand, we then trivially have (𝛾+𝜔²_𝑛)⟨𝑢_𝑛∣𝑀𝑢_𝑚⟩_𝐻1=𝜆_𝑛⟨𝑢_𝑛∣𝑀𝑢_𝑚⟩_𝐻1=𝛿_𝑛,𝑚, (4.24) and, on the other hand, we further conclude that 𝛿_𝑛,𝑚=⟨𝑢_𝑛∣𝑆_𝛾𝑢_𝑚⟩_𝐻1=𝛾⟨𝑢_𝑛∣𝑀𝑢_𝑚⟩_𝐻1+⟨𝑢_𝑛∣𝑆_𝑘𝑢_𝑚⟩_𝐻1=𝛾𝛾+𝜔²_𝑛𝛿_𝑛,𝑚+⟨𝑢_𝑛∣𝑆_𝑘𝑢_𝑚⟩_𝐻1. (4.25) Finally each 𝑢∈𝐻¹ can be uniquely represented by an 𝐻¹-convergent series 𝑢=^∞_𝑛=0𝛼_𝑛𝑢_𝑛, (4.26) with an ℓ₂-sequence 𝛼_𝑛=⟨𝑢∣𝑆_𝛾𝑢_𝑛⟩_𝐻1. Evaluating the scalar product of both sides of the above equation with 𝑆_𝑘𝑢_𝑚, we get ⟨𝑢∣𝑆_𝑘𝑢_𝑚⟩_𝐻1=^∞_𝑛=0𝛼_𝑛⟨𝑢_𝑛∣𝑆_𝑘𝑢_𝑚⟩_𝐻1=𝛼_𝑚𝜔²_𝑚𝛾+𝜔²_𝑚. (4.27)

The next proposition together with Proposition 4.1 finally shows that all these weak 𝐻¹-solutions are indeed strong solutions.

Proposition 4.3. If 𝑢∈𝐻¹ fulfills (4.2) for all 𝑣∈𝐻¹, then 𝑢 is continuous and layerwise 𝒞^∞.

Proof. We inductively show that the distributional derivatives of 𝑢 can be identified with smooth functions. Expanding the operator 𝐿 according to its definition and rearranging (4.2) give ^+𝐻_−𝐻𝜕_𝑧𝑣^𝑇𝐿^𝑇_𝑧𝐶𝐿_𝑧𝜕_𝑧𝑢𝑑𝑧=^+𝐻_−𝐻𝑣^𝑇(𝜔²𝜌+(𝑖𝑘_𝑥𝐿_𝑥+𝑖𝑘_𝑦𝐿_𝑦)^𝑇𝐶𝐿_𝑘)−𝑢𝑑𝑧_+𝐻−𝐻𝜕_𝑧𝑣^𝑇𝐿^𝑇_𝑧𝐶(𝑖𝑘_𝑥𝐿_𝑥+𝑖𝑘_𝑦𝐿_𝑦)𝑢𝑑𝑧. (4.28) We fix a layer 𝑙 and take an arbitrary 𝒞^∞-function 𝑣 with compact support in (𝐻_𝑙−1,𝐻_𝑙). As 𝐶_𝑙 is symmetric, for those 𝑣 the left-hand side can be written as ^𝐻𝑙_𝐻𝑙−1𝜕_𝑧(𝐿^𝑇_𝑧𝐶_𝑙𝐿_𝑧)𝑣^𝑇𝜕_𝑧𝑢𝑑𝑧. (4.29) We integrate by parts on the right-hand side to get ^𝐻𝑙_𝐻𝑙−1𝜕_𝑧(𝐿^𝑇_𝑧𝐶_𝑙𝐿_𝑧)𝑣^𝑇𝜕_𝑧𝑢𝑑𝑧=^𝐻𝑙_𝐻𝑙−1𝑣^𝑇(𝜔²𝜌+(𝑖𝑘_𝑥𝐿_𝑥+𝑖𝑘_𝑦𝐿_𝑦)^𝑇𝐶_𝑙𝐿_𝑘)+𝑢𝑑𝑧_{𝐻𝑙𝐻𝑙−1}𝑣^𝑇𝐿^𝑇_𝑧𝐶_𝑙(𝑖𝑘_𝑥𝐿_𝑥+𝑖𝑘_𝑦𝐿_𝑦)𝜕_𝑧𝑢𝑑𝑧. (4.30) The operator 𝐿^𝑇_𝑧𝐶_𝑙𝐿_𝑧 is invertible since 𝐶_𝑙 is positive definite and 𝐿_𝑧 is injective; hence the previous equation for all 𝑣 is equivalent to ^𝐻𝑙_𝐻𝑙−1𝜕_𝑧𝑣^𝑇𝜕_𝑧𝑢𝑑𝑧=−^𝐻𝑙_𝐻𝑙−1𝑣^𝑇𝐴_𝑙𝑢𝑑𝑧, (4.31) where 𝐴_𝑙𝑢∶=−(𝐿^𝑇_𝑧𝐶_𝑙𝐿_𝑧)⁻¹(𝜔²𝜌+(𝑖𝑘_𝑥𝐿_𝑥+𝑖𝑘_𝑦𝐿_𝑦)^𝑇𝐶_𝑙𝐿_𝑘+𝐿^𝑇_𝑧𝐶_𝑙(𝑖𝑘_𝑥𝐿_𝑥+𝑖𝑘_𝑦𝐿_𝑦)𝜕_𝑧)𝑢. (4.32) Equation (4.31) states that the second distributional derivative 𝜕²_𝑧𝑢|_{(𝐻𝑙−1,𝐻𝑙)} of 𝑢 restricted to (𝐻_𝑙−1,𝐻_𝑙) can be identified with 𝐴_𝑙𝑢|_{(𝐻𝑙−1,𝐻𝑙)}. Since 𝑢∈𝐻¹, we have 𝐴_𝑙𝑢|_{(𝐻𝑙−1,𝐻𝑙)}∈𝐿²(𝐻_𝑙−1,𝐻_𝑙). Consequently, we may assume that 𝑢|_{(𝐻𝑙−1,𝐻𝑙)}∈𝐻²(𝐻_𝑙−1,𝐻_𝑙) and hence 𝑢|_{(𝐻𝑙−1,𝐻𝑙)}∈𝒞¹([𝐻_𝑙−1,𝐻_𝑙]). From this we in turn infer that 𝐴_𝑙𝑢|_{(𝐻𝑙−1,𝐻𝑙)}∈𝒞([𝐻_𝑙−1,𝐻_𝑙]), and (4.31) ensures that 𝜕²_𝑧𝑢|_{(𝐻𝑙−1,𝐻𝑙)}∈𝒞([𝐻_𝑙−1,𝐻_𝑙]); that is, 𝑢|_{(𝐻𝑙−1,𝐻𝑙)}∈𝒞²([𝐻_𝑙−1,𝐻_𝑙]). Repeating this argument proves the assertion.

We conclude this section by deriving a mode decomposition of general waves travelling in the plate. For a function 𝑢=𝑢(𝑥,𝑦), we denote by 𝑢=𝑢(𝑘_𝑥,𝑘_𝑦) its Fourier transform in the wavenumber domain 𝑢(𝑘_𝑥,𝑘_𝑦1)=(2𝜋)^2ℝ2𝑢(𝑥,𝑦)𝑒^{−𝑖(𝑘𝑥𝑥+𝑘𝑦𝑦)}𝑑𝑥𝑑𝑦. (4.33) Let again 𝑢=(𝑢_𝑥,𝑢_𝑦,𝑢_𝑧)^𝑇=𝑢(𝑥,𝑦,𝑧,𝑡) be the displacement vector of a wave travelling in the plate in the absence of external forces (e.g., as a result of an initial excitation which now has stopped) such that it fulfills (3.9), (3.10), and (3.11). By formally applying the Fourier transform with respect to (𝑥,𝑦) to these equations, we see that 𝑢=𝑢(𝑘,𝑧,𝑡)=𝑢(𝑘_𝑥,𝑘_𝑦,𝑧,𝑡) then fulfills 𝜌_𝑙𝑢̈=𝐿^𝑇_𝑘𝐶_𝑙𝐿_𝑘𝑢,∀𝑧∈(𝐻_𝑙−1,𝐻_𝑙),𝑙=1,…,𝑁 (4.34) and the boundary and interface conditions (3.14) and (3.15). Suppose that for all 𝑘 and 𝑡 (or at least for the ones of interest) the functions 𝑧↦𝑢(𝑘,𝑧,𝑡) and 𝑧↦𝑢̈(𝑘,𝑧,𝑡) are in 𝐻¹. Analogously to the previous considerations, the weak form of (4.34), (3.14), and (3.15) is then equivalent to 𝑀𝑢̈=−𝑆_𝑘𝑢. (4.35) According to Proposition 4.2 and (2.11), we can represent 𝑢(𝑘,𝑧,𝑡) by a series with the eigenvectors 𝑢_𝑛(𝑘,𝑧) as basis functions 𝑢(𝑘,𝑧,𝑡)=^∞_𝑛=0𝛼_𝑛(𝑘,𝑡)𝑢_𝑛(𝑘,𝑧),𝑢̈(𝑘,𝑧,𝑡)=^∞_𝑛=0𝛼̈_𝑛(𝑘,𝑡)𝑢_𝑛(𝑘,𝑧). (4.36) We remark that for different 𝑘 we have different basis functions 𝑢_𝑛=𝑢_𝑛(𝑘), hence the dependence 𝛼_𝑛=𝛼_𝑛(𝑘); but for fixed 𝑘 we fix this basis for all times 𝑡, hence the dependence 𝛼_𝑛=𝛼_𝑛(𝑘,𝑡). Furthermore we assume that 𝑢̈ can indeed be represented with 𝛼̈_𝑛; that is, summation and time derivation commute. Applying the orthogonality relations in (4.21) to (4.35), we get for fixed 𝑘 the decoupled system of ordinary differential equations in 𝑡: 𝛼̈_𝑛1(𝑘,𝑡)𝛾+𝜔_𝑛(𝑘)²=−𝛼_𝑛𝜔(𝑘,𝑡)^𝑛(𝑘)²𝛾+𝜔_𝑛(𝑘)²⇔𝛼̈_𝑛(𝑘,𝑡)+𝜔_𝑛(𝑘)²𝛼_𝑛(𝑘,𝑡)=0, (4.37)

with general solutions 𝛼_𝑛(𝑘,𝑡)=𝑎_𝑛(𝑘)𝑒^{𝑖𝜔𝑛(𝑘)𝑡}+𝑏_𝑛(𝑘)𝑒^{−𝑖𝜔𝑛(𝑘)𝑡}. (4.38) Assuming that the following operations are valid, we get the mode decomposition 𝑢(𝑥,𝑦,𝑧,𝑡)=_ℝ2𝑢(𝑘_𝑥,𝑘_𝑦,𝑧,𝑡)𝑒^{𝑖(𝑘𝑥𝑥+𝑘𝑦𝑦)}𝑑𝑘_𝑥𝑑𝑘_𝑦=^∞_𝑛=0ℝ2(𝑎_𝑛(𝑘)𝑢_𝑛(𝑘,𝑧)𝑒^{𝑖𝜔𝑛(𝑘)𝑡}+𝑏_𝑛(𝑘)𝑢_𝑛(𝑘,𝑧)𝑒^{−𝑖𝜔𝑛(𝑘)𝑡})×𝑒^{𝑖(𝑘𝑥𝑥+𝑘𝑦𝑦)}𝑑𝑘_𝑥𝑑𝑘_𝑦. (4.39)

5. Numerical Results

The generalized eigenvalue problem can efficiently be solved by a suitable finite-element discretization, leading to a finite-dimensional generalized eigenvalue problem of the same form of (4.10) with mass matrix 𝑀 and stiffness matrix 𝑆_𝑘; see [7–9]. A comprehensive treatment of the mathematical theory and practice of FEM can be found in the study by Ern and Guermond in [21]. For theoretical results concerning order-preserving convergence of the eigenvalues of the discretized problem to the eigenvalues of the infinite-dimensional problem, we refer to the studies by Yang and Chen in [22].

As an application we want to examine how dispersion curves vary with an increasing degree of anisotropy. To do so, we consider a transversely isotropic material with 𝑦-axis as the axis of rotational symmetry. We note that unidirectional fibre-reinforced composites can often be approximately modelled by transversely isotropic material so that symmetry-axis and fibre direction coincide. For a description of the elastic material constants in this case we refer to the studies by Ambartsumyan in [23] and Altenbach et al. in [24]. We used as Young's moduli 𝐸₁=𝐸₃=70 GPa, 𝐸₂=𝜃𝐸₁, as Poisson’s ratios 𝜈₃₂=𝜈₃₁=𝜈₁₃=𝜈₁₂=0.3, 𝜈₂₃=𝜈₂₁=𝜈₃₂/𝜃, and as shear moduli 𝐺₂₃=𝐺₁₂=𝐺₁₃=𝐸₁/(2(1+𝜈₂₃)), where the ratio 𝜃=𝐸₂/𝐸₁ determines the degree of anisotropy. The material behaves isotropic for 𝜃=1 while for 𝜃>1 it becomes transversely isotropic with rotational symmetry around the 𝑦-axis. In this case the matrix 𝐶 of material constants can be written as 𝐶=(1−𝜈²₃₂/𝜃)𝐸₁Δ(𝜈₃₂+𝜈²₃₂)𝐸₁Δ(𝜈₃₂+𝜈²₃₂/𝜃)𝐸₁Δ000(1−𝜈²₃₂)𝜃𝐸₁Δ(𝜈₃₂+𝜈²₃₂)𝐸₁Δ000(1−𝜈²₃₂/𝜃)𝐸₁Δ𝐸000₁2+2𝜈₃₂𝐸00symmetric¹2+2𝜈₃₂0𝐸₁2+2𝜈₃₂, (5.1) with 2Δ∶=1−(𝜃+1)𝜈²₃₂−2𝜃𝜈³₃₂. (5.2)

The plate is assumed to have a thickness of 𝑑=1mm and a mass density of 𝜌=2700kg⋅m⁻³. We computed the circular frequency 𝜔(𝑘,𝜙) for a fixed circular wavenumber in dependence of the propagation direction 𝜙 and for different ratios 𝜃=𝐸₂/𝐸₁. Some results for the first four modes are displayed in Figure 2 which illustrate the increasingly anisotropic behaviour. In the strictly transversely isotropic case 𝜃>1, we labeled the modes according to their isotropic counterparts, but it becomes clear from Figure 3 that the distinction between pure 𝑆𝐻- and 𝑆-modes is growing less obvious as both kinds of modes get components in all directions.

Figure 2 

Direction dependence of circular frequency $\omega (k,\varphi )$$[{\text{s}}^{-1}]$ at circular wavenumber $k=1000$ ${\text{m}}^{-1}$ for different ratios $E_2/E_1$. As one can see, each mode has its own individual anisotropic behaviour; even the sequence of the wave modes ordered by their frequencies is angle dependant.

Figure 3 

Components of the displacement vector $u=(u_x,u_y,u_z)$ over the thickness of the plate at $\varphi =70^{\circ }$ of $S_0$- and $S{H}_1$-modes for different ratios $E_2/E_1$. Here $u_x$ and $u_y$ are the in-plane components in longitudinal and transversal directions, respectively. Obviously the distinction between shear-horizontal and in-plane Lamb-waves becomes less clear in an anisotropic material. This means especially that common measuring techniques quantifying the $z$ component will detect also the $SH$-modes.

Acknowledgment

The work of the authors is being supported by Deutsche Forschungsgemeinschaft (DFG) under Schu 1978/4-1.