vibration_absorber

# Absorvedor Dinâmico de Vibrações¶

Uma forma eficiente de promover o controle passivo de um sistema mecânico primário, consiste no acoplamento de um sistema secundário com o objetivo de absorver energia do sistema primário visando diminuir as amplitudes de vibração.

A Figura mostra um esquemático do sistema primário acoplado com o sistema secundário. Esse sistema pode ser modelado como um sistema de dois graus de liberdade, em que as equações de movimento podem ser escritas na forma matricial, conforme segue:

\begin{equation} \begin{bmatrix} m & 0 \\ 0 & m_a \\ \end{bmatrix} \begin{Bmatrix} \ddot{u} \\ \ddot{u}_a \end{Bmatrix} + \begin{bmatrix} k + k_a & -k_a \\ - k_a & k_a \\ \end{bmatrix} \begin{Bmatrix} u \\ u_a \end{Bmatrix} = \begin{Bmatrix} F_0 sen( \Omega t) \\ 0 \end{Bmatrix} \end{equation}

Considerando o sistema em regime permanente, as respostas $u(t)$ e $u_a(t)$ são dada por:

\begin{equation} u = U sen(\Omega t) \\ u_a = U sen(\Omega t) \end{equation}

Substituindo essas soluções nas equações de movimento, obtém-se:

\begin{equation} \begin{bmatrix} k + k_a - m \Omega^2 & -k_a \\ - k_a & k_a -m_a \Omega^2 \\ \end{bmatrix} \begin{Bmatrix} U \\ U_a \end{Bmatrix} = \begin{Bmatrix} F_0 \\ 0 \end{Bmatrix} \end{equation}

ou, na forma compacta:

\begin{equation} ([k] - \Omega^2 [m])\{ U \} = \{f(t)\} \end{equation}

As amplitudes do movimento, podem ser obtidas pré-mulitplicando a equação anterior pela inversa de $([k] - \Omega^2 [m])$, resultando em:

\begin{equation} \begin{Bmatrix} U \\ U_a \end{Bmatrix} = \frac{1}{(k+k_a-m\Omega^2)(k_a-m_a\Omega^2)-k_a^2} \begin{bmatrix} k_a - m_a \Omega^2 && k_a \\ k_a && k + k_a -m \Omega^2 \end{bmatrix} \begin{Bmatrix} F_0 \\ 0 \end{Bmatrix} \end{equation}

Desta forma, a amplitude do movimento do sistema primário é dada por:

\begin{equation} U = \frac{(k_a - m_a \Omega^2) F_0} {(k+k_a - m\Omega^2)(k_a - m_a \Omega^2) -k_a^2} \end{equation}

Enquanto a amplitude do movimento do sistema absorvedor é:

\begin{equation} U_a = \frac{ k_a F_0} {(k+k_a - m\Omega^2)(k_a - m_a \Omega^2) -k_a^2} \end{equation}

Definindo as seguintes variáveis:

\begin{align} \omega_p &= \sqrt{\frac{k}{m}} \rightarrow \text{ Frequência natural do sistema primário} \\ \\ \omega_a &= \sqrt{\frac{k_a}{m_a}} \rightarrow \text{ Frequência natural do absorvedor isolado} \\ \\ \mu &= \frac{m_a}{m} \rightarrow \text{ Razão entre a massa do absorvedor e a massa do sistema primário}\\ \\ u_E &= \frac{F_0}{k} \rightarrow \text{ Deslocamento estático do sistema primário} \\ \end{align}

As amplitudes do movimento podem ser reescritas da seguinte forma:

\begin{equation} U = \frac{[1 - (\Omega/\omega_a)^2] u_E} {[1 + \mu (\omega_a/\omega_p)^2 - (\Omega/\omega_p)^2][1 - (\Omega/\omega_a)^2] -\mu (\omega_a/\omega_p)^2} \\ \,\,\,\,\, \\ U_a = \frac{u_E} {[1 + \mu (\omega_a/\omega_p)^2 - (\Omega/\omega_p)^2][1 - (\Omega/\omega_a)^2] -\mu (\omega_a/\omega_p)^2} \end{equation}

