Dispersion curves for multilayer elastic materials using Trasfer Matrix Method and Bloch theorem.

example_layered_elastic.py

# -*- coding: utf-8 -*-
"""
Wave propagation in layered elastic media.

@author: Nicolas Guarin-Zapata
@email: [email protected]
"""
import matplotlib.pyplot as plt
from matplotlib import rcParams
from layered_elastic import *

rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
rcParams['legend.fontsize'] = 15

params = [[70e9, 0.35, 2770],[92e9, 0.33, 8270], [200e9, 0.285, 7800]]
L_list = [1, 1, 1]
nparams = 3

nfre = 10000
omega_vec = np.linspace(1e-6, 25000, nfre)
k_vec = disp_multilay_iso(nparams, params, L_list, omega_vec)
thick = 2
vel = np.sqrt(params[0][0]/params[0][2])

k_vec[np.abs(np.abs(k_vec) - 1) > 3e-1] = np.NaN  # Get rid of the complex freqs

# S-waves
plt.plot(np.abs(np.angle(k_vec[:,0]))/np.pi, omega_vec*thick/vel,'k',lw=2)
plt.plot(np.abs(np.angle(k_vec[:,1]))/np.pi, omega_vec*thick/vel,'b', lw=2)
# P-wave
plt.plot(np.abs(np.angle(k_vec[:,2]))/np.pi, omega_vec*thick/vel,'r--',lw=2)

plt.xlim(0,1)
plt.xlabel(r"$2dk/\pi$", size=20)
loc, labels = plt.xticks()
plt.xticks(loc, size=16)
plt.ylabel(r"$2d\omega\sqrt{\rho/E_1}$", size=20)
loc, labels = plt.yticks()
plt.yticks(loc, size=16)
plt.grid('on')
plt.show()

layered_elastic.py

# -*- coding: utf-8 -*-
"""
Wave propagation in layered elastic media.

@author: Nicolas Guarin-Zapata
@email: [email protected]
"""
from __future__ import division
import numpy as np

def elast_lay_mat(E, nu, rho, omega, L):
    r"""Propagation matrix for a layer of thickness L
    
    Parameters
    ----------
    E : float
        Young modulus.
    nu : float
        Poisson ratio.
    omega : array-like
        Circular Frequency.
    L : float
        Thickness of the layer.
        
    Returns
    -------
    T : (6,6) ndarray
        Propagation matrix.
        
        
    Examples
    --------
    If we propagate the wave a distance `L` and then a distance `-L`
    the product of the two propagation matrices should give the
    identity matrix
    
    >>> T1 = elast_lay_mat(70e9, 0.3, 2700, 100, 1)
    >>> T2 = elast_lay_mat(70e9, 0.3, 2700, 100, -1)
    >>> print np.dot(T1, T2)
    [[ 1.  0.  0.  0.  0.  0.]
     [ 0.  1.  0.  0.  0.  0.]
     [ 0.  0.  1.  0.  0.  0.]
     [ 0.  0.  0.  1.  0.  0.]
     [ 0.  0.  0.  0.  1.  0.]
     [ 0.  0.  0.  0.  0.  1.]]
       
    The square of the propagation matrix for `L` should be the same
    that a propagation for `2L`
      
    >>> T3 = elast_lay_mat(70e9, 0.3, 2700, 100, 2) 
    >>> print np.round(np.dot(T1, T1) - T3, 4)
    [[ 0.  0.  0.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.]
     [ 0.  0.  0.  0.  0.  0.]]
       
    """
   
    alpha = np.sqrt(E*(1 - nu)/((1 + nu)*(1 - 2*nu)*rho))
    beta = np.sqrt(E/(2*(1 + nu)*rho))
    k_P = omega/alpha
    k_S = omega/beta
    
    T = np.zeros((6,6))
    T[0,0] = T[3,3] = np.cos(k_P*L)
    T[1,1] = T[2,2] = np.cos(k_S*L)
    T[4,4] = T[5,5] = np.cos(k_S*L)
    T[0,3] = np.sin(k_P*L)/(k_P*alpha**2*rho)
    T[1,4] = T[2,5] = np.sin(k_S*L)/(k_S*beta**2*rho)
    T[3,0] = -k_P*alpha**2*rho*np.sin(k_P*L)
    T[4,1] = T[5,2] = -k_S*beta**2*rho*np.sin(k_S*L)
    
    return T

def elast_multilay_mat(n, param_list, omega, L_list):
    r"""Propagation matrix for the n-layers material
    
    Parameters
    ----------
    n : int
        Number of layers.
    param_list : (n) list
        List of parameters. Each entry is a list with `[E, nu, rho]`,
        being `E` the Young modulus, `nu` the Poisson coefficient and
        `rho` the mass density.
    omega : float
        Angular frequency.
    L_list : (n) list
        List with the width of each layer.
        
    Returns
    -------
    T : (6,6) ndarray
        Propagation matrix for the multilayer material.
    
    """
    assert len(param_list) == n and len(L_list) == n
    
    T = np.eye(6)
    for k in range(n):
        E = param_list[k][0]
        nu = param_list[k][1]
        rho = param_list[k][2]
        L = L_list[k]
        T_aux = elast_lay_mat(E, nu, rho, omega, L)
        T = np.dot(T, T_aux)
        
    return T

