3.1 Axisymmetric two-dimensional flow assumption

The two-dimensional flow condition sets $v_{\unicode[STIX]{x1D703}}=0$ and axisymmetry sets $\unicode[STIX]{x2202}f/\unicode[STIX]{x2202}\unicode[STIX]{x1D703}=0$, where $f$ is any scalar flow quantity or component of a vector. In cylindrical coordinates, the velocity gradient tensor reduces to

(3.1)$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{v}^{\prime }=\left[\begin{array}{@{}ccc@{}}{\displaystyle \frac{\unicode[STIX]{x2202}v_{r}^{\prime }}{\unicode[STIX]{x2202}r}} & 0 & {\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}r}}\\ 0 & {\displaystyle \frac{v_{r}^{\prime }}{r}} & 0\\ {\displaystyle \frac{\unicode[STIX]{x2202}v_{r}^{\prime }}{\unicode[STIX]{x2202}z}} & 0 & {\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}z}}\end{array}\right].\end{eqnarray}$$

Note that the gradient $(\unicode[STIX]{x1D735}\boldsymbol{v}^{\prime })_{\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}}$ is not zero because of the rotation of radial unit vector in $\boldsymbol{e}_{\unicode[STIX]{x1D703}}(1/r)(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D703})(v_{r}\boldsymbol{e}_{r})$, although $\unicode[STIX]{x2202}v_{r}/\unicode[STIX]{x2202}\unicode[STIX]{x1D703}=0$. Accordingly the viscous stress tensor, applying the Stokes hypothesis, is

(3.2)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D749}^{\prime } & = & \displaystyle \unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}\boldsymbol{v}^{\prime }+(\unicode[STIX]{x1D735}\boldsymbol{v}^{\prime })^{T})-\frac{2}{3}\unicode[STIX]{x1D707}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}^{\prime })\unicode[STIX]{x1D644}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D707}\left[\begin{array}{@{}ccc@{}}{\displaystyle \frac{4}{3}}{\displaystyle \frac{\unicode[STIX]{x2202}v_{r}^{\prime }}{\unicode[STIX]{x2202}r}}-{\displaystyle \frac{2}{3}}{\displaystyle \frac{v_{r}^{\prime }}{r}}-{\displaystyle \frac{2}{3}}{\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}z}} & 0 & {\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}r}}+{\displaystyle \frac{\unicode[STIX]{x2202}v_{r}^{\prime }}{\unicode[STIX]{x2202}z}}\\ 0 & {\displaystyle \frac{4}{3}}{\displaystyle \frac{v_{r}^{\prime }}{r}}-{\displaystyle \frac{2}{3}}{\displaystyle \frac{\unicode[STIX]{x2202}v_{r}^{\prime }}{\unicode[STIX]{x2202}r}}-{\displaystyle \frac{2}{3}}{\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}z}} & 0\\ {\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}r}}+{\displaystyle \frac{\unicode[STIX]{x2202}v_{r}^{\prime }}{\unicode[STIX]{x2202}z}} & 0 & {\displaystyle \frac{4}{3}}{\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}z}}-{\displaystyle \frac{2}{3}}{\displaystyle \frac{1}{r}}{\displaystyle \frac{\unicode[STIX]{x2202}(rv_{r}^{\prime })}{\unicode[STIX]{x2202}r}}\end{array}\right].\qquad\end{eqnarray}$$

Hence the acoustic viscous dissipation function ${\mathcal{D}}_{\unicode[STIX]{x1D708}}$ identified in (2.21) takes the reduced form of

(3.3)$$\begin{eqnarray}\displaystyle {\mathcal{D}}_{\unicode[STIX]{x1D708}} & = & \displaystyle \unicode[STIX]{x1D749}^{\prime }\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{v}^{\prime }\nonumber\\ \displaystyle & = & \displaystyle 2\unicode[STIX]{x1D707}\left[\left({\displaystyle \frac{\unicode[STIX]{x2202}v_{r}^{\prime }}{\unicode[STIX]{x2202}r}}\right)^{2}+\left({\displaystyle \frac{v_{r}^{\prime }}{r}}\right)^{2}+\left({\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}z}}\right)^{2}+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}r}}+{\displaystyle \frac{\unicode[STIX]{x2202}v_{r}^{\prime }}{\unicode[STIX]{x2202}z}}\right)^{2}-\frac{1}{3}\left({\displaystyle \frac{1}{r}}\frac{\unicode[STIX]{x2202}(rv_{r}^{\prime })}{\unicode[STIX]{x2202}r}+{\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}z}}\right)^{2}\right]\!.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

This concludes all the simplification obtained with the axisymmetric cylindrical assumption.

3.2 Rott’s assumptions

Rott (Reference Rott1969) proposed three key assumptions to further simplify the linearized governing equations (2.6)–(2.8). The first two assumptions can be mathematically summarized as

(3.4)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}r}=0\end{eqnarray}$$

and

(3.5)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}T_{0}}{\unicode[STIX]{x2202}r}=0.\end{eqnarray}$$

The third one, ‘axial heat conduction in the acoustic wave and friction due to axial gradients are ignored’, is more subtle. The first part states that heat conduction happens in the radial direction only:

(3.6)$$\begin{eqnarray}q^{\prime }=k\unicode[STIX]{x1D735}T^{\prime }=k\frac{\unicode[STIX]{x2202}T^{\prime }}{\unicode[STIX]{x2202}r}\boldsymbol{e}_{r}.\end{eqnarray}$$

The second part states that the friction terms with $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}z$ are ignored, as indicated by the inclined cross line:

(3.7)

Despite its simplicity, Rott’s theory has proven to be in good agreement with experimental data when the linear approximation holds (see Backhaus & Swift Reference Backhaus and Swift2000; Wu et al. Reference Wu, Zhang, Dai and Luo2014). Therefore these underlying assumptions were used in the present study to simplify (2.6)–(2.8), yielding

(3.8)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime }}{\unicode[STIX]{x2202}t}+v_{z}^{\prime }\frac{\text{d}\unicode[STIX]{x1D70C}_{0}}{\text{d}z}+\unicode[STIX]{x1D70C}_{0}\frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}z}+\unicode[STIX]{x1D70C}_{0}{\displaystyle \frac{1}{r}}\frac{\unicode[STIX]{x2202}(rv_{r}^{\prime })}{\unicode[STIX]{x2202}r}=0, & \displaystyle\end{eqnarray}$$

(3.9)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{0}\frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}t}=-\frac{\unicode[STIX]{x2202}p^{\prime }}{\unicode[STIX]{x2202}z}+\unicode[STIX]{x1D707}\frac{\unicode[STIX]{x2202}^{2}v_{z}^{\prime }}{\unicode[STIX]{x2202}r^{2}}+\unicode[STIX]{x1D707}{\displaystyle \frac{1}{r}}{\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}r}}, & \displaystyle\end{eqnarray}$$

(3.10)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{0}T_{0}\frac{\unicode[STIX]{x2202}s^{\prime }}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D70C}_{0}T_{0}v_{z}^{\prime }\frac{\text{d}s_{0}}{\text{d}z}=k\frac{\unicode[STIX]{x2202}^{2}T^{\prime }}{\unicode[STIX]{x2202}r^{2}}+k{\displaystyle \frac{1}{r}}\frac{\unicode[STIX]{x2202}T^{\prime }}{\unicode[STIX]{x2202}r}, & \displaystyle\end{eqnarray}$$

which in this study are subject to the following boundary conditions:

(3.11a-d)$$\begin{eqnarray}T_{0}(R,z)=T_{w}(z),\quad T^{\prime }(R,z)=0,\quad \boldsymbol{v}^{\prime }(R,z)=0,\quad v_{r}^{\prime }(0,z)=0,\quad \forall z\in [0,L],\end{eqnarray}$$

where function $T_{w}(z)$ is given. The momentum equation in the radial direction is reduced to $\unicode[STIX]{x2202}p^{\prime }/\unicode[STIX]{x2202}r=0$ and mean temperature

(3.12)$$\begin{eqnarray}T_{0}(r,z)=T_{w}(z).\end{eqnarray}$$

This concludes the simplification obtained with Rott’s assumptions.

3.3 Dimensional analysis

To further simplify the disturbance energy equation (2.27) and reduce the number of independent parameters, the problem is non-dimensionalized. Scaling procedures follow the ones used by In’t Panhuis et al. (Reference In’t Panhuis, Rienstra, Molenaar and Slot2009). The dimensionless variables, denoted by an asterisk, are

(3.13)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle r=Rr^{\ast },\quad z=Lz^{\ast },\quad t=\frac{1}{\unicode[STIX]{x1D714}}t^{\ast },\quad v_{z}^{\prime }=c_{ref}v_{z}^{\prime \ast },\quad v_{r}^{\prime }=\unicode[STIX]{x1D700}c_{ref}v_{r}^{\prime \ast },\\ \displaystyle \unicode[STIX]{x1D70C}^{\prime }=\unicode[STIX]{x1D70C}_{ref}\unicode[STIX]{x1D70C}^{\prime \ast },\quad \unicode[STIX]{x1D70C}_{0}=\unicode[STIX]{x1D70C}_{ref}\unicode[STIX]{x1D70C}_{0}^{\ast },\quad p^{\prime }=\unicode[STIX]{x1D70C}_{ref}c_{ref}^{2}p^{\prime \ast },\quad p_{0}=\unicode[STIX]{x1D70C}_{ref}c_{ref}^{2}p_{0}^{\ast },\quad T^{\prime }=\frac{c_{ref}^{2}}{c_{p}}T^{\prime \ast },\\ \displaystyle T_{0}=\frac{c_{ref}^{2}}{c_{p}}T_{0}^{\ast },\quad s=c_{p}s^{\prime \ast },\quad s_{0}=c_{p}s_{0}^{\ast },\quad \unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{ref}\unicode[STIX]{x1D707}^{\ast },\quad k=k_{ref}k^{\ast },\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D700}=R/L$ is the aspect ratio and the reference point for $\unicode[STIX]{x1D70C}_{ref}$, $c_{ref}$, $\unicode[STIX]{x1D707}_{ref}$ and $k_{ref}$ was set at the mid-point ($z_{ref}=L/2$) of the regenerator. The governing equations in dimensionless form are

