https://web.archive.org/web/20090623041338/http://mshades.free.fr/thermoacoustique/thermoacoustique.html

Voici à gauche le montage thermoacoustique standard et à droite ce qu'on voudrait réaliser. Le stack (la pile de plaques) avec son fort gradient thermique amplifie les ondes sonores. Dans le montage de droite ce rôle est assumé par les alvéoles de l'échangeur. En fonctionnement réfrigérateur le montage de gauche fournit un écart de 16°C entre source froide et source chaude pour 1000Pa de pression acoustique. Le montage de droite vise un fonctionnement moteur. Les différences principales: à gauche une onde stationnaire et un stack de longueur quart d'onde, à droite une onde progressive avec flux de circulation et une longueur d'échangeur égale à plusieurs dizaines de longueurs d'onde (fréquence acoustique élevée, rendement exergétique proche de 100%). Pour analyser correctement le fonctionnement de ces systèmes il faut descendre au niveau de la particule fluide. c'est ce que font les 2 schémas ci-dessous.

cycles de la particule en diagramme (P,v)

cycles de la particule en diagramme (T,x)

Considérons le cycle thermodynamique de la particule fluide de la manière un peu caricaturale représentée ci-contre en diagramme (T,x), c'est à dire température-position. On a 4 évolutions: 0-1. compression adiabatique (sans échange de chaleur) sans frottement, donc isentropique. 1-2. réchauffe à pression constante. 2-3. détente adiabatique isentropique. 3-0. refroidissement à pression constante. Si les compression et détente adiabatiques sont bien synchronisées avec le déplacement de la particule, celle-ci peut suivre exactement le gradient thermique de la paroi. Ensuite, lors de la réchauffe ou du refroidissement isobare, on continue à longer le gradient avec un léger DT pour assurer le transfert thermique, DT qui peut être aussi petit qu'on veut, typiquement qq dixièmes de degré, en jouant sur l'écartement entre plaques ou le diamètre d'alvéole: on l'a prouvé dans l'optimisation de l'échangeur. En diagramme (T,x) le cycle peut donc être quasiment plat et donc les pertes exergétiques par saut de température presque nulles. Cela ne signifie pas que la puissance motrice du cycle est nulle. Pour le voir il faut passer en diagramme (T,S) ci-dessous. Notons toutefois avant de quitter ce diagramme que comme les phases d'échange thermique ne sont pas exactement (et même pas du tout, voir plus bas) superposées dans l'espace, le flux de chaleur moyen sur un cycle n'est pas constant le long de la paroi, ce qui doit pourtant être le cas dans une configuration d'échangeur, il faut le compenser par de multiples cycles décalés en abcisse x correspondant à autant de particules voisines, autrement dit la longueur de déplacement x2-x0 doit être faible devant la longueur totale L de l'échangeur.