def disp_multilay_iso(nlayers, params, L_list, omega_vec):
    """Compute the dispersion curves for a multilayered material
    
    Parameters
    ----------
    nlayers : int
        Number of layers in the cell.
    params : (nlayers) list
        List of parameters. Each entry is a list with `[E, nu, rho]`,
        being `E` the Young modulus, `nu` the Poisson coefficient and
        `rho` the mass density.
    L_list : list
        List of widths for different layers.
    omega_vec : (n) array
        Array of frequencies.
        
    Returns
    -------
    k_vec : (n) array
        Array of vector numbers.       
        
    
    """
    assert len(params) == nlayers and len(L_list) == nlayers
    
    nfre = len(omega_vec)
    k_vec = np.zeros((nfre,3), dtype=complex)
    for j, omega in enumerate(omega_vec):
        T = elast_multilay_mat(nlayers, params, omega, L_list)
        vals, vecs = np.linalg.eig(T)
        vals, vecs = order_vec(vals, vecs)
        k_vec[j,:] = vals
        
    return k_vec

   
def order_vec(vals, vecs):
    """
    Order the eigenvalues according to their imaginary part and then,
    order the eigenvectors according to their projection respect
    to a canonical coordinate system :math:`[u_x, u_y, u_z]`.
    
    This procedure is done for transverse isotropic media, with
    incident perpendicular incidence. Where we can assure that one of
    the polarizations is longitudinal and the other two are
    transversal. For other incidences (or more general anisotropies)
    this procedure can be cumbersome.[1]
    
    Parameters
    ----------
    vals : (6) ndarray
        Eigenvalues of the propagation matrix.
    vecs : (6,6) ndarray
        Eigenvectors of the propagation matrix.
        
    Returns
    -------
    vals : (3) ndarray
        Ordered eigenvalues of the propagation matrix. 
    vecs : (6,3) ndarray
        Ordered eigenvectors of the propagation matrix.  
    
    References
    ----------
      [1] Carcione, J. Jose M., ed. Wave fields in real media. Elsevier 
      Science, 2007.
      
    """
    
    ivals = np.array([k.imag if (abs(k) - 1.0)<=1e-6 else k.real for k in vals])
    order = np.argsort(ivals) # Order the eigenvalues
    vals = vals[order]
    vecs = vecs[:,order]    
    
    proj1 = abs(np.dot(np.array([0,0,0,1,0,0]),vecs[:,0:3]))
    proj2 = abs(np.dot(np.array([0,0,0,0,1,0]),vecs[:,0:3]))
    proj3 = abs(np.dot(np.array([0,0,0,0,0,1]),vecs[:,0:3]))
    max1 = proj1[0]    
    max2 = proj2[0]
    max3 = proj3[0]
    k1 = 0    
    k2 = 0
    k3 = 0
    for k in range(0,3):
        if max1<proj1[k]:
            max1 = proj1[k]
            k1 = k
        if max2<proj2[k]:
            max2 = proj2[k]          
            k2 = k
        if max3<proj3[k]:
            max3 = proj3[k]
            k3 = k
            
    return vals[[k1,k2,k3]], vecs[:,[k1, k2, k3]]


if __name__=="__main__":
    import doctest
    doctest.testmod()

    import matplotlib.pyplot as plt
    from matplotlib import rcParams
    
    rcParams['font.family'] = 'serif'
    rcParams['font.size'] = 16
    rcParams['legend.fontsize'] = 15
    
    params = [[70e9, 0.35, 2770],[92e9, 0.33, 8270]]
    L_list = [1, 1]
    nparams = 2
    
    nfre = 1000
    omega_vec = np.linspace(1e-6, 15000, nfre)
    k_vec = disp_multilay_iso(nparams, params, L_list, omega_vec)
    thick = 2
    vel = np.sqrt(params[0][0]/params[0][2])
    
    k_vec[np.abs(np.abs(k_vec) - 1) > 3e-1] = np.NaN  # Get rid of the complex freqs
    
    # S-waves
    plt.plot(np.abs(np.angle(k_vec[:,0]))/np.pi, omega_vec*thick/vel,'k',lw=2)
    plt.plot(np.abs(np.angle(k_vec[:,1]))/np.pi, omega_vec*thick/vel,'b', lw=2)
    # P-wave
    plt.plot(np.abs(np.angle(k_vec[:,2]))/np.pi, omega_vec*thick/vel,'r--',lw=2)
    
    plt.xlim(0,1)
    plt.xlabel(r"$2dk/\pi$", size=20)
    loc, labels = plt.xticks()
    plt.xticks(loc, size=16)
    plt.ylabel(r"$2d\omega\sqrt{\rho/E_1}$", size=20)
    loc, labels = plt.yticks()
    plt.yticks(loc, size=16)
    plt.grid('on')
    plt.show()