High-fidelity simulation of an ultrasonic standing-wave thermoacoustic engine with bulk viscosity effects

TL;DR

This study uses high-fidelity simulations to analyze an ultrasonic thermoacoustic engine, incorporating bulk viscosity effects to better understand energy conversion and optimize small-scale engine performance.

Contribution

It introduces a boundary-layer-resolved simulation framework that includes bulk viscosity effects from molecular relaxation, providing new insights into ultrasonic thermoacoustic engine behavior.

Findings

Bulk viscosity varies with pressure, temperature, humidity, and frequency.

Limit cycle amplitudes are maximized at 5% humidity.

Linear models overpredict growth rates but accurately estimate limit cycles.

Abstract

We have carried out boundary-layer-resolved, unstructured fully-compressible Navier--Stokes simulations of an ultrasonic standing-wave thermoacoustic engine (TAE) model. The model is constructed as a quarter-wavelength engine, approximately 4 mm by 4 mm in size and operating at 25 kHz, and comprises a thermoacoustic stack and a coin-shaped cavity, a design inspired by Flitcroft and Symko (2013). Thermal and viscous boundary layers (order of 10 $\mathrm{\mu}$m) are resolved. Vibrational and rotational molecular relaxation are modeled with an effective bulk viscosity coefficient modifying the viscous stress tensor. The effective bulk viscosity coefficient is estimated from the difference between theoretical and semi-empirical attenuation curves. Contributions to the effective bulk viscosity coefficient can be identified as from vibrational and rotational molecular relaxation. The inclusion of the coefficient captures acoustic absorption from infrasonic ($\sim$10 Hz) to ultrasonic ($\sim$100 kHz) frequencies. The value of bulk viscosity depends on pressure, temperature, and frequency, as well as the relative humidity of the working fluid. Simulations of the TAE are carried out to the limit cycle, with growth rates and limit-cycle amplitudes varying non-monotonically with the magnitude of bulk viscosity, reaching a maximum for a relative humidity level of 5%. A corresponding linear model with minor losses was developed; the linear model overpredicts transient growth rate but gives an accurate estimate of limit cycle behavior. An improved understanding of thermoacoustic energy conversion in the ultrasonic regime based on a high-fidelity computational framework will help to further improve the power density advantages of small-scale thermoacoustic engines.

Equations

$$\frac{\partial }{\partial t}\left(\rho \right)$$

$$\frac{\partial }{\partial t}\left(\rho u_i\right)$$

$$\frac{\partial }{\partial t}\left(\rho E\right)$$

$${\tau }_{ij}$$

$$q_j$$

$${\mu }_B\equiv \lambda +\frac{2}{3}\mu ,$$

$${\alpha }_{cl}$$

$${\alpha }_a$$

$${\alpha }_k$$

$${\left({\alpha }_k\lambda \right)}_m$$

$$c_{v,k}$$

$$\lambda$$

$${\mu }_{B,vib}$$

$$=\sum _k\left[\frac{p_0}{2\pi }{\left(\gamma -1\right)}^2\left(\frac{n_k}{n}{\left(\frac{T_k^{\ast }}{T_k}\right)}^2\exp (-T_k^{\ast }/T_k)\right)\right]\frac{f_k}{f_k^2+f^2},$$

$$\frac{1}{2\pi {\tau }_O}\equiv f_O$$

$$\frac{1}{2\pi {\tau }_N}\equiv f_N$$

$${\log }_{10}\left(p_{\text{sat}}/p_{s0}\right)$$

$$i\sigma \hat{p}$$

$$i\sigma \hat{U}$$

$$f_{\nu }=\frac{2}{i{\eta }_w}\frac{J_1(i{\eta }_w)}{J_0(i{\eta }_w)},f_{\kappa }=\frac{2}{i{\eta }_w\sqrt{\text{mbox}Pr}}\frac{J_1(i{\eta }_w\sqrt{\text{mbox}Pr})}{J_0(i{\eta }_w\sqrt{\text{mbox}Pr})}$$

