# Modeling of thermoacoustic systems using the nonlinear frequency domain method

A A de Jong, A H Wijnant, A A Wilcox, A A de Boer

► To cite this version:

A A de Jong, A H Wijnant, A A Wilcox, A A de Boer. Modeling of thermoacoustic systems using the nonlinear frequency domain method. Journal of the Acoustical Society of America, 2015, 138, pp.1241-1252. 10.1121/1.4928301. hal-01443188

HAL Id: hal-01443188

https://hal.science/hal-01443188v1

Submitted on 22 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Modeling of thermoacoustic systems using the nonlinear frequency domain method

J. A. de Jong,¹ Y. H. Wijnant,¹ D. Wilcox,² and A. de Boer¹

¹Structural Dynamics and Acoustics, Faculty of Engineering Technology, University of Twente, Enschede, The Netherlands

²Chart Incorporated—Qdrive, 302 Tenth Street, Troy, New York 12180, USA

When modeling thermoacoustic (TA) devices at high amplitude, nonlinear effects such as time-average mass flows, and the generation of higher harmonics can no longer be neglected. Thus far, modeling these effects in TA devices required a generally computationally costly time integration of the nonlinear governing equations. In this paper, a fast one-dimensional nonlinear model for TA devices is presented, which omits this costly time integration by directly solving the periodic steady state. The model is defined in the frequency domain, which eases the implementation of phase delays due to viscous resistance and thermoacoustic heat exchange. As a demonstration, the model is used to solve an experimental standing wave thermoacoustic engine. The obtained results agree with experimental results, as well as with results from a nonlinear time domain model from the literature. The low computational cost of this model opens the possibility to do optimization studies using a nonlinear TA model.

## I. INTRODUCTION

When designing thermoacoustic (TA) systems, fast computational models are required for quick feedback on design decisions. However, fast computational models are often based on some crude simplifications of the underlying physics. Hence, generally a trade-off has to be made between the accuracy of the physics and the computational cost. On the one hand, we have well-established one-dimensional linear TA theory,¹,² which is fast to solve, but only valid for low amplitude, while on the other hand, multi-dimensional computational fluid dynamics (CFD) simulations offer deep insight into the physical processes of a TA system, but these simulations demand significantly more time to solve.

From a computational time perspective, simulating a complete TA system with CFD is considered too costly and therefore not suitable for optimization purposes. Hence, additional assumptions are often applied. In the work of Hamilton *et al.*,³ the Courant number restriction on the time step is loosened by neglecting the momentum equation in the direction perpendicular to the wave propagation direction. In other multi-dimensional modeling efforts, a low Mach number analysis⁴,⁵ has been applied and/or only a subdomain of the full system⁶ is taken into account. For the interaction with the part of the system not modeled in detail, typically a linear acoustic model is used, or linear acoustic boundary conditions are applied. Consequently, nonlinear effects are neglected in that part of the system which is not modeled.

The main focus of this paper is the development of a nonlinear one-dimensional TA model. The research on such

a model has been started by Watanabe *et al.*⁷ and Yuan *et al.*⁸ With this time domain model, they were able to simulate a thermoacoustic system starting from initial conditions, up to the periodic steady-state. The onset of a self-excited system was triggered with an initial pressure disturbance. Initially, their model suffered from spurious instabilities, caused by unstable linear eigenmodes which are introduced by the viscous resistance and heat transfer models. To suppress these instabilities, artificial diffusion was added. Later, Karpov and Prosperetti⁹ solved this problem with better models for resistance and heat transfer, which they call the “exchange terms” as they are responsible for exchanging mass and momentum with the surrounding solid. Still, their exchange terms are only in agreement with linear theory for the first and second harmonic of the solution’s spectrum. The frequencies of these harmonics needed to be determined before starting the simulation, which makes their model, in some sense, dependent on global (tunable) parameters.

A more recent work discusses the application of the method of characteristics¹⁰ to model the experimental standing wave thermoacoustic engine of Swift.¹¹ Here, heat transfer and viscous resistance is computed using semi-empirical models.

The model presented in this work is based on the quasi-one-dimensional model for thermoacoustic devices used by Karpov and Prosperetti.⁹ Computational time is diminished by directly solving the periodic steady-state. This is facilitated by a reformulation in the frequency domain using the nonlinear frequency domain (NLFD) method.¹²,¹³ Earlier work has shown that the NLFD method is able to solve the periodic steady-state of a one-dimensional acoustic system with a computational time of ~1 s. With the NLFD method, all physical quantities are assumed to be periodic in time,
hence, they can be described by a Fourier series. Then, instead of doing a time integration, all Fourier coefficients are simultaneously solved. Using the Fourier coefficients, the periodic time domain response can be reconstructed by applying an inverse discrete Fourier transform. An additional advantage of a frequency domain formulation is that the frequency domain exchange models of linear TA theory can be used in unmodified form, which makes the model in exact agreement with linear thermoacoustic theory. A drawback of the method is that the initial transients of a system cannot be modeled. However, in most cases of engine and refrigerator design, the interest lies primarily on the performance of a TA system in steady-state.

In Sec. II, the model is explained, including the application of the NLFD method. Furthermore, a frequency domain form is derived for the exchange terms. Then, in Sec. III, the numerical method used to solve the model of Sec. II will be described. Notably, several problems accompanying the NLFD method, which are related to solving for closed self-excited systems (TA prime movers), are treated. Then, in Sec. IV, the model will be applied to solve the steady-state of an experimental standing wave TA engine.¹⁴ The results are discussed, and finally in Sec. V conclusions are drawn.

## II. MODEL

### A. Governing equations

As shown by, for example, Watanabe et al.,⁷ a quasi-one-dimensional nonlinear continuity, momentum, and energy equation can be derived from the governing equations. For convenience, the model is restated here,

$$ \frac{\partial \rho}{\partial t} + \frac{1}{S_f} \frac{\partial}{\partial x} (S_f \rho u) = 0, \quad (1a) $$

$$ \frac{\partial}{\partial t}(\rho u) + \frac{1}{S_f} \frac{\partial}{\partial x}(S_f \rho u^2) + \frac{\partial p}{\partial x} = -\mathcal{R}_u, \quad (1b) $$

$$ \frac{\partial}{\partial t} \left( \frac{p}{\gamma-1} + \frac{1}{2} \rho u^2 \right) + \frac{1}{S_f} \frac{\partial}{\partial x} \left[ u S_f \left( \frac{\gamma}{\gamma-1} p + \frac{1}{2} \rho u^2 \right) \right] = \mathcal{H}_u - \mathcal{Q}_u \frac{dT_w}{dx}, \quad (1c) $$

accompanied by the perfect gas equation of state

$$ p = \rho R_s T, \quad (1d) $$

where $x$ is the wave propagation direction, $t$ is time, $\rho$ density, $u$ (axial) velocity, $p$ pressure, $T$ temperature, $\gamma$ the ratio of specific heats, $R_s$ the specific gas constant, and $S_f$ is the fluid cross-sectional area. So $S_f = \phi S$, where $\phi$ is the porosity and $S$ the total cross-sectional area. $\mathcal{R}_u$ is the resistive drag force [N m⁻³] in the wave propagation direction due to shear stress from the fluid acting on the wall. It is by definition,

$$ \mathcal{R}_u \equiv -\frac{1}{S_f} \oint_{\Pi} (\boldsymbol{\tau} \cdot \mathbf{n})_x d\Pi, \quad (2) $$

where $\boldsymbol{\tau}$ is the viscous stress tensor and $\mathbf{n}$ is the unit normal pointing transversely from the fluid in the direction of the

solid. The transverse direction is perpendicular to the wave propagation direction. $\Pi$ is the “wetted perimeter,” the contact line between solid and fluid. The two terms on the right side of Eq. (1c) together describe heat flow from the solid to the fluid,