(3.14)$$\begin{eqnarray}\displaystyle & \displaystyle He\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{\prime \ast }}{\unicode[STIX]{x2202}t^{\ast }}+v_{z}^{\prime \ast }\frac{\text{d}\unicode[STIX]{x1D70C}_{0}^{\ast }}{\text{d}z^{\ast }}+\unicode[STIX]{x1D70C}_{0}^{\ast }\frac{\unicode[STIX]{x2202}v_{z}^{\prime \ast }}{\unicode[STIX]{x2202}z^{\ast }}+\unicode[STIX]{x1D70C}_{0}^{\ast }\frac{1}{r^{\ast }}\frac{\unicode[STIX]{x2202}(r^{\ast }v_{r}^{\prime \ast })}{\unicode[STIX]{x2202}r^{\ast }}=0, & \displaystyle\end{eqnarray}$$

(3.15)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{0}^{\ast }\frac{\unicode[STIX]{x2202}v_{z}^{\prime \ast }}{\unicode[STIX]{x2202}t^{\ast }}=-\frac{1}{He}\frac{\unicode[STIX]{x2202}p^{\prime \ast }}{\unicode[STIX]{x2202}z^{\ast }}+\frac{1}{Wo^{2}}\unicode[STIX]{x1D707}^{\ast }\left(\frac{\unicode[STIX]{x2202}^{2}v_{z}^{\prime \ast }}{\unicode[STIX]{x2202}r^{\ast 2}}+\frac{1}{r^{\ast }}\frac{\unicode[STIX]{x2202}v_{z}^{\prime \ast }}{\unicode[STIX]{x2202}r^{\ast }}\right), & \displaystyle\end{eqnarray}$$

(3.16)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{0}^{\ast }T_{0}^{\ast }\frac{\unicode[STIX]{x2202}s^{\prime \ast }}{\unicode[STIX]{x2202}t^{\ast }}+\frac{1}{He}\unicode[STIX]{x1D70C}_{0}^{\ast }T_{0}^{\ast }v_{z}^{\prime \ast }\frac{\text{d}s_{0}^{\ast }}{\text{d}z^{\ast }}=\frac{1}{Pr\cdot Wo^{2}}\frac{k\unicode[STIX]{x1D714}}{p_{0}}\left(\frac{\unicode[STIX]{x2202}^{2}T^{\prime \ast }}{\unicode[STIX]{x2202}r^{\ast 2}}+\frac{1}{r^{\ast }}\frac{\unicode[STIX]{x2202}T^{\prime \ast }}{\unicode[STIX]{x2202}r^{\ast }}\right), & \displaystyle\end{eqnarray}$$

where

(3.17a,b)$$\begin{eqnarray}He=\frac{\unicode[STIX]{x1D714}L}{c_{ref}},\quad Wo=R\left(\frac{\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{ref}}{\unicode[STIX]{x1D707}_{ref}}\right)^{1/2}\end{eqnarray}$$

are the Helmholtz number and the Womersley number, respectively, with $Pr$ being the Prandtl number.

Parameter $Wo^{2}$ is the ratio of transient inertial force ($\unicode[STIX]{x1D70C}_{ref}(c_{ref}/(1/\unicode[STIX]{x1D714}))$) to viscous force ($\unicode[STIX]{x1D707}_{ref}(c_{ref}/R)(1/R)$). In fact $Wo^{2}$ is also the Reynolds number when the length and time scales, instead of velocity scale, are imposed (which leads to the Strouhal number $St=1$), namely

(3.18)$$\begin{eqnarray}Wo^{2}=Re=\frac{\unicode[STIX]{x1D70C}_{ref}(\unicode[STIX]{x1D714}R)R}{\unicode[STIX]{x1D707}_{ref}}.\end{eqnarray}$$

Correspondingly, $Pr\cdot Wo^{2}$ would be the Péclet number, i.e.

(3.19)$$\begin{eqnarray}Pr\cdot Wo^{2}=\frac{c_{p}\unicode[STIX]{x1D707}_{ref}}{k_{ref}}\frac{\unicode[STIX]{x1D70C}_{ref}(\unicode[STIX]{x1D714}R)R}{\unicode[STIX]{x1D707}_{ref}}=Pe.\end{eqnarray}$$

Some authors (e.g. Worlikar & Knio Reference Worlikar and Knio1996) prefer the notion of $Re$ and $Pe$, but in this paper $Wo^{2}$ is used because its physical meaning can be interpreted in various ways. Firstly, $Wo^{2}$ and $Pr\cdot Wo^{2}$ can be interpreted as the ratio of two distances:

(3.20a,b)$$\begin{eqnarray}Wo=\sqrt{2}\frac{R}{\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}},\quad \sqrt{Pr}\cdot Wo=\sqrt{2}\frac{R}{\unicode[STIX]{x1D6FF}_{k}},\end{eqnarray}$$

where

(3.21a,b)$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}=\sqrt{\frac{2\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0}}},\quad \unicode[STIX]{x1D6FF}_{k}=\sqrt{\frac{2k}{\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0}c_{p}}}\end{eqnarray}$$

are the depth of the Stokes layer and its thermal counterpart, respectively. They are the characteristic length of momentum and thermal diffusion in one period. Therefore if $Wo\sim O(1)$ the whole channel is filled with oscillatory boundary layer as in the case of regenerators considered here, and $Wo\gg 1$ is found in ordinary acoustic resonators.

Another interpretation of $Wo^{2}$ is proposed by Tominaga (Reference Tominaga1995), as the ratio of two time scales:

(3.22a,b)$$\begin{eqnarray}Wo^{2}=\unicode[STIX]{x1D714}t_{\unicode[STIX]{x1D708}},\quad Pr\cdot Wo^{2}=\unicode[STIX]{x1D714}t_{k},\end{eqnarray}$$

and

(3.23a,b)$$\begin{eqnarray}t_{\unicode[STIX]{x1D708}}=\frac{R^{2}}{\unicode[STIX]{x1D708}},\quad t_{k}=\frac{R^{2}}{\unicode[STIX]{x1D6FC}},\end{eqnarray}$$

where $\unicode[STIX]{x1D708}$ is the kinematic viscosity and $\unicode[STIX]{x1D6FC}$ the thermal diffusivity, while $t_{\unicode[STIX]{x1D708}}$ and $t_{k}$ characterize the time of momentum and thermal diffusion across the channel area. Therefore $Wo^{2}$ is the time scale of diffusion measured in the period of oscillation.

3.3.1 Simplification of disturbance energy terms

In the following, the scales obtained from dimensional analysis are used to achieve the simplification of disturbance energy terms. First consider the velocity gradients:

(3.24)$$\begin{eqnarray}\left.\begin{array}{@{}cccccc@{}} & {\displaystyle \frac{\unicode[STIX]{x2202}v_{r}^{\prime }}{\unicode[STIX]{x2202}z}} & {\displaystyle \frac{\unicode[STIX]{x2202}v_{r}^{\prime }}{\unicode[STIX]{x2202}r}} & {\displaystyle \frac{v_{r}^{\prime }}{r}} & {\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}z}} & {\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}r}}\\ {\displaystyle \frac{c_{ref}}{L}}\cdot & \bigg(\unicode[STIX]{x1D700} & 1 & 1 & 1 & {\displaystyle \frac{1}{\unicode[STIX]{x1D700}}}\bigg).\end{array}\right\}\end{eqnarray}$$

The scales are written beneath as the product of a common factor (in this case $c_{ref}/L$) and relative orders defined as $\unicode[STIX]{x1D700}^{2}$, $\unicode[STIX]{x1D700}$, 1, $1/\unicode[STIX]{x1D700}$, $1/\unicode[STIX]{x1D700}^{2}$. Note that the aspect ratio $\unicode[STIX]{x1D700}$ is small, $\unicode[STIX]{x1D700}\sim O(0.01)$.

The rationale behind Rott’s third assumption, namely (3.6) and (3.7), becomes clearer with this dimensional analysis, because

(3.25)$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}^{\prime })\boldsymbol{\cdot }\boldsymbol{e}_{z}=\unicode[STIX]{x1D707}\left[{\displaystyle \frac{\unicode[STIX]{x2202}^{2}v_{z}^{\prime }}{\unicode[STIX]{x2202}r^{2}}}+{\displaystyle \frac{\unicode[STIX]{x2202}^{2}v_{r}^{\prime }}{\unicode[STIX]{x2202}r\unicode[STIX]{x2202}z}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}r}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{\unicode[STIX]{x2202}v_{r}^{\prime }}{\unicode[STIX]{x2202}z}}+{\displaystyle \frac{1}{\unicode[STIX]{x2202}z}}\left({\displaystyle \frac{4}{3}}{\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}z}}-{\displaystyle \frac{2}{3}}{\displaystyle \frac{1}{r}}{\displaystyle \frac{\unicode[STIX]{x2202}(rv_{r}^{\prime })}{\unicode[STIX]{x2202}r}}\right)\right]\!\\ \hspace{18.0pt}{\displaystyle \frac{\unicode[STIX]{x1D707}_{ref}c_{ref}}{L^{2}}}\cdot \hspace{15.60004pt}\left({\displaystyle \frac{1}{\unicode[STIX]{x1D700}^{2}}}\hspace{24.0pt}1\hspace{26.39996pt}{\displaystyle \frac{1}{\unicode[STIX]{x1D700}^{2}}}\hspace{28.79993pt}1\hspace{92.39996pt}1\hspace{40.79993pt}\right).\end{array}\right\}\end{eqnarray}$$

The leading terms $\unicode[STIX]{x2202}^{2}v_{z}^{\prime }/\unicode[STIX]{x2202}r^{2}$ and $(1/r)(\unicode[STIX]{x2202}v_{z}^{\prime }/\unicode[STIX]{x2202}r)$, which are kept in the assumption, are two orders of magnitude larger than the rest. Similarly,

