Numerical analysis of wave propagating in a periodic layered structure

When a periodic external vibration load is applied to a periodic structure in which the same layered structures are stacked as shown in Fig. 1, the wave field in the structure must become a periodic motion due to the periodic boundary condition. This section describes formulation of the wave propagation governed by the periodic boundary condition and wave equation for isotropic media.

Zoom In Zoom Out Reset image size

Consider periodic external forces and a periodic layered structure with fundamental layers of thickness l that consist of isotropic homogeneous regions and intermediate scatterer as shown in Fig. 1. It is assumed that the cross-sectional structure in the x–y plane is periodic at the distance of l in the y direction, and the external vibration distribution also maintains the periodicity of the distance l. Now, the external force and the structure are uniform in the z direction, which is perpendicular to the x–y plane, and plane strain is assumed. Then vibration on the x–y plane is considered. At this time, using the scalar potential $\phi$ and the z component of the vector potential $\psi ,$ the wave equations for longitudinal and transverse waves in the isotropic homogeneous regions can be written as ³⁰⁾

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}\phi }{\partial {t}^{2}}={c}_{L}^{2}\left(\displaystyle \frac{{\partial }^{2}\phi }{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}\phi }{\partial {y}^{2}}\right),\,\displaystyle \frac{{\partial }^{2}\psi }{\partial {t}^{2}}={c}_{T}^{2}\left(\displaystyle \frac{{\partial }^{2}\psi }{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}\psi }{\partial {y}^{2}}\right),\end{eqnarray} \tag{ 1 }$$

where t is time, and ${c}_{L}$ and ${c}_{T}$ are the longitudinal wave velocity and transverse wave velocity in the isotropic homogeneous regions, respectively. Assuming a harmonic vibration field with an angular frequency ω propagating with the wave number ${k}_{\alpha }$ of the longitudinal wave component and the wave number ${k}_{\beta }$ of the transverse wave component in the x direction, the solution of Eq. (1) can be written as

$$\begin{eqnarray*}&&\phi ={C}_{1}{e}^{-i\alpha y}{e}^{i{k}_{\alpha }x}{e}^{-i\omega t}+{C}_{2}{e}^{i\alpha y}{e}^{i{k}_{\alpha }x}{e}^{-i\omega t},\end{eqnarray*}$$

$$\begin{eqnarray}&&\psi ={D}_{1}{e}^{-i\beta y}{e}^{i{k}_{\beta }x}{e}^{-i\omega t}+{D}_{2}{e}^{i\beta y}{e}^{i{k}_{\beta }x}{e}^{-i\omega t},\end{eqnarray} \tag{ 2 }$$

where ${C}_{1},\,{C}_{2},\,{D}_{1},$ and ${D}_{2}$ are arbitrary constants, and i is an imaginary unit. The first term of the first equation in Eq. (2) represents a plane wave having a wavenumber component ${k}_{\alpha }$ and $-\alpha$ in the x and y directions, respectively, and the second term represents a plane wave having a wavenumber component ${k}_{\alpha }$ and $+\alpha$ in the x and y directions. Because the phases must match at every y position in a distance of l in the periodic structure as shown in Fig. 1 even in the intermediate scatterer region, the following equations hold

$$\begin{eqnarray}&&{\alpha }_{n}l=2n\pi ,\,{\beta }_{n}l=2n\pi ,\end{eqnarray} \tag{ 3 }$$

where $n$ is an arbitrary integer, ${\alpha }_{n}$ and ${\beta }_{n}$ are the y component of wave number for the order of n. The wavenumber component in the x-direction with respect to the order n is also written as,

$$\begin{eqnarray*}&&{k}_{{\alpha }_{n}}=\pm \sqrt{{\omega }^{2}/{c}_{L}^{2}-{\alpha }_{n}^{2}}=\pm \sqrt{{\omega }^{2}/{c}_{L}^{2}-{(2n\pi /l)}^{2}},\end{eqnarray*}$$

$$\begin{eqnarray}&&{k}_{{\beta }_{n}}=\pm \sqrt{{\omega }^{2}/{c}_{T}^{2}-{\beta }_{n}^{2}}=\pm \sqrt{{\omega }^{2}/{c}_{T}^{2}-{(2n\pi /l)}^{2}}.\,\end{eqnarray} \tag{ 4 }$$

