In [1]:
%matplotlib inline
import sys
package_directory = r'C:\Users\ntotaro\Documents\TOOLBOX\PYTHON' #remplacer ici par l'adresse du dossier SiRoCo_Project dans votre arbrescnce
sys.path.append(package_directory)
import numpy as np #import des fonctions mathématiques du package numpy
import matplotlib.pyplot as plt
import os
os.chdir(package_directory + r'\SiRoCo_Project\CHECK\Vol1')
Soit une cavité d'air $\Omega$, à parois rigides, dont les caractéristiques sont :
La réponse de la cavité peut se décomposer sur la base de ses modes propres (parois rigides) telle que : \begin{equation} p(M) = \sum_{n} a_n(\omega)\phi_n(M) \label{eq3} \end{equation} où $\phi_n(M)$ est le mode $n$ au point $M$ et $a_n(\omega)$ est l'amplitude modale du mode $n$ à la pulsation $\omega$. Les modes respectent l'équation de Helmholtz sans second membre : \begin{equation} \Delta \phi_n(M) + k_n^2 \phi_n(M) = 0 \end{equation} Les modes de cavité respectent aussi la condition de Neumann homogène sur les conditions aux limites : \begin{equation} \frac{\partial \phi_n}{\partial n} (M) = 0 \end{equation}
Soit une source ponctuelle $S$, placée au point $M_0$, de débit volumique $q$ ($kg/s$). L'équation de Helmholtz s'écrit alors : \begin{equation} \Delta p(M) + k^2 p(M) = -j\omega q \delta(M-M_0) \label{eq1} \end{equation} où l'opérateur Laplacien $\Delta(\bullet) = \nabla^2(\bullet) = \frac{\partial^2}{\partial x^2}(\bullet) + \frac{\partial^2}{\partial y^2}(\bullet) + \frac{\partial^2}{\partial z^2}(\bullet)$. $k=\frac{\omega}{c}$ est le nombre d'onde à la pulsation $\omega$. La pression $p(M)$ respecte les conditions aux limites de la cavité : \begin{equation} \frac{\partial p}{\partial n} (M) = 0 \end{equation} ce qui représente une condition de Neumann homogène (vitesse nulle) sur toute la surface $\Sigma$ de la cavité. L'identité de Green, appliquée à la pression $p(M)$ et aux modes $\phi_n(M)$, peut s'écrire : \begin{equation} \int_\Omega (\Delta p(M) + k^2 p(M))\phi_n(M) dM = \int_\Omega (\Delta \phi_n(M) + k^2 \phi_n(M))p(M) dM + \int_\Sigma p(Q)\frac{\partial \phi_n}{\partial n}(Q)-\phi_n(Q)\frac{\partial p}{\partial n}(Q) dQ \label{eqgreen} \end{equation} Compte tenu des équations de Helmholtz (pour $p(M)$ et $\phi_n(M)$) et des conditions aux limites associées, il vient : \begin{equation} -j\omega q\phi_n(M_0) = \int_V (k^2-k_n^2)\phi_n(M)p(M) dM \end{equation} A l'aide de la décomposition modale de la pression $p(M)$ (Eq. \ref{eq3}) et de la propriété d'orthogonalité des modes propres, on peut exprimer l'amplitude modale comme : \begin{equation} a_n(\omega) = \frac{j\omega q\phi_n(M_0)}{\Lambda_n(k_n^2-k^2)} \end{equation} où $\Lambda_n$ est la norme du mode $n$. Finalement, la pression au point $M$ lorsqu'un monopole est placé en $M_0$, s'écrit : \begin{equation} p(M) = \sum_n \frac{j\omega q\phi_n(M_0) \phi_n(M)}{\Lambda_n(k_n^2-k^2)} \end{equation}
Soit un élément de surface $S_p$ ($S_p \in \Sigma$, surface enveloppe de la cavité) vibrant à la vitesse V. Considérons qu'aucune source n'est présente dans le volume $\Omega$ de la cavité. L'équation de Helmholtz s'écrit alors : \begin{equation} \Delta p(M) + k^2 p(M) = 0 \label{eq_helm_0} \end{equation} Sur la surface de la cavité ($\Sigma-S_p$), la condition de Neumann homogène s'applique : \begin{equation} \frac{\partial p}{\partial n} (M) = 0 \quad \forall M \in (\Sigma-S_p) \end{equation} Sur la surface $S_p$, la condition de vitesse est imposée par les vibrations du piston : \begin{equation} \frac{\partial p}{\partial n} (M) = -j\rho\omega V \quad \forall M \in S_p \end{equation} En utilisant les modes de la cavité à parois rigides dans ce cas sur l'Eq. (\ref{eqgreen}), il vient : \begin{equation} a_n(\omega) = \frac{j\omega \rho V S_p \left\langle\phi_n\right\rangle_{S_p}}{\Lambda_n(k_n^2-k^2)} \end{equation} où \begin{equation} \left\langle\phi_n\right\rangle_{S_p} = \frac{1}{S_p} \int_{S_p} \phi_n(M) dM \end{equation} Finalement, la pression au point $M$ lorsqu'un piston de surface $S_p$ vibre, s'écrit : \begin{equation} p(M) = \sum_n \frac{j\omega \rho V S_p\left\langle\phi_n\right\rangle_{S_p} \phi_n(M)}{\Lambda_n(k_n^2-k^2)} \end{equation}
L'amortissement est généralement introduit par l'intermédiaire d'une célérité complexe de la forme $c^*=c(1+j\eta)$. Cette célérité complexe s'applique uniquement à $k$ et non à $k_n$. Ainsi, on obtient : \begin{equation} p(M) = \sum_n \frac{j\omega q\phi_n (M_0) \phi_n (M)}{\Lambda_n (k_n^2-{k^*}^2)} = \sum_n \frac{j\omega q.c^2.\phi_n(M_0) \phi_n(M)}{\Lambda_n( \omega_n^2-\omega^2/(1+j\eta)^2)} = \sum_n \frac{j\omega q.c^2.\phi_n(M_0) \phi_n(M)}{\Lambda_n( \omega_n^2-\omega^2(1-2j\eta))} \label{eqp} \end{equation} pour le monopole et \begin{equation} p(M) = \sum_n \frac{j\omega \rho c^2 V S_p \left\langle\phi_n\right\rangle_{S_p} \phi_n(M)}{\Lambda_n( \omega_n^2-\omega^2(1-2j\eta))} \label{eqp2} \end{equation} pour le piston vibrant
L'intensité acoustique se déduit de la pression et la vitesse des particules d'air par la relation suivante : \begin{equation} \vec{I}(M) = \frac{1}{2}p(M)\vec{v}^*(M) \end{equation} où $\bullet^*$ représente le complexe conjugué. La partie imaginaire de l'intensité est l'intensité réactive et la partie réelle de l'intensité est l'intensité active.
La densité d'énergie cinétique est liée à la norme du vecteur vitesse par l'équation : \begin{equation} e_c(M) = \frac{1}{4}\rho |\vec{v}(M)|^2 \end{equation} La densité d'énergie potentielle est liée à la pression quadratique locale : \begin{equation} e_p(M) = \frac{1}{4} \frac{p(M)p^*(M)}{\rho c^2} \end{equation}
Pour déclarer une telle cavité, on utilise la classe classAnalyticCavity :
In [2]:
from SiRoCo_Project.AnalyticalStructuralAcoustics.Cavity.AnalyticCavity import *
monVolume = AnalyticCavity(Lx=0.4, Ly=0.6, Lz=0.25, density=1.29, celerity=340, damping=0.01)
Les paramètres density, celerity et damping sont optionnels et ont comme valeurs par défaut 1.29, 340 et 0.01. Les labels des paramètres Lx, Ly et Lz peuvent être ommis si l'ordre est respecté. Ainsi, la déclaration précédente est en tout point identique à :
In [3]:
monVolume = AnalyticCavity(0.4,0.6,0.25)
Dans le cas d'une cavité parallélépipédique, les pulsations propres sont données par : \begin{equation} \omega_{mnp} = c \sqrt{ \left(\frac{m\pi}{L_x}\right)^2 + \left(\frac{n\pi}{L_y}\right)^2 + \left(\frac{p\pi}{L_z}\right)^2 } \end{equation} les déformées propres par : \begin{equation} \phi_{mnp}(x,y,z) = cos\left(\frac{m\pi}{L_x}x\right)cos\left(\frac{n\pi}{L_y}y\right)cos\left(\frac{p\pi}{L_z}z\right) \label{eqphi} \end{equation} et les normes des modes par : \begin{equation} \Lambda_{mnp} = \frac{L_xL_yL_z}{8} 2^d \end{equation} où (m,n,p) est le triplet d'indices du mode et
Pour calculer les modes de \textit{monVolume}, il faut utiliser la méthode \textit{setFmax} qui permet de définir la fréquence maximum de calcul des modes :
In [4]:
monVolume.setFmax(2000)
La méthode calcule les fréquences propres et les triplets d'indices de modes correspondants. On peut ensuite accéder à ces fréquences propres et ces indices par les méthodes \textit{getEigenFreq} et \textit{getEigenIndices}:
In [5]:
fmnp = monVolume.getEigenFreq()
mnp = monVolume.getEigenIndices()
# Affichage des données (par exemple)
for f, n in zip(fmnp, mnp):
print('Mode ' + str(n) + ', Fréquence = ' + str(f) + 'Hz.')
In [6]:
fid = open('freqp_vol1.txt')
data = fid.read()
f_nastran = eval('[' + data + ']')
Comparaison
In [7]:
plt.plot(fmnp,'s')
plt.hold(True)
plt.plot(f_nastran,'o')
plt.xlabel('Mode index')
plt.ylabel('Frequency (Hz)')
plt.legend(('Calcul analyticCavity','Calcul Nastran'),loc=4)
Out[7]:
Compte tenu de la forme des modes (Eq. \ref{eqphi}) et des Eqs. (\ref{eqp}) et (\ref{eqv}), on peut calculer la réponse en fréquence de la cavité en un point M défini par les coordonnées $(x_M, y_M, z_M)$ (voir Fig. 1) lorsqu'un monopole de débit volumique q est placé en un point S de coordonnées $(x_S, y_S, z_S)$.
In [8]:
s1 = Source(monopoleQ = ([0.091,0.18,0.1471],1))
m1 = Sensor(point = ([0.1062,0.3,0.235],))
monVolume.addSource(s1)
monVolume.addSensor(m1)
freq = np.linspace(10.,2000,1991)
res = monVolume.computeAtFreq(freq)
La variable res est un dictionnaire contenant les résultats de calculs pour la pression (scalaire), la vitesse (vecteur) et l'intensité (vecteur) et les densité d'énergie potientielle et cinétique (scalaires) au point de mesure m1.
In [9]:
for type_data in res.keys():
print(type_data)
Par exemple, pour accéder aux données de pression il faut utiliser la syntaxe suivante :
In [10]:
p = res['pressure']
print(len(p))
print(len(p[0]))
#equivalent à
print(len(res['pressure']))
print(len(res['pressure'][0]))
Ce sont des listes de résultats. Dans le cas où N "sensors" sont définis, les listes sont de taille N (ici N=1 donc len(p)=1). Ainsi, dans le cas d'un seul sensor, on accède à la réponse fréquentielle en pression par res['pressure'][0]. C'est une donnée scalaire pour chaque fréquence. res['pressure'][0] est donc de type np.ndarray et de dimensions 1xnb_freq (ici 1991). On peut ainsi par exemple tracer l'amplitude de la pression en fonction de la fréquence par :
In [11]:
plt.plot(freq, np.abs(res['pressure'][0]))
Out[11]:
Concernant les données de vitesse et d'intensité, res['velocity'][0] et res['intensity'][0] sont de listes de dimensions 1x3, correspondant aux données sur les axes x, y, et z. Ainsi on peut accéder à la réponse fréquentielle de l'intensité sur l'axe y par res['intensity'][0][1]. La partie réelle de l'intensité est l'intensité active et la partie imaginaire est l'intensité réactive.
Le calcul Actran porte sur la réponse d'une cavité de dimensions 0.4 x 0.6 x 0.25m remplie d'air (masse volumique $\rho_0 = 1.29$kg/m^3, célérité $c=340$m/s, amortissement $\eta=0.01$). L'amortissement est introduit par l'intermédiaire d'une célérité complexe $c^*=c(1+j\eta)$. Un monopole de débit volumique $q=1$ (voir Eq. \ref{eq1}) est placé à la position (0.091;0.18;0.1471)m. Sous Actran, cette source correspond à un monopole de type Q. Le point d'écoute est en (0.1062;0.3;0.235)m. La bande de fréquence de calcul va de 10Hz à 2000Hz par pas de 1Hz.
Chargement des fichiers de résultats Actran :
In [12]:
from SiRoCo_Project.tools.readingFilesTools import readASCII
#vecteur fréquence fréquences de calcul
Preal_actran = readASCII('Preal.txt', ' \t ')
freq_actran = Preal_actran[:,0]
#chargement de la pression (partie réelle, partie imaginaire)
Pimag_actran = readASCII('Pimag.txt', ' \t ')
P_actran = Preal_actran[:,1] + 1j*Pimag_actran[:,1]
#Chargement de la vitesse suivant x, y et z (partie réelle, partie imaginaire)
Vxreal_actran = readASCII('Vxreal.txt', ' \t ')
Vximag_actran = readASCII('Vximag.txt', ' \t ')
Vx_actran = Vxreal_actran[:,1] + 1j*Vximag_actran[:,1]
Vyreal_actran = readASCII('Vyreal.txt', ' \t ')
Vyimag_actran = readASCII('Vyimag.txt', ' \t ')
Vy_actran = Vyreal_actran[:,1] + 1j*Vyimag_actran[:,1]
Vzreal_actran = readASCII('Vzreal.txt', ' \t ')
Vzimag_actran = readASCII('Vzimag.txt', ' \t ')
Vz_actran = Vzreal_actran[:,1] + 1j*Vzimag_actran[:,1]
#Calcul de l'intensité suivant x, y et z
Ix_actran = 0.5*P_actran*np.conjugate(Vx_actran)
Iy_actran = 0.5*P_actran*np.conjugate(Vy_actran)
Iz_actran = 0.5*P_actran*np.conjugate(Vz_actran)
Ces résultats du calcul Actran peuvent ainsi être comparés aux calculs analyticCavity :
In [13]:
plt.figure()
f, (ax1, ax2) = plt.subplots(2, sharex=True)
ax1.plot(freq_actran, np.real(P_actran),label="Actran")
ax1.plot(freq, np.real(res['pressure'][0]),label="analyticCavity")
ax2.plot(freq_actran, np.imag(P_actran))
ax2.plot(freq, np.imag(res['pressure'][0]))
ax2.set_xlabel('Frequency (Hz)')
ax1.set_ylabel('real part (Pa)')
ax2.set_ylabel('imag part (Pa)')
ax1.set_title('Pressure')
ax1.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.figure()
f, (ax1, ax2) = plt.subplots(2, sharex=True)
ax1.plot(freq_actran, np.real(Vx_actran), label="Actran")
ax1.plot(freq, np.real(res['velocity'][0][0]), label="analyticCavity")
ax2.plot(freq_actran, np.imag(Vx_actran))
ax2.plot(freq, np.imag(res['velocity'][0][0]))
ax2.set_xlabel('Frequency (Hz)')
ax1.set_ylabel('real part (Pa)')
ax2.set_ylabel('imag part (Pa)')
ax1.set_title('Velocity /x')
ax1.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.figure()
f, (ax1, ax2) = plt.subplots(2, sharex=True, sharey=True)
ax1.plot(freq_actran, np.real(Vy_actran), label="Actran")
ax1.plot(freq, np.real(res['velocity'][0][1]), label="analyticCavity")
ax2.plot(freq_actran, np.imag(Vy_actran))
ax2.plot(freq, np.imag(res['velocity'][0][1]))
ax2.set_xlabel('Frequency (Hz)')
ax1.set_ylabel('real part (Pa)')
ax2.set_ylabel('imag part (Pa)')
ax1.set_title('Velocity /y')
ax1.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.figure()
f, (ax1, ax2) = plt.subplots(2, sharex=True)
ax1.plot(freq_actran, np.real(Vz_actran), label="Actran")
ax1.plot(freq, np.real(res['velocity'][0][2]), label="analyticCavity")
ax2.plot(freq_actran, np.imag(Vz_actran))
ax2.plot(freq, np.imag(res['velocity'][0][2]))
ax2.set_xlabel('Frequency (Hz)')
ax1.set_ylabel('real part (Pa)')
ax2.set_ylabel('imag part (Pa)')
ax1.set_title('Velocity /z')
ax1.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.figure()
f, (ax1, ax2) = plt.subplots(2, sharex=True)
ax1.plot(freq_actran, np.real(Ix_actran), label="Actran")
ax1.plot(freq, np.real(res['intensity'][0][0]), label="analyticCavity")
ax2.plot(freq_actran, np.imag(Ix_actran))
ax2.plot(freq, np.imag(res['intensity'][0][0]))
ax2.set_xlabel('Frequency (Hz)')
ax1.set_ylabel('real part (Pa)')
ax2.set_ylabel('imag part (Pa)')
ax1.set_title('Intensity /x')
ax1.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.figure()
f, (ax1, ax2) = plt.subplots(2, sharex=True)
ax1.plot(freq_actran, np.real(Iy_actran), label="Actran")
ax1.plot(freq, np.real(res['intensity'][0][1]), label="analyticCavity")
ax2.plot(freq_actran, np.imag(Iy_actran))
ax2.plot(freq, np.imag(res['intensity'][0][1]))
ax2.set_xlabel('Frequency (Hz)')
ax1.set_ylabel('real part (Pa)')
ax2.set_ylabel('imag part (Pa)')
ax1.set_title('Intensity /y')
ax1.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
plt.figure()
f, (ax1, ax2) = plt.subplots(2, sharex=True)
ax1.plot(freq_actran, np.real(Iz_actran), label="Actran")
ax1.plot(freq, np.real(res['intensity'][0][2]), label="analyticCavity")
ax2.plot(freq_actran, np.imag(Iz_actran))
ax2.plot(freq, np.imag(res['intensity'][0][2]))
ax2.set_xlabel('Frequency (Hz)')
ax1.set_ylabel('real part (Pa)')
ax2.set_ylabel('imag part (Pa)')
ax1.set_title('Intensity /z')
ax1.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
Out[13]:
Plusieurs types de sources et de "capteurs" sont disponibles. Ils sont détaillés ci-dessous
Les Sensors disponibles sont :
In [14]:
from SiRoCo_Project.tools.classCartesianLine import * #import du package cartesianLine (contient les objets gridLine)
from SiRoCo_Project.tools.classCartesianMap import * #import du package cartesianMap (contient les objets gridMap)
m1 = Sensor(point = ([0.2,0.3,0.4],)) #déclaration d'un point de mesure de coordonnées (0.2,0.3,0.4). Attention, point prend comme paramètre un tuple de dimension 1x2 (le deuxième est vide)
m2 = Sensor(elSurf = ([0.2,0.3,0.4],[0.1,0.2,0.])) # déclaration d'une elSurf de point de coordonnées (0.2,0.3,0.4) et de longueur (0.1,0.2). La longueur suivant z étant nulle, la surface est perpendiculaire à z
line1 = gridLine3D(x0=0., y0=0.1, z0=0.0, L=0.5, d=0.01, theta=0., phi=0.) # déclaration d'une gridLine3D (voir package cartesianLine)
m3 = Sensor(gridLine = line1) # déclaration d'une gridLine3D comme Sensor
map1 = gridMap3D(x0=0., y0=0., z0=0.1, Lx=0.5, Ly=0.6, dx=0.02, dy=0.02, theta=0., phi=0.) # déclaration d'une gridMap3D (voir package cartesianMap)
m4 = Sensor(gridMap = map1) # déclaration d'une gridMap3D comme Sensor</div><i class="fa fa-lightbulb-o "></i>
Il suffit ensuite d'ajouter successivement tous les "Sensor" déclarés (m1...m4) pour qu'ils soient pris en compte dans le calcul. Les "Source" et "Sensor" ajoutés précédemment peuvent être supprimés en utilisant respectivement les méthodes "deleteSources" et "deleteSensors" de l'objet analyticCavity. Nous supprimerons uniquement les "Sensor" dans l'exemple suivant :
In [15]:
monVolume.deleteSensors()
monVolume.addSensor(m1)
monVolume.addSensor(m2)
monVolume.addSensor(m3)
monVolume.addSensor(m4)
print(monVolume)
On peut ensuite lancer le calcul de réponse en fréquence sur ces données, comme précédemment :
In [16]:
res = monVolume.computeAtFreq(freq)
Les résultats resP, resV et resI sont des listes (ici de dimension 4 car nous avons déclaré 4 "Sensor"). Les données de sortie sont de même type que les "Sensor" déclarés. Ainsi, pour la pression, resP[0] est une liste de valeurs scalaires correspondant au calcul au "point" pour chaque fréquence, resP[1] est une liste de valeurs scalaires correspondant au calcul sur la "elSurf" pour chaque fréquence, resP[2] est une liste de gridLine3D calculées pour chaque fréquence, resP[3] est une liste de gridMap3D calculées pour chaque fréquence. On peut ainsi tracer la réponse en fréquence au "point" ou la cartographie sur la gridMap. Par exemple :
In [17]:
#Tracé de la réponse en fréquence sur le "point"
plt.figure()
plt.semilogy(freq, np.abs(res['pressure'][0]))# trace la réponse en fréquence au "point"
plt.xlabel('Frequency (Hz)')
plt.ylabel('FRF')
#Tracé de la cartographie sur la "gridMap" à la fréquence d'indice 526
res['pressure'][3][526].plot()
Pour les données vectorielles (vitesse, intensité) et dans le cas de "Sensor" de type "point", resV[i] est de taille 3 correspondant aux 3 directions. Ainsi, resV[0][1] est la liste des vitesses calculées suivant l'axe y pour toutes les fréquences et resI[0][0] est la liste des intensités calculées suivant l'axe x pour toutes les fréquences. Pour les "Sensor" de type "elSurf", c'est le même principe mais seule la composante normale à la surface est non-nulle. Dans le cas des "Sensor" de type "gridLine" ou "gridMap", resV[i] est la liste des "gridLine" (ou "gridMap") calculées pour chaque fréquence. Les objets "gridLine" et "gridMap" permettent de gérer les données vectorielles. Par exemple, resV[3][526] est la "gridMap" calculée pour la fréquence d'indice 526. Comme c'est une données vectorielle, on a deux méthodes d'affichage :
In [18]:
res['velocity'][526].plot(dataDirection = 'x', displayMode = 'ReIm', contour = 'on')
ou :
In [ ]:
res['velocity'][3][526].plotArrows(displayMode = 'Re')
Les sources disponibles sont :
In [ ]:
s1 = Source(monopoleQ = ([0.2,0.3,0.15],1)) #monopole de type Q, placé en (0.2,0.3,0.15) et dont q=1
s2 = Source(monopoleP = ([0.1,0.2,0.1],-1)) #monopole de type P, placé en (0.1,0.2,0.1) et dont q'=-1
s3 = Source(elSurf = ([0.,0.,0.1],[0.1,0.1,0.],0.5)) #piston dont le coin inférieur gauche est en (0.,0.,0.1), dont les dimensions sur x et y sont 0.1 et 0.1 (il est donc perpendiculaire à z) et dont V=0.5
Il suffit ensuite d'ajouter successivement tous les \textit{Source} déclarées (s1...s4) pour qu'elles soient pris en compte dans le calcul. Attention la réponse calculée est la réponse à l'ensemble des sources ajoutées dans \textit{monVolume}. Les \textit{Source} et ajoutées précédemment peuvent être supprimés en utilisant la méthode \textit{deleteSources} de l'objet \textit{analyticCavity}. Par exemple :
In [ ]:
monVolume.deleteSources()
monVolume.addSource(s1)
monVolume.addSource(s2)
monVolume.addSource(s3)
print(monVolume)
Enfin, les réponses sont calculées, toujours de la même manière par :
In [ ]:
res = monVolume.computeAtFreq(freq)
La densité modale (en nombre de modes par Hertz) d'une cavité parallélépipédique est donnée par la relation suivante (LESUEUR, p. 81) : \begin{equation} n(f) = \frac{4\pi f^2 V}{c^3} + \frac{\pi f S}{2c^2} + \frac{L}{8c} \end{equation} avec S la surface enveloppe de la cavité et L la somme des arêtes de la cavité. La méthode \textit{getModalDensity} permet de calculer la densité modale de la cavité en fonction de la fréquence :
In [ ]:
n = monVolume.getModalDensity(freq)
plt.figure()
plt.plot(freq,n)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Modal Density (modes/Hz)')
In [ ]:
nbf = monVolume.nbEigenFreq() # renvoie le nombre de modes calculés
print("{0} modes calculés".format(nbf))
nbm = monVolume.nbEigenFreq(fMin = 500, fMax = 1000) # renvoie le nombre de modes calculés dont les fréquences propres sont supérieures à 500Hz et inférieures ou égales à 1000Hz
print("{0} modes calculés ont une fréquence supérieure à 500Hz et inférieure ou égale à 1000Hz".format(nbm))
nbs = monVolume.nbSources() # renvoie le nombre de sources associées à monVolume
print("{0} sources ont été affectées à monVolume".format(nbs))
nbp = monVolume.nbSensors() # renvoie le nombre de sensors associés à monVolume
print("{0} sensors ont été affectés à monVolume".format(nbp))
In [ ]:
listEigenFreq = monVolume.getEigenFreq() #liste totale des fréquences propres
listEigenFreq_partial = monVolume.getEigenFreq(fMin = 500, fMax = 1000) #liste des fréquences propres supérieures à 500Hz (exclu) et inférieures à 1000Hz (inclus)
listEigenIndices = monVolume.getEigenIndices(fMax = 800)# liste des triplets d'indices dont les modes ont une fréquence inférieure à 800Hz (inclus)
In [ ]:
V = monVolume.getVolume() # renvoie le volume de la cavité
print(V)