two_block

# Vibrações Livres de Sistemas Discretos¶

Considerando um protótipo de um sistema de 2 graus de liberdade apresentado na Figura.

As equações de movimento podem ser escritas na forma matricial, conforme segue:

\begin{equation} \begin{bmatrix} m_1 & 0 \\ 0 & m_2 \\ \end{bmatrix} \begin{Bmatrix} \ddot{u}_1 \\ \ddot{u}_2 \end{Bmatrix} + \begin{bmatrix} k_1 + k_2 & -k_2 \\ - k_2 & k_2 + k_3 \\ \end{bmatrix} \begin{Bmatrix} u_1 \\ u_2 \end{Bmatrix} = \begin{Bmatrix} 0 \\ 0 \end{Bmatrix} \end{equation}

Ou ainda, na forma compacta:

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

Utilizando uma abordagem modal, a solução recai em um problema de autovalores. Desta forma, assumindo os seguintes dados:

\begin{align} m_1 = \, 2 \, kg && m_2 = \, 1 \, kg && k_1 &= k_2 = k_3 = \, 1 \, N/m \end{align}

As frequências naturais podem ser obtidas utilizando a rotina desenvolvida anteriormente para a solução de um problema de autovalores, conforme mostrado a seguir.

In [2]:

    %matplotlib inline
    import numpy as np
    from eigen_problem import eigen
    import matplotlib.pyplot as plt
    import matplotlib.animation as animation

    m1 = 2.0      # mass of block 1
    m2 = 1.0      # mass of block 2
    k1 = 1.0      # stiffness of system
    k2 = 1.0      # stiffness of system
    k3 = 1.0      # stiffness of system
    #========================================================================#
    # Natural frequencies and modes                                          #
    #========================================================================#
    print "#"+40*"="+"#\n","#"+" Frequências e modos de vibração"+8*" "+"#\n",\
          "#"+40*"="+"#"
    M = np.array( [[m1, 0.], [0., m2]] )            # inertia matrix
    K = np.array( [[k1+k2, -k2], [-k2, k2+k3]] )    # stiffness matrix
    l, L = eigen(M, K)
    w = np.sqrt(l)
    print "frequências naturais: ", w
    print "modos               : ", L[:,0], L[:,1],"\n"

    #========================================#
    # Frequências e modos de vibração        #
    #========================================#
    frequências naturais:  [ 0.79622522  1.538189  ]
    modos               :  [-0.62796303 -0.45970084] [ 0.32505758 -0.88807383]

Em seguida, considera-se a mudança de coordenadas $\{u\} = [\Gamma]\{ \eta \}$ em que $[\Gamma]$ é a matrix modal normalizada em relação a matriz de massa, $[m]$, e $\{ \eta \}$ é denominado vetor de coordenadas normais. Essa mudança de variáveis, permite desacoplar o sistema de equações de movimento, obtendo:

\begin{equation} \begin{Bmatrix} \ddot{ \eta }_1 \\ \ddot{ \eta }_2 \end{Bmatrix} + \begin{bmatrix} \omega_1^2 & 0 \\ 0 & \omega_2^2 \end{bmatrix} \begin{Bmatrix} \eta_1 \\ \eta_2 \end{Bmatrix} = \begin{Bmatrix} 0 \\ 0 \end{Bmatrix} \end{equation}

ou na forma compacta:

\begin{equation} \ddot{ \eta }_r + \omega_r^2 \eta_r = 0 \,\,\,\,\,\, (r=1, 2) \end{equation}

em que $\omega_1$ e $\omega_2$ são as frequências naturais calculadas anteriormente. Esse procedimento resulta em um sistema de EDOs, cuja solução é discutida em vibrações livres sem amortecimento e tem a forma:

\begin{equation} \eta_r = H_r cos(\omega_r t - \phi_r) \end{equation}

onde as constantes $H_r$ e $\phi_r$ são obtidas aplicando as condições iniciais.

Como tem-se um sistema de duas EDOs de segunda ordem, necessita-se de quatro condições iniciais:

\begin{align} \begin{cases} u_1(0) = u_{1,0} \\ u_2(0) = u_{2,0} \\ \dot{u}_1(0) = v_{1,0} \\ \dot{u}_2(0) = v_{2,0} \end{cases} \end{align}

Essas condições iniciais podem ser transformadas para o domínio das coordenadas naturais, conforme segue:

\begin{align} \{ \eta(0) \} &= [\Gamma]^{-1} \{ u(0) \} = [\Gamma]^{T} [m] \{ u(0) \} \\ \{ \dot{\eta}(0) \} &= [\Gamma]^{-1} \{ \dot{u}(0) \} = [\Gamma]^{T} [m] \{ \dot{u}(0) \} \end{align}

Desta forma, definindo:

\begin{equation} \{ p_0 \} = \begin{Bmatrix} \eta_{1,0} \\ \eta_{2,0} \end{Bmatrix} = \{ \eta (0) \} \end{equation}\begin{equation} \{ q_0 \} = \begin{Bmatrix} \dot{\eta}_{1,0} \\ \dot{\eta}_{2,0} \end{Bmatrix} = \{ \dot{\eta} (0) \} \end{equation}