That is, due to the above constraints under the periodic boundary condition as shown in Fig. 1, the general solution of Eq. (1) can be expressed by the sum with respect to n as follows

$$\begin{eqnarray}&&\begin{array}{l}\phi =\displaystyle \sum _{n=-\infty }^{+\infty }{A}_{n}{e}^{i{\alpha }_{n}y}{e}^{i{k}_{{\alpha }_{n}}x}{e}^{-i\omega t}+\displaystyle \sum _{n=-\infty }^{+\infty }A{{\prime} }_{n}{e}^{i{\alpha }_{n}y}{e}^{-i{k}_{{\alpha }_{n}}x}{e}^{-i\omega t},\\ \psi =\displaystyle \sum _{n=-\infty }^{+\infty }{B}_{n}{e}^{i{\beta }_{n}y}{e}^{i{k}_{{\beta }_{n}}x}{e}^{-i\omega t}+\displaystyle \sum _{n=-\infty }^{+\infty }{{B}^{^{\prime} }}_{n}{e}^{i{\beta }_{n}y}{e}^{-i{k}_{{\beta }_{n}}x}{e}^{-i\omega t},\end{array}\end{eqnarray} \tag{ 5 }$$

where ${A}_{n},$ $A{{\prime} }_{n},$ ${B}_{n},$ and ${{B}^{{\prime} }}_{n}$ are arbitrary constants.

For an isotropic periodic layered structure, the scalar potential $\phi$ and vector potential $\psi$ can be expressed as the summation of plane wave modes with different propagation directions. Because the expressions of $\phi$ and $\psi$ are only different in the velocities ${c}_{L}$ and ${c}_{T},$ now we consider only the transverse wave-derived $\psi$ that is mainly discussed in the later numerical analyses. When considering only the wave in the +x direction, the first term remains in Eq. (5) as follows

$$\begin{eqnarray}&&\psi =\displaystyle \sum _{n=-\infty }^{+\infty }{\psi }_{n},\,{\psi }_{n}={B}_{n}{e}^{i{\beta }_{n}y}{e}^{i{k}_{{\beta }_{n}}x}{e}^{-i\omega t},\end{eqnarray} \tag{ 6 }$$

where the wave number ${k}_{{\beta }_{n}}$ in the x direction takes a real number or a pure imaginary number depending on the sign of radicand of Eq. (4). Figure 2 is a dispersion relation of the second equation in Eq. (4) showing with the vertical axis representing the normalized frequency $fl/{c}_{T}$ and the horizontal axis representing the real and imaginary parts of the normalized wavenumber ${k}_{{\beta }_{n}}l/2\pi ,$ where $f$ is the frequency and $f=\omega /2\pi .$ If the frequency satisfies $fl/{c}_{T}\gt n,$ the wave number in the x direction ${k}_{{\beta }_{n}}$ becomes a real number, forming a wave front towards the x direction. Then the angle of wave direction with respect to the x axis is expressed as,

$$\begin{eqnarray}&&{\theta }_{n}={\tan }^{-1}\left(\displaystyle \frac{{\beta }_{n}}{{k}_{{\beta }_{n}}}\right).\end{eqnarray} \tag{ 7 }$$

The right sketches of Fig. 2 show the wave front for $n=0,\,+1,\,+2.$ All of these plane wave modes propagate at the transverse wave speed of ${c}_{T}$ (wavelength $\lambda ={c}_{T}/f$). On the other hand, when $fl/{c}_{T}\lt n,$ ${k}_{{\beta }_{n}}$ is a pure imaginary number and the wave field becomes attenuated as x increases. For example, if square or elliptical periodic obstacles are located as shown in Fig. 1, such non-propagation modes are indispensable for representing the wave field scattered from the obstacles. Strictly speaking, a large number of modes of the longitudinal wave component $\phi$ are also superposed in the scattering wave field as well as the propagating and non-propagating modes of transverse waves.

Zoom In Zoom Out Reset image size

Considering the wave field propagating the periodic structure as a superposition of modes propagating (or attenuating) in the x direction as represented by Eq. (5), all modes are orthogonal. Therefore, the mode conversion characteristics of the waves transmitted through the scattering region can be evaluated by using the energy flux of each mode propagating through the cross section in the y direction.