Vamos visualizar gráficamente a amplitude do sistema primário normalizada pelo deslocamento estático, $|U/u_E|$.

In [2]:

    %matplotlib inline
    from numpy import linspace,sqrt,logical_and  # import functions from library
    import matplotlib.pyplot as plt              # plotting library

    def amplification_factor_1(w, wp):
        """ Amplification factor of the primary system (1 gdl) """
        r = w/wp
        return abs(1./(1. - r**2))

    def amplification_factor_2(w, wa, wp, mi):
        """ Amplification factor of system with absorber (2 gdl) """
        w1 = w/wa
        w2 = wa/wp
        w3 = w/wp
        return abs((1.-(w1)**2)/((1.+mi*(w2**2)-(w3**2))*(1.-(w1**2))-mi*(w2**2)))
    #========================================================================#
    # Imputs                                                                 #
    #========================================================================#
    mi = 0.25                           # mass ratio
    wa = 1.0                            # absorber frequency
    wp = 1.0                            # primary system
    w = linspace(0., 2., 1000)          # frequency data
    G1 = amplification_factor_1(w, wp)
    G2 = amplification_factor_2(w, wa, wp, mi)
    #========================================================================#
    # Plot figure                                                            #
    #========================================================================#
    plt.close('all')
    plt.figure(figsize=(8,5))
    plt.plot(w, G1, '-r', lw=1.5, label='1 gdl')
    plt.plot(w, G2, '-b', lw=1.5, label='2 gdl')
    plt.xlabel('$\\Omega / \omega_{a}$', fontsize=14)
    plt.ylabel('$|(U/u_{E})|$', fontsize=14)
    plt.ylim(0., 5.)
    plt.legend()
    plt.show()

É possível observar que existe uma faixa de frequências de funcionamento em que o absorvedor apresenta um desempenho adequado, que pode ser definida quando $|U/u_E|$ é menor que essa quantidade para o sistema de 1 gdl isolado. Matematicamente, essa região pode ser obtida fazendo:

\begin{equation} \bigg|\frac{1 - (\Omega/\omega_a)^2}{[1+\mu(\omega_a/\omega_p)^2-(\Omega/\omega_p)^2][1-(\Omega/\omega_a)^2]-\mu(\omega_a/\omega_p)^2}\bigg| = \bigg|\frac{1}{1 - (\Omega/\omega_p)} \bigg| \end{equation}

Resolvendo, obtém-se:

\begin{equation} r_{3,4}^2 = \bigg( \frac{\Omega}{\omega_a} \bigg)^2 = \frac{3+\beta^2(1+\mu) \pm \sqrt{[3+\beta^2(1+\mu)]^2 - 8 (\beta^2+1)}}{2 (\beta^2 +1)} \end{equation}

Um ponto importante é que o sistema primário isolado possui apenas 1 gdl e, portanto, apenas uma frequência de ressonância. Ao adicionar o absorvedor, introduz-se um grau de liberdade adicional, fazendo com que o novo sistema tenha duas frequências de ressonância, uma maior e outra menor que a frequência de ressoância do sistema primário isolado. As frequências naturais do sistema de 2 gdl podem ser obtidas fazendo $det([k] - \Omega^2 [m]) = 0$. Assim, interpretando $\Omega$ como $\omega_n$, obtém-se a seguinte equação característica:

\begin{equation} \beta^2 \bigg(\frac{\omega_n^2}{\omega_a^2} \bigg)^2 - [1 + \beta^2(1 + \mu)] \bigg(\frac{\omega_n^2}{\omega_a^2} \bigg) + 1 = 0 \end{equation}

em que $\beta = \omega_a/\omega_p$. Desta forma, resolvendo o polinômio característico, tem-se:

\begin{equation} r_{1,2}^2 = \bigg( \frac{\omega_n}{\omega_a} \bigg)^2 = \frac{1+\beta^2(1+\mu)}{2 \beta^2} \pm \frac{1}{2 \beta^2} \sqrt{[1 + \beta^2 (a+\mu)]^2 - 4 \beta^2} \end{equation}

em que mostra como as frequências naturais variam com a razão de massa, $\mu$, e a razão de frequências, $\beta$. A função parameters calcula as frequências naturais descritas anteriormente, assim como a região em que o absorvedor tem um desempenho adequado.