$$\eta \equiv \sqrt{\frac{i\omega }{{\nu }_0}}r=\sqrt{2i}\frac{r}{{\delta }_{\nu }}$$

$${\delta }_{\nu }=\sqrt{\frac{2{\nu }_0}{\omega }},{\delta }_{\kappa }=\sqrt{\frac{2k}{\omega {\rho }_0{c}_p}}$$

$$i\sigma \hat{p}$$

$$i\sigma \hat{U}$$

$$A_g=\sum _{m=1}^{n_s+1}{A}_g^{(m)},A_g^{(m)}={\int }_{r_{\text{bot}}^{(m)}}^{r_{\text{top}}^{(m)}}2\pi rdr$$

$$\hat{U}=\sum _{m=1}^{n_s+1}{\hat{U}}^{(m)},{\hat{U}}^{(m)}={\int }_{r_{\text{bot}}^{(m)}}^{r_{\text{top}}^{(m)}}2\pi r\hat{u}(r)dr$$

$$\begin{split}f^{(m)}_{\nu}=-\frac{\pi\,\delta_{\nu}^{2}}{A_{g}^{(m)}}\Big{\{}\frac{1}{J_{0}(i\,\eta^{(m)}_{\textrm{top}})}\left[\eta^{(m)}_{\textrm{top}}J_{1}(i\eta^{(m)}_{\textrm{top}})-\eta^{(m)}_{\textrm{bot}}J_{1}(i\eta^{(m)}_{\textrm{bot}})\right]+\\ \frac{1}{H^{(1)}_{0}(i\,\eta^{(m)}_{\textrm{bot}})}\left[\eta^{(m)}_{\textrm{top}}H^{(1)}_{1}(i\eta^{(m)}_{\textrm{top}})-\eta^{(m)}_{\textrm{bot}}H^{(1)}_{1}(i\eta^{(m)}_{\textrm{bot}})\right]\Big{\}}\end{split}$$

$$\begin{split}f^{(m)}_{\kappa}=-\frac{\pi\,\delta_{\kappa}^{2}\,\sqrt{\mbox{{Pr}}}}{A_{g}^{(m)}}\Big{\{}\frac{1}{J_{0}(i\,\eta^{(m)}_{\textrm{top}}\sqrt{\mbox{{Pr}}})}\left[\eta^{(m)}_{\textrm{top}}J_{1}(i\eta^{(m)}_{\textrm{top}}\sqrt{\mbox{{Pr}}})-\eta^{(m)}_{\textrm{bot}}J_{1}(i\eta^{(m)}_{\textrm{bot}}\sqrt{\mbox{{Pr}}})\right]+\\ \frac{1}{H^{(1)}_{0}(i\,\eta^{(m)}_{\textrm{bot}}\sqrt{\mbox{{Pr}}})}\left[\eta^{(m)}_{\textrm{top}}H^{(1)}_{1}(i\eta^{(m)}_{\textrm{top}}\sqrt{\mbox{{Pr}}})-\eta^{(m)}_{\textrm{bot}}H^{(1)}_{1}(i\eta^{(m)}_{\textrm{bot}}\sqrt{\mbox{{Pr}}})\right]\Big{\}}\;,\end{split}$$

$$i\sigma \hat{p}$$

$$i\sigma \hat{U}$$

$$f$$

$$i\sigma {\hat{p}}_J$$

$$\Delta {\hat{p}}_{ml}$$

$${\zeta }_e$$

$${\zeta }_c$$

$$p_{\text{eff}}$$

$$\frac{1}{{\rho }_0{a}_0^2}\frac{\partial p^{\prime }}{\partial t}+\nabla \cdot \overrightarrow{u^{\prime }}$$

$$\frac{\omega }{{\rho }_0{a}_0^2}\left|p^{\prime }\right|+\left|\nabla \cdot \overrightarrow{u^{\prime }}\right|$$

Full text

1 Introduction

1.1 Background

