# Semianalytical Approach for Single- and Multilayered Anisotropic Plates

If one refers books, research articles, and Section 5.3.1 on analytical methods to solve fundamental wave modes in multilayered anisotropic wave guide, it will immediately appear that the solution of guided wave modes in multilayered solid is a daunting task. Hence, many researchers including author Banerjee  have pursued a semianalytical approach to understand the guided wave dispersion in multilayered plate-like structure.

All multilayered situations could be handled by this semianalytical approach. By virtue of its name, it is self-explanatory that the semianalytical method uses partially analytical and partially numerical formulations. When the numerical approach is introduced. Finite Element Method (FEM) is most common to use and adopt for a problem. This method is no exception. Hence, finite element discretization is used in this method across the depth of the plate. Thus, this method is called Semi Analytical Finite Element (SAFE) method, or simply the SAFE. In SAFE approach, to model guided waves in composite structures assumes analytical guided wave propagation along the plate and performs a finite element discretization across the multilayer thickness of the plate. The SAFE formulation forms a stable eigenvalue problem for solving wave numbers and mode shapes in a given frequency range. In this book, numerical computation of wave field for NDE problems is discussed. Although there are many drawbacks discussed in Chapter 1, FEM is one of the prominent methods used in wave field modeling. Despite FEM being not discussed in this book explicitly, readers are recommended to refer any FEM text books for basic understanding of FEM, such as Reference .

Let us assume a multilayered anisotropic plate shown in Fig. 5.16. In Fig. 5.16, the conventional right-hand coordinate system used in Fig. 5.15 is flipped to define the thickness of the plate along the positive .Гз-axis from the top surface of the plate. ,vraxis is the direction of guided wave propagation designated with a wave vector к at frequency to. The cross section of the plate is on the x2-x} plane. Across depth, along x3 -axis, multinodded (say n noded, where n could be any number between 3 and infinity) one-dimensional (ID) elements that are dedicated for each layer are used to discretize the cross section of the plate as shown in Fig. 5.16 (Fig. 5.16 shows FIGURE 5.16 Semi-analytical modeling approach for simulating wave propagation in composite layered plate.

three nodded element for clarity). The discretization is also to describe the mode shapes of the guided wave modes in the plate. The harmonic displacement (и, ), stress (a,,), and strain (e,j) field components at any point (jc* ) of the waveguide in Cartesian coordinate system are expressed by where standard index notation described in Chapter 2 is 3 and is valid. Please note that the factor 2 with each shear strain component is added to satisfy the constitutive relation described between Eqs. (3.78) and (3.81). Hence, recalling Eq. (3.75), the constitutive equation of the material at a point is given by where C,yW is the stiffness or constitutive matrix, described in Section 3.9. Recalling Eq. (3.72) the stress-displacement relation at any material point can be written as However, let us expand Eq. (5.81) like in Appendix in Chapter 3 and express each strain component in Eq. (5.79), with respect to the displacements in Eq. (5.77) at any arbitrary material point as follows: By scrutinizing Eqs. (5.82.1) and (5.82.6), it is realized that all of them have the derivative terms with factors either 1 or 0. For example, in Eq. (5.82.1),

has factor 1, but Eq. (5.82.4) has factor 0 for Similarly, in Eq. (5.82.6), a|r has factor 0, but both Eqs. (5.82.4) and (5.82.5) have factor 1 for Hence, it is wise to construct three matrices as follows and to express the strain filed in a bit different form. Simplifying Eq. (5.83), one can write and further where T indicates transpose of the vector u, and

## Hamiltons Principle and the Governing Equation