Então:

\begin{align} p_{r,0} &= H_r cos(\phi_r) \\ q_{r,0} &= -H_r \omega_r sen(\phi_r) \end{align}

Resolvendo para $H_r$ e $\phi_r$, obtém-se:

\begin{align} H_r &= \sqrt{ \bigg( \frac{q_{r,0}}{\omega_r} \bigg)^2 + \big(p_{r,0} \big)^2 } \\ \phi_r &= \arctan \bigg( \frac{q_{r,0}}{\omega_r p_{r,0}} \bigg) \end{align}

Com base na discussão apresentada, vamos implementar uma função denominada coeff_homog em Python para calcular os coeficientes $H_r$ e $\phi_r$ da solução homogênea.

In [3]:

    def coeff_homog(p0, q0, w):
        """
        Function to calculate coefficients (H, phi) of homogeneous solution.
        """
        n = len(w)
        H = np.zeros(n); phi = np.zeros(n)
        for r in range(n):
            H[r] = np.sqrt( (q0[r]/w[r])**2 + (p0[r])**2 )
            phi[r] = np.arctan( q0[r]/(w[r]*p0[r]) )
        return H, phi

Por fim, precisamos escrever uma função solution que receba as constantes $H_r$ e $\phi_r$, converta a solução em coordenadas normais, $\{\eta \}$, para o domínio físico, através da seguinte relação $\{u\} = [\Gamma] \{ \eta \}$.

In [4]:

    def solution(t, H, phi, w, L):
        """
        Function to calculate the solution.
        """
        nr = len(w); nc = len(t)
        p = np.zeros((nr, nc))
        q = np.zeros((nr, nc))
        for r in range(nr):
            p[r,:] = H[r]*np.cos(w[r]*t - phi[r])
            q[r,:] = -H[r]*w[r]*np.sin(w[r]*t - phi[r])

        return np.dot(L, p), np.dot(L, q)

A fim de avaliar a dinâmica do sistema, vamos considerar três casos com condições iniciais diferentes.

## 1) Vibração no 1º modo:¶

Assumindo condições iniciais de velocidade nulas, isto é, $\dot{u}_{1,0}=\dot{u}_{2,0} = 0$ e condições iniciais de deslocamento iguais ao autovetor do primeiro modo, obtém-se os resultados apresentados a seguir.

In [8]:

    u0 = L[:,0]                 # position initial condition (physical domain)
    v0 = np.array([0.0, 0.0])   # velocity initial condition (physical domain)
    #========================================================================#
    # Coordinate transformation                                              #
    #========================================================================#
    LtM = np.dot( L.T, M)       # (L^T)*M
    p0 = np.dot(LtM, u0)        # normal coord., initial position
    q0 = np.dot(LtM, v0)        # normal coord., initial velocity

    print "#"+40*"="+"#\n","#"+" Coeficientes solução homegênea "+8*" "+"#\n",\
          "#"+40*"="+"#"

    H, phi  = coeff_homog(p0, q0, w)

    print "H   = ", H.round(8)
    print "phi = ", phi.round(8)

    t = np.linspace(0., 160., 1000);        dt = t[1] - t[0]
    u, v = solution(t, H, phi, w, L)
    #========================================================================#
    # Plot Results                                                           #
    #========================================================================#
    plt.close('all')
    fig = plt.figure(figsize=(14,5))
    ax = fig.add_subplot(121)
    ax.plot(t, u[0], '-b', lw=1, label=u"M1",   zorder=1)
    ax.plot(t, u[1], '-r', lw=1, label=u"M2",zorder=0)
    ax.set_xlabel(" Tempo (s)", fontsize=12); plt.ylabel(" Deslocamento (m)", fontsize=12)
    plt.legend(ncol=2)
    ymin = min([min(u[0]), min(u[1])])
    ymax = max([max(u[0]), max(u[1])])
    ax.set_ylim(ymin*1.4, ymax*1.4)

    ax = fig.add_subplot(122)
    ax.plot(u[0], v[0], '-b', lw=1, label=u"M1",   zorder=1)
    ax.plot(u[1], v[1], '-r', lw=1, label=u"M2",zorder=0)
    ax.set_xlabel(" Deslocamento (m)", fontsize=12); plt.ylabel(" Velocidade (m/s)", fontsize=12)
    plt.legend(ncol=2)
    plt.show()

    #========================================#
    # Coeficientes solução homegênea         #
    #========================================#
    H   =  [ 1.  0.]
    phi =  [ 0.  0.]

Avaliando os coeficientes da solução homogênea, verifica-se que $H_2 = 0$, ou seja, $\eta_2 = 0$ e o sistema vibra no 1º modo com uma frequência $\omega_1$.

## 2) Vibração no 2º modo:¶

Novamente, adotando condições iniciais de velocidade $\dot{u}_{1,0}=\dot{u}_{2,0} = 0$ e assumindo condições iniciais de deslocamento iguais ao autovetor do segundo modo, obtém-se os resultados apresentados a seguir.

