The direct approach is used to compute the wave solutions of the unit cell. The dynamic of the unit cell is described with the dynamic stiffness matrix \mathbf{D} (we assume no damping, \mathbf{C}=0)

\underbrace{\left(\mathbf{K}-\omega^2 \mathbf{M}\right)}_{\mathbf{D}}\mathbf{q}=\mathbf{f}

The matrices \mathbf{K} and \mathbf{M} are the stiffness and mass matrix, \mathbf{f} is the load vector (see article Equation of Motion). The entries in the dynamic stiffness matrix of the unit cell can be ordered with respect to the left, interior and right degrees of freedom

(1) \begin{equation} \begin{bmatrix} \mathbf{D}_{LL} & \mathbf{D}_{LI} & \mathbf{D}_{LR} \\ \mathbf{D}_{IL} & \mathbf{D}_{II} & \mathbf{D}_{IR} \\ \mathbf{D}_{RL} & \mathbf{D}_{RI} & \mathbf{D}_{RR} \end{bmatrix} \begin{bmatrix} \mathbf{q}_{L}\\ \mathbf{q}_{I}\\ \mathbf{q}_{R} \end{bmatrix} = \begin{bmatrix} \mathbf{f}_{L}\\ \mathbf{f}_{I}\\ \mathbf{f}_{R} \end{bmatrix} \label{eq:orderDynamicStiffness} \end{equation}

The idea behind this approach is to directly derive the same relation that is described by the Floquet theroem (Periodic Structures:1) using the FE-model of the unit cell.

Theory

First, we condense displacements of the internal degrees of freedom assuming that there is no external load (\mathbf{f}_I=0). Thus the second line of \mathbf{Dq=f} is

\mathbf{D}_{IL} \mathbf{q}_L + \mathbf{D}_{II} \mathbf{q}_I + \mathbf{D}_{IR} \mathbf{q}_R = \mathbf{0}

Rearranging yields the displacements \mathbf{q}_I depending on the displacements \mathbf{q}_L and \mathbf{q}_R

(2) \mathbf{q}_I = -\mathbf{D}^{-1}_{II} (\mathbf{D}_{IL} \mathbf{q}_L + \mathbf{D}_{IR} \mathbf{q}_R )

Inserting (2) in the first and third line of (1) yields

\begin{align} \mathbf{f}_L &= \left(\mathbf{D}_{LL}-\mathbf{D}_{LI} \mathbf{D}^{-1}_{II} \mathbf{D}_{IL}\right) \mathbf{q}_L + \left(\mathbf{D}_{LR} - \mathbf{D}_{LI} \mathbf{D}^{-1}_{II} \mathbf{D}_{IR}\right) \mathbf{q}_R \\ \mathbf{f}_R &= \left(\mathbf{D}_{RL}-\mathbf{D}_{RI} \mathbf{D}^{-1}_{II} \mathbf{D}_{IL}\right) \mathbf{q}_L + \left(\mathbf{D}_{RR} - \mathbf{D}_{RI} \mathbf{D}^{-1}_{II} \mathbf{D}_{IR}\right) \mathbf{q}_R \end{align}

Summarizing, we get the reduced dynamic stiffness \tilde{\mathbf{D}}

(3) \begin{align} \begin{bmatrix} \tilde{\mathbf{D}}_{LL} & \tilde{\mathbf{D}}_{LR} \\ \tilde{\mathbf{D}}_{RL} & \tilde{\mathbf{D}}_{RR} \end{bmatrix} \begin{bmatrix} \mathbf{q}_L \\ \mathbf{q}_R \end{bmatrix} = \begin{bmatrix} \mathbf{f}_L \\ \mathbf{f}_R \end{bmatrix} \end{align}
\begin{align} \tilde{\mathbf{D}}_{LL} = \mathbf{D}_{LL} - \mathbf{D}_{LI} \mathbf{D}^{-1}_{II} \mathbf{D}_{IL} \\ \tilde{\mathbf{D}}_{RL} = \mathbf{D}_{RL} - \mathbf{D}_{RI} \mathbf{D}^{-1}_{II} \mathbf{D}_{IL} \\ \tilde{\mathbf{D}}_{LR} = \mathbf{D}_{LR} - \mathbf{D}_{LI} \mathbf{D}^{-1}_{II} \mathbf{D}_{IR} \\ \tilde{\mathbf{D}}_{LL} = \mathbf{D}_{RR} - \mathbf{D}_{RI} \mathbf{D}^{-1}_{II} \mathbf{D}_{IR} \\ \label{eq:condensedS} \end{align}