Before going into any further detail, readers are recommended to read the section on classical mechanics in Section 3.10.4 in Chapter 3. Wave propagation in a composite plate can be considered as a dynamic system that is evolving over time and must abide by the fundamental physics and principles of a classical system. According to the Hamilton’s principle [62, 63], the system must satisfy the 8£ —> 0, which is the principle of least action in classical mechanics. That means that the system could have evolved only in one way my minimizing the action and the action £ in minimum. £ is called Lagrangian of a system, £ = £(«,, гл, м,.Д as described in Eq. (A4.5). Lagrangian of a system in the light of Section 3.8.3 in Chapter 3 is the difference of the kinetic energy density and the strain energy density (C = K-£) as described in Eq. (A4.15). Hence, the Hamilton’s principle can be written as where К is the kinetic energy density and £ is the strain energy density of the system as described in Eq. (A4.16) using index notation. In light of strain and stress vectors at a point described at the beginning of this section, the energy densities can be written in vector form as follows: where p is the material density, and the dot represents a time derivative. After integrating by parts, Eq. (5.88) can be written as Eq. (5.89) is the fundamental governing equation, which should be solved using SAFE approach. By now, readers are expected to see a pattern in solving wave propagation in solid and fluid, which should not change even if SAFE is used. That pattern is used in this chapter over and over again. The displacement filed is assumed to be harmonic along the propagation direction k. Similarly, here in this problem, the harmonic displacement function is assumed. However, here it is done in a bit different way for later convenience with discretization. As the thickness direction of the plate belongs to the jc2 — х.ч plane, and wave propagation is along the direction X, the phase term is specifically assigned for the X direction only and kept a generic two dimensional function for the x2 -x2 plane as written below. The displacement field over the x2-x3 plane, £/,•(x2,x2) a spatial function can be further expressed using discretized function along the .v3-axis with the values described at the nodal points (refer Fig. 5.16) only.

## Discretization of Plate Thickness

Across thickness, i.e. along дгз-axis, each lamina of the plate is discretized  with multinodded element. In Fig. 5.16, a three nodded element in any arbitrary layer is shown. The t/,(at2,jc3) function (i.e., U(x2,x2), U2(and U)(x2,x2) along X|, x2, and x}, respectively) in Eq. (5.90) is described by its nodal values at the node points, e.g., Uj at the y'-th node and respective discretization shape function Afj(x2,x2) for the y'-th node. Hence, displacements at any point on the element can be expressed in terms of the shape functions Jfj(x2,x3) and the nodal unknown displacements Uf(x2,x2). Hence, displacement at any point on the element (e) (Uj{exi,x2,xi,t)) is expressed by where and

## Element Strain Equation

Substituting the newly found expression for u(t,) written in Eq. (5.91) in the strain equation expressed in Eq. (5.85), strain vector for each element can take the following mathematical shape. Applying differential operator on the shape function matrix Af(x2,xi), the strain vector for each element will be In some cases, the above matrix in Eq. (5.95) is written as where

## Governing Wave Equation

Let us assume along the thickness direction of the plate has total Ne numbers of layers. Hence, there are total Ne numbers of n nodded elements in the plate. Assuming dvie) being the elemental volume of the e-th element and assembling all the elements in the light of Eq. (5.96), the Hamilton’s principle or the governing equation of the system described in Eq. (5.89) can be further modified into where C(<,) and pu'’ are the constitutive matrix and the density of respective e-th element. In this case, where multilayered plate is analyzed, material properties for each layer or each element (i.e., the e-th element) must be separately specified. The symbol U symbolizes that the elemental equations for each element are required to be assembled in a bigger matrix by matching the interfacial conditions. The assembly procedure is not discussed in this chapter and basic understanding of FEM  is assumed throughout this book. However, while presenting Spectral Element Method (SEM) in Chapter 10, the assembly procedure is briefly discussed.

Further substituting the strain vector in Eq. (5.96), the strain energy part of the governing equation takes the following form. Please note that being a nonzero term.

the exponential phase and time harmonic part is removed without the loss of the generality from the equation intentionally. Similarly, substituting the displacement vector in Eq. (5.91), the kinetic energy part of the governing equation takes the following form. Please note that like Eq. (5.99), being a nonzero term, the exponential phase and time harmonic part is removed from the equation intentionally, without the loss of the generality. Eq. (5.98), the governing equation of the system, has two parts and those two parts are now expanded in Eqs. (5.99) and (5.100). The discretized governing equation while satisfying the Hamilton’s principle can be further obtained by substituting Eqs. (5.99) and (5.100) in Eq. (5.98). Collecting the terms with integrals and rearranging the equation, governing equation will be where elemental expressions of the newly introduced terms are Next applying the assembling procedure and joining the nodes of two adjacent elements, the above equation will become where U is the global displacement vector of unknown nodal displacements across the depth of the plate. Like q(<,) written in Eq. (5.93) consists of the nodal displacements from each n nodded element, the U consists of nodal displacements for all the nodes in assembled system, where nodes at the element interface are shared and has only one unknown value assigned to each displacement (iq, m2, m3). With this understanding, the new terms introduced in Eq. (5.102) cane be expressed by From Eq. (5.102), it is easily realized that to satisfy the Hamilton’s principle, the integral should be equal to zero, and moreover, the integrand must be equal to zero. This should be valid for any arbitrary virtual displacement 5U. Hence, the governing equation of the system with wave propagating along .vraxis follows the following homogeneous equation Further, Eq. (5.104) can be written as follows :

## Eigen Value Problem: Wave Dispersion Solution and Phase Velocity