$$ \mathcal{H}_u - \mathcal{Q}_u \frac{dT_w}{dx} \equiv -\frac{1}{S_f} \oint_{\Pi} (\mathbf{q} \cdot \mathbf{n}) d\Pi, \quad (3) $$

where $\mathbf{q}$ is the heat flux vector.

Unless otherwise stated, all quantities in Eq. (1) are area-averaged. With the emphasis on modeling of TA systems, Watanabe et al.⁷ also provide an overview of the major assumptions behind this model. In the energy equation, two sources are defined for heat exchange with the wall. The source $\mathcal{H}_u$, is due to transverse conduction as a result of a difference between the wall temperature and the fluid area-averaged temperature. The source $\mathcal{Q}_u(dT_w/dx)$ is due to convective enthalpy flow, originating from the velocity and temperature boundary layer close to the wall. In Sec. II C, the momentum and energy exchange terms $\mathcal{R}_u$, $\mathcal{H}_u$, and $\mathcal{Q}_u$ are treated in detail. We assume that the temperature of the wall is constant in time, but a known function of space: $T_w = T_w(x)$.

### B. Application of the NLFD method

We assume a periodic steady state, so all dependent variables can be described by a Fourier series. The basis of the NLFD method¹³ is to substitute a truncated Fourier series for each dependent variable in Eqs. (1) and to solve for the Fourier coefficients. To explain this method, we take the density as the example dependent variable. The truncated Fourier series for the density is

$$ \rho(x, t) = \hat{\rho}_{c,0}(x) + \Re \left[ \sum_{n=1}^{N_f} \hat{\rho}_{c,n}(x) e^{in\omega t} \right], \quad (4) $$

where $\omega$ is the fundamental radial frequency, $n$ the harmonic index, $N_f$ is the total number of harmonics taken into account, $\Re[\cdots]$ denotes the real part of its argument, and $\hat{\rho}_{c,i}$ with $0 \le i \le N_f$ is the Fourier coefficient (phasor) for harmonic $i$. We store all Fourier coefficients in a real-valued vector $\hat{\rho}$ as

$$ \hat{\rho} = \{ \hat{\rho}_{c,0} \Re[\hat{\rho}_{c,1}] \Im[\hat{\rho}_{c,1}] \Re[\hat{\rho}_{c,2}] \Im[\hat{\rho}_{c,2}] \cdots \Im[\hat{\rho}_{c,N_f}] \}^T, \quad (5) $$

where $\Im[\cdots]$ denotes the imaginary part of its argument. The imaginary part of the time-averaged Fourier component ($\hat{\rho}_{c,0}$) is zero and is therefore not taken into account in the storage vector. From the vector of density Fourier coefficients, a vector of time instances $\rho$ can be constructed by applying the inverse discrete Fourier transform (DFT),

$$ \rho_m \equiv \rho(t_m) = \hat{\rho}_0 + \Re \left[ \sum_{n=1}^{N_f} (\hat{\rho}_{2n-1} + i\hat{\rho}_{2n}) e^{inm\omega\Delta t} \right], \quad 0 \le m < 2N_f + 1, \quad (6) $$
where $t_m = m\Delta t$ and $\Delta t$ is the time step between each time instance, which can be computed from the fundamental frequency and the number of harmonics with

$$ \Delta t = \frac{2\pi}{\omega(2N_f + 1)}. \qquad (7) $$

The operation of Eq. (6) applied to all time instance simultaneously is written as

$$ \rho = \mathcal{F}^{-1} \cdot \hat{\rho}, \qquad (8) $$

where $\mathcal{F}$ is referred to as the discrete Fourier transform matrix, and hence $\mathcal{F}^{-1}$ the backward discrete Fourier transform matrix. For the continuity equation, we also need the time derivative of the density. For the periodic function of Eq. (4), the time derivative can be decomposed in terms of the time derivatives of the individual complex exponentials. So the vector of the time derivatives of the density can be evaluated as

$$ \frac{\partial}{\partial t} \rho = \mathcal{F}^{-1} \cdot \omega \check{\mathbf{D}} \cdot \hat{\rho}, \qquad (9) $$

where $\check{\mathbf{D}}$ is the “time-derivative matrix,” scaled with $\omega$,

$$ \omega \check{\mathbf{D}} = \omega \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & \cdots \\ 0 & 0 & -1 & 0 & 0 & \cdots \\ 0 & 1 & 0 & 0 & 0 & \cdots \\ 0 & 0 & 0 & 0 & -2 & \cdots \\ 0 & 0 & 0 & 2 & 0 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots \end{bmatrix}. \qquad (10) $$

Thus, the time derivative of the density in the frequency domain is

$$ \frac{\partial \rho}{\partial t} = \mathcal{F} \cdot \omega \mathcal{F}^{-1} \cdot \check{\mathbf{D}} \cdot \hat{\rho} = \omega \check{\mathbf{D}} \cdot \hat{\rho}. \qquad (11) $$

To compute the spatial derivative of the mass flux, the nonlinear product $m = \rho u$ is required. This is done by first applying the backward discrete Fourier transform to the density and velocity vectors of Fourier coefficients, then computing the product of density and velocity at each time instance, and finally transforming back into the frequency domain. Accordingly, the Fourier coefficients of the mass flux are computed as

$$ \hat{\mathbf{m}} = \mathcal{F} \cdot ((\mathcal{F}^{-1} \cdot \hat{\rho}) \circ (\mathcal{F}^{-1} \cdot \hat{\mathbf{u}})), \qquad (12) $$

where the operator $\circ$ denotes the elementwise product of two vectors$^{15}$ yielding a vector.

We apply this method to Eqs. (1a)–(1d). As dependent variables the frequency domain storage vectors for cross-sectional area averaged density ($\hat{\rho}(x)$), volume flow ($\hat{\mathbf{U}}(x)$), temperature ($\hat{\mathbf{T}}(x)$), and pressure ($\hat{\mathbf{p}}(x)$) are chosen. The continuity, momentum, and energy equation are multiplied with the fluid cross-sectional area. For the continuity equation, we find

$$ S_f \omega \check{\mathbf{D}} \cdot \hat{\rho} + \frac{d}{dx} \mathcal{F} \cdot (\rho \circ \mathbf{U}) = \mathbf{0}. \qquad (13) $$

For brevity, the inverse DFT inside the parentheses is not shown. Therefore, for the time domain vectors $\rho$ and $\mathbf{U}$, one should read $\mathcal{F}^{-1} \cdot \hat{\rho}$ and $\mathcal{F}^{-1} \cdot \hat{\mathbf{U}}$, respectively. The volume flow can be related to the velocity by

$$ U \equiv uS_f. \qquad (14) $$

For the NLFD momentum equation we find

$$ \omega \check{\mathbf{D}} \cdot \mathcal{F} \cdot (\rho \circ \mathbf{U}) + \frac{d}{dx} \mathcal{F} \cdot \left( \frac{\rho \circ \mathbf{U}^2}{S_f} \right) + S_f \frac{d\hat{\mathbf{p}}}{dx} = -\hat{\mathcal{R}}_U, \qquad (15) $$

and for the NLFD form of the energy equation

$$ \omega \check{\mathbf{D}} \cdot \frac{S_f \hat{\mathbf{p}}}{\gamma - 1} + \mathcal{F} \cdot \frac{1}{2S_f} (\rho \circ \mathbf{U}^2) - \frac{d}{dx} \mathcal{F} \cdot \left( S_f \kappa \circ \frac{d\mathbf{T}}{dx} \right) \\ + \frac{d}{dx} \mathcal{F} \cdot \left( \mathbf{U} \circ \frac{\gamma}{\gamma - 1} \mathbf{p} + \frac{1}{2S_f^2} (\rho \circ \mathbf{U}^2) \right) = \hat{\mathcal{H}}_U - \hat{\mathcal{Q}}_U \frac{dT_w}{dx}, \qquad (16) $$