(3.26)$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}\hspace{18.0pt}\boldsymbol{q}^{\prime }=k\left({\displaystyle \frac{\unicode[STIX]{x2202}T^{\prime }}{\unicode[STIX]{x2202}r}}\boldsymbol{e}_{r}+{\displaystyle \frac{\unicode[STIX]{x2202}T^{\prime }}{\unicode[STIX]{x2202}z}}\boldsymbol{e}_{z}\right)\!\\ {\displaystyle \frac{k_{ref}T_{ref}}{L}}\cdot \hspace{6.0pt}\left(~{\displaystyle \frac{1}{\unicode[STIX]{x1D700}}}\hspace{48.0pt}1\right).\end{array}\right\}\end{eqnarray}$$

The radial heat flux is one order of magnitude larger than the axial heat flux.

Now one applies the order-of-magnitude analysis to the terms in the disturbance energy equation. The scales of squared terms in viscous dissipation are

(3.27)$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}\hspace{42.0pt}\left({\displaystyle \frac{\unicode[STIX]{x2202}v_{r}^{\prime }}{\unicode[STIX]{x2202}r}}\right)^{2}+\left({\displaystyle \frac{v_{r}^{\prime }}{r}}\right)^{2}+\left({\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}z}}\right)^{2}+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}r}}+{\displaystyle \frac{\unicode[STIX]{x2202}v_{r}^{\prime }}{\unicode[STIX]{x2202}z}}\right)^{2}-{\displaystyle \frac{1}{3}}\left({\displaystyle \frac{1}{r}}{\displaystyle \frac{\unicode[STIX]{x2202}(rv_{r}^{\prime })}{\unicode[STIX]{x2202}r}}+{\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}z}}\right)^{2}\\ \left({\displaystyle \frac{c_{ref}}{L}}\right)^{2}\cdot \hspace{6.0pt}\left(1\hspace{42.0pt}1\hspace{42.0pt}1\hspace{48.0pt}{\displaystyle \frac{1}{\unicode[STIX]{x1D700}^{2}}}\hspace{20.39996pt}\unicode[STIX]{x1D700}^{2}\hspace{66.0pt}1\hspace{30.0pt}1\right).\end{array}\right\}\end{eqnarray}$$

The leading term $(\unicode[STIX]{x2202}v_{z}^{\prime }/\unicode[STIX]{x2202}r)^{2}$ is two orders of magnitude larger than the rest, so the viscous dissipation can be assumed to be

(3.28)$$\begin{eqnarray}{\mathcal{D}}_{\unicode[STIX]{x1D708}}=\unicode[STIX]{x1D707}\left({\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}r}}\right)^{2},\end{eqnarray}$$

without losing much accuracy.

For thermal dissipation, simply applying Rott’s assumption equation (3.6) yields

(3.29)$$\begin{eqnarray}{\mathcal{D}}_{k}=k\frac{\left(\unicode[STIX]{x1D735}T^{\prime }\right)^{2}}{T_{0}}=\frac{k}{T_{0}}\left(\frac{\unicode[STIX]{x2202}T^{\prime }}{\unicode[STIX]{x2202}r}\right)^{2}.\end{eqnarray}$$

For the production term, given that the direction of $\unicode[STIX]{x1D735}T_{0}$ as shown in figure 2 is orthogonal to the direction of simplified heat flux, it takes the form

(3.30)$$\begin{eqnarray}\displaystyle {\mathcal{P}} & = & \displaystyle \left(k\frac{T^{\prime }}{T_{0}}\unicode[STIX]{x1D735}T^{\prime }-T_{0}\unicode[STIX]{x1D70C}_{0}s^{\prime }\boldsymbol{v}^{\prime }\right)\boldsymbol{\cdot }\frac{\unicode[STIX]{x1D735}T_{0}}{T_{0}}\nonumber\\ \displaystyle & = & \displaystyle k\frac{T^{\prime }}{T_{0}^{2}}\frac{\unicode[STIX]{x2202}T^{\prime }}{\unicode[STIX]{x2202}r}\boldsymbol{e}_{r}\frac{\text{d}T_{0}}{\text{d}z}\boldsymbol{e}_{z}-\unicode[STIX]{x1D70C}_{0}s^{\prime }v_{z}^{\prime }\boldsymbol{e}_{z}\frac{\text{d}T_{0}}{\text{d}z}\boldsymbol{e}_{z}-\unicode[STIX]{x1D70C}_{0}s^{\prime }v_{r}^{\prime }\boldsymbol{e}_{r}\frac{\text{d}T_{0}}{\text{d}z}\boldsymbol{e}_{z}\nonumber\\ \displaystyle & = & \displaystyle -\unicode[STIX]{x1D70C}_{0}s^{\prime }v_{z}^{\prime }\frac{\text{d}T_{0}}{\text{d}z}.\end{eqnarray}$$

Finally, one considers the disturbance energy flux term by term:

(3.31)$$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }I=\unicode[STIX]{x1D735}\boldsymbol{\cdot }(p^{\prime }\boldsymbol{v}^{\prime })-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D749}^{\prime }\boldsymbol{\cdot }\boldsymbol{v}^{\prime })-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{T^{\prime }}{T_{0}}k\unicode[STIX]{x1D735}T^{\prime }\right).\end{eqnarray}$$

With $\unicode[STIX]{x1D749}^{\prime }:\unicode[STIX]{x1D735}\boldsymbol{v}^{\prime }$ and $\boldsymbol{v}^{\prime }\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}^{\prime })$ already individually simplified, given the divergence identity, it must satisfy

(3.32)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D749}^{\prime }\boldsymbol{\cdot }\boldsymbol{v}^{\prime }) & = & \displaystyle \boldsymbol{v}^{\prime }\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}^{\prime })+\unicode[STIX]{x1D749}^{\prime }\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{v}^{\prime }\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D707}\left[v_{z}\frac{\unicode[STIX]{x2202}^{2}v_{z}}{\unicode[STIX]{x2202}r^{2}}+v_{z}{\displaystyle \frac{1}{r}}{\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}r}}+\left({\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}r}}\right)^{2}\right].\end{eqnarray}$$

From the simplified form of $\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D749}^{\prime }\boldsymbol{\cdot }\boldsymbol{v}^{\prime })$, $\unicode[STIX]{x1D749}^{\prime }\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{v}^{\prime }$ and $\boldsymbol{v}^{\prime }\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}^{\prime })$, an equivalent stress tensor can be defined as

(3.33)$$\begin{eqnarray}\unicode[STIX]{x1D749}_{e}^{\prime }=\left[\begin{array}{@{}ccc@{}}0 & 0 & \unicode[STIX]{x1D70F}_{rz}\\ 0 & 0 & 0\\ 0 & 0 & 0\end{array}\right]=\unicode[STIX]{x1D707}\left[\begin{array}{@{}ccc@{}}0 & 0 & {\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}r}}\\ 0 & 0 & 0\\ 0 & 0 & 0\end{array}\right].\end{eqnarray}$$

Note that the stress tensor is said to be ‘equivalent’ because it leads to the correct simplified form (3.7), (3.28), (3.32) but it does not obey symmetry and the Stokes hypothesis. The equivalent tensor is useful in determining $\unicode[STIX]{x1D749}^{\prime }\boldsymbol{\cdot }\boldsymbol{v}^{\prime }$ in the following paragraphs.

Similarly,

(3.34)$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{T^{\prime }}{T_{0}}k\unicode[STIX]{x1D735}T^{\prime }\right) & = & \displaystyle \frac{T^{\prime }}{T_{0}}k\unicode[STIX]{x1D6FB}^{2}T^{\prime }+k\frac{(\unicode[STIX]{x1D735}T^{\prime })^{2}}{T_{0}}-k\frac{T^{\prime }}{T_{0}^{2}}\unicode[STIX]{x1D735}T_{0}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T^{\prime }\nonumber\\ \displaystyle & = & \displaystyle \frac{T^{\prime }}{T_{0}}k\left(\frac{\unicode[STIX]{x2202}^{2}T^{\prime }}{\unicode[STIX]{x2202}r^{2}}+{\displaystyle \frac{1}{r}}\frac{\unicode[STIX]{x2202}T^{\prime }}{\unicode[STIX]{x2202}r}\right)+\frac{1}{T_{0}}k\left(\frac{\unicode[STIX]{x2202}T^{\prime }}{\unicode[STIX]{x2202}r}\right)^{2}.\end{eqnarray}$$

Taking advantage of the boundary conditions, $\unicode[STIX]{x1D735}\boldsymbol{\cdot }I$ can be further simplified in the integral form. The control volume in this study was a cylinder bounded by three surfaces $\unicode[STIX]{x1D6F4}_{1}$, $\unicode[STIX]{x1D6F4}_{2}$ and $\unicode[STIX]{x1D6F4}_{3}$ as shown in figure 2. Integrating the time-averaged disturbance energy balance relation (2.27) over the control volume and applying Green’s theorem yields

(3.35)$$\begin{eqnarray}\iint _{\unicode[STIX]{x1D6F4}}\boldsymbol{n}\boldsymbol{\cdot }\overline{I}\,\text{d}S=-\iiint _{\unicode[STIX]{x1D6FA}}\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}}\,\text{d}V-\iiint _{\unicode[STIX]{x1D6FA}}\overline{{\mathcal{D}}_{k}}\,\text{d}V+\iiint _{\unicode[STIX]{x1D6FA}}\overline{{\mathcal{P}}}\,\text{d}V,\end{eqnarray}$$

where $\boldsymbol{n}$ is the outward-pointing unit normal vector on $\unicode[STIX]{x1D6F4}$. Equation (3.35) is the integral form of the time-averaged disturbance energy balance relation.

The instantaneous disturbance energy flux is

(3.36)$$\begin{eqnarray}\iint _{\unicode[STIX]{x1D6F4}}\boldsymbol{n}\boldsymbol{\cdot }I\,\text{d}S=\iint _{\unicode[STIX]{x1D6F4}}\boldsymbol{n}\boldsymbol{\cdot }\left(p^{\prime }\boldsymbol{v}^{\prime }-\unicode[STIX]{x1D749}^{\prime }\boldsymbol{\cdot }\boldsymbol{v}^{\prime }-\frac{T^{\prime }}{T_{0}}k\unicode[STIX]{x1D735}T^{\prime }\right)\,\text{d}S.\end{eqnarray}$$