The energy flux that passes through one layer of the periodic structure in the +x direction can be expressed as

$$\begin{eqnarray}&&\begin{array}{c}E=-\frac{1}{2}\left[{\int }_{0}^{l}{\boldsymbol{\sigma }}{\left(\dot{{\boldsymbol{u}}}\right)}^{\ast }\,\cdot \,{\boldsymbol{n}}dy\right]\\ \,=-\frac{1}{2}\left[{\int }_{0}^{l}\left\{{\sigma }_{x}{\left(\dot{{u}_{x}}\right)}^{\ast }+{\tau }_{xy}{\left(\dot{{u}_{y}}\right)}^{\ast }\right\}dy\right],\end{array}\end{eqnarray} \tag{ 8 }$$

where ${\boldsymbol{\sigma }}$ is the stress tensor, $\dot{{\boldsymbol{u}}}$ is the first derivative of the displacement vector ${\boldsymbol{u}}$ with respect to time, ${\boldsymbol{n}}$ is the unit vector in the +x direction, and * is the complex conjugate. ${\sigma }_{x}$ and ${\tau }_{xy}$ are normal stress and shear stress in the x direction, and ${u}_{x}$ and ${u}_{y}$ are displacements in the x and y directions. Using the first term in the +x direction of Eq. (5), each of these components required for the calculation of the energy flux can be written as

$$\begin{eqnarray*}&&{u}_{x}=\displaystyle \frac{\partial \phi }{\partial x}+\displaystyle \frac{\partial \psi }{\partial y}=i\displaystyle \sum _{n=-\infty }^{+\infty }{k}_{{\alpha }_{n}}{\phi }_{n}+i\displaystyle \sum _{n=-\infty }^{+\infty }{\beta }_{n}{\psi }_{n},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\phi }_{n}={A}_{n}{e}^{i{\alpha }_{n}y}{e}^{i{k}_{{\alpha }_{n}}x}{e}^{-i\omega t},\,{\psi }_{n}={B}_{n}{e}^{i{\beta }_{n}y}{e}^{i{k}_{{\beta }_{n}}x}{e}^{-i\omega t},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{u}_{y}=\displaystyle \frac{\partial \phi }{\partial y}-\displaystyle \frac{\partial \psi }{\partial x}=i\displaystyle \sum _{n=-\infty }^{+\infty }{\alpha }_{n}{\phi }_{n}-i\displaystyle \sum _{n=-\infty }^{+\infty }{k}_{{\beta }_{n}}{\psi }_{n},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\sigma }_{x}=\rho {c}_{L}^{2}\left(\displaystyle \frac{{\partial }^{2}\phi }{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}\phi }{\partial {y}^{2}}\right)-2\rho {c}_{T}^{2}\left(\displaystyle \frac{{\partial }^{2}\phi }{\partial {y}^{2}}-\displaystyle \frac{{\partial }^{2}{\psi }_{z}}{\partial x\partial y}\right)\end{eqnarray*}$$

$$\begin{eqnarray*}&&=\displaystyle \sum _{n=-\infty }^{+\infty }\left\{-\rho {c}_{L}^{2}\left({k}_{{\alpha }_{n}}^{2}+{\alpha }_{n}^{2}\right)+2\rho {c}_{T}^{2}{\alpha }_{n}^{2}\right\}{\phi }_{n}-2\rho {c}_{T}^{2}\displaystyle \sum _{n=-\infty }^{+\infty }{k}_{{\beta }_{n}}{\beta }_{n}{\psi }_{n},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\tau }_{xy}=\rho {c}_{T}^{2}\left(2\frac{{\partial }^{2}\phi }{\partial x\partial y}+\frac{{\partial }^{2}{\psi }_{}}{\partial {y}^{2}}-\frac{{\partial }^{2}{\psi }_{}}{\partial {x}^{2}}\right)\end{eqnarray*}$$

$$\begin{eqnarray}&&=-2\rho {c}_{T}^{2}\displaystyle \sum _{n=-\infty }^{+\infty }{k}_{{\alpha }_{n}}{\alpha }_{n}{\phi }_{n}+\rho {c}_{T}^{2}\displaystyle \sum _{n=-\infty }^{+\infty }\left({k}_{{\beta }_{n}}^{2}-{\beta }_{n}^{2}\right){\psi }_{n}.\end{eqnarray} \tag{ 9 }$$