where $\kappa$ is the thermal conductivity, and hence $\mathbf{\kappa}$ the vector of time instances for the thermal conductivity. Equation (16) is not the full NLFD-equivalent variant of Eq. (1c), because we have added an extra term corresponding to axial conduction. As stated by Karpov et al.,$^9$ this term is only important near walls. The reason for adding this term is explained in Appendix A. The last equation to solve for is the NLFD perfect gas equation of state

$$ \hat{\mathbf{p}} = R_s \mathcal{F} \cdot (\rho \circ \mathbf{T}). \qquad (17) $$

In Eqs. (15) and (16), viscous resistance and heat transfer terms are redefined such that $\hat{\mathcal{R}}_U = S_f \hat{\mathcal{R}}_u$, $\hat{\mathcal{H}}_U = S_f \hat{\mathcal{H}}_u$, and $\hat{\mathcal{Q}}_U = S_f \hat{\mathcal{Q}}_u$, where the exchange terms with lower case $u$ are defined as the Fourier transform of the resistance and heat transfer operators of Eqs. (1a)–(1d).

### C. The exchange terms

We assume laminar, fully developed—oscillating—flow, so the resistance $\hat{\mathcal{R}}$ is a linear function of the volume flow. A closed form laminar resistance coefficient $\hat{\mathbf{R}}$ in frequency domain is derived from the linearized momentum equation, as is explained by Watanabe et al.$^7$ Applying the NLFD method to this model results in

$$ \hat{\mathcal{R}}_U = \hat{\mathbf{R}} \diamond \hat{\mathbf{U}}, \qquad (18) $$

where $\hat{\mathbf{R}}$ is the frequency domain resistance coefficient. The symbol $\diamond$ denotes the product in the frequency domain. Since the real and imaginary part of the value at each frequency is stacked in the form of Eq. (5), an equivalent matrix-vector product is generated. For example, for $\hat{\mathbf{R}}$ this matrix is
$$
M(\hat{\mathbf{R}}) = \begin{bmatrix}
\hat{R}_0 & 0 & 0 & 0 & 0 & \cdots \\
0 & \Re(\hat{R}_1) & \Im(-\hat{R}_1) & 0 & 0 & \cdots \\
0 & \Im(\hat{R}_1) & \Re(\hat{R}_1) & 0 & 0 & \cdots \\
0 & 0 & 0 & \Re(\hat{R}_2) & \Im(-\hat{R}_2) & \cdots \\
0 & 0 & 0 & \Im(\hat{R}_2) & \Re(\hat{R}_2) & \cdots \\
\vdots & \vdots & \vdots & \vdots & \vdots & \ddots
\end{bmatrix} \tag{19}
$$

Then, the frequency domain product $\hat{\mathbf{R}} \diamond \hat{\mathbf{U}}$ is computed as $M(\hat{\mathbf{R}}) \cdot \hat{\mathbf{U}}$, where the operator $\cdot$ is the matrix-vector product. The components of $\hat{\mathbf{R}}$ are defined such that

$$
\hat{R}_{2n-1} + i\hat{R}_{2n} = \frac{\hat{\mu}_0}{r_h^2} \frac{is_n^2 f_{\nu,n}}{1 - f_{\nu,n}} \quad \text{for } n > 0, \qquad (20)
$$

where $\hat{\mu}_0$ is the time-averaged dynamic viscosity, $f_{\nu,n} = f_{\nu,n}(s_n)$ is the cross-sectional geometry-dependent viscous Rott function$^{1,2}$ evaluated at shear wave number $s_n$. $s_n$ is the shear wave number for harmonic $n$, given as

$$
s_n = r_h \sqrt{\frac{\hat{\rho}_0 \omega n}{\hat{\mu}_0}}, \qquad (21)
$$

where $r_h$ is the hydraulic radius. It is defined as the cross-sectional area divided by the wetted perimeter

$$
r_h = \frac{S_f}{\Pi}. \tag{22}
$$

The shear wave number can be related to the viscous penetration depth for harmonic $n$ as

$$
s_n \equiv \sqrt{2} \frac{r_h}{\delta_{\nu,n}}. \qquad (23)
$$

If one would use this definition of the viscous resistance in the time domain, computing the viscous drag would require the computation of the convolution integral of the inverse Fourier transform of the frequency domain resistance coefficient with the volume flow. However, performing this integration over all cycles necessary to reach steady-state is completely impractical.$^7$ Moreover, for simple geometries, such as circular tubes and parallel plates, no closed form solution exists for the inverse Fourier transform of the resistance coefficient. These are important reasons for choosing a frequency domain formulation rather than a time domain formulation.

In the drag model it is assumed that nonlinear effects in the viscous boundary layer are small, such that the boundary layer flow can be regarded as a superposition of linear boundary layer velocity profiles. Hence, the drag coefficient is based on the following transverse velocity distribution:

$$
u(x, \mathbf{r}, t) = \Re \left[ \sum_{n=0}^{N_f} \langle \hat{u}_n \rangle(x) \frac{1 - h_{\nu,n}(x, \mathbf{r})}{1 - f_{\nu,n}(x)} e^{in\omega t} \right], \quad (24)
$$

where $\mathbf{r}$ is the transverse coordinate, $\langle \cdots \rangle$ denotes area-averaging, and $h_{\nu,n}$ is the transverse velocity shape function, from which the Rott functions are derived:

$$
\langle h_{\nu,n} \rangle = f_{\nu,n}. \qquad (25)
$$

For example, for parallel-plate pores, $h_{\nu,n}$ is

$$
h_{\nu,n}(s_n, y) = \frac{\cosh\left(s_n \sqrt{i} \frac{y}{r_h}\right)}{\cosh\left(s_n \sqrt{i}\right)}, \qquad (26)
$$

where $y$ is the distance from the center line between two plates.

Equation (24) implies that the transverse “modes” of the axial velocity distribution do not mix, hence in the boundary layer the linear (viscous) effects dominate.

As $\hat{R}_n$ becomes singular for $n=0$, we evaluate $\hat{R}_0$ by computing

$$
\hat{R}_0 = \lim_{n\omega \to 0} \frac{\hat{\mu}_0 is_n^2 f_{\nu,n}}{r_h^2 (1 - f_{\nu,n})}. \qquad (27)
$$

This procedure is done for the viscous Rott functions for the earlier mentioned parallel-plate pores, for cylindrical tubes, for which the viscothermal Rott functions are defined in terms of Bessel functions,$^{16,17}$ and for the boundary layer approximation. In Table I, the computed limits are listed.

For the heat exchange terms, $\mathcal{H}_U$ and $\hat{Q}_U$, the same procedure is followed as for the viscous resistance. For $\mathcal{H}_U$ we obtain

$$
\hat{\mathcal{H}}_U = S_f \hat{\mathbf{H}} \diamond (T_w \hat{\mathbf{I}} - \hat{\mathbf{T}}), \qquad (28)
$$

where

$$
\hat{H}_{2n-1} + i\hat{H}_{2n} = \frac{\mathrm{Pr}\hat{\kappa}_0}{r_h^2} \frac{is_n^2 f_{\kappa,n}}{1 - f_{\kappa,n}}, \qquad (29)
$$

where $\mathrm{Pr}$ is the Prandtl number and $f_{\kappa,n}$ is the thermal Rott function. The thermal Rott function can be evaluated as the viscous Rott function at shear wave number $\sqrt{\mathrm{Pr}s_n}$:

$$
f_{\kappa,n} = f_{\kappa}(s_n, \mathrm{Pr}) = f_{\nu}(\sqrt{\mathrm{Pr}s_n}). \qquad (30)
$$

The vector $\hat{\mathbf{I}}$ is