Reordering the first line of (3) to resolve for \mathbf{q}_R and inserting this into the second line, one receives the following relation

\begin{align} \begin{bmatrix} -\tilde{\mathbf{D}}^{-1}_{LR} \tilde{\mathbf{D}}_{LL} & \tilde{\mathbf{D}}^{-1}_{LR} \\ \tilde{\mathbf{D}}_{RL} - \tilde{\mathbf{D}}_{RR} \tilde{\mathbf{D}}^{-1}_{LR} \tilde{\mathbf{D}}_{LL} & \tilde{\mathbf {D}}_{RR} \tilde{\mathbf{D}}^{-1}_{LR} \end{bmatrix} \begin{bmatrix} \mathbf{q}_L \\ \mathbf{f}_L \end{bmatrix} = \begin{bmatrix} \mathbf{q}_R \\ \mathbf{f}_R \end{bmatrix} \end{align}

To get the same relation as with the Floquet theorem, the second line is multiplied with a minus sign (remember \mathbf{q}_R^n=\mathbf{q}_L^{n+1}, while \mathbf{f}_R^n=-\mathbf{f}_L^{n+1} from Floquet Theorem)

\begin{align} \underbrace{ \begin{bmatrix} -\tilde{\mathbf{D}}^{-1}_{LR} \tilde{\mathbf{D}}_{LL} & \tilde{\mathbf{D}}^{-1}_{LR} \\ -\tilde{\mathbf{D}}_{RL} + \tilde{\mathbf{D}}_{RR} \tilde{\mathbf{D}}^{-1}_{LR} \tilde{\mathbf{D}}_{LL} & -\tilde{\mathbf{D}}_{RR} \tilde{\mathbf{D}}^{-1}_{LR} \end{bmatrix} }_{\mathbf{T}} \begin{bmatrix} \mathbf{q}_L^n \\ \mathbf{f}_L^n \end{bmatrix} = \begin{bmatrix} \mathbf{q}_R^n \\ -\mathbf{f}_R^n \end{bmatrix} = \begin{bmatrix} \mathbf{q}_L^{n+1} \\ \mathbf{f}_L^{n+1} \end{bmatrix} \end{align}

\mathbf{T} is called transfer matrix. It describes the relation of the state variables of two neighboring cells. With the transfer matrix \mathbf{T} and the Floquet theorem (Periodic Structures:1), there are two expressions describing the relation between the displacements and the forces of to neighboring unit cells

\begin{align} \begin{bmatrix} \mathbf{q}_L^{n+1} \\ \mathbf{f}_L^{n+1} \end{bmatrix} = \mathbf{T} \begin{bmatrix} \mathbf{q}_L^{n} \\ \mathbf{f}_L^{n} \end{bmatrix} = e^{-ikL} \begin{bmatrix} \mathbf{q}_L^{n} \\ \mathbf{f}_L^{n} \end{bmatrix} \end{align}

Computing the eigenvalues \mu_i and eigenvectors \Phi_i, the transfer matrix \mathbf{T} can be diagonalized

\begin{align} \left[\Phi_1, \dots ,\Phi_k\right]^T \mathbf{T} \left[\Phi_1, \dots ,\Phi_k\right] = diag\{\left[\mu_1, \dots ,\mu_k\right]\} \end{align}

The eigenvalues of the transfer matrix \mathbf{T} are the wave solutions that in general can appear in the periodic structure. For each frequency the wave solutions can be computed, determing the eigenvalues \mu=e^{-ikL}.

Visualization

The left plot shows the eigenvalues in the complex plane. To visualize at which frequencies waves can propagate the eigenvalues are visualized with a frequency dependent decay and a phase shift (right plots). The decay and the phase are determined by extracting k_{im} and k_{re} L from the eigenvalues

\begin{align} k_{im}&=\frac{\ln\left(\|\mu\|\right)}{L}\\ k_{re}L&=\tan^{-1}\left(\frac{\mu_{im}}{\mu_{re}}\right) \end{align}