eigen_problem

# Abordagem Modal¶

Considerando o problema de vibrações livres de um sistema não dissipativo. O sistema de equações deferenciais de movimento tem a seguinte forma:

\begin{equation} [m]\{\ddot{u}\} + [k] \{u\} = \{0\} \end{equation}

Uma abordagem de solução para esse sistema consiste em considerar uma separação de variáveis do tipo:

\begin{equation} \{u\} = \{U\} g \end{equation}

em que o vetor $\{U\}$ representa a característica espacial do problema, enquanto $g=g(t)$ é uma função do tempo e igual para todas componentes de $\{u\}$. Substituindo essa relação na equação de movimento, obtém-se:

\begin{equation} [m]\{U\}\ddot{g} + [k]\{U\}g = \{0\} \end{equation}

Multiplicando ambos os lado da equação anterior por $\{U\}^T$, podemos escrever o sistema na seguinte forma:

\begin{equation} \{U\}^T [m] \{U\} \ddot{g} + \{U\}^T [k] \{U\} g = 0 \end{equation}

Essa equação pode ser reescrita como segue:

\begin{equation} \frac{ \{U\}^T [k] \{U\} }{ \{U\}^T [m] \{U\} } = -\frac{\ddot{g}}{g} \end{equation}

Avaliando a equação, nota-se que o lado direito depende do tempo, enquanto o lado, por sua vez, é uma característica espacial que não depende do tempo. Desta forma, essa expressão deve ser uma constante, digamos $\omega^2$. Assim, o problema pode ser descrito conforme segue:

\begin{equation} \begin{cases} \ddot{g} + \omega^2 g = 0 \\ \big( [k] - \omega^2 [m] \big) \{U\} = \{0\} \end{cases} \end{equation}

A primeira equação é uma EDO no tempo e corresponde a um problema de vibrações livre não dissipativo visto anteriormente. A segunda equação é espacial definindo um problema de autovalores e autovetores.

Como o sistema é homogêneo, sempre possui uma solução trivial, $\{U\} = \{0\}$. Entretanto, buscando soluções diferentes da trivial, necessita-se que o sistema tenha infinitas soluções, isto é imposto da seguinte forma:

\begin{equation} det \big( [k] - \omega^2 [m] \big) = 0 \end{equation}

Essa condição está associada ao polinomio característico do sistema. Portanto, existem $r$ autovalores $\lambda_r$ ou $r$ raízes do polinômico característico, tal que $\lambda_r = \omega_r^2$. Desta forma, os autovalores estão associados as frequências naturais do sistema, $\omega_r$.

Cada autovalor possui um autovetor associado, $\{U^{(r)}\}$, denominado de modo natural. Cada autovetor $\{U^{(r)}\}$, pode ser obtido resolvendo:

\begin{equation} \big( [k] - \omega_r^2 [m] \big) \{U^{(r)}\} = \{0\} \end{equation}

Considerando ainda, uma normalização em relação a matriz de massa $[m]$:

\begin{equation} \{ U^{(r)} \}^T [m] \{ U^{(r)} \} = 1 \end{equation}

Tem-se, a partir do problema de autovalores:

\begin{equation} [k] \{ U^{(r)} \} = \omega_r^2 [m] \{ U^{(r)} \} \end{equation}

em que, multiplicando por $\{ U^{(r)} \}^T$, resulta:

\begin{equation} \{ U^{(r)} \} [k] \{ U^{(r)} \} = \omega_r^2 \{ U^{(r)} \} [m] \{ U^{(r)} \} \end{equation}

Portanto, da normalização, obtém-se que:

\begin{equation} \{ U^{(r)} \} [k] \{ U^{(r)} \} = \omega_r^2 \end{equation}

## Problema de autovalores¶

A solução do problema de autovalores possui uma forma-padrão do tipo:

\begin{equation} \big( [A] - \lambda [I] \big) \{U\} = \{0\} \end{equation}

onde $[I]$ é a matriz identidade. O uso dessa forma-padrão é interessante, pois existem muitos algoritmos disponíveis para solucionar esse problema específico. Faremos uso de um algoritmo implementado em Python. Para isso, vamos converter o problema de autovalores em questão para a forma usual. Isto pode ser feito considerando a seguinte mudança de coordenadas:

\begin{equation} \{ U \} = [m]^{-1} \{q\} \end{equation}

Usando essa transformação no problema de autovalores, tem-se que:

\begin{equation} \big( [k] [m]^{-1} - \lambda [m] [m]^{-1} \big) \{q\} = \{0\} \end{equation}