From the above expression in Eq. (5.105), it can be immediately realized that the governing equation has now become an eigenvalue problem. Here, the eigen values are also solved using a root-finding approach at a given frequency. First, lower and upper bounds of the frequency axis are selected and the axis is finely discretized to input different frequency values to the system to find the roots. If the dimension of vector U is assumed to be D (D = nNe-{L-1), where L is the number of layers in the plate), then at each frequency со, 2D eigenvalues and ID eigenvectors are obtained. The eigen values are obtained as pair of wave numbers that are both positive and negative (±A, = +(£,'+/£,'), where k[ is the real part of the wave number and A'[ is the imaginary part of the wave number), representing wave propagating along positive or negative xraxis as shown in Fig. 5.16. The wave numbers could be only real, pure, imaginary, and complex conjugates. Real w'ave numbers are the ones that represent the real propagating wave modes with no decaying term associated.

Complex conjugates always have both real and imaginary parts, w'hich signifies propagating wave modes w'ith decaying amplitude. This could be explained as follow's: In the exponential term in Eq. (5.91), ikX part contributes to the phase part of the propagating wave mode. When an eigen value for the wave number is found to be (&i = (k[ + ikl), ikX will result in ikX = (ik[x -k[x,). where the imaginary oscillatory part is contributed by the k[ real part of the wave number and the exponentially decaying part is contributed by the -k[ imaginary part of the wave number. Hence, complex conjugate wave numbers are representative of damped propagative waves decaying in the ±X directions. With similar explanation, if the eigen wave numbers are found to be pure imaginary, then they will only contribute to the decaying part of the wave mode and will not propagate; such wave modes are called evanescent wave modes. Based on our previous discussion in Chapter 4 and beginning of Chapter 5, the phase velocity of the propagating wave modes can be found as ## Dispersion Behavior

Here in this section, few sample dispersion plots are presented for multilayered plates. Fig. 5.17 shows the typical dispersion curves in a (a) eight-layered transversely isotropic composite w'ith all 0° layers, (b) (0/90/0/90)s eight-layered composite plate, (c) (45/-45/0/90)s composite laminates, and (d) (0/0/90/90)s composite laminates with the material properties for transversely isotropic material presented in Eq. (4.102). In order to substitute right material properties for 0° and 90°, the axis of symmetry was changed from Araxis to xvaxis, w'hich became their dominant fiber direction, respectively. All the dispersion curves presented are for the propagation direction along the Araxis as showm in Fig. 5.16. Please note that with the change in the direction of propagation of the wave (k vector) the dispersion behavior changes. The dispersion curves in Fig. 5.17 are obtained using LAMSS-Composites software developed by Laboratory for Active Materials and Smart Structures (LAMSS) laboratory at the University of South Carolina, Columbia, SC, USA. developed under the project funded by NASA Langley Research Center, contract no. NNL15AA16C.

## Group Velocity of Propagating Wave Modes

As discussed before in Chapter 4, the group velocity of the propagating wave modes can be found by taking derivatives of the frequency-wave number dispersion relations. Conventionally, group velocity is found based on the differences calculated in the wave number values at two adjacent points on the same mode but at two different frequencies. For example, from a dispersion solution, let us take two adjacent points, p and q, from a same mode. The difference in frequency and wave numbers is, say со,, -со,, and k4 -kr. Then the group velocity is cg = cov -оэгч -kp. Thus, it is apparent that the accuracy of the group velocity calculation heavily depends on the frequency discretization or the frequency resolution of the dispersion solution. Hence, to solve the eigen values from Eq. (105), the frequency axis must be finely discretized. An alternative method to find the group velocity was proposed by FIGURE 5.17 Typical dispersion curves in a) 8 layered transversely isotropic composite with all 0° layers, b) (0/90/0/90), 8 layered composite plate, с) (45/-45/0/90), composite laminates, d) (0/0/90/90), composite laminates with the material properties for transversely isotropic material presented in Eq. (4.102).

Svante Finnveden. . By virtue of this new method, without any contribution from adjacent points, the group velocity can be directly calculated at each (kt, со) point anywhere in the frequency-wave number solution space. The approach starts from the governing equation written in Eq. (5.104). Derivative of Eq. (5.104) with respect to the wave number along the propagation direction can be written as Upon simplification and rearrangement, the governing equation is Next, multiplying Eq. (5.108) by the transpose of the left eigenvectors U[, one can write From Eq. (5.104), &i2K| + ik{K} + K2 - orM = 0, left with only one choice to satisfy the following equation In Eq. (5.110), the derivative term is a scalar quantity. Hence, the group velocity can be found as Previously in Chapter 4, how the group velocity direction can be found from the known phase velocity profile (in fact from slowness profile) is discussed in detail. Readers are recommended to follow the discussion on finding the steering angle direction in Section 4.5.4 in Chapter 4.