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 : \begin{itemize} \item une célérité $c$ ($m/s$) des ondes \item une densité $\rho$ ($kg/m^3$) du fluide \item un volume $V$ \end{itemize} La pression dans le volume est notée $p(M)$, elle respecte l'équation de Helmholtz.
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} La partie imaginaire de l'intensité est l'intensité active et la partie réelle de l'intensité est l'intensité réactive.
Pour déclarer une telle cavité, on utilise la classe classAnalyticCavity :
In [2]:
from SiRoCo_Project.AnalyticalStructuralAcoustics.Cavity.classAnalyticCavity 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 \begin{itemize} \item d=3 si les 3 indices sont nuls \item d=2 si 2 indices sont nuls \item d=1 si 1 indice est nuls \item d=0 si aucun indice n'est nul \end{itemize}
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)
Les variables resP, resV et resI sont les résultats de calculs pour respectivement la pression (scalaire), la vitesse (vecteur) et l'intensité (vecteur) au point de mesure m1. Ce sont des listes de résultats. Dans le cas où N "sensors" sont définis, les listes sont de taille N. Ainsi, dans le cas d'un seul sensor, on accède à la réponse fréquentielle en pression par resP[0]. C'est une donnée scalaire pour chaque fréquence. resP[0] est donc de type np.ndarray et de dimensions 1xnb_freq. On peut ainsi par exemple tracer l'amplitude de la pression en fonction de la fréquence par :
In [9]:
plt.semilogy(freq, np.abs(resP[0]))
Concernant les données de vitesse et d'intensité, resV[0] et resI[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 resI[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 [ ]:
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 [ ]:
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.)
Plusieurs types de sources et de "capteurs" sont disponibles. Ils sont détaillés ci-dessous
Les Sensors disponibles sont : \begin{itemize} \item "point" : un point particulier défini par ses coordonnées \item "elSurf" : une surface élémentaire sur laquelle est intégrée la solution (intégration exacte). Pour la vitesse et l'intensité, seules les quantités normales à la surface sont calculées. Une "elSurf" est définie par un point dans l'espace (coin inférieur gauche du quadrangle) et par 3 longueurs dont une doit être nulle (surface élémentaire perpendiculaire à la direction de la longueur nulle). \item "gridline" : définit une ligne de points régulièrement espacés. Cette ligne doit être de type gridLine3D (package cartesianLine) \item "gridmap" : définit une grille de points régulièrement espacés. Cette grille doit être de type gridMap3D (package cartesianMap) \end{itemize} Ci-dessous, un exemple de déclaration de chacun des types de Sensors :
In [10]:
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 [11]:
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 [12]:
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 [13]:
#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 [ ]:
res['velocity'][526].plot(dataDirection = 'x', displayMode = 'ReIm', contour = 'on')
ou :
In [ ]:
res['velocity'][3][526].plotArrows(displayMode = 'Re')
Les sources disponibles sont : \begin{itemize} \item "monopoleQ" : monopole de type débit volumique (modélisation monopole type Q sous Actran). Correspond à la modélisation de l'Eq. (\ref{eq16}) (on impose la valeur de $q$ dans l'équation). Un "monopoleQ" est défini par ses coordonnées et sa valeur. \begin{equation} \Delta p(M) + k^2 p(M) = -j\omega q \delta(M-M_0) \label{eq16} \end{equation} \item "monopoleP" : monopole de type variation de débit volumique (modélisation monopole de type P sous Actran. Correspond à la modélisation de l'Eq (\ref{eq17})(on impose la valeur de $\dot{q}$ dans l'équation). Un "monopoleP" est défini par ses coordonnées et sa valeur. \begin{equation} \Delta p(M) + k^2 p(M) = -\dot{q} \delta(M-M_0) \label{eq17} \end{equation} \item "elSurf" : vitesse imposée sur une surface élémentaire. Une "elSurf" est définie par un point dans l'espace (coin inférieur gauche du quadrangle) et par 3 longueurs dont une doit être nulle (surface élémentaire perpendiculaire à la direction de la longueur nulle) \end{itemize} Ci-dessous, un exemple de déclaration de chacun des types de Sensors :
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)