$$
\hat{\mathbf{I}} \equiv \left\{ \begin{matrix} 1 \\ 0 \\ 0 \\ \vdots \end{matrix} \right\} = \mathcal{F} \cdot 1 \equiv \mathcal{F} \cdot \left\{ \begin{matrix} 1 \\ 1 \\ 1 \\ \vdots \end{matrix} \right\}, \qquad (31)
$$

so it only acts on the time-averaged part of the heat transfer. Similarly to Eq. (28), for $\hat{Q}$ we find

<table><caption>TABLE I. Derived low-frequency limits for the thermoacoustic exchange terms.</caption>
  <thead>
    <tr>
      <th></th>
      <th><i>f</i><sub>v</sub></th>
      <th><i>R̂</i><sub>0</sub></th>
      <th><i>H̃</i><sub>0</sub></th>
      <th><i>Q̂</i><sub>0</sub></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>Parallel-plate pores</td>
      <td>
        <div style="text-align: center;">
          <div style="border-bottom: 1px solid black; padding-bottom: 2px;">tanh(√is<sub>n</sub>)</div>
          <div style="padding-top: 2px;">√is<sub>n</sub></div>
        </div>
      </td>
      <td>
        <div style="text-align: center;">
          <div style="border-bottom: 1px solid black; padding-bottom: 2px;">3μ̃<sub>0</sub></div>
          <div style="padding-top: 2px;">r<sup>2</sup><sub>h</sub></div>
        </div>
      </td>
      <td>
        <div style="text-align: center;">
          <div style="border-bottom: 1px solid black; padding-bottom: 2px;">3κ̃<sub>0</sub></div>
          <div style="padding-top: 2px;">r<sup>2</sup><sub>h</sub></div>
        </div>
      </td>
      <td>
        <div style="text-align: center;">
          <div style="border-bottom: 1px solid black; padding-bottom: 2px;">1/5</div>
          <div style="padding-top: 2px;">ρ̃<sub>0</sub>c<sub>p</sub></div>
        </div>
      </td>
    </tr>
    <tr>
      <td>Cylindrical pores</td>
      <td>
        <div style="text-align: center;">
          <div style="border-bottom: 1px solid black; padding-bottom: 2px;">J<sub>1</sub>(2i√is<sub>n</sub>)</div>
          <div style="padding-top: 2px;">i√iJ<sub>0</sub>(2i√is<sub>n</sub>)</div>
        </div>
      </td>
      <td>
        <div style="text-align: center;">
          <div style="border-bottom: 1px solid black; padding-bottom: 2px;">2μ̃<sub>0</sub></div>
          <div style="padding-top: 2px;">r<sup>2</sup><sub>h</sub></div>
        </div>
      </td>
      <td>
        <div style="text-align: center;">
          <div style="border-bottom: 1px solid black; padding-bottom: 2px;">2κ̃<sub>0</sub></div>
          <div style="padding-top: 2px;">r<sup>2</sup><sub>h</sub></div>
        </div>
      </td>
      <td>
        <div style="text-align: center;">
          <div style="border-bottom: 1px solid black; padding-bottom: 2px;">1/3</div>
          <div style="padding-top: 2px;">ρ̃<sub>0</sub>c<sub>p</sub></div>
        </div>
      </td>
    </tr>
    <tr>
      <td>Boundary layer approximation</td>
      <td>
        <div style="text-align: center;">
          <div style="border-bottom: 1px solid black; padding-bottom: 2px;">-i√i</div>
          <div style="padding-top: 2px;">s<sub>n</sub></div>
        </div>
      </td>
      <td>0</td>
      <td>0</td>
      <td>-ρ̃<sub>0</sub>c<sub>p</sub></td>
    </tr>
  </tbody>
</table>
$$ \hat{\mathcal{Q}}_U = \hat{\mathcal{Q}} \diamond \hat{\mathbf{U}}, \qquad (32) $$

where

$$ \hat{\mathcal{Q}}_{2n-1} + i\hat{\mathcal{Q}}_{2n} = \frac{\hat{\rho}_0 c_p}{1 - \mathrm{Pr}} \left( \frac{f_{\nu,n}}{1 - f_{\nu,n}} - \mathrm{Pr} \frac{f_{\kappa,n}}{1 - f_{\kappa,n}} \right). \quad (33) $$

In Eq. (33), $c_p$ is the specific heat at constant pressure. The results for $\hat{H}_0$ and $\hat{Q}_0$ are also listed in Table I. Note that if we linearize Eqs. (13) and (15), (16) with Eqs. (18), (28), and (33) for the thermoacoustic exchange terms and neglect axial conduction, the model reduces to the exact linearized model of Swift and Rott, except that multiple non-interacting harmonics can be computed simultaneously. Likewise, for an inviscid non-conducting gas, $\hat{\mathbf{R}}_U$, $\hat{\mathbf{H}}_U$, and $\hat{\mathcal{Q}}_U$ all reduce to the null vector and we obtain the quasi-one-dimensional Euler equations for a calorically perfect gas.

It is striking that the boundary layer approximation of $\hat{Q}_0$ in Table I has a different sign than the limit for parallel-plate and cylindrical pores. Nevertheless, this solution is in contradiction with its underlying assumption. This assumption states that the boundary layer thickness is small compared to the typical cross-sectional length scale, yet the boundary layer thickness for fully developed steady flow is infinite by definition, because as $\omega \to 0$, $\delta_{\nu,\kappa} \to \infty$. Therefore, this boundary layer approximation limit of Table I is nonphysical and should not be used to model time-averaged interaction with the wall. Note that this does not mean that in general the boundary layer approximation is invalid, but that the boundary layer approximation is inconsistent for the steady (time-average) frequency component. For the case of Sec. IV, the cylindrical resonator and hot end are modeled using the boundary layer approximation for $n > 0$, while for $n = 0$, the cylindrical pore limit is used.

A problem with the heat exchange terms is that the model for wall heat exchange significantly underestimates the time-averaged heat transfer between the solid and the fluid. This is because as a nonlinear effect, the presence of an acoustic wave affects the time-averaged transverse temperature profile, which in effect influences the time-averaged heat transfer. In previous work, this effect is modeled in the framework of second order thermoacoustic theory.18,19 Additionally, empirical models have been developed to compute time-averaged heat transfer.20 The latter uses a clear distinction between oscillating heat transfer by the acoustic wave and induced time-averaged heat transfer. The oscillating heat transfer is modeled using linear thermoacoustic theory, and time-averaged heat transfer with an empirical model. The empirical model needs statistical input from the acoustic field, such as the root mean square (RMS) acoustic Reynolds number. A fully nonlinear model does not allow such a distinction, and the model presented in this section does not take such nonlinear heat transfer effects into account. Inclusion of these effects requires a complete revision of the current thermoacoustic heat exchange terms. At this point, this is left for future research.

### III. NUMERICAL IMPLEMENTATION

In this section, the numerical method used to solve the model of Sec. II will be explained. For the spatial discretization, a staggered cell-centered finite volume method is used.21 For constant grid spacing this results in second order spatial discretization. The volume flow is evaluated at the cell walls and density, temperature, and pressure at the vertices.

For self-excited thermoacoustic systems, the fundamental frequency is unknown and part of the solution. This introduces an extra degree of freedom for the numerical system, which requires solving an extra equation. We add the additional equation that fixes the phase of the first harmonic at one of the walls of the system. For example, for the left wall we add

$$ \Im(\hat{p}_1(x=0)) = 0 \qquad (34) $$

to the system of equations. The details of this implementation are explained in Appendix B.

In the system of equations, one time-averaged continuity equation is replaced by the global conservation of mass equation

$$ \sum_{j=0}^{N_{\text{cells}}-1} \hat{p}_{0,j} V_{f,j} = m_p, \qquad (35) $$

where $j$ is the cell index, $N_{\text{cells}}$ is the total number of cells in the system, $V_{f,j}$ is the local volume of a cell and $m_p$ the prescribed (fluid) mass in the system, which is computed as