On $\unicode[STIX]{x1D6F4}_{2}$, the boundary conditions (3.11) set the integration (3.36) to zero and on $\unicode[STIX]{x1D6F4}_{1}$ and $\unicode[STIX]{x1D6F4}_{3}$, $\boldsymbol{n}$ was parallel to $\boldsymbol{e}_{z}$, so (3.36) becomes

(3.37)$$\begin{eqnarray}\displaystyle \iint _{\unicode[STIX]{x1D6F4}}\boldsymbol{n}\boldsymbol{\cdot }I\,\text{d}S & = & \displaystyle \iint _{\unicode[STIX]{x1D6F4}_{3}}\boldsymbol{e}_{z}\boldsymbol{\cdot }\left(p^{\prime }\boldsymbol{v}^{\prime }-\unicode[STIX]{x1D749}^{\prime }\boldsymbol{\cdot }\boldsymbol{v}^{\prime }-\frac{T^{\prime }}{T_{0}}k\unicode[STIX]{x1D735}T^{\prime }\right)\,\text{d}S\nonumber\\ \displaystyle & & \displaystyle -\iint _{\unicode[STIX]{x1D6F4}_{1}}\boldsymbol{e}_{z}\boldsymbol{\cdot }\left(p^{\prime }\boldsymbol{v}^{\prime }-\unicode[STIX]{x1D749}^{\prime }\boldsymbol{\cdot }\boldsymbol{v}^{\prime }-\frac{T^{\prime }}{T_{0}}k\unicode[STIX]{x1D735}T^{\prime }\right)\,\text{d}S.\end{eqnarray}$$

With the equivalent stress tensor already defined, it is found that

(3.38)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D749}^{\prime }\boldsymbol{\cdot }\boldsymbol{v}^{\prime }=\unicode[STIX]{x1D749}_{e}^{\prime }\boldsymbol{\cdot }\boldsymbol{v}^{\prime }=\unicode[STIX]{x1D70F}_{rz}^{\prime }v_{z}^{\prime }\boldsymbol{e}_{r}=\unicode[STIX]{x1D707}{\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}r}}v_{z}^{\prime }\boldsymbol{e}_{r}, & \displaystyle\end{eqnarray}$$

(3.39)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{T^{\prime }}{T_{0}}k\unicode[STIX]{x1D735}T^{\prime }=\frac{T^{\prime }}{T_{0}}k\frac{\unicode[STIX]{x2202}T^{\prime }}{\unicode[STIX]{x2202}r}\boldsymbol{e}_{r}, & \displaystyle\end{eqnarray}$$

which are both orthogonal to $\boldsymbol{e}_{z}$.

Hence the final simplified expression of $\iint _{\unicode[STIX]{x1D6F4}}\boldsymbol{n}\boldsymbol{\cdot }I\,\text{d}S$ in a cylindrical domain is

(3.40)$$\begin{eqnarray}\iint _{\unicode[STIX]{x1D6F4}}\boldsymbol{n}\boldsymbol{\cdot }I\,\text{d}S=\iint _{\unicode[STIX]{x1D6F4}_{3}}p^{\prime }v_{z}^{\prime }\,\text{d}S-\iint _{\unicode[STIX]{x1D6F4}_{1}}p^{\prime }v_{z}^{\prime }\,\text{d}S.\end{eqnarray}$$

Therefore it is shown that, after the simplification, the disturbance energy flux in a cylindrical domain reduces to the classic acoustic flux.

This section has reported in detail the simplification of ${\mathcal{D}}_{\unicode[STIX]{x1D708}}$, ${\mathcal{D}}_{k}$, ${\mathcal{P}}$ and $\iint \boldsymbol{n}\boldsymbol{\cdot }I\,\text{d}S$, taking the form of (3.28), (3.29), (3.30) and (3.40), respectively.

3.3.2 Dimensionless disturbance energy balance

The integral form of the disturbance energy balance is non-dimensionalized using integrated inlet disturbance energy flux $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{\unicode[STIX]{x1D6F4}_{1}}$. Shorthand notations $(\cdot )_{cv}$ and $(\cdot )_{cs}$ are used for $\iiint _{\unicode[STIX]{x1D6FA}}\cdot \text{d}V$ and $\iint _{\unicode[STIX]{x1D6F4}}\cdot \text{d}S$, respectively. The equation takes the form

(3.41)$$\begin{eqnarray}(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }=-(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }-(\overline{{\mathcal{D}}_{k}})_{cv}^{\ast }+(\overline{{\mathcal{P}}})_{cv}^{\ast },\end{eqnarray}$$

where each term

(3.42)$$\begin{eqnarray}(\cdot )_{cs|cv}^{\ast }=\frac{(\cdot )_{cs|cv}}{(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{\unicode[STIX]{x1D6F4}_{1}}}.\end{eqnarray}$$

The term $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$ then represents the relative wave growth; $(\overline{{\mathcal{P}}})^{\ast }$, $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})^{\ast }$ and $(\overline{{\mathcal{D}}_{k}})^{\ast }$ are production, viscous dissipation and thermal dissipation relative to the energy transported by the wave, respectively. The benefit of scaling in this way is to offer better intuition since terms are expressed as fractions or multiples of inlet energy and are hence $O(1)$. Also the equality of left- and right-hand sides of the equation can be easily shown by taking the difference $\unicode[STIX]{x1D6E5}^{\ast }=(\overline{{\mathcal{P}}})_{cv}^{\ast }-(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }-(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }-(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$.

It is worth noting that another possible scaling procedures is to non-dimensionalize the energy terms by the respective flow variables. For example,

(3.43)$$\begin{eqnarray}(\overline{{\mathcal{D}}_{k}})_{cv}=\int _{0}^{L}\int _{0}^{R}\frac{k}{T_{0}}\left(\frac{\unicode[STIX]{x2202}T^{\prime }}{\unicode[STIX]{x2202}r}\right)^{2}r\,\text{d}r\,\text{d}l=k_{ref}T_{ref}L\int _{0}^{1}\int _{0}^{1}\frac{k^{\ast }}{T_{0}^{\ast }}\left(\frac{\unicode[STIX]{x2202}T^{\prime \ast }}{\unicode[STIX]{x2202}r^{\ast }}\right)^{2}r^{\ast }\,\text{d}r^{\ast }\,\text{d}l^{\ast }.\end{eqnarray}$$

Define

(3.44a,b)$$\begin{eqnarray}(\overline{{\mathcal{D}}_{k}})^{\dagger }=\frac{k^{\ast }}{T_{0}^{\ast }}\left(\frac{\unicode[STIX]{x2202}T^{\prime \ast }}{\unicode[STIX]{x2202}r^{\ast }}\right)^{2},\quad (\overline{{\mathcal{D}}_{k}})_{cv}^{\dagger }=\int _{0}^{1}\int _{0}^{1}(\overline{{\mathcal{D}}_{k}})^{\dagger }r^{\ast }\,\text{d}r^{\ast }\,\text{d}l^{\ast },\end{eqnarray}$$

where the dagger ($\dagger$) superscript denotes non-dimensionalization using corresponding flow variables. It can be shown that the two non-dimensionalization methods are equivalent, because

(3.45)$$\begin{eqnarray}(\overline{{\mathcal{D}}_{k}})_{cv}^{\ast }=\frac{(\overline{{\mathcal{D}}_{k}})_{cv}}{(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{\unicode[STIX]{x1D6F4}_{1}}}=\frac{k_{ref}{\displaystyle \frac{c_{ref}^{2}}{c_{p}}}L(\overline{{\mathcal{D}}_{k}})^{\dagger }}{\unicode[STIX]{x1D70C}_{ref}c_{ref}^{3}R^{2}(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{\unicode[STIX]{x1D6F4}_{1}}^{\dagger }}=\frac{He}{Wo^{2}}\frac{(\overline{{\mathcal{D}}_{k}})^{\dagger }}{(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{\unicode[STIX]{x1D6F4}_{1}}^{\dagger }}.\end{eqnarray}$$

For two dynamically similar flows, $He$, $Wo^{2}$, $(\overline{{\mathcal{D}}_{k}})^{\dagger }$ and $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{\unicode[STIX]{x1D6F4}_{1}}^{\dagger }$ are identical and hence so are $(\overline{{\mathcal{D}}_{k}})_{cv}^{\ast }$. Similarly,

(3.46a,b)$$\begin{eqnarray}(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})^{\dagger }=\unicode[STIX]{x1D707}^{\ast }\left({\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime \ast }}{\unicode[STIX]{x2202}r^{\ast }}}\right)^{2},\quad (\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\dagger }=\int _{0}^{1}\int _{0}^{1}(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})^{\dagger }r^{\ast }\,\text{d}r^{\ast }\,\text{d}l^{\ast }.\end{eqnarray}$$

(3.47a,b)$$\begin{eqnarray}(\overline{{\mathcal{P}}})^{\dagger }=-\unicode[STIX]{x1D70C}_{0}^{\ast }s^{\prime \ast }v_{z}^{\prime \ast }\frac{\text{d}T_{0}^{\ast }}{\text{d}z^{\ast }},\quad (\overline{{\mathcal{P}}})_{cv}^{\dagger }=\int _{0}^{1}\int _{0}^{1}(\overline{{\mathcal{P}}})^{\dagger }r^{\ast }\,\text{d}r^{\ast }\,\text{d}l^{\ast }.\end{eqnarray}$$

The equivalence of these two methods allows us to use the flow variables and (3.44), (3.46) and (3.47) to explain the trends of disturbance energy terms defined in (3.41).

The differential disturbance energy balance is non-dimensionalized using inlet flux density $(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{I})_{\unicode[STIX]{x1D6F4}_{1}}$ on the axis $r=0$:

(3.48)$$\begin{eqnarray}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{I})^{\ast }=-(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})^{\ast }-(\overline{{\mathcal{D}}_{k}})^{\ast }+(\overline{{\mathcal{P}}})^{\ast },\end{eqnarray}$$

where each term

(3.49)$$\begin{eqnarray}(\cdot )^{\ast }=\frac{(\cdot )}{(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{I})_{\unicode[STIX]{x1D6F4}_{1},r=0}}.\end{eqnarray}$$

