différentes à gauche et à droite du stack, correspondant aux températures moyennes $T_h$ et $T_c$, et les dissipations concentrées à l'extrémité froide du résonateur. De nouveau, les conditions limites aux extrémités du stack sont obtenues par raccord avec l'acoustique dans le résonateur (continuité des vitesses et de la pression acoustiques). En imposant que le système global ait une solution non-nulle on obtient une condition d'existence qui représente une relation de dispersion faisant intervenir la pulsation complexe $\tilde{\omega}$. Cette équation est résolue numériquement par une méthode de tir qui nécessite une valeur initiale proche de la solution, valeur qui sera fournie par les résultats numériques précédents. Les diagrammes de stabilité obtenus dans le plan $(T_h, Z)$ ainsi que le calcul des fréquences des modes les plus instables sont en très bon accord avec les résultats de stabilité sur les simulations numériques.

### 3.2.4 Comparaison avec les résultats expérimentaux

Les résultats numériques sont également comparés aux résultats expérimentaux de Atchley et al. [7], sur les conditions de déclenchement de l'instabilité. Le dispositif expérimental consiste en un moteur à onde stationnaire, pour lequel la charge est représentée par les dissipations thermo-visqueues dans le résonateur. La figure 3.8 est une représentation, dans le plan (température chaude $\tilde{T}_h$, pression moyenne $\tilde{P}_m$), des régions qui présentent des régimes de fonctionnement différents : dans la région I, le mode fondamental $\lambda/2$ est observé ; dans la région II, c'est le premier harmonique ; dans la région III, les deux modes sont présents. Les frontières de ces régions sont données par les seuils d'instabilité. On note de plus que si la pression moyenne $\tilde{P}_m$ est très faible, il faut imposer un gradient de température très important pour démarrer le moteur, et c'est le premier harmonique qui apparaît alors. La figure 3.8 montre que pour $\tilde{P}_m = 2.4$

<figure><img src="image_33.png" /><figcaption>FIGURE 3.8 – Diagramme de stabilité de [7]. Courbes marginales des modes 1 et 2. Résultats des calculs numériques représentés par points</figcaption></figure>

bar, lorsque $\tilde{T}_h = 653$ K ($\Delta\tilde{T} = 360$ K), seul le mode fondamental $\lambda/2$ est amplifié. Et lorsque $\tilde{T}_h = 730$ K ($\Delta\tilde{T} = 437$ K), le premier harmonique devient à son tour instable.

Pour comparer aves les expériences, plusieurs simulations ont été faites pour différentes pressions moyennes. Comme les dissipations ne se font pas de la même manière entre simulations et expériences, il est nécessaire de caler les deux. La démarche suivie a été de fixer la pression moyenne et de trouver par simulation la valeur de la charge $Z$ qui correspond à la température critique du diagramme de stabilité. Par exemple, pour $\tilde{P}_m = 4.4$ bar on trouve que la valeur de la charge $f = 45$, soit $Z = -0.928$ correspond à la température critique $\tilde{T}_h = 660$ K pour laquelle le mode fondamental $\lambda/2$ est le plus instable ($\Delta\tilde{T} = 367$ K). Cette valeur correspond à une résistance dimensionnelle $\tilde{R} = f M \tilde{P}_m / (A U_{ref}) = 17.1$ MPa.s/m$^3$, $A$ étant l'aire de la section du résonateur ($\tilde{R} = \gamma \frac{1 - Z M \tilde{P}_m}{1 + Z A \tilde{U}_{ref}}$). Si on augmente la température jusqu'à
$\tilde{T}_h = 821$ K ($\Delta\tilde{T} = 528$ K), avec la même charge, le premier harmonique $\lambda$ devient également instable, en accord avec la valeur $\tilde{T}_h = 800$ K de l'expérience. Les résultats numériques sont représentés avec des points bleus et verts sur la diagramme d'Atchley. Les signaux de pression associés aux calculs sont montrés à droite du diagramme. On a également vérifié que pour des faibles valeurs de pression moyenne, $\tilde{P}_m = 0.15$ bar on trouve que le premier harmonique est le mode le plus instable et la valeur de la charge $f = 16$, soit $Z = -0.839$, correspond à la température critique $\tilde{T}_h = 718$ K ($\Delta\tilde{T} = 367$ K). Cette valeur correspond à une résistance dimensionnelle $\tilde{R} = 2.1$ MPa.s/m$^3$. Si on essaye de faire une interprétation physique de cette charge, on trouve qu'elle pourrait jouer le rôle d'une résistance thermique. Une estimation de la résistance thermique de l'expérience d'Atchley, basée sur l'expression classique ([78])

$\tilde{R}_{\kappa} = \frac{2\gamma\tilde{P}_m}{\tilde{\omega}(\gamma - 1)\tilde{S}_s\tilde{\delta}_{\kappa}}$, où $\tilde{S}_s = \pi\tilde{d}\tilde{L}_{res}$ est la surface latérale du résonateur, conduit à un bon accord qualitatif avec les expériences (cf. figure 12 dans [45], article en annexe). Ceci donne une indication sur la manière dont on peut rendre compte du fonctionnement d'une machine réelle où les dissipations sont réparties dans tout le résonateur, à l'aide du modèle hybride avec une charge résistive localisée.

## 3.3 Analyse du régime périodique

Le modèle numérique décrit précédemment a été utilisé pour étudier le comportement transitoire d'un moteur chargé : influence de la charge, de l'espacement entre le stack et les échangeurs, de la pression moyenne ainsi que de la différence de température entre les échangeurs $\Delta T$. Pour des différences de température $\Delta T > \Delta T_{crit}$ on observe une onde avec un comportement multi fréquentiel, résultat confirmé par les études de stabilité de Atchley [7]. Un comportement similaire a été observé aussi dans les expériences sur le dispositif expérimental du LIMSI pour des drive ratio supérieurs à 1% (voir mesures expérimentales chapitre 2). Ceci avait été attribué à la propagation acoustique non-linéaire dans le résonateur. Il y a probablement une combinaison des effets non linéaires : ceux intrinséques à la cellule active et ceux liés au résonateur. L'analyse du transitoire montre aussi qu'il y a une redistribution des températures dans le stack : le gradient de température diminue lentement pendant le transitoire ayant comme conséquence un réajustement de l'amplitude de saturation après la phase rapide d'amplification, traduit dans les évolution temporelles du signal de pression par un "overshoot". La saturation se fait sur une échelle de temps bien plus longue que l'amplification initiale (10000 périodes = 100hCPU pour la saturation et 400 périodes = 2hCPU pour l'amplification), et elle a lieu à des vitesses du fluide très élevées, ce qui nécessite le réajustement du pas de temps, par exemple de $10^{-2}$ à $10^{-4}$. Pour des valeurs importantes de $\Delta T$ il est difficile d'atteindre le régime stable (sans évolution de l'amplitude de pression) car les simulations deviennent très coûteuses en temps de calcul. De ce fait, une étude paramétrique n'a pas pu être envisagée et probablement ce modèle n'est pas adapté pour cet objectif. Pour avoir une saturation à des niveaux raisonnables du drive ratio il faudrait pouvoir prendre en compte la propagation non linéaire dans le résonateur et donc changer le modèle concernant l'acoustique du résonateur. Néanmoins, le modèle qu'on a développé montre le caractère non linéaire de l'écoulement dans la cellule active, responsable de la saturation. La charge quand à elle va avoir une influence sur le $\Delta T_{crit}$ et sur la sélection des modes instables. En régime périodique, on a pu calculer le drive ratio du moteur, la puissance acoustique délivrée dans le résonateur froid, et on a pu faire un bilan énergétique pour caractériser le moteur simulé. Une analyse de la dynamique de l'écoulement au voisinage des extrémités du stack et des échangeurs pour différents drive ratio a été également faite en régime périodique, mettant en évidence les non-linéarités hydrodynamiques dans la cellule ac-
tive. Le calcul des champs instantanés de vorticité montre que pendant une période acoustique, une paire de tourbillons se forme en sortie de l'échangeur froid, qui sont ensuite aspirés entre les plaques du stack, suivi d'une réformation de deux tourbillons contrarotatifs en fin de cycle. Le champ de température subi des fluctuations importantes.

<figure><img src="image_34.png" /><figcaption>FIGURE 3.9 – Champ instantané de vorticité à la sortie de l'échangeur froid pendant une période acoustique à 5 instants de la période. Config. Atchley. $Dr = 3.4\%$..</figcaption></figure>

## 3.4 Application : thermophones

Les compétences développées dans la modélisation et simulation des moteurs thermoacoustiques ont nourri une collaboration avec J. Rémus (Ipotam Mécamusique) sur une application musicale dans le cadre d'un projet Art&Sciences, "Thermophonia", en 2015/2016 (collab. J. Rémus - Ipotam Mécamusique, C. Weisman, V. Daru, F. Jebali). Jacques Rémus, musicien et plasticien, constructeur de sculptures et machines musicales, a créé depuis plusieurs années des tuyaux sonores thermo-acoustiques, qu'il a baptisés "thermophones" (Figure 3.10). Le projet Thermophonia avait comme objectif de réaliser un ensemble visuel et sonore de jeux de tuyaux de thermophones électriques, afin de créer des installations sonores spatialisées.
<figure><img src="image_35.png" /><figcaption>FIGURE 3.10 – Thermophones aux support mobiles. Installation Curiositas 2015, J. Rémus</figcaption></figure>

Le thermophone est un instrument de musique particulier, où le son est produit par effet thermoacoustique. Il est constitué d'un tuyau en acier doux, verre ou aluminium, de section circulaire ou carrée, ouvert aux deux extrémités (flûte) ou bien ouvert à une extrémité et fermé à l'autre extrémité (bourdon). À l'intérieur du tuyau se trouve un stack solide. Le stack est chauffé par une résistance électrique à une extrémité. Lorsque la température de l'extrémité chaude dévient suffisamment élevée, le thermophone se met à chanter. Le démarrage du son est le résultat d'une instabilité thermoacoustique. Un thermophone s'apparente donc à un générateur d'ondes thermoacoustique, à onde stationnaire, le gaz utilisé étant l'air atmosphérique. Pour des raisons pratiques, les thermophones ne sont pas pourvus d'un échangeur froid permettant de contrôler la température de l'autre extrémité du stack, comme les moteurs thermoacoustiques, et celle-ci peut donc varier le long du tube (voir figure 3.11) : Le son du thermophone est dominé par la

<figure><img src="image_36.png" /><figcaption>FIGURE 3.11 – Schéma : comparaison moteur thermoacoustique/thermophone</figcaption></figure>

fréquence de résonance fondamentale du tuyau, liée à sa longueur, et il dépend bien évidem-
ment de la chaleur apportée et de la géométrie du stack et du résonateur. Il est donc important
de connaître les paramètres de fonctionnement du thermophone pour avoir le son désiré pour un
ensemble des tuyaux. Nous avons utilisé un modèle 1D d'un thermophone, basé sur la théorie
linéaire de Rott que l'on a adapté pour prendre en compte les particularités du thermophone
(voir section 3.2). L'équation de dispersion a été résolue numériquement sur Matlab pour une
variété de thermophones existants (tuyaux de section carrée ou circulaire, ouvert/ouvert ou ou-
vert/fermé). Pour chaque instrument on a determiné par une étude de stabilité le seuil de l'écart
de température entre les extrémités du stack correspondant au démarrage, ainsi que la fréquence
du mode le plus instable. Plusieurs études paramétriques ont été menées sur la solution numé-
rique afin de guider la conception de futurs thermophones : position du stack, longueur du stack,
température de la résistance chauffante (stage M1 T.Hung, 2015). Les résultats ont fait l'objet
d'une communication au CFA 2016.
# Chapitre 4

## Couplages fluide-poreux

La modélisation de la cellule active où se réalise la conversion d'énergie est une étape importante dans la modélisation de la machine complète, comme on a pu le voir dans le chapitre précédent. L'élément clé de la cellule active est un milieu poreux, constitué d'un empilement de plaques ou de grilles métalliques, ou autre géométries [78]. L'approche actuelle de modélisation pour la prise en compte de la géométrie des pores est à travers des fonctions de forme et des approximations de couche limite. Ces modèles sont adaptés à des géométries simples de milieux poreux comme un empilement de plaques (stack) ou un réseau de capillaires parallèles de section circulaire, rectangulaire ou arbitraire. Dans le cas d'un empilement de grilles, utilisé pour des machines plus efficaces, il s'agit d'un milieu poreux topologiquement complexe et difficile à modéliser. Plusieurs études sur des modèles de milieux poreux adaptés aux régénérateurs des systèmes thermoacoustiques ont été menés récemment [42, 68, 79].

### Travaux de thèse : modélisation du transport d'un traceur, par convection et diffusion, dans les milieux poreux et multiporeux

Dans mes travaux de thèse, faite sous la direction de T. Lévy, j'ai utilisé la méthode d'homogénéisation des milieux périodiques (développée par Sanchez-Palencia [71]) pour étudier le transport de masse à travers un milieu poreux périodique. L'homogénéisation est une méthode qui emploie des techniques asymptotiques conjointes aux techniques d'équations aux dérivées partielles pour modéliser le comportement des milieux comportant beaucoup d'hétérogénéités de très petites dimensions, réparties avec une certaine régularité, de façon qu'à une échelle bien plus grande que celles des hétérogénéités, le milieu puisse être perçu comme homogène. En partant des équations écrites à l'échelle du pore et en développant les inconnues en échelles spatiales doubles (échelle du pore et échelle du milieu), avec des approximations successives et un processus de moyenne spatiale sur un volume élémentaire représentatif (ici la période du milieu) on déduit les équations macroscopiques ou moyennes. Ces équations font intervenir des quantités moyennes et des coefficients "macroscopiques" qui dépendent de la géométrie du pore.

Le premier problème étudié a été celui du transport d'un traceur passif par convection-diffusion, dans un milieu poreux périodique, dans deux cas : le cas où la convection et la diffusion sont du même ordre de grandeur au niveau macroscopique et le cas où les deux phénomènes sont du même ordre de grandeur au niveau microscopique (au niveau du pore). L'équation macroscopique vérifiée par la concentration moyenne du traceur a été obtenue par homogénéisation, ainsi que les coefficients macroscopiques de diffusion. Quand les deux modes de transport sont du même ordre de grandeur au niveau microscopique, l'analyse asymptotique a été faite avec deux échelles de temps différentes et on a montré que l'utilisation d'une seule échelle temporelle est suffisante, le temps caractéristique important étant le temps de convection
macroscopique.

Le deuxième problème étudié a été le problème de convection-diffusion dans des milieux multiporeux, dans le cas où la convection et la diffusion sont du même ordre de grandeur au niveau macroscopique. Dans le cas des milieux poreux contenant des inclusions imperméables, une nouvelle équation macroscopique de convection-diffusion a été obtenue, avec un tenseur de diffusion ne dépendant pas de la vitesse de filtration, pour une concentration volumique en inclusions rigides d'ordre 1. Dans le cas des faibles concentrations, les coefficients homogénéisés ont été approximés avec un développement en fonction de la concentration des inclusions solides pour des géométries spécifiques des inclusions [10]. Pour finir, le problème de convection-diffusion dans des milieux poreux comportant des cavités saturées de fluide a été étudié dans deux cas : cavités non connexes et cavités connexes. Cette partie de modélisation a été complétée par le calcul numérique des coefficients de diffusion homogénéisés macroscopiques. Pour cela il faut résoudre numériquement les problèmes locaux sur la cellule de base. La résolution numérique a été faite avec la méthode des éléments finis, en utilisant deux codes de calcul : MODULEF en 2D et CASTEM 2000 en 3D, dans le cas isotrope, pour plusieurs géométries. À partir de ces solutions on a calculé numériquement les coefficients homogénéisés en fonction de la porosité du milieu. Les résultats numériques ont été validés avec les résultats analytiques trouvés pour des inclusions sphériques en faible concentration.

Un des objectifs que j'avais en arrivant au LIMSI était d'établir un nouveau modèle du transport ayant lieu au sein du régénérateur, modèle qui devait prendre en compte explicitement le caractère poreux de celui-ci et les non linéarités locales qui y surviennent. L'investissement important que j'ai fourni pour la recherche expérimentale lors de mes premières années au LIMSI m'ont éloignée de cet objectif. Néanmoins, mon expérience de modélisation en milieu poreux a été utilisée pour les études de convection dans le stack soumis à un fort gradient de température.

## 4.1 Étude numérique de la convection naturelle autour et dans le stack moteur

Dans l'étude des machines thermacoustiques, les écoulements de convection naturelle ont été généralement négligés dans les modèles 1D des systèmes thermoacoustiques. Ils existent pourtant à la fois dans la cellule active et dans la résonateur. Les essais infructueux de démarrage du moteur thermoacoustique d'un grand démonstrateur construit par les équipes LIMSI/INPO en 2005-2006 ont motivé un travail plus fondamental sur les écoulements de convection naturelle présents dans le stack moteur. Le grand démonstrateur thermoacoustique a été rempli avec de l'air, azote ou hélium, sous une pression moyenne allant de 1 à 30 bars. Malgré une différence de température aux échangeurs $\Delta T > \Delta T_{crit}$, le "bon champ" de température nécessaire au déclenchement de l'instabilité n'était pas obtenu. En effet, l'analyse des mesures de température à l'intérieur du stack montrait une inhomogénéité dans le plan de la gravité, avec des températures plus chaudes en haut et plus froides en bas, à de grandes valeurs de pression moyenne. Ces expériences indiquaient la présence d'écoulements de convection naturelle, comme en cavité différentiellement chauffée.

Je me suis donc intéressée à la convection naturelle dans la cellule active, modélisée comme une cavité 2D différentiellement chauffée et partiellement remplie d'un empilement de plaques planes conductrices horizontales régulièrement espacées, représentant un stack moteur placé entre deux échangeurs de chaleur. Ce travail a été effectué en collaboration avec P. Duthil, C. Weisman et P. Le Quéré. La géométrie considérée est évidemment une simplification par rapport à la géométrie réelle du stack qui remplit un tube cylindrique : cavité rectangulaire 2D, remplie partiellement par un empilement de plaques parallèles (4.1). Le rapport de forme de
<figure><img src="image_37.png" /><figcaption>FIGURE 4.1 – Schéma de la cavité rectangulaire, partiellement remplie d'un empilement de plaques planes horizontales.</figcaption></figure>

la cavité L/H (longueur/hauteur) est fixé pour cette application par le dimensionnement du moteur thermoacoustique (= 2.5 dans notre cas d'étude). De même, les échangeurs sont ici représentés par une simple condition limite en température sur les parois verticales de la cavité, ce qui représente également une simplification (on appelle $\Delta T = T_h - T_c$ la différence de température imposée entre les parois, $T_h$ = température de l'échangeur chaud, $T_c$ = température de l'échangeur froid). L'empilement de plaques est situé à une distance l de chaque paroi verticale représentant l'échangeur de chaleur. L'espacement inter plaques d et l'épaisseur des plaques e définissent la porosité $\varepsilon$ du stack, $\varepsilon = d/(d + \varepsilon)$. Les conditions d'adiabaticité sur les parois horizontales sont une bonne représentation de ce que l'on souhaite obtenir dans la pratique en isolant le tube du mieux possible. Aux interfaces solide/fluide on suppose la continuité des températures et des flux de chaleur ainsi que l'adhérence pour le champ de vitesses. Les stacks utilisés dans l'application sont constitués d'une centaine de plaques, et les nombres de Rayleigh $Ra = g\beta\Delta T H^3/(\nu\alpha)$ (basés sur la hauteur et les propriétés physiques du fluide) peuvent varier entre $10^5$ et $10^8$. Les paramètres $\beta = 1/T_m$ (avec $T_m = (T_h + T_c)/2$), $\alpha$, $\nu$ et $\rho$ tous pris constants. Afin d'éviter les études paramétriques systématiques sur la géométrie finale, ce qui nécessiterait un maillage très fin et des calculs lourds, l'étude est réalisée pour une configuration fixée, et pour un nombre de plaques variable (mais inférieur à la réalité) : 17, 34, 68, 100. Les calculs ont été réalisés en conservant la porosité $\varepsilon$ constante. L'objectif était de détermier le comportement global de l'écoulement, décrit comme un écoulement en milieu poxe, pour prédire ensuite la solution dans un cas réaliste. Dans le domaine fluide on considère les équations de Navier-Stokes, écrites dans l'approximation de Boussinesq :

$$ 
\begin{aligned}
\nabla \cdot \mathbf{V} &= 0 \\
\frac{\partial (\mathbf{V})}{\partial t} + \nabla \cdot (\mathbf{V} \otimes \mathbf{V}) &= -\nabla p + \nu \Delta \mathbf{V} + g\beta(T - T_m) \\
\frac{\partial T}{\partial t} + (\mathbf{V} \cdot \nabla)T &= \kappa \Delta T
\end{aligned}
$$

Dans les plaques du stack la température vérifie l'équation de conduction stationnaire :

$$ \frac{\partial T}{\partial t} = \alpha_s \Delta T \qquad (4.1) $$

avec $\alpha_s$ la diffusivité thermique dans le solide. Les équations sont résolues numériquement avec un code volumes finis (développé par P. Le Quéré), utilisant une discrétisation spatiale de second ordre. Les termes de diffusion sont traités implicitement et les termes de convection explicitement. Le code inclut un algorithme de projection avec un solveur ADI et des méthodes multi-grilles. Les équations sont résolues simultanément dans le fluide et dans le solide
<figure><img src="image_38.png" /><figcaption>FIGURE 4.2 – Isovaleurs des températures dans le stack avec 17 plaques, pour a) $Ra = 10^5$ et b) $Ra = 10^7$.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>FIGURE 4.3 – Champ de vitesses horizontales pour $Ra = 10^7$.</figcaption></figure>

en utilisant une variable de phase. Les équations sont intégrées jusqu'à l'obtention d'un état stationnaire. Les grilles spatiales sont uniformes dans les deux directions. On a montré que la convection modifie considérablement la distribution de température à l'intérieur du stack.

Les calculs réalisés ont montré qu'en fonction de la distance entre les plaques, le phénomène global de convection s'accompagne éventuellement de transferts convectifs de plus petite échelle localisés entre deux plaques (figure 4.3). Plusieurs modèles de type milieu-poreux ont été testés pour proposer une formulation analytique approchée de l'écoulement du cœur de la cavité, montrant la nécessité d'utiliser une description anisotrope. Le modèle retenu est celui de la convection naturelle stationnaire en milieu poreux vérifiant la loi de Darcy, dans l'hypothèse de Boussinesq [13] :

$$ \nabla \cdot \vec{V} = 0 \qquad (4.2) $$

$$ \mu \vec{V} = -\overline{\overline{K}}(\nabla p - \rho_0 g \beta (T - T_m) \vec{e}_z) \qquad (4.3) $$

$$ \vec{V} \cdot \nabla T = \nabla \cdot (\overline{\overline{\alpha}}_m \nabla T) \qquad (4.4) $$

$\overline{\overline{K}}$ étant le tenseur de perméabilité et $\overline{\overline{\alpha}}_m$ est le tenseur de diffusivité thermique effective du milieu poreux. Pour un milieu poreux homogène et anisotrope, on a :

$$ \overline{\overline{K}} = \begin{bmatrix} K_{xx} & 0 \\ 0 & K_{zz} \end{bmatrix} \qquad \overline{\overline{\alpha}}_m = \begin{bmatrix} \alpha_{xx} & 0 \\ 0 & \alpha_{zz} \end{bmatrix} \qquad (4.5) $$
avec

$$K_{xx} = \frac{d^3}{12(d+e)}, \quad K_{zz} = 0, \quad \alpha_{xx} = \frac{\lambda_{xx}}{\rho_0 c_p}, \quad \alpha_{zz} = \frac{\lambda_{zz}}{\rho_0 c_p} \qquad (4.6)$$

avec $c_p$ la chaleur spécifique du fluide, et

$$\lambda_{xx} = \epsilon\lambda + (1 - \epsilon)\lambda_s \quad \text{et} \quad \frac{1}{\lambda_{zz}} = \frac{\epsilon}{\lambda} + \frac{(1 - \epsilon)}{\lambda_s} \qquad (4.7)$$

où $\lambda$ est la conductivité thermique du fluide, $\lambda_s$ est la conductivité thermique du solide.

Avec le nombre de Rayleigh modifié $Ra_K = g\beta\Delta THK_{xx}/\nu\alpha_{zz}$ et une vitesse de référence $U_{refK} = \frac{\alpha_{zz}}{L}$, une solution approchée 1D adimensionnée ($u^*, T^*$) dans le cœur de la cavité est retenue :