Sabendo que:

\begin{equation} [m] [m]^{-1} = [I] \end{equation}

E definindo:

\begin{equation} [k] [m]^{-1} = [\tilde{k}] \end{equation}

Resulta:

\begin{equation} \big( [\tilde{k}] - \lambda [I] \big) \{q\} = \{0\} \end{equation}

onde $\lambda = \omega^2$.

Segue uma implementação em Python de uma função para resolver o problema de autovalores.

In [1]:

    #============================================================================#
    #                       Eigenvalues and Eigenvectors                         #
    #============================================================================#
    import numpy as np
    def eigen(M, K):
        """
        Function that receives [M] and [K] matrices, return eigenvalues and
        normalized eigenvectors with respect to [M].

        Parameters:
        -----------

        M : ndarray
            inertia matrix (n x n)
        K : ndarray
            stiffness matrix (n x n)
        Mi : ndarray
            inverse matrix of [M]
        Kt : ndarray
            computes the following matrix operation: [K]*[Mi]

        Return:
        -------
        l : ndarray
            eigenvalues
        Ut : ndarray
            dimensionless eigenvectors as: ( [Kt] - l*[I] )*{Ut} = {0}

        Reference:
        G. Strang, Linear Algebra and Its Applications, 2nd Ed., Orlando, FL,
        Academic Press, Inc., 1980, Various pp.
        """
        if (len(M) != len(K)):
            raise ValueError("[M] and [K] matrices must have same size.")
        else:
            Mi = np.linalg.inv(M)          # compute the inverse of matrix [M].
            Kt = np.dot(K, Mi)             # performe matrix multiplication.
            l, Ut = np.linalg.eig(Kt)      # solve the eigenvalue problem.
            U = np.dot(Mi,Ut)
            #====================================================================#
            # Normalization with respect to [M]                                  #
            #====================================================================#
            for i in range(len(l)):
                alpha = np.dot(U[:,i].T,M)
                alpha = np.dot(alpha, U[:,i])
                U[:,i] = U[:,i]/np.sqrt(alpha)
            return l, U

In [6]:

    #============================================================================#
    #                               MAIN PROGRAM                                 #
    #============================================================================#
    if __name__=='__main__':

        print "#================================================================#"
        print "# Exemplo 1                                                      #"
        print "#================================================================#"
        M = np.array([[9., 0.],[0.,1.]])        # inertia matrix
        K = np.array([[27.,-3.],[-3.,3.]])      # stiffness matrix
        #========================================================================#
        # Eigenvalue and eigenvectors                                            #
        #========================================================================#
        l, U = eigen(M,K)
        w = np.sqrt(l)
        #========================================================================#
        # Print result on screen                                                 #
        #========================================================================#
        for i in range(len(w),0,-1):
            print "Frequência Natural %i : "%i, w[i-1]
            print "Modo Natural       %i : "%i, U[:,i-1],"\n"
        print "#================================================================#"
        print "# Exemplo 2                                                      #"
        print "#================================================================#"
        M = np.array([[1., 0.,0.],[0.,1.,0.],[0.,0.,1.]])       # inertia matrix
        K = np.array([[2.,-1.,0.],[-1.,2.,-1],[0.,-1.,1.]])     # stiffness matrix
        #========================================================================#
        # Eigenvalue and eigenvectors                                            #
        #========================================================================#
        l, U = eigen(M,K)
        w = np.sqrt(l)
        #========================================================================#
        # Print result on screen                                                 #
        #========================================================================#
        for i in range(len(w),0,-1):
            print "Frequência Natural %i : "%i, w[i-1]
            print "Modo Natural       %i : "%i, U[:,i-1],"\n"
    #============================================================================#

    #================================================================#
    # Exemplo 1                                                      #
    #================================================================#
    Frequência Natural 2 :  1.41421356237
    Modo Natural       2 :  [ 0.23570226  0.70710678]

    Frequência Natural 1 :  2.0
    Modo Natural       1 :  [ 0.23570226 -0.70710678]

    #================================================================#
    # Exemplo 2                                                      #
    #================================================================#
    Frequência Natural 3 :  0.445041867913
    Modo Natural       3 :  [ 0.32798528  0.59100905  0.73697623]

    Frequência Natural 2 :  1.24697960372
    Modo Natural       2 :  [-0.73697623 -0.32798528  0.59100905]

    Frequência Natural 1 :  1.8019377358
    Modo Natural       1 :  [-0.59100905  0.73697623 -0.32798528]