Thermoacoustic engines (TAEs) are devices capable of converting external heat sources into acoustic power, which in turn can be converted to mechanical or electrical power. TAEs do not require moving parts and are thermoacoustically unstable past a critical heat input; given this condition, an initial perturbation is sufficient to generate acoustic power amplification. The acoustic nature of wave energy propagation in TAEs guarantees close-to-isentropic stages in the thermoacoustic energy conversion process, suggesting the possibility for high efficiency external heat engine designs. For example, advanced TAEs have achieved thermal-to-acoustic energy conversion efficiency of 32%, corresponding to 49% of Carnot’s theoretical limit. [2] There are a variety of TAEs in use for energy production and heat pumping, with varying sizes, and heat sources and energy extraction strategies. [3]

In a TAE, understanding the energy conversion process from thermal to acoustic is crucial. The latter is fluid dynamic in nature and understood and predictable at various levels of fidelity, from quasi one-dimensional linear acoustics [4] to fully compressible three-dimensional Navier–Stokes models. [5] However, energy production and dissipation in high-amplitude or high-frequency devices are difficult to model without high-fidelity simulations. Similitude offers some answers, but may break down as certain assumptions are violated. [6] In the following, we demonstrate a computational modeling strategy, building upon a Navier–Stokes solver, which can simulate high-frequency thermoacoustic amplification with high fidelity.

Modern research has been focused on achieving conversion efficiencies comparable to theoretical expectations. Ceperley realized that the thermodynamic cycle induced by engines with traveling wave phasing comprises discrete stages of compression, heating, expansion, and cooling; in standing-wave engines, these stages are partly overlapped and lead to lower energy conversion efficiencies. [7] However, Ceperley was unsuccessful in developing a working traveling-wave TAE; the first practical realization is attributed to Yazaki et al. [8] TAEs can therefore largely be classified into standing-wave and traveling-wave configurations, the latter being often more efficient but more complicated to build. Hybrid configurations are also possible, with the two concepts combined in a cascaded system. [9]

Modern thermoacoustic engines are generally large and commonly have operating frequencies between 60 and 400 Hz. This design space is constrained by design considerations such as construction limitations and electroacoustic transducer efficiency. However, there are significant benefits to engines operating at higher frequencies. Power density for standing-wave engines scales favorably with frequency and pressure amplitude, holding constant the operating temperature range.[10]. Further, acoustic-to-electric conversion via piezoelectric transducers can be more efficient at higher frequencies.[11, 12, 13]

Acoustic simulations in the ultrasonic regime present some unique challenges. Loss mechanisms from thermodynamic non-equilibrium, which can be neglected at lower frequencies, can become dominant[14], making Stokes’s hypothesis invalid as a result. Moreover, acoustic streaming may no longer be assumed to be a second-order flow quantity.[15] Oscillating viscous boundary layers may no longer evolve in the continuum regime, and the onset of slip-flow near the wall may be possible. As a result, wall-heat transfer may have longer timescales, limiting the intensity of the thermoacoustic response.[1]

Previous high-fidelity efforts by Scalo et al.[5] demonstrated a full-scale three-dimensional simulation of a large $\sim$60 Hz TAE, revealing the presence of transitional turbulence and providing support for direct low-order modeling of acoustic nonlinearities such as Gedeon streaming. More recently, Lin et al.[16] carried out high-fidelity modeling of a $\sim$390 Hz mid-size piezoelectric thermoacoustic energy harvester with particular emphasis on time-domain modeling of electromechanical transmittance functions. The goal of the current work is to tackle the modeling of a miniature, realistic ultrasonic $\sim$25 kHz engine by accurately capturing ultrasonic attenuation concurrently with thermoacoustic instability, with the support of both high-fidelity and low-fidelity numerical prediction tools.

1.2 Research Aims

In this paper we present a high-fidelity fully compressible Navier–Stokes simulation of a thermoacoustic engine operating in the ultrasonic frequency regime. A model for bulk viscosity, incorporating rotational and vibrational relaxation effects, has been developed; these effects are not insignificant compared with Stokesian viscous and thermal dissipation[17] and thermoacoustic energy production. The engine design is based on a standing-wave engine construction first presented by Flitcroft and Symko [1] and was chosen due to its simple design. To the authors’ knowledge, this design is the sole experimental example of an ultrasonic TAE in literature. Being able to capture thermoacoustic onset and nonlinear effects are preliminary steps towards the development of computational tools to better predict and optimize energy generation for miniaturized and high-frequency thermoacoustic engines.