This section has shown the non-dimensionalization of the governing equations. It has also discussed key dimensionless groups: $He$ and $Wo$, the simplification of disturbance energy terms and the non-dimensionalization of the disturbance energy balance.

3.4 Solution

Assuming that all the disturbances are time-harmonic signals with angular frequency $\unicode[STIX]{x1D714}$, equations (3.8)–(3.10) can be solved in the frequency domain $f^{\prime }(\boldsymbol{x},t)=\text{Re}(\hat{f}(\boldsymbol{x})\text{e}^{\text{i}\unicode[STIX]{x1D714}t})$, where $\hat{f}(\boldsymbol{x})$ is the complex amplitude. Since instantaneous values are of no interest in this study, no distinction is made between frequency-domain and time-domain variables. Complex and real disturbances are uniformly denoted as $f^{\prime }$, whose meaning should be clear in the context. The equations were solved in two dimensions instead of using area average performed by Rott in the original paper and in most subsequent work.

Note that the solutions are given in dimensional form on purpose because this means that it is easier for readers to compare the solutions reported here with those in the previous literature. Dimensional form is also easier to interpret in practice for implementation. Dimensionless solutions are given on an ad hoc basis when they are needed in the results section.

Velocity and temperature fluctuations were solved in terms of $p^{\prime }$ and $\text{d}p^{\prime }/\text{d}z$, with $\text{i}$ being the complex unit:

(3.50)$$\begin{eqnarray}\displaystyle & \displaystyle v_{z}^{\prime }=\frac{1-{\mathcal{H}}_{\unicode[STIX]{x1D708}}}{-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0}}\frac{\text{d}p^{\prime }}{\text{d}z}, & \displaystyle\end{eqnarray}$$

(3.51)$$\begin{eqnarray}\displaystyle & \displaystyle v_{r}^{\prime }=\frac{1}{2}r\left[-\frac{\text{i}\unicode[STIX]{x1D714}}{p_{0}}\frac{1+(\unicode[STIX]{x1D6FE}-1){\mathcal{G}}_{k}}{\unicode[STIX]{x1D6FE}}p^{\prime }+\frac{\text{i}R_{g}}{\unicode[STIX]{x1D714}p_{0}}\frac{\text{d}T_{w}}{\text{d}z}\frac{{\mathcal{G}}_{k}-{\mathcal{G}}_{\unicode[STIX]{x1D708}}}{1-Pr}\frac{\text{d}p^{\prime }}{\text{d}z}-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\left(\frac{1-{\mathcal{G}}_{\unicode[STIX]{x1D708}}}{-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0}}\frac{\text{d}p^{\prime }}{\text{d}z}\right)\right],\qquad & \displaystyle\end{eqnarray}$$

(3.52)$$\begin{eqnarray}\displaystyle & \displaystyle T^{\prime }=\frac{1-{\mathcal{H}}_{k}}{\unicode[STIX]{x1D70C}_{0}c_{p}}p^{\prime }-\frac{1}{\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D714}^{2}}\frac{\text{d}T_{w}}{\text{d}z}\frac{(1-{\mathcal{H}}_{k})-Pr(1-{\mathcal{H}}_{\unicode[STIX]{x1D708}})}{1-Pr}\frac{\text{d}p^{\prime }}{\text{d}z}, & \displaystyle\end{eqnarray}$$

(3.53)$$\begin{eqnarray}\displaystyle & \displaystyle s^{\prime }=\frac{-{\mathcal{H}}_{k}}{\unicode[STIX]{x1D70C}_{0}T_{0}}p^{\prime }-\frac{c_{p}}{\unicode[STIX]{x1D70C}_{0}T_{0}\unicode[STIX]{x1D714}^{2}}\frac{\text{d}T_{w}}{\text{d}z}\frac{(1-{\mathcal{H}}_{k})-Pr(1-{\mathcal{H}}_{\unicode[STIX]{x1D708}})}{1-Pr}\frac{\text{d}p^{\prime }}{\text{d}z}, & \displaystyle\end{eqnarray}$$

where the radial distribution is given by

(3.54)$$\begin{eqnarray}{\mathcal{H}}(r,z)=\frac{\text{J}_{0}[(\text{i}-1)r/\unicode[STIX]{x1D6FF}]}{\text{J}_{0}[(\text{i}-1)R/\unicode[STIX]{x1D6FF}]}.\end{eqnarray}$$

Here $\text{J}_{0}$ and $\text{J}_{1}$ are Bessel functions of the first kind of order 0 and 1, respectively; ${\mathcal{H}}_{\unicode[STIX]{x1D708}}$ uses viscous boundary layer thickness $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$, and ${\mathcal{H}}_{k}$ uses thermal boundary layer thickness $\unicode[STIX]{x1D6FF}_{k}$, as defined in (3.21). The solution structure was adapted from Swift (Reference Swift2002). Notice that ${\mathcal{H}}$ contains $R/\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ and $R/\unicode[STIX]{x1D6FF}_{k}$ which are the local $(1/\sqrt{2})Wo$ and $(1/\sqrt{2})\sqrt{Pr}Wo$, underlying again the importance of the Womersley number. The radial distribution function ${\mathcal{G}}$ is defined as the area average of ${\mathcal{H}}$ in a circle of radius $r$:

(3.55)$$\begin{eqnarray}{\mathcal{G}}(r,z)=\frac{2}{r^{2}}\int _{0}^{r}{\mathcal{H}}x\,\text{d}x=\frac{2\unicode[STIX]{x1D6FF}}{(\text{i}-1)r}\frac{\text{J}_{1}[(\text{i}-1)r/\unicode[STIX]{x1D6FF}]}{\text{J}_{0}[(\text{i}-1)R/\unicode[STIX]{x1D6FF}]},\end{eqnarray}$$

which also has two forms ${\mathcal{G}}_{\unicode[STIX]{x1D708}}$ and ${\mathcal{G}}_{k}$. In the integral $x$ is a dummy variable. At $r=0$,

(3.56)$$\begin{eqnarray}{\mathcal{G}}(0,z):=\lim _{r\rightarrow 0^{+}}{\mathcal{G}}(r,z)=1/\text{J}_{0}[(\text{i}-1)R/\unicode[STIX]{x1D6FF}]={\mathcal{H}}(0,z).\end{eqnarray}$$

The derivatives needed to evaluate the simplified ${\mathcal{P}}$, ${\mathcal{D}}_{\unicode[STIX]{x1D708}}$, ${\mathcal{D}}_{k}$ were also derived:

(3.57)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}r}=\frac{{\mathcal{I}}_{\unicode[STIX]{x1D710}}}{-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0}}\frac{\text{d}p^{\prime }}{\text{d}z}, & \displaystyle\end{eqnarray}$$

(3.58)$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}v_{z}^{\prime }}{\unicode[STIX]{x2202}z}}=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\left(\frac{1-{\mathcal{H}}_{\unicode[STIX]{x1D710}}}{-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0}}\frac{\text{d}p^{\prime }}{\text{d}z}\right), & \displaystyle\end{eqnarray}$$

(3.59)$$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}v_{r}^{\prime }}{\unicode[STIX]{x2202}r}} & = & \displaystyle -\frac{\text{i}\unicode[STIX]{x1D714}}{2p_{0}\unicode[STIX]{x1D6FE}}p^{\prime }-\frac{\text{i}}{2\unicode[STIX]{x1D714}}\frac{\text{d}}{\text{d}z}\left(\frac{1}{\unicode[STIX]{x1D70C}_{0}}\frac{\text{d}p^{\prime }}{\text{d}z}\right)-\frac{\text{i}\unicode[STIX]{x1D714}}{p_{0}}\frac{(\unicode[STIX]{x1D6FE}-1)\left({\mathcal{H}}_{k}-\frac{1}{2}{\mathcal{G}}_{k}\right)}{\unicode[STIX]{x1D6FE}}p^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\text{i}R_{g}}{\unicode[STIX]{x1D714}p_{0}}\frac{\text{d}T_{w}}{\text{d}z}\frac{({\mathcal{H}}_{k}-\frac{1}{2}{\mathcal{G}}_{k})-({\mathcal{H}}_{\unicode[STIX]{x1D708}}-\frac{1}{2}{\mathcal{G}}_{\unicode[STIX]{x1D708}})}{1-Pr}\frac{\text{d}p^{\prime }}{\text{d}z}-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}z}\left(\frac{{\mathcal{H}}_{\unicode[STIX]{x1D708}}-\frac{1}{2}{\mathcal{G}}_{\unicode[STIX]{x1D708}}}{-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0}}\frac{\text{d}p^{\prime }}{\text{d}z}\right),\qquad\end{eqnarray}$$

(3.60)$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}T^{\prime }}{\unicode[STIX]{x2202}r}=\frac{1}{\unicode[STIX]{x1D70C}_{0}c_{p}}{\mathcal{I}}_{k}p^{\prime }-\frac{1}{\unicode[STIX]{x1D70C}_{0}\unicode[STIX]{x1D714}^{2}}\frac{\text{d}T_{w}}{\text{d}z}\frac{{\mathcal{I}}_{k}-Pr{\mathcal{I}}_{\unicode[STIX]{x1D708}}}{1-Pr}\frac{\text{d}p^{\prime }}{\text{d}z},\end{eqnarray}$$

with

(3.61)$$\begin{eqnarray}{\mathcal{I}}(r,z)=-\frac{\unicode[STIX]{x2202}{\mathcal{H}}}{\unicode[STIX]{x2202}r}=\frac{\text{i}-1}{\unicode[STIX]{x1D6FF}}\frac{\text{J}_{1}[(\text{i}-1)r/\unicode[STIX]{x1D6FF}]}{\text{J}_{0}[(\text{i}-1)R/\unicode[STIX]{x1D6FF}]}.\end{eqnarray}$$

This concludes the definition of radial distribution functions ${\mathcal{H}}$, ${\mathcal{G}}$, ${\mathcal{I}}$ and the analytical solutions in the radial direction in terms of $p^{\prime }$ and $\text{d}p^{\prime }/\text{d}z$.

3.5 Numerical scheme