$$ m_p = \frac{p_{\text{ref}}}{R_s T_{\text{ref}}} \sum_{j=0}^{N_{\text{cells}}-1} V_{f,j}, \qquad (36) $$

where $p_{\text{ref}}$ is the reference pressure (which is set when the system is initially filled with gas), and $T_{\text{ref}}$ is the reference temperature. Equation (35) guarantees that the total mass in the system is conserved. The necessity of this replacement is explained in Appendix A.

The inclusion of thermal conduction in the model results in thermal boundary layers of thickness $\sim \delta_k$ at the tube ends where the temperature is prescribed. One such temperature boundary condition is always present, as at least one Dirichlet boundary condition for the temperature is needed to make the solution unique. The small thermal boundary layer is responsible for so-called thermal hysteresis (thermal relaxation),22 which dissipates acoustic power. To resolve the acoustic temperature in the boundary layers, the grid spacing needs to be $< \delta_k$. Therefore, a local grid refinement is made at the expected location of boundary layers. This boundary grid stretches exponentially to the inner segment (standard) grid spacing. This way, boundary layer effects can be resolved while the number of grid-points stays manageable.

The computations are done on a standard desktop computer with Intel i7 860 cores running at a clock frequency of 2.80 GHz with 8192 kB of cache size. A C++ program has been written, which runs serially and uses the ARMADILLO
<figure><img src="image_2.png" /><figcaption>FIG. 1. (Color online) Schematic overview of the standing wave engine of Atchley. The picture is not to scale. The wall temperature profile is shown in the graph below. On the right side of the stack, $T_w$ is set to $T_H$.</figcaption></figure>

linear algebra package²³ for computing the Jacobian terms. For solving the linear system at each Newton iteration, the direct SPARSELU solver of the EIGEN linear algebra package²⁴ is used. Typical solution times range from 1 to 5 min, which is fast compared to computation times of ~10 h when using the method of characteristics.¹⁰

## IV. APPLICATION: STANDING WAVE TA ENGINE

In this section we introduce the application of the NLFD method to solve the periodic steady-state of a standing wave thermoacoustic engine. As the application test case, the standing wave engine of Atchley¹⁴ is chosen, because we have both experimental data of the linear and nonlinear aspects of this device. Moreover, using the results of this test case allows us to compare them with earlier numerical results, presented by Karpov et al.⁹ and Yuan et al.⁸

Figure 1 shows a schematic overview of the geometry of this standing wave engine. The geometry consists of discrete “segments”: the resonator, cold heat exchanger, stack, hot heat exchanger, and the hot duct. The geometrical and modeling details for each segment are listed in Table II. The total cross-sectional area of each segment is the same and is based on the inner diameter of the tube, which is 38.2 mm. One geometrical detail has been left out of the model. In the experimental setup, the cold heat exchanger consists of two identical structures. However, in our case this is modeled as a single heat exchanger with twice the original length. The wall temperature distribution $T_w$ is prescribed. It is constant everywhere except in the stack, where $T_w$ is linear in $x$. For the fluid temperature, the resonator segment has an adiabatic wall boundary condition, while the hot end has an isothermal wall at $T_H$. On both sides, the volume flow is set to zero as a

result of the no-slip condition. The working gas is helium and the system is filled to an elevated pressure of 3.76 atm. Due to heating and the presence of thermoacoustic oscillations, the mean pressure increases to a steady-state value which is 0.1–0.2 atm higher.

At the segment interfaces, the geometry has discrete jumps in cross-sectional area and hydraulic radius. These jumps are softened using a smooth transition of all geometric parameters. However, the viscothermal Rott functions are not smoothed. Smoothing the geometries together has several disadvantages. First, the solution becomes dependent on the chosen transition length. Second, flow “minor losses”²⁵ due to sudden expansions and contractions cannot be modeled. An improved approach (beyond the current scope) would be to impose conservative jump conditions at the segment interfaces. This would also allow for “branching” different segments, as is required in (traveling wave) TA systems with side resonators.

The amplitude of the acoustic oscillations in the system is determined by the strength of the thermoacoustic effect, which is controlled by the wall temperature gradient in the stack. This temperature gradient is in turn determined by the temperature difference between the cold and the hot side: $\Delta T = T_H - T_0$, where $T_H$ is the wall temperature of the hot end, and $T_0$ is the ambient temperature, which is set to 293.15 K.

As a measure of grid convergence, the amplitude of the first harmonic of the pressure is used. For all shown results, this pressure amplitude does not change more than 1.5% with a doubling of the number of grid-points. For the low-amplitude cases, the number of harmonics in the simulations is set to $N_f=6$, which is sufficient to find the amplitude of the first harmonic to have negligible dependence on $N_f$. With $N_f=6$, the number of grid-points is 1059, which results in a total of 55 068 degrees of freedom for the solution.

### A. Results

Using this nonlinear model, the onset temperature difference can be determined approximately. Below onset, the numerical method does not converge because, as discussed in Appendix B, the fundamental frequency becomes undefined when the oscillation amplitude goes to zero. We use the following procedure to determine the onset temperature difference. First, a temperature difference is chosen which is higher than the onset temperature. Then, the temperature difference is lowered in fixed steps of 1 K, until the solution

<table><caption>TABLE II. Geometrical and modeling details of the standing wave engine of Atchley. The column “BL” denotes if for that segment the boundary layer approximation for the viscothermal Rott functions is used. For circular (circ.) tube geometries, $r_h = 1/2R$, where R is the tube radius. For parallel (par.) plates, $r_h = y_0$, which is half the plate spacing.</caption><thead><tr><th>Segment</th><th>Length (mm)</th><th>r<sub>h</sub> (mm)</th><th>φ</th><th>Shape</th><th>BL</th></tr></thead><tbody><tr><td>Resonator</td><td>879.7</td><td>9.55</td><td>1.00</td><td>Circ. tube</td><td>Yes</td></tr><tr><td>Cold HX</td><td>20.4</td><td>0.51</td><td>0.70</td><td>Par. plates</td><td>No</td></tr><tr><td>Stack</td><td>35.0</td><td>0.385</td><td>0.73</td><td>Par. plates</td><td>No</td></tr><tr><td>Hot HX</td><td>7.62</td><td>0.51</td><td>0.70</td><td>Par. plates</td><td>No</td></tr><tr><td>Hot end</td><td>473.0</td><td>9.55</td><td>1.00</td><td>Circ. tube</td><td>Yes</td></tr></tbody></table>
does not converge. The last value for which the solution does converge is the onset temperature.

A close agreement has been found in the onset temperature difference between our computational results and the experimental results. For the onset temperature difference we find a value of 324 K. This value is only 1 K off from the experimentally obtained value of 325 K. Close to the onset temperature, the amplitude dependence on ΔT is strong. This was also found by Karpov et al.⁹

Figure 2 shows the spatial distribution of the pressure harmonics at ΔT = 325 K. From these distributions, it is clear that the resulting time response of the pressure strongly depends on the axial position. The shapes of the harmonic profiles closely resemble the linear acoustic eigenmodes of a tube closed at both sides. The relative phases are determined by the nonlinear coupling. The resulting amplitudes of the higher harmonics are determined by the eigenfrequencies and quality factors of the linear acoustic modes. Specifically, the amount of dissonance²⁶ is of importance: the difference between the acoustic eigenfrequencies and the higher harmonics of the fundamental (driven) tone. To suppress this cascade energy transfer to higher harmonics, more dissonance could be created by shaping the cross-sectional area of the resonator.²⁷,²⁸ For x ≥ 0.9 m, the influence of the heat exchangers and stack causes these profiles to deviate from “pure” cosines.

For higher temperature differences—and hence higher drive ratios—the shape of the pressure phasors does not change significantly. However, the higher harmonics become more prominent in relation to the fundamental tone, hence nonlinear distortion increases.