In the following, the adopted theoretical TAE model is first introduced, together with the governing equations and computational setup (section 2). A technique to capture bulk viscosity and a setup for absorption verification is presented (section 3.2.2). A linear thermoacoustic model predicting the onset and growth of oscillations in the TAE model is presented (section 4). Finally, results for the thermoacoustic engine model for both the linear and Navier–Stokes models are shown and discussion follows (section 6).

2 Engine Model

2.1 Engine Model Design and Computational Setup

The chosen computational setup (fig. 1) is a simplified axisymmetric model of the miniature standing-wave thermoacoustic engine first presented by Flitcroft and Symko[1]. Our model is a quarter-wavelength resonator with a thermoacoustic stack in a straight circular tube, closed on one end and connected to a coin-shaped cavity on the other. The cavity lowers the critical temperature necessary for onset and also provides for the possibility of pressurizing the engine, and the cavity geometry is estimated from a presentation by Flitcroft and Symko[18]. The stack temperature profile varies from $T_{c}$ of 300 K to $T_{h}$ of 600 K in our simulations.

The referenced literature reports minimal geometrical information, with only the diameter and length of the tube being explicitly provided. As a result, the stack position and the diameter and height/width of the cavity were redesigned such that the engine achieves onset of thermoacoustic instability. The stack is constructed as radially-concentric plates, as in Lin et al.,[16] with porosity $\phi_{s}=0.6$. The number of concentric stack elements ($n_{s}$) was chosen to be 7, including the centered cylindrical rod, resulting in a stack gap width $h_{g}$ of 0.0706 mm. Linear approximations, using Rott’s wave equations, suggest an operating frequency of approximately 21 kHz, while the fully nonlinear Navier–Stokes simulations suggest an operating frequency of approximately 25 kHz. Because the defined cavity volume may be different from the experimental setup, the simulation-derived frequency is not expected to match the reported experimental operating frequency of 21 kHz.

The computational grid, as also shown in fig. 1, is axisymmetric and designed to resolve thermoviscous boundary layers. Rotational extrusion of five layers, each of one degree, along the $x$ axis is used to construct the three-dimensional computational grid. Adiabatic slipwall conditions are used to impose axial symmetry.

The high-fidelity model was run both without bulk viscosity (reference) and varying levels of bulk viscosity, as tuned by relative humidity. For presented results, several cases were run, differing by gas attenuation magnitude. As discussed in detail in section 3.2, the relative humidity for atmospheric air varies the effective bulk viscosity of the fluid significantly.

3 Governing Equations and Bulk Viscosity Model

3.1 Fully compressible Navier–Stokes equations

The conservation equations for mass, momentum, and energy, solved in the fully compressible Navier–Stokes simulations of the presented TAE model are, respectively,

[TABLE]

where $x_{1}$, $x_{2}$, and $x_{3}$ (equivalently, $x$, $y$, and $z$) are axial and cross-sectional coordinates, $u_{i}$ are the velocity components in each of those directions, and $p$, $\rho$, and $E$ are respectively pressure, density, and total energy per unit mass. The gas is assumed to be ideal, with equation of state $p=\rho\,R_{gas}\,T$ and a constant ratio of specific heats, $\gamma$. The gas constant is fixed and calculated as $R_{gas}=p_{\textrm{ref}}\left(T_{\textrm{ref}}\,\rho_{\textrm{ref}}\right)^{-1}$, based on the reference thermodynamic density $\rho_{\textrm{ref}}$, pressure $p_{\textrm{ref}}$, and temperature $T_{\textrm{ref}}$. The viscous and conductive heat fluxes are, respectively,

[TABLE]

