<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Análise de um oscilador com N graus de liberdade sujeito a uma excitação dinâmica aplicada nalguns graus de liberdade
Formulação do problema
Equação de equilíbrio dinâmico em coordenadas reais de um oscilador com N graus de liberdade sujeito a uma excitação dinâmica
Step1:  Dados
Oscilador com três graus de liberdade
Step2:  Análise modal
Step3:  Ordenação dos modos
Step4:  Normalização dos modos de vibração relativamente à matriz de massa
Step5:  Massa e rigidez modais
Step6:  Excitação dinâmica
Step8:  Função auxiliar para o cálculo da resposta modal
Step9:  Solução clássica
Step10:  Solução alternativa
<ASSISTANT_TASK:>
Python Code:
import sys
import math
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
%matplotlib inline
print('System: {}'.format(sys.version))
for package in (np, mpl):
    print('Package: {} {}'.format(package.__name__, package.__version__))
Explanation: Análise de um oscilador com N graus de liberdade sujeito a uma excitação dinâmica aplicada nalguns graus de liberdade
Formulação do problema
Equação de equilíbrio dinâmico em coordenadas reais de um oscilador com N graus de liberdade sujeito a uma excitação dinâmica:
$$[M] \times {\ddot u} + [C] \times {\dot u} + [K] \times {u} = {s} \cdot p(t)$$
O vetor ${s}$ contém a descrição espacial da excitação cuja variação ao longo do tempo é definida através do termo $p(t)$. Admite-se no presente caso, sem perda de generalidade, que a excitação é uma função harmónica ($p(t) = \sin (\omega \cdot t)$).
Adotando a seguinte transformação de coordenadas reais em coordenadas modais:
$${u} = [\Phi] \times {q} = \sum_{n=1}^{N} {\phi_n} \cdot q_n$$
é possível converter as N equações de equilíbrio dinâmico para coordenadas modais. Para o modo n obtém-se:
$$M_n \cdot {\ddot q_n} + C_n \cdot {\dot q_n} + K_n \cdot {q_n} = {\phi_n}^T \times {s} \cdot p(t)$$
em que:
$M_n = {\phi_n}^T \times [M] \times {\phi_n}$ - Massa modal
$C_n = {\phi_n}^T \times [C] \times {\phi_n} = 2 \cdot \zeta_n \cdot M_n \cdot \omega_n$ - Amortecimento modal clássico
$K_n = {\phi_n}^T \times [K] \times {\phi_n}$ - Rigidez modal
$\zeta_n$ é o coeficiente de amortecimento modal e $\omega_n$ é a frequência de vibração modal.
Solução claśsica
Expansão modal do vetor de excitação:
$${s} = \sum_{n=1}^{N} {s_n} = \sum_{n=1}^{N} \Gamma_n \cdot [M] \times {\phi_n}$$
em que $\Gamma_n$ é o fator de participação do modo n dado por:
$$\Gamma_n = \frac{{\phi_n}^T \times {s}}{M_n}$$
A equação de equilíbrio dinâmico em coordenadas modais passa a escrever-se:
$$\ddot q_n + 2 \cdot \zeta_n \cdot \omega_n \cdot \dot q_n + \omega_n^2 \cdot q_n = \Gamma_n \cdot p(t)$$
Para a excitação harmónica acima definida, a resposta modal em regime estacionário é dada por:
$$q_n(t) = \frac{\Gamma_n}{K_n} \cdot \frac{(1 - \beta_n^2) \cdot \sin (\omega \cdot t) - 2 \cdot \zeta_n \cdot \beta_n \cdot \cos (\omega \cdot t)}{(1- \beta_n^2)^2 + (2 \cdot \zeta_n \cdot \beta_n)^2}$$
em que $\beta_n = \frac{\omega}{\omega_n}$ é a relação entre a frequência de excitação e a frequência de vibração modal. Naquela expressão o termo $\frac{\Gamma_n}{K_n} = q_{0n}$ corresponde à amplitude em deslocamento da resposta modal, dado que o outro fator é adimensional. De facto:
$$\lim_{\beta \to 0} q_n(t) = \frac{\Gamma_n}{K_n} \cdot \sin (\omega \cdot t) = q_{0n} \cdot \sin (\omega \cdot t)$$
A resposta em deslocamento do oscilador nas coordenadas reais é dada por:
$${u(t)} = [\Phi] \times {q(t)} = \sum_{n=1}^{N} {\phi_n} \cdot q_n(t)$$
Solução alternativa
A excitação aplicada produz deslocamentos estáticos dados por:
$${u_{st}} = [K^{-1}] \times {s}$$
Realizando a expansão modal deste vetor de deslocamentos, obtêm-se as coordenadas modais estáticas correspondentes:
$$q_{n,st} = \frac{{\phi_n}^T \times [M] \times {u_{st}}}{M_n}$$
As respostas modais são então dadas em função destas coordenadas:
$$q_n(t) = q_{n,st} \cdot \frac{(1 - \beta_n^2) \cdot \sin (\omega \cdot t) - 2 \cdot \zeta_n \cdot \beta_n \cdot \cos (\omega \cdot t)}{(1- \beta_n^2)^2 + (2 \cdot \zeta_n \cdot \beta_n)^2}$$
e a resposta em deslocamento do oscilador nas coordenadas reais, à semelhança do caso anterior, é dada por:
$${u} = [\Phi] \times {q(t)} = \sum_{n=1}^{N} {\phi_n} \cdot q_n(t)$$
Comparação de soluções
Conforme se observa, as duas soluções diferem apenas no termo multiplicativo da resposta modal pelo que serão equivalentes se e só se $q_{0n} = q_{n,st}$, isto é:
$$\frac{{\phi_n}^T \times {s}}{M_n \cdot K_n} = \frac{{\phi_n}^T \times [M] \times [K^{-1}] \times {s}}{M_n}$$
Eliminando os fatores comuns, tem-se:
$$\frac{1}{K_n} = [M] \times [K^{-1}]$$
$$\frac{1}{K_n} \cdot [K] = [M]$$
Pré-multiplicando por ${\phi_n}^T$ e pós-multiplicando por ${\phi_n}$ ambos os lados da igualdade:
$$\frac{1}{K_n} \cdot {\phi_n}^T \times [K] \times {\phi_n} = {\phi_n}^T \times [M] \times {\phi_n}$$
$$\frac{K_n}{K_n} = M_n$$
$$1 = M_n$$
Conclui-se assim que as duas soluções serão equivalentes se os modos de vibração forem normalizados relativamente à massa. Este é um aspeto que habitualmente não é tratado nos livros de dinâmica de estruturas.
Exemplo de aplicação
Laboratório computacional
End of explanation
MM = np.matrix(np.diag([1.,2.,3.]))
KK = np.matrix([[2,-1,0],[-1,2,-1],[0,-1,1]])*1000.
print(MM)
print(KK)
Explanation: Dados
Oscilador com três graus de liberdade:
End of explanation
W2, F1 = np.linalg.eig(KK.I@MM)
print(W2)
print(F1)
Explanation: Análise modal:
End of explanation
ix = np.argsort(W2)[::-1]
W2 = 1./W2[ix]
F1 = F1[:,ix]
Wn = np.sqrt(W2)
print(W2)
print(Wn)
print(F1)
Explanation: Ordenação dos modos:
End of explanation
Fn = F1/np.sqrt(np.diag(F1.T@MM@F1))
print(Fn)
Explanation: Normalização dos modos de vibração relativamente à matriz de massa:
End of explanation
Mn = np.diag(Fn.T@MM@Fn)
Kn = np.diag(Fn.T@KK@Fn)
print(Mn)
print(Kn)
Explanation: Massa e rigidez modais:
End of explanation
sp = np.matrix([[0.], [1.], [0.]])
wp = 2.*np.pi*4.
print(sp)
print(wp)
tt = np.arange(1000)*0.005
pt = np.sin(wp*tt)
plt.figure()
plt.plot(tt, pt)
plt.xlabel('Tempo (s)')
plt.ylabel('Força (kN)')
plt.show()
Explanation: Excitação dinâmica:
End of explanation
def qn(amplitude, wp, beta, zn, tt):
    Calcula a resposta modal do modo n.
    qn_t = amplitude * ((1.-beta**2)*np.sin(wp*tt)-2.*zn*beta*np.cos(wp*tt))/((1.-beta**2)**2+(2.*zn*beta)**2)
    return qn_t