Figures 3 and 4 show the acoustic pressure as a function of time on the left side of the resonator for a stack temperature difference of ΔT = 325 and 368 K. The first temperature difference is slightly above onset and the profile is nearly

sinusoidal. In the second case, the drive ratio has been increased to approximately 10%. The fundamental frequency of the ΔT = 368 K is 520 Hz, which is a bit higher than the fundamental of the ΔT = 325 K case of 516 Hz. This is due to the higher average speed of sound in the system. Figure 4 can directly be compared with Fig. 11 of Karpov et al. Some small wrinkles are found in the wave valleys; these are also found in the experimental results,²⁹ and in the results of the time domain model of Karpov et al.⁹ Increasing N<sub>f</sub> does not visibly change the wave profile. Moreover, the wrinkles in the wave valleys do not disappear.

A typical quantity of interest in a thermoacoustic system is the acoustic power flow, or the cycle-average of the work flow

$$E = \overline{pU}, \qquad (37)$$

where the overline means time-averaging. This quantity can be computed from the frequency domain data by

$$E = p_0 U_0 + \sum_{n=1}^{N_f} E_n = p_0 U_0 + \frac{1}{2} \Re \left[ \sum_{n=1}^{N_f} p_{c,n} U_{c,n}^* \right]. \quad (38)$$

In Eq. (38), E<sub>n</sub> is the acoustic power flow of the nth harmonic. Figure 5 shows the acoustic power in the first three harmonics at ΔT = 325 K. For all i, the power is zero at both ends of the tube. It can be seen that both the fundamental and second harmonic gain acoustic power in the stack. Yet the second harmonic loses nearly all its gained power in the cold heat exchanger by viscothermal losses. All harmonics higher than the second transport net power from the resonator to the stack and are therefore dissipative. For this situation, linear losses are dominant, since the power in the higher harmonics is negligible compared to the fundamental, as would be expected close to the onset temperature

<figure><img src="image_3.png" /><figcaption>FIG. 2. Spatial distribution of the real (solid line) and imaginary (dashed line) part of the acoustic pressure harmonics in the engine running at the lowest temperature difference of 325 K. The influence of the engine core can be seen at x ≈ 0.9 m, where the curvature in the profiles is relatively high.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3. Solution of the pressure waveform at $x=0$ for $\Delta T = 325$ K and $N_f = 6$.</figcaption></figure>

difference. This is no longer true in the high drive ratio/high temperature case of Fig. 6, where the higher harmonics—especially the second—have a visible effect on the fundamental.

At the high drive ratio, the heat load in the heat exchangers is relatively low, because the transverse temperature difference between solid and fluid in the heat exchangers and stack is small compared to the axial temperature difference of 368 K. Therefore, the underestimation of the time-averaged heat transfer with the model of Sec. II C does not have a significant influence on the computed saturation amplitude. For a traveling wave engine running at similar pressure amplitude, the wave would have 20 to 100 times more acoustic power, and therefore the heat load would be approximately that same factor higher. In that case, the overestimation time-average temperature difference in the heat exchangers would have a significant effect on the saturation amplitude.

For the high drive ratio, the peak RMS Reynolds number is $\sim 3 \times 10^4$. With a fundamental shear wave number of 273, this results in “weakly turbulent flow” according to Fig. 7.4 of the book by Swift.³⁰ Thus, the laminar flow

<figure><img src="image_5.png" /><figcaption>FIG. 4. Solution of the pressure waveform at $x=0$ for $\Delta T = 368$ K and $N_f = 6$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5. Acoustic power of the first three harmonics as a function of position in the tube slightly above onset ($\Delta T = 325$ K). From the left, the first vertical line is located at the left side of the cold HX. The second line is at the right side of the hot HX.</figcaption></figure>

drag model fails. In order to obtain results closer to the experimental results, a turbulent flow drag model is required, which is beyond the scope of this research.

The computed amplitude for the high temperature case is about 40% higher than the result from the time domain model of Karpov *et al.*,⁹ due to several differences between the models. First, since the wall temperature is fixed in this model, the wall temperature gradient is not decreased by thermoacoustic streaming, i.e., the wall acts as a perfect heat source/sink. Second, in the model of Karpov *et al.*, heat transfer is lowered by a factor $K_h = 25/35$, to compensate for the fact that not all plates “participate in the heat transport along the stack.” We have not included this factor in our model. Finally, their approximate drag and heat transfer

<figure><img src="image_7.png" /><figcaption>FIG. 6. Acoustic power of the first three harmonics as a function of position in the tube at $\Delta T = 368$ K). The vertical lines are the same as in Fig. 5.</figcaption></figure>
model overestimates damping for the higher frequencies, as can be seen from the drag curve (Fig. 1 in the paper by Karpov et al.⁹). The combined differences bring their model closer to the experimental results and this model closer to linear thermoacoustic theory.

## B. Time-average pressure

The system modeled in this paper is of the standing wave type. As such, the walls at the ends of the resonator do not allow a time-average mass flow ($\rho U$) through the system. Because of this, there is no Gedeon streaming.³¹ This restriction of mass streaming by the boundary conditions induces a time-average pressure distribution in the system, which also known as the acoustic radiation pressure. A correct prediction of this radiation pressure verifies the implementation and method of the current model.

Using Langevin's first relation, for a single frequency standing wave, the Eulerian excess pressure is found to be³²

$$ \overline{p_E - p_0} = \frac{|\hat{p}_{c,1}|_{x=0}^2}{4\rho_0 c_0^2} \cos(2kx) + C, \quad (39) $$

where C is a constant, that depends on external boundary conditions, $k = \omega/c_0$ is the wave number of the fundamental frequency, and $c_0$ is the undisturbed speed of sound evaluated at the ambient temperature in the resonator. To compare the numerically obtained result with the analytical result, C is set to zero and for the numerical result the spatial average over the domain $(1/L) \int_0^L \hat{p}_0 dx$ is subtracted. Figure 7 shows a comparison of the time-average pressure profile of the simulations with the analytical result of Eq. (39). For the low amplitude case, the result matches well with the analytical result in the resonator. At the position of the stack, a deviation is

<figure><img src="image_8.png" /><figcaption>FIG. 7. Analytical (solid line) vs numerical (dashed gray line) result for the time-average pressure profile in the system. Above: low amplitude case ($\Delta T = 325$ K), below: high amplitude case ($\Delta T = 368$ K). The vertical lines are the same as in Fig. 5.</figcaption></figure>

visible. This is due to the change in the cross-sectional area, as well as the non-constant time-average temperature in the stack. In case of a higher amplitude, the agreement is poorer. This is because in the high amplitude case, the influence of the higher harmonics plays a larger role.

From this we conclude that the time-averaged pressure profile shows a good agreement with the theoretical results. This gives confidence that Gedeon streaming in systems with a “looped tube” will correctly be modeled.

## V. CONCLUSIONS

The periodic steady-state of a standing wave thermoacoustic engine has been numerically solved using a quasi-one-dimensional thermoacoustic model. The periodic steady-state is found in the frequency domain by applying the nonlinear frequency domain method to this model. Several difficulties are treated related to solving self-excited thermoacoustic systems with the nonlinear frequency domain method. Among them, a method to obtain the unknown fundamental frequency of the periodic steady-state. The low-amplitude limit of this model is in exact agreement with linear thermoacoustic theory. Nonlinear effects appear in the results in the form of the generation of higher harmonics. It has been found that taking only six harmonics into account is enough to accurately describe nonlinear effects in the experimental setup and to reach global convergence. The results are in good agreement with a time domain model from literature, which shows a similar qualitative agreement with numerical experiments. Typical computation times of the developed computer code are in the order of minutes, which is small enough so that the model can be used in future optimization studies.

## ACKNOWLEDGMENTS