where $S_{ij}$ is the strain-rate tensor, given by $S_{ij}=(1/2)\left(\partial u_{j}/\partial x_{i}+\partial u_{i}/\partial x_{j}\right)$; Pr is the Prandtl number; and $\mu$ is the dynamic viscosity, given by $\mu=\mu_{\textrm{ref}}\left(T/T_{\textrm{ref}}\right)^{n}$, where $n$ is the viscosity power-law exponent and $\mu_{\textrm{ref}}$ is the reference viscosity. $\lambda$ is the second viscosity defined by

[TABLE]

where $\mu_{B}$ is the bulk viscosity. One significant advancement in this work is the adoption of a newly developed bulk viscosity model, accounting for both rotational and vibrational molecular relaxation, as outlined in the following section. Simulations have been carried out with the following gas properties: $\gamma=1.4$, $\rho_{\textrm{ref}}=1.2\,\textrm{kg m}^{-3}$, $p_{\textrm{ref}}=101\,325\,\textrm{Pa}$, $T_{\textrm{ref}}=300\,\textrm{K}$, $\mu_{\textrm{ref}}=1.98\times 10^{-5}\,\textrm{kg}\,\textrm{m}^{-1}\textrm{s}^{-1}$, $\mbox{{Pr}}=0.72$, and $n=0.76$, valid for air.[19]

The governing equations are solved using CharLESX, a control-volume-based, finite-volume solver for the fully compressible Navier–Stokes equations on unstructured grids, developed as a joint-effort among researchers at Stanford University. CharLESX employs a three-stage, third-order Runge-Kutta time discretization and a grid-adaptive reconstruction strategy, blending a high-order polynomial interpolation with low-order upwind fluxes.[20] The code is parallelized using the Message Passing Interface (MPI) protocol and highly scalable on a large number of processors.[21]

3.2 Bulk Viscosity Model

3.2.1 Model formulation

In the following, we outline a novel procedure for estimating the bulk viscosity, $\mu_{B}$, from the absorption coefficient. The latter is the measure of wave attenuation over a given traveled distance and has traditionally been of particular interest for atmospheric acoustics. The classical absorption coefficient $\alpha_{cl}$ is

[TABLE]

where $\omega$ is angular frequency, and $\kappa$ is the heat conductivity. However, multispecies interactions and rotational and vibrational relaxation result in deviations from classic absorption characteristics at various frequency regimes. [17] When both bulk viscosity and vibrational relaxation contributions are considered, the absorption coefficient $\alpha_{a}$ is

[TABLE]

where the subscript $k$ indicates the contribution from the $k$-th species, $\tau_{k}$ is the associated relaxation time given by the semi-empirical relationships for relaxation frequencies, as in eq. 12, $\left(\alpha_{k}\lambda\right)_{m}$ denotes the maximum absorption per wavelength for the $k$-th species, $n_{k}/n$ is the mole fraction for the $k$-th species, and $T_{k}^{*}$ is the characteristic molecular vibration temperature for the $k$-th species. Species in air, for example, are that of nitrogen and gas, with corresponding $n_{k}/n$ of 0.21 and 0.78, respectively.

In the present work, the bulk viscosity and absorption contributions from rotational and vibrational relaxation are collapsed into one equation, such that

[TABLE]

where at atmospheric conditions, $\mu_{B,rot}\approx 0.6\mu$, [17] the functional form of $f_{k}$ is dependent only on temperature, and thus the form of effective bulk viscosity $\mu_{B}$ is both frequency and temperature-dependent.

The relaxation frequencies of species in air is given by Bass et al.[22]:

[TABLE]

where $p_{s0}$ is the reference atmospheric pressure, $p_{\textrm{sat}}$ is the saturation vapor pressure, and $T_{01}$ is the triple-point isotherm pressure. The contribution of water vapor to the relaxation frequency of each species in air is determined via an adjustment using the relative humidity $h_{r}$, which defines the absolute humidity $h=h_{r}\frac{p_{\textrm{sat}}}{p_{0}}$.