Voici le même cycle en diagramme (T,S), c'est à dire température-entropie. Toute quantité de chaleur infinitésimale dQ reçue par la particule entraîne une variation d'entropie dS=dQ/T, d'où l'on tire dQ=TdS, puis Q=integralecycle (TdS). La surface à l'intérieur du cycle est donc égale à la chaleur Q reçue par la particule pendant le cycle. Soit W le travail reçu, on a W+Q=0 sur un tour (par conservation de l'énergie). Le travail fourni -W est donc égal à la surface du cycle, qu'on cherche à maximiser. Comme T0 et T2 sont déterminés par les positions x0 et x2 sur la paroi qu'on ne veut pas trop écarter pour éviter les frottements, reste à trouver T1 qui maximise Q à T0 et T2 données. Comme dS=CpdT/T-RdP/P, on a dS=Cp*dT/T sur l'évolution 1-2, donc S2-S1=Cpln(T2/T1)=Cpln(T3/T0). De plus Q1-2 =Cp(T2-T1)* et *Q3-0 =Cp(T0-T3). Comme Q=Q1-2 +Q3-0, Q=Cp(T2-T1+T0-T3). Or T2/T1=T3/T0=K. Alors Q=Cp(T1-T0)(K-1), Q=Cp(T1-T0)(T2/T1-1). Dérivons Q par rapport à T1: dQ/d(T1)=Cp[(T2/T1-1)+(T0-T1)T2/T1²]=Cp[T0*T2/T1²-1]. Le maximum de Q est atteint quand T1²=T0T2*, c'est à dire quand T1 est moyenne géométrique de T0 et T2. Alors Q=Cp*T01/2 (T21/2 -T01/2 )(T21/2 /T01/2 *-1), *travail fourni=Q=Cp(T21/2 -T01/2 )2 ** (en J/kg). Quand T0 et T2 sont voisines, posons T0=T et T2=T+dT. On obtient la formule approximative *Q=CpdT²/(4T)*. On voit donc qu'on a intérêt à éviter les dT trop petits. Au passage calculons T3: K=T3/T0=T2/T1=T1/T0 donc *T3=T1, x3=x1*; les zones d'échange thermique [x1,x2] et [x3,x0] ne sont pas du tout superposées. La particule qui dépose la chaleur n'est pas la même que celle qui la prend. Corrigez (mentalement) les 2 diagrammes ci-dessus. Ce travail Q est fourni aux particules voisines sous forme de compression, et il faut absolument l'évacuer du système le plus rapidement possible, sans quoi des phénomènes dissipatifs (ondes de choc,frottements intenses,écarts fluide-paroi…) vont se former, réduisant à néant notre optimisation. On peut récupérer ce travail en organisant dans le stack ou l'alvéole d'échangeur thermoacoustique un gradient de pression *dP/dx*, joint à une circulation moyenne *<u>* dans le même sens, et des clapets anti-retour (à temps de réponse rapide) aux bornes du stack. De cette façon le travail produit est emmagasiné au plus près de la source, dans la particule elle-même, sous forme de compression. Alors fQrf²=f²<u>*dP/dx (en W/m), où *f* est le diamètre de l'alvéole, *f* la fréquence du cycle (en Hz), *r* la masse volumique (en kg/m3). On simplifie cette formule en *fQr=<u>dP/dx (en W/m3). <u> et dP/dx étant déterminés par la chaleur de combustion et l'optimisation de l'échangeur, reste à déterminer la fréquence f. En première approche, le travail fourni à l'abcisse x pendant la phase 0-1, dW/dx=-d(Pu)/dx (en W/m3), est égal à la quantité de chaleur dQ/dx échangée dans la phase isobare 1-2 (valeur appelée compacité à la base du dimensionnement optimisé de l'alvéole, disons 1MW/m3), et la vitesse est de l'ordre de 50cm/s, disons 1m/s. Quand la pression est disons 1bar(=1e5 Pa), Pu varie entre 0 (au point 0) et 1e5 (au point 1). On a donc d(Pu)=1e5, puis dx=1e5/1e6=0.1m=10cm. La phase de compression dure 10cm, la phase de réchauffe 10cm aussi, et pour un gradient de 1500K/m, le dT=T2-T0 est de 300K, donc Q=256 J/kg pour de l'air,Qr=308 J/m3=308Pa. Si dP/dx est de l'ordre de 2000 Pa/m (du même ordre de grandeur que les frottements sans effet thermoacoustique) et <u>=1m/s, fQr=2000W/m3 , d'où f=2000/308=6.5Hz et la longueur d'onde 300/6.5=46m bien trop grande. Transformons la formule fQr=<u>dP/dx. <u>dP/dx=fCpdT²/(4T)r=f(g/4(g-1))P(dT/T)² puis <u>dP/Pdx=f(g/4(g-1))(dT/T)². Il y a une relation entre f,u et dT: dT=dxTu/f à un facteur près. Donc dP/Pdx=udxT²/T²f. Il faut donc maximiser u/f (en m), c'est à dire le déplacement, quand dxT est donné. | |Ceux qui ne connaissent pas très bien la thermodynamique préfèrent utiliser le diagramme (P,V) ci-contre, c'est à dire pression-volume, dans lequel le travail W=integralecycle (-PdV) apparaît directement. Mais comme on travaille à T0 et T2 imposées, son utilisation est moins pratique dans le cas présent, et d'ailleurs maintenant inutile. On peut quand même y lire que le travail est fourni pendant les phases 1-2 et 2-3.

D'autre part le gradient de paroi est du/dy=-(du/dr)_R =-Umax2R(2pbR²-(p+b))=Umax2R(p-b)=-4(f/n)dP/dx, où n est la viscosité cinématique.
d²u/dr²=Umax(-2(p+b)+12pbr²), 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(2prdru®)/pR²=integrale(u®d(r²))/R², U=Umax(R²-(p+b)R⁴ /2+pbR⁶ /3)/R², U=Umax(1-(p+b)R²/2+pbR⁴ /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²>=Umaxintegrale¹⁾/R², <u²>=Umaxintegrale((1-2(p+b)r²+((p+b)²+2pb)r⁴ -2(p+b)pbr⁶ +p²b²r⁸ )d(r²))/R², <u²>=Umax(R²-(p+b)R⁴ +((p+b)²+2pb)R⁶ /3-(p+b)pbR⁸ /2+p²b²R¹⁰ /5)/R², <u²>=Umax(1-(p+b)R² +²⁾, <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²>=Umaxintegrale³⁾/R², <d²u/dr²>=Umax(-2(p+b)R²+6pbR⁴ )/R², <d²u/dr²>=Umax(-2(p+b)+6pbR² ), *<d²u/dr²>=2Umax(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²>=2pUmax(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)=nd²u/dy²-sTdP/Pdx, où n=µ/rest la viscosité cinématique (m²/s). En approx. isotherme dP=sTdret dr/rdx=dP/Pdx, on obtient du/dt-u(du/dx)=nd²u/dy²+(u²-sT)*dP/Pdx. u² est petit devant sT, ce qui conduit à *du/dt-u(du/dx)=nd²u/dr²-sTdP/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²>=-sTdP/Pdx en tout point (x,t). Puis
Umax[-dz/6dt-dz/12dx+zdz/30dx-2pn(2z-1)]+(d(Umax)/dt)(1-z/3)/2-(d(Umax)/dx)(1-z/2+z²/10)/6=-sTdP/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)+12pbR²), *dP/dx=2pµUmax*(5z-1) *puis -sTdP/Pdx=2pnUmax(5z-1). En notant L=ln(Umax), il vient [-dz/6dt-dz/12dx+zdz/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-pr²)(1-br²))/dt=(1-pr²)[(1-br²)d(Umax)/dt-r²Umaxdb/dt].
du/dx=d(Umax(1-pr²)(1-br²))/dx=(1-pr²)[(1-br²)d(Umax)/dx-r²Umaxdb/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).

**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/41350=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)/748296/(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.