In [3]:

    def parameters(mi, beta):
        """
        r1 : first dimensionless ressonance frequency of 2 gdl system
        r2 : second dimensionless ressonance frequency of 2 gdl system
        r3 <= r <= r4 : range of frequency to passive absorver
        y3 = f(r3) : dimensionless amplitude at r3
        y4 = f(r4) : dimensionless amplitude at r4
        """
        r1 = sqrt(((1.+(beta**2)*(1.+mi)) - sqrt((1.+(beta**2)*(1.+mi))**2 -     \
                   4.*beta**2))/(2.*beta**2))
        r2 = sqrt(((1.+(beta**2)*(1.+mi)) + sqrt((1.+(beta**2)*(1.+mi))**2 -     \
                   4.*beta**2))/(2.*beta**2))
        r3 = sqrt(((3.+(beta**2)*(1.+mi)) - sqrt((3.+(beta**2)*(1.+mi))**2 -     \
                   8.*(beta**2 +1.)))/(2.*(beta**2 +1.)))
        y3 = amplification_factor_1(r3, 1.)
        r4 = sqrt(((3.+(beta**2)*(1.+mi)) + sqrt((3.+(beta**2)*(1.+mi))**2 -     \
                   8.*(beta**2 +1.)))/(2.*(beta**2 +1.)))
        y4 = amplification_factor_1(r4, 1.)
        return r1,r2,r3,y3,r4,y4

Vamos visualizar graficamente a região de frequência de funcionamento em que o absorvedor tem desempenho adequado, e as frequências de ressonância do sistema de 2 gdl.

In [4]:

    #========================================================================#
    # Imputs                                                                 #
    #========================================================================#
    mi = 0.25                           # mass ratio
    wa = 1.0                            # absorber frequency
    wp = 1.0                            # primary system
    w = linspace(0., 2., 1000)          # frequency data
    y = linspace(0., 1., 1000)          # dummy vector
    G1 = amplification_factor_1(w, wp)
    G2 = amplification_factor_2(w, wa, wp, mi)
    r1,r2,r3,y3,r4,y4 = parameters(mi,wa/wp)
    #========================================================================#
    # Plot figure                                                            #
    #========================================================================#
    plt.figure(figsize=(8,5))
    plt.plot(w/wa, G1, '-r', lw=1.5, label='1 gdl')
    plt.plot(w/wa, G2, '-b', lw=1.5, label='2 gdl')
    plt.plot(0*w/wa + r1, 4*w, '--k')
    plt.plot(0*w/wa + r2, 4*w, '--k')
    plt.plot(0*w/wa + r3, y*y3, '-k')
    plt.plot(0*w/wa + r4, y*y4, '-k')
    plt.fill_between(w/wa,0, G2, where=logical_and(r3<=w/wa,w/wa<=r4),
                     interpolate=True, alpha=0.3)
    plt.xlabel('$\\Omega / \omega_{a}$', fontsize=14)
    plt.ylabel('$|(U/u_{E})|$', fontsize=14)
    plt.legend()
    plt.text(r1*0.84, 0.1, '%.2f'%r1)
    plt.text(r2*1.02, 0.1, '%.2f'%r2)
    x_ticks=[0.,0.5, r3, 1.0, r4, 1.5, 2.]
    x_ticks = [round(elem, 2) for elem in x_ticks]
    plt.xticks(x_ticks)
    plt.ylim(0., 5.)
    plt.show()

## Simulações Numéricas¶

Em seguida, vamos analisar o comportamento do absorvedor dinâmico projetado anteriormente, no domínio do tempo.

### 1) Sistema sem amortecimento¶

Conforme visto, as equações de movimento do sistema primário acoplado com o absorvedor são dadas por:

\begin{equation} \begin{cases} m \ddot{u} + (k + k_a) u - k_a u_a = F_0 sen(\Omega t) \\ m_a \ddot{u}_a - k_a u + k_a u_a = 0 \end{cases} \end{equation}

Na função func1 é descrito as equações de movimento na sua forma canônica.

