Nonlinear theory for thermoacoustic waves in a narrow channel and pore subject to a temperature gradient

A nonlinear theory for thermoacoustic waves in a gas-filled, narrow channel and pore subject to an axial temperature gradient is developed based on the fluid dynamical equations for an ideal gas. Under the narrow-tube approximation, three small parameters are introduced as asymptotic parameters, one being the ratio of a span length to a typical thickness of the thermoviscous diffusion layer, another the ratio of the typical propagation speed of thermoacoustic waves to an adiabatic sound speed and the final parameter is the ratio of the typical magnitude of a pressure disturbance to uniform pressure in a quiescent state. No thermal interaction between the gas and the solid wall is taken into account on assuming that the wall has a large heat capacity. Using the three small parameters, the fluid dynamical equations are approximated asymptotically to be reduced to a single nonlinear diffusion wave (advection) equation for an excess pressure. All field variables are determined consistently in terms of the excess pressure so as to satisfy the boundary conditions on the wall. Supposing a time-periodic solution to the equation derived, the mean value of the excess pressure over one period is examined. It is shown that while the mean vanishes in the linear theory, it decreases monotonically due to nonlinearity. It is also shown that mean values of the shear stress and the heat flux at the wall, as well as those of the vector fields of the mass and energy fluxes representing, respectively, acoustic and thermoacoustic streaming, are expressed in terms of the mean values of the products of the spatial and/or temporal pressure gradients, which are reduced to the spatial derivatives of the mean pressure.


766
N. Sugimoto diffusion prevails quickly and spanwise over a period of the acoustic wave. In the motion of gas, a shear stress due to viscous diffusion dominates over fluid inertia and a pressure gradient balances with the wall friction. The resultant flow is similar to Poiseuille flow of incompressible fluids driven by a pressure difference. On the other hand, the thermal diffusion sets the gas to be in thermally perfect contact with a solid wall, yielding an isothermal process locally.
The above physical consideration suggests that the axial velocity of the gas is proportional to the pressure gradient and directed downward, while a density variation in the isothermal process is proportional to a pressure variation. Using these relations in the linearised equation of continuity averaged over a cross-section of a flow passage, it is reduced to a diffusion equation for the pressure variation. This is the basis of the statement made at the beginning. Such a picture is valid when the channel or the pore is free from a temperature gradient.
If a temperature gradient is present, the situation is dramatically changed. Sugimoto (2010) derived the thermoacoustic wave equation which describes generally the spatiotemporal behaviour of thermoacoustic waves in a channel and a pore. This is the linear equation valid in the time domain and corresponds to an equation transformed from the Rott equation derived in the frequency domain (Rott 1969;Swift 2002). Taking a limit that the frequency of the thermoacoustic waves is so low that the thermoviscous diffusion layer is much thicker than the span length, the equation is reduced to a diffusion wave (advection) equation. It is thus revealed that the gradient gives rise to wave propagation (or advection) directed in the sense of the temperature gradient.
The diffusion wave equation may be derived not only in the case of a simple geometry such as the channel and the pore, but also in complicated cases such as in the regenerative heat exchangers (called simply regenerators) used in thermoacoustic devices. Flow passages in the regenerators are tortuous, but because of the slow motion of the gas, the mean velocity over a cross-section of the regenerator is considered to be proportional to the pressure gradient. This is the very case to which the empirical Darcy's law is applicable. Although the coefficient of permeability in Darcy's law depends on temperature and also on porosity, the diffusion wave equation is expected to be derived from the equation of continuity. Such a picture, based on the diffusion wave equation, would be applicable as long as the magnitude of the pressure disturbances remains very small and the span remains very narrow.
In recent thermoacoustic devices, however, the pressure has attained the level of 10 % of the mean pressure so that finite effects in the magnitude are no longer negligibly small, while the span length is not always small enough. Then, an appreciable temperature variation arises so that no isothermal process can be maintained. While the latter effects are identified in the thermoacoustic wave equation, the former finite effects are unknown. One effect due to nonlinearity is the emergence of slow and steady streaming of mass and energy, which affect crucially both the power and efficiency of devices. Regarding acoustic streaming by a standing wave in a closed resonator, the literature is extensive. For example, see Nyborg (2008) and a review by Boluriaan & Morris (2003) in which many key papers and reviews so far published are cited. In particular, the latter review may be the first in referring to a new type of acoustic streaming such as Gedeon streaming (Gedeon 1997) in the thermoacoustic context.
Impacts of a temperature gradient on acoustic streaming were first studied by Rott (1974) for the case of a thin diffusion layer, in which the temperature dependence of viscosity and heat conductivity is also taken into account. This is Rayleigh streaming by a standing wave driven by the boundary layer in the presence of a temperature gradient. With historical perspective introduced, Bailliet et al. (2001) developed the analysis for Rayleigh streaming due to a thermoacoustic standing wave in a closed resonator, and also performed experiments to measure the acoustic streaming. It was shown that the impacts of the temperature gradient are not significant despite expectations, and the analytical results agree with the experimental ones. Later Moreau, Bailliet & Valière (2009) studied the influence of a stack on Rayleigh streaming and found new streaming vortices just outside the stack end, while the same authors (2008) developed the experimental techniques to measure streaming by laser Doppler velocimetry (LDV) in view of application to thermoacoustic devices, though in the case without a temperature gradient. Thompson, Atchley & Maccarone (2005) investigated the influences of a thermoacoustically induced temperature gradient and fluid inertia on the acoustic streaming in a standing wave. They reported that the streaming cell and its magnitude depend sensitively on the thermal condition imposed on the tube wall and therefore the axial temperature gradient induced along it, except for the isothermal (fixed) case. Although the influence of the temperature gradient is noticeable, it appears to suppress the streaming in some places and deforms significantly streaming cells. They also reported that the streaming velocities do not agree with any available theories. In the case of a looped tube, Gusev et al. (2000) made an analysis of the streaming and derived an equation which determines the streaming velocity. It is shown, however, that the sense of the acoustic streaming is sensitive to the thickness of the viscous diffusion layer.
For the energy flux, a similar steady streaming occurs, especially of the enthalpy flux, called 'thermoacoustic streaming', which occurs generally (Merkli & Thomann 1975;Rott 1980). In the isothermal process, the enthalpy variation vanishes locally and so does the enthalpy flux. By virtue of the thermodynamic relation, a pressure variation is then proportional to an entropy variation with the sign reversed so that the acoustic energy flux vector (or so-called acoustic intensity) is set equal to the entropy flux vector in the opposite sense. According to the terminology in thermoacoustics, the former is called the work flow, while the latter is called the heat flow. Although these fluxes are of quadratic order in amplitude, they can be evaluated in terms of linear solutions only and therefore no nonlinear effects are included.
Effects of nonlinearity are not restricted to acoustic and thermoacoustic streaming. Recently, it is shown by Biwa, Takahashi & Yazaki (2011) that the shock wave emerges in a looped tube. When the shock wave hits the regenerator, it is unknown what effect will be brought about in it. The ends of the regenerators play a very important role in efficiency and power because vortical motions are generated by sudden contraction and expansion of flow passages and both heat conduction and heat transfer are enhanced near the edges. Jet-driven streaming also occurs. These lead to so-called minor losses (Swift 2002), though they are major in their effects. However, such a complicated field would be limited in the vicinity of the end, and the end effects would disappear towards the inside of the regenerator. While the far-field behaviour of the gas from the end and in the outside of the regenerator is known to some extent even in a nonlinear regime, its behaviour in the inside of the regenerator is unknown. Little is known about the nonlinear behaviour of the gas inside so analytical investigations are desirable. Hence the purpose of this paper is to develop a nonlinear theory for thermoacoustic waves in such a simple case in the geometry of a narrow channel or pore, and to demonstrate its application to analysis of acoustic and thermoacoustic streaming.
The theory is developed based on a system of fluid dynamical equations for an ideal gas. Non-uniformity in the temperature distribution is taken into account through the narrow-tube approximation. This rests on an assumption of a gentle temperature gradient along the wall axially so that the first-order derivative of the wall temperature is taken into account but the second-order one is neglected (Sugimoto 2010). In addition, introducing three small parameters, one being the ratio of the span length to the thickness of the thermoviscous diffusion layer, another the ratio of the propagation speed of the thermoacoustic waves to an adiabatic sound speed and the final one being the ratio of the magnitude of the pressure disturbance to the uniform pressure in a quiescent state, the system of equations is approximated asymptotically to be reduced to a single equation for the excess pressure. In this process, all field variables are sought consistently so as to satisfy the boundary conditions on the wall. No thermal interaction between the gas and the solid is taken into account on assuming that the wall has so large a heat capacity that no temperature fluctuations occur on the wall. The linear diffusion wave equation is derived to the lowest-order approximation in these three parameters. Proceeding to a higher-order approximation, a nonlinear diffusion wave equation is derived.
In what follows, a set of fluid dynamical equations for an ideal gas is presented in § 2. On the basis of these equations, a weakly nonlinear theory is developed in § 3 for the case of a channel. It is an asymptotic theory with respect to the three small parameters based on the narrow-tube approximation. The fluid dynamical equations are approximated asymptotically so as to satisfy the boundary conditions on the wall and all field variables are expressed in terms of the excess pressure. The nonlinear diffusion wave equation is derived finally from consistency with the boundary conditions on the wall. For the case of a pore, the results are given in § 4 because this analysis is made in parallel to the one for the channel. Supposing a time-periodic solution to the equation derived, general relations among the mean values over the period are examined in § 5. Using these results, the mean values of the shear stress, the heat flux on the wall and the vector fields of the mass and energy fluxes are evaluated in terms of the mean values of products of the spatial and/or temporal gradients of the excess pressure. In § 6, some discussions are given on the evaluation of these mean values.