Substituting Eq. (9) into (8), it can be written as

$$\begin{eqnarray*}&&E=\frac{\rho \omega l}{2}\displaystyle \sum _{\begin{array}{c}n=-\infty \\ {k}_{{\alpha }_{n}}:\mathrm{real}\end{array}}^{+\infty }\left\{{c}_{L}^{2}\left({k}_{{\alpha }_{n}}^{2}+{\alpha }_{n}^{2}\right)-4{c}_{T}^{2}{\alpha }_{n}^{2}\right\}{k}_{{\alpha }_{n}}{A}_{n}{A}_{n}^{* }\end{eqnarray*}$$

$$\begin{eqnarray*}&&+i\frac{\rho \omega l}{2}\displaystyle \sum _{\begin{array}{c}n=-\infty \\ {k}_{{\alpha }_{n}}:\mathrm{imaginary}\end{array}}^{+\infty }{c}_{L}^{2}\left({k}_{{\alpha }_{n}}^{{\prime} 2}-{\alpha }_{n}^{2}\right){k}_{{\alpha }_{n}}^{{\prime} }{A}_{n}{A}_{n}^{* }\end{eqnarray*}$$

$$\begin{eqnarray*}&&+\frac{\rho \omega l}{2}\displaystyle \sum _{\begin{array}{c}n=-\infty \\ {k}_{{\beta }_{n}}:\mathrm{real}\end{array}}^{+\infty }{c}_{T}^{2}\left({k}_{{\beta }_{n}}^{2}+{\beta }_{n}^{2}\right){k}_{{\beta }_{n}}{B}_{n}{B}_{n}^{* }\end{eqnarray*}$$

$$\begin{eqnarray}&&+i\frac{\rho \omega l}{2}\displaystyle \sum _{\begin{array}{c}n=-\infty \\ {k}_{{\beta }_{n}}:\mathrm{imaginary}\end{array}}^{+\infty }{c}_{T}^{2}\left({k}_{{\beta }_{n}}^{{\prime} 2}+3{\beta }_{n}^{2}\right){k}_{{\beta }_{n}}^{{\prime} }{B}_{n}{B}_{n}^{* }.\end{eqnarray} \tag{ 10 }$$

The energy flux can be written separately as the sum of the modes with real wave numbers ${k}_{{\alpha }_{n}}$ and ${k}_{{\beta }_{n}}$ (propagating mode) and the sum of the modes with pure imaginary wave numbers ${k}_{{\alpha }_{n}}$ and ${k}_{{\beta }_{n}}$ (evanescent mode). Now, because ${k}_{{\alpha }_{n}}^{{\prime} }$ and ${k}_{{\beta }_{n}}^{{\prime} }$ are real numbers represented by ${k}_{{\alpha }_{n}}=i{k}_{{\alpha }_{n}}^{{\prime} }\,\,$and ${k}_{{\beta }_{n}}=i{k}_{{\beta }_{n}}^{{\prime} },$ the four summation terms are real, respectively. Since the time average of the energy flux passing through the cross section ${E}^{AVE}$ is the real part of Eq. (10), it is represented by the sum of only the propagation modes as,

$$\begin{eqnarray*}&&{E}^{AVE}=\mathrm{Re}\left(E\right)=\displaystyle \sum _{\begin{array}{c}n=-\infty \\ {k}_{{\alpha }_{n}}:\mathrm{real}\end{array}}^{+\infty }{E}_{n}^{L}+\displaystyle \sum _{\begin{array}{c}n=-\infty \\ {k}_{{\beta }_{n}}:\mathrm{real}\end{array}}^{+\infty }{E}_{n}^{T}\,,\end{eqnarray*}$$

$$\begin{eqnarray*}&&{E}_{n}^{L}=\displaystyle \frac{\rho \omega l}{2}\left\{{c}_{L}^{2}\left({k}_{{\alpha }_{n}}^{2}+{\alpha }_{n}^{2}\right)-4{c}_{T}^{2}{\alpha }_{n}^{2}\right\}{k}_{{\alpha }_{n}}{A}_{n}{A}_{n}^{* },\end{eqnarray*}$$