Explanation: Função auxiliar para o cálculo da resposta modal:
End of explanation
Gn = np.diag(1./Mn)@Fn.T@sp
print(Gn)
qn_t = []
plt.figure()
for n in range(3):
    an = Gn[n]/Kn[n]
    bn = wp/Wn[n]
    q = qn(an[0,0], wp, bn, 0.05, tt)
    qn_t.append(q)
    plt.plot(tt, q, label='an={:.2e},bn={:.2f}'.format(an[0,0], bn))
plt.legend()
plt.xlabel('Tempo (s)')
plt.ylabel('Deslocamento modal (m)')
plt.show()
plt.figure()
u_t = Fn@qn_t
for n in range(3):
    plt.plot(tt, u_t[n].T, label='{:.2e}'.format(np.max(u_t[n])))
plt.legend()
plt.xlabel('Tempo (s)')
plt.ylabel('Deslocamento (m)')
plt.show()
Explanation: Solução clássica
End of explanation
ust = KK.I@sp
print(ust)
qn_st_t = []
plt.figure()
for n in range(3):
    an = Fn.T[n]@MM@ust/Mn[n]
    bn = wp/Wn[n]
    qst = qn(an[0,0], wp, bn, 0.05, tt)
    qn_st_t.append(qst)
    plt.plot(tt, qst, label='an={:.2e},bn={:.2f}'.format(an[0,0], bn))
plt.legend()
plt.xlabel('Tempo (s)')
plt.ylabel('Deslocamento modal (m)')
plt.show()
plt.figure()
u_t = Fn@qn_st_t
for n in range(3):
    plt.plot(tt, u_t[n].T, label='{:.2e}'.format(np.max(u_t[n])))
plt.legend()
plt.xlabel('Tempo (s)')
plt.ylabel('Deslocamento (m)')
plt.show()
Explanation: Solução alternativa
End of explanation
<END_TASK>