Basic equations and boundary conditions
At the outset, we summarize the basic equations for compressible Newtonian fluids obeying the ideal gas law, gravity being ignored. They consist of the equations of continuity, of motion, i.e. Navier-Stokes equation and of energy, which are given, respectively, by and where D/Dt stands for the Lagrangian derivative defined by ∂/∂t + v · ∇, t being the time; ρ, v, p and T denote, respectively, density, velocity vector, pressure and absolute temperature, while µ, µ v and k denote, respectively, a shear viscosity, a bulk viscosity and a thermal conductivity, c p being the specific heat at constant pressure; e and I denote, respectively, the rate of stain tensor and the unit tensor, and the components of e will be given later in terms of those of v; Φ denotes the viscous dissipation function defined by with tr(e) = ∇ · v, where tr(· · ·) designates taking the trace of the tensor, and the viscous stress tensor σ is related to e through σ = 2µ[e − tr(e)I/3] + µ v tr(e)I. (2.5) The temperature dependence of µ, µ v and k is taken into account by a power law of the form where quantities with a subscript 0 imply values in a reference state, and β is a positive constant between 0.5 and 0.6 for air. In addition to these equations, the ideal gas is stipulated by the equation of state given by where R is the gas constant. Using this, the equation of continuity is rewritten as On the other hand, boundary conditions on the surface of solid wall require that v = 0 and T = T w on wall surface, (2.9a,b) where T w represents the temperature on the wall surface. It is assumed that the heat capacity of the solid wall is so large that no temperature variations take place in the solid (for the effects of heat conduction, see Sugimoto & Hyodo (2012)). Hence T w is given and fixed. Boundary conditions at both ends of a channel and a pore are left open, assuming that they are long axially.

Nonlinear theory for a gas in a channel
3.1. Dimensionless parameters and approximation We consider a gas-filled channel of width 2H confined between two flat and parallel walls, as shown in figure 1. The x-axis is taken along the channel, while the y-axis is taken normal to it with the origin at a midpoint between the plates. This illustration depicts not only the channel but also a circular pore of radius R, to be treated later.
Supposing that the span is so narrow that the diffusion layer is thick enough to fill the channel spanwise, we develop a weakly nonlinear theory for disturbances to the quiescent state. Since gravity is neglected, the pressure in this state is uniform throughout the gas, and is denoted by p 0 . In the first place, we remark that the 770 N. Sugimoto