$$\begin{eqnarray}&&{E}_{n}^{T}=\displaystyle \frac{\rho \omega l}{2}{c}_{T}^{2}\left({k}_{{\beta }_{n}}^{2}+{\beta }_{n}^{2}\right){k}_{{\beta }_{n}}{B}_{n}{B}_{n}^{* }.\end{eqnarray} \tag{ 11 }$$

Here, since ${E}_{n}^{L}\,\,$and ${E}_{n}^{T}$ depend only on the amplitudes ${A}_{n}$ and ${B}_{n}$ of the nth mode, respectively, and are not affected by the other modes, they represent the time average of the energy flux of the nth mode derived from longitudinal and transverse waves, respectively.

The squares of the amplitudes of each mode, ${A}_{n}{A}_{n}^{* }$ and ${B}_{n}{B}_{n}^{* },$ can be extracted from the displacement distribution in the y direction as follows by using the orthogonality of the modes. Let the integral of the product of the displacement ${u}_{x},\,{u}_{y}$ and ${e}^{-i{\alpha }_{n}y}$ be ${u}_{x\alpha n},\,{u}_{y\alpha n},$ respectively, as

$$\begin{eqnarray*}&&{u}_{x\alpha n}\equiv {\int }_{0}^{l}{e}^{-i{\alpha }_{n}y}{u}_{x}dy=\,i{k}_{{\alpha }_{n}}l{A}_{n}{e}^{i{k}_{n}x}{e}^{-i\omega t}\,,\end{eqnarray*}$$

$$\begin{eqnarray}&&{u}_{y\alpha n}\equiv {\int }_{0}^{l}{e}^{-i{\alpha }_{n}y}{u}_{y}dy=\,i{\alpha }_{n}l{A}_{n}{e}^{i{k}_{\beta n}x}{e}^{-i\omega t}.\end{eqnarray} \tag{ 12 }$$

The squares of the complex values are

$$\begin{eqnarray}&&{u}_{x\alpha n}{u}_{x\alpha n}^{* }={k}_{\alpha n}^{2}{l}^{2}{A}_{n}{A}_{n}^{* }\,,\,{u}_{y\alpha n}{u}_{y\alpha n}^{* }={\alpha }_{n}^{2}{l}^{2}{A}_{n}{A}_{n}^{* }\,,\end{eqnarray} \tag{ 13 }$$

then, the squares of the amplitudes can be obtained as,

$$\begin{eqnarray}&&{A}_{n}{A}_{n}^{* }={u}_{x\alpha n}{u}_{x\alpha n}^{* }/{k}_{\alpha n}^{2}{l}^{2}\,{\rm{or}}\,{A}_{n}{A}_{n}^{* }={u}_{y\alpha n}{u}_{y\alpha n}^{* }/{\alpha }_{n}^{2}{l}^{2}\,.\end{eqnarray} \tag{ 14 }$$

Similarly, ${B}_{n}{B}_{n}^{* }$ can be obtained as,

$$\begin{eqnarray}&&{B}_{n}{B}_{n}^{* }={u}_{x\beta n}{u}_{x\beta n}^{* }/{\beta }_{n}^{2}{l}^{2}\,{\rm{or}}\,{B}_{n}{B}_{n}^{* }={u}_{y\beta n}{u}_{y\beta n}^{* }/{k}_{\beta n}^{2}{l}^{2}\,.\,\end{eqnarray} \tag{ 15 }$$

Inserting these values into the second and third equations of Eq. (11), the time average of the energy flux of the nth mode can be calculated.

As described above, all modes propagating in the $\pm x$ directions in a periodic layered structure satisfy the condition of Eq. (4). In the vicinity of the scattering obstacles, the wave field is expressed as a superposition of propagation modes and an infinite number of evanescent modes as shown in Eq. (5). On the other hand, in the regions apart from the scattering obstacles, the wave field can be expressed as the sum of a finite number of propagating modes. Therefore, if the scattering region is designed so that a certain propagation mode becomes dominant, direction-controlled transmitted waves and reflected waves can be obtained. This section describes the direction control of ultrasonic wave in a periodic structure that consists of homogeneous elastic regions and an intermediate plate stacked region.