In the vibrational relaxation term, the relaxation times $\tau_{k}$ are assigned according to relaxation frequencies $f_{k}=\left(2\pi\tau_{k}\right)^{-1}$, which have been semi-empirically determined by Bass et al.[22] Sample curves along various relative humidity levels, using the appropriate effective bulk viscosity developed in the preceding equations, are shown in fig. 3; these curves accurately replicate experimental measurements of absorption in air.

The form of $\mu_{B}$ is interpreted in some literature as a frequency-dependent bulk viscosity, as absorption of acoustic power is the primary technique for experimentally measuring the bulk viscosity of a gas. [23] The dependence of $\mu_{B,vib}$ on temperature and pressure is qualitatively depicted in fig. 2. The current modeling framework and the ideal gas assumption can break down in dense gas situations; for example, in dense gases, even monatomic gases can exhibit bulk viscosity. [24] It is also important to note that an implementation of $\mu_{B}$ does not allow a solver to fully capture dispersion effects. However, because most thermoacoustic engines do not change frequencies significantly in transient and in limit cycle operation, especially if care is taken to minimize thermoacoustic streaming, the predictive ability of this setup is expected to hold in most cases.

3.2.2 Time-domain ultrasonic acoustic absorption verification

The modifications made to the governing equations are tested against semi-empirical and analytical expressions for the absorption of sound in air.

Single-frequency traveling wave one-dimensional simulations with 4096 points per wavelength have been performed in a periodic domain. Frequencies in the range $f=10^{3}-10^{5}$ Hz have been tested with relative humidity levels of 20%. Due to the monochromatic nature of the wave propagation, the frequency input to the effective bulk viscosity is fixed throughout a single simulation. Initial conditions were set as a pure traveling isentropic wave. Numerical experiments yielded $\alpha_{a}$ as annotated in fig. 3, extracted from time-series decaying pressure amplitudes (fig. 4).

4 Linear Thermoacoustic Eigenvalue Model

A system-wide linear model has been developed based on Rott’s theory, to support the Navier–Stokes calculations both in the start-up phase and the limit cycle. The engine is divided into five Eulerian control volumes, as shown in fig. 5: a pre-stack heated duct, the gas-filled volume of the stack, an after-stack duct, an acoustic junction, and the disk-shaped cavity. The governing equations have been linearized about the thermodynamic state $\{\rho_{0},T_{0},P_{0}\}$. The base pressure, $P_{0}$, is assumed to be uniform, and the mean density and temperature vary with the axial coordinate according to $P_{0}=\rho_{0}(x)\,R_{gas}\,T_{0}(x)$. The base speed of sound is calculated as $a_{0}=\sqrt{\gamma R_{gas}T_{0}}$. All fluctuating quantities are assumed to be harmonic. The $e^{+i\sigma\,t}$ convention is adopted where $\sigma=-i\alpha+\omega$, with $\alpha$ and $\omega$ being the growth rate and angular frequency, respectively. The thermoacoustic growth rate $\alpha$, measured as $\textrm{s}^{-1}$, can be related to the absorption constant $\alpha_{a}=\alpha_{cl}+\sum_{k}\alpha_{k}$ by $-\alpha_{a}\,a_{0}$.

4.1 Ducts in the $x$ direction

In the ducts, a constant axial mean temperature is assumed, yielding the linearized equations

[TABLE]

which enforce the conservation of mass and energy (15a) and momentum (15b), respectively. The total cross-sectional area is the area available to the gas, $A=A_{g}$. The complex thermoviscous functions $f_{\nu}$ and $f_{\kappa}$ in (15) are

[TABLE]

where $J_{n}(\cdot)$ are Bessel functions of the first kind and $\eta$ is the dimensionless complex radial coordinate

[TABLE]

where $\nu_{0}=\mu(T_{0})/\rho_{0}$ is the kinematic viscosity based on mean values of density and temperature, and $\eta_{w}$ in (16) is the dimensionless coordinate (17) calculated at the radial location of the isothermal, no-slip wall. The viscous, $\delta_{\nu}$, and thermal, $\delta_{\kappa}$, Stokes thicknesses are

[TABLE]

and are related via the Prandtl number, $\delta_{\nu}=\sqrt{\mbox{{Pr}}}\;\delta_{\kappa}$.

