# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Acoplamento mecânico - acústico

# +
import numpy as np
from scipy import special
from matplotlib import pyplot as plt
plt.rcParams.update({'font.size': 14})

def zrad_piston_infbaffle(a, freq, rho0 = 1.21, c0 = 343):
    # Equivalent circuit
    Ra1 = (0.1404 * rho0 * c0) / (a**2)
    Ra2 = (rho0 * c0)/(np.pi * a**2)
    Ca1 = 5.94 *a**3/(rho0 * c0**2)
    Ma1 = (0.27 * rho0) / a
    # Equivalent impedance
    omega = 2*np.pi*freq
    Zrad_eq=(-(omega**2)*Ma1*Ra1*Ra2*Ca1+1j*omega*Ma1*(Ra1+Ra2))/\
        (Ra1+Ra2-(omega**2)*Ma1*Ra1*Ca1+1j*omega*(Ra1*Ra2*Ca1+Ma1))

    # Radiation impedance 
    k0 = omega/c0; 
    #s = pi*a^2;
    Jka = special.jv(1,2*k0*a)
    Ska = special.struve(1,2*k0*a)

    R1 = 1 - 2*Jka/(2*k0*a);
    X1 = 2*Ska/(2*k0*a);

    Zrad = (rho0*c0/(np.pi*a**2))*(R1 + 1j*X1)
    return Ra1, Ra2, Ca1, Ma1, Zrad_eq, Zrad


# -

# # Exemplo 1: Pistão no baffle infinito
#
# Seja um pistão de raio $a$ e massa $M_m$, acoplado a um sistema de suspensão com compliância $C_m$ e amortecimento $R_m$. O pistão sofre ação de uma força $f(t)$ e está montado em um baffle infinito e, por isso sofre ação da impedância de radiação acústica.
#
# ![image.png](attachment:image.png)
#
# Os circuitos do sistema mecânico-acústico acoplado por um transformador e com o lado acústico referido ao lado mecânico são dados na figura.
#
# Se conseguirmos mensurar a velocidade de vibração do pistão, $\tilde{\dot{X}}$, podemos usar a Eq. (1) para mensurar a pressão sonora radiada na direção ***normal*** e a uma distância $r$ do pistão.
#
# \begin{equation}
# \tilde{P}(r,0^{\circ})=\text{j} 2 \rho_{0} c_{0} \tilde{\dot{X}} \ \sin\left[\frac{k_{0}}{2}\left(\sqrt{r^{2}+a^{2}}-r\right)\right] \ \text{e}^{-\text{j} (k_{0}/2) (\sqrt{r^{2}+a^{2}}+r)}
# \tag{1}
# \end{equation}

def pres_spkbaffle_odeg(freq, X, r, a, rho0 = 1.21, c0 = 343):
    w = 2 * np.pi * freq
    k0 = w/c0
    pres = 1j*2*rho0*c0 * X *\
        (np.sin((k0/2)*(np.sqrt(r**2+a**2)-r))) *\
        np.exp(-1j*(k0/2)*(np.sqrt(r**2+a**2)+r))
    return pres


# +
# dados acústicos
rho0 = 1.21
c0 = 343
freq = np.arange(1, 10000, 1)
w = 2 * np.pi * freq
k0 = w/c0

# Dados do pistão
a = (6/2) * 2.54e-2
Sp = np.pi*a**2
Mm = 0.05
Rm = 1.4
Cm = 0.001
Ra1, Ra2, Ca1, Ma1, Zrad_eq, Zrad = zrad_piston_infbaffle(a, freq, 
    rho0 = 1.21, c0 = 343)

# Frequência de ressonância mecânica
fs = (1/(2*np.pi)) * np.sqrt(1/(Cm*Mm))
print("A frequência de ressonância mecânica é {:.2f} [Hz]".format(fs))
# -

# ### Velocidade do pistão
#
# \begin{equation}
# \frac{F}{\tilde{\dot{X}}} = \text{j}\omega M_{m}+\frac{1}{\text{j}\omega C_{m}}+R_{m}+S_p^{2}Z_{\text{rad}}
# \tag{2}
# \end{equation}
#