$$u^* = -C(z^* - \frac{1}{2}) \qquad (4.8)$$

$$T^* = Cx^* + D + \left(\frac{H}{L}\right)^2 C^2 Ra_K \left(\frac{1}{4}z^{*2} - \frac{1}{6}z^{*3}\right) \qquad (4.9)$$

avec $C$ et $D$ deux constantes qui dépendent du nombre de Rayleigh et du rapport $L/H$ ([13]). Dans notre étude on a déterminé ces constantes à partir des solutions numériques. L'analyse a montré que l'approche milieu poreux fournit un ordre de grandeur correct pour l'échelle de vitesse et permet une description qualitative des solutions. La solution approche bien la solution numérique lorsque le rapport de forme de la cavité est suffisamment grand $L/H = 10$. Dans notre application, le rapport de forme n'étant pas très grand, il n'y a pas de solution analytique de type milieu poreux qui puisse décrire correctement le phénomène. Néanmoins, l'approche milieu poreux est utile pour donner une idée de l'écoulement et de la distribution de tempéramenture dans le cœur de la cellule. On utilisera la représentation du stack comme milieu poreux homogène anisotrope dans l'étude de la section suivante.

## 4.2 Étude numérique de la convection naturelle dans un système thermoacoustique

Les études de convection naturelle, abandonnées en 2008, ont été repris en 2018 sous l'impuissance de l'ANR TACOT. Le travail a été mené en collaboration avec O. Hireche (post-doc pendant 1 an), C. Weisman, V. Daru et Y. Fraigneau. Les résultats numériques ont montré un très bon accord avec les expériences effectuées au laboratoire Pprime, par H. Bailliet et I. Radaman.

L'objectif de ce travail a été de modéliser et simuler en 3D les écoulements de convection naturelle qui peuvent se développer dans des systèmes thermoacoustiques constitués d'une cellule active soumise à une différence de température via deux échangeurs de chaleur, le tout placé dans un résonateur rempli de gaz. Les simulations ont été faites sur des géométries simplifiées spécifiques à deux configurations réalistes construites au laboratoire P Prime et sur lesquelles l'équipe de H. Bailliet a fait des mesures de température et de vitesse par PIV et LDV. Les expériences et les simulations sont effectuées en l'absence d'onde acoustique. Les deux approches ont été élaborés en même temps, en faisant des ajustements au fur et à mesure, dans le but de rapprocher le plus possible les conditions de la simulation et des expériences. Il est important de préciser que la partie chaude du résonateur et celle de la cellule active ne sont pas accessibles aux mesures de vitesse. Les mesures de température ne sont pas nombreuses non plus. De ce fait, la simulation numérique seule peut nous donner des informations sur les écoulements dans cette partie importante du système, dès lors que modèle est validé. L'objectif principal de ce
travail a été de montrer la capacité de prédiction des modèles classiques de convection naturelle en cavité fluide, partiellement remplie de milieu poreux dans des configurations représentatives des systèmes thermoacoustiques.

Deux configurations disposées horizontalement ont été modélisées. La première consiste en un résonateur rempli d'air, contenant un stack en céramique, chauffé à une extrémité par une résistance électrique jouant le rôle d'échangeur chaud. La deuxième consiste en un résonateur de section rectangulaire, rempli d'air, contenant une cellule active constituée d'un régénérateur (empilement de grilles métalliques) placé entre une résistance électrique (échangeur chaud) et un échangeur à température ambiante (échangeur froid). Celui-ci utilisait l'eau comme fluide secondaire refroidissant et des ailettes pour un meilleur échange entre la partie solide de l'échangeur et le gaz (voir [40]).

Le modèle numérique consiste à résoudre les équations de Navier-Stokes avec l'approximation de Boussinesq dans les cavités fluides et une équation de type Darcy-Brinkman dans la cellule active [12, 13, 34], dans chacune des configurations. À noter que dans la cellule active, dans le cas d'un moteur, la différence de températures peut être très grande et le critère d'application pour l'approximation de Boussinesq n'est plus vraiment respecté ($\Delta T/T_m = O(10^{-1})$ ce qui n'est pas très petit). Néanmois, ce modèle peut servir comme premier choix d'étude. Dans les applications thermoacoustiques, la cellule active est un milieu poreux avec une porosité très grande (typiquement 80%) et l'écoulement qui la traverse est très lent (le nombre de Reynolds du pore est très inférieur à 1). Par conséquent, l'écoulement à travers le milieu poreux peut être modélisé avec l'équation de Darcy-Brinkman où la viscosité effective dans le terme de Brinkman $\mu_{eff}\Delta V$ est prise $\mu_{eff} \approx \mu$.

$$ 
\begin{gathered}
\nabla \cdot \mathbf{V} = 0 \\
\frac{1}{\epsilon} \frac{\partial (\rho \mathbf{V})}{\partial t} + X \frac{1}{\epsilon^2} \nabla \cdot (\rho \mathbf{V} \otimes \mathbf{V}) = -\nabla p + \mu \Delta \mathbf{V} - \rho \beta (T - T_{ref}) \mathbf{g} - (1-X) \mu \mathbf{K}^{-1} \mathbf{V} \\
((1-X)(\rho c_p)_{eff} + X \rho c_p) \frac{\partial T}{\partial t} + \rho_0 c_p (\mathbf{V} \cdot \nabla) T = \nabla \cdot [(X k + (1-X) k_p) \nabla T]
\end{gathered}
$$

Dans les équations précédentes la masse volumique $\rho$, la conductivité thermique $k$, la capacité calorifique $c_p$, la viscosité dynamique $\mu$ sont constantes et égales aux valeurs prises à $T_{ref}$, $\beta = 1/T_{ref}$ est le coefficient de dilatation isobare pour un gaz parfait, $(\rho c_p)_{eff} = \epsilon(\rho c_p)_{fluide} + (1-\epsilon)(\rho c_p)_{solide}$ est la valeur effective pour le milieu poreux, $K$ est le tenseur de perméabilité, $X$ est une fonction de phase ($X = 0$ dans le poreux et $X = 1$ dans le fluide) et $k_p$ la conductivité effective du milieux poreux. Les caractéristiques du milieu poreux sont basées sur [13] : la perméabilité axiale $K_{xx} = \frac{a^3 \epsilon}{8}$, $K_{yy} = K_{zz} = 0$ pour le stack anisotrope ($a$= taille du pore, $\epsilon$= porosité), respectivement $K_{xx} = K_{yy} = K_{zz} = 10^{-8}$ pour le régénérateur considéré isotrope et la conductivité thermique effective $k_p = \epsilon k_{fluide} + (1-\epsilon)k_{solide}$.

Dans la première configuration le stack en céramique est modélisé comme un milieu poreux homogène anisotrope. L'échangeur chaud a été modélisé avec plusieurs corps solides à température fixée. L'échangeur froid étant inexistant dans cette expérience, l'air du résonateur joue ce rôle. Dans la simulation il n'y a donc pas d'échangeur froid et la température est imposée' au long de la paroi du résonateur en utilisant les mesures expérimentales. Dans la deuxième configuration, l'échangeur chaud est le même, le régénérateur est modélisé avec un milieu poreux homogène isotrope tandis que l'échangeur froid est modélisé avec des barres solides à température fixée (température de l'eau), imergées dans un milieux poreux anisotrope (voir figure 4.4). Les conditions en température sur les frontières du résonateur ont été imposées de la façon suivante : adiabatique à gauche de l'échangeur chaud, température constante égale à la température ambiante à droite de l'échangeur froid et variation linéaire de température entre les deux échangeurs.
<figure><img src="image_40.png" /><figcaption>FIGURE 4.4 – Gauche : Configuration numérique pour l'expérience 2. Droite : échangeur froid du dispositif expérimental</figcaption></figure>

Les simulations numériques ont été effectuées avec le code SunfluidH, développé au LIMSI par Y. Fraigneau pour des simulations d'écoulements incompressibles. Les équations (4.10) sont résolues par une méthode volumes finis, avec une discrétisation de second ordre en espace et en temps. Les termes de diffusion sont traités implicitement et les termes de convection explicitement. Le code inclut un algorithme de projection avec un solveur ADI et des méthodes multi-grilles. Le maillage est uniforme dans toutes les directions. Le pas de temps est ajusté durant la simulation pour satisfaire le critère de stabilité CFL < 0.5.

## 4.2.1 Résultats sur la configuration 1

Les résultats des simulations sont en très bon accord avec les mesures expérimentales des vitesses (PIV/LDV). L'analyse des résultats montre la présence d'un écoulement de convection dans la cavité fluide à droite du milieu poreux. Cette cavité est différentiellement chauffée (chaude à gauche, en sortie du stack poreux $T = 310K$) et froide à droite (l'extrémité droite du résonateur $T = 290K$). L'air monte du côté chaud et descend du côté froid avec un écoulement global du milieu poreux vers l'extrémité droite du résonateur dans la partie haute et de l'extrémité du résonateur vers le milieu poreux dans la partie basse. La vitesse près du milieu poreux est quasi nulle avec un profil de couche limite pour la composante verticale. Les différents profils axiaux, à différentes positions radiales montrent le caractère tridimensionnel des écoulements. L'ordre de grandeur du champ de vitesses est de 0.04 m/s. A l'intérieur du stack, un poreux fortement anisotrope, l'écoulement est très faible (vitesse de 0.04 m/s) malgré la différence de température très grande, $\Delta T = 220K$, et suit la géométrie de l'échangeur chaud. Le faible écoulement obtenu est en accord avec l'estimation du nombre de Rayleigh pour un milieu poreux, $Ra_p = \frac{\rho g \beta \Delta T H_{xx}}{\mu \alpha_p} = 1.1$ dans ce cas ($\alpha_p = \frac{k_p}{\rho c_p}$).

Dans la cavité fluide à gauche de l'échangeur chaud, le sens du $\Delta T$ est inversé. La cavité est froide à gauche et chaude à droite. La différence de température est très importante, $\Delta T = 110K$ et comme c'est une cavité fluide, l'écoulement de convection est très fort, avec des vitesses de 0.15m/s. Le calcul du nombre de Rayleigh fluide donne ici $Ra = \frac{\rho g \beta \Delta T_l H^3}{\mu \alpha} = 3.36 \times 10^5$, $\alpha = \frac{k}{\rho c_p}$. Ce calcul montre la nécessité d'isolation de cette partie chaude dans le but de réduire la différence de température et diminuer donc l'écoulement de convection.

Les mesures effectuées sur ce système thermoacoustique avec une onde entretenue par la même puissance électrique fournie à l'exchangeur chaud, ont montré l'existence d'un écoulement de streaming dans la partie froide du résonateur de 1cm/s, superposé au champ acoustique d'amplitude 4m/s. L'écoulement de streaming est donc du même ordre de grandeur que l'écoulement de convection en absence de l'onde. Il sera donc important de prendre en compte la convection naturelle pour comprendre l'origine et le comportement des écoulements continus
présents dans les machines thermoacoustiques.

<figure><img src="image_41.png" /><figcaption>FIGURE 4.5 – Vecteurs vitesse (flèches) et isolignes de la composante axiale de l’écoulement pour (a) les expériences et (b) les simulations numériques, sur un plan vertical central, à droite du milieu poreux.</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>FIGURE 4.6 – (a) Profil vertical de la composante axiale de vitesse, à droite du milieu poreux, dans le plan vertical à $x = 40mm$. (b) Profil de la composante axiale de vitesse en fonction de la coordonnée radiale, dans le plan horizontal à $x = 60mm$.</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>FIGURE 4.7 – Vecteurs vitesse (gauche) et champ de température (droite) dans le plan vertical central, à gauche de l’échangeur chaud</figcaption></figure>
4.2.2 Résultats sur la configuration 2

Dans cette configuration, la présence de l'échangeur froid change complètement la topologie des écoulements convectifs : tout d'abord, dans la partie froide de l'échangeur nous avons observé un écoulement de convection inversé par rapport à la configuration 1 sans échangeur : l'air monte du côté de l'extrémité du résonateur et descend du côté du régénérateur. Une analyse plus attentive du champ de température montre que la différence de température est inversée $\Delta T = T_{\text{cold}} - T_{\text{amb}} = -3K$. L'ordre de grandeur du champ de vitesse convective est 0.025 m/s, en accord avec l'estimation du nombre de Rayleigh de la cavité fluide $Ra = 1.63 \times 10^4$. L'échangeur froid empêche donc le développement d'un fort écoulement convectif. Dans le régénérateur, l'écoulement qui se développe est similaire à un écoulement de convection naturelle dans un milieu poreux, soumis à une différence de température $\Delta T = T_{\text{hot}} - T_{\text{cold}} = 147K$. La vitesse de 0.002 m/s est en accord avec l'estimation du nombre de Rayleigh poreux $Ra_p = 1.26 \times 10^{-1}$. La convection dans le régénérateur est freinée par sa faible perméabilité. Par conséquent, l'écoulement de convection est bien faible et le champ de température a une stratification presque horizontale, comme dans un solide conducteur. Du côté gauche de l'échangeur chaud, le résonateur a été isolé et la température est presque uniforme dans la cavité fluide, comme on peut le voir sur la figure 4.7. L'écoulement résultant n'est qu'un écoulement résiduel de convection, issu de l'écoulement dans le régénérateur qui traverse les grilles métalliques.

<figure><img src="image_44.png" /><figcaption>FIGURE 4.8 – Vecteurs vitesse et isovaleurs de la composante axiale de vitesse à droite de AHX, dans le plan vertical $y = -0$ (gauche expériences, droite : simulations numériques).</figcaption></figure>

Les résultats des simulations numériques et des mesures expérimentales par PIV et LDV sont en très bon accord pour les deux configurations de point de vue qualitatif et quantitatif. Les échelles de vitesse obtenues dans les différentes régions des cavités thermoacoustiques sont en accord avec les échelles classiques des écoulements de convection en cavités fluide et poreuse. Il sera donc possible de prédire l'ordre de grandeur et la topologie des écoulements de convection à l'intérieur des machines thermoacoustiques avec une approche similaire à celle présentée ici. Le résultat principal est que l'écoulement de convection et le champ de température sont tridimensionnels. Ils dépendent des propriétés thermophysiques du gaz et du milieux poreux dont la cellule active est constituée ou les échangeurs ainsi que de la géométrie et des conditions thermiques extérieures.

La question à laquelle on devra répondre par la suite est : quel sera le comportement de ces écoulements en présence de l'acoustique, après le démarrage du moteur ? Jusqu'à maintenant, dans la littérature on a considéré que la convection naturelle serait effacée par l'acoustique, mais est-ce réellement le cas ? Certaines expériences dédiées aux mesures de streaming montrent une présence possible de la convection naturelle.
# Projet de recherche

Les échanges scientifiques et humains apportés par un travail en équipe ont été le moteur dans mon activité de recherche. Il y une dizaine d'années, j'ai démarré une collaboration étroite avec C. Weisman sur la modélisation asymptotique et numérique d'un moteur thermoacoustique et la convection en milieu poreux, puis avec V. Daru sur l'étude numérique du streaming et plus récemment avec H. Bailliet qui apporte ses compétences expérimentales et sa spécificité acoustique. Nous avons ainsi constitué une équipe solide et solidaire, qui me motive fortement pour continuer dans les mêmes directions. Le contexte du projet TACOT a donné une impulsion nouvelle aux thèmes de recherches que j'ai abordés jusqu'à maintenant : le streaming acoustique et la simulation des effets multidimensionnels à l'intérieur d'une machine thermoacoustique. Dans les résonateurs des machines thermoacoustiques il y a des éléments maintenus sous un fort gradient de température comme les tubes à gaz pulsé, les stacks ou les régénérateurs. Le streaming créé par l'onde sera donc influencé par ce champ de température fortement inhomogène. La présence des plaques solides (le cas du stack) ou d'un régénérateur, les deux étant deux types de milieu poreux très différents, va modifier également le streaming. Les simulations récentes, validées par les expériences ont montré l'influence de la gravité qui engendre des écoulements caractéristiques de convection naturelle. Une question à laquelle on devra répondre sera : quel rôle joue la gravité sur le streaming ? Y a-t-il aussi un couplage entre convection naturelle et streaming ? Pendant toute cette période de recherche en thermoacoustique, je me suis un peu éloignée des techniques apprises pendant ma thèse de doctorat. Je souhaite y revenir pour essayer d'appliquer la méthode d'homogéisation des milieux périodiques pour obtenir des modèles de type milieux poreux pour les cellules actives, modèles qui pourraient être utilisés par la suite dans les simulations.

Je souhaite donc continuer le travail de recherche sur le streaming acoustique en guide d'onde à fort niveau acoustique, qui est le thème que j'ai abordé à mon arrivée au LIMSI et qui s'avère être très riche.

## Streaming acoustique

Dans l'analyse des effets thermiques sur le streaming à faible niveau nous avons trouvé un critère pour le changement de configuration des cellules de streaming extérieures, la mutation du streaming, qui permet de passer de une à deux cellules de streaming extérieur sur une étendue axiale de quart de longueur d'onde acoustique.

L'explication de cette dynamique observée numériquement et expérimentalement était la suivante : dans le contexte d'un guide d'onde d'une machine thermoacoustique, un gradient de température axial est créé par l'onde, dû au pompage thermoacoustique de chaleur. Ce gradient dépend de l'amplitude de l'onde, de la paroi solide et des échanges thermiques avec l'extérieur du résonateur. Le streaming extérieur qui se met en place par entraînement du streaming intérieur, est responsable ensuite du transport convectif de chaleur dans le résonateur quand $Re_{NL} > 1$. Le champ de température moyen devient bidimensionnel et un gradient vertical est créé (voir figure 4.9). L'acoustique sera donc modifiée et avec elle le streaming extérieur (
<figure><img src="image_45.png" /></figure>

FIGURE 4.9 – Schèma des isothermes dans un guide d'onde à faible niveau acoustique (a) et à fort niveau acoustique (b). Les rectangles blancs indiquent les régions où la température est stratifiée dans la direction transverse.

[27]). Cette explication phénoménologique de la dynamique du streaming nécessite une validation numérique car les couplages streaming-acoustique sont non-linéaires. De plus, nous avons fait deux hypothèses importantes :

- premièrement, nous avons supposé que les effets peuvent être découplés : l'acoustique crée le streaming en première phase, puis le streaming modifie l'acoustique qui modifie le streaming et on ferme la boucle. Il y a donc deux questions auxquelles on doit répondre à ce stade : notre hypothèse de séparation des effets est-elle correcte ? Jusqu'où peut-on continuer cette boucle et si elle s'arrête quel sera le mécanisme de saturation ?

-deuxièmement, nous avons fait l'hypothèse que la différence de température moyenne transversale est du même ordre de grandeur que la différence de température moyenne longitudinale, en nous basant sur une analyse dimensionnelle sur l'équation de transport de chaleur moyenne. Ceci n'a pas pu être vérifié avec les premières simulations puisque le streaming à fort niveau est obtenu pour une onde de forte intensité qui présente des chocs. Dans ce cas la température moyenne dans le résonateur augmente de manière excessive dépassant largement le petit pompage de chaleur thermoacoustique en bord de couche limite.

Il faudra donc vérifier toutes ces hypothèses avec une simulation complète du streaming en régime rapide, crée par une onde monofréquentielle, de forte intensité. V. Daru travaille actuellement sur la parallélisation du code DNS dans le but d'effectuer des simulations sur une configuration plus proche de l'expérience [67], dans laquelle l'onde est monofréquentielle et l'influence de la paroi du résonateur est prise en compte. Bien évidemment, de nouvelles expériences où on a accès au champ de température bidimensionnel à l'intérieur du résonateur seraient plus que bienvenues. Mais il est difficile de faire ce genre de mesures sans affecter l'écoulement. De plus, les précédentes meures de température en paroi ont montré de très petites valeurs de $\Delta T$ tout au long de la paroi, entre 0.3 K et 3 K, ce qui est difficilement mesurable. L'équipe d'Hélène Bailliet a fait récemment une cartographie du champ de température à différents endroits dans le résonateur avec une faisceau de thermocouples. Nous espérons que ces mesures seront utiles pour une confrontation numérique / expérimental.

En étroite liaison avec l'étude précédente, le streaming de Rayleigh pourra être analysé par la suite en résonateurs maintenus sous un important gradient de température axial, ce qui est le cas des pulse tubes dans les machines thermoacoustiques. Les pulse tubes sont des portions de résonateur, courtes par rapport à la longueur d'onde, maintenues sous une différence de température très importante (dizaines ou centaines de degrés), grâce à des échangeurs de chaleur. Ils ont le rôle d'isoler thermiquement la cellule active du reste du résonateur et jouent également un rôle déphaseur entre vitesse et pression acoustique. Olson et Swif [55] ont utilisé un pulse tube conique dans le but de réduire, voire annuler le streaming de Rayleigh. Ils ont développé un modèle pour estimer le streaming massique dans la couche limite en fonction de l'angle entre la paroi conique et l'axe du résonateur, et en fonction du gradient de température axial. Le modèle leur a permis de déduire une valeur de l'angle qui annule le streaming dans le pulse tube, prédictions confirmées indirectement par des mesures de température. Il serait alors intéressant de
pouvoir simuler des écoulements de streaming dans des résonateurs à section variable. L'idée serait, une fois la dynamique du streaming comprise, de trouver des solutions de contrôle du streaming : suppression ou augmentation du streaming, si celle-ci présente un intérêt, en jouant sur son couplage avec les effets de température ou sur la géométrie.

## Streaming acoustique : influence de la gravité

Les observations expérimentales des premiers travaux sur le streaming durant la thèse de P. Debesse avaient montré des déformations du profil de streaming, non expliquées par les modèles linéaires ou les simulations de la littérature. A l'époque nous avions pensé à une superposition de l'écoulement de convection naturelle et le streaming de Rayleigh, mais nous n'avions pas encore les moyens de calcul pour estimer l'écoulement de convection. La possibilité que le streaming soit modifié par la convection naturelle et les résultats récents concernant les effets thermiques sur le streaming ont relancé les études sur la simulation des écoulements de convection en cavité fluide/poreux, avec des besoins de modélisation.

Récemment, les études de convection naturelle démarrées dans le cadre du projet TACOT ont montré l'existence d'un écoulement de convection dans la cavité fluide située entre l'extrémité du résonateur et l'échangeur froid en sortie du moteur thermoacoustique, même pour de faibles différences de température axiale. Les simulations sont en très bon accord avec les mesures réalisées sur un système avec la même géométrie [40]. L'écoulement observé est du même ordre de grandeur que le streaming de Rayleigh prédit par un modèle linéaire de streaming lent, associé à une onde dont l'amplitude est celle de l'onde générée par le moteur thermoacoustique, pour la même différence de température. L'écoulement est tridimensionnel, comme dans la thèse de P. Debesse. Les études envisagées par la suite ont comme objectif l'analyse du couplage entre la convection naturelle et le streaming acoustique. L'équipe d'Hélène Bailliet a effectué des mesures d'écoulements moyens en présence et en absence de l'onde pour le moteur à ondes stationnaires présenté dans le chapitre 4. Notre contribution sera du point de vue de la simulation : il s'agit d'étudier et de simuler le couplage des deux écoulements.

Il y a plusieurs manières d'aborder le problème : on peut considérer comme état de base l'écoulement de convection naturelle et perturber cet état par une oscillation acoustique modélisée soit comme une source volumique, soit comme un piston mobile. On pourra étudier si l'acoustique qui en résulte est modifiée par le champ de température tridimensionnel issu de la convection naturelle, ce que semble indiquer le modèle présenté dans le chapitre précédent [27]. Des mesures précises pour la composante verticale acoustique pourraient compléter ce travail. Certains auteurs avaient conclu récemment que le mécanisme de convection naturelle est inhibé par l'acoustique et remplacé par un écoulement de convection forcée [57, 56]. La question de la convection naturelle reste donc ouverte.

Une autre stratégie va consister à considérer les équations de Navier-Stokes compressibles, modifiées par la présence de la gravité, et à effectuer des simulations directes sur les équations complètes pour calculer à la fois l'acoustique et le streaming. Ces calculs pourront bien sûr se faire initialement sur des configurations géométriques simplifiées et de taille réduite, mais qui pourront être choisies en similitude avec l'expérience, à la fois du point de vue de l'écoulement et des conditions thermiques. Cette étude sera plus une étude fondamentale, sur un guide d'onde représenté par un canal plan bidimensionnel, contenant un fluide et soumis à une différence de température entre les deux frontières verticales.

## Application à la réfrigération thermoacoustique : projet TACOT

Enfin, depuis 2017 nous faisons partie d'un projet ANR TACOT dont le porteur est H. Bailliet du laboratoire Pprime et qui regroupe plusieurs équipes de recherches françaises (LAUM,
LIMSI, ECL) et un partenaire industriel PSA. L'objectif est la construction et l'étude d'une pompe à chaleur thermoacoustique compacte pour le transport terrestre. Un prototype instrumenté et compact a été construit à l'Institut Pprime suite à la collaboration étroite entre G. Poignand (LAUM) et I. Ramadan et H. Bailliet (Pprime), avec les conseils précieux de D. Gardner (ancien collaborateur de G.Swift). Le prototype est constitué de deux sources acoustiques qui permettent de régler le déphasage entre pression et vitesse et réaliser ainsi une géométrie compacte sans avoir besoin d'un résonateur [60]. La pompe à chaleur est une cellule active de type thermoacoustique, constituée d'un régénérateur placé entre deux échangeurs de chaleur (un échangeur froid et un échangeur à température ambiante) (Figure 4.2.2). Entre le cœur thermoacoustique et les parois de l'enceinte il y a un petit espacement pour permettre le passage du gaz et réaliser une boucle de retour. Le gaz utilisé est un mélange d'Helium et d'Argon à pression moyenne de 40 bar (valeur maximale testée). Pour 210 W de puissance électrique fournie, 200 W sont pompés à la source froide avec une différence de température de 36 K entre les échangeurs (pour un déphasage de 180° entre les deux sources), c'est à dire un gradient de température de 900 K/m. Le COP global de la machine est 0.95, et celui de la conversion d'énergie thermoacoustique seule est de 2.37. Ces résultats montrent l'intérêt du projet et des applications possibles pour le secteur automobile. Le prototype a été instrumenté avec des thermocouples et des capteurs de pression pour pouvoir suivre l'évolution de la température dans la cellule active et contrôler les pressions moyenne et acoustique. Cependant l'instrumentation du prototype est délicate et les grandeurs mesurées disponibles sont insuffisantes pour suivre le mouvement du gaz et les échanges de chaleur à l'intérieur de la machine. La simulation numérique doit donc pouvoir apporter des informations précieuses.