Lamb waves propagating in the longitudinal direction in a thin plate have dispersion characteristics in which the phase velocity and group velocity differ depending on the thickness and frequency. Figure 3(a) is the phase velocity dispersion curve of A0 mode for an aluminum alloy plate with longitudinal and transverse wave velocities of 6400 m s⁻¹ and 3170 m s⁻¹, respectively. The lower horizontal axis is the product of frequency $f$ and plate thickness $d$ divided by the transverse wave velocity ${c}_{T},$ and the left vertical axis is the normalized phase velocity $c/{c}_{T},$ where the phase velocity of A0 mode of Lamb wave $c$ is divided by the transverse wave velocity ${c}_{T}.$ As seen in the figure, in the low $fd$ region $fd/{c}_{T}\lt 1,$ the phase velocity changes significantly with respect to the frequency and plate thickness. Considering this property and setting the structure in which flat plates with different thicknesses are stacked as shown in Fig. 3 in the middle, this section discusses control of the direction of transmitted waves by adjusting the thickness.

Zoom In Zoom Out Reset image size

Intensive analyses of Lamb waves in Ref. 31 showed that reflection and scattering are relatively small at the entrance of a plate stacked region when a transverse plane wave is incident from the left region. Therefore, each plate vibrates in the vertical direction as an A0 mode of Lamb wave. Because the A0 modes have different phase velocities in the plates with different thickness, the phase shifts at the exit, the right edge of the thin plate stacked region as shown in Fig. 3(b). Since this phase shift can be predicted from the dispersion curve, the phase shift at the exit can be controlled by adjusting the plate thickness ${d}_{i}$ and the length of stacked plate region L for the frequency ${f}_{0}$ (or angular frequency ${\omega }_{0}$) used. For example, by adjusting the phase at the exit to be $\phi \left(y\right)=2n\pi y/l+{\phi }_{0}$ (${\phi }_{0}$ is an arbitrary constant), the vibration distribution in the y direction becomes ${e}^{i\phi }={e}^{i{\phi }_{0}}{e}^{2in\pi y/l}={e}^{i{\phi }_{0}}{e}^{i{\beta }_{n}y},$ and then the nth-order plane wave mode can dominantly form in the right transmission region.

Table I shows the combination of plates in which the thickness and length of each stacked plate are set, considering that the transverse wave mode with $n=+1$ dominantly forms in the transmitted region at the frequency of 4.5 MHz. It was assumed that there are four plates per periodic layer. Figure 3(c) shows the phase difference between the entrance and exit of the stacked plate region represented by the following equation.

$$\begin{eqnarray}&&{\phi }_{i}=2\pi fL/{c}_{i}\,,\end{eqnarray} \tag{ 16 }$$

where ${c}_{i}$ ($i=1,2,3,4$) is the phase velocity of the Lamb wave A0 mode propagating in the ith plate as shown in Figs. 3(a) and 3(b) that depends on the plate thickness ${d}_{i}.$ $L$ is the length of the stacked plate region, 3.5 mm, and f is the frequency of 4.5 MHz.

Table I. Parameters for the thin plate stacked structure.

${d}_{1}({\rm{mm}})$	${d}_{2}({\rm{mm}})$	${d}_{3}({\rm{mm}})$	${d}_{4}({\rm{mm}})$	$L({\rm{mm}})$	$f({\rm{MHz}})$
0.32	0.28	0.22	0.18	3.5	4.5

The circles in Fig. 3(c) plots the phase ${\phi }_{i}$ against the center positions of the four plates. As shown in Fig. 3(b), $y=0$ and $y=l\,(=1.0\,{\rm{mm}})$ are set on the periodic boundary, and flat plates with thicknesses ${d}_{1},{d}_{2},\,{d}_{3}$ and ${d}_{4}$ are stacked in ascending order of y position. Since the Lamb wave of A0 mode vibrates in approximately uniform distribution in the plate thickness direction, it can be considered that the A0 mode has the phase shown in the bar graph in Fig. 3(c) at the exit of the stacked plate region. In Fig. 3(c), a straight line with a slope of $2n\pi /l\,$(l = 1 mm, n = +1) is added. This line means the ideal distribution for the formation of the transverse wave mode with $n=+1.$ Therefore, the single plane wave mode is expected to appear largely when the bar chart becomes closer to the line.