This research has been carried out as part of the Agentschap NL EOS-KTO (KTOT03009) research program. The financial support is gratefully acknowledged.

## APPENDIX A: MASS BALANCE

When a time domain model for a closed system is solved using a conservative integration scheme, global quantities such as mass and momentum are conserved as the solution progresses in time. If the NLFD equivalent of this time domain model is solved, this is no longer true as the solution does not follow physical states during the solution process. It results in a numerical problem, which can be explained by assuming a steady-state system. This can be done using the NLFD method by setting $N_f=0$. For $N_f=0$ all storage vectors have length one, so they contain only the time-average of the state. The model we consider is a piece of prismatic tube of length L with an isothermal wall on both sides, as shown in Fig. 8. On the left side, the temperature is kept at the reference temperature $T_0$. On the right side, the tube is heated to a temperature $T_H$. In the domain, we assume adiabatic tube walls and we neglect viscothermal effects. The problem is one-dimensional and we consider a calorically perfect gas with constant thermal conductivity. Then, the continuity equation states
<figure><img src="image_9.png" /><figcaption>FIG. 8. Example situation of a closed system with adiabatic walls on the side. At $x=0$, and $x=L$, two isothermal walls are present. A cell-centered grid is shown, where the discrete density variables lie in the middle of the cells.</figcaption></figure>

$$u \frac{d\rho}{dx} + \rho \frac{du}{dx} = 0. \qquad (A1)$$

This equation can be integrated to obtain

$$u(x) = \frac{\rho_{\text{ref}}}{\rho(x)} u_{\text{ref}}. \qquad (A2)$$

So if we set $u_{\text{ref}} = u(0) \equiv 0$, the solution we find is that the velocity is zero at all positions in the tube: the tube is quiescent. If we use this result for the momentum equation, we find

$$\frac{dp}{dx} = 0, \qquad (A3)$$

for which we find the solution $p = \text{constant}$ in the domain. Next, the energy equation reduces to

$$\frac{d^2 T}{dx^2} = 0. \qquad (A4)$$

Using the left and right boundary condition, we find the linear temperature distribution $T = T_0 + (T_H - T_0)x/L$. Finally, using the perfect gas equation of state, the density can be determined up to the same constant: $\rho = \text{constant}/R_s T$. Since we do not have a boundary condition where the pressure or the density is prescribed, the value of the pressure and the density cannot be determined. Moreover, we have too many boundary conditions for the velocity: when the velocity is zero at one boundary, it has to be zero at the other boundary due to Eq. (A2). For the numerical method, this situation results in a singular Jacobian matrix.

Nevertheless, the pressure in the system can uniquely be determined by applying knowledge of the global state of the system. For example, if we know the amount of mass in the tube, we obtain the global equation

$$\int_{V_f} \rho dV = m_p, \qquad (A5)$$

where $V_f$ is the total fluid volume in the tube and $m_p$ the amount of fluid mass in the tube. Then we can obtain the pressure as

$$p = \frac{(T_H - T_0) R_s m_p}{\ln\left(\frac{T_H}{T_0}\right) V_f}. \qquad (A6)$$

For the numerical method, we apply a discretization of the governing equations, which results in a nonlinear system of algebraic equations. Adding Eq. (A5) to this system yields

an over-determined system. Therefore, the time-averaged continuity equation for one grid-point is replaced by the equation describing the global conservation of mass. For a finite volume scheme, we replace one local DC continuity equation with Eq. (35). Since closed thermoacoustic systems are often filled at room temperature, a beneficial effect of this global equation is the fact that we obtain the DC pressure increase in the system under operating conditions. For systems with large resonators, this increase in mean pressure is small. Yet for thermoacoustic systems which have a high heated volume with respect to the total volume, such as engines with mechanical resonators, the mean pressure can increase significantly.

At this point, we are in the position to say something about the inclusion of axial conduction in Eq. (16). If we had neglected conduction, Eq. (A4) would reduce to an identity and the system would be completely under-determined. Inclusion of the conduction term solves this problem. Moreover, it enables us to resolve the thermal boundary layer adjacent to an isothermal wall.

## APPENDIX B: THE FUNDAMENTAL FREQUENCY

A TA engine is a self-excited system, therefore the fundamental frequency is unknown and part of the solution. This Appendix describes the numerical implementation used to solve for unknown fundamental frequency. A global residual operator $L$ is introduced, which combines the governing equations for all degrees of freedom. A state vector $x$ is the solution if

$$L(x) = 0. \qquad (B1)$$

This solution is searched for using the Newton-Raphson method, where we analytically evaluate the Jacobian matrix $dL/dx$.

Since for a self-excited system no boundary conditions are given for the phase, this residual operator is invariant to a shift in time. So if a global time shift is applied to one solution, this yields another valid solution. A unique solution is found by fixing the phase. Therefore, an extra equation is added to the system, following the original idea of Nakhla and Vlach,³³ as developed by Gilbert et al.³⁴ To obtain one extra equation, the phase of the solution is fixed by adding an extra equation. We call this extra equation the timing constraint. This timing constraint can be set to any dependent variable. As explained in Sec. III, we add the equation setting the imaginary part of the amplitude of the first harmonic of the pressure at the left wall equal to zero, as represented by Eq. (34). Note that the relative phase of higher harmonics with the fundamental tone is determined by the nonlinear coupling.

We define a new system of equations in which the fundamental frequency is added as unknown and the timing constraint as an equation. So the augmented solution vector is

$$y = (x, \omega), \qquad (B2)$$

and the augmented residual is
$$ M = (L(x), \Im(\hat{p}_1(x=0))). \quad (B3) $$

To search for the solution, Newton iterations are done using the Jacobian of the augmented residual operator **M**:

$$ \delta \mathbf{y} = -\frac{d\mathbf{M}^{-1}}{dy} \cdot \mathbf{M}, \quad (B4) $$

where

$$ \frac{d\mathbf{M}}{dy} = \begin{bmatrix} \frac{\partial \mathbf{L}}{\partial \mathbf{x}} & \frac{\partial \mathbf{L}}{\partial \omega} \\ \frac{\partial \Im(\hat{p}_1(x=0))}{\partial \mathbf{x}} & 0 \end{bmatrix}. \quad (B5) $$

In Eq. (B5), $\partial\Im(\hat{p}_1(x=0))/\partial\mathbf{x}$ is a single row which is zero everywhere, but is one at the global degree of freedom number corresponding to $\Im(\hat{p}_1(x=0))$. The column $\partial\mathbf{L}/\partial\omega$, is the sensitivity of the residual to a change in frequency. For brevity we only show the semi-discrete form of these sensitivities. These sensitivities can be derived from the governing equations. For the continuity equation this sensitivity is

$$ \frac{\partial \mathbf{L}_c}{\partial \omega} = S_f \tilde{\mathbf{D}} \cdot \hat{\mathbf{p}}. \quad (B6) $$

For the momentum and energy equation, however, we neglect the sensitivity to the operators $\tilde{\mathbf{D}}$, $\tilde{\mathbf{H}}$, and $\tilde{\mathbf{Q}}$, so for the momentum equation we use

$$ \frac{\partial \mathbf{L}_m}{\partial \omega} \approx \tilde{\mathbf{D}} \cdot \mathcal{F} \cdot (\rho \circ \mathbf{U}), \quad (B7) $$

and for the energy equation

$$ \frac{\partial \mathbf{L}_e}{\partial \omega} \approx \tilde{\mathbf{D}} \cdot \left( \frac{S_f}{\gamma - 1} \hat{\mathbf{p}} + \frac{1}{2S_f} \mathcal{F} \cdot (\rho \circ \mathbf{U}^2) \right), \quad (B8) $$

and finally for the equation of state

$$ \frac{\partial \mathbf{L}_s}{\partial \omega} = 0. \quad (B9) $$