<figure><img src="image_46.png" /><figcaption>FIGURE 4.10 – Projet ANR TACOT, schéma du prototype construit à l'institut Pprime.</figcaption></figure>

La contribution de notre équipe (C. Weisman, V. Daru, Y. Fraigneau) au projet concerne la modélisation et la simulation des phénomènes non linéaires présents dans ce type de machines : convection naturelle, effets multidimensionnels et streaming.

La deuxième configuration expérimentale sur laquelle nous avons fait les calculs de convection avec O. Hireche contient une cellule thermoacoustique similaire à celle du prototype avec quelques différences : le gradient thermique plus petit dans le prototype et sa géométrie compacte. Une estimation du nombre de Rayleigh basée pour la configuration du prototype TACOT nous indique que le nombre de Rayleigh poreux est du même ordre de grandeur que dans le cas d'étude précédent [40]. Il sera nécessaire donc de faire des calculs de convection naturelle dans cette configuration, d'autant plus que les premières mesures de température dans le prototype montrent une différence de température importante (de l'ordre de 30K) sur une ligne verticale, entre le haut et le bas du régénérateur.

L'étude des effets multidimensionnels et de streaming en machine compacte sera réalisée avec deux outils numériques différents : un code compressible utilisable à des hautes fréquences
et un code Faible Mach adapté pour des basses fréquences (donc des géométries courtes par rapport à la longueur d'onde). L'approximation Faible Mach peut être particulièrement intéressante pour limiter les temps de calcul, et son utilisation se justifie par la compacité des prototypes du projet TACOT par rapport à la longueur d'onde. Nous souhaitons appliquer cette approximation à la modélisation du pompage de chaleur thermoacoustique dans le prototype, en géométrie axisymétrique tout d'abord. Le code SunfluidH développé au LIMSI par Y. Fraigneau est un code 3D instationnaire qui résout les equations de Navier-Stokes en approximation faible Mach (modèle présenté dans le chapitre 3).

Le premier travail envisagé est de comparer les résultats des deux codes et estimer la gamme de fréquences dans laquelle leurs résultats sont similaires pour déterminer ensuite une stratégie de leur utilisation (collab. Y. Fraigneau, V. Daru, C. Weisman). Une première étude de comparaison a démarré l'an dernier durant le stage de A. Krishakumar, sur un cas académique d'une cavité remplie de fluide, dans laquelle l'onde est créé par un piston. Les premiers résultats ont montré la validité de l'approximation Faible Mach concernant les caractéristiques de l'écoulement oscillant obtenu. Le travail sera poursuivi d'une étude comparative quant au pompage de chaleur thermoacoustique à proximité d'une plaque solide à l'intérieur de la cavité et d'un empilement de plaques, ainsi que l'écoulement moyen de streaming.

Elle sera suivie de simulations en présence d'un milieu poreux qui constitue la cellule active d'une machine thermoacoustique. Les cellules actives les plus efficaces sont celles constituées d'un milieu poreux avec une géométrie de la microstructure bien plus compliquée que celle d'un stack, avec une taille caractéristique des pores plus petite que l'épaisseur des couches limites thermique et visqueuse acoustiques. Les modèles de milieux poreux utilisés dans la littérature thermoacoustique sont peu nombreux. Le modèle généralement utilisé dans la construction des machines thermoacoustiques est celui des milieux poreux tortueux soumis aux oscillations harmoniques, basé sur des coefficients de friction et de transfert de chaleur spécifiques des écoulements quasi-stationnaires [78]. Plus récemment, la théorie linéaire de la thermoacoustique a été étendue pour un milieu poreux où l'écoulement est régi par la loi de Darcy, ou l'extension semi-empirique de Brinkman-Forchheimer, en effectuant une approximation acoustique de ces équations moyennées en espace [79]. Dans ces modèles on suppose que le déplacement acoustique est assez grand pour qu'on puisse négliger l'influence du pore sur la propagation acoustique et le transfert de chaleur.

Le but du travail futur sera d'établir un modèle de propagation thermo-acoustique en milieu poreux, en partant de l'approximation acoustique au niveau du pore pour arriver à des équations moyennes (en espace) caractérisant le milieu poreux à l'échelle macroscopique, avec une technique basée sur des développements asymptotiques en échelles multiples. Ceci nous permettra de développer des équations de propagation acoustique et de transfert de chaleur en milieux poreux pour n'importe quelle géométrie du pore et pour différents ordres de grandeur de l'écoulement à travers le milieu poreux. Ces équations vont nourrir ensuite les simulations numériques dans le code compressible, où on peut envisager de séparer les opérateurs pour l'acoustique et pour le transport convectif (séparation d'échelles de temps), permettant ainsi d'analyser l'influence des caractéristiques du milieu poreux sur les performances des machines.
Bibliographie

[1] ADEFF J.A., HOFLER T.J, 2005 Design and construction of a solar-powered, thermoacoustically driven, thermoacoustic refrigerator *J. Acoust. Soc. Am.* **107**(6), L37-L42.

[2] AKTAS M., FAROUK B. 2004 Numerical simulation of acoustic streaming produced by finite-amplitude resonant oscillations in an enclosure. *J. Acoust. Soc. Am.* **116**, 2822-2831.

[3] AKTAS, M. K. & OZGUMUS, T. 2010 The effects of acoustic streaming on thermal convection in an enclosure with differentially heated horizontal walls. *Int. J. Heat and Mass Transfer* **53** (23-24), 5289-5297.

[4] ARNOTT W. P., BASS H. E., RASPET R. 1991 General formulation of thermoacoustics for stacks having arbitrarily shaped pore cross section, J. Acoust. Soc. Am. **90**, 3228-3237.

[5] ARROYO M.P., GREATED C.A. 1991 Stereoscopic particle image velocimetry. *Meas. Sci. Technol.* **2**,1181-1186.

[6] ATCHLEY A. A., BASS H. E., HOFLER T. J., LIN H. T. 1992 Study of a thermoacoustic prime mover below onset of self-oscillation, *J. Acoust. Soc. Am.* **91**, 734-743.

[7] ATCHLEY A. A. AND F. M. KUO 1994 Stability curves for a thermoacoustic prime mover. *J. Acoust. Soc. Am.* **95**, 1401-1404.

[8] S. BACKHAUS, E. TWARD, M. PETACH 2004 Traveling wave thermoacoustic electric generator. *Appl. Phys. Lett.* **85** 1085.

[9] BAILLIET H., GUSEV V., RASPET R. & HILLER R. A. 2001 Acoustic streaming in closed thermoacoustic devices. *J. Acoust. Soc. Am.* **110**, 1808-1821.

[10] BALTEAN D., LÉVY T., BALINT S. 2003 Diffusion-convection in a porous medium with impervious inclusions at low flow rates. *Transport in Porous Media* **51**, 19-39.

[11] BALTEAN-CARLÈS D., DARU V., WEISMAN C., TABAKOVA S. & BAILLIET H. 2019 An unexpected balance between outer Rayleigh streaming sources. *J. Fluid Mech.* **867**, 985-1011.

[12] A. BEJAN 1985 Convection Heat Transfer, *Mac Graw Hill*.

[13] NIELD D.A., BEJAN A. 2006 Convection in Porous Media, Springer and Business Media, New York

[14] BESNOIN E. 2001 Numerical study of thermoacoustic heat exchangers. Thèse de doctorat, Johns Hopkins University.

[15] BOLURIAAN S. & MORRIS P. J. 2003 Numerical simulation of acoustic streaming in high amplitude standing wave. *9th AIAA/CEAS Aeroacoustics Conference and Exhibit 12-14 May 2003, Hilton Head, South Carolina AIAA 2003-3152*

[16] BROOMHEAD D.S., KING G.P. 1986 Extracting qualitative dynamics from experimental data. *Physica D: Nonlinear Phenomena* **20**, 217-236.

[17] CAMPBELL M., COSGROVE J.A., GREATED C.A., JACK S., ROCKLIFF D. 2000 Review of LDA and PIV applied to the measurement of sound and acoustic streaming. *Opt. Laser Technol* **32**, 629-639.
[18] CERVENKA M. & BEDNARRIK M. 2017 Effect of inhomogeneous temperature fields on acoustic streaming structures in resonators. *J. Acoust. Soc. Am.* **141**, 4418–4426.

[19] CERVENKA M. & BEDNARRIK M. 2018 Numerical study of the influence of the convective heat transport on acoustic streaming in a standing wave. *J. Acoust. Soc. Am.* **143**, 727.

[20] CHINI G. P., MALECHA Z. & DREEBEN T. B. 2014 Large-amplitude acoustic streaming. *J. Fluid Mech.* **744**, 329–351.

[21] CHURCHILL S. W. & CHU H. H. S. 1975 Correlating equations for laminar and turbulent free convection from a horizontal cylinder. *Int. J. Heat Mass Transfer* **18**(9), 1323–1329.

[22] DAI W., LUO E., HU J., LING H. 2005 A heat driven thermoacoustic cooler capable of reaching liquid nitrogen temperature. *Appl. Phys. Lett.* **86**, 224103.

[23] V. DARU, C. TENAUD 2009 Numerical simulation of the viscous shock tube problem by using a high resolution monotonicity preserving scheme, *Computers and Fluids* **38**(3), 664-676.

[24] DARU V., BALTEAN-CARLÈS D., WEISMAN C., DEBESSE P. & GANDIKOTA G. 2013 Two-dimensional numerical simulations of nonlinear acoustic streaming in standing waves. *Wave Motion* **50**, 955-96

[25] DARU V., REYT I., BAILLIET H., WEISMAN C. & BALTEAN-CARLÈS D. 2017a Acoustic and streaming velocity components in a resonant wave guide at high acoustic levels. *J. Acoust. Soc. Am.* **141** (1), 563–574

[26] DARU V., WEISMAN C., BALTEAN-CARLÈS D., REYT I., BAILLIET H. 2017b Inertial effects on acoustic Rayleigh streaming flow : Transient and established regimes. *Wave Motion* **74**, 1–17.

[27] DARU V., WEISMAN C., BALTEAN-CARLÈS D., & BAILLIET H. 2020 Thermal effects on fast acoustic streaming inside a resonator. soumis au *J. Fluid Mech.*

[28] DEBESSE P., BALTEAN-CARLES D., LUSSEYRAN F., FRANÇOIS M.-X., 2014 Oscillating and streaming flow identification in a thermoacoustic resonator from under-sampled PIV measurements. *Meas. Sci. Technol.* **25**, 025005.

[29] DUFFOURD S., 2001 Réfrigérateur thermoacoustique : Études analytiques et expérimentales en vue d'une miniaturisation. *Dissertation, Ecole Centrale de Lyon*.

[30] DUTHIL P., BALTEAN CARLÈS D., BÉTRANCOURT A., FRANÇOIS M.-X., YU Z.B., THERMEAU J.P, 2006 Experiments and valve modelling in Thermoacoustic device. in *Advances in Cryogenic Engineering* **51 B**, 1095-1102.

[31] DUTHIL P., WEISMAN C., BRETAGNE E., FRANÇOIS M-X. 2004 Experimental and numerical investigation of heat transfer and flow within a thermoacoustic cell, *International Journal of Transport Phenomena* **6 (4)**, 265-272.

[32] FAURE T.M., LUSSEYRAN F., GOUGAT P., LAUNAY F. 2006 Experimental and numerical investigation of the flow distribution inside a tubular heat exchanger. *J. Fluids Eng.* **128**,1218–1227.

[33] GARDNER D. L. AND SWIFT G. W. 2003 A cascade thermoacoustic engine *J. Acoust. Soc. Am.* **114**(4), 1905-1919.

[34] D. GOBIN AND B. GOYEAU 2008 Natural convection in partially porous media : a brief overview, *Int. J. Num. Meth. Heat Fluid Flow, Emerald* **18 (3-4)**, 465–490.

[35] GOPINATH A. & MILLS A. F. 1994 Convective heat transfer due to acoustic streaming across the ends of a Kundt tube, *J. Heat Transfer* **102**(3) 47-53
[36] GOPINATH A., TAIT N. L. & GARRETT S. L. 1998 Thermoacoustic streaming in a resonant channel : The time averaged emperature distribution. *J. Acoust. Soc. Am.* **102**(3), 1388-1405.

[37] HAMILTON M. F., ILINSKII Y. A. & ZABOLOTSKAYA E. A. 2003a Acoustic strea-
ming generated by standing waves in two-dimensional channels of arbitrary width.
*J. Acoust. Soc. Am.* **113**(1), 153–160.

[38] HAMILTON M. F., ILINSKII Y. A. & ZABOLOTSKAYA E. A. 2003 Thermal effects on acoustic streaming in standing waves. *J. Acoust. Soc. Am.* **114**(6), 3092–3101.

[39] HIRECHE O., WEISMAN C., BALTEAN-CARLÈS D., LE QUÉRÉ P., BAUWENS L. 2010 Low Mach number analysis of idealized thermoacoustic engines with numerical solution, *J. Acoust. Soc. Am.* **128** (6), 3438–3448.

[40] HIRECHE O., RAMADAN I., WEISMAN C., BAILLIET H., FRAIGNEAU Y., BALTEAN-CARLES D., DARU V. 2020 Experimental and numerical investigation of natural convection flows in two horizontal thermoacoustic cavities. *Int. J. Heat Mass Transfer* **149**, 119195.

[41] INCROPERA F. P., DEWIT D. P., BERGMAN T. L. & LAVINE A. S. 2013 Principles of heat and mass transfer. Seventh Edition. *John Wiley & Sons*.

[42] JENSEN C., AND RASPET R. 2010 Thermoacoustic properties of fibrous materials. *The Journal of the Acoustical Society of America* **127**, 3470–3484.

[43] LIGHTHILL, J. 1978 Acoustic streaming. *J. Sound Vib.* **61**, 391–418.

[44] LIN Y. & FAROUK B. 2008 Heat transfer in a rectangular chamber with differentially heated horizontal walls : Effects of a vibrating sidewall. *Int. J. Heat Mass Transfer* **51**, 3179–3189.

[45] MA L., WEISMAN C., BALTEAN-CARLÈS D., DELBENDE I., BAUWENS L. 2015 Effect of a resistive load on the starting performance of a standing wave thermoacoustic engine : A numerical study. *J. Acoust. Soc. Am.* **138 (2)**, 847–857.

[46] WOLFRAM RESEARCH, INC. 2018 Mathematica, Version 11.3, Champaign, IL.

[47] MENGUY L. & GILBERT J. 2000 Non-linear Acoustic Streaming Accompanying a Plane Stationary Wave in a Guide. *Acta Acustica united with Acustica* **86**, 249–259.

[48] MERKLI P. & THOMANN H. 1975 Thermoacoustic effects in a resonance tube. *J. Fluid Mech.* **70**, 161–177.

[49] MICHEL G. & CHINI G. P. 2019 Strong wave-mean-flow coupling in baroclinic acoustic streaming. *J. Fluid Mech.* **858**, 536–564.

[50] MOREAU S., BAILLIET H. & VALIÈRE J.-C. 2008 Measurements of inner and outer streaming vortices in a standing waveguide using laser doppler velocimetry. *J. Acoust. Soc. Am.* **123**(2), 640–647.

[51] MORZUKEWICH, G. 2002 Heat transport by acoustic streaming within a cylindrical resonator. *Applied Acoustics* **63**, 713–735.

[52] NABAVI M., SIDDIQUI K. & DARGAHI, J. 2008 Influence of differentially heated hori-
zontal walls on the streaming shape and velocity in a standing wave resonator. *International
Communications in Heat and Mass Transfer* **35**, 106–1064.

[53] NYBORG, W. L. 1953 Acoustic streaming due to attenuated plane waves. *J. Acoust. Soc. Am.* **25**, 68–75.

[54] NYBORG, W. L. 1965 Acoustic streaming. in *Physical Acoustics, W. P. Mason (ed)*, Academic Press, New York **2B**, 265–331.
[55] OLSON J. R. & SWIFT G. W. 1997 Acoustic streaming in pulse tube refrigerators : Tape-red pulse tubes. *Cryogenics*. **37**, 769-776.

[56] PAN N., WANG S., SHEN C. 2012 Visualization investigation of the flow and heat transfer in thermoacoustic engine driven by loudspeaker, *Int. J. Heat and Mass Transf.* **55**, 7737-7746.

[57] PAN N., WANG S., SHEN C. 2014 A fundamental study on characteristic of thermoacoustic engine with different tilt angles, *Int. J. Heat and Mass Transf.* **74**, 228-237.

[58] PAOLUCCI S. 1982 On filtering of sound from the Navier-Stokes equations. *Sandia National lab. Report SAND 82-8257*, (non publié).

[59] PASTUR L.R., LUSSEYRAN F., FAURE T.M., FRAIGNEAU Y., PETHIEU R., DEBESSE P. 2008 Quantifying the non-linear mode competition in the flow over an open cavity at medium Reynolds number. *Exp. Fluids* **44**, 597-608

[60] POIGNAND G., LIHOREAU B., LOTTONP., GAVIOT E., BRUNEAU M. AND GUSEV V. 2007 Optimal acoustic fields in compact thermoacoustic refrigerators, *Applied Acoustics* **68** 642-659.

[61] QUÉNOT G, PAKLEZA J, KOWALEWSKI TA 1998 Particle image velocimetry with optical flow. *Exp. Fluids* **25**,177-189

[62] RAMADAN I., BAILLIET H. AND VALIÈRE J.-C. 2018 Experimental investigation of the influence of natural convection and end-effects on Rayleigh streaming in a thermoacoustic engine, *J. Acoust. Soc. Am.* **143**(1), 361–372.

[63] RASPET R., BRUSTER J., BASS H. E 1998 New Approximation Method for Thermoacoustic Calculations, *J. Acoust. Soc. Am.* **103**(5), 2395-2402

[64] LORD RAYLEIGH 1884 On the circulation of air observed in Kundt tubes, and on some allied acoustical problems. *Philos. Trans. R. Soc. London* **175**, 1-21

[65] RÉMUS J. http://jacques-remus.fr/

[66] REYT I., DARU V., BAILLIET H., MOREAU S., VALIÈRE J.-C., BALTEAN-CARLÈS D. & WEISMAN C. 2013 Fast acoustic streaming in standing waves : Generation of an additional outer streaming cell. *J. Acoust. Soc. Am.* **134**, 1791-1801.

[67] REYT I., BAILLIET H. & VALIÈRE J.-C. 2014 Experimental investigation of acoustic streaming in a cylindrical wave guide up to high streaming Reynolds numbers . *J. Acoust. Soc. Am.* **135** (1), 27-37.

[68] ROH H., RASPET R., AND BASS H. E. 2007 Parallel capillary-tube-based extension of thermoacoustic theory for random porous media. *The Journal of the Acoustical Society of America* **121**(3), 1443-1421.

[69] ROTT, N. 1974 The influence of heat conduction on acoustic streaming. *Zeitschrift für Angewandte Mathematik und Physik (ZAMP)* **25**, 417–421.

[70] N. ROTT 1980 Thermoacoustics. *Adv. Appl. Mech.* **20**, 135-243.

[71] SANCHEZ-PALENCIA, E. 1980 Non-Homogeneous Media and Vibration Theory. *Lecture Notes in Physics* **127**, Springer Verlag Berlin

[72] SCHLICHTING H 1932 Berechnung ebener periodischer Grenzschichtströmungen (calculation of plane periodic boundary layer streaming). *Phys. Zeit* **33**, 327-335.

[73] SCHUSTER V. K. & MATZ W. 1940 Über stationare Strömungen im Kundtschen Rohr (On Stationary Streaming in Kundt Tubes). *Akust. Zeitschrift* **5**, 349-352.

[74] SPOELSTRA S., TIJANI. M.E.H. 2005 Thermoacoustic heat pumps for energy savings, in seminar "Boundary crossing acoustics" of the Acoustical Society of the Netherlands
[75] STANISLAS M., OKAMOTO K., KÄHLER C.J., WESTERWEELN J. 2005 Main results of the Second International PIV Challenge. *Exp. Fluids* **39**,170–191.

[76] STANISLAS M., OKAMOTO K., KÄHLER C.J., WESTERWEEL J., SCARANO F. 2008 Main results of the third international PIV Challenge. *Exp. Fluids* **45**, 27–71.

[77] STUART J.T. 1966 Double boundary layer in oscillatory viscous flow. *J. Fluid Mech.* **24**, 673–987

[78] SWIFT, G. W. 1988 Thermoacoustic engines. *J. Acoust. Soc. Am.* **84**, 1145-1180.

[79] S. H. TASNIM, 2011 Porous Media Thermoacoustic Stacks : Measurements and Models. *Thesis, University of Waterloo, Ontario, Canada*

[80] THOMPSON M. W., ATCHLEY A. A. & MACCARONE M. J. 2005 Influences of a temperature gradient and fluid inertia on acoustic streaming in a standing wave. *J. Acoust. Soc. Am.* **117**(4) 1839–1849.

[81] WARD B., CLARK J. AND SWIFT G. W. 2008 Design Environment for Low-Amplitude Thermoacoustic Energy Conversion (DeltaEC) User's Guide (version 6.2). Los Alamos National Laboratory.

[82] WEISMAN C., BALTEAN-CARLÈS D., DUTHIL P., LE QUÉRÉ P. 2008 Natural convection in a stack of horizontal plates in a differentially heated cavity, in : *Proceedings of the 19th International Symposium on Transport Phenomena (ISTP-19), Reykjavik, Island 2008*.

[83] WESTERVELT P. J. 1953 The Theory of Steady Rotational Flow Generated by a Sound Field. *J. Acoust. Soc. Am.* **25**, 60–67.

[84] WORLIKAR A. S., AND KNIO O. M., 1996 Numerical simulation of a thermoacoustic refrigerator. I : Unsteady adiabatic flow around a thermoacoustic refrigerator. *J. Comput. Phys.* **127**, 424-451

[85] WORLIKAR A. S., KNIO O. M., AND KLEIN R., 1998 Numerical simulation of a thermoacoustic refrigerator. II : stratified flow around the stack. *J. of Comput. Fluid.* **144**, 299–324.

[86] ZOONTJENS L., HOWARD C., ZANDER A. AND B. CAZZOLATO 2005. Feasibility Study of an Automotive Thermo acoustics Refrigerator in : *Proceedings of Acoustics, Busselton, Australia*
P. Debesse, D. Baltean Carlès, F. Lusseyran and M-X François.

Oscillating and streaming flow identification in a thermoacoustic resonator, from
undersampled PIV measurements.

Meas. Sci. Technol. 25, 025005 (16pp), (2014).
# Oscillating and streaming flow identification in a thermoacoustic resonator, from undersampled PIV measurements

P Debesse<fnref n="1" />,<fnref n="2" />, D Baltean Carlès<fnref n="1" />,<fnref n="2" />, F Lusseyran<fnref n="2" /> and M-X François<fnref n="1" />,<fnref n="2" />

E-mail: baltean@limsi.fr

Received 15 July 2013, revised 7 October 2013

Accepted for publication 12 November 2013

Published 3 January 2014

## Abstract

The nonlinear phenomena that appear in thermoacoustic systems generate continuous secondary flows (streaming), superimposed on the dominant acoustic oscillations, reducing system efficiency. In this paper, an experimental study is presented for evaluating both acoustics and streaming in such a system. The velocity measurements, based on particle image velocimetry, are undersampled with respect to the acoustic period and without external phase reference. A new method is proposed for treating the undersampled data, in which the acoustic field is reconstructed via an embedding technique. The application of this method and use of different data processing tools provide results in good agreement with linear acoustic theory. The harmonic content of the velocity as well as the wavelength are obtained. The estimate of the mean velocity over the time period indicates the existence of a continuous flow different from the Rayleigh-Schlichting streaming.

Keywords: acoustic streaming, thermoacoustics, particle image velocimetry, phase averaging, singular value decomposition

(Some figures may appear in colour only in the online journal)

## 1. Introduction

Thermoacoustic systems (Swift 2002) convert thermal energy into acoustic energy (thermoacoustic engines) and vice versa (thermoacoustic refrigerators or heat pumps). Such devices consist of a long resonator tube filled with gas, containing an active cell, which is composed of a porous solid matrix placed between a heater and a cooler. In the past twenty years, great effort has been made by researchers to develop such systems, because of the large field of industrial applications, ranging from cryogenics to air conditioning and domestic refrigeration. Despite these efforts, a lot of phenomena limiting the efficiency of thermoacoustic systems are still not well understood. Acoustic streaming and in general continuous flows, source of convective heat transport, are nonlinear phenomena responsible for thermal losses. Acoustic streaming is a continuous secondary flow, superimposed on the first-order acoustic motion, generated in a resonator as an interaction

between an acoustic wave and a solid boundary, caused by Reynolds stresses. In the case of a wide standing wave resonator (half-wavelength), it consists of four steady vortices outside the viscous boundary layer, referred to as Rayleigh streaming (outer streaming). Inside the boundary layer, there are other four vortices, referred to as Schlichting streaming (inner streaming), rotating in directions opposite to those of outer vortices. The streaming patterns are represented in figure 1 and are characteristic of slow (linear) streaming.

These types of continuous flows have been modelled and precisely measured for resonators with no active thermoacoustic cell (Rayleigh 1884, Stuart 1966, Schlichting 1932, Hamilton et al 2003, Arroyo and Greated 1991, Campbell et al 2000, Bailliet et al 2001, Thompson and Atchley 2005, Nabavi et al 2007). In the case of a realistic thermoacoustic system, the acoustic streaming can be influenced by the presence of the active cell maintained under a large temperature gradient, and even more by the

<footnote n="1">¹ Université Pierre et Marie Curie, 4 Place Jussieu, F-75252 Paris Cedex 05, France</footnote>

<footnote n="2">² LIMSI-UPR CNRS 3251, BP133, F-91403 Orsay Cedex, France</footnote>
<figure><img src="image_47.png" /><figcaption>**Figure 1.** Streaming patterns (streamlines) for a standing wave in a half-wavelength resonator obtained with the model of Hamilton *et al* (2003).</figcaption></figure>

nonlinear regime attained for high intensity waves with large harmonic content. For high intensity waves, the streaming itself becomes nonlinear and the streaming cells are deformed (experimental studies: Arroyo and Greated 1991, Campbell *et al* 2000, Thompson *et al* 2005, Moreau *et al* 2008, Nabavi *et al* 2009, Reyt *et al* 2012; theoretical studies: Menguy and Gilbert 2000, Aktas and Farouk 2004, Daru *et al* 2013). In this nonlinear regime, few experiments were performed in the case when the wave is generated with a thermoacoustic engine, under high mean pressure, driven at very high acoustic pressures.

In this type of system, the measurement of both acoustic velocity field and acoustic streaming velocity is a challenge in metrology and signal processing, and requires high precision. Indeed, both phenomena have two very different time scales, some m s⁻¹ for the acoustic component and some cm s⁻¹ or even less for the continuous component. However, non-intrusive optical velocimetry, laser Doppler velocimetry (LDV) and particle image velocimetry (PIV), mainly applied for turbulent flows, can be used in the acoustic context if associated with adapted post-treatments.

Thompson and Atchley (2005) used LDV in a cylindrical standing-wave resonator to study the streaming onset and also the influence of temperature gradient and fluid inertia on the developed streaming. Moreau *et al* (2008) measured acoustic streaming at different points in a standing-wave flow as the nonlinear Reynolds number increases, and highlighted the inner vortices. Desjouy *et al* (2009) compared successfully experimental results obtained in a closed-loop waveguide with theoretical results and revealed that acoustic velocity and streaming are impacted by the resonator curvature. In works employing PIV measurements (Duffourd 2001, Michard *et al* 2004, Nabavi *et al* 2007, 2009), acoustic and streaming fields are always extracted using phase averaging based on well-sampled external signal. Indeed, the phase averaging method is necessarily linked to the identification of a periodic component of the phenomenon under study. The general framework of the conditional averaging allows structures identification even for non-periodic phenomena, but requires a specific criterion (see pioneering works on turbulence, for example Adrian 1979, Hussain and Haykawa 1987). In the case of thermoacoustics, the average process has to be linked to the dominant acoustic oscillation, and therefore a phase reference is needed.

Classical phase averaging based on an external reference (Kiya and Matsumara 1985), presents a risk of bias induced by phase jitter between the trigger and the averaged signal. It is of crucial importance to avoid this effect when the structure of interest is several orders of magnitude smaller than the leading velocity, as in our case. A more favourable case can be mentioned, when the external reference is the triggering force also. In this case, the phase jitter may not affect the averaged signal, as shown by Widmann *et al* (2012), who identified very small structures (Tollmien-Schlichting waves) having the same magnitude as the turbulent noise. To overcome a potential loss of phase coherence, a few procedures are proposed in the literature, which are not based on an external reference signal. One way consists in using local information taken inside the measured field itself (PIV time resolved velocity records), as in the case of the average of the wake behind a cylinder, either using as phase reference the circulation computed at a precise location inside the recirculation domain (Lourenco *et al* 1997) or using band-pass filtered time series (Negri *et al* 2011). Still, the results appear to be very dependent on details of the procedure.

More recently, instead of using local information taken inside the velocity field, global approaches are proposed. Konstantinidis *et al* (2005) used the cross-correlation between velocity fields as a phase identifier. This procedure allows randomly and not resolved time sampling, but cannot identify the value of each phase, the phase ordering being done visually. With another global approach, Perrin *et al* (2007) obtained the phase angle of the vortex shedding in the wake behind a cylinder using the two first modes obtained with the proper orthogonal decomposition (POD) method applied to the velocity snapshots. This phase averaging gives a more precise description than those using an external pressure signal, but the quality of the phase reference depends on the spectral content of the chosen POD modes.

The present study proposes an original method of signal processing (Debesse 2008), within the frame of PIV, based on the embedment technique of the underlying dynamics proposed by Broomhead and King (1986). The embedment technique has been successfully used in the case of self-sustained oscillations developing in open cavity flows in order to quantify the rate of intermittency of nonlinearly competing modes (Pastur *et al* 2008). In the present work,
<figure><img src="image_48.png" /><figcaption>Figure 2. Simplified diagram of the apparatus: in green, the PIV measurement zone; M, the thermoacoustic engine; $P_{1k}$, the dynamical pressure transducer; $P_m$, the mean pressure transducer; the proportions are not respected.</figcaption></figure>

the embedment technique was adapted to a different aim: the reconstruction of both acoustic part and time-independent part from experimental undersampled PIV velocity fields. No reference signal is needed to reconstruct the velocity field through phase-averaged series. Singular value decomposition (SVD) projects, on their principal components, the velocity fields strongly undersampled with respect to the fundamental acoustic period. This allows reconstruction of a pseudo phase space and then defining a self-referenced phase inside the acoustic fundamental period. The velocity fields can be then rearranged with respect to their phase, initially randomly distributed in time. Hence, we obtain a signal with enough information in space and time through phase sampling. A phase average over the fundamental period is then possible, yielding the mean velocity corresponding to second-order steady effect. The embedment technique together with the estimation of the mean velocity (second-order) from undersampled velocity PIV series represents an original method for measuring second-order steady flow in thermoacoustic waves.

This approach is developed and applied to experiments in a cylindrical standing-wave resonator, with a thermoacoustic engine acting as a wave generator. The following sections of this paper present the experimental setup (section 2), the acoustic velocity field estimate (section 3), the self-referenced phase arrangement (section 4), the streaming flow estimate and analysis (section 5) and an overall summary (section 6).

## 2. Experimental setup and metrology

The experimental device contains two major parts: the thermoacoustic system and the PIV measurement zone (figure 2).

The thermoacoustic system consists of a cylindrical steel tube closed at both ends, with length $L = 7.88$ m and inner diameter $D = 56.3$ mm, suspended by chains and shock absorbers in order to reduce structure vibration losses.

The thermoacoustic engine (active cell) consists of a stack of plates, placed between two heat exchangers. The wave is generated and maintained inside the system with a temperature gradient imposed between the hot heat exchanger and the cold heat exchanger. The hot heat exchanger is heated by an electrical power supply, while the cold heat exchanger is maintained at constant temperature with cooling water. This allows hot temperatures to be obtained in the range of 423 K (150 °C) to 873 K (600 °C); for our experiments, the temperature differences between hot and cold exchangers were 458 K (185 °C) to 573 K (300 °C). The wave generated inside

<figure><img src="image_49.png" /><figcaption>Figure 3. Power spectrum of the pressure recorded by transducer $P_1$.</figcaption></figure>

the resonator is a standing wave. The corresponding acoustic pressure amplitudes vary within the range 0.1 to 0.5 bar.

The resonator is filled with nitrogen, confined under a mean pressure of $P_m = 7.26$ bar ($7.26 \times 10^5$ Pa). The experiments were conducted at the mean temperature $T = 289$ K. The kinematic viscosity is $\nu = 2.09 \times 10^{-6}$ m$^2$ s$^{-1}$ and the mean density is $\rho = 8.47$ kg m$^{-3}$.

The pressure signal is measured near the engine at $z_{1k} = L/2 - 1.27 = 2.67$ m from the middle of the resonator (figure 2). At the measurement location, the pressure fundamental amplitude is estimated by fitting the entire pressure signal with a Fourier series. Then, maximum pressure amplitude $P_1$ is calculated at the end of the resonator (standing wave in perfect fluid hypothesis):

$$P_1 = \frac{P_{1k}}{\cos\left(\frac{2\pi f_1}{c_0} z_{1k}\right)} \quad (1)$$

with $f_1$ the fundamental frequency, $\omega = 2\pi f$ the angular frequency and $c_0$ the sound speed; here, we have $P_1 = 127$ mbar. This corresponds to a drive ratio $P_1/P_m = 1.75\%$. The measured frequency of the acoustic wave is $f_1 = 21.92$ Hz (resonance frequency of the tube with stack and heat exchangers) and $c_0 = 347.28$ m s$^{-1}$ (adiabatic sound speed). Up to the eighth harmonic of $f_1$ are present in the pressure power spectrum (figure 3).

The velocity fields are measured by PIV in a small pipe section located at 1.16 m from the left end of the resonator. It is constituted of a cylindrical tube of glass 0.1 m in
<figure><img src="image_50.png" /><figcaption>**Figure 4.** Diagram of PIV and pressure acquiring system.</figcaption></figure>

length, with an anti-reflection coating. Seeding particles of magnesium oxide (MgO), less than 1 µm in diameter, are injected into the pipe close to the measurement section using a Dantec spreading device, limited to a maximum pressure of 10 bar. A buffer tank allows seeding of the resonator up to 20 bar. The particles have a sedimentation velocity between 1 and 6.6 mm min⁻¹ (respectively for diameter of 0.4 µm and 1 µm). The light source is a 10 Hz double cavity pulsed YAG of 532 nm wavelength and 30 mJ/pulse. Particle images are recorded with a 20 Hz CCD camera (1024 × 768 pixels). The sampling double frame rate is $f_e = 10$ Hz with a delay of 2.000 ms between the two images of one pair. The time record of one series is close to 100 s. The data acquisition system (SigLab) has two parallel synchronized channels, one for the acoustic pressure and the other for the cavities laser Q-switch control, with the signal sample frequency at 5120 Hz. An anti-aliasing filter is activated only for the pressure channel (see figure 4).

Seven different measurement series were performed with the same control parameters, allowing verification of the reproducibility of the results. Two data sets of measurements are carried out with two days lag, respectively labelled A (A1–A3) and B (B1–B4). The first record A1 starts 12.5 min after

a first seeding injection, followed by A2 and A3 carried out respectively 2.5 min and 7.5 min after a second injection. The same procedure is applied for the four records of set B with only one seeding injection.

The velocity field ensues from image pairs processing in two steps: image pre-processing and particle displacement computation. Firstly, the optical reflections on the pipe walls are removed and the ground noise reduced (figure 5(a)). Usually this is obtained by subtracting from each frame a reference one, acquired before seeding. In our case, this solution is not well adapted because the seeding particles induce a strong additional light source that modifies reflections. Consequently, subtraction of an image averaged over the entire series is preferred (figure 5(b)), suppressing most of the reflections. The image contrast is finally improved by applying a gamma factor to the gray level (see figure 5(c)).

Then, the nonhomogeneous intensity of images is corrected by applying a high-pass filter and the intensity mean levels between the images are equalized.

The second step addresses the particle displacement fields computed from a pair of images as usual, but with an original optical flow algorithm using orthogonal dynamic programming (Quénot *et al* 1998). Dynamic programming algorithms are generally used for optimization (Horowitz *et al* 1996). This advanced PIV method is based on an iterative search for a global alignment of particle images, with regularity and continuity constraints (minimization of a Minkowski norm). The particle image displacement problem inside a plane is reduced to a one-dimension problem, dividing the images in strips along one direction and considering the displacement in the same direction. The two-dimensional (2D) solution is obtained by successive iterative resolution within smaller and smaller orthogonal strips. When the algorithm converges a displacement vector for each pixel is obtained. This algorithm presents some decided advantages (Stanislas *et al* 2005) when looking for a wide velocity range, compared with classical or even advanced correlation methods of PIV:

- a vector field without outliers (<0.02%), saving their identification and the interpolation artefacts for their replacement,

- the velocity field is by construction continuous and differentiable yielding an accurate and reliable estimate where strong velocity gradients occur (Stanislas *et al* 2008), in particular inside the boundary layer if the reflection problem can be overcome,

<figure><img src="image_51.png" /><figcaption>**Figure 5.** (a) Unprocessed image with particles; (b) reference average image; (c) processed image after subtraction.</figcaption></figure>
<figure><img src="image_52.png" /><figcaption>Figure 6. Velocity vector field $\vec{U}_{ac}(x, y, t)$ and the limited region.</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>Figure 7. Evolution of axial component of space average acoustic velocity $u_x(t_n)$ versus acquisition time $t$.</figcaption></figure>

- a high spatial resolution of one vector per pixel and an accuracy comparable to the best advanced algorithms based on cross-correlation using adaptive meshes.

- very good stability, especially in noisy images processing with low spatial definition.

The algorithm presents several disadvantages:

- sensitivity to local out-of-plane velocity component in the two-component version of the algorithm (Stanislas et al 2005);

- adjustment of many parameters and their sensitivity to flow configuration.

The optical flow resolution with the image analysis setting used in the present work is 1/181 pixel. That resolution corresponds to a velocity accuracy of $4 \times 10^{-4} \text{ m s}^{-1}$.

The cylindrical curvature of the measurement region deforms vertical displacements and positions of the particles. This dioptre effect is corrected with an optical model, detailed in appendix B.

The image processing presented above gives instantaneous velocity fields, strongly undersampled in time. In the following section, a first estimate of the acoustic field is obtained from this set of undersampled data.

### 3. Characterization of the first-order acoustic field

The height of the domain of measurements is the entire diameter of the waveguide and the axial length represents 0.63% of the wavelength and 2.5% of the axial length of the expected Rayleigh streaming cells. The velocity fields obtained in the previous section are referred to as $\vec{U}(x, y, t)$, $x \in [0, 10 \text{ cm}]$ and $y \in [-28.15 \text{ mm}, 28.15 \text{ mm}]$ being the axial and radial local coordinates of the PIV measurement zone and $t$ the acquisition time, respectively.

Consider the decomposition $\vec{U}(x, y, t) = \vec{U}_{ac}(x, y, t) + \vec{U}_{ac}(x, y, t)$. The time mean velocity is first estimated over a long integration time $t_{tot}$, $\vec{U}(x, y) = \frac{1}{t_{tot}} \int_{0}^{t_{tot}} \vec{U}(x, y, t) dt$. The acoustic velocity is obtained by subtraction:

$$ \vec{U}_{ac}(x, y, t) = \vec{U}(x, y, t) - \vec{U}(x, y). $$

Then, a limited region is selected, situated far from the boundary layers (figure 6), over which space variations are considered negligible. In this region, spatial average of the axial acoustic velocity component $\vec{U}_{ac}$ is performed yielding a new discrete time set $u_x(t_n)$.

The time dependent set $u_x(t_n)$ is represented in figure 7 for one measurement series. This time behavior is characteristic of an undersampled periodic signal. The class distribution of that series (figure 8) is also typical of a mono-periodical signal
<figure><img src="image_54.png" /><figcaption>**Figure 8.** Velocity class distribution.</figcaption></figure>

<table><caption>**Table 1.** Values of velocity extrema and amplitudes (in m s⁻¹) for the experimental data of the seven measurement series and their mean values over all the series from A1 to B4 (in bold type).</caption><thead><tr><th rowspan="2">Velocities (m s<sup>-1</sup>)</th><th colspan="8">Name of experimental series</th></tr><tr><th>A1</th><th>A2</th><th>A3</th><th>B1</th><th>B2</th><th>B3</th><th>B4</th><th>Mean</th></tr></thead><tbody><tr><td>Experimental</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Maximum</td><td>1.647</td><td>1.663</td><td>1.649</td><td>1.699</td><td>1.684</td><td>1.690</td><td>1.685</td><td><strong>1.672</strong></td></tr><tr><td>Minimum</td><td>-2.064</td><td>-2.103</td><td>-2.083</td><td>-2.142</td><td>-2.124</td><td>-2.128</td><td>-2.131</td><td><strong>-2.107</strong></td></tr><tr><td>Amplitude</td><td>1.856</td><td>1.883</td><td>1.866</td><td>1.921</td><td>1.904</td><td>1.909</td><td>1.908</td><td><strong>1.890</strong></td></tr><tr><td>Analytical amplitude</td><td>1.867</td><td>1.892</td><td>1.870</td><td>1.923</td><td>1.901</td><td>1.908</td><td>1.904</td><td><strong>1.894</strong></td></tr><tr><td>Amplitude difference<br>(% of the analytical)</td><td>-0.012</td><td>-0.009</td><td>-0.004</td><td>-0.002</td><td>0.003</td><td>0.001</td><td>0.004</td><td><strong>-0.004</strong></td></tr></tbody></table>

with some noise. The two extrema of the set are distributed over several classes and not in the most probable class, as for a pure cosine. Therefore, their minimum $u_x(t_n)_{min}$ and maximum $u_x(t_n)_{max}$ values must be averaged at least on the two last classes.

The amplitude of $u_x(t_n)$ is then roughly given by $(u_x(t_n)_{max} - u_x(t_n)_{min})/2$.

The estimated amplitudes of the seven series are compared with the corresponding analytical values obtained from the theoretical acoustic field (Hamilton *et al* 2003):

$$u_z(z, y) = v_0 \left( \frac{\cosh(\alpha z)}{\cosh(\alpha L/2)} - 1 \right) R_z(y). \quad (2)$$

In the previous equation, $u_z$ represents the axial component of acoustic velocity inside the resonator, $z \in [-L/2, L/2]$ represents the axial coordinate in the resonator, $v_0$ is the amplitude of the harmonic excitation of the resonator (control parameter), $\alpha$ and $R_z(y)$ are intermediate functions depending on the thermal-viscous properties of the gas, the characteristics of the wave and the Bessel functions $J_0$ and $J_1$:

$$\alpha^2 = - \left( \frac{1 + (\gamma - 1) f_{\kappa}}{1 - f_{\nu}} \right) \frac{\omega^2}{c^2} \\ f_{\nu,\kappa} = \frac{2J_1(\beta_{\nu,\kappa} r_0)}{\beta_{\nu,\kappa} r_0 J_0(\beta_{\nu,\kappa} r_0)}, \quad \beta_{\nu,\kappa} = \frac{i-1}{\delta_{\nu,\kappa}}$$

$$\delta_{\nu} = \sqrt{\frac{2\nu}{\omega}}, \quad \delta_{\kappa} = \sqrt{\frac{2\kappa}{\omega}}$$

$$R_z(r) = 1 - \frac{J_0(\beta_{\nu} r)}{J_0(\beta_{\nu} r_0)}$$

$$R_r(r) = \frac{r}{r_0} - \frac{J_1(\beta_{\nu} r)}{J_1(\beta_{\nu} r_0)} + \frac{(\gamma - 1)(1 - f_{\nu}) f_{\kappa}}{[1 + (\gamma - 1) f_{\kappa}] f_{\nu}} \left( \frac{r}{r_0} - \frac{J_1(\beta_{\kappa} r)}{J_1(\beta_{\kappa} r_0)} \right).$$

The control parameter $v_0$ is not directly accessible by experiment but can be calculated from the maximal pressure defined in equation (1), knowing the pressure $P_{1k}$ at the end of the resonator ($P_{1k}$ is measured with a pressure transducer as shown in figure 2), $v_0 = \frac{P_{1k}\alpha \cosh(\alpha \frac{\gamma}{2})}{i\omega\rho_0 \sinh(\alpha z_{1k})}$. Table 1 summarizes the measured acoustic wave characteristics and the analytical values obtained with equations (2), for the seven records.

The reproducibility of the results is good whatever the series, validating the measurement system and the data treatments. In particular, these results confirm that the initial flow conditions, induced by the seeding injection process, vanish rapidly, as concluded from direct visualization. Moreover, experimental and analytical values are in close agreement, their difference lying between 1 to 4 mm s⁻¹. The velocity extrema point out a large asymmetry. The difference between maximum and minimum velocity is about 20.7% of the minimum velocity. This means that, within an acoustic cycle, the flow coming from the thermoacoustic engine is faster than the flow going back to the thermoacoustic engine. This asymmetry is generated by the harmonics presence as confirmed by the power spectrum for one series (B1), plotted in figure 9.

The velocity temporal set being undersampled, the spectrum is aliased. The fundamental frequency of the wave $f_1 = 21.89$ Hz appears in this spectrum at $f_1 = 1.892$ Hz (one can verify that $1.892$ Hz + $2 \times 10$ Hz = $21.89$ Hz). Similar expressions can be found for the harmonics $f_2 = 43.78$ Hz to $f_6 = 131.34$ Hz. The same frequencies can be observed in the pressure power spectrum (figure 3). The local rate of the first harmonic (ratio of the amplitude of the first harmonic to the fundamental) is evaluated with the Fourier spectrum to reach 19.3%. Asymmetrical wave extrema indicate a phase lag of harmonics with respect to the fundamental mode, which cannot be evaluated at this stage. The wave appears as a mix
<figure><img src="image_55.png" /><figcaption>Figure 9. Aliased power spectrum PS (fourth root) for the axial velocity component; $f_1$ is the fundamental mode frequency and $f_2$ to $f_6$ are the frequencies of the five first harmonics.</figcaption></figure>

of a standing-wave component and a low-rate travelling wave component.

To go further in the wave analysis, a precise wave cycle reconstruction is needed.

## 4. Phase rearrangement of the undersampled data

Let us recall that quantifying the continuous secondary flow $\vec{U}(x, y)$ is this paper's main goal, with $\vec{U}(x, y) = \vec{U}(x, y, t) - \vec{U}_{ac}(x, y, t)$. As previously seen, the PIV fields are strongly undersampled (approximately one snapshot every two periods). Moreover, the periodic component is at least two orders of magnitude greater than the searched steady component $\vec{U}(x, y)$. The separation can therefore be achieved only by knowing $\vec{U}_{ac}(x, y, t)$ very precisely over one cycle of the lowest fundamental frequency $f_1$ (see figure 9). Usually such a rearrangement is done by phase averaging based on an external clock, which can be obtained by filtering a well-sampled pressure signal around the fundamental frequency. But in the nonlinear range this procedure presents some limitations, particularly when the pressure reference has to be taken far from the velocity field measurement. The shift lag between the two measurements can induce correlation losses or phase shift not necessarily frequency independent. We propose here a direct identification of the torus associated with the fundamental frequency, from the undersampled time series itself. Therefore, a phase portrait is built with an embedding technique: the signal is projected on its principal components using SVD (Broomhead and King 1986, Mees et al 1987). The time set is reorganized over one cycle using the phase order.

### 4.1. Self-referenced acoustic velocity cycle

The results presented above show that the time set $u_x(t_n)$ is mainly a periodic signal with a fundamental frequency $f_1$ and its harmonics (figure 9), which built the acoustic wave. In order to analyse the acoustic wave, the set $u_x(t_n)$ is rearranged over one cycle of $f_1$ using self-referenced phase ordering. The periodic part of the dynamics is obtained by projecting the most energetic mode on the first two principal components of

<figure><img src="image_56.png" /><figcaption>Figure 10. Phase portrait ($\xi_{svd}^{(1)}(t)$; $\xi_{svd}^{(2)}(t)$) obtained from SVD on 99 modes; $\varphi_{svd}$ represents the current phase.</figcaption></figure>

the signal. To this end, the SVD is applied to the delay matrix (called 'trajectory matrix' by Broomhead and King (1986)), denoted by $S$ and constructed with the discrete time set $u_x(t_n)$:

$$S = \begin{pmatrix} u_x(t_0) & u_x(t_1) & \dots & u_x(t_{M_{SVD}-1}) \\ u_x(t_1) & \dots & \dots & \dots \\ \vdots & \vdots & \ddots & \vdots \\ u_x(t_{r-1}) & u_x(t_r) & \dots & u_x(t_{N_{max}-1}) \end{pmatrix}$$

In the previous definition, $M_{SVD}$ is the number of SVD decomposition modes, $r = N_{max} - (M_{SVD} - 1)$, $N_{max}$ being the total number of velocity data and $t_n = \frac{n}{f_e}$, $1 \le n \le N_{max} - 1$.

$S$ is decomposed according to the relation:

$$S = UVD^t,$$

$D$ being the diagonal matrix of the singular values and $V$ the unitary matrix of eigenvectors. The phase space is generated with a basis provided with the columns of $\Phi_{SVD}$, defined as:

$$\Phi_{SVD} = UD = SV.$$

The $j$th column of $\Phi_{SVD}$ represents the $j$th mode $\xi_{svd}^{(j)}(t)$, ranked from the most energetic mode to the lowest. The first two modes $\xi_{svd}^{(1)}(t)$ and $\xi_{svd}^{(2)}(t)$ are used to construct a phase portrait, which is the projection of the dynamic on its fundamental periodic part. But this phase portrait is not sufficiently circular, because SVD decomposition modes 1 and 2 contain a harmonic residue. If the phase portrait is not a perfect circle, the velocity reconstruction will yield a deformed profile. This will have a consequence on the mean velocity estimate, as detailed in appendix A. Therefore, a band-pass filtering centered on the aliased fundamental acoustic frequency is applied to $\xi_{svd}^{(1)}(t)$ and $\xi_{svd}^{(2)}(t)$, in order to filter harmonics. Then the phase portrait is constructed (figure 10). Let us remark that the band-pass filter cannot be applied directly on the time set, because the aliasing associated with the undersampling is uncontrolled.
<figure><img src="image_57.png" /><figcaption>Figure 11. Phase rearrangements of velocity $u_x(\varphi_{svd})$ (·) and of pressure $P_1(\varphi_{svd})$ (+) measured at $P_{1k}$, using SVD decomposition based on velocity field.</figcaption></figure>

A phase $\varphi_{svd}$, defined as the phase of the complex number $\xi_{svd}^{(1)} + i\xi_{svd}^{(2)}$ (figure 10), can be assigned to each instantaneous field velocity $u_x(t_n)$ and therefore to the whole corresponding velocity field $\bar{U}(x, y, t_n)$. Velocity data are then rearranged in phase on the interval $[-\pi; \pi]$ according to $\varphi_{svd}$ values. Pressure data measured near the engine (position $P_{1k}$ in figure 2) are also rearranged on the interval $[-\pi; \pi]$ according to $\varphi_{svd}$ values. In figure 11 are represented both the velocity $u_x(\varphi_{svd})$ and the pressure $P_1(\varphi_{svd})$. The phase origin cannot be determined from the phase portrait and is arbitrarily taken in figure 11 as being $-\pi$. This phase rearrangement shows a sinusoidal signal, slightly deformed by the harmonics. The asymmetry values of the extrema are close to those found in section 3.

A validation of the embedding procedure is done using the oversampled pressure record $P_1(t)$. An undersampled pressure set is constructed linking the pressure data to their corresponding velocity one. Hence, we obtain an undersampled pressure signal, like the velocity data. The phase order obtained previously is used to rearrange this undersampled pressure series. In figure 12 are plotted the acoustic pressure cycle $P_1(\varphi_{svd})$ together with a representative period $[-\pi; \pi]$ of the original pressure signal. The good agreement validates the phase ordering.

## 4.2. Standing-wave properties

The undersampled data being rearranged over one period, information on acoustics and streaming can be extracted. In this section, local and global properties of acoustic wave are successively obtained.

Harmonic content of velocity cycle is estimated identifying the signal with a sinusoidal function, having one fundamental oscillating mode and its first harmonic:

$$v_{\text{fit}}(\varphi_{\text{svd}}) = a_{0}^{(\nu)} + a_{1}^{(\nu)} \cos(d_{1}^{(\nu)}\varphi_{\text{svd}} + \psi_{1}^{(\nu)}) \\ + a_{2}^{(\nu)} \cos(d_{2}^{(\nu)}\varphi_{\text{svd}} + \psi_{2}^{(\nu)}).$$

<figure><img src="image_58.png" /><figcaption>Figure 12. Comparison between the oversampled pressure data for one period $P_1(t)$ (—) and the undersampled data rearranged by SVD $P_1(\varphi_{svd})$ (+).</figcaption></figure>

In this equation, $a_0^{(\nu)}$ represents the time averaged value, $a_i^{(\nu)}$ and $\psi_i^{(\nu)}$, respectively, the amplitude and the phase of mode 'i'. The parameters $d_1^{(\nu)}$ and $d_2^{(\nu)}$ are considered unknown in order not to impose the wave shape. The time average being already removed from the signal before the phase rearrangement, the constant $a_0^{(\nu)}$ should be useless. Nevertheless, many calculations have established that the identification does not run properly when this term is left out. The main reason is that the exact mean value will only be accessible at the end of the procedure (see section 5.1). The results averaged over all series A1–A3, B1–B4 are summarized in table 2.

The fundamental mode is identified through the coefficient $d_1^{(\nu)} = 1.0$ and its first harmonic with $d_2^{(\nu)} = 2.0$, with an accuracy of 5%. The fundamental mode amplitude is $1.816 \text{ m s}^{-1}$, slightly smaller than the amplitude obtained with the histograms method ($1.890 \text{ m s}^{-1}$, cf table 1). One
<table><caption>**Table 2.** Averaged coefficients obtained by identification of the velocity cycle rearranged in phase.</caption><thead><tr><th colspan="6">Coefficients of fit</th></tr><tr><th>a<sub>0</sub><sup>(v)</sup> (m s<sup>-1</sup>)</th><th>a<sub>1</sub><sup>(v)</sup> (m s<sup>-1</sup>)</th><th>a<sub>2</sub><sup>(v)</sup> (m s<sup>-1</sup>)</th><th>|a<sub>2</sub><sup>(v)</sup>/a<sub>1</sub><sup>(v)</sup>| (%)</th><th>d<sub>1</sub><sup>(v)</sup></th><th>d<sub>2</sub><sup>(v)</sup></th></tr></thead><tbody><tr><td>0.030</td><td>1.816</td><td>0.339</td><td>18.7</td><td>1.0</td><td>2.0</td></tr></tbody></table>

<table><caption>**Table 3.** Coefficients obtained by identification of the pressure cycle rearranged in phase.</caption><thead><tr><th colspan="6">Coefficients of fit</th></tr><tr><th>a<sub>0</sub><sup>(p)</sup> (bar)</th><th>a<sub>1</sub><sup>(p)</sup> (bar)</th><th>a<sub>2</sub><sup>(p)</sup> (bar)</th><th>|a<sub>2</sub><sup>(p)</sup>/a<sub>1</sub><sup>(p)</sup>| (%)</th><th>d<sub>1</sub><sup>(p)</sup></th><th>d<sub>2</sub><sup>(p)</sup></th></tr></thead><tbody><tr><td>0.001</td><td>0.110</td><td>0.006</td><td>5.8</td><td>1.0</td><td>2.0</td></tr></tbody></table>

explanation of this difference is that the histogram method overestimates the velocity amplitude, by integrating all modes. The harmonic amplitude is about 34 cm s⁻¹ and its rate with respect to the fundamental is 18.7%, in good agreement with the result obtained from the folded spectrum (19.3% harmonic rate, figure 9). This high rate is then explained by the measurement zone location, situated near the harmonic maximum amplitude. The phase lag between the fundamental and harmonic |ψ₁⁽⁽⁾ - ψ₂⁽⁽⁾| varies significantly from one series to another and is responsible for the velocity extrema asymmetry.

The same approach is used for the pressure cycle, which is identified with the function

$$P_{\text{fit}}(\varphi_{\text{svd}}) = a_{0}^{(P)} + a_{1}^{(P)} \cos(d_{1}^{(P)} \varphi_{\text{svd}} + \psi_{1}^{(P)}) + a_{2}^{(P)} \cos(d_{2}^{(P)} \varphi_{\text{svd}} + \psi_{2}^{(P)}).$$

The results averaged over all series A1–A3, B1–B4 are summarized in table 3.

Again the identified values of the parameters $d_1^{(v)}$ and $d_2^{(v)}$ indicate that the two first modes match the fundamental and the first harmonic modes. The first harmonic rate with respect to the fundamental is 5.8%, much smaller than the velocity harmonic rate as a consequence of the pressure sensor location in the standing wave (0.5 m from the wave generator, close to fundamental pressure antinodes). The cycle is then less deformed.

Now, one can estimate the phase-lag between fundamental modes of pressure and velocity. For all the measurement series A and B, the same value is obtained, $|\psi_1^{(P)} - \psi_1^{(v)}| = 0.48\pi$, which is close to the theoretical value of $0.50\pi$, specific for a standing wave inside a resonator.

In the following, the global space evolution of the wave (wavelength, maximum amplitude, etc) is deduced from the undersampled velocity data, measured on a very small fraction (0.63%) of the wavelength.

The initial axial velocity data are used again, but averaged only in the radial direction y. A new space-time set is obtained, depending on x and t, denoted by $u_x(x, t) = \langle U_x(x, y, t) \rangle_y$. For each location x in the measurement zone, the time average component is subtracted from the previous set $u_x(x, t)$. Then, SVD phase rearrangement is done on the new velocity set after subtraction, leading to the velocity profiles $u_x(\varphi_{svd})$ over one

<figure><img src="image_59.png" /><figcaption>**Figure 13.** y-averaged velocity profiles evolution according to x and phases $\varphi_{svd}$.</figcaption></figure>

<table><caption>**Table 4.** Identification coefficients for data $u_x(\varphi_{svd})$.</caption><thead><tr><th colspan="6">Coefficients of fit</th></tr><tr><th>u<sub>m0</sub> (m s<sup>-1</sup>)</th><th>u<sub>m1</sub> (m s<sup>-1</sup>)</th><th>u<sub>m2</sub> (m s<sup>-1</sup>)</th><th>|u<sub>m2</sub>/u<sub>m1</sub>| (%)</th><th>λ<sub>exp</sub> (m)</th><th>λ<sub>th</sub> (m)</th></tr></thead><tbody><tr><td>0.004</td><td>4.046</td><td>0.438</td><td>10.9</td><td>15.56</td><td>15.84</td></tr></tbody></table>

cycle. In figure 13 are represented a large number of those profiles according to their phases, on about 8.5 cm length, corresponding to the measurement area.

For each phase lag, the x-evolution is typical of a standing wave. Despite the length of the measurement window being about 1% of the resonator length, a proper identification can be achieved using the whole space {$x, \varphi_{svd}$}x∈[0,0.1], $\varphi_{svd}∈[-\pi, \pi]$ information. Therefore, supposing a two-frequency plane standing wave, the fitted velocity depending on the global axial coordinate $z \in [-\frac{L}{2}; \frac{L}{2}]$, is given by:

$$u_{\text{fit}}(z, \varphi_{\text{svd}}) = u_{m0} + u_{m1} \sin \left[ k \left( z + \frac{L}{2} \right) \right] \cos(\varphi_{\text{svd}} + \psi_1) + u_{m2} \sin \left[ 2k \left( z + \frac{L}{2} \right) \right] \cos(2\varphi_{\text{svd}} + \psi_2). \quad (3)$$

The unknowns of the fit procedure are the amplitudes ($u_{m0}, u_{m1}, u_{m2}$), the phases ($\psi_1, \psi_2$) and the wave number $k$. The identification is done for each measurement series A and B and the coefficients are averaged over all series (table 4).

The global harmonic rate $u_{m2}/u_{m1}$ is about 11%. The wavelength estimate of the fundamental mode is $\lambda_{exp} = 15.56$ m. This experimental estimation can be compared to $2L$ (L the tube length) which yields 15.76 m or with $\lambda_{th} = c_0/f_1 = 15.84$ m ($f_1$ measured fundamental frequency and $c_0$ the sound wave velocity). The latter is preferred, closer to reality of the system (resonator with thermoacoustic stacks and heat exchangers) in which the temperature varies significantly (Bétrancourt 2008). The difference with the experience is then 1.77% (0.28 m) of the theoretical wavelength, which represents a slight difference when knowing that $\lambda_{exp}$ is estimated on a pipe segment which represents only 0.63% of the wavelength.

Finally, the maximum velocity amplitude measured by PIV at the location $z_{piv} = L/2+1.16$ is compared with the
estimated one at the same location using equation (3). One finds 2.172 m s⁻¹, in quite good agreement with the maximum velocity amplitude 1.894 m s⁻¹ estimated by the histogram method (section 3). The phase values $\psi_1 = 1.12\pi$ and $\psi_2 = 0.95\pi$ are obtained through the fit.

After this careful characterization of the acoustic velocity, it becomes possible to evaluate the mean velocity field.

# 5. Streaming estimation

## 5.1. Phase average method

As already done in section 3 to extract the oscillations of the acoustic part, the mean velocity component calculation could be performed by a simple time average, approximated in the discrete form by:

$$ \bar{U}(\vec{x}) \approx \frac{1}{N_{\max}} \sum_{n=0}^{N_{\max}} \vec{U}(x_i, y_j, t_n). $$

This zero order discrete averaging form presents several limitations that penalize an exact partition between the strong periodic oscillation and the tiny steady component. Practically, the most challenging point comes from the slow convergence to zero of the periodic residue when the number $N_{\max}$ of velocity fields increases. With the considered time series, the discrete average approximation approaches the exact time average with an error of the same magnitude as the term which has to be estimated yielding high imprecision. To illustrate that, let us estimate the residue of the average calculation for the velocity axial component. Consider the Hussain and Reynolds (1970) decomposition:

$$ v_x(\vec{x}, t) = \tilde{u}_x(\vec{x}, t) + u'_x(\vec{x}, t) + \tilde{u}_x(\vec{x}). $$

In this equation, $\tilde{u}_x(\vec{x}, t)$ is the periodic part (including the fundamental acoustic mode and harmonics of the wave) such as $\tilde{u}_x(\vec{x}, t + T) = \tilde{u}_x(\vec{x}, t)$, T being the fundamental period of the oscillation and $\vec{x} = (x, y)$. $\tilde{u}_x(\vec{x})$ is the time-independent component (time average) and $u'_x(\vec{x}, t)$ the non-periodic velocity component (turbulence, noise, etc). Over an infinite time, the mean value of the whole time dependent part of the velocity field $\tilde{u}_x(\vec{x}, t) = \tilde{u}_x(\vec{x}, t) + u'_x(\vec{x}, t)$ is by definition equal to zero:

$$ \bar{u}(\vec{x}, t) = \lim_{t_{\text{tot}} \to \infty} \frac{1}{t_{\text{tot}}} \int_{0}^{t_{\text{tot}}} \tilde{u}(\vec{x}, t) \, dt = 0 $$

$$ \overline{u'_x(\vec{x}, t)} = 0. $$

Let us estimate the maximum error related to an average over a finite time interval [0, $t_{\text{tot}}$]. Consider $t_{\text{tot}} \gg 1$ and $N_T$ the number of intervals [$kT$, $(k+1)T$] included in [0, $t_{\text{tot}}$], with T the fundamental period of the oscillation. We shall leave out the spatial dependence in $\vec{x}$:

$$ \begin{aligned} \frac{1}{t_{\text{tot}}} \int_{0}^{t_{\text{tot}}} v_x(t) \, dt &= \frac{1}{t_{\text{tot}}} \left( \int_{0}^{t_{\text{tot}}} [\tilde{u}_x(t) + u'_x(t) + \tilde{u}_x] \, dt \right) \\ &= \frac{1}{t_{\text{tot}}} \left( \int_{0}^{t_{\text{tot}}} u'_x(t) \, dt + \int_{0}^{N_T T} \tilde{u}_x(t) \, dt \right. \\ & \qquad \left. + \int_{N_T T}^{t_{\text{tot}}} \tilde{u}_x(t) \, dt \right) + \tilde{u}_x. \end{aligned} $$

The periodicity of $\tilde{u}_x$ implies $\int_{0}^{N_T T} \tilde{u}_x(t) \, dt = 0$. Furthermore, the time mean value of $u'(t)$ converges to zero faster than the residual integral $\frac{1}{t_{\text{tot}}} \int_{N_T T}^{t_{\text{tot}}} \tilde{u}_x(t) \, dt$, so that one obtains at the first order of approximation:

$$ \frac{1}{t_{\text{tot}}} \int_{0}^{t_{\text{tot}}} v_x(t) \, dt = \bar{u}_x + \underbrace{\frac{1}{t_{\text{tot}}} \int_{N_T T}^{t_{\text{tot}}} \tilde{u}_x(t) \, dt}_{(1)} $$

The residual integral (1) can be estimated by changing the time variable $t$ into the phase variable $\varphi$. Consider the following transformation:

$$ t = \frac{T}{2\pi}\varphi + N_T T, \quad \varphi \in [0, 2\pi]. $$

The residue becomes then

$$ \frac{t}{t_{\text{tot}}} \int_{N_T T}^{t_{\text{tot}}} \tilde{u}_x(t) \, dt = \frac{T}{2\pi t_{\text{tot}}} \int_{0}^{2\pi \left( \frac{t_{\text{tot}}}{T} - N_T \right)} \tilde{u}_x(\varphi) \, d\varphi. $$

Suppose now that $\tilde{u}_x(\varphi)$ has the simplified form $2 \sin(\varphi)$ (one frequency wave with an amplitude close to the experiment), and consider the case the most unfavourable where the integral has a maximal contribution, corresponding to $2\pi (\frac{t_{\text{tot}}}{T} - N_T) = \pi$. In our experiment, $T = 4.56 \times 10^{-2}$ s, $t_{\text{tot}} = 100$ s, $N_T = 0(000)$ and the residual part estimate is:

$$ \frac{T}{2\pi t_{\text{tot}}} \int_{0}^{\pi} 2 \cos(\varphi) \, d\varphi = 2.9 \times 10^{-4} \, \text{m s}^{-1}. $$

Knowing that the expected mean velocity can be around $10^3$ smaller than the acoustic component, this residue can represent 30% of the mean velocity. Consequently, the averaging of the velocity field has to be performed over an integer number of reconstructed cycles using the SVD phase rearrangement discussed in section 4.1.

This being established, two other sources of errors must be reduced. Firstly, the average integral is done using a first-order scheme, because the adaptive Simpson quadrature does not improve the result here. The velocity field averaged with respect to the phase is calculated for the 2D velocity field at each point of the measurement zone:

$$ \bar{U}(\vec{x}) = \frac{1}{2\pi} \int_{-\pi}^{\pi} \bar{U}(\vec{x} \cdot \varphi_{\text{svd}}) \, d\varphi_{\text{svd}}. $$

Phases are deduced from the phase portrait $(\xi_{\text{svd}}^{(1)}; \xi_{\text{svd}}^{(2)})$ (figure 10). As mentioned in section 4.1, both SVD modes are filtered around the fundamental frequency to avoid harmonic pollution. Appendix A highlights this problem with a synthetic signal, composed of a standing wave in $(x, t)$ (containing a fundamental mode and its first harmonic) superimposed on a parabolic y-profile representing a streaming mean flow. We observed that filtering the harmonic into the two SVD modes before rearrangement gives a better approximation.

## 5.2. Mean velocity field

In order to test the reproducibility of the phase mean velocity estimation, we compare this average result for all measurement series A and B. The comparison is performed in a limited region of the measurement zone, to avoid the loss of image quality coming from the boundaries. We notice quite good stability of this mean flow, whatever the series, despite distinct
<figure><img src="image_60.png" /><figcaption>**Figure 14.** Mean velocity profiles averaged over y axis $\langle \bar{u}_x(x, y) \rangle_y$ (left), $\langle \bar{u}_y(x, y) \rangle_y$ (right) for all six series.</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>**Figure 15.** Phase mean velocity field (limited zone) for the series: A2 (left) and B1 (right).</figcaption></figure>

measurement times after particle injection according to series (see section 2). In order to compare the time averaged velocity field, the variance between the six records A1, A2, A3, B1, B2, B3 is computed at each spatial location and then averaged over the whole space. Finally, the square root is computed. This deviation is normalized with the mean value, i.e. $\sqrt{\frac{\langle \text{var}(\bar{u}_{x,i}(x, y)) \rangle_{x,y}}{\langle \text{mean}(\bar{u}_{x,i}(x, y)) \rangle_{x,y}}}$, yielding 3% maximum variation for the axial velocity component $\bar{u}_x(x, y)$ and 7% for the radial component $\bar{u}_y(x, y)$. It is also worthwhile to compare velocity profiles versus $x$ averaged over $y$. For the axial velocity component (figure 14 left), the agreement is very good, keeping in mind that the acoustic velocity field is 100 times stronger than the mean flow. For the radial component (figure 14 right), the discrepancy between the six series is still acceptable because the same trend is observed even if the two first series A1 and B1 are less similar. As an example, figure 15 presents two mean velocity fields, obtained in the limited region (defined in figure 6) for the series A2 and B1.

The reproducibility being quantified, let us consider now the whole measured field, excepted for the very near wall region corrupted by reflection. Figure 16 presents the mean

<figure><img src="image_62.png" /><figcaption>**Figure 16.** Phase mean velocity field (PIV zone) for the series B1.</figcaption></figure>

field for the whole PIV measurement zone obtained for the series B1. In the following, only this series is treated, the fields from the other series being similar.

The flow has a complicated structure that can be separated into three regions:
(i) One central region with small velocities probably due to the presence of deformed vortices, corresponding to a stagnation boundary between two streaming cells, eventually accompanied by a transverse flow.

(ii) A right side region in which the flow coming from the engine moves towards the left side near the upper edge; in the middle, it goes down again to head for the engine.

(iii) A left side region, roughly symmetric to the right part, but with the presence of a little vortex at the bottom corner.

Globally, it seems there are two large and two small counter-
rotating three-dimensional cells exactly at the position of
the measurement zone. Arroyo and Greated (1991) noticed
the existence of a slight transverse velocity in the case of the
Rayleigh streaming that stands for 11% of the maximum
Rayleigh streaming velocity. This could explain the transverse
drift of the seeding particles in the laser sheet. In addition, in
the present work, cells are strongly asymmetric. This flow
is different from that expected for a regular axisymmetric
Rayleigh streaming in terms of patterns. Nevertheless, the
velocity values of this field are between 0 cm s⁻¹ and
5 cm s⁻¹, which is of the same order of magnitude as the
Rayleigh theoretical estimation: 2.2 cm s⁻¹ in our experimental
conditions. The observed cell asymmetry could be due to a
nonlinear regime of streaming. This regime is characterized
by the critical values of the streaming Reynolds number.
Two definitions of the Reynolds number are generally used.
According to Stuart (1966), $Re_s = \frac{u_s \lambda}{v}$ where $u_s$ is the
streaming velocity, $\lambda$ the wavelength and $v$ the kinematic
viscosity. For $Re_s > 1$, inertial effects deform the streaming
cells. In our case, $u_s = 2.2$ cm s⁻¹ and $Re_s$ is of the
order of 166 000. According to Menguy and Gilbert (2000),
the nonlinear Reynolds number $Re_{nl} = \frac{U_{\max}^2}{\frac{1}{2}\left(\frac{\delta u}{R_c}c_0\right)^2}$, based
on the maximal amplitude of acoustic velocity in the center of
the resonator, separates the linear and nonlinear regime for a
critical value $Re_{nl} = 1$. In the present case, it reaches the value
$Re_{nl} = 9.7$, which means once again that streaming cells would
be deformed by the inertial effects. As $Re_{nl}$ increases above
one, some authors indicate the possible presence of more than
four outer streaming cells in a half-wavelength resonator, still
due to inertial effects (Nabavi et al 2009, Reyt et al 2012, Daru
et al 2013). The 2D velocity field in our measurements seems
to correspond to this case.

Let us examine the transverse evolution of axial velocity
profiles $\bar{u}_x(y)$, associated with the vector field in figure 16.
These profiles are represented in figure 17 for different local
axial positions, from $x = 0.4$ mm to 9.1 cm (in steps of 1.1 mm).
Three values of $x$ are displayed in the graph to indicate the axial
evolution of the velocity.

One can observe three different flow regions, delimited
by $y = -2$ mm and $y = 1.5$ mm, with a strong asymmetry,
as already observed from the 2D velocity field. The flow near
boundaries ($y > 1.5$ cm and $y < -2$ cm) is roughly opposite
to that existing in the zone $y \in [-2; 1.5]$ cm. The same
observation can be made for the left/right side regions. All
this indicates the presence of four counter-rotating convection
cells between the core of the resonator and the wall boundaries,
and also between left and right, confirming the conclusions
of the 2D flow analysis. The flow close to the tube walls

<figure><img src="image_63.png" /><figcaption>**Figure 17.** Velocity profiles $\bar{u}_x(y)$ in m s⁻¹ plotted for different x values (in steps of 1.1 mm) for the series B1.</figcaption></figure>

does not approach zero as expected for a boundary layer
region (of 0.17 mm) because the measurements are polluted by
reflections as already explained previously. Nevertheless, the
global evolution of velocity profiles near the walls indicates
the existence of strong mean flows in this region.

In conclusion, analysis of experimental results shows that
the measurement area is located at the boundary between
four irregular streaming cells, indicating the existence of
multiple streaming vortices inside the resonator. The flow
regime observed is then far above the linear case of Rayleigh
streaming cells.

**6. Conclusion**

This work is dedicated to the study of a thermoacoustic system reproducing the working conditions of a real machine: the low frequency wave is generated with a thermoacoustic engine, which allows high intensity waves, under large mean pressure (above 7 bars). This makes our configuration different from other academic studies on acoustic streaming. Thus, one has to deal with the metrological challenge of optical measurements in realistic systems. Moreover, first-order (acoustic oscillation) and second-order quantities (streaming motion) have to be estimated at different time scales (acoustic fast scale and slow streaming scale).

All measurements are performed at fixed drive-ratio in
order to validate the method. The velocity field measurements
are provided by PIV in a small segment of a long resonator. The
data, recorded on a long time scale, are strongly undersampled.
The acquisition rate of particle images is low (10 Hz) in
relation to the frequency wave (~22 Hz). The velocity fields
are computed with an optical flow technique, which allows
accuracy of $3 \times 10^{-4}$ m s⁻¹. The undersampled velocity data
are rearranged with respect to their phase over one period
[-$\pi$; $\pi$], with Bromhead and King's embedding technique
based on SVD. Then, information on acoustics and streaming
<figure><img src="image_64.png" /><figcaption>**Figure 18.** Diagram of the processing methodology in ten steps: 1, velocity field snapshots input (cf section 2, figure 6); 2, space averaged field reduced to a temporal set (cf section 3); 3, temporal set of velocity fluctuation (cf section 3); 4, embedding of the dynamics into a phase-space (cf section 4.1); 5, band pass filtering centered on the fundamental frequency (cf section 4.1); 6, 2D phase portrait projection and phase ordering (cf section 4.1); 7, rearrangement of the velocity field input according to phase ordering (cf section 4.1); 8, 9, acoustic wave properties (cf section 4); 10, phase-averaged velocity field (cf section 5).</figcaption></figure>

are extracted. The methodology analysis developed in this paper is summarized in the diagram below (figure 18).

First, PIV velocity measurements allow characterization of the first-order flow. The fundamental acoustic component (of 1.890 m s⁻¹ at the PIV location) is identified in agreement with the analytic amplitude (1.894 m s⁻¹). The first harmonic rate amplitude is 19% of the fundamental amplitude at the PIV location. The phase lag between the first harmonic and the fundamental mode is measured, confirming the existence of a travelling wave component. The comparison with the acoustic pressure cycle shows a phase lag velocity–pressure of 0.48 π, close but not equal to the theoretical value specific for standing waves. The wavelength is estimated from measurements taken on a small distance representing 0.63% of this wavelength, yielding a difference of 1.77% with respect to the theoretical value.

Finally, the time-mean velocity is estimated, using a phase average over the period [-π; π]. This velocity of the order of 10⁻² m s⁻¹ corresponds to the streaming motion. Several calculations over different sets of data show good reproducibility of the mean measurement. Results show the existence of irregular streaming, of the same order of magnitude as Rayleigh streaming but with different space patterns. We can see a complex flow with several counter-rotating vortices at the measurement location, indicating the existence of multiple streaming cells inside the resonator, with an axial periodicity higher than quarter of the wavelength. This flow differs from the linear Rayleigh–Schlichting streaming. This is because the inertia effects of the oscillating flow

cannot be neglected, the Re$_{nl}$ characteristic of this experiment being higher than one. In the literature, there is no model available in the parameter range corresponding to our chosen configuration. However, recent experiments and numerical simulations in the nonlinear regime of streaming show the existence of multiple distorted streaming cells (Nabavi et al 2009, Daru et al 2013, Reyt et al 2012).

From the methodology point of view, a new procedure is proposed, that overcomes the undersampling of acoustic oscillation for evaluating both acoustics and streaming. It is a general method that can be used to study the mean flow evolution when varying parameters in a range beyond existing theoretical models. In particular, when increasing acoustic wave intensity, the irregular streaming regime appears, under the influence of inertia effects as well as thermal effects, which are not yet entirely understood.

## Appendix A. Harmonic influence on SVD embedding

In order to show the influence of the harmonic in SVD mode decomposition, a synthetic sinusoidal standing wave is used, which contains a fundamental frequency of 21.95 Hz and a first harmonic with an amplitude rate of 18% from the fundamental mode, each mode having its own phase. The characteristic parameters of the wave are chosen to fit those of experiment. The Rayleigh–Nyborg 2D analytical solution for streaming in large channels bounded by two parallel solid
<figure><img src="image_65.png" /><figcaption>Figure A1. Evolution of signal s(x_c, y_c, t) versus phase for a period; (x_c, y_c) is the central point of the observation zone.</figcaption></figure>

walls is superimposed on the standing wave:

$$s(x; y; t) = u_m(x, y) + u_A \cos(2\pi ft + \psi_1) \sin\left(\frac{2\pi}{\lambda}x\right) + u_A 0.18 \cos(4\pi ft + \psi_2) \sin\left(\frac{4\pi}{\lambda}x\right).$$

Here $u_m(x, y)$ is the Rayleigh–Nyborg streaming velocity, given by

$$u_m(x, y) = \frac{u_A^2}{8c} \sin(2kx)[F_1(y) + 3 - 18F_2(y)],$$

with

$$F_1(y) = -\left[e^{-2y/\delta_v} + 2e^{-y/\delta_n} \cos\left(\frac{y}{\delta_v}\right) + 6e^{-y/\delta_v} \sin\left(\frac{y}{\delta_v}\right)\right], \\ F_2(y) = \frac{y}{2y_0}\left(1 - \frac{y}{2y_0}\right).$$

Here, $u_A = 3.8 \text{ m s}^{-1}$, $x \in [1.115 \text{ m}; 1.205 \text{ m}]$, $y \in [0; 2y_0]$, $y_0 = 56.5 \times 10^{-3} \text{ m}$, $f = 21.95 \text{ Hz}$, $\psi_1 = 1.12\pi$ and $\psi_2 = 0.95\pi$.

The synthetic signal $s(x, y, t)$ is oversampled at the frequency 21950 Hz. In figure A1 the phase/time evolution of $s(x_c, y_c, t)$ is represented over one acoustic period, ($x_c, y_c$) being the center of the observation zone. In figure A2 the time mean profile of $s(x, y, t)$ is shown at $x = x_c$.

With the previous function, a synthetic signal is created, undersampled in time, at the frequency 10 Hz over 100 s, like in PIV measurements. The SVD is applied on temporal set $s(t) = \langle s(x, y, t) - \tilde{s}(x, y, t) \rangle_{x,y}$, as explained in section 4.1. The average with respect to the phase is then calculated for the entire velocity field. Figure A3 presents: (i) the phase average evolution with respect to y-coordinate of the axial velocity $s_m(y) = s(x_c, y, t)$, at $x = x_c$, (ii) the streaming Rayleigh–Nyborg velocity $u_m(x, y)$ and (iii) the simple time average (over $t_{tot} = 100 \text{ s}$) of $s(x_c, y, t)$. A difference of 3.8 mm s$^{-1}$ between the phase average velocity and the Rayleigh–Nyborg solution can be observed at the maximum value of the profile, representing an error of 62% of the maximum Rayleigh–Nyborg solution. The simple time

<figure><img src="image_66.png" /><figcaption>Figure A2. The axial streaming Rayleigh–Nyborg velocity profile versus y-coordinate.</figcaption></figure>

<figure><img src="image_67.png" /><figcaption>Figure A3. Axial mean velocity y-profile for the synthetic signal with harmonic: (i) —, $s_m(y)$ (phase average), (ii) ---, time average and (iii) +, Nyborg solution.</figcaption></figure>

<figure><img src="image_68.png" /><figcaption>Figure A4. Aliased power spectrum PS of $\xi_{svd}^{(1)}$ and $\xi_{svd}^{(1)}$ SVD modes in semi-log scale.</figcaption></figure>

average velocity also underestimates the exact value, the error being smaller (34% of the maximum).

Analysis of the aliased power spectrum of the two first SVD modes, plotted in figure A4, shows the fundamental acoustic mode (at $f_1 = 1.949 \text{ Hz}$) and a residue of the
<figure><img src="image_69.png" /><figcaption>**Figure A5.** Axial mean velocity y-profile for the synthetic signal without harmonic: (i) —, $s_m(y)$ (phase average), (ii) ---, time average and (iii) +, Nyborg solution.</figcaption></figure>

<figure><img src="image_70.png" /><figcaption>**Figure A6.** Axial mean velocity y-profile for the synthetic signal with harmonic: (i) —, $s_m(y)$ (phase average with filtered harmonic), (ii) ---, time average and (iii) +, Nyborg solution.</figcaption></figure>

harmonic (at $f_2 = 3.898$ Hz). Even though its intensity is small (0.002% of fundamental), this harmonic residue is responsible for deforming the phase portrait so that the average calculation is biased. This is confirmed performing the same procedure on the previous undersampled signal but without the harmonic part. The results show very good agreement between the phase average velocity and the Rayleigh–Nyborg solution (see figure A5). The maximum difference is 0.02 mm s⁻¹, representing 0.3% of the maximum Rayleigh–Nyborg solution, while the time average velocity still underestimates the exact value (error estimation is 25%).

In the case of a signal with harmonic content, a great part of the harmonic is contained in the third and fourth SVD modes. However, increasing the number of SVD modes does not allow suppression of the harmonic pollution of the two first modes. Instead, a simple band-pass filter applied to the two first SVD modes $\xi_{svd}^{(1)}$ and $\xi_{svd}^{(2)}$ (three-order Butterworth filter) significantly reduces the estimate error (figure A6). The maximum difference between phase-averaged velocity and the Rayleigh–Nyborg solution is 0.7 mm s⁻¹ (representing 11% of the maximum Nyborg’s solution).

<figure><img src="image_71.png" /><figcaption>**Figure B1.** Schema of different trajectories of the diffracted light through the glass tube of inner radius $R_i$ and glass thickness $e$.</figcaption></figure>

## Appendix B. Optical correction for dioptre effect

The cylindrical curvature of the measurement region deforms vertical displacements and positions of the particles. This dioptre effect is corrected with an optical model (Faure et al 2006) applied to the present configuration. In figure B1 is represented a schema of a plane orthogonal to the resonator axis and to the laser sheet plane. The light diffracted by particles crosses three different media before reaching the camera lens: the working gas (nitrogen), the glass wall (circular cylinder) and the air, outside the tube.

Here, $r$ represents the current distance from the center, and $r_a$ its apparent distance on the camera image. The refractive indices of the three media are respectively $n_N$, $n_{glass}$, $n_{air}$. The beam deviations at the nitrogen/glass and glass/air interfaces are characterized by the angles $\alpha_{1,2,3}$ and $\beta_{1,2}$.

Using Snell–Descartes equations, the dioptre effect correction is found to be:

$$r = \frac{n_{\text{air}}}{n_{N2}} \frac{r_a}{\cos(\alpha_3 - \beta_1 + \beta_2 - \alpha_2)}$$

## References

* Adrian R J 1979 Conditional eddies in isotropic turbulence *Phys. Fluids* **22** 2065–70

* Aktas M and Farouk B 2004 Numerical simulation of acoustic streaming produced by finite-amplitude resonant oscillations in an enclosure *J. Acoust. Soc. Am.* **116** 2822–31

* Arroyo M P and Greated C A 1991 Stereoscopic particle image velocimetry *Meas. Sci. Technol.* **2** 1181–6

* Bailliet H, Gusev V, Raspet R and Hiller R A 2001 Acoustic streaming in closed thermoacoustic devices *J. Acoust. Soc. Am.* **110** 1808–21

* Bétrecourt A 2008 Simulations numériques et validations expérimentales de machines de réfrigération thermoacoustiques *Dissertation Université Pierre et Marie Curie (Paris)*

* Broomhead D S and King G P 1986 Extracting qualitative dynamics from experimental data *Physica D: Nonlinear Phenom.* **20** 217–36

* Campbell M, Cosgrove J A, Greated C A, Jack S and Rockliff D 2000 Review of LDA and PIV applied to the measurement of sound and acoustic streaming *Opt. Laser Technol.* **32** 629–39
Daru V, Balteanu Carlès D, Weisman C, Debesse P and Gandikota G
2013 Two-dimensional numerical simulations of nonlinear
acoustic streaming in standing waves *Wave Motion*
**50** 955–63

Debesse P 2008 Vers une mesure du vent thermoacoustique
*Dissertation* Université Pierre et Marie Curie (Paris)

Desjouy C, Pénelet G, Lotton P and Blondeau J 2009 Measurement of acoustic streaming in a closed-loop traveling wave resonator using laser Doppler velocimetry *J. Acoust. Soc. Am.* **126** 2176–83

Duffourd S 2001 Réfrigérateur thermoacoustique: études analytiques et expérimentales en vue d'une miniaturisation
*Dissertation* Ecole Centrale de Lyon, Lyon

Faure T M, Lusseyran F, Gougat P and Launay F 2006 Experimental and numerical investigation of the flow distribution inside a tubular heat exchanger *J. Fluids Eng.* **128** 1218–27

Hamilton M F, Ilinskii Y A and Zabolotskaya E A 2003 Thermal effects on acoustic streaming in standing waves *J. Acoust. Soc. Am.* **114** 3092–101

Horowitz E, Sahni S and Rajasekaran S 1996 *Computer Algorithms in C++* (New York: Freeman)

Hussain A K M F and Hayakawa M 1987 Eduction of large-scale organized structures in a turbulent plane wake *J. Fluid. Mech.* **180** 193–229

Hussain A K M F and Reynolds W C 1970 The mechanics of an organized wave in turbulent shear flow *J. Fluid. Mech.* **41** 241–58

Kiya M and Matsumara M 1985 Turbulence structure in intermediate wake of a circular cylinder *Bull. Japan Soc. Mech. Eng.* **28** 2617–24

Konstantinidis E, Balabani S and Yianneskis M 2005 Conditional averaging of PIV plane wake data using a cross-correlation approach *Exp. Fluids* **39** 38–47

Lourenco L, Subramanian S and Ding Z 1997 Time series velocity field reconstruction from PIV data *Meas. Sci. Technol.* **8** 1533–8

Mees A I, Rapp P E and Jennings L S 1987 Singular-value decomposition and embedding dimension *Phys. Rev. A* **36** 340–6

Menguy L and Gilbert J 2000 Non-linear acoustic streaming accompanying a plane stationary wave in a guide *Acustica* **86** 249–59

Michard M, Blanc-Benon P, Grosjean N and Nicot C 2004 Apport de la véclocimétrie par images de particules pour la caractérisation du champ de vitesse acoustique dans une maquette de réfrigérateur thermoacoustique 9ème Congrès Francophone de Vélocimétrie Laser (Brussels, 14–17 Sept.)

Moreau S, Bailliet H and Valière J-C 2008 Measurements of inner and outer streaming vortices in a standing waveguide using laser Doppler velocimetry *J. Acoust. Soc. Am.* **123** 640–7

Nabavi M, Siddiqui M H K and Dargahi J 2007 Simultaneous measurement of acoustic and streaming velocities using

synchronized PIV technique *Meas. Sci. Technol.*
**123** 1811–7

Nabavi M, Siddiqui M H K and Dargahi J 2009 Analysis of regular and irregular acoustic streaming patterns in a rectangular enclosure *Wave Motion* **46** 312–22

Negri M, Cozzi F and Malavasi S 2011 Self-synchronized phase averaging of PIV measurements in the base region of a rectangular cylinder *Meccanica* **46** 423–35

Pastur L R, Lusseyran F, Faure T M, Fraigneau Y, Pethieu R and Debesse P 2008 Quantifying the non-linear mode competition in the flow over an open cavity at medium reynolds number *Exp. Fluids* **44** 597–608

Perrin R, Braza M, Cid E, Cazin S, Barthet A, Sevrain A and Mockett C 2007 Phase-averaged measurements of the turbulence properties in the near wake of a circular cylinder at high Reynolds number by 2C-PIV and 3C-PIV *Exp. Fluids* **43** 341–55

Quénot G, Pakleza J and Kowalewski T A 1998 Particle image velocimetry with optical flow *Exp. Fluids* **25** 177–89

Rayleigh L 1884 On the circulations of air observed in Kundt's tubes, and on some allied acoustical problems *Phil. Trans. R. Soc. Lond.* **175** 1–21

Reyt I, Moreau S, Bailliet H and Valière J C 2012 Experimental study of acoustic streaming in a high level standing wave guide: Influence of mean temperature and higher harmonics distribution *Nonlinear Acoustics: State-of-the-Art and Perspectives (ISNA 19)*, AIP Conf. Proc. **1474** 83–6

Schlichting H 1932 Berechnung ebener periodischer grenzschichtströmungen (calculation of plane periodic boundary layer streaming) *Z. Phys.* **33** 327–35

Stanislas M, Okamoto K, Kähler C J, Westerweel J and Scarano F 2008 Main results of the third international PIV challenge *Exp. Fluids* **45** 27–71

Stanislas M, Okamoto K, Kähler C J and Westerweel J 2005 Main results of the second international PIV challenge *Exp. Fluids* **39** 170–91

Stuart J T 1966 Double boundary layer in oscillatory viscous flow *J. Fluid Mech.* **24** 673–987

Swift G W 2002 *Thermoacoustics: A Unifying Perspective for Some Engines and Refrigerators* (Melville, NY: Acoustical Society of America)

Thompson M W and Atchley A A 2005 Simultaneous measurement of acoustic and streaming velocities in a standing wave using laser doppler anemometry *J. Acoust. Soc. Am.*
**117** 1828–38

Thompson M W, Atchley A A and Maccarone M J 2005 Influences of a temperature gradient and fluid inertia on acoustic streaming in a standing wave *J. Acoust. Soc. Am.* **117** 1839–49

Widmann A, Duchmann A, Kurz A, Grundmann S and Tropea C 2012 Measuring Tollmien–Schlichting waves using phase-averaged particle image velocimetry *Exp. Fluids* **53** 707–15
O. Hireche, C. Weisman, D. Baltean-Carlès, P. Le Quéré and L. Bauwens

Low Mach number analysis of idealized thermoacoustic engines with numerical solution.
J. Acoust. Soc. Am. 128(6), 3438-3448, (2010).
# Low Mach number analysis of idealized thermoacoustic engines with numerical solution

Omar Hireche

Laboratoire de Mécanique Appliquée, Université des Sciences et de la Technologie d'Oran Mohamed Boudiaf,
BP 1523, EL-M'naouer Oran 31000, Algeria

Catherine Weisman<sup>a)</sup> and Diana Baltean-Carlès<sup>b)</sup>

Université Pierre et Marie Curie, 4 Place Jussieu 75252 Paris Cedex 05, France

Patrick Le Quéré

Laboratoire d'Informatique pour la Mécanique et les Sciences de l'Ingénieur, Centre National de la
Recherche Scientifique UPR 3251, 91403 Orsay Cedex, France

Luc Bauwens

Department of Mechanical and Manufacturing Engineering, University of Calgary, Calgary, Alberta T2N
1N4, Canada

(Received 4 June 2009; revised 29 September 2010; accepted 4 October 2010)

A model of an idealized thermoacoustic engine is formulated, coupling nonlinear flow and heat exchange in the heat exchangers and stack with a simple linear acoustic model of the resonator and load. Correct coupling results in an asymptotically consistent global model, in the small Mach number approximation. A well-resolved numerical solution is obtained for two-dimensional heat exchangers and stack. The model assumes that the heat exchangers and stack are shorter than the overall length by a factor of the order of a representative Mach number. The model is well-suited for simulation of the entire startup process, whereby as a result of some excitation, an initially specified temperature profile in the stack evolves toward a near-steady profile, eventually reaching stationary operation. A validation analysis is presented, together with results showing the early amplitude growth and approach of a stationary regime. Two types of initial excitation are used: Random noise and a small periodic wave. The set of assumptions made leads to a heat-exchanger section that acts as a source of volume but is transparent to pressure and to a local heat-exchanger model characterized by a dynamically incompressible flow to which a locally spatially uniform acoustic pressure fluctuation is superimposed. © 2010 Acoustical Society of America. [DOI: 10.1121/1.3506348]

PACS number(s): 43.35.Ud [RR]

Pages: 3438–3448

## I. INTRODUCTION

The nonlinear nature of the energy exchanges within the stack plays a crucial role in thermoacoustic engines, affecting power generation, development of thermoacoustic instability, frequency selection, and saturation. Classical linear theory<sup>1–3</sup> has been successful in predicting onset temperatures and stability curves for each acoustic mode that are in good agreement with experiments.<sup>4–7</sup> Quasi-one-dimensional linear codes<sup>8</sup> yield parameters such as temperature profile, frequency, power of the saturated regime for a fixed geometry, working gas, and mean pressure. Subsequent work based upon linear theory includes, for instance, approximate relations for the thermoacoustic gain and efficiency, damping coefficient, oscillation frequency, and onset temperature.<sup>9,10</sup> So far, extension to nonlinear regimes has either used simplified quasi-one-dimensional models<sup>11,12</sup> or if based upon multidimensional models,<sup>13,14</sup> these rely upon restrictive assumptions, or a rather coarse spatial resolution in the

stack and heat exchangers. The current work presents a well-resolved model, in which the flow in the stack and heat exchangers remains fully nonlinear and multidimensional. Resonator and load are described by a standard linear one-dimensional acoustic model. Proper coupling between these two submodels results in a global model that is asymptotically consistent in the small Mach number approximation.

The classical thermoacoustic stability analysis by Rott<sup>1</sup> was based upon the standard assumptions leading to linear acoustics, i.e., a formulation in the small Mach number limit, with length and time in a ratio of the order of the speed of sound, which he also coupled with a boundary layer model. Rott proceeded by considering the limit case whereby temperature changes are reduced to a jump with zero thickness between two uniform values and introducing *ad-hoc* continuity conditions across that singularity. In order to obtain this limit solution, he thus considered a temperature gradient taking place over a length that becomes negligible compared with the wavelength, as the Mach number approaches zero.

In the “short engine” model, Swift<sup>2</sup> also assumed that the stack is short compared with the wavelength but long compared with the sweep by the fluid motion. Indeed, in the resulting model, pressure is time-dependent but with negligible phase shift over the stack length, while a linear

<footnote><sup>a)</sup>Author to whom correspondence should be addressed. Present address: LIMSI, CNRS UPR 3251, 91403 Orsay Cedex, France. Electronic mail: catherine.weisman@upmc.fr</footnote>

<footnote><sup>b)</sup>Current address: LIMSI, CNRS UPR 3251, 91403 Orsay Cedex, France.</footnote>
formulation is retained, hence the implicit assumption
whereby the stack is longer than the sweep.

In the current model, the system is split into two different
regions. The length of the stack and heat exchangers assem-
bly is assumed comparable with the sweep by the fluid
motion. While often realistic, this scaling is interesting
because the resulting local model for the heat exchanger as-
sembly retains all the physics included in Swift's short engine
model to which it adds the physics of nonlinear advection (at
the cost of additional complications). In the resonators, the
usual assumption is made of a length of the order of the
wavelength. With the local flow in the heat exchangers thus
adequately modeled, consistent jump conditions across Rott's
singularity are obtained, which, coupled with an acoustic for-
mulation in the resonator, result in an asymptotically consist-
ent global model. There is no *ad-hoc* assumption in this
approach and no phenomenological imposed law.

The model can deal with complex and realistic heat
exchanger and stack configurations, and it retains the essential
adptive nonlinearity. Although that nonlinearity remains
fully reversible, it does result in the stack and heat exchangers
behaving as a source in the acoustics, with a strength that will
be reduced as the amplitude grows. It is thus suitable for tran-
sient operation and if the computation is pursued long enough
it should result in a saturation and hence stationary operation.

Different scalings characterize different physics, which
opens the door to a variety of singular perturbation or multi-
ple scales studies. There are situations that require using a
single length scale and multiple time scales, such as slowly
varying acoustics.<sup>11,15</sup> Here, the presence of two distinct
regions described by different length scales while sharing a
single characteristic time lead to a singular perturbation prob-
lem (i.e., matched asymptotics or boundary layer methods).

The goal here is to demonstrate the effectiveness of the
singular perturbation analysis. In doing so, the emphasis was
placed on the flow in the stack and heat exchangers, which
was coupled with the simplest meaningful consistent resona-
tor model. On that basis, an exploration of the fundamental
consequences of the current set of scaling assumptions is
presented, but no attempt is made to obtain realistic results
that could be used for practical purposes. To illustrate the
power of the method, results were obtained focusing mainly
upon startup and stability and also saturation due to nonlin-
ear effects in the heat exchangers and stack, in a much ideal-
ized resonator and load configuration.

The two assumptions proposed for the stack and heat
exchangers, of a small Mach number and a sweep of the
order of the stack and heat-exchangers length, lead to a local
model in which spatially uniform pressure fluctuations can
still actually occur at up to leading order. Large temperature
and density variations can also be dealt with, both in time
and in space.¹⁶ (If additionally the solution domain would
be connected to a large pressure reservoir such as the
atmosphere, and if the boundary conditions would restrict
temperature variations to some small values, the standard
incompressible Navier-Stokes equations would be obtained;
however, these additional assumptions are not applicable in
the current context.) The resulting model is more or less
standard in Stirling engine and pulse-tube cryorefrigerator

models, in which leading order pressure varies,¹⁷ in contrast
to the current situation, where due to the presence of the res-
onators, pressure fluctuations only occur at acoustic order. It
has been used in an analysis in straight ducts,¹⁸ and simula-
tions motivated by thermoacoustic devices, either using a
vortex method¹⁹ or second-order accurate finite volume
method,²⁰ that dealt with the heat exchangers alone, with
boundary conditions meant to represent resonator acoustics.

In the current formulation, the computational fluid dy-
namics (CFD) solution for the local flow in the stack and
heat exchangers is coupled with a consistent acoustic model
in the resonators. Other recent CFD studies have been based
upon a fully viscous, conductive, and compressible solu-
tion.²¹ Some have relied upon commercial codes to simulate
the thermoacoustic engine.²² Given the broad range of scales
present in the fully compressible problem, these approach
leads to extremely large simulations²³ or simulation of devi-
ces of very small (miniature) size.²⁴ In contrast, the current
approach allows focusing a smaller computational effort
upon the heat exchangers and stack, where it is important
because this is where small scales truly matter.

II. PHYSICAL MODEL

A. Geometry

The geometry of the device consists of a duct with rec-
tangular cross-section within which a set of heat exchangers
is placed. The set of heat exchangers is made up of a stack
placed between a cold and a hot heat exchanger. This assem-
bly is inserted within the resonator tube, which is closed at
one end and terminates with a load at the other end, which
thus undergoes a small amplitude oscillation. Dead space
separates the stack from the heat exchangers. In this article,
the expression "heat-exchanger section" is consistently used
to refer to the entire region, including heater, stack, and
cooler. The geometry is shown in Fig. 1 (top).

B. Formulation

The problem that is being addressed consists of an initial
boundary value problem, on the spatial domain shown in
Fig. 1 (top), characterized by the conservation laws of vis-
cous, conducting, compressible fluid dynamics, together
with conservation of energy in the stack walls. The working
fluid is a known ideal gas with constant specific heats.

Boundary conditions include no-slip conditions on the
heat exchanger and stack walls. Resonator walls are free slip
and adiabatic. A temperature boundary condition is provided
on the heat-exchanger walls. Both temperature and heat flux
are continuous on the stack walls. The load is characterized
by a known relationship between the local acoustic pressure
and velocity, which is approximately the same at the mean
load location and the actual load location, assuming its dis-
placement is negligible compared with the wavelength.
These initial conditions are discussed in detail in Sec. II G.

The problem thus described can be solved numerically
without further simplification. However, this entails very
large computations, while, clearly, as long as the Mach num-
bers remain small, flow in the resonators will differ very lit-
tle from linear acoustics. In contrast, in the heat-exchanger
<figure><img src="image_72.png" /><figcaption>FIG. 1. Geometry of entire system (top) and of numerical simulation domain (bottom).</figcaption></figure>

section, there is a need for a simulation that deals with transient and spatially resolved velocities and temperatures. Next, an appropriate simplified model is derived from asymptotic analysis, reducing the flow in the resonator to an acoustic model, eliminating the need for an acoustic Couplant-Friedrichs-Lewy condition (time steps smaller than mesh size divided by the speed of sound) and focusing on the heat-exchanger section.

### C. Scaling assumptions

Since temperatures must be properly resolved throughout the device, temperatures, hence the speed of sound, must be taken to be of the same order throughout. Although the presence of a load requires that one of the extremities of the resonator must be moving, the displacement will be small as long as fluid velocities are small. Thus one can introduce a reference resonator length $L_{res}$, a reference temperature, which here, is set equal to the cold heat-exchanger temperature $T_{ref}$, and consequently, a reference speed of sound $c_{ref}$. Given that the period will be determined by an acoustic resonance, time is scaled by $\tilde{t}_{ref} = L_{res}/c_{ref}$, which equals one half of the resonant period in an empty tube closed at both ends and at the cold end temperature. This provides a single time scale for the entire problem.

It is reasonable to expect the heat-exchanger section to be acoustically compact, i.e., its characteristic length is small compared with the resonator length. If the heat exchangers are not overly large compared to the stack, the length of the heat-exchanger section is of the same order as the stack length, which provides a suitable reference, and $L_{stack} \ll L_{res}$.

As to fluid velocities, very low initial values associated with initial conditions are not particularly representative and do not provide a meaningful scale. However, close to the stationary regime, the fluid motion should sweep a length comparable with the stack length, over which the temperature gradient takes place, hence velocities of order $L_{stack}/\tilde{t}_{ref}$. Taking into account the definition of $\tilde{t}_{ref}$ above results in $L_{stack}/L_{res} = U_{ref}/c_{ref}$, so that

$$M = \frac{\tilde{U}_{ref}}{c_{ref}} = \frac{\tilde{L}_{stack}}{\tilde{L}_{res}} \ll 1. \qquad (1)$$

Pressure is scaled by its value at rest $\tilde{P}_{ref}$ and density by its value at reference pressure and temperature $\tilde{\rho}_{ref} = \tilde{P}_{ref}/RT_{ref}$. All variables are scaled uniformly throughout the device except for the longitudinal coordinate $\tilde{x}$, which is scaled by $\tilde{L}_{stack}$ in the heat-exchangers section but by $L_{res}$ in the resonators. Perturbation series in the reference Mach number are introduced for all dependent variables, using the notation, for example, for pressure

$$p = p^{(0)} + Mp^{(1)} + M^2 p^{(2)} + O(M^3). \qquad (2)$$

For velocity, however, the solution below only includes the leading order contribution so there is no need to introduce extra symbols. Below, $u$ is used for the dimensionless longitudinal velocity and $v$ for the transverse velocity; $\mathbf{u}$ represents the dimensionless velocity vector.

If the length of the heat-exchanger section is smaller than the resonator length by a factor of order $M$, then it only contributes to $O(M)$ on the overall length. The symbol $\tilde{x}$ is used for the longitudinal coordinate as scaled in the resonators and setting its origin at the heat-exchangers location. Then at the closed end on the hot side, $\tilde{x} = -l_L = -\tilde{l}_L/\tilde{L}_{res}$. Likewise, at the mean load location, $\tilde{x} = l_R = \tilde{l}_R/\tilde{L}_{res}$. Both the stack plate thickness and the distance between plates are of the order of the viscous and heat penetration depth. In the heat-exchanger section, the dimensionless spatial coordinates scaled by the stack length are denoted as $x$ and $y$, with the origin set at the stack center.

Next, the formulation is written in the dimensionless variables (without $\tilde{\cdot}$) introduced above, respectively, in the resonators and then in the heat-exchanger section.

### D. Resonators

If the resonator cross-section is much wider than the viscous penetration depth, the wall boundary layer thickness is negligible compared to the resonator width. Furthermore, if the thermal boundary conditions are adiabatic, the flow is one-dimensional. Likewise, longitudinal diffusion results in negligible irreversibility. Thus, the flow is approximately isentropic, and as long as the Mach number remains small, temperature, density, and pressure remain constant at leading
order, hence equal to their spatially uniform initial values: $p^{(0)} = 1$ and $\rho^{(0)} = 1/T^{(0)}$. The temperature is $T^{(0)} = T_{\text{cold}} = 1$ on the cold side while it equals $T^{(0)} = T_{\text{hot}} = \tilde{T}_{\text{hot}}/\tilde{T}_{\text{cold}}$ on the hot side. At leading order, these assumptions lead to linear acoustics. Indeed, the mass and momentum conservation equations become

$$ \frac{\partial \rho^{(1)}}{\partial t} + \rho^{(0)} \frac{\partial u}{\partial \hat{x}} = 0, \qquad \rho^{(0)} \frac{\partial u}{\partial t} + \frac{1}{\gamma} \frac{\partial p^{(1)}}{\partial \hat{x}} = 0, \qquad (3) $$

in which $\rho^{(1)}$ and $p^{(1)}$ are related by the isentropic relation: $p^{(1)} = \gamma T^{(0)} \rho^{(1)}$, where $\gamma$ is the ratio of specific heats.

Boundary conditions at the two resonator ends correspond, respectively, to a closed end at $\hat{x} = -l_L$ and a load at the other extremity, instantaneously located at $\hat{x}_R = l_R + O(M)$, and characterized by some differential relationship between velocity and pressure, which becomes linear in the small Mach number limit. A purely dissipative load is then represented by a reflection coefficient, applied at the fixed location $\hat{x} = l_R$ and denoted by $f$ in dimensionless form. (The current formulation can readily be extended to more realistic loads.) Thus,

$$ u(-l_L, t) = 0, \qquad p^{(1)}(l_R, t) = fu(l_R, t). \qquad (4) $$

The acoustic problem depends upon continuity or jump conditions at the heat-exchanger location, yet to be determined. However, based upon the current information, the problem can be reduced to a single relationship between acoustic pressure and velocity at the location $\hat{x} = 0$, on each side, based upon a d'Alembert solution on characteristics.

Eliminating density,

$$ \frac{\partial p^{(1)}}{\partial t} + \gamma \frac{\partial u}{\partial \hat{x}} = 0, \qquad \gamma \frac{\partial u}{\partial t} + T^{(0)} \frac{\partial p^{(1)}}{\partial \hat{x}} = 0. \qquad (5) $$

This system is readily diagonalized, adding and subtracting the second equation to the first one multiplied by $\sqrt{T^{(0)}}$, yielding the Riemann variables, which remain constant on characteristics moving at the speed of sound $\sqrt{T^{(0)}}$, respectively, left and right:

$$ \mathcal{L} = \gamma u - \sqrt{T^{(0)}} p^{(1)}, \qquad \mathcal{R} = \gamma u + \sqrt{T^{(0)}} p^{(1)}. \qquad (6) $$

At the hot (left) side, the boundary condition is $u = 0$, which, since $\mathcal{L} + \mathcal{R} = 2\gamma u$ determines a reflection whereby $\mathcal{R}_L(-l_L, t) = -\mathcal{L}_L(-l_L, t)$. Since acoustic waves move at the speed of sound, $\mathcal{L}(\hat{x}, t) = \mathcal{L}(0, t + \hat{x}/\sqrt{T^{(0)}})$ and $\mathcal{R}(\hat{x}, t) = \mathcal{R}(0, t - \hat{x}/\sqrt{T^{(0)}})$, so that $\mathcal{R}_L(0, t + l_L/\sqrt{T^{(0)}}) = -\mathcal{L}_L(0, t - l_L/\sqrt{T^{(0)}})$. Shifting time by $2l_L/\sqrt{T^{(0)}}$, at the heat-exchanger location, $\hat{x} = 0^{-}$, taking into account that $T^{(0)} = T_{\text{hot}}$, one finds

$$ \mathcal{R}_L(0^{-}, t) = -\mathcal{L}_L\left(0^{-}, \frac{t - 2l_L}{\sqrt{T_{\text{hot}}}}\right). \qquad (7) $$

Likewise, on the cold side, taking into account that $\mathcal{L} + \mathcal{R} = 2\gamma u$ and $\mathcal{R} - \mathcal{L} = 2\sqrt{T^{(0)}}p^{(1)}$, and in view of the boundary condition of Eq. (6), at the heat-exchangers location $\hat{x} = 0^{+}$, and since $T^{(0)} = 1$,

$$ \mathcal{L}_R(0^{+}, t) = Z \mathcal{R}_R(0^{+}, t - 2l_R), \qquad (8) $$

in which

$$ Z = \frac{\gamma - f}{\gamma + f}. \qquad (9) $$

An open end corresponds to $f = 0$ hence $Z = 1$ while at a closed end, $f$ is infinite and $Z = -1$. Both cases result in a zero load since the power developed is the time-averaged value of $up^{(1)}$.

## E. Heat-exchanger section

Within the heat-exchanger section, the problem is not amenable to analytical solution. A numerical simulation is required.

The stack is made up of a large number of identical uniformly spaced plates. It is thus reasonable to assume vertical periodicity, i.e., flow and heat transfer are identical in all gaps. The flow in only one slot between two contiguous plates needs being simulated, as shown in Fig. 1 (bottom). The same scaling as described above is used but, in contrast with the resonator model, length is now scaled by $\tilde{L}_{\text{stack}}$. The perturbation series equations, Eq. (2), are then introduced in the compressible, viscous, and conductive Navier-Stokes equations, respectively, conservation of mass, momentum, and energy. The system equations become

$$ \frac{\partial \rho^{(0)}}{\partial t} + \nabla \cdot (\rho^{(0)} \mathbf{u}) = 0, \qquad (10) $$

$$ \left\{ \begin{array}{l} \nabla p^{(0)} = 0, \\ \nabla p^{(1)} = 0, \\ \frac{\partial(\rho^{(0)} \mathbf{u})}{\partial t} + \nabla \cdot (\rho^{(0)} \mathbf{u} \otimes \mathbf{u}) = -\nabla p^{(2)} + \frac{1}{\text{Re}} \nabla \cdot \tau, \end{array} \right. \qquad (11) $$

$$ \rho^{(0)} \left[ \frac{\partial T^{(0)}}{\partial t} + (\mathbf{u} \cdot \nabla) T^{(0)} \right] = \frac{\gamma - 1}{\gamma} \left[ \frac{\partial p^{(0)}}{\partial t} + (\mathbf{u} \cdot \nabla) p^{(0)} \right] + \frac{1}{\text{Pe}} \nabla \cdot (k \nabla T^{(0)}), \qquad (12) $$

where $\tau = \mu[\nabla\mathbf{u} + (\nabla\mathbf{u})^t - \frac{2}{3}(\nabla\cdot\mathbf{u})\mathbf{I}]$, Re is the reference Reynolds number $Re = \tilde{\rho}_{\text{ref}}\tilde{U}_{\text{ref}}\tilde{L}_{\text{stack}}/\tilde{\mu}_{\text{ref}}$, Pe is the reference Péclet number $Pe = \tilde{U}_{\text{ref}}\tilde{L}_{\text{stack}}/\tilde{\alpha}_{\text{ref}}$, and $\mu$, $k$, and $\alpha$ are the dynamic viscosity, thermal conductivity, thermal diffusivity, respectively, scaled by their value at the reference state $\tilde{\mu}_{\text{ref}}$, $\tilde{k}_{\text{ref}}$, and $\tilde{\alpha}_{\text{ref}}$. Gravity has been neglected. $k$ and $\mu$ are assumed to be independent of temperature.

From momentum at leading order, $p^{(0)}$ is a function of time only. However, this value has to equal the pressure $p^{(0)} = 1$ in the resonator. At order $M$, likewise, $p^{(1)}$ is a function of time only, but given that in the resonator $p^{(1)}$ fluctuates, it remains a function of time. Taking the leading order equation of state into account, $\rho^{(0)} = 1/T^{(0)}$, the problem becomes

$$ \frac{\partial \rho^{(0)}}{\partial t} + \nabla \cdot (\rho^{(0)} \mathbf{u}) = 0, \qquad (13) $$
$$ \frac{\partial(\rho^{(0)}\mathbf{u})}{\partial t} + \nabla \cdot (\rho^{(0)}\mathbf{u} \otimes \mathbf{u}) = -\nabla p^{(2)} + \frac{H}{\mathrm{Re}_{\tilde{H}}} \nabla \cdot \tau, \quad (14) $$

$$ \rho^{(0)} \left[ \frac{\partial T^{(0)}}{\partial t} + (\mathbf{u} \cdot \nabla) T^{(0)} \right] = \frac{H}{\mathrm{Pe}_{\tilde{H}}} \nabla^2 T^{(0)}, \quad (15) $$

$$ \rho^{(0)} T^{(0)} = 1, \quad (16) $$

introducing $\mathrm{Re}_{\tilde{H}}$, $\mathrm{Pe}_{\tilde{H}}$, based on the distance $\tilde{H}$ separating two successive plates centerlines, and $H$ being the corresponding dimensionless distance, $H = \tilde{H}/\tilde{L}_{\mathrm{stack}}$.

In the solid plates, the heat conduction equation written in dimensionless form is

$$ \frac{\partial T}{\partial t} = \frac{1}{\mathrm{Pe}_s} \nabla^2 T, \quad (17) $$

where $\mathrm{Pe}_s = \mathrm{Pe}_{\tilde{H}} \alpha_{\mathrm{ref}} / \alpha_s$ is the Péclet number in the wall.

Boundary conditions on symmetry lines (due to vertical periodicity of the stack) include zero normal velocities within the fluid domain and zero normal temperature gradients on all symmetry lines. Fixed temperatures are imposed at the heat-exchanger walls, while continuity of temperature and normal heat flux are imposed at the interface between fluid and stack. Both velocity components vanish on all solid walls. Boundary conditions at the left and right ends are dealt with below and so are the initial conditions.

## F. Matching

In the resonator scaling, the heat-exchanger section shrinks to zero thickness as $M = \tilde{L}_{\mathrm{stack}}/\tilde{L}_{\mathrm{res}} \to 0$. In the heat-exchanger scaling, the resonator length becomes infinity. In other words, in the limit of $M \to 0$, the location $\hat{x} = 0^-$, i.e., the right end of the left resonator section corresponds to $x \to -\infty$ in the heat-exchanger section. Likewise, $\hat{x} = 0^+$, the left end of the right resonator section corresponds to $x \to +\infty$ in the scaling used in the heat-exchanger section. Thus, velocity, pressure, and temperature values must approach the same values at these locations in both solutions.

First, according to the formulation in the heat-exchanger section, pressure is spatially uniform to order $M$; thus, for acoustic pressures,

$$ p^{(1)}(0^{-}) = p^{(1)}(0^{+}) = p^{(1)}(0). \quad (18) $$

Leading order temperatures are constant and uniform in each resonator side. Thus, leading order temperatures in the heat exchangers must approach these fixed values for both $x \to -\infty$ and $x \to +\infty$. The acoustic solution in the resonators provides one equation on each side relating the acoustic pressure $p^{(1)}(0)$, respectively, to longitudinal velocity for $x \to -\infty$, Eq. (7), and $x \to +\infty$, Eq. (8). Transverse velocities vanish in the resonators, hence also for $x \to \pm\infty$.

Finally, integrating the energy equation, Eq. (15), over the entire heat-exchanger section and making use of continuity, a global energy balance is obtained, relating velocities at the two ends to heat transfer (with the end boundaries set to

be adiabatic $\partial T / \partial n = 0$). Taking into account that $u(0^-) = \lim_{x \to -\infty} u$, $u(0^+) = \lim_{x \to +\infty} u$,

$$ [u(0^{-}) - u(0^{+})] + \frac{1}{\mathrm{Pe}_{\tilde{H}}} \int \nabla T^{(0)} \cdot \mathbf{n} \, ds = 0. \quad (19) $$

Finally, solving Eqs. (7), (8), and (19) determines $u(0^-)$, $u(0^+)$, and $p^{(1)}(0)$, providing boundary conditions for both resonators and for the low Mach number problem in the heat-exchanger section. These results show that the heat-exchanger section is transparent to acoustic pressure but appears in the acoustics as a point source of volume.

## G. Initial conditions

The onset problem is a stability problem, in which the issue is the growth of initial perturbations. When using numerical simulation, it is often the case that numerical errors provide an initial perturbation which triggers the instability. However, here, for initial conditions in which the fluid is initially at rest, the initial field is uniform. Since our code is, appropriately, symmetrically implemented, subtracting, however, approximate but identical binary representations at neighboring spatial locations yields differences with a binary representation that is exactly zero; likewise, the value of the initial velocities is zero, which likewise has an exact binary representation. As a result, numerical errors do not provide an initial perturbation. For initial conditions in which the fluid is initially at rest, it stays at rest, as long as the temperature in the heat exchangers and the fluid corresponds to a steady solution. However, if the initial temperature field is not uniform lengthwise, round-off errors result in an initial numerical value of the integral in Eq. (19) that is not exactly zero.

The current model being nonlinear, different initial conditions might trigger different results, and it might be worthwhile testing different initial models. Since the acoustic solution in the resonators, Eqs. (7) and (8), relates the current solution to solutions at earlier times, defining initial conditions is equivalent to defining the previous acoustic history. Likewise, specific temperature distributions can be tested. These issues are discussed in detail below.

## III. NUMERICAL SOLUTION

A finite volume solver suitable for low Mach number flows with large temperature gradients and large pressure fluctuations²⁵ was adapted to deal with the present configuration. The solver has been validated on a benchmark problem of non-Boussinesq natural convection. Time discretization combines an implicit treatment of the viscous and diffusive terms with an explicit discretization for the convective terms. The algorithm implements a second-order spatial and time discretization on a uniform mesh.

The simulation domain is extended beyond the heat-exchanger section by a length equal to twice the stack length. In the matching equations, values on the left and right boundaries are used, where the flow is forced to be one-dimensional. An additional phase variable is used, specifying whether a given cell is to be treated as solid or fluid. The equations are then solved through the whole domain, using appropriate tests
of the phase variable, ensuring continuity of both temperature and heat flux through the plate-fluid interfaces.

Spatial discretization is based upon a staggered mesh. Velocities are defined on cell faces, while density, temperature, and phase are defined on cell centers. In the continuity equation, the divergence of the velocity field is adjusted at each time step using a fractional time step method derived from the projection method for incompressible flows.

The sequence of steps for advancing in time from time $t^n$ to time $t^{n+1}$ is as follows:

(1) The energy equation, Eq. (15), is solved for the new temperature field $T^{n+1}$, using density at time $t^n$, by either a nonincremental alternating direction implicit (ADI) algorithm or the generalized minimal residual method (GMRES).²⁶

(2) The linear system of Eqs. (7), (8), and (19) is solved, determining $u(0^{-})^{n+1}$, $u(0^{+})^{n+1}$, and $p^{(1)}(0)^{n+1}$, based upon values at previous times. The velocity values provide boundary conditions in the calculation of the velocity field. These values and the new value of the acoustic pressure are stored since they will be needed at future time steps.

(3) The new density is determined from the equation of state, Eq. (16).

(4) The new divergence, at time $t^{n+1}$, is obtained by combining the energy equation and the continuity equation, yielding $\nabla \cdot \mathbf{u}^{n+1} = (H/\mathrm{Pe}_{\tilde{H}})\Delta T^{n+1}$.

(5) An intermediate velocity vector is calculated which does not satisfy continuity, from momentum, Eq. (14), using ADI or GMRES, with dynamic pressure at $t_n$.

(6) The dynamic pressure correction is determined using the multigrid method.²⁷ Numerical parameters in the multigrid algorithm are adjusted for optimal convergence.

(7) Respective arrays are updated with the new velocity field $\mathbf{u}^{n+1}$ and dynamic pressure.

# IV. RESULTS

## A. Device and configurations simulated

The various cases that were simulated correspond to a single engine geometry, based upon an existing experimental set-up.²⁸ The overall length is $\tilde{L}_{\mathrm{res}} = 7.57$ m, with a stack $\tilde{L}_{\mathrm{stack}} = 150$ mm long, made up of stainless steel plates 0.20 mm thick separated by a spacing $\tilde{h} = 0.77$ mm wide. The stack porosity is $\tilde{h}/\tilde{H} = 0.794$, where $\tilde{H}$ is the distance separating two successive plate centerlines, and $\tilde{h}$ is the gap or flow passage between two successive stack plates. The distance between the heat exchangers and the ends of the stack plates equals $\tilde{h}$, and the heat-exchanger length is taken to be $\tilde{L}_{\mathrm{stack}}/20$. The length of the left resonator section is $\tilde{L}_L = 676$ mm. The device is filled with helium at 1 MPa and a temperature of 293 K on the cold side. The heat exchangers are taken to be made of the same material as the stack plates. The resonant frequency of the device is 66.6 Hz when filled with helium at 293 K in both cold and hot sections, no thermoacoustic heat pumping, and no load. These data result in a ratio of 0.38 between thermal boundary layer and spacing between two successive

stack plates, which is appropriate for thermoacoustic energy transfer.²

Using the scalings previously defined, the reference Mach number is equal to $M = 0.02$, so that $\tilde{L}_{\mathrm{res}} = 50 \tilde{L}_{\mathrm{stack}}$ [Eq. (1)]. The Péclet and Reynolds numbers, based on $\tilde{H}$ are $\mathrm{Pe}_{\tilde{H}} = 1093$ and $\mathrm{Re}_{\tilde{H}} = 1584$. The solid Péclet number for the stack plates and the heat exchangers is $\mathrm{Pe}_s = 4.485 \times \mathrm{Pe}_{\tilde{H}} = 31.7$.

A combination of values of the dimensionless temperature in the hot heat exchanger ($T_{\mathrm{hot}} = 1.2, 1.3$, and 2) and of the load ($f = \infty, 90$) were simulated. Two initial conditions were tested, with initial velocities and acoustic pressures in a resonator excited by either: (1) A random noise or (2) a standing wave at the cold resonator frequency, with no load and no thermoacoustics. In the heat-exchanger section, both initial dynamic pressure and velocity were set to zero. Two sets of initial temperature conditions in the stack and the fluid in the gap were also considered: A longitudinally linear distribution between heat exchangers, in both solid and fluid, or a distribution corresponding to steady-state heat conduction, with the fluid at rest. The latter was calculated numerically, taking into account the geometry of the heat-exchanger section, including the gas layer between heat exchangers and stack plates, and the thermal conductivities in the wall and the gas.

## B. Convergence analysis

A convergence analysis was performed for $T_{\mathrm{hot}} = 1.2$, with a closed cold end, starting from random noise, and assuming balanced initial temperatures due to conduction in the heat-exchanger section, with spatial resolutions of 1024 × 32, 2048 × 32, and 2048 × 64 grid points. The coarse mesh allows for three grid points inside each half-plate width and three grid points in the fluid gap between heat exchanger and plate end. The time step was chosen as large as possible, maximizing the Courant number while ensuring stability. The coarsest resolution used was about 600 time steps per dimensionless unit of time, which, since initially velocities are small, is appropriate. As velocities increased, the time step was decreased by a factor of 10, to a value of 1/6000. The CPU (Central Processing Unit) time was about 1 µs per time step and per mesh point on a NEC SX8. The numerical parameters for the multigrid algorithm dealing with pressure correction and the divergence of the velocity field are adjusted so that for the coarsest grid, an average of only three or four iterations is needed. On the coarsest grid, the average CPU time per run is 2 h when only simulating the early phase of the process, showing amplification of the initial perturbation. It took about 50 h until saturation was observed, when the wave amplitude no longer increases. These are good performance figures when compared with current literature reports.²³

While qualitatively the results were observed to be independent of the grid size, quantitative differences of about 10% in the growth rates were noted between coarse and fine mesh, but only of 1% between the two finest meshes. For each spatial grid, solutions were obtained using three different time steps, from 1/600th to 1/6000th of the unit time. In order to
assess any effect of the truncation error inherent to the ADI solver, while GMRES yields an exact solution of the discretized equations, results using both solvers were compared. For all the time steps that were used, they were nearly identical. Since GMRES is more costly, all results below were obtained using the ADI solver, and for the coarsest mesh.

# V. DISCUSSION

## A. Initial amplification

Figures 2 and 3 show results obtained for an initially linear temperature profile between the heat exchangers, with $T_{\text{hot}} = 1.3$, in a tube closed at both ends ($f = \infty$). The initial excitation of the resonator was a standing wave of small amplitude at the resonator fundamental frequency in cold fluid. The results shown correspond to the early part of the process, while wave amplification is taking place, for $0 \le t \le 410$. Figures 2(a) and 2(b) show respectively, the evolution of acoustic pressure $Mp^{(1)}(0)$ and velocity $u(0^{-})$ (on the left side of the heat exchangers). The results show the thermoacoustic instability resulting in an approximately exponential amplitude increases. Indeed when acoustic pressure is plotted on a logarithmic scale, Fig. 3(a), the envelope of the signal appears to be linear. However, in Fig. 3(b), the positive and negative parts of the acoustic pressure and velocity signals are noticeably asymmetric, as predicted from nonlinear theoretical mod-

<figure><img src="image_73.png" /><figcaption>FIG. 2. Acoustic pressure $Mp^{(1)}(0)$ (a) and left velocity $u(0^{-})$ (b) history, $f = \infty$, $T_{\text{hot}} = 1.3$, with linear initial temperature distribution and initial standing wave.</figcaption></figure>

<figure><img src="image_74.png" /><figcaption>FIG. 3. Acoustic pressure $Mp^{(1)}(0)$ history, in semi-log scale (a) and detail (b) of acoustic pressure $p^{(1)}(0)$, left velocity $u(0^{-})$, right velocity $u(0^{+})$ history, $f = \infty$, $T_{\text{hot}} = 1.3$, with linear initial temperature distribution and standing wave as initial condition.</figcaption></figure>

els and observed in experiments.<sup>11,14,28</sup> Predictably, the dominant frequency is very close to 0.5, corresponding to the resonance in a tube at a uniformly cold temperature, closed at both ends, with no thermoacoustic power generation and no load. For the relatively modest temperature ratio of 1.3, the Lees/Rott dispersion relation,<sup>1,29</sup> which accounts for a temperature discontinuity at the location of the heat-exchangers section, yields a frequency that only differs by 0.053% from the resonant frequency in a tube entirely filled with cold fluid, which was used as the initial excitation. The effect of the thermoacoustic load resulting in a nonzero growth rate, which also shifts frequencies, is small too, as readily verified in a linear model accounting for a source of volume at the location of the temperature jump. The fundamental is the most energetic mode, which is consistent since the initial excitation starts by feeding energy in the fundamental exclusively. However, closer analysis of the signal shows that the first harmonic, oscillating twice as fast, is also unstable, as reported in both theoretical studies and experiments.<sup>4,6,28</sup>

The difference between the left and the right velocity amplitudes is shown in Fig. 3, which presents a close-up view of acoustic pressure and axial velocities at the end of the computation. The phase shift between pressure and left velocity is close to 90°. Since the acoustic energy flux is proportional to the integral of the product pressure times velocity, which on the left, becomes zero in the stationary regime, a phase shift of 90° corresponds to a standing wave. Here,
<figure><img src="image_75.png" /><figcaption>FIG. 4. Acoustic pressure $Mp^{(1)}(0)$ history, $f = \infty$, $T_{\text{hot}} = 1.3$, with linear initial temperature distribution and random noise as initial perturbation.</figcaption></figure>

however, the amplitude is still increasing, which requires a nonzero energy flux, smaller on the left than on the right since the left side represents only about 9% of the resonator length.

In Fig. 3(b), corresponding to the end of the run in Fig. 2, the drive ratio at the heat-exchanger section location [defined as $\max|Mp^{(1)}(0)|$] is 3.5% and dimensionless velocities are approximately 0.15. Pursuing the simulation with such high velocities requires the time step to be significantly reduced. With no load (closed end, infinite $f$), saturation will only occur on a time scale of $1/M$. Indeed, the irreversibility due to damping in the heat-exchanger section is smaller than the acoustic energy fluxes by order $M$, while the resonator model is isentropic.

Figure 4 shows the growth of pressure for the same configuration and parameters ($T_{\text{hot}} = 1.3$, $f = \infty$, and linear initial temperature distribution), except for the initial excitation, which is a random noise. Since, in this case, there is no significant energy initially deposited in any mode in particular, the effect of the strong growth rate of the first harmonic becomes clear. Both the fundamental and the first harmonic appear to be unstable.

When increasing the hot temperature to $T_{\text{hot}} = 2$ in Fig. 5, starting from a standing wave, the first harmonic (second

<figure><img src="image_76.png" /><figcaption>FIG. 5. Acoustic pressure $Mp^{(1)}(0)$ history, $f = \infty$, $T_{\text{hot}} = 2$, with linear initial temperature distribution and standing wave as initial perturbation.</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>FIG. 6. Acoustic pressure $Mp^{(1)}(0)$ history, $t > 2500$, $f = 90$, $T_{\text{hot}} = 1.2$, with linear initial temperature distribution (a) or steady-state heat conduction initial temperature (b) and initial random noise.</figcaption></figure>

mode) is clearly unstable and can be seen already very early on. Similar observations are made when starting with noise.

## B. Periodic regime

In order to show the evolution of the wave toward a stationary regime, a purely resistive load was introduced at the right end. For $f$ ranging between 80 and infinity, increasing

<figure><img src="image_78.png" /><figcaption>FIG. 7. Time history (detail) of acoustic pressure $p^{(1)}(0)$, left velocity $u(0^{-})$, right velocity $u(0^{+})$, difference $u(0^{+}) - u(0^{-})$ (scale on right vertical axis), $t = \tilde{t}/(2\tilde{r}_{\text{ref}})$, $f = 90$, $T_{\text{hot}} = 1.2$ with linear initial temperature distribution and initial random noise.</figcaption></figure>
the impedance $f$ leads to an increase in the growth rate, and saturation comes sooner. Results below were obtained for $f = 90$, which leads to saturation with a reasonable pressure amplitude level and in a reasonable time. Two separate calculations were done for $f = 90$, both for a temperature ratio of 1.2. The first one started from an initially linear temperature distribution between heat exchangers, with results shown in Figs. 6(a) and 7. In the second simulation, temperature initial conditions considered a steady-state heat conduction field. Results are shown in Fig. 6(b). In the former case, saturation occurs after $t \approx 4500$ (which corresponds to about 35 s in real time). Figure 7 shows acoustic pressure, axial velocities, and the difference between axial velocities for a narrow time window at the end of the simulation. The phase shift between acoustic pressure and left velocity is very close to $90^\circ$ and the period is approximately 0.998, which is again very close to the first harmonic, and even more so if taking the presence of a temperature jump into account. Indeed, its period is 0.9986 in a tube closed at both ends (i.e., with no load) according to the Lees/Rott dispersion relation.¹,²⁹ At the stationary regime, acoustics on the left side is a monochromatic standing wave since there is no longer any energy transfer. On the right side, however, the power developed by the heat-exchanger section acting as a source of volume is carried to the load by the acoustics, which consequently do no longer behave as a pure standing wave (Fig. 7). The final drive ratio ($\max|Mp^{(1)}(0)|$) is about 3.2% and the velocity amplitudes are 0.59. The signal in Fig. 7 is now much more symmetric than for infinite $f$. An overshoot is clearly noticeable in the pressure amplitudes, before the stationary oscillatory regime is reached. The same observation has been reported in previous theoretical, numerical,¹⁴,²³ and experimental studies. This overshoot may be due to interaction between the evolution of the stack temperature and the thermoacoustic wave amplification, which are taking place simultaneously, each of which is characterized by its own time scale. The initially linear temperature gradient evolves mostly by conduction in the solid plates and in the fluid, toward a stationary solution which, for heavy stack plates, will be nearly steady. The time scale characterizing conduction in the stack is of the order of 1–2 h, i.e., approximately $3.5 \times 10^6$ dimensionless time units (based on the stack length and the thermal diffusivity of the stack plates $\tau \approx \overline{L}_{stack}^2/\alpha_s$). Meanwhile, the thermoacoustic engine provides for an increase in amplitude and delivers power to the load, until reaching an equilibrium. Here equilibrium, hence a nearly periodic regime, appears to be reached while the stack temperature is still adjusting due to its high thermal mass. In contrast, in the case of steady-state heat conduction field, the initial temperature profile has already been adjusted for conduction (but not for the effect of oscillating flow), and it is much closer to its final profile.

In the case of steady-state heat conduction, as a result of relatively poor conduction between heat exchangers and stack, the actual temperature gradient in the stack is much lower already when the device starts, the engine is clearly weaker, and wave amplification is slower, as observed by comparing times in Fig. 6(b) with those in Fig. 6(a). Figure 8 shows the initial temperature field, zooming between heater and stack plate, due to conduction but no motion, which is

<figure><img src="image_79.png" /><figcaption>FIG. 8. Temperature field detail between the hot heat exchanger and the stack plates due to conduction (no motion), $T_{hot} = 1.2$, same view as Fig. 1 (bottom).</figcaption></figure>

clearly two-dimensional. The initial temperature gradient in the stack plates is significantly smaller than the linear profile used as initial conditions in the previous case and much closer to the final profile, as shown in Fig. 9. This is due primarily to the gap between the heat exchangers and the stack ends. The gas layer in this gap is much less conductive than the solid plates and, therefore, it absorbs much of the temperature difference, resulting in a smaller gradient in the stack. The same drive ratio is reached eventually but, as Fig. 6(b) shows, there is almost no overshoot.

The unavoidable presence of a fluid layer between heat exchangers and plates has a definite impact on the growth rate. Thermal edge effects and gaps between the stack ends and the heat exchangers were further analyzed through linear and nonlinear studies in order to derive improved heat exchange laws and to optimize the systems.³⁰⁻³⁴ The current model will allow for optimization of its width and exploration of the optimum values of other design parameters, such as stack and heat-exchanger location.

## VI. CONCLUSION

A globally consistent low Mach number model of the thermoacoustic engine has been developed, based upon the assumption that the stack and heat-exchangers length is

<figure><img src="image_80.png" /><figcaption>FIG. 9. Temperature distribution along the stack symmetry plane and in the gas from left to right heat exchanger, showing both initial conditions, linear distribution and heat conduction solution, and the mean stationary periodic solution for $f = 90$, $T_{hot} = 1.2$, and initial random noise.</figcaption></figure>
comparable with the fluid particle displacement, resulting in a heat-exchanger section that is transparent to acoustic pressure but acts as a source of volume.

Transient processes whereby amplitude grows and eventually reaches a steady value are predicted, including frequencies and power fluxes. The key components of the thermoacoustic device are included. Yet by virtue of not needing a numerical solution in the resonator, and because the flow in the stack and heat exchangers is reduced to a low Mach number flow, implementation is orders of magnitude faster than compressible models of the entire device. Indeed, the current model entails solving numerically for only a small fraction of the overall volume, and the time step is not restricted by an acoustic Courant-Friedrichs-Lewy condition. The current approach is a good example whereby the singular perturbation technique can be used successfully, allowing the computational effort to be concentrated where it truly matters. Even so, large computations are still required, since the direct simulation of a transient period flow in the heat-exchanger section is still required.

The model will allow for exploring and understanding crucial trade-offs on stack position, stack wall thickness and thermal mass, heat-exchanger size and location, and load impedance including the trade-off between optimizing for device performance and optimizing for its ability to start, which may lead to conflicting requirements.

The resonator and the load models that were used are admittedly much simplified. However, the current approach can readily be extended to more realistic resonators as long as they can be described by linear acoustics, and to more complex load models, including multidimensional resonators with losses, complex losses, and loops, or in which the long acoustic resonator is replaced by a free piston, hence allowing for arbitrarily large pressure amplitudes, and a much higher power density. In particular, losses due to heat and momentum diffusion can readily be included, since these are linear in nature. Wave steepening cannot be dealt with; however, until shocks appear, wave steepening remains a fully reversible phenomenon, which as such does not directly contribute to losses.

## ACKNOWLEDGMENTS

Support of the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged (L.B.). Computations were performed on the NEC-SX8 of the IDRIS-CNRS computing center, under Project No. i2009020599. The authors thank Professor Maurice-Xavier François for his interest and numerous fruitful exchanges on this work and for initiating the collaboration that made this study possible. Finally, support of the Algerian and French governments through a doctoral fellowship 2007–2008, “Programme Boursier Enseignants du Ministère Algérien de l’Enseignement Supérieur et de la Recherche Scientifique” is gratefully acknowledged (O.H.).

<footnote n="1">¹N. Rott, “Damped and thermally driven acoustic oscillations in wide and narrow tubes,” Z. Angew. Math. Phys. **20**, 230–243 (1969).</footnote>

<footnote n="2">²G. W. Swift, “Thermoacoustic engines,” J. Acoust. Soc. Am. **84**, 1145–1180 (1988).</footnote>

<footnote n="3">³W. P. Arnott, H. E. Bass, and R. Raspet, “General formulation of thermoacoustics for stacks having arbitrarily shaped pore cross-sections,” J. Acoust. Soc. Am. **90**, 3228–3237 (1991).</footnote>

<footnote n="4">⁴A. A. Atchley, H. E. Bass, T. J. Hofler, and H. T. Lin, “Study of a thermoacoustic prime mover below onset of self-oscillation,” J. Acoust. Soc. Am. **91**, 734–743 (1992).</footnote>

<footnote n="5">⁵A. A. Atchley, “Standing wave analysis of a thermoacoustic prime mover below the onset of self-oscillation,” J. Acoust. Soc. Am. **92**, 2907–2914 (1992).</footnote>

<footnote n="6">⁶A. A. Atchley and F. M. Kuo, “Stability curves for a thermoacoustic prime mover,” J. Acoust. Soc. Am. **95**, 1401–1404 (1994).</footnote>

<footnote n="7">⁷W. P. Arnott, J. R. Belcher, R. Raspet, and H. E. Bass, “Stability analysis of a helium-filled thermoacoustic engine,” J. Acoust. Soc. Am. **96**, 370–375 (1994).</footnote>

<footnote n="8">⁸B. Ward, J. Clark, and G. W. Swift, *Design Environment for Low-Amplitude Thermoacoustic Energy Conversion (DeltaEC) User's Guide* (version 6.2) (Los Alamos National Laboratory, 2008), http://www.lanl.gov/thermoacoustics/UsersGuide.pdf.</footnote>

<footnote n="9">⁹R. Raspet, J. Brewster, and H. E. Bass, “A new approximation method for thermoacoustic calculations,” J. Acoust. Soc. Am. **103**, 2395–2402 (1998).</footnote>

<footnote n="10">¹⁰A. T. A. M. de Waele, “Basic treatment of onset conditions and transient effects in thermoacoustic Stirling engines,” J. Sound Vib. **325**, 974–988 (2009).</footnote>

<footnote n="11">¹¹S. Karpov and A. Prosperetti, “Nonlinear saturation of the thermoacoustic instability,” J. Acoust. Soc. Am. **107**, 3130–3147 (2000).</footnote>

<footnote n="12">¹²S. Karpov and A. Prosperetti, “A nonlinear model of thermoacoustic devices,” J. Acoust. Soc. Am. **112**, 1431–1444 (2002).</footnote>

<footnote n="13">¹³N. Cao, J. R. Olson, G. W. Swift, and S. Chen, “Energy flux density in a thermoacoustic couple,” J. Acoust. Soc. Am. **99**, 3456–3464 (1996).</footnote>

<footnote n="14">¹⁴M. F. Hamilton, Y. A. Ilinski, and E. Zabolotskaya, “Nonlinear two-dimensional model for thermoacoustic engines,” J. Acoust. Soc. Am. **111**, 2076–2086 (2002).</footnote>

<footnote n="15">¹⁵L. Bauwens, C. R. L. Bauwens, and I. Wierzba, “Oscillating flames: Multiple-scale analysis,” Proc. R. Soc. London, Ser. A **465**, 2089–2110 (2009).</footnote>

<footnote n="16">¹⁶S. Paolucci, “On the filtering of sound from the Navier-Stokes equations,” Report No. SAND82-8257, Sandia National Laboratories (1982).</footnote>

<footnote n="17">¹⁷L. Bauwens, “Near-isothermal regenerator: A perturbation analysis,” AIAA J. Thermophys. Heat Transfer **9**, 749–756 (1995).</footnote>

<footnote n="18">¹⁸L. Bauwens, “Oscillating flow of a heat-conducting fluid in a narrow tube,” J. Fluid Mech. **324**, 135–161 (1996).</footnote>

<footnote n="19">¹⁹A. Worlikar and O. Knio, “Numerical study of oscillatory flow and heat transfer in a loaded thermoacoustic stack,” Numer. Heat Transfer, Part A **35**, 49–65 (1999).</footnote>

<footnote n="20">²⁰P. Duthil, C. Weisman, E. Bretagne, and M.-X. François, “Experimental and numerical investigation of heat transfer and flow within a thermoacoustic cell,” Int. J. Transp. Phenom. **6**, 265–272 (2004).</footnote>

<footnote n="21">²¹D. Marx and P. Blanc-Benon, “Computation of the temperature distortion in the stack of a standing-wave thermoacoustic refrigerator,” J. Acoust. Soc. Am. **118**, 2993–2999 (2005).</footnote>

<footnote n="22">²²J. A. Lycklama à Nijeholt, M. E. H. Tijani, and S. Spoelstra, “Simulation of a traveling-wave thermoacoustic engine using computational fluid dynamics,” J. Acoust. Soc. Am. **118**, 2265–2270 (2005).</footnote>

<footnote n="23">²³G. Y. Yu, E. C. Luo, W. Dai, and J. Y. Hu, “Study of nonlinear processes of a large experimental thermoacoustic-Stirling heat engine by using computational fluid dynamics,” J. Appl. Phys. **102**, 074901–074907 (2007).</footnote>

<footnote n="24">²⁴D. Marx and P. Blanc-Benon, “Numerical simulation of stack-heat exchangers coupling in a thermoacoustic refrigerator,” AIAA J. **42**, 1338–1347 (2004).</footnote>

<footnote n="25">²⁵P. Le Quéré, C. Weisman, H. Paillère, J. Vierendeels, E. Dick, R. Becker, M. Braack, and J. Locke, “Modelling of natural convection flows with large temperature differences: A benchmark problem for low Mach number solvers. Part 1. Reference Solutions,” M2AN, Math. Modell. Numer. Anal. **39**, 609–616 (2005).</footnote>

<footnote n="26">²⁶Y. Saad and M. Schultz, “GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems,” SIAM (Soc. Ind. Appl. Math.) J. Sci. Stat. Comput. **7**, 856–869 (1986).</footnote>

<footnote n="27">²⁷W. Briggs, *A Multigrid Tutorial*, 2nd ed. (SIAM Publications, Philadelphia, 1987).</footnote>

<footnote n="28">²⁸P. Debesse, D. Baltean-Carlès, and F. Lusseyran, “Experimental analysis of nonlinear phenomena in a thermoacoustic system, nonlinear acoustics fundamentals and applications,” in *Proceedings of the 18th International Symposium on Nonlinear Acoustics (ISNA)*, Stockholm, Sweden (July 7–10, 2008), Vol. 1022, pp. 355–358.</footnote>