In order to verify the above-mentioned transverse wave propagation control, we calculated the wave propagation for such a periodic structure using the finite element software COMSOL Multiphysics. Figure 4 is a schematic figure of the calculation region. The upper and lower boundaries are periodic boundaries, and the left and right ends are absorbing regions so that reflected waves from the both ends can be sufficiently suppressed. The upper and lower boundaries of each plate in the intermediate scattering region where four plates are stacked are set as traction free boundaries force. The displacement vibrates in the cross section, and a plane strain state is assumed, which is the same assumptions as in Sect. 2. In addition, the calculation was performed in the frequency domain, and a steady-state solution for a certain frequency f was obtained. The calculation area is one cycle of the periodic structure (length l = 1.0 mm in the y direction) whose dimensions are as shown in Fig. 4. The absorbing region is divided into 148 912 on each side by a triangular element with a maximum side length of 0.01 mm, and the remaining area is divided into 1050 in the x direction and 100 in the y direction by a square element with a side length of 0.01 mm. The incident wave was a transverse plane wave towards the stacked plate region from the left by applying a body force in the y direction to an element located 1.0 mm away from the left end, "entrance" of the intermediate stacked region. Assuming that the material is an aluminum alloy, the longitudinal wave velocity and transverse wave velocity are the same as above, 6400 and 3170 m s⁻¹, and the density is 2700 kg m⁻³. The generated transverse wave passes through the entrance of the stacked region, propagates in the four plates as A0 mode of Lamb wave, and arrives at the exit (right end of the intermediate stacked region). These Lamb waves have different phases at the exit as described above.

Zoom In Zoom Out Reset image size

Figure 5 shows the wave propagation at 4.5 MHz under the above calculation conditions. The color is the vertical component of the rotation of the displacement vector ${\left[\text{rot}\,{\boldsymbol{u}}\right]}_{z},$ which represents the transverse wave component in the wave field. Four periodic calculation regions are displayed, stacking vertically in the y direction so that the wave propagation can be clearly seen. The vibration was attenuated in the absorbing regions at the left and right ends, and the reflected waves from the both ends was sufficiently suppressed. In addition, the wave propagation seems to be continuous at the periodic boundaries, showing the calculation with the periodic boundary conditions was correctly executed.

Zoom In Zoom Out Reset image size

In the right region where the transmitted wave is shown, the plane wave seems to propagate in the upper right direction. However, because modes other than $n=+1$ are also included, the broken striped pattern appeared in whole area. To analyze the transmitted waves quantitatively, the energy ratio of each mode in the transmitted wave to the total energy of the transmitted wave was calculated.

The energy of the transmitted wave ${E}^{Trans}$ can be calculated as the time average of the vibrating energy flux transmitted through the cross section by taking the real part of Eq. (8) in the transmission region. The time average of the energy flux in each mode can also be calculated according to the second and third equations of Eq. (11). In the frequencies and a structure used in this calculation, four propagation modes exist in the periodic structure: longitudinal wave with $n=0$ and transverse wave with $n=0,\pm 1.$ The time averaged energy fluxes ${E}_{0}^{L},\,{E}_{0}^{T}$ and ${E}_{\pm 1}^{T}$ were calculated using the displacement and stress in the transmission region.

Table II shows the energy of each mode as a ratio to the energy of the entire transmitted wave. ${E}^{Trans}$ and ${E}_{0}^{L}+{E}_{0}^{T}+{E}_{+1}^{T}+{E}_{-1}^{T}$ match with an error of 1% or less, showing that the energy calculation method for each mode in the previous section is correct. In addition, ${E}_{0}^{L}/{E}^{Trans}$ is very small, which means that the longitudinal wave component hardly appears due to the mode conversion when a transverse plane wave is incident. The transverse wave with $n=0$ appears largely in the transmission region, and ${E}_{0}^{T}/{E}^{Trans}$ exceeds 70%. On the other hand, since the energy ratio ${E}_{+1}^{T}/{E}^{Trans}$ of the desired transverse wave with $n=+1$ is about 30%, the conversion to the mode with $n=+1$ was not well performed in this stacked plate region. In the transmission region of Fig. 5, the transmitted wave appears to propagate diagonally upward (large arrow in the Fig. 5). However, the white lines connecting the same phase points differ from the wave front of the large arrow. In other words, as a result of the mixture of the mode with $n=0$ that propagates in the horizontal direction of the solid arrow and the mode with $n=+1$ that propagates in the diagonal direction, it propagates in the direction of the thick arrow. This is supported by the fact, shown in Table II, that the mode with $n=0$ and the mode with $n=+1$ are relatively large in the energy, and the mode with $n=-1$ is almost nonexistent.