# +
F = 1
X = (1j*w*Mm + 1/(1j*w*Cm) + Rm + (Sp**2)*Zrad)**(-1)
r = 1
P = pres_spkbaffle_odeg(freq, X, r, a, rho0 = rho0, c0 = c0)


# plot
fig, axs = plt.subplots(2,1, figsize = (7,5))
axs[0].semilogx(freq, 20*np.log10(np.abs(X)), linewidth = 2)
axs[0].grid(linestyle = '--', which='both')
#axs[0].set_ylim((-40, 5))
axs[0].set_xlim((1, 10000))
axs[0].set_ylabel(r'$|\tilde{\dot{X}}|$ [dB]')
axs[1].semilogx(freq, 20*np.log10(np.abs(P)), linewidth = 2)
axs[1].grid(linestyle = '--', which='both')
axs[1].set_ylabel(r'$|\tilde{P}|$ [dB]')
axs[1].set_xlabel(r'Frequência [Hz]')
axs[1].set_xlim((1, 10000))
fig.tight_layout()
# -

# # Exemplo 2: Pistão numa caixa acústica e baffle infinito
#
# Seja um pistão de raio $a$ e massa $M_m$, acoplado a um sistema de suspensão com compliância $C_m$ e amortecimento $R_m$. O pistão sofre ação de uma força $f(t)$ e está montado em um baffle infinito e, por isso sofre ação da impedância de radiação acústica. Além disso, o pistão a parte traseira do pistão é enclausurada por uma caixa selada com volume $V$ e compliância $C_a$.
#
# ![image.png](attachment:image.png)
#
# Os circuitos do sistema mecânico-acústico acoplado por um transformador e com o lado acústico referido ao lado mecânico são dados na figura.
#
# Da mesma forma, vamos mensurar a velocidade de vibração do pistão, $\tilde{\dot{X}}$ e usar a Eq. (1) para mensurar a pressão sonora radiada na direção ***normal*** e a uma distância $r$ do pistão.

# +
# Dados da caixa
V = 18/1000
Ca = V/(rho0*c0**2)
print("A compliância mecânica é {:.4f} [m/N]".format(Cm))
print("A compliância acústica referida ao lado mecânico é {:.4f} [m/N]".format(Ca/Sp**2))

Ceq = (1/Cm + 1/(Ca/Sp**2))**(-1)
print("A compliância mecânica equivalente é {:.5f} [m/N]".format(Ceq))

fs = (1/(2*np.pi)) * np.sqrt(1/(Ceq*Mm))
print("A frequência de ressonância aprox. é {:.2f} [Hz]".format(fs))

X_ca = (1j*w*Mm + 1/(1j*w*Cm) + Rm + (Sp**2)/(1j*w*Ca) + (Sp**2)*Zrad)**(-1)
P_ca = pres_spkbaffle_odeg(freq, X_ca, r, a, rho0 = rho0, c0 = c0)

# plot
fig, axs = plt.subplots(2,1, figsize = (7,5))
axs[0].semilogx(freq, 20*np.log10(np.abs(X)), 
                linewidth = 2, label = 'livre'),
axs[0].semilogx(freq, 20*np.log10(np.abs(X_ca)), 'r', 
                linewidth = 2, label = 'caixa')
axs[0].grid(linestyle = '--', which='both')
axs[0].legend()
axs[0].set_xlim((1, 10000))
axs[0].set_ylabel(r'$|\tilde{\dot{X}}|$ [dB]')
axs[1].semilogx(freq, 20*np.log10(np.abs(P)), 
                linewidth = 2, label = 'livre')
axs[1].semilogx(freq, 20*np.log10(np.abs(P_ca)), 'r',
                linewidth = 2, label = 'caixa')
axs[1].grid(linestyle = '--', which='both')
axs[1].set_ylabel(r'$|\tilde{P}|$ [dB]')
axs[1].set_xlabel(r'Frequência [Hz]')
axs[1].set_xlim((1, 10000))
fig.tight_layout()

StarcoderdataJupyterScriptsDedupFiltered