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

Dispersion curves of elastic guided waves in plates can be efficiently computed by the StripElement 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.


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 finiteelement-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 infinitedimensional 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 C ∞ , 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.

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 H, H 1 , H 2 be real or complex Hilbert spaces with scalar products

2.1
Range and nullspace are denoted by R T and N T , respectively.We say that T ∈ L H is self-adjoint if T * T , and positive if If T is injective and has closed range, then it is bijective and hence continuously invertible.
Proof.For X ⊂ H, let X ⊥ {x ∈ H | x | y 0, for all y ∈ H} be the orthogonal complement of X.The assertion follows from which means that T is also surjective.
Proof.See, for example, the study by Schr öder in 17 .
We obtain the following.
Corollary 2.3.An operator T ∈ L H is injective and has closed range if there exists some constant c > 0 such that Proof.Together with the Cauchy-Schwartz inequality, we obtain from 2.5 v | Su H indeed defines a scalar product on H which induces an equivalent norm on H.With respect to this scalar product, the injective and compact operator M : S −1 M is self-adjoint and positive.Therefore 2.9 is equivalent to the standard eigenvalue problem Mu n λ n u n for the injective, compact, self-adjoint, and positive operator M ∈ L{H, • | • S }.

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 L is given as, with the constant matrices ii Displacement vector u is given as u u x, y, z, t u x , u y , u z T .

iv Strain-displacement relation is given as
Lu. 3.5 v Stress vector σ is given as vi And the generalized Hooke Law is given as with C C ij i,j 1,...,6 being the matrix of material constants.In the following, C is supposed to be real, symmetric, and positive definite spd .For an isotropic material we have, for example, Now consider a laminated plate of thickness 2H in direction of the z-axis middle z 0, top z H, bottom z −H and infinite in the x-y plane; see Figure 1.The plate consists of N layers.Layer l 1, . . ., N has thickness H and is supposed to consist of homogenous, anisotropic, elastic material with density ρ l and material constants C l .
Let u u x , u y , u z T u x, y, z, t 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 3.9 The traction-free boundary conditions on the top and bottom surfaces are Besides continuity of the displacement vector, the following interface conditions concerning continuity of the stresses are supposed to hold:

3.11
As ansatz for a wave mode, we take a plane harmonic wave travelling in the x-y plane u x, y, z, t u z e i k x x k y y ± ωt , 3.12 with z-dependent amplitude vector u z u x z , u y z , u z z T , real circular frequency ω, and real wave vector k k x , k y T .For such a wave mode, the wave equation 3.9 , boundary 3.10 and interface conditions 3.11 reduce to with the differential operator matrix L k given as

3.16
Obviously for complex conjugation we have 3.17 The question is "for which combinations of circular frequencies ω and wave vectors k do nontrivial solutions u of 3.13 , 3.14 , and 3.15 exist?"The answer leads to the dispersion relations ω k and is given in the next section.For better readability we will write ω instead of ω k in the following.

Weak Form of the Reduced Wave Equation and Related Eigenvalue Problem
We define the piecewise constant functions of density ρ and matrix of material constants C as For suitable virtual displacements v v z , which we define later, the weak form of 3.13 , 3.14 , and 3.15 can then be written as Mathematical Problems in Engineering 7 The next proposition shows that strong and weak solutions coincide for smooth enough functions.Proposition 4.1.Suppose that u u z is continuous and layerwise C

4.3
Let u fulfill this equation for all v ∈ C 1 −H, H .At first we fix a layer l 1, . . ., N and choose arbitrary functions v that have compact support in the interior H l−1 , H l of the layer.For such functions, 4.3 reduces to from which we infer that u 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

4.5
Now by successively choosing functions v that equal 1 in a vicinity of one of the points −H H 0 , H 1 , . . ., H N H and that are 0 everywhere else, we see that u 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 H 1 H 1 −H, H be the Sobolev space of all complex-valued square integrable functions on −H, H 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 H 1 also for H 1 3  H 1 × H 1 × H 1 , H 1 6 , . . ., and likewise for other spaces like the spaces of square integrable functions L 2 L 2 −H, H 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 the space H 1 is a complex Hilbert space.The space C 1 is continuously embedded and dense in H 1 and the inclusions H 1 → C and J : H 1 → L 2 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 C, respectively, define bijective, positive, self-adjoint, continuous linear operators on the respective L 2 -spaces.Furthermore, the differential operator matrix L k 3.16 defines a continuous linear operator L k : H 1 → L 2 .For functions u, v ∈ H 1 , which we assume in the following, the left-hand side of 4.2 can then be written as with the compact, injective, positive, and self-adjoint operator M : J * ρJ, and the right-hand side can be written as with the positive and self-adjoint operator S k : L * k CL k .Thus, a function u ∈ H 1 fulfills 4.2 for all v ∈ H 1 if and only if which is equivalent to the generalized eigenvalue problem with λ : γ ω 2 and the positive and self-adjoint operator S γ : γM S k for an arbitrary γ > 0. We will need this disturbance with γM to prove that the operator S γ is also bijective and hence can apply Proposition 2.4 to guarantee nontrivial solutions.By Lemma 2.1 and Corollary 2.3, bijectivity of S γ follows by showing that there exists some constant c > 0 such that

4.11
To show this, we write

4.12
Since layerwise C is an spd-matrix and ρ > 0, there is a constant c 1 > 0 such that pointwise a.e.
and thus we have and the Hermitian matrix where 4.17 are strictly positive.Hence there is a constant c 2 > 0 such that pointwise a.e.
and by the definition of the scalar product 4.6 we arrive at