In [5]:

    #============================================================================#
    # Importing libraries                                                        #
    #============================================================================#
    %matplotlib inline
    import numpy as np
    import matplotlib.pyplot as plt
    import matplotlib.animation as animation
    from rk4 import integrate
    #============================================================================#
    # Right-hand side of differential equations                                  #
    #============================================================================#
    def func1(x, t):
        ff = np.zeros_like(x)
        ff[0] = x[1]
        ff[1] = (f0*np.sin(omega*t)+Ka*x[2]-(K+Ka)*x[0])/(M)
        ff[2] = x[3]
        ff[3] = (Ka*x[0]-Ka*x[2])/(Ma)
        return ff

A função ShowResults, a seguir, é utilizada para plotar os deslocamentos do sistema primário e do absorvedor, assim como os respectivos espaços de estado.

In [6]:

    def ShowResults(matrix):
        #========================================================================#
        # Plot Displacements                                                     #
        #========================================================================#
        t = matrix[:,0]          # time data
        fig = plt.figure(figsize=(14,5))
        ax = fig.add_subplot(121)
        ax.plot(t, matrix[:,1],'-b', lw=1.5, label='$M_1$', zorder=20)
        ax.plot(t, matrix[:,3],'-r' ,lw=1.5, label='$M_2$', zorder=18)
        ax.set_xlabel(' $ t $ ', fontsize=14)
        ax.set_ylabel(' $ u(t) $ ', fontsize=14)
        ax.set_xlim(0., t[-1])
        ax.set_ylim(-1.1, 1.1)
        plt.legend(ncol=2)
        #========================================================================#
        # State space                                                            #
        #========================================================================#
        ax = fig.add_subplot(122)
        ax.plot(matrix[:,1],matrix[:,2],'-b', lw=1.5, label='$M_1$')
        ax.plot(matrix[:,3],matrix[:,4],'-r', lw=1.5, label='$M_2$')
        ax.set_xlabel('$ u(t) $', fontsize=14)
        ax.set_ylabel('$ \dot{u}(t) $', fontsize=14)
        ax.set_ylim(-2.5, 2.5)
        plt.legend(ncol=2)
        plt.show()

Em seguida, definimos os parâmetros utilizados nas simulações numéricas. Lembrando que o absorvedor dinâmico tem um funcionamento adequado quando $\omega_a \approx \omega_p$.

In [7]:

    #========================================================================#
    # Inputs                                                                 #
    #========================================================================#
    M  = 100.0              # mass of block
    Ma = 4.0                # mass of absorber
    K  = 600.0              # stiffness
    Ka = 25.0               # stiffness of absorber system
    f0 = 10.                # excitation amplitude
    omega = np.sqrt(K/M)    # excitation frequency

    t0 = 0.                 # initial time
    dt = 0.05               # time step
    N  = 3200               # number of time steps

    u1 = 0.0                # initial displacement of block
    v1 = 0.0                # initial velocity of block
    u2 = 0.0                # initial displacement of absorber
    v2 = 0.0                # initial velocity of absorber
    init_conditions = np.array([u1, v1, u2, v2])
    #========================================================================#
    # Numerical Solution using RK4                                           #
    #========================================================================#
    matrix = integrate(t0, dt, N, init_conditions, func1)# matrix with results
    ShowResults(matrix)

Vale a pena ressaltar que a análise realizada no domínio da frequência assume a condição de regime permanente. Nesta simulação, não é esse o caso, uma vez que oscilações transientes ocorrem quando o sistema sai do repouso para o estado forçado de acordo com $f(t)$. Como nesse caso, o sistema não tem amortecimento, a resposta transiente nunca decai. No entanto, note que a amplitude de $u(t)$ é limitida. Isto não ocorreria sem o absorvedor, uma vez que o sistema está sendo excitado na frequência de ressonância do sistema primário isolado e, consequentemente, as amplitudes de resposta tenderiam ao infinito.

### 2) Sistema com amortecimento viscoso¶

Neste momento vamos acrescentar um amortecimento viscoso, $c$, em paralelo com a mola de rigidez $k$, conforme mostrado na Figura. Nesse caso, as equações de movimento podem ser escritas na forma matricial, conforme segue:

\begin{equation} \begin{bmatrix} m & 0 \\ 0 & m_a \\ \end{bmatrix} \begin{Bmatrix} \ddot{u} \\ \ddot{u}_a \end{Bmatrix} + \begin{bmatrix} c & 0 \\ 0 & 0 \\ \end{bmatrix} \begin{Bmatrix} \dot{u} \\ \dot{u}_a \end{Bmatrix} + \begin{bmatrix} k + k_a & -k_a \\ - k_a & k_a \\ \end{bmatrix} \begin{Bmatrix} u \\ u_a \end{Bmatrix} = \begin{Bmatrix} F_0 sen(\Omega t) \\ 0 \end{Bmatrix} \end{equation}

ou ainda:

\begin{equation} \begin{cases} m \ddot{u} + c \dot{u}+ (k + k_a) u - k_a u_a = F_0 sen(\Omega t) \\ m_a \ddot{u}_a - k_a u + k_a u_a = 0 \end{cases} \end{equation}

Na função func2 é descrito as equações de movimento na sua forma canônica.

In [8]:

    #============================================================================#
    # Right-hand side of differential equations                                  #
    #============================================================================#
    def func2(x, t):
        ff = np.zeros_like(x)
        ff[0] = x[1]
        ff[1] = (f0*np.sin(omega*t)+Ka*x[2]-(K+Ka)*x[0]- C*x[1])/(M)
        ff[2] = x[3]
        ff[3] = (Ka*x[0]-Ka*x[2])/(Ma)
        return ff

Inicialmente, vamos analisar o sistema primário sem o absorvedor acoplado. Isto pode ser feito atribuido $k_a = 0$. Os resultados obtidos são apresentados a seguir, onde é possivel observar que, apesar do sistema está sendo excitado na frequência de ressonância, o amortecimento limita as amplitudes em regime permanente. Note que o absorvedor tem deslocamento nulo, conforme desejado.

In [9]:

    #========================================================================#
    # Inputs                                                                 #
    #========================================================================#
    M  = 100.0              # mass of block
    Ma = 4.0                # mass of absorber
    K  = 600.0              # stiffness
    Ka = 0.0                # stiffness of absorber system
    C  = 10.0               # damping
    f0 = 10.                # excitation amplitude
    omega = np.sqrt(K/M)    # excitation frequency

    t0 = 0.                 # initial time
    dt = 0.05               # time step
    N  = 3200               # number of time steps

    u1 = 0.0                # initial displacement of block
    v1 = 0.0                # initial velocity of block
    u2 = 0.0                # initial displacement of absorber
    v2 = 0.0                # initial velocity of absorber
    init_conditions = np.array([u1, v1, u2, v2])
    #========================================================================#
    # Numerical Solution using RK4                                           #
    #========================================================================#
    matrix = integrate(t0, dt, N, init_conditions, func2)# matrix with results
    ShowResults(matrix)

Agora, vamos acoplar o absorvedor dinâmico ao sistema primário, fazendo novamente $k_a = 25$.

In [10]:

    #========================================================================#
    # Inputs                                                                 #
    #========================================================================#
    M  = 100.0              # mass of block
    Ma = 4.0                # mass of absorber
    K  = 600.0              # stiffness
    Ka = 25.0               # stiffness of absorber system
    C  = 10.0               # damping
    f0 = 10.                # excitation amplitude
    omega = np.sqrt(K/M)    # excitation frequency

    t0 = 0.                 # initial time
    dt = 0.05               # time step
    N  = 3200               # number of time steps

    u1 = 0.0                # initial displacement of block
    v1 = 0.0                # initial velocity of block
    u2 = 0.0                # initial displacement of absorber
    v2 = 0.0                # initial velocity of absorber
    init_conditions = np.array([u1, v1, u2, v2])
    #========================================================================#
    # Numerical Solution using RK4                                           #
    #========================================================================#
    matrix = integrate(t0, dt, N, init_conditions, func2)# matrix with results
    ShowResults(matrix)

É possivel visualizar que o absorvedor dinâmico faz o papel para qual foi projeto, retirando energia do sistema primário, limitando as amplitudes de resposta. Note que, devido ao amortecimento introduzido, as oscilações transientes decaem ao longo do tempo até atingir regime permanente.