In most cases it is only the fundamental mode which is unstable.³⁵ The frequency of the instability can be found using linear theory. For the application of Sec. IV, a close enough guess for the numerical solution is

$$ \omega \approx \frac{\pi c_0}{L}, \quad (B10) $$

where $c_0$ is the ambient undisturbed speed of sound and $L$ is the length of the engine.

It should be noted that for self-excited systems the quiescent situation cannot be used as an initial guess for the solver, since a quiescent solution makes the augmented Jacobian matrix singular. This is because at quiescent conditions, the solution is no longer dependent on the fundamental frequency. Therefore, as a starting guess, a “driven” solution is used. For the application of Sec. IV, the starting

guess is the solution of this system driven at the left boundary with frequency as defined in Eq. (B10) and unit amplitude.

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

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

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

<footnote n="4">⁴O. Hireche, C. Weisman, D. Baltean-Carles, P. L. Quere, and L. Bauwens, “Low Mach number analysis of idealized thermoacoustic engines with numerical solution,” J. Acoust. Soc. Am. **128**(6), 3438–3448 (2010).</footnote>

<footnote n="5">⁵A. S. Worlikar, O. M. Knio, and R. Klein, “Numerical simulation of a thermoacoustic refrigerator: II. Stratified flow around the stack,” J. Comput. Phys. **144**(2), 299–324 (1998).</footnote>

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

<footnote n="7">⁷M. Watanabe, A. Prosperetti, and H. Yuan, “A simplified model for linear and nonlinear processes in thermoacoustic prime movers. Part I. Model and linear theory,” J. Acoust. Soc. Am. **102**, 3484–3496 (1997).</footnote>

<footnote n="8">⁸H. Yuan, S. Karpov, and A. Prosperetti, “A simplified model for linear and nonlinear processes in thermoacoustic prime movers. Part II. Nonlinear oscillations,” J. Acoust. Soc. Am. **102**(6), 3497–3506 (1997).</footnote>

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

<footnote n="10">¹⁰A. I. A. El-Rahman and E. Abdel-Rahman, “Characteristic-based nonlinear simulation of large-scale standing-wave thermoacoustic engine,” J. Acoust. Soc. Am. **136**(2), 649–658 (2014).</footnote>

<footnote n="11">¹¹G. W. Swift, “Analysis and performance of a large thermoacoustic engine,” J. Acoust. Soc. Am. **92**(3), 1551–1563 (1992).</footnote>

<footnote n="12">¹²K. C. Hall, K. Ekici, J. P. Thomas, and E. H. Dowell, “Harmonic balance methods applied to computational fluid dynamics problems,” Int. J. Comput. Fluid Dyn. **27**(2), 52–67 (2013).</footnote>

<footnote n="13">¹³M. Mcmullen, A. Jameson, and J. Alonso, “Demonstration of nonlinear frequency domain methods,” Am. Inst. Aero. Astron. J. **44**(7), 1428–1435 (2006).</footnote>

<footnote n="14">¹⁴A. A. Atchley, “Study of a thermoacoustic prime mover below onset of self-oscillation,” J. Acoust. Soc. Am. **91**, 734–743 (1992).</footnote>

<footnote n="15">¹⁵Also called the “Hadamard” product, or the “Schur” product. For MATLAB users: the “*” operation.”</footnote>

<footnote n="16">¹⁶W. P. Arnott, “General formulation of thermoacoustics for stacks having arbitrarily shaped pore cross sections,” J. Acoust. Soc. Am. **90**(6), 3228–3237 (1991).</footnote>

<footnote n="17">¹⁷G. W. Swift, *Thermoacoustics: A Unifying Perspective for some Engines and Refrigerators* (Acoustical Society of America, Melville, NY, 2003), p. 88.</footnote>

<footnote n="18">¹⁸J. A. de Jong, Y. H. Wijnant, and A. de Boer, “A one-dimensional heat transfer model for parallel-plate thermoacoustic heat exchangers,” J. Acoust. Soc. Am. **135**(3), 1149–1158 (2014).</footnote>

<footnote n="19">¹⁹G. Mozurkewich, “Time-average temperature distribution in a thermoacoustic stack,” J. Acoust. Soc. Am. **103**(1), 380–388 (1998).</footnote>

<footnote n="20">²⁰I. Paek, J. E. Braun, and L. Mongeau, “Characterizing heat transfer coefficients for heat exchangers in standing wave thermoacoustic coolers,” J. Acoust. Soc. Am. **118**(4), 2271–2280 (2005).</footnote>

<footnote n="21">²¹K. W. Morton and D. F. Mayers, *Numerical Solution of Partial Differential Equations: An Introduction* (Cambridge University Press, Cambridge, UK, 2005), p. 110.</footnote>

<footnote n="22">²²B. Ward, J. Clark, and G. W. Swift, “DeltaEC users guide,” version 6.3b11, http://lanl.gov/thermoacoustics/UsersGuide.pdf (Last viewed 9/1/2015).</footnote>

<footnote n="23">²³C. Sanderson, “Armadillo: An open source C++ linear algebra library for fast prototyping and computationally intensive experiments,” technical report, NICTA, Australia (October, 2010), http://arma.sourceforge.net (last viewed 9/1/2015).</footnote>

<footnote n="24">²⁴G. Guennebaud, B. Jacob, P. Avery, A. Bachrach, and S. Barthelemy, “Eigen v3”, 2010, http://eigen.tuxfamily.org (Last viewed 1/9/2015).</footnote>

<footnote n="25">²⁵G. W. Swift, D. L. Gardner, and S. Backhaus, “Acoustic recovery of lost power in pulse tube refrigerators,” J. Acoust. Soc. Am. **105**(2), 711–724 (1999).</footnote>

<footnote n="26">²⁶Y. A. Ilinskii, B. Lipkens, T. S. Lucas, T. W. V. Doren, and E. A. Zabolotskaya, “Nonlinear standing waves in an acoustical resonator,” J. Acoust. Soc. Am. **104**(5), 2664–2674 (1998).</footnote>
²⁷D. F. Gaitan and A. A. Atchley, "Finite amplitude standing waves in harmonic and anharmonic tubes," J. Acoust. Soc. Am. **93**(5), 2489-2495 (1993).

²⁸M. F. Hamilton, Y. A. Ilinskii, and E. A. Zabolotskaya, "Linear and non-linear frequency shifts in acoustical resonators with varying cross sections," J. Acoust. Soc. Am. **110**(1), 109-119 (2001).

²⁹A. Atchley, H. Bass, and T. Hofler, "Development of nonlinear waves in a thermoacoustic prime mover," in *Frontiers of Nonlinear Acoustics; 12th International Symposium on Nonlinear Acoustics* (Elsevier, London, 1990), pp. 603-608.

³⁰G. W. Swift. *Thermoacoustics: A Unifying Perspective for some Engines and Refrigerators* (Acoustical Society of America, Melville, NY, 2003), p. 161.

³¹D. Gedeon, “DC gas flows in Stirling and pulse tube refrigerators,” *Cryocoolers* **9**, 385–392 (1997).

³²C. P. Lee and T. G. Wang, “Acoustic radiation pressure,” J. Acoust. Soc. Am. **94**(2), 1099–1109 (1993).

³³M. Nakhla and J. Vlach, “A piecewise harmonic balance technique for determination of periodic response of nonlinear systems,” *IEEE Trans. Circuits Syst.* **23**(2), 85–91 (1976).

³⁴J. Gilbert, J. Kergomard, and E. Ngoya, “Calculation of the steady-state oscillations of a clarinet using the harmonic balance technique,” J. Acoust. Soc. Am. **86**(1), 35–41 (1989).

³⁵S. Karpov and A. Prosperetti, “Nonlinear saturation of the thermoacoustic instability,” J. Acoust. Soc. Am. **107**(6), 3130–3147 (2000).