4.19
Now we can apply Proposition 2.4 to the pair M, S γ and get the following.

Proposition 4.2.
There exists an increasing sequence 0 ≤ ω 2 n → ∞ and a sequence u n ∈ H 1 such that and each u ∈ H 1 can be uniquely represented by an 22 Proof.Let the assertions of Proposition 2.4 appropriately hold for the pair M, S γ .At first we observe that by 2.10 we have and hence ω 2 n : λ n −γ ≥ 0 and ω 2 n → ∞ increasingly.By 2.10 , on one hand, we then trivially have and, on the other hand, we further conclude that

4.25
Finally each u ∈ H 1 can be uniquely represented by an H 1 -convergent series Evaluating the scalar product of both sides of the above equation with S k u m , we get

4.27
The next proposition together with Proposition 4.1 finally shows that all these weak H 1 -solutions are indeed strong solutions.
Proof.We inductively show that the distributional derivatives of u can be identified with smooth functions.Expanding the operator L according to its definition and rearranging 4.2 give

4.28
We fix a layer l and take an arbitrary C ∞ -function v with compact support in H l−1 , H l .As C l is symmetric, for those v the left-hand side can be written as

4.29
We integrate by parts on the right-hand side to get

4.30
The operator L T z C l L z is invertible since C l is positive definite and L z is injective; hence the previous equation for all v is equivalent to where    We conclude this section by deriving a mode decomposition of general waves travelling in the plate.For a function u u x, y , we denote by u u k x , k y its Fourier transform in the wavenumber domain u k x , k y 1 2π 2 R 2 u x, y e −i k x x k y y dx dy.

4.33
Let again u u x , u y , u z T u x, y, z, t 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 x, y to these equations, we see that u u k, z, t u k x , k y , z, t then fulfills × e i k x x k y y dk x dk y . 4.39

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 M and stiffness matrix S k ; 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 y-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 E 1 E 3 70 GPa, E 2 θE 1 , as Poisson's ratios ν 32 ν 31 ν 13 ν 12 0.3, ν 23 ν 21 ν 32 /θ, and as shear moduli G 23 G 12 G 13 E 1 / 2 1 ν 23 , where the ratio θ E 2 /E 1 determines the degree of anisotropy.The material behaves isotropic for θ 1 while for θ > 1 it becomes transversely isotropic with rotational symmetry around the y-axis.In this case the matrix C of material constants can be written as

Mathematical Problems in Engineering
The plate is assumed to have a thickness of d 1 mm and a mass density of ρ 2700 kg • m −3 .We computed the circular frequency ω k, φ for a fixed circular wavenumber in dependence of the propagation direction φ and for different ratios θ E 2 /E 1 .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 SH-and S-modes is growing less obvious as both kinds of modes get components in all directions.

1 Figure 2 :
Figure 2: Direction dependence of circular frequency ω k, φ s −1 at circular wavenumber k 1000 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.

E 2 /E 1 3 Figure 3 :
Figure3: Components of the displacement vector u u x , u y , u z over the thickness of the plate at φ 70 • of S 0 -and SH 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.
• | • .By L H 1 , H 2 we denote the space of continuous linear operators from H 1 to H 2 and set L H : L H, H .For T ∈ L H 1 , H 2 we denote by T * ∈ L H 2 , H 1 the adjoint operator defined by Su n H . Especially, v | u S : v | Su H defines a scalar product on H which induces an equivalent norm on H, and u n is an orthonormal basis with respect to this scalar product.
2H , ∀u ∈ H.2.6With Lemma 2.2, we can conclude the assertion.For M, S ∈ L H , consider the generalized eigenvalue problem l 1, . . ., N. Then u fulfills 4.2 for all v ∈ C 1 −H, H if and only if u is a solution of 3.13 , 3.14 , and 3.15 .Proof.Note that we shortly write u| H l−1 ,H l ∈ C 2 H l−1 , H l meaning that there is a function u ∈ C 2 H l−1 , H l such that u| H l−1 ,H l u| H l−1 ,H l .By the smoothness assumptions on u, relation 3.17 , and layerwise partial integration on the right-hand side of 4.2 , we see that for all v ∈ C 1 −H, H 4.2 is equivalent to Equation 4.31 states that the second distributional derivative∂ 2 z u| H l−1 ,H l of u restricted to H l−1 , H l can be identified with A l u| H l−1 ,H l .Since u ∈ H 1 , we have A l u| H l−1 ,H l ∈ L 2 H l−1 , H l .Consequently, we may assume that u| H l−1 ,H l ∈ H 2 H l−1 , H l and hence u| H l−1 ,H l ∈ C 1 H l−1 , H l .From this we in turn infer that A l u| H l−1 ,H l ∈ C H l−1 , H l , and 4.31 ensures that ∂ 2 z u| H l−1 ,H l ∈ C H l−1 , H l ; that is, u| H l−1 ,H l ∈ C 2 H l−1 , H l .Repeating this argument proves the assertion.
Analogously to the previous considerations, the weak form of 4.34 , 3.14 , and 3.15 is then equivalent to We remark that for different k we have different basis functions u n u n k , hence the dependence α n α n k ; but for fixed k we fix this basis for all times t, hence the dependence α n α n k, t .Furthermore we assume that ¨ u can indeed be represented with αn ; that is, summation and time derivation commute.Applying the orthogonality relations in 4.21 to 4.35 , we get for fixed k the decoupled system of ordinary differential equations in t: l 1, ...,N  4.34and the boundary and interface conditions 3.14 and 3.15 .Suppose that for all k and t or at least for the ones of interest the functions z → u k, z, t and z → ¨ u k, z, t are in H 1 .