Since all the flow quantities were expressed in terms of $p^{\prime }$ and $\text{d}p^{\prime }/\text{d}z$, it is important to find an efficient numerical scheme to solve them in the axial direction. Now let

(3.62)$$\begin{eqnarray}{\mathcal{F}}(z)={\mathcal{G}}(R,z).\end{eqnarray}$$

Then evaluating (3.51) at $R$ and applying the non-slip boundary condition $v_{r}(R,z)=0$, yields

(3.63)$$\begin{eqnarray}\displaystyle & & \displaystyle -\frac{\text{i}\unicode[STIX]{x1D714}}{p_{0}}\frac{1+(\unicode[STIX]{x1D6FE}-1){\mathcal{F}}_{k}}{\unicode[STIX]{x1D6FE}}p^{\prime }+\frac{1}{T_{w}}\frac{\text{d}T_{w}}{\text{d}z}\frac{{\mathcal{F}}_{k}-{\mathcal{F}}_{\unicode[STIX]{x1D708}}}{(1-{\mathcal{F}}_{\unicode[STIX]{x1D708}})(1-Pr)}\left(\frac{1-{\mathcal{F}}_{\unicode[STIX]{x1D708}}}{-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0}}\frac{\text{d}p^{\prime }}{\text{d}z}\right)\nonumber\\ \displaystyle & & \displaystyle \quad -\,\frac{\text{d}}{\text{d}z}\left(\frac{1-{\mathcal{F}}_{\unicode[STIX]{x1D708}}}{-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0}}\frac{\text{d}p^{\prime }}{\text{d}z}\right)=0,\end{eqnarray}$$

which is Rott’s thermoacoustic wave equation. Function ${\mathcal{F}}(z)$ is the area average of ${\mathcal{H}}$ on the entire circular cross-section and hence (3.63) only depends on $z$. The physical meaning of the term $((1-{\mathcal{F}}_{\unicode[STIX]{x1D708}})/-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0})(\text{d}p^{\prime }/\text{d}z)$ is area-averaged axial velocity $(v_{z}^{\prime })_{av}$ because

(3.64)$$\begin{eqnarray}(v_{z}^{\prime })_{av}=\frac{1}{\unicode[STIX]{x03C0}R^{2}}\int _{0}^{R}\frac{1-{\mathcal{H}}_{\unicode[STIX]{x1D708}}}{-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0}}\frac{\text{d}p^{\prime }}{\text{d}z}2\unicode[STIX]{x03C0}r\,\text{d}r=\frac{1-{\mathcal{G}}_{\unicode[STIX]{x1D708}}(R,z)}{-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0}}\frac{\text{d}p^{\prime }}{\text{d}z}=\frac{1-{\mathcal{F}}_{\unicode[STIX]{x1D708}}(z)}{-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0}}\frac{\text{d}p^{\prime }}{\text{d}z}.\end{eqnarray}$$

Therefore the second-order ordinary differential equation (3.63) can be conveniently converted to two simultaneous first-order ordinary differential equations:

(3.65a)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}p^{\prime }}{\text{d}z}=\frac{-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0}}{1-{\mathcal{F}}_{\unicode[STIX]{x1D708}}(z)}(v_{z}^{\prime })_{av}, & \displaystyle\end{eqnarray}$$

(3.65b)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}(v_{z}^{\prime })_{av}}{\text{d}z}=-\frac{\text{i}\unicode[STIX]{x1D714}}{p_{0}}\frac{1+(\unicode[STIX]{x1D6FE}-1){\mathcal{F}}_{k}}{\unicode[STIX]{x1D6FE}}p^{\prime }+\frac{1}{T_{w}}\frac{\text{d}T_{w}}{\text{d}z}\frac{{\mathcal{F}}_{k}-{\mathcal{F}}_{\unicode[STIX]{x1D708}}}{(1-{\mathcal{F}}_{\unicode[STIX]{x1D708}})(1-Pr)}(v_{z}^{\prime })_{av}, & \displaystyle\end{eqnarray}$$

$(v_{z}^{\prime })_{av}$

$((1-{\mathcal{F}}_{\unicode[STIX]{x1D708}})/-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0})(\text{d}p^{\prime }/\text{d}z)$

$R$

$\unicode[STIX]{x1D714}$

$p_{0}$

$p^{\prime }$

$\text{d}p^{\prime }/\text{d}z$

$(v_{z}^{\prime })_{av}$

$p^{\prime }$

$\unicode[STIX]{x1D6F4}_{1}$

$10^{-6}$

$p^{\prime }$

$\text{d}p^{\prime }/\text{d}z$

$T^{\prime }$

$v_{z}^{\prime }$

${\mathcal{H}}$

${\mathcal{G}}$

$(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}z)(((1-{\mathcal{G}}_{\unicode[STIX]{x1D708}})/-\text{i}\unicode[STIX]{x1D714}\unicode[STIX]{x1D70C}_{0})(\text{d}p^{\prime }/\text{d}z))$

$v_{r}^{\prime }$

$z$

$\unicode[STIX]{x1D6FF}$

${\mathcal{H}}$

${\mathcal{G}}$

Figure 3.

Illustration of computational grid in this study, including main output points (○), auxiliary points (♢) and interim points ($\times$). Here $h_{z}$ and $h_{r}$ are the axial and radial output step sizes, respectively, and $h$ is the self-adaptive integration step size.

A global error is needed to validate the simplification made and ensure overall accuracy of the numerical procedures. In the linear regime, equation (3.35) holds exactly. The equality is compromised by simplification and use of numerical methods. Therefore in this study the global relative error was defined as

(3.66)$$\begin{eqnarray}\unicode[STIX]{x1D716}=\frac{|(\overline{{\mathcal{P}}})_{cv}-(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}-(\overline{{\mathcal{D}}_{k}})_{cv}-(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}|}{\max (|(\overline{{\mathcal{P}}})_{cv}|,|(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}|,|(\overline{{\mathcal{D}}_{k}})_{cv}|,|(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}|)}.\end{eqnarray}$$

The mesh numbers were adjusted during the computation to reach a satisfactory error level (here of the order of $10^{-5}$) to balance computation speed and accuracy.

4.1 Parametric study: Womersley number

The Womersley number is the most widely discussed dimensionless number in the research of thermoacoustic devices. This is reasonable since it appears in both the dimensionless momentum equation as $Wo^{2}$ and the energy equation as $Pr\cdot Wo^{2}$. Although $\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}(T_{0})$ and $k=k(T_{0})$, the Prandtl number is very weakly dependent on temperature, varying only by 0.01 in the temperature range [325, 525] K. Therefore the Prandtl number is treated as a constant for helium in this study. Parameter $Wo^{2}$ being a function of radius, density, frequency and viscosity, can be varied in many ways. In this study, it was varied by changing $p_{0}$ and, ultimately, $\unicode[STIX]{x1D70C}_{ref}$, given that $\unicode[STIX]{x1D735}p_{0}=0$ and $p_{0}=\unicode[STIX]{x1D70C}_{ref}R_{g}T_{ref}$. Note that $|p^{\prime }|_{\unicode[STIX]{x1D6F4}_{1}}$ needed to be changed to keep $p^{\prime \ast }$ constant. It follows that the inlet disturbance energy flux $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{\unicode[STIX]{x1D6F4}_{1}}$ used to scale the disturbance energy terms was a linear function of $Wo^{2}$.

The terms in the integral disturbance energy balance (3.41) are shown in figure 4 against $Wo^{2}$. The vertical heights of light, medium and dark grey areas represent, respectively, dimensionless values of the disturbance energy flux $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$, thermal dissipation $(\overline{{\mathcal{D}}_{k}})_{cv}^{\ast }$ and viscous dissipation $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }$, of which the sum should equal the thermoacoustic production $(\overline{{\mathcal{P}}})_{cv}^{\ast }$ (dashed line in figure 4a). The deviation from the equality, defined as $\unicode[STIX]{x1D6E5}^{\ast }=(\overline{{\mathcal{P}}})_{cv}^{\ast }-(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }-(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }-(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$, is plotted in figure 4(b). It can be seen that the difference is four orders of magnitude smaller and hence validates the simplification and assumption made previously. In figure 5, $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }$, $(\overline{{\mathcal{D}}_{k}})_{cv}^{\ast }$ and total dissipation $(\overline{{\mathcal{D}}}_{k+\unicode[STIX]{x1D708}})_{cv}^{\ast }$ are shown again for better visualization.

Figure 4.

The balance of the dimensionless disturbance energy terms against $Wo^{2}$. (a) The values of the dimensionless thermoacoustic production in the control volume $(\overline{{\mathcal{P}}})_{cv}^{\ast }$ (– – –), viscous dissipation $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }$ (▪), thermal dissipation $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }$ (

) and disturbance energy flux $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$ (

). (b) The deviation from equality defined as $\unicode[STIX]{x1D6E5}^{\ast }=(\overline{{\mathcal{P}}})_{cv}^{\ast }-(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }-(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }-(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$.

Figure 5.

The changes of (a) the dimensionless thermal dissipation $(\overline{{\mathcal{D}}_{k}})_{cv}^{\ast }$ (——), viscous dissipation $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }$ (– – –) and total dissipation $(\overline{{\mathcal{D}}}_{k+\unicode[STIX]{x1D708}})_{cv}^{\ast }$ (— ⋅ —) and (b) the viscous dissipation (▪) and thermal dissipation (

) as a percentage of the total dissipation against $Wo^{2}$.

A straightforward yet important conclusion from figure 5(b) is that the relative importance of viscous and thermal dissipation, which can be defined as the percentages $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}/(\overline{{\mathcal{D}}}_{k+\unicode[STIX]{x1D708}})_{cv}$ and $(\overline{{\mathcal{D}}_{k}})_{cv}/(\overline{{\mathcal{D}}}_{k+\unicode[STIX]{x1D708}})_{cv}$, respectively, shifts with $Wo^{2}$. When other conditions are maintained, viscous dissipation is more significant at low $Wo^{2}$ while thermal dissipation is the main loss mechanism at high $Wo^{2}$. This is because $Wo$, according to (3.20), represents the channel radius relative to boundary-layer thickness; as the relative radius becomes larger, the portion of main flow, which has little velocity gradient, also grows. Hence the relative viscous dissipation decreases. For the temperature profile, no oscillation is allowed (i.e. $T^{\prime }=0$) at the wall due to the implicit assumption of its infinite heat capacity originated from the boundary condition $T(R,z)=T_{w}(z)$. As the relative radius becomes larger compared to the thermal boundary layer, the wall’s ability to ‘smooth out’ the temperature profile decreases as shown in figure 6: the magnitude of temperature oscillation becomes larger and so does its radial gradient. This leads to a surge of thermal dissipation. In summary, figure 5(b) points out the main loss mechanism to tackle when building thermoacoustic devices.

