Ceci est une ancienne révision du document !
Thermoacoustique
équation | unité | coord. euleriennes (x,t) | définitions |
conservation de la masse | kg/m³ | (dr/dt)x +(d(ru)/dx)t =0 | u=(dx/dt)a =vitesse (faible: <1m/s) |
conservation de l'impulsion | N/m³ | d(ru)/dt+dP/dx=0 (oublie les frottements) | Ti =T+u²/(2*cp) |
travail reçu | W/m³ | dW/dx=-d(Pu)/dx | Pu est le work flow (en W/m²) |
transfert thermique | W/m³ | dQ/dx=8plDT/f² | DT=Tparoi-T (K) |
conservation de l'énergie | W/m³ | dH/dt=d(rcp*T)/dt=(dP/dt)*g/(g-1)=(dW/dx)+(dQ/dx) | |
équation d'état du gaz parfait | P=rrT | ||
gradient de paroi | K/m | d(Tp)/dx=dx T constant, d(Tp)/dt=0 | h=viscosité(Pa.s) |
frottement (parabolique) | (dP/dx)f =-32uh/f² | grave;s pratiques car les particules successives ne sont pas en début de cycle au même moment. Il vaut mieux remplacer la coordonnée lagrangienne a par une coordonnée de phase j qui varie de 0 à 1 quand la particule parcourt son cycle. Pour une particule donnée j est adimensionnelle et varie linéairement dans le temps: j=(t-t0)/T, où t est le temps, t0 le temps au début de cycle (point 0) pour la particule considérée et T la période (en secondes). Il faut faire attention à ce que le temps de période T varie éventuellement avec la particule. | |
Sur cette figure on a repéré en rouge les 4 points clés 0,1,2,3 du cycle dans la partie réchauffante de l'échangeur (dans la partie refroidissante le sens de propagation de l'onde sonore est inversé: en fonctionnement moteur, la partie haute pression du cycle doit toujours être aux hautes températures). Au point 0, la vitesse u est nulle, donc dP/dx>0 pour avoir une accélération. Ceci nous donne une relation entre la pression acoustique p=A*sin(w(t-x/a)) où a est la vitesse du son, et le gradient de pression temporel moyen dxP. -dp/dx=A*(w/a)*cos(w(t-x/a))>dxP au point 0, donc A*(w/a)>dxP. Comme w=2p*f, plus la fréquence est élevée et courte la longueur d'onde, plus on peut se contenter d'une faible amplitude. La vitesse u aussi est sinusoidale, à une constante près: u=B*sin(w(t-x/a)+D)+<u>, où D est le déphasage. On en tire du/dt=B*w*cos(w(t-x/a)+D) et du/dx=-B*(w/a)*cos(w(t-x/a)+D). En considérant que les efforts visqueux à la paroi compensent en moyenne le gradient de pression dxP, on a l'équation du/dt-u*(du/dx)=-dp/rdx (voir plus bas). Avec des ondes sinusoidales, ça donne:
B*w*cos(w(t-x/a)+D) +(B*sin(w(t-x/a)+D)+<u>)*B*(w/a)*cos(w(t-x/a)+D) =A*(w/a)*cos(w(t-x/a))/r.
B*w*cos(w(t-x/a)+D) +B²*(w/2a)*sin(2w(t-x/a)+D) +<u>*B*(w/a)*cos(w(t-x/a)+D) =A*(w/a)*cos(w(t-x/a))/r.
B*w*(1+<u>/a)*cos(w(t-x/a)+D) +B²*(w/2a)*sin(2w(t-x/a)+D) =A*(w/a)*cos(w(t-x/a))/r.
B*w*(1+<u>/a)*[cos(w(t-x/a))cos(D)-sin(w(t-x/a))sin(D)] +B²*(w/2a)*sin(2w(t-x/a)+D) =A*(w/a)*cos(w(t-x/a))/r.
w*[B*(1+<u>/a)*cos(D)-A/ar]*cos(w(t-x/a)) -w*B*(1+<u>/a)*sin(D)*sin(w(t-x/a)) +B²*(w/2a)*sin(2w(t-x/a)+D)=0.
<u> est petit devant a:
w*[B*cos(D)-A/ar]*cos(w(t-x/a)) -w*B*sin(D)*sin(w(t-x/a)) +B²*(w/2a)*sin(2w(t-x/a)+D)=0.
B est petit devant a:
w*[B*cos(D)-A/ar]*cos(w(t-x/a)) -w*B*sin(D)*sin(w(t-x/a))=0.
Ceci est valable en tout (x,t) donc [B*cos(D)-A/ar]=0 et B*sin(D)=0, puis sin(D)=0, D=0 ou p, cos(D)=-1 ou 1. Comme A et B sont positifs, cos(D)=0, D=0, puis B=A/ar, B/a=A/a²r. a est la vitesse isotherme du son, a²=sT=P/r, donc B/a=A/P, relation qui lie les amplitudes de vitesse et de pression.
Calculons maintenant le déplacement Dx d'une particule pendant la phase 0-2. Comme u=(dx/dt) particulaire, dx/dt=B*sin(w(t-x/a))+<u>.
x est des 2 cotés de l'équation, c'est ennuyeux mais on peut considérer que pour une particule donnée x/a varie beaucoup moins que t, donc sur un demi-cycle Dx=integrale[0,p/w] 1)=Umax*2R*(p-b)=-4(f/n)*dP/dx, où n est la viscosité cinématique.
d²u/dr²=Umax*(-2(p+b)+12pb*r²), est la dérivée seconde au point r, à t et x constants.
d²u/dr² s'annule quand r²=(p+b)/6pb=(R²+R1²)/6, r(inflexion)=0.58*r(rebroussement).
Calculons la vitesse débitante U. U=integrale(2pr*dr*u®)/pR²=integrale(u®*d(r²))/R², U=Umax*(R²-(p+b)R4 /2+pbR6 /3)/R², U=Umax*(1-(p+b)R²/2+pbR4 /3), U=Umax*(1/2-b/2p+b/3p), U=Umax*(1-b/3p)/2.
Quant au carré moyen de la vitesse <u²> =integrale(u²®*d(r²))/R², <u²>=Umax*integrale2)/R², <u²>=Umax*integrale3)/R², <u²>=Umax*(R²-(p+b)R4 +4), <u²>=Umax*(1/3+(b/p)*(-1/6)+(b/p)²*(1/30)), <u²>=Umax*(1-(b/p)/2+(b/p)²/10)/3.
Le d²u/dr² moyen est <d²u/dr²>=integrale(d²u/dr²*d(r²))/R², <d²u/dr²>=Umax*integrale5)/R², <d²u/dr²>=Umax*(-2(p+b)R²+6pb*R4 )/R², <d²u/dr²>=Umax*(-2(p+b)+6pb*R2 ), <d²u/dr²>=2*Umax*(2b-p).
En notant z=b/p pour simplifier les expressions, il vient: U=Umax*(1-z/3)/2, <u²>=Umax*(1-z/2+z²/10)/3, <d²u/dr²>=2p*Umax*(2z-1).
Voyons maintenant l'équation différentielle à laquelle répond notre profil.
Par conservation de l'impulsion, en tout point du profil d(ru)/dt+dP/dx-µ*d²u/dy²=0, r*(du/dt)+u*(dr/dt)=µ*d²u/dy²-dP/dx, r*(du/dt)-u*(d(ru)/dx)=µ*d²u/dy²-dP/dx, r*(du/dt)-ru*(du/dx)-u²*(dr/dx)=µ*d²u/dy²-dP/dx. En divisant par s et en utilisant P=rsT (où s est la constante massique du gaz parfait (s=8.314/massemolaire J/kg/K), habituellement notée r, interdite dans le cas présent pour cause de double signification), il vient du/dt-u*(du/dx)-u²*(dr/rdx)=n*d²u/dy²-sT*dP/Pdx, où n=µ/rest la viscosité cinématique (m²/s). En approx. isotherme dP=sT*dret dr/rdx=dP/Pdx, on obtient du/dt-u*(du/dx)=n*d²u/dy²+(u²-sT)*dP/Pdx. u² est petit devant sT, ce qui conduit à du/dt-u*(du/dx)=n*d²u/dr²-sT*dP/Pdx en tout point (x,r,t). Ceci est donc valable en moyenne sur toute une section: dU/dt-d<u²>/2dx-n*<d²u/dr²>=-sT*dP/Pdx en tout point (x,t). Puis
Umax*[-dz/6dt-dz/12dx+z*dz/30dx-2pn*(2z-1)]+(d(Umax)/dt)*(1-z/3)/2-(d(Umax)/dx)*(1-z/2+z²/10)/6=-sT*dP/Pdx.
Rappelons nous la relation entre z et dP/dx. Sur la paroi u=0, donc µ*d²u/dr²=dP/dx=µ*Umax*(-2(p+b)+12pb*R²), dP/dx=2pµ*Umax*(5z-1) puis -sT*dP/Pdx=2pn*Umax*(5z-1). En notant L=ln(Umax), il vient [-dz/6dt-dz/12dx+z*dz/30dx-2pn*(2z-1)]+(dL/dt)*(1-z/3)/2-(dL/dx)*(1-z/2+z²/10)/6=2pn*(5z-1)
C'est maintenant que nous allons faire intervenir notre équation de profil; Umax et b dépendent de x et de t, pas de r.
du/dt=d(Umax*(1-p*r²)*(1-b*r²))/dt=(1-p*r²)*[(1-b*r²)*d(Umax)/dt-r²Umax*db/dt].
du/dx=d(Umax*(1-p*r²)*(1-b*r²))/dx=(1-p*r²)*[(1-b*r²)*d(Umax)/dx-r²Umax*db/dx].
d²u/dr²=Umax*(-2(p+b)+12pb*r²).
On ne peut dériver le modèle beaucoup plus car on tombe sur des incohérences, dûes au fait que ce n'est qu'une approximation d'ordre 4.
En prenant un peu de recul sur le modèle on s'aperçoit que les particules centrales et périphériques ne suivent pas le même cycle. Les dT sont différents, et donc le travail produit. Elles échangent du travail par compression et pas seulement par le biais des forces visqueuses, ce que notre modèle ne laisse pas percevoir (le r particulaire restant constant, les vitesses radiales sont nulles).
Le problème est alors trop complexe pour être résolu à la main, il faut utiliser un logiciel de CFD.
Lorsqu'on tente de miniaturiser un clapet antiretour à ressort de rappel, on bute -à l'échelle millimétrique- sur une difficulté de fabrication et une grande sensibilité au colmatage par des poussières.
La présence de hautes températures (>200°C) dégrade les propriétés élastiques des aciers et donc celle d'un éventuel ressort.
Pour ces 3 raisons, on a parfois intérêt à utiliser un clapet statique (sans pièce mobile) dont le montage ci-contre, à effet de réflexion, constitue la solution la plus simple (l'autre solution, dûe à M.Barry, consiste à utiliser un vortex à échappement central et entrée périphérique tangentielle à forte perte de charge pour bloquer un éventuel contre-courant -vortex court-circuité dans le sens de passage normal-, dans un montage appelé diode hydraulique, est difficile à inclure dans une alvéole d'échangeur).
Mode de fonctionnement: en régime stationnaire et quand la pression dynamique est largement inférieure à la pression ambiante, la perte de charge d'élargissement brusque est légèrement inférieure à la perte de charge de rétrécissement brusque (effet de veina contracta) mais surtout
-une onde de compression rebroussante se réfléchit sur la paroi de brusque changement de section, ce qui affaiblit sa propagation.
-si la pression dynamique approche l'ordre de grandeur de la pression ambiante (et donc la vitesse est proche de celle du son) le ralentissement périphérique occasionne une forte compression qui écrase la veine de circulation centrale et amplifie considérablement l'effet de veina contracta.
réalisation, applications: la section de passage doit varier au moins du simple au double. Le montage se prête bien à une réalisation emboutie qui peut être incluse en n'importe quel point d'un échangeur thermique alvéolaire à plaques (figure du bas).
Pour conclure ce paragraphe, nous avons là un accessoire permettant de bloquer d'éventuels contre-courants. On ne l'a pas inclus dans la modélisation thermoacoustique.|
pages francophones
centre acoustique LMFA Thermoacoustique des systèmes frigorifiques miniatures
Générateur thermoacoustique annulaire à ondes progressives (PDF, page 4 du cv de Stéphane Job)
Machines thermoacoustiques (université du Mans)
Simulation numérique d'effets thermoacoustiques (PS, travail de Jean-Yves Tinevez)
Simulation numérique en propagation acoustique (PS, travail de Stanislas Antczak)
Un prototype de compresseur thermoacoustique et modélisation du régénérateur
Prototype de réfrigérateur Stirling thermoacoustique utilisant un agent fluide critique
Instabilité thermoacoustique : modélisation numérique d'un compresseur thermoacoustique (travail de Y Delbende)
pages anglophones
Thermoacoustics une bonne introduction de Ralph Muehleisen: définition, fonctionnement, état de l'art. On y apprend que le rendement exergétique a déjà été poussé à 40% du rendement de Carnot (qui est le maximum théorique). Ici, à isentropics.org, on vise, par modélisation et optimisation systématique, à effleurer les 100%.
Solar powered thermoacoustic refrigerator (STADTAR) montage incluant un moteur et un réfrigérateur thermoacoustiques en série. Performances: côté chaud .457²p/4*1350=221Watts à 475°C(=748K), côté froid 2.5W à 5°C(=278K), côté ambiance à 23°C(=296K). Rendement de carnot correspondant rc=(748-296)/748*296/(296-278)=0.604*16.44=9.93 . Rendement réel r=2.5/221=1.13%. Rendement exergétique re=r/rc=1.14e-3, soit 3.37% de chaque côté.
POWER AND PROPULSION IN THE NEW MILLENIUM Thermo-Acoustic Cycle (TAC) Engines. Tiré du site de Fellows Research Group, Inc .
Tout cela est très ambitieux et prometteur.
Performance measurements on a thermoacoustic refrigerator “Frankenfridge” contient aussi de la théorie. p.29
logiciels
DSTAR (Design Simulation for ThermoAcoustic Research). Le fichier zippé contient entre autres les équations de la thermoacoustique, issues du modèle de Rott.
Simtube un programme de simulation de tube pulsé de Pierre Neveu.
Le <u> sert avec le dxP, à récupérer le travail du cycle. Oublions le pour l'instant. On a alors Dx=B/f, puis dT=dxT*B/f. Le travail du cycle est donc Q=rCp*dT²/(4*T) (en W/m3, voir plus haut) puis Q=rCp*(dxT*B/f)²/(4*T)=<u>*dxP. En utilisant B=Aa/P et a²=sT, on obtient Cp*(dxT*A/f)²/4(g-1)P=<u>*dxP. En posant A=K*a*dxP/f (avec K>1/2p), on trouve Cp*(dxT*K/f²)²/4(g-1)r=<u>/dxP.
Si <u> est fixé (par l'optimisation de l'échangeur), dxP est donc proportionnel à f4 . Bon point pour les hautes fréquences. Mais en fait la fréquence f n'est pas choisie librement. Pour un échangeur donné, elle est conditionnée par l'apport de chaleur en bout d'échangeur, apport qu'à l'idéal on souhaiterait variable à loisir.|
Il serait tentant de modéliser le profil de vitesses par une double parabole (ce qui donne un modèle à 2 paramètres: S1 et (du/dy)paroi), mais un polynôme unique sur toute la section (forcément une quadrique, car la vitesse s'annule 4 fois, ce qui donne un modèle à 5 paramètres) est plus facile à manipuler et offre plus de possibilités. Soient R=f/2 et r=R-y; u® est une fonction paire qui s'annule en r=R et r=R1, qui possède un maximum local en r=0. u®=A*(1-r/R)*(1+r/R)*(1-r/R1)*(1+r/R1), u®=A*(1-r²/R²)*(1-r²/R1²). Incidemment le modèle n'a plus que 2 paramètres.
Le profil parabolique est un cas particulier correspondant à R1=infini.
Si R1=R, (du/dy)paroi =0. C'est la famille de profils qu'on obtient par dégénerescence du profil parabolique lorsque la pression motrice dP/dx a disparu.
Si R1=i*R, u®=A*(1-r4 /R4 ) et (du/dy)paroi est maximal pour A=Umax donné. C'est la configuration d'entrée dans l'alvéole, la plus proche de la distribution uniforme des vitesses.
On peut reformuler ainsi le profil: u®=Umax*(1-p*r²)*(1-b*r²), où p=1/R² est le coefficient de paroi (invariable), b=1/R1² le coefficient de rebroussement et Umax la vitesse centrale. u/Umax=1-(p+b)*r²+pb*r4 , (du/dr)/Umax=-2(p+b)*r+4pb*r3 , (du/dr)=Umax*2r*(2pb*r²-(p+b)) qui s'annule pour r=0 et r²=(p+b)/(2pb)=R²*R1²*(1/R²+1/R1²)/2=(R²+R1²)/2 donc r(rebroussement)~=(R+R1)/2 ce qui correspond à peu près à la figure.
D'autre part le gradient de paroi est du/dy=-(du/dr)R =-Umax*2R*(2pb*R²-(p+b