In [7]:

    u0 = L[:,1]                 # position initial condition (physical domain)
    v0 = np.array([0.0, 0.0])   # velocity initial condition (physical domain)
    #========================================================================#
    # Coordinate transformation                                              #
    #========================================================================#
    LtM = np.dot( L.T, M)       # (L^T)*M
    p0 = np.dot(LtM, u0)        # normal coord., initial position
    q0 = np.dot(LtM, v0)        # normal coord., initial velocity

    print "#"+40*"="+"#\n","#"+" Coeficientes solução homegênea "+8*" "+"#\n",\
          "#"+40*"="+"#"

    H, phi  = coeff_homog(p0, q0, w)

    print "H   = ", H.round(8)
    print "phi = ", phi.round(8)

    t = np.linspace(0., 160., 1000);        dt = t[1] - t[0]
    u, v = solution(t, H, phi, w, L)
    #========================================================================#
    # Plot Results                                                           #
    #========================================================================#
    plt.close('all')
    fig = plt.figure(figsize=(14,5))
    ax = fig.add_subplot(121)
    ax.plot(t, u[0], '-b', lw=1, label=u"M1",   zorder=1)
    ax.plot(t, u[1], '-r', lw=1, label=u"M2",zorder=0)
    ax.set_xlabel(" Tempo (s)", fontsize=12); plt.ylabel(" Deslocamento (m)", fontsize=12)
    plt.legend(ncol=2)
    ymin = min([min(u[0]), min(u[1])])
    ymax = max([max(u[0]), max(u[1])])
    ax.set_ylim(ymin*1.4, ymax*1.4)

    ax = fig.add_subplot(122)
    ax.plot(u[0], v[0], '-b', lw=1, label=u"M1",   zorder=1)
    ax.plot(u[1], v[1], '-r', lw=1, label=u"M2",zorder=0)
    ax.set_xlabel(" Deslocamento (m)", fontsize=12); plt.ylabel(" Velocidade (m/s)", fontsize=12)
    plt.legend(ncol=2)
    plt.show()

    #========================================#
    # Coeficientes solução homegênea         #
    #========================================#
    H   =  [ 0.  1.]
    phi =  [ 0.  0.]

Nesse caso, uma avaliação dos coeficientes da solução homogênea, motra que $H_1 = 0$, ou seja, $\eta_1 = 0$ e o sistema vibra no 2º modo com uma frequência $\omega_2$.

## 3) Vibração com composição dos dois modos:¶

Novamente, adotando condições iniciais de velocidade $\dot{u}_{1,0}=\dot{u}_{2,0} = 0$ e assumindo condições iniciais de deslocamento $u_{1,0} = 1.0$ e $u_{2,0} = -1.2$, obtém-se os resultados apresentados a seguir.

In [9]:

    u0 = np.array([0.1,-0.3])   # position initial condition (physical domain)
    v0 = np.array([0.1, 0.2])   # velocity initial condition (physical domain)
    #========================================================================#
    # Coordinate transformation                                              #
    #========================================================================#
    LtM = np.dot( L.T, M)       # (L^T)*M
    p0 = np.dot(LtM, u0)        # normal coord., initial position
    q0 = np.dot(LtM, v0)        # normal coord., initial velocity

    print "#"+40*"="+"#\n","#"+" Coeficientes solução homegênea "+8*" "+"#\n",\
          "#"+40*"="+"#"

    H, phi  = coeff_homog(p0, q0, w)

    print "H   = ", H.round(8)
    print "phi = ", phi.round(8)

    t = np.linspace(0., 160., 1000);        dt = t[1] - t[0]
    u, v = solution(t, H, phi, w, L)
    #========================================================================#
    # Plot Results                                                           #
    #========================================================================#
    plt.close('all')
    fig = plt.figure(figsize=(14,5))
    ax = fig.add_subplot(121)
    ax.plot(t, u[0], '-b', lw=1, label=u"M1",   zorder=1)
    ax.plot(t, u[1], '-r', lw=1, label=u"M2",zorder=0)
    ax.set_xlabel(" Tempo (s)", fontsize=12); plt.ylabel(" Deslocamento (m)", fontsize=12)
    plt.legend(ncol=2)
    ymin = min([min(u[0]), min(u[1])])
    ymax = max([max(u[0]), max(u[1])])
    ax.set_ylim(ymin*1.4, ymax*1.4)

    ax = fig.add_subplot(122)
    ax.plot(u[0], v[0], '-b', lw=1, label=u"M1",   zorder=1)
    ax.plot(u[1], v[1], '-r', lw=1, label=u"M2",zorder=0)
    ax.set_xlabel(" Deslocamento (m)", fontsize=12); plt.ylabel(" Velocidade (m/s)", fontsize=12)
    plt.legend(ncol=2)
    plt.show()

    #========================================#
    # Coeficientes solução homegênea         #
    #========================================#
    H   =  [ 0.27348261  0.33942195]
    phi =  [-1.52574112 -0.21738381]

Note que os coeficientes da solução homogênea são ambos diferentes de zero, com isso, a resposta do sistema é uma composição dos dois modos de vibrar.