Table II. Transmission energy rate for all propagating modes (%).

${E}_{0}^{L}/{E}^{Trans}\,$	${E}_{0}^{T}/{E}^{Trans}$	${E}_{+1}^{T}/{E}^{Trans}$	${E}_{-1}^{T}/{E}^{Trans}$
0.08	71.22	28.59	0.11

Although the thickness and length of each plate were determined so as to match the transverse wave mode with $n=+1$ by calculating the phase shift from the dispersion curve, the transmitted wave could not be made into a wave field dominated by the mode with $n=+1$ as described above. Therefore, in this section, the energy of each mode in the transmitted wave is calculated for the stacked plate region determined in the previous section when the frequency and a stacked plate length are changed.

Figure 6 shows the energy flux rate of each mode propagating in the transmission region when the frequency is changed in 0.01 MHz increments from 4 to 5 MHz. 4.50 MHz is the frequency used in the previous section. The ratio of each mode in the transmission region changes significantly due to subtle changes in frequency. At all frequencies, the longitudinal wave mode is close to 0, which indicates that the ratio of mode conversion to longitudinal waves by the stacked plate region is very small when a transverse plane wave is incident. Moreover, although the mode with $n=+1$ is larger than the $n=-1$ mode over the frequency range, the mode with $n=0$ is the largest in some frequency ranges, which denotes that frequency tuning is necessary for forming the transmitted wave of a single $n=+1$ mode. For example, the wave motion shown in Table II and Fig. 5 is given at the frequency where the transverse wave mode with $n=0$ is dominant.

Zoom In Zoom Out Reset image size

Figures 7(a) and 7(b) show the wave propagation at 4.20 MHz and 4.63 MHz, respectively, where the energy rate of the transverse wave mode with $n=+1$ is the largest and the second largest. As Fig. 5, the vertical component of the rotation of the displacement vector ${\left[\text{rot}\,{\boldsymbol{u}}\right]}_{z}$ is shown in color. Unlike in Fig. 5, the transverse wave mode with $n=0$ becomes smaller, and the mode with $n=+1$ exceeds 70%, then the transmitted wave with the dominant mode of $n=+1$ can be observed. However, because small amount of the $n=-1$ mode is superposed, the plane wave propagation towards the right lower direction can also be seen.

Zoom In Zoom Out Reset image size

Next, the energy of the transmitted mode was calculated for different lengths of the intermediate plate region. Figure 8 shows the energy flux rate of all propagating modes in the transmitted wave at the frequency of 4.5 MHz for different lengths of the stacked plate region in 0.1 mm increments from L = 3.0 mm to 4.0 mm. Since the mode with $n=+1$ is generally larger than the $n=-1$ mode in this length range, the wave propagating to the lower right is small.

Zoom In Zoom Out Reset image size

Figures 9(a) and 9(b) show the wave motion for L = 3.8 and 3.4 mm, where the ratio of the mode with $n=+1$ is the largest and the second largest in the length range of Fig. 8. The color represents ${\left[\text{rot}\,{\boldsymbol{u}}\right]}_{z}$ as in Figs. 5 and 7. Since both the $n=0$ mode and the $n=-1$ mode are small, a clear plane wave front propagating to the upper right is formed as a result of the fact that the $n=+1$ mode is dominant.

Zoom In Zoom Out Reset image size

The results shown above indicates that it is possible to realize a transmitted wave with the dominant mode of $n=+1$ by adjusting frequency and length of a stacked plate region from the stacked plate intermediate region originally designed with dispersion curves.

In this way, when a periodic external load is applied to a periodic layered structure, the transmitted wave field is limited to the propagation modes in which the wave number in the propagation direction expressed by Eq. (4) is a real value, and each mode propagates in its own angular direction. Therefore, each single mode can be measured by installing a receiving transducer that matches the unique direction.