Figure 6.

The changes of the magnitude of (a) the dimensionless temperature fluctuation $|T^{\prime \ast }|$ and (b) the dimensionless radial derivative of temperature fluctuation $|\unicode[STIX]{x2202}T^{\prime \ast }/\unicode[STIX]{x2202}r^{\ast }|$ against dimensionless radial position $r/R$ at mid-axial position $z=L/2$ for four different $Wo^{2}$: $Wo^{2}=0.06$ (— ⋅ —), $Wo^{2}=0.19$ (——), $Wo^{2}=0.32$ (– – –), $Wo^{2}=1.9$ ($\cdots \cdots$).

The changes of $(\overline{{\mathcal{P}}})_{cv}^{\ast }$, $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }$, $(\overline{{\mathcal{D}}_{k}})_{cv}^{\ast }$ and $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$ with $Wo^{2}$ can be divided into four regions by three critical values of $Wo^{2}$ as labelled in figure 4. With $Wo^{2}$ increasing from 0, $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$ grows faster than $(\overline{{\mathcal{P}}})_{cv}^{\ast }$. Given $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }=(\overline{{\mathcal{P}}})_{cv}^{\ast }-(\overline{{\mathcal{D}}}_{k+\unicode[STIX]{x1D708}})_{cv}^{\ast }$, this faster increase of acoustic output clearly benefits from the initial drop of total dissipation as shown in figure 5(a) from 0 to $C_{1}$; $C_{1}$ is $Wo^{2}$ where total dissipation reaches a minimum. It follows that $C_{1}$ is also the position where $(\overline{{\mathcal{P}}})_{cv}^{\ast }$ and $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$ have parallel tangent lines, i.e. $\text{d}(\overline{{\mathcal{P}}})_{cv}^{\ast }=\text{d}(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$. Formally, $C_{1}$ is defined as

(4.3)$$\begin{eqnarray}C_{1}=\min (Wo^{2}:\text{d}(\overline{{\mathcal{D}}}_{k+\unicode[STIX]{x1D708}})_{cv}^{\ast }=0),\end{eqnarray}$$

where the $\min (\,)$ function ensures the first such $Wo^{2}$. The physical meaning of $C_{1}$ is $Wo^{2}$ where thermoacoustic production starts to grow faster than acoustic output and where the system has a minimum loss.

On continuing to increase $Wo^{2}$, the total dissipation starts to pick up and the rate of change of the production exceeds that of the disturbance energy flux and the second critical value $C_{2}$ can be defined as

(4.4)$$\begin{eqnarray}C_{2}=\min \left(Wo^{2}:\text{d}\frac{(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }}{(\overline{{\mathcal{P}}})_{cv}^{\ast }}=0\right).\end{eqnarray}$$

Equivalently, $C_{2}$ is also where $\text{d}\text{ln}(\overline{{\mathcal{P}}})_{cv}^{\ast }=\text{d}\text{ln}(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$. It is clear that this ratio,

(4.5)$$\begin{eqnarray}\unicode[STIX]{x1D702}:=\frac{(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }}{(\overline{{\mathcal{P}}})_{cv}^{\ast }},\end{eqnarray}$$

is the efficiency of the conversion from thermal energy to acoustic energy. Then $C_{2}$ is $Wo^{2}$ where maximum efficiency is reached. This new definition of ‘true’ thermoacoustic efficiency reflects only the energy conversion between different aspects of the disturbance and does not concern any steady-state physical process. The changes of thermoacoustic efficiency are better visualized in figure 7.

Figure 7.

The changes of the thermoacoustic efficiency $\unicode[STIX]{x1D702}$ against $Wo^{2}$.

On further increasing $Wo^{2}$, it can be seen that both thermoacoustic production and dissipation (mainly thermal) rise sharply. It is clear in figure 4(a) that $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$ has a turning point $C_{3}$, i.e. $\text{d}(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }=0$. It follows that $C_{3}$ is also the position where $(\overline{{\mathcal{P}}})_{cv}^{\ast }$ and $(\overline{{\mathcal{D}}}_{k+\unicode[STIX]{x1D708}})_{cv}^{\ast }$ have the same rate of change, namely $\text{d}(\overline{{\mathcal{P}}})_{cv}^{\ast }=\text{d}(\overline{{\mathcal{D}}}_{k+\unicode[STIX]{x1D708}})_{cv}^{\ast }$. Formally,

(4.6)$$\begin{eqnarray}C_{3}=\min (Wo^{2}:\text{d}(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }=0).\end{eqnarray}$$

The physical meaning of $C_{3}$ is $Wo^{2}$ where $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$ reaches a maximum and also where dissipation starts to grow faster than production.

It follows that the four distinct regions are

1. (i) $Wo^{2}<C_{1}$: total dissipation decreases, the efficiency of the device rises and energy flux grows faster than production;

2. (ii) $C_{1}<Wo^{2}<C_{2}$: efficiency rises, total dissipation rises and production grows faster than energy flux;

3. (iii) $C_{2}<Wo^{2}<C_{3}$: efficiency drops but production grows faster than dissipation; and

4. (iv) $Wo^{2}>C_{3}$: dissipation rises faster than production.

4.2 Parametric study: Helmholtz number

The Helmholtz number represents the length of the channel relative to the wavelength. A thermoacoustic regenerator is generally very short compared with the wavelength: $He\sim O(0.01)$. In this section, the effects of $He$ on disturbance energy budget are studied by varying the channel length $L$ while maintaining other conditions.

The terms in the disturbance energy balance equation are shown in figure 8(a) against the Helmholtz number. Similar to figure 4, the heights of the light, medium and dark grey areas represent dimensionless values of disturbance energy flux $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$, thermal dissipation $(\overline{{\mathcal{D}}_{k}})_{cv}^{\ast }$ and viscous dissipation $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }$, respectively, of which the sum should equal the thermoacoustic production $(\overline{{\mathcal{P}}})_{cv}^{\ast }$ (dashed line in figure 8). The drop of thermoacoustic production with $He$ is expected because a larger length reduces the temperature gradient. Compared to figure 4, no turning point of $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$ is observed in the computed range; $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$ decreases monotonically. This result seems to suggest that shorter regenerators are always preferable; however, in practice, the conductivity of the material needs to be taken into account. If the length is too short, most of the heat input would be conducted through the solid wall instead of being transferred to the fluid. The deviation from equality $\unicode[STIX]{x1D6E5}^{\ast }$ against $He$ is shown in figure 8(b); again excellent agreement is seen.

Figure 8.

The balance of dimensionless disturbance energy terms against $He$. (a) The values of the dimensionless thermoacoustic production in the control volume $(\overline{{\mathcal{P}}})_{cv}^{\ast }$ (– – –), viscous dissipation $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }$ (▪), thermal dissipation $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }$ (

) and disturbance energy flux $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$ (

). (b) The deviation from equality defined as $\unicode[STIX]{x1D6E5}^{\ast }=(\overline{{\mathcal{P}}})_{cv}^{\ast }-(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }-(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }-(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$.

The minimum total dissipation $C_{1}$ was obtained at $He=0.014$, a combined effect of decreasing thermal dissipation and rising viscous dissipation with $He$ as shown in figure 9(a). Hence the main loss mechanism for a long regenerator is friction while for short ones it is irreversible heat conduction as visualized in figure 9(b).

Figure 9.

The changes of (a) the dimensionless thermal dissipation $(\overline{{\mathcal{D}}_{k}})_{cv}^{\ast }$ (——), viscous dissipation $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }$ (– – –) and total dissipation $(\overline{{\mathcal{D}}}_{k+\unicode[STIX]{x1D708}})_{cv}^{\ast }$ (— ⋅ —) and (b) the viscous dissipation (▪) and thermal dissipation (

) as a percentage of total dissipation against $He$.

Figure 10.

The changes of the magnitude of (a) the dimensionless temperature fluctuation $|T^{\prime \ast }|$ and (b) the dimensionless entropy fluctuation $|s^{\prime \ast }|$ against dimensionless axial position $z/L$ on the axis $r=0$ for four different $He$: $He=0.008$ (— ⋅ —), $He=0.012$ (——), $He=0.019$ (– – –), $He=0.031$ ($\cdots \cdots$).

The viscous dissipation is expected to rise with increasing $He$ because at high $He$, the boundary layer extends further in the axial direction, consequently leading to more dissipation. The decline of thermal dissipation with increasing $He$ can be traced back to $T^{\prime \ast }$ as shown in figure 10(a). It is clear that the magnitude of $T^{\prime \ast }$ diminishes as $He$ increases and, given that $R$ is constant, so does $|\unicode[STIX]{x2202}T^{\prime \ast }/\unicode[STIX]{x2202}r^{\ast }|$. According to (3.44), the thermal dissipation also diminishes.

The thermoacoustic production, by (3.47), depends on $|s^{\prime \ast }|$, $|v_{z}^{\prime \ast }|$ and $\unicode[STIX]{x1D719}_{vs}$. The main contributor in this case was $|s^{\prime \ast }|$, which also decreases with increasing $He$ as shown in figure 10(b). This is because the solutions of $T^{\prime \ast }$ and $s^{\prime \ast }$ are

(4.7)$$\begin{eqnarray}\displaystyle & \displaystyle T^{\prime \ast }=\frac{1-{\mathcal{H}}_{k}}{\unicode[STIX]{x1D70C}_{0}^{\ast }}p^{\prime \ast }-\frac{1}{He^{2}}\frac{1}{\unicode[STIX]{x1D70C}_{0}^{\ast }}\frac{\text{d}T_{w}^{\ast }}{\text{d}z^{\ast }}\frac{(1-{\mathcal{H}}_{k})-Pr(1-{\mathcal{H}}_{\unicode[STIX]{x1D708}})}{1-Pr}\frac{\text{d}p^{\prime \ast }}{\text{d}z^{\ast }}, & \displaystyle\end{eqnarray}$$

(4.8)$$\begin{eqnarray}\displaystyle & \displaystyle s^{\prime \ast }=\frac{-{\mathcal{H}}_{k}}{\unicode[STIX]{x1D70C}_{0}^{\ast }T_{0}^{\ast }}p^{\prime \ast }-\frac{1}{He^{2}}\frac{1}{\unicode[STIX]{x1D70C}_{0}^{\ast }T_{0}^{\ast }}\frac{\text{d}T_{w}^{\ast }}{\text{d}z^{\ast }}\frac{(1-{\mathcal{H}}_{k})-Pr(1-{\mathcal{H}}_{\unicode[STIX]{x1D708}})}{1-Pr}\frac{\text{d}p^{\prime \ast }}{\text{d}z^{\ast }}. & \displaystyle\end{eqnarray}$$

It can be seen the second parts of the solutions contain $\text{d}T_{w}^{\ast }/\text{d}z^{\ast }$; they are the temperature and entropy fluctuation originated from the temperature gradient. In the parametric study, $\text{d}T_{w}^{\ast }/\text{d}z^{\ast }$ are identical across $He$: the same cold and hot temperature varying linearly from 0 to 1. It is $He$ that accounts for the temperature gradient changes in the dimensionless domain. Since both solutions contain the inverse of $He$, the magnitudes drop with increasing $He$.

Figure 11.

The changes of the thermoacoustic efficiency against $He$.

The thermoacoustic efficiency $\unicode[STIX]{x1D702}$ is shown against $He$ in figure 11. Although there exists the highest efficiency $C_{2}$ at $He=0.012$, the change of the efficiency is rather insensitive to $He$ compared with that against $Wo^{2}$ in figure 7. This indicates that as long as $He$ is small (${\sim}O(0.01)$), the effect of the channel length (and consequently temperature gradient) on the efficiency is less important than the effects of shear and heat transfer in the radial direction.

As mentioned previously, the wave growth decreases monotonically with larger $He$; in practice the choice of $\unicode[STIX]{x1D714}$ and $L$ should ideally be such that the resulting $He$ falls to the left of $C_{1}$, i.e. $He<C_{1}$, taking advantage of both high efficiency and high output.

4.3 Two-dimensional variation

Up to this point, the global values (volume integral $(\overline{{\mathcal{P}}})_{cv}^{\ast }$, $(\overline{{\mathcal{D}}_{k}})_{cv}^{\ast }$, $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})_{cv}^{\ast }$ and surface integral $(\boldsymbol{n}\boldsymbol{\cdot }\overline{I})_{cs}^{\ast }$) and the integral form of the disturbance energy balance relation (3.35) have been used. In this section, the differential form of disturbance energy (2.27) is employed to analyse the two-dimensional variation of the point-based values $(\overline{{\mathcal{P}}})^{\ast }$, $(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{I})^{\ast }$, $(\overline{{\mathcal{D}}_{k}})^{\ast }$ and $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})^{\ast }$.

Compared with previous one-dimensional models, one of the highlights of this work is that the problem is considered in an axisymmetric cylindrical domain. Hence, this work is capable of showing the two-dimensional variation of the disturbance energy terms against radial and axial location. The profiles of $(\overline{{\mathcal{P}}})^{\ast }$, $(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{I})^{\ast }$, $(\overline{{\mathcal{D}}_{k}})^{\ast }$ and $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})^{\ast }$ are given in figure 12. The vertical axis is the dimensionless radius $r^{\ast }=r/R$, where $r^{\ast }=0$ represents the axis and $r^{\ast }=1$ is the wall.

Figure 12.

The profiles of the cycle-averaged (a) dimensionless thermoacoustic production $(\overline{{\mathcal{P}}})^{\ast }$, (b) viscous dissipation $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})^{\ast }$, (c) disturbance energy flux density $(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{I})^{\ast }$ and (d) thermal dissipation $(\overline{{\mathcal{D}}_{k}})^{\ast }$ at six axial locations: $z/L=0.02$, 0.2, 0.4, 0.6, 0.8, 0.98 which are denoted by ○, ▫, ♢, ▵, ▹, ◃, respectively, at the baseline point.

The focus first is on the radial variation. At the wall, $(\overline{{\mathcal{P}}})^{\ast }$ is zero because $v_{z}^{\prime \ast }$ is zero (no-slip boundary condition). Both $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})^{\ast }$ and $(\overline{{\mathcal{D}}_{k}})^{\ast }$ are significant because $|\unicode[STIX]{x2202}T^{\prime \ast }/\unicode[STIX]{x2202}r^{\ast }|$ and $|\unicode[STIX]{x2202}v_{z}^{\prime \ast }/\unicode[STIX]{x2202}r^{\ast }|$ are both large at the wall. Consequently, the near-wall region actually sees a negative energy flux density $(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{I})^{\ast }=(\overline{{\mathcal{P}}})^{\ast }-(\overline{{\mathcal{D}}}_{k+\unicode[STIX]{x1D708}})^{\ast }<0$, representing a wave decay. The $(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{I})^{\ast }=0$ level mark is shown in figure 12(c) for clearer visualization. It is interesting to realize that even for thermoacoustic devices designed to amplify waves, a small region of wave decay always exists. The near-wall wave decay also provides support to the conclusion in § 4.1 that overly small $Wo^{2}$ is not optimal. On the axis, $(\overline{{\mathcal{P}}})^{\ast }$ reaches the maximum while dissipation is zero (axisymmetric condition) and hence $(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{I})^{\ast }=(\overline{{\mathcal{P}}})^{\ast }$ on the axis. The fact that wave growth is uneven in the radial direction implies the presence of radial mixing loss associated with fluid near both ends entering and leaving the channels during oscillation.

With regard to the axial variation, it is too simplistic to assume the shapes of the profiles are self-similar. Take $\overline{{\mathcal{P}}}$ for example. Each profile is normalized using respective values on the axis as shown in figure 13. It can be seen that the profiles are almost self-similar but not quite. This is because radial distribution functions contain the local viscous and thermal boundary-layer thicknesses $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D708}}$ and $\unicode[STIX]{x1D6FF}_{k}$, which ultimately depend on $T(z)$.

Figure 13.

The profiles of the thermoacoustic production normalized by its value on the axis $\overline{{\mathcal{P}}}/(\overline{{\mathcal{P}}})_{r\,=\,0}$ against dimensionless radial position $r/R$ at six axial locations: $z/L=0.02$, 0.2, 0.4, 0.6, 0.8, 0.98 which are denoted by ○, ▫, ♢, ▵, ▹, ◃, respectively, at the baseline point.

The axial changes were monotonic: $(\overline{{\mathcal{P}}})^{\ast }$, $(\overline{{\mathcal{D}}_{k}})^{\ast }$ and $(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{I})^{\ast }$ decreased as the axial location moves towards the hot end; on the contrary, $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})^{\ast }$ increased. Recalling equations (3.46) and (3.44), the changes of $(\overline{{\mathcal{D}}_{k}})^{\ast }$ and $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})^{\ast }$ are readily explained by the profiles of $|\unicode[STIX]{x2202}T^{\prime \ast }/\unicode[STIX]{x2202}r^{\ast }|$ and $|\unicode[STIX]{x2202}v_{z}^{\prime \ast }/\unicode[STIX]{x2202}r^{\ast }|$ as shown in figure 14. Note that although $k(T_{w}(z))$ increased with axial location (actually $k(L)=1.41k(0)$), the factor $k^{\ast }(1/T_{0}^{\ast })$ in (3.44) actually decreased with axial location given the power law $k\propto T^{0.716}$ as mentioned previously.

Figure 14.

The profiles of the magnitude of (a) the dimensionless radial derivative of temperature fluctuation $|\unicode[STIX]{x2202}T^{\prime \ast }/\unicode[STIX]{x2202}r^{\ast }|$ and (b) the dimensionless radial derivative of axial velocity fluctuation $|\unicode[STIX]{x2202}v_{z}^{\prime \ast }/\unicode[STIX]{x2202}r^{\ast }|$ against dimensionless radial position $r/R$ at six axial locations: $z/L=0.02$, 0.2, 0.4, 0.6, 0.8, 0.98, which are denoted by ○, ▫, ♢, ▵, ▹, ◃, respectively, at the baseline point.

By referring to the thermoacoustic production (3.47), the effect of $\text{d}T_{0}^{\ast }/\text{d}z^{\ast }$ vanished here by the choice of linear temperature distribution. The factor $-\overline{s^{\prime \ast }v_{z}^{\prime \ast }}$ is the specific entropy transported by the volumetric flow rate per unit area per unit time. The profiles of $-\overline{s^{\prime \ast }v_{z}^{\prime \ast }}$ are shown in figure 15. It is clear that $-\overline{s^{\prime \ast }v_{z}^{\prime \ast }}$ increased as axial position moves towards the hot end. However, this trend was reversed by $\unicode[STIX]{x1D70C}_{0}^{\ast }$, which decreases with $z$.

Figure 15.

The profiles of the cycle-averaged product of dimensionless fluctuating entropy and axial velocity $-\overline{s^{\prime \ast }v_{z}^{\prime \ast }}$ against dimensionless radial position $r/R$ at six axial locations: $z/L=0.02$, 0.2, 0.4, 0.6, 0.8, 0.98, which are denoted by ○, ▫, ♢, ▵, ▹, ◃, respectively, at the baseline point.

In summary, it has been shown that there is a near-wall wave decay region and the two-dimensional variation of $(\overline{{\mathcal{P}}})^{\ast }$, $(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\overline{I})^{\ast }$, $(\overline{{\mathcal{D}}_{k}})^{\ast }$ and $(\overline{{\mathcal{D}}_{\unicode[STIX]{x1D708}}})^{\ast }$ can be summarized as similar radial profiles and monotonic axial changes.