4.2 Thermoacoustic stack in the $x$ direction

In the thermoacoustic stack, assuming that all annular flow passages share the same pressure field and collapsing via area-weighted averaging leads to a set of approximate linearized equations,

[TABLE]

where the total cross-sectional area available to the gas, $A_{g}$, and flow rate, $\hat{U}$, are

[TABLE]

[TABLE]

and an area-weighted equipartitioning of the flow rates, $\hat{U}^{(m)}=A_{g}^{(m)}/A_{g}\;\hat{U}$, has been assumed. The accompanying thermoviscous functions are

[TABLE]

where $J_{n}$ and $H_{n}^{(1)}$ are the Bessel functions of the first kind and Hankel functions of the first kind, respectively. A detailed mathematical derivation is found in Lin et al.[16]

4.3 Cavity in the $r$ direction

In the radial disk cavity, the governing equations for $\hat{p}$ and $\hat{U}$ vary with perimeter and area as a function of the radius:

[TABLE]

The thermoviscous functions used in the radial disk are that of parallel plates,

[TABLE]

where $L$ is the width of the cavity.

4.4 Junction and losses

In the compliance junction, a relationship on the pressure, input, and output volume flow rates is imposed as

[TABLE]

where $\hat{U}_{x}$ and $\hat{U}_{r}$ are the volume flow rates into and out of the junction, and $V_{J}$ is the cylindrical volume of the junction.

Three minor losses—between the duct and stack, stack and duct, and duct and radial cavity—are also accounted for in the linear model in the form of pressure jumps

[TABLE]

where we have adopted the expansion and contraction formulas of Borda-Carnot and Idelchik minor losses:

[TABLE]

where $A_{0}$ and $A_{1}$ are the smaller and larger areas, respectively. The pressure drop eq. 26 is linearized about a given acoustic velocity amplitude, $U$, which is updated in time based on growth rate information, $U=U_{0}\exp\left(\alpha t\right)$, where $U_{0}$ is an initial velocity amplitude.

Further implementation details can be found in Lin et al.[16] Linear stability analysis, in the present model, does not incorporate bulk viscosity effects.

5 Results

The thermoacoustic engine model (see fig. 1) without bulk viscosity effects (the reference case) and with various levels of bulk viscosity, as tuned by working fluid humidity, were simulated through transient growth and thermoacoustic amplification into the first limit cycle. Since the effective bulk viscosity implies a direct effect on the pressure gradient and has a non-trivial relationship with pressure, it is also expected that Navier–Stokes computational results will be highly dependent on humidity levels.

We first present a pressure amplitude timeseries of the results without a bulk viscosity model, as shown in fig. 5. We also have overlaid transient results from the linear model, which reaches limit cycle due to nonlinearities from minor losses. The time series for the linear model is integrated along its evolution over $\exp\left(\left(\alpha+i\omega\right)t\right)$. While the limit cycle is well-predicted, as minor losses increase with larger flow rate amplitude, the transient growth rate is overpredicted by the linear model.

Acoustic amplitudes along the engine, as defined through the tubes, stack, and radial cavity in fig. 6, were also computed for the reference model. These are shown in fig. 6 and qualitatively compare well.

Three additional cases were run, corresponding to a working fluid of atmospheric air with 1%, 5%, and 20% relative humidity. The transients for each of these are shown in fig. 7.

Notably, the reference case has the greatest growth rate and limit cycle amplitude, approximately 3503 Pa. This is to be expected, as the attenuation from bulk viscosity is not accounted for in the reference case. The limit cycle amplitude of the 5% relative humidity case is calculated to be approximately 3304 Pa, which is significantly higher than that of the 20% humidity case (3208 Pa) and that of the 1% humidity case (3054 Pa).

At the limit cycle, the Navier–Stokes model also reveals the presence of imperfect Helmholtz resonator behavior. Compressibility within the system is not neglected, and acoustic energy fluctuates through not only the quarter-wavelength engine sections but also through the radial cavity.