Wall
Gas q x o s FIGURE 1. (Colour online) Illustration of a gas-filled channel of width 2H between two flat and parallel walls subjected axially to a non-uniform temperature distribution, where the x-axis is taken along the channel and the y-axis is taken normal to it. This illustration includes a gas-filled pore of radius R if y and ±H are replaced, respectively, with the radial coordinate r and the radius R. Here T w denotes the wall temperature, while T e and ρ e denote, respectively, the gas temperature and density in a quiescent state at uniform pressure p 0 ; p , ρ and T denote, respectively, excess quantities in pressure, density and temperature, while u and v denote, respectively, axial and spanwise (radial) components of the velocity vector, s and q designating, respectively, the shear stress acting on the gas in contact with the wall surface and the heat flux flowing into it. present theory is based on the narrow-tube approximation stipulated by H/L 1, where L represents the typical length in the x-direction. The theory takes account of quantities of order H/L but neglects those of (H/L) 2 . Consider the wall temperature T w (x), which varies in the x direction so gradually that the following inequalities are satisfied: Because T w |dT w /dx| −1 specifies a typical axial length of the temperature gradient, the middle term is small, of order H/L, while the term on the left-hand side is of order (H/L) 2 and is smaller. When the second-order derivative of T w is ignored, it is shown that the gas temperature in the quiescent state, denoted by T e , may be regarded as being spanwise uniform and may be identified with T w (Sugimoto 2010). Then distinction between T w and T e is unnecessary and T e is used instead of T w hereafter. When T e is spanwise uniform, so also is the gas density ρ e because ρ e T e must be equal to ρ 0 T 0 by Charles' law. Hereafter, the quantities with subscript 0 represent those at temperature T 0 in the quiescent state under pressure p 0 , while quantities with the subscript e represent those in the quiescent state determined in terms of T e so that they are functions of x.
Besides the narrow-tube approximation, we are concerned with the slow dynamical process of a small, but finite magnitude, pressure disturbance. Letting a typical angular frequency of the disturbance be ω, the process is characterized by the following three parameters: Nonlinear theory for thermoacoustic waves in a narrow channel and pore 771 where ν, a and p represent, respectively, typical values of kinematic viscosity µ/ρ, adiabatic sound speed and the magnitude of the excess pressure p (= p − p 0 ). In the case of the circular pore, H should be read as R. The first parameter δ measures the thickness of the viscous diffusion layer √ ν/ω relative to the half-width of the channel. There also appears the thermal diffusion layer of thickness √ κ/ω, κ (= k/ρc p ) being a typical thermal diffusivity. Because the Prandtl number Pr (= ν/κ) is of order unity, for example 0.72 for air, both diffusion layers may be regarded as being comparable in thickness, though the thermal one is a little thicker.
The parameter δ 2 may also be regarded as the ratio of the typical period of slow temporal variations ω −1 relative to the time for the viscous diffusion to spread spanwise over H, i.e. H 2 /ν. When the diffusion layer is thick, disturbances are propagated neither with the adiabatic sound speed nor with the isothermal one, though the isothermal process takes place. Hence both the adiabatic and isothermal sound speeds have no physical significance in this situation. When the temperature gradient is present, as was already shown in Sugimoto (2010), disturbances are not only diffused but also propagated axially in the sense of the positive gradient. The second parameter χ measures this slow speed of propagation relative to the adiabatic sound speed. The last parameter ε measures the magnitude of the disturbances. This is the parameter measuring the order of nonlinearity. Relationships of these parameters with the Reynolds number and the Mach number will be discussed later.
Because ν and a depend on temperature, δ and χ may be defined based on local values of ν e (= µ e /ρ e ) and a e [= √ where µ e denotes the viscosity at the temperature T e at a position of x and γ (≡ c p /c v ) denotes the ratio of specific heats, c v being the specific heat at constant volume. These local values of δ and χ may be designated by attaching the subscript e. When δ and χ are defined in the reference state at p 0 and T 0 , their local values are related through Although δ e and χ e change with position, they may be regarded as being of order δ and χ, unless T e /T 0 is extremely large or small. In the following, T e /T 0 is assumed to be of order unity so that no subscript e is attached to δ and χ.

Lowest-order approximation
We consider the lowest-order approximation in the limits of δ → ∞, χ → ∞ and ε → 0. In these limits, fluid inertia is negligible and the wall friction dominates to counteract the pressure gradient. On the other hand, since the solid wall has a large heat capacity, the heat conduction suppresses gas temperature fluctuations from the local wall temperature so that T (= T − T e ) vanishes. Then the viscosity and thermal conductivity at position x take the values µ e , µ ve and k e at T = T e .
Designating a disturbance to the quiescent state by attaching a prime to a variable, and assuming the disturbance to depend on x, y and t, the equation of continuity (2.1) is linearised around the quiescent state as where v is thus set to v [= (u , v )] and use has been made of the relation ρ −1 e dρ e /dx = −T −1 e dT e /dx, since ρ e T e is constant in x. From this equation, v is estimated to be smaller than u by H/L.
When the temperature dependence of the viscosity is taken into account, the x-and y-components of the equation of motion (2.2) are given, respectively, by (A 1) and (A 2) in appendix A. Equation (A 1) is approximated by neglecting the inertia terms to be where the inertia and the other linear terms neglected are small, of order 1/δ 2 and (H/L) 2 , respectively, in comparison with the terms retained, and µ is taken to be the value at T = T e . It is estimated from this that u is comparable with H 2 p /µL. Using this, (A 2) for the y-component is similarly approximated by neglecting (H/L) 2 to be Thus p is uniform over the channel width so that u is obtained from (3.5) as This parabolic distribution in y is simply the velocity profile due to Poiseuille flow. When approximating the equation of energy (2.3), T /T 0 is assumed momentarily to be comparable with p /p 0 by the linearised equation of state (2.7), i.e. p /p 0 = ρ /ρ e + T /T e , though T is anticipated to be negligible physically because the gas is always in thermally perfect contact with the wall. Then (2.3) is approximated to be 0 = ∂ ∂y k e ∂T ∂y . (3.8) In deriving this, ρc p ∂T/∂t and ∂p/∂t are negligibly small, of order 1/δ 2 , and the other linear terms are of order (H/L) 2 relative to (3.8). The advective term u dT e /dx is of χ 2 /δ 2 relative to ∂T /∂t, but its magnitude is unknown at the present stage because the relation between δ and χ is unknown. As shown later, however, this is negligible. It then follows from the boundary conditions at x = ±H that T = 0, as expected. Thus the equation of state yields the isothermal relation locally It is found from this and (3.4) that u is of order εωL. Because u is stipulated by (3.5) to be of order H 2 p /µL and p /p 0 is regarded as being of order ε, it follows that where ∼ implies equality in order of magnitude and the numerical factor γ (and Pr later) are regarded as being of order unity. Thus χ turns out to be comparable with δ. If (3.7) and (3.9) are substituted into (3.4), v will be obtained. Before doing this, integration of (3.4) over the channel width with use of the boundary condition v = 0 at y = ±H yields where the overline stands for the average over the channel width defined, for example, for u by Using (3.7) and (3.9), p is governed by the following diffusion wave equation: where α e is the diffusivity of the gas confined in the channel and is given by (3.14) It is found that because χ ∼ δ, the first term in (3.13) balances with the second term. For a moderate temperature gradient, the third term is also of the same order as the second term. If the temperature dependence of the shear viscosity is neglected, then α e is a constant independent of x. Equation (3.13) is the same as the one derived from the thermoacoustic wave equation for a thick diffusion layer (Sugimoto 2010). This equation means that p is not only diffused by the second term but also is propagated by the third term with a velocity (α e /T e ) dT e /dx. The explicit form of v is obtained by integrating (3.4) with respect to y as v = − ∂ ∂t where (3.13) has been used. Thanks to this, it is confirmed that v satisfies the boundary conditions v = 0 at y = ±H. Furthermore v is expressed in a compact form as where v is odd in y and |v | has peaks at y = ±H/ √ 3. Although u has the same profile as that of Poiseuille flow of incompressible fluids, note that v exists due to the compressibility of the gas. It is verified by (3.15) that v is of order εωH and is smaller than u by the first order of H/L. What is also to be noted is that u is determined by the spatial pressure gradient as expected, whereas v is determined by the temporal pressure gradient.
With all lowest-order terms available, we check the relations among them. The acoustic Mach number Ma of the axial velocity relative to the local adiabatic sound speed a e is defined by Thus the acoustic Mach number is much smaller than ε. The quantity p /u defines an acoustic impedance, and this is found to be of order χρ e a e . Hence the characteristic acoustic impedance is much larger than ρ e a e for a travelling wave in a 774 N. Sugimoto lossless medium. Next we check the Reynolds number Re. Taking a typical speed and length to be εωL and H, respectively, it is given by Hence the Reynolds number is determined by the competition among the three small parameters 1/δ, ε and H/L. Finally, we evaluate the components of the rate of strain tensor e. These are calculated straightforwardly as It is confirmed from these that the shear components e xy and e yx are much larger than the other components and that the divergence of the velocity vector is smaller than them by the factor H/L because Since this is equivalent to the last two terms of (3.4), (3.20) is consistent with (3.4) for an isothermal process with ρ /ρ e = p /p 0 .

Higher-order approximation
We now proceed to take account of the finiteness of the small parameters δ −1 , χ −1 and ε. We start with the y-component of the equation of motion and see how the uniformity of pressure over the channel width is affected. The equation is given by reviving the neglected terms in (A 2) in appendix A as where µ and µ v are expanded around T e as µ = µ e 1 + βT T e + · · · and µ v = µ ve 1 + βT T e + · · · , (3.22a,b) but since T vanishes to the lowest order, µ and µ v have been set equal to µ e and µ ve , respectively. The term (∂/∂x)(µ∂v/∂x) in (A 2) is small, of order (H/L) 2 , relative to the first term on the right-hand side of (3.21), while the inertial term ρ e ∂v /∂t is smaller by 1/δ 2 and the nonlinear advection terms are even smaller by ε/δ 2 .
With the right-hand side of (3.21) evaluated using the lowest-order expressions available, it is solved for p . Its solution, denoted by p 2 , is given by where p (x, t), which is an arbitrary integration constant and uniform over a cross-section, corresponds to the lowest-order solution. It is found from (3.23) that the pressure becomes non-uniform spanwise and exhibits a parabolic distribution in y. The second term is estimated relative to p to be 3µ e p 0 H 2 p ∂p ∂t (3.24) The third term is obviously of order (H/L) 2 . Hence, as long as the ratio (H/L) 2 is negligible, p 2 is regarded as being uniform over the channel. In the following analysis, we still maintain this assumption as the narrow-tube approximation that the pressure is uniform over the cross-section. Next we look into the x-component of (A 1) to revive neglected terms in (3.5). Among these, terms next to be included are given as follows: (3.25) With the lowest-order expressions of u and v and with T = 0, the right-hand side is evaluated in terms of these and (3.25) is solved for u . On the first line of the right-hand side of (3.25), the second term is estimated to be of order 1/δ 2 relative to the first pressure gradient, while the three nonlinear terms due to inertia are evaluated to be small, of order of ε/δ 2 or ε/χ 2 , relative to it. All terms on the second line are neglected because they are of order (H/L) 2 relative to the left-hand side. It is noted that, although the bulk viscosity is often neglected from the outset, it does not come into play even if included. When the nonlinear terms are included, whereas the second line is ignored, this means that Since (H/L) 2 is neglected, the relation is satisfied. In actual situations this will be fulfilled because ε is of 10 −2 -10 −1 , while H/L is of 10 −2 -10 −3 and δ of 10-10 2 . Under such a condition, (3.25) is solved for u with the boundary conditions u = 0 at y = ±H. It then follows straightforwardly that with φ i (i = 2, 4, 6) defined, respectively, by which is determined by the solution of d 2 φ 64 /dy 2 = −(H 2 − y 2 )y 2 /2 with φ 64 (±H) = 0. Using (3.27), we estimate the order of each term relative to the leading term. The second term on the first line is due to inertia and is of order 1/δ 2 . This reflects the finite effect of the span length, The terms on the second line are due to nonlinear advection and small, of ε/δ 2 or εχ 2 /δ 4 . Thus it is found that the nonlinear effect remains secondary in (3.27) in comparison with the finite effect of span length. From the standpoint of the asymptotic expansion, this implies that the expressions of the nonlinear terms may be inaccurate. If the lowest-order approximation is regarded as the zeroth-order approximation, the second term in (3.27) belongs to the first-order one but the nonlinear terms belong to the second-order one. To validate them, therefore, we must proceed to the second-order approximation. This will be discussed later.
Next we examine the equation of energy to specify a small temperature variation by taking higher-order terms into account. Noting that T vanishes at the lowest order, quadratic terms neglected in (3.8) are revived as On the right-hand side, the leading terms in magnitude are ρ e c p u dT e /dx and −∂p /∂t, while ρ e c p ∂T /∂t is neglected because it is smaller than the left-hand side by 1/δ 2 and T vanishes to the lowest order. Also neglected are (d/dx)(k e dT e /dx) and (∂/∂x)(k e ∂T /∂x), which are smaller than the terms retained by (H/L) 2 . The nonlinear terms arise from c p ρ u dT e /dx due to the density change advected under the temperature gradient, the advective term of pressure u ∂p /∂x and the viscous dissipation function Φ. The latter two terms are free from the temperature gradient. The function Φ is evaluated in terms of the components (2.4) and is contributed to mainly by e xy and e yx . It is given by (3.19c) as to the lowest order. Substituting the lowest-order relations in (3.30) and solving for T , it follows that where ρ e c p /k e µ e = Pr/ρ e ν 2 e with Pr = ν e /κ e , κ e being k e /ρ e c p and where Pr is independent of T e . Here we check the order of T relative to T e . The first term is estimated to be For δ ∼ χ, the second term in (3.32) is comparable with the first term. In a similar fashion, both the third and fourth terms are estimated to be of order ε 2 χ 2 /δ 4 . Hence it is found that the first two terms are comparably small, of order ε/δ 2 , relative to T e , which confirms that T vanishes at the lowest order, and that the temperature shows a parabolic distribution plus a quartic one in y. The third term in (3.32) results from the second term on the right-hand side of (3.30), while the last term results from the third and fourth terms. The last term in (3.32) is independent of the temperature gradient and is always negative. Because Φ is positive and T vanishes at y = ±H, −Φ contributes in (3.30) to making T positive, whereas because the advective term u ∂p /∂x is negative, it cancels Φ to make T negative in the form of the last term.
With the sign of the third term in (3.32) unknown yet, these nonlinear terms are of higher order than the first two terms by ε. As in the case of u , the nonlinear terms must be validated by evaluating higher-order terms.
To facilitate further reduction of equations in proceeding to the second-order approximation, the order of magnitude of all variables is summarized. On the basis of p /p 0 of order ε, it follows that Now that T has been specified, we reconsider (3.25) to the second order and check by using (3.27) and (3.32) (up to the first two terms) whether or not there are any other terms missing, but comparable with the terms taken. Because of the finite temperature disturbance, the temperature dependence of the viscosity should be included so that µ e on the left-hand side of (3.25) should be modified. From (A 1) with (3.22), it is taken into account as ∂ ∂y µ e ∂u ∂y + βµ e T T e ∂u ∂y . (3.36) The second term is comparable with the advection terms. This term is evaluated by using u in (3.7) and the lowest-order terms of T in (3.32). In addition, there arises a linear term of higher order emerging from the second term of (3.25) by taking u in ( where φ 60 is defined by a solution to d 2 φ 60 /dy 2 = −φ 4 with φ 60 (±H) = 0 and given by (5H 2 φ 4 − φ 64 )/12, while φ 62 is defined through d 2 φ 62 /dy 2 = (d/dy)(φ 4 dφ 2 /dy) with φ 62 (±H) = 0, and they are given, respectively, by The third line in (3.37) reflects the thermomechanical coupling due to the temperature dependence of the viscosity. Thus u is specified correctly up to the order of 1/δ 4 , ε/δ 2 and εχ 2 /δ 4 relative to the first term. Next to be checked is the equation of energy (3.30). Because the heat conductivity also changes with the temperature, k is also expanded in the form of (3.22). Then there comes a term βk e (T /T e )∂T /∂y into the parentheses on the left-hand side of (3.30). Yet it is found that this term is smaller than the nonlinear terms taken and is negligible. On the right-hand side, however, ρ e c p ∂T /∂t must be revived, and u in the first term should be taken up to the second term in (3.27). Then (3.32) is corrected as The second line gives the second-order corrections to the two linear terms on the first line. Thus T is specified correctly up to the order of ε/δ 4 , εχ 2 /δ 6 and ε 2 χ 2 /δ 4 relative to T e . With u and T specified, we finally consider the equation of continuity expressed in the form of (2.8). Noting that p is spanwise uniform and the nonlinear advection terms of T are smaller than u dT e /dx, (2.8) is approximated to include the nonlinear terms of the lowest order as Expanding 1/(p 0 + p ) and 1/(T e + T ) with respect to the disturbances, and taking the mean of (3.40) over the channel width, it follows that to quadratic order in the disturbances. Here (u /p 0 )∂p /∂x and (p /p 2 0 )∂p /∂t are comparable in magnitude of the lowest order, whereas (T /T 2 e )∂T /∂t and (T /T 2 e )u dT e /dx have been neglected because they are smaller than the lowest-order ones by 1/δ 2 and 1/δ 4 , respectively.
To close this equation for p , u and T are necessary. They are calculated by (3.37) and (3.39) and using the averages: φ 2 = H 2 /3, φ 4 = 2H 4 /15, φ 2 2 = 2H 4 /15, φ 6 = 2H 6 /35, φ 64 = 2H 6 /105, φ 60 = 17H 6 /315 and φ 62 = 8H 6 /315 as (3.44) This is the desired higher-order equation next to (3.13). Here it is noted that, as far as (3.44) is concerned, the higher-order expressions of (3.37) and (3.39) are unnecessary. However, it will turn out that they are required when mean values over one period of oscillations are discussed for the acoustic and thermoacoustic streaming in § 5. As the first three terms of (3.44) represent the lowest-order equation, the fourth and fifth terms represent the finite effect of the span length of 1/δ 2 relative to the lowest-order terms, while the last two terms represent the finite (nonlinear) effect of the pressure disturbance of ε or (χ 2 /δ 2 )ε. Using the lowest-order balance, the first term in the square bracket may be replaced by It is noted that the nonlinear terms are independent of the temperature gradient dT e /dx, and also of the temperature, if α e were taken as constant for β = 0. Hence it is revealed that the temperature and its gradient do not impact nonlinear terms in evolution of the pressure, and that the temperature dependence of the viscosity and heat conductivity appear only through α e and the coefficient β.
With the higher-order relations available, v is specified to the next higher order than (3.16) by integrating (3.40) with respect to y from y = 0. Executing the integrals of the functions in y involved, i.e. 1, φ 2 and φ 4 , it follows that where v vanishes at y = ±H owing to (3.44). Here it should be noted that for ∂u /∂x in (3.40), only the first two terms in (3.37) have been used because v is smaller than u by H/L. Using (3.44), (3.47) may be expressed in an alternative form as We have demonstrated a way to carry the theory up to higher order in ε. As was seen, the leading nonlinear terms result from p −1 Dp/Dt in (2.8) so that the first term and the nonlinear terms in (3.41) may be expressed without expansion with respect to p /p 0 as This form may be useful if a higher-order equation than (3.46) is required. In taking the average of v · ∇ log(1 + p /p 0 ) over the channel width, u should be replaced by (3.42), and p is assumed to be independent of y. On top of these, nonlinear terms arise from T −1 DT/Dt written as ∂ ∂t + v · ∇ log(T e + T ). (3.50) On taking the average, (3.43) is useful but spanwise variations in T need to be included together with v . The calculations are expected to be complicated but higher-order theory could be developed step-by-step without any difficulties in principle.

Case of a circular pore
For a circular pore of radius R, the analysis can be made in parallel with that presented for the channel. Thus its description is discussed briefly. Assuming axisymmetry, the x-axis is taken along the pore, while the radial coordinate r is used instead of y ( figure 1). The x-and r-components of the velocity vector v are denoted by u and v , respectively, which depend on x, r and t. The results for the channel are basically applicable if y and H are replaced by r and R, respectively, but there occur differences in the coefficients.

General relations for means of time-periodic oscillations
This section is devoted to discussion of general relations derived from (3.46) and (4.19). In application to thermoacoustic devices, they describe the nonlinear behaviour of a gas in regenerators and should be solved together with other equations valid in the domains outside of the regenerators. In this sense, it is not relevant to seek specific solutions to the equations alone. Because self-excited oscillations are generated in real devices, it is expected that the nonlinear equations (3.46) and (4.19) will be able to support time-periodic solutions. In the linear case, in fact, it is shown by Hyodo & Sugimoto (2014) that the linearised versions of those equations describe a marginal state of instability and yield a periodic oscillation, though infinitesimally small in amplitude. On the basis of an a priori assumption based on the existence of periodic solutions, general relations among mean values over one period are considered.
To discuss the two cases for the channel and the pore simultaneously, (3.46) and (4.19) are written as where the pairs of fractions (m, n) are (2/5, 6/5) and (1/6, 8/6) for the channel and for the pore, respectively and α e is defined by ma 2 e H 2 /nγ ν e (= mp 0 H 2 /nµ e ), H being replaced with R for the pore.

Relations for mean values
Denoting a time period of the oscillations in p by τ , i.e. p (x, t) [= p (x, t + τ )], a mean value of p is defined by attaching a tilde as It is then straightforward to show that Taking the mean of (5.1), it follows that where the terms associated with differentiation with respect to t drop out because of the periodicity. This relates the spatial derivatives of the mean pressure p to the mean of the pressure gradient squared. Note that the finite effect of the span length does not contribute to (5.4). Because (5.4) is written as the mean of the pressure gradient squared and multiplied with α e /T e on the left-hand side is expressed in terms of the spatial derivative of the mean pressure gradient multiplied with −α e p 0 /T e . Alternatively this implies that the gradient of the mean pressure is small, of quadratic order. In fact, ∂ p /∂x is expressed as x α e p 0 T e ∂p ∂x 2 dx. (5.6) When the nonlinear terms are neglected, the mean pressure gradient vanishes so that p is taken to be zero. However, when the nonlinear terms are taken into account, it is obvious that the mean pressure decreases with x. Next multiplying (5.1) with p and taking the mean, it follows that where P(x) stands for the higher-order terms given by where the tilde over the last term is effective with p inclusive. The first three terms in (5.7) are of the same order and each term in P is small, of the order of 1/χ 2 , 1/δ 2 and ε relative to the leading terms. Elimination of the second term in (5.7) by use of (5.4) and rearrangement of the terms yield When the right-hand side is neglected, it follows that where c 1 and c 2 are arbitrary constants in x. Because the right-hand side is independent of p , and is present even in the linear case, it is reasonable to take c 1 and c 2 to vanish. Then it is found that the mean pressure p is negative and small, of order ε 2 p 0 . The terms neglected in (5.10) (due to P) are even smaller, of order ε 2 /χ 2 , ε 2 /δ 2 and ε 3 relative to p 0 . When (5.1) is multiplied with α e ∂p /∂x and the mean is taken, the mean of the product of the temporal gradient of p and its spatial gradient is expressed in terms of the spatial derivative of the spatial gradient squared and multiplied with α 2 e /2T 2 e as Here it should be remembered that the relations (5.11) to (5.13) hold to the lowest order only, whereas (5.5) is derived rigorously from (5.1). It is found from (5.11) that the mean of the product of the spatial and temporal gradients of p is expressed, besides the factors α e and T e , in terms of the spatial derivative of the mean of the pressure gradient squared, which is related to the gradient of the mean excess pressure through (5.5). Similarly, it is found from (5.13) that the mean of the temporal gradient squared is expressed in terms of the second-order derivative of the mean of the pressure gradient squared. Hence the means of the products of the spatial and/or temporal gradients of p are related to each other through spatial differentiation and ultimately they are reduced to the spatial derivatives of the mean pressure p by (5.5).

Mean shear stress and mean heat flux on walls
Knowing the relations among the mean values, we now consider mean values of the shear stress s acting on the gas and of the heat flux q flowing into the gas on the walls. Noting that the shear viscosity and the heat conductivity on the wall are, respectively, µ e and k e because T vanishes there, s and q are defined for the channel by s = 2µ e ∂u ∂y y=H and q = 2k e ∂T ∂y y=H , (5.14a,b) where the factor 2 reflects the contribution from the lower wall by using even functions of u and T in y, and defined for the pore by s = µ e ∂u ∂r r=R and q = k e ∂T ∂r r=R .
(5.15) Substituting into these expressions u and T obtained in § § 3 and 4, s and q are expressed in terms of p . For the channel, u and T are given, respectively, by (3.37) and (3.39), while for the pore, they are given, respectively, by (4.10) and (4.15). By doing this and taking the mean, s and q are obtained for the respective cases as follows: where note that mPrk e H 2 /nρ e ν 2 e = [γ /(γ − 1)]α e /T e and γ /(γ − 1) = c p /R = c p ρ e T e /p 0 . For the case of the pore, 2H on the left-hand sides of (5.16) and (5.17) is replaced by R/2, whereas H on the right-hand side is replaced by R. The quantities on the left-hand sides of (5.16) and (5.17) represent, respectively, the total shear force and the heat flux along the wetted perimeter divided by the cross-sectional area. For a channel of unit depth perpendicular to the (x, y)-plane, the wetted perimeter in contact with the gas is 2 (which has already been included in (5.14)) and the cross-sectional area is 2H.
It is found that the mean shear stress is contributed to not only by the gradient of the mean pressure but also by the means of the products of the spatial and/or temporal pressure gradients. However, using (5.11) and (3.7) (or (4.3) in the case of the pore), (5.16) is rewritten as where u 2 = (m 2 H 4 /nµ 2 e )(∂p /∂x) 2 to the lowest-order approximation. The first and last terms represent, to quadratic order, the balance of the x component of the momentum flux density vector averaged over the span length. This form is rather more obvious physically than (5.16). Because of (5.10), the gradient of the mean pressure is small and comparable with (p ∂p /∂x)/p 0 , yet the remaining nonlinear terms are smaller than this term by 1/δ 2 and χ 2 /δ 4 . This means that the shear stress is mainly caused by the pressure gradient and the contribution from ρ e u 2 remains small.
On the other hand, (5.17) is derived directly by averaging (2.3) over the span length. The terms on the right-hand side stem from the advection under the temperature gradient, c p (ρ e + ρ )u dT e /dx, while the sum u ∂p /∂x + Φ vanishes on averaging. The relation (5.17) shows that the mean heat flux flows only when the gradient is present. To examine this, we start with the lowest-order approximation. Because p vanishes and so does the mean pressure gradient, there is no mean heat flux. This result is reasonable because T vanishes everywhere. Hence the mean heat flux is brought about by the finite effects of 1/δ, 1/χ and ε in the presence of a temperature gradient.
When the higher-order terms in (5.9) are taken into account, it is integrated as (5.19) Hence q in (5.17) is rewritten as (5.20) The mean heat flux is thus available after integration of the right-hand side. Taking a reference heat flux to be p 0 ωH, it is estimated from (5.20) that q is of order ε 2 /δ 2 , (χ 2 /δ 4 )ε 2 and (χ 2 /δ 2 )ε 3 relative to the reference value. No further reduction appears to be possible at this stage, but integration of P/T e is possible in special cases. This will be discussed later.
5.3. Mean mass flux Next we consider a mean mass flux due to the periodic oscillations, which gives rise to the acoustic streaming. (Strictly speaking, the mass flux in the present context should accompany with the term 'density vector' as the mass flux density vector, but this will be omitted below for brevity not only for the mass flux but also for the energy flux.) The mass flux ρv is split into (5.21) The excess density ρ is expressed by using the equation of state (2.7), i.e. p /p 0 = ρ /ρ e + T /T e + ρ T /ρ e T e , as with a specific enthalpy designated by h (= c p T), where h e (= c p T e ) and h (= h − h e = c p T ) denote, respectively, a local specific enthalpy in the quiescent state and an excess specific enthalpy from it. Note that since h /h e (= T /T e ) is smaller than p /p 0 , the first term dominates on the right-hand side of (5.22). At any rate, it is found from this that ρ v is related to the acoustic energy flux p v and the excess enthalpy flux ρh v .
On the other hand, however, the y-component is not calculated by (3.48) accurately, because v has not been specified up to the approximation consistent with (3.37). To do this, it is necessary to derive the higher-order equation than (3.44) by taking account of the nonlinear advection terms (the second and third lines in (3.42)) in (3.41) and to obtain the expression of v higher than (3.48) from (3.40). The procedure is the same as was used and is simple in principle but becomes very awkward. To circumvent this, it is helpful to note that the mean mass flux satisfies the divergence-free condition i.e. ∇ · (ρv ) = 0. Noting that the integral of ρu over the span length is constant in x, (5.24) ρv is calculated by the following integral of (5.23) as (5.25) where, since ρu is even in y, while ρv is odd so that it vanishes at y = 0, (5.25) satisfies the boundary conditions ρv = 0 at y = ±H owing to (5.24). For the case of the pore, ρu and ρv are replaced, respectively, by r ρu and r ρv , while y and H are replaced, respectively, by r and R. In any case, it is found that ρv contains the spatial derivatives of the terms with the tilde in (5.23).
Here we estimate the order of the mass flux. It is noted that the sum in the first parentheses in (5.23) which is proportional to φ 2 vanishes to the lowest order as in the case of the mean heat flux, and therefore it must be evaluated by using (5.9). Taking a reference magnitude of the mass flux to be ρ e ωL, and using (5.19) with (5.8), it is found that since φ 46 ∼ φ 2 2 ∼ H 4 , and the second terms are of the same order as the first terms, i.e. χ 2 /δ 2 relative to the first ones. Hence ρu is found to be of order For ρv given by (5.25), on the other hand, the integration of ρu with respect to y and then differentiation with respect to x lowers the order by H/L. Thus ρv is estimated to be ρv ∼ ρu H L .
( 5.29) It is noted that ρv includes the spatial derivatives of the factors with the tilde in (5.23) by differentiation with respect to x. Hence when the three parameters 1/δ 2 , 1/χ 2 and ε are regarded as being comparably small, the magnitude of the x-component of the mass flux is of order ε 3 ρ 0 ωL, while the y-component is smaller than it by H/L. If the mean pressure gradient is known by any means, it is straightforward to evaluate the means of the spatial and/or temporal gradients by (5.11) and (5.13). Then the vector field of the mean of the mass flux can be analysed by (5.23) and (5.25).

Mean energy flux
It has been revealed that the mean acoustic energy flux and excess enthalpy flux are engaged in the mean mass flux. Here we consider the mean energy flux. While the energy equation is given by (2.3), it is originally derived from the following conservation form: where E stands for the total energy density ρ(e + v · v/2), e being a specific internal energy c v T and stands for the energy flux density vector given by Taking the mean of (5.30) over one period, it follows that ∇ · = 0. As consists of four terms, we look at the magnitude of each term.
The first term ρ(v · v)v/2 is due to transfer of the kinetic energy. This may be approximated by ρ e u 2 v /2 because the contribution from v 2 /2 is smaller by the order H 2 /L 2 relative to u 2 /2. Taking a typical magnitude of the energy flux to be ρ 0 h 0 ωL and noting that h 0 = c p T 0 = a 2 0 /(γ − 1), the transfer of kinetic energy is estimated to be of order ε 3 χ 4 /δ 6 relative to the typical one. Because we are concerned with the energy flux at the same level of approximation as the mass flux relative to ρ 0 ωL, the transfer of the kinetic energy is neglected.
The second term represents the enthalpy flux ρhv , which is transported substantially with the mass flux. It may be viewed as the sum of the flux carrying the local value h e and the excess enthalpy h . Thus it is broken into (5.32) Using (5.22) to eliminate ρ v , the enthalpy flux is written as Comparing this with the mass flux (5.21) (with (5.22)) multiplied with h e , the only difference is due to the term h ρv . The third term in (5.31) represents the heat conduction and the leading term is given by the spanwise component. Using (3.39) and noticing that the linear terms associated with differentiation with respect to t drop out on taking the mean and the sum of the second and third terms becomes smaller than the fourth term by (5.26), it follows that where for the case of the pore, y is replaced with r. The term (5.34) is estimated to be of order ε 2 δ 2 H/χ 2 L relative to ρ 0 h 0 ωL and is found to be smaller. The fourth term represents the flux due to the viscous power. Because the leading term is due to the y-or r-component, it is given for the case of the channel by precise measurement of them. In this sense, the measurement of the spatial gradient would be difficult. As suggested by the theory, however, we know that the spatial gradient gives rise to the axial velocity. Because velocity measurements may be available relatively easily by techniques of hot-wire anemometry and LDV, the spatial gradient is available by measuring the axial velocity. Furthermore, if the maximum spanwise component of the velocity can be measured, it bears relation to the temporal gradient of the pressure. However, because the spanwise component is much smaller, it seems to be less practical. Rather the temporal gradient is available if the temperature can be measured. Therefore if the axial velocity and the temperature are measured pointwise so that the means of the products of the spatial and/or temporal pressure gradients are available, then the mean values may be obtained according to the formulae given in § § 5.2-5.5. This will be described briefly below. For the measurement of the pressure, in passing, it may be available via the density measured optically. Since the density is uniform spanwise thanks to the isothermal process to the lowest-order approximation, the measurement of the density may provide a complementary means to that of the pressure. Let temporal fluctuations of the axial velocity and the gas temperature averaged over the span be available pointwise in the axial direction. Then the averaged value of u is evaluated approximately by (3.42)  where H is replaced by R for the case of pore. Making use of these, the temporal mean values of the products of the spatial and/or temporal pressure gradients may be given by where k e /µ e = c p /Pr. Hence if the temporal mean values of the averaged axial velocity and temperature over the span were measured, then it would be possible to determine approximately all the vector fields of the mean mass and energy fluxes. Here it is emphasized that because the temperature fluctuation is associated with the mean fluxes, the lowest isothermal approximation fails to determine them.
In the case that no spanwise averaged values were available, but maximum values over the cross-section were available experimentally, (6.10) is modified. Designating the maximum value by attaching an asterisk, the maximum value of u is given for the case of the channel by (6.11) to the lowest order, while the one of T is given by (6.12) For the case of the pore, H 2 is replaced by R 2 /2. Then the relations (6.10) are replaced by the following ones: If the means of the products of the spatial and/or temporal gradients of p were evaluated via the relations above by measured data of the axial velocity and temperature, all the mean values of not only the shear stress and heat flux on the wall but also those of the vector fields of the mass and energy fluxes would be available.
When the temperature gradient is present, the above relations are modified because (6.9) and (6.12) are given, respectively by in the case of the channel. Since the respective second terms are written in terms of u or u * , the means of the products of the spatial and/or temporal pressure gradients are expressed in terms of the means of the axial velocity and/or the temperature.

Conclusion
A nonlinear theory for thermoacoustic waves in a gas-filled, narrow channel and pore subjected to a temperature gradient has been developed on the basis of the fluid dynamical equations for an ideal gas. No thermal interaction between the gas and the wall is taken into account by assuming that the heat capacity of the wall is so large that the wall temperature is fixed. The theory is asymptotic in the three small parameters 1/δ, 1/χ and ε defined by (3.2) on the basis of the narrow-tube approximation such that the terms of order (H/L) 2 are negligible so that the pressure may be regarded as being uniform spanwise. Starting with the lowest-order approximation that the pressure gradient balances with the wall friction and the isothermal process takes place everywhere locally, the fluid dynamical equations are approximated asymptotically and successively by the three parameters so as to satisfy the boundary conditions on the wall.
In the limiting case that 1/δ, 1/χ and ε tend to zero, the lowest-order approximation yields the linear diffusion wave equation for the excess pressure. The axial velocity of the gas is simply equal to that of Poiseuille flow as in incompressible fluids, though the spanwise velocity is present because of the compressibility. When the mean velocity averaged over the flow passage is concerned, it is proportional to the pressure gradient with the sign reversed. This result endorses the validity of the empirical Darcy's law for the mean velocity. On the other hand, the local temperature disturbance remains zero so that the gas temperature is equal locally to the wall temperature. Such a picture would be applicable to complicated cases of tortuous flow passage in regenerators, provided that the span length and the magnitude of the pressure disturbances are extremely small. In these cases, the diffusivity α e must be determined experimentally.
For finite values of the parameters, all field variables have been determined consistently on the basis of the lowest, i.e. zeroth-order approximation, in terms of the excess pressure up to the higher-order terms. When the three parameters satisfy the relation 1/δ 2 ∼ 1/χ 2 ∼ ε, the finite effect of span length first manifests in the axial velocity and the temperature. Because of the latter disturbance, the local isothermal process is no longer valid. To first-order approximation, the finite effect of the magnitude arises from the density change ρ −1 Dρ/Dt in the equation of continuity, especially through the pressure change p −1 Dp/Dt. Taking the lowest-order nonlinear terms, the single nonlinear diffusion wave equations (3.46) and (4.19) are derived from the consistency condition with the boundary conditions on the wall.
The linear terms in the nonlinear diffusion wave equation agree with those derived from the approximation of the thermoacoustic wave equation for a thick diffusion layer at low frequency (Sugimoto 2010). The other nonlinear terms in the velocity and temperature due to the advection terms in inertia, the temperature dependence of the viscosity and the heat conductivity in the equation of motion and the viscous dissipation function in the equation of energy remain secondary in the equations derived. To discuss the temporal mean values for time-periodic oscillations, however, the axial velocity and the temperature have been specified to second-order approximation. To the extent of the present approximation, it is found that the bulk viscosity does not come in to play even if included, though usually neglected without any reasonable arguments in spite of its being comparable with the shear viscosity.
When the theory is applied to the analysis of thermoacoustic devices, the equation derived is to be solved in conjunction with other wave equations describing nonlinear propagation in a buffer tube and a tube with uniform temperature. As a consequence of this, it is expected that self-excited, periodic oscillations would be obtained. Supposing the existence of time-periodic solutions to the nonlinear diffusion wave equation, the temporal means of the shear stress and the heat flux on the wall, and the vector fields of the mass and energy fluxes over one period of oscillation, have been expressed in terms of the means of the products of the spatial and/or temporal gradients of the excess pressure.
As the mean pressure gradient vanishes in the linear case, it is found to be small, of quadratic order, and to decrease monotonically in the axial direction. Because the order of the mean pressure gradient is thus lowered, the nonlinear terms which are secondary in the nonlinear diffusion wave equation come up in to play in the mean mass and energy fluxes. It has been revealed that the means of the products of the spatial and/or temporal pressure gradients are reduced to the spatial derivatives of the mean pressure gradient. This is an interesting outcome because acoustic and thermoacoustic streaming are determined in terms of the mean pressure gradient and its spatial derivatives only. This result is expected to provide a method to determine the streaming by measuring the means of the products of the axial velocity and/or temperature rapidly changing, and further to illuminate a clue to suppress the undesired streaming.