Streaming patterns appear similar to other standing-wave engines thus far studied,[16] and thermoacoustic heat transport away from the stack qualitatively is lower than expected (not shown).

6 Discussion

6.1 Effective acoustic pressure

Bulk viscosity has the largest effect near vibrational energy “peaks,” where the resonance frequency approaches the natural frequency of a vibrational mode. In order to evaluate the effect of bulk viscosity under different conditions, we consider its effect on the acoustic effective pressure. The physical manifestation of the bulk viscosity within the momentum and energy equations, eq. 1, is as an adjustment to the thermodynamic pressure, i.e.:

[TABLE]

Per linear acoustics and assuming simple oscillations for which pressure and velocity are in standing-wave or traveling-wave phasing,

[TABLE]

hence resulting in

[TABLE]

This suggests that the dimensionless group $\frac{\mu_{B}\omega}{\gamma p_{0}}$ is a measure of the relative importance of bulk viscosity effects on the effective acoustic pressure.

Traditional attenuation curves, measured relative to attenuation per meter, tend to belie the effect of bulk viscosity at high frequencies. When measured relative to acoustic wavelength, the bulk viscosity contribution to attenuation is shown to be as large as 1% of pressure amplitude and has a magnitude peak varying with gas temperature, pressure, and humidity, as seen in fig. 8.

6.2 Contribution of bulk viscosity

In the previous section, it was shown that the growth rate and limit cycle amplitude vary non-monotonically with the relative working fluid humidity. Figure 8b suggests one possible explanation of the phenomenon. The bulk viscosity contribution to attenuation at the base temperature of 300 K increases monotonically with humidity (for $h_{r}<20\%$); as a result, the 20% relative humidity case has the highest attenuation in the regions of the engine which are at ambient temperature. The bulk viscosity contribution to attenuation at the peak temperature of 600 K decreases with humidity; as a result, the 1% relative humidity case has the highest attenuation in the heated portion of the engine. For the optimal case, with 5% relative humidity, attenuation in the heated and ambient temperature portions of the engine are reduced.

The phase difference between pressure fluctuations and the velocity field divergence (dilatation) term contributes directly to both acoustic wave attenuation and thermoacoustic energy production. For both purely standing- and traveling-waves, pressure fluctuation and dilatation are 90 degrees out of phase. In the case of the presented standing-wave engine, it is found that such a phase difference holds everywhere but in the thermoacoustic stack, where it approaches zero (fig. 9). This is consistent with interpretations of standing-wave thermoacoustic energy production, in which, due to the imposed temperature gradient, a fluid parcel oscillates between high pressure and expansion and low pressure and contraction.

It is yet unclear whether bulk viscosity directly affects thermoacoustic energy conversion. However, these results suggest that there is an opportunity for miniature thermoacoustic engines to be designed with bulk viscosity effects in mind. Tracing the Lagrangian fluid parcel within the stack would suggest that bulk viscosity attenuation is minimal within the stack channel and has a much larger effect in the rest of the engine. The limit cycle pressure amplitude differs by over 7% between cases which incorporate bulk viscosity, suggesting a non-negligible effect on the acoustic power output.

6.3 Summary and Future Work

Several results are currently presented. A Navier–Stokes code has been shown to accurately capture acoustic absorption under a range of frequencies relevant to both acoustic and thermoacoustic applications. A method for evaluating the attenuation strength of bulk viscosity in different conditions is presented. The use of bulk viscosity adjustments via temperature, pressure, and humidity provides a computational baseline. Further investigation will reveal how relaxation at high-frequency affects thermoacoustic onset, relative to the baseline. Results for thermoacoustic amplification have been presented, using a modified version of the Flitcroft and Symko engine. As suggested in the paper, limit cycle results also differ significantly between numerics which account for bulk viscosity and for solvers which do not account for bulk viscosity.

Continuing work include the computation and analysis of additional cases, a generalized predictive model for the effect of $\mu_{B}$ on thermoacoustic onset and limit cycle behavior, and visual intuition for how an effective bulk viscosity can affect thermoacoustic behavior. These results are expected to conclude in techniques for optimization of high-frequency engines.