Evolution of Flow and Streaming in Exponential Variable Cross-Section Resonators

A gas-kinetic scheme (GKS) based on an unstructured grid is applied to simulate the evolution of the fluid motions in exponential variable cross-section resonators. The effects of the acoustic field intensity on the oscillatory pressure, velocity, temperature, and flow streaming structure were investigated numerically, and the model was validated. The results demonstrate that the geometry and driving strength are the main factors affecting the final performance of the system. For the quasi-linear and moderate non-linear cases in optimum exponential tube, the periodic generation, evolution, and shedding of vortices in flow fields are associated with the storage and release of energy, which is the transmission mode of the third type of direct current (DC) flow, and its driving mechanism is attributed to the asymmetrical pressure and temperature. Meanwhile, some new physical characteristics were also discovered for the highly non-linear case, e.g., the disorder and unsteadiness of the flow direction accomplished with turbulent flow streaming structures. The secondary flow is manifested as multiscale, irregular and unsteady vortices throughout the tube. The smallest increment of pressure and velocity amplitude occurs concurrently with the biggest increment of temperature amplitude. These evidences suggest that there is an optimal driving strength, even for a good configuration tube, with which the maximum efficiency can be obtained.


Introduction
Owing to its high reliability and environmental friendliness, thermoacoustic engines without moving parts have attracted increasing attention as a major technology to develop more efficient energy conversing or generating systems [1,2]. The working gas in thermoacoustic engines is characteristic of the alternating current (AC) gas flow, which alternately moves backward and forward with pulsating pressure at certain frequencies. Until recently, the thermoacoustic conversion efficiency had been continually improving with both novel configurations and insight into the physics of flow and heat transfer in resonators [3,4].
The geometry of the resonator is a crucial factor in determining the non-linear standing pressure waveform and amplitude generated within it. If the shape of a resonator is a straight cylinder, a shock wave accompanied by acoustic saturation will appear when the inner gas oscillates at its resonance frequency [5]. In recent years, researchers have been focusing on developing measures and approaches to suppress these non-linearities [3,4]. It is well-known that one can obtain high-amplitude and shock-free acoustic oscillations by varying the shape of resonators. This method was first presented by Lawrenson et al. [6] and is now called resonant macrosonic synthesis (RMS). Luo et al. [7] in exponential varying cross-section resonators. They found that the exponential resonators have the ability to suppress the shock waves, acoustic streaming, and vortices and gave an optimum index of winding for exponential resonators.
Stimulated by the successful applications of GKS in turbulent flow [30] and acoustic and shock-wave capturing [31], we applied the GKS method to study the problem of high-amplitude acoustic fields in shaped resonators and investigated the influence of driving strength on the flow and streaming. The purpose of this work was to obtain an optimal shape parameter and driving strength for a highly efficient acoustic resonator that has the highest pressure fluctuation amplitude, ratio of output work to input work, and the capacity of suppressing dissipation effects (including waveform distortion, streaming, secondary flow, etc.).

Methodology of GKS-BGK
The presented method involves a two-dimension GKS-BGK model. This section gives the main idea of this scheme. Interested readers can refer to References [23,24] for more details about this scheme.
The fully compressible NS equations can be recovered from the Boltzmann equation [23]. From this point, the present GKS is based on the analytical solution of the following two-dimension BGK-Boltzmann equation: where f is the gas distribution function, and g is the equilibrium state function approached by f. Both are functions of time (t), space (x, y), particle velocity (ξ 1 , ξ 2 ), and internal variables (ζ ∈ R K ); τ is the particle collision time, and the function g is supposed to be a Maxwellian distribution.
Here, λ = 1/2RT and ζ 2 = ζ 2 1 + ζ 2 2 + · · · + ζ 2 K relate with the internal energy of the particles. For the two-dimension ideal gas, the total number of degrees of freedom is K = (4 − 2γ)/(γ − 1), where γ is the ratio of specific heats; u and v are the flow velocities in directions x and y, respectively. During the process of particle collisions, the conservation laws are satisfied. Then, one can impose the conservation constraint on functions f and g: where dΞ = dξ 1 dξ 2 dζ, dζ = dζ 1 dζ 2 · · · dζ K and ψ = (ψ 1 , ψ 2 , ψ 3 , ψ 4 ) T = 1, ξ 1 , ξ 2 , ξ 1 2 + ξ 2 2 + ζ 2 /2 T (4) After taking the moment, ψ, in Equation (1) and integrating it with respect to dΞ in the velocity space, dS on a cell S i,j , and dt in a time step [t n , t n+1 ], we obtain the following: where W n i,j are the average macroscopic conservative quantities over S i,j with boundary ∂Ω i,j at time t n . F is the flux across the cell interface, which is related with the function f as follows: where ρε = p/(γ − 1) + ρ u 2 + v 2 /2 is the total energy. Therefore, Equations (5) and (6) mean that it is critical for the GKS (Equation (5)) to obtain the function (f ) at the cell interface.
Appl. Sci. 2020, 10, 1694 4 of 18 The BGK Equation (1) has the analytical solution: Here, and f 0 is the gas distribution function at the beginning of each time step (t = 0). The key point of the GKS-BGK is to construct the functions f 0 and g around the cell interface by Chapman-Enskog expansion. More details about the construction process for functions f 0 and g can be found in References [23,24].
It must be stressed that the inclusion of the discontinuity at a cell interface in functions f 0 and g is a quality of the GKS, which enhances its performance in capturing flow discontinuities, such as shocks.

GKS on the Triangle Mesh
The numerical scheme introduced in Reference [23] can be modified such that it can be implemented on unconstructed grids. The unconstructed mesh must be applied when the computation region is not regular, such as a closed resonator with variable cross-section. Figure 1 shows a schematic of unconstructed triangle mesh.
Appl. Sci. 2020, 10, x 4 of 17 )  1  2  1  2  0   0  1  2  1  2   1  , , , , ,  , , , , ,   , , , , Here, , and f0 is the gas distribution function at the beginning of each time step (t = 0). The key point of the GKS-BGK is to construct the functions f0 and g around the cell interface by Chapman-Enskog expansion. More details about the construction process for functions f0 and g can be found in References [23,24]. It must be stressed that the inclusion of the discontinuity at a cell interface in functions f0 and g is a quality of the GKS, which enhances its performance in capturing flow discontinuities, such as shocks.

GKS on the Triangle Mesh
The numerical scheme introduced in Reference [23] can be modified such that it can be implemented on unconstructed grids. The unconstructed mesh must be applied when the computation region is not regular, such as a closed resonator with variable cross-section. Figure 1 shows a schematic of unconstructed triangle mesh. A crucial step in the application of the GKS on a triangle mesh is to compute the fluxes through the grid sides. Taking side e as an example, the computational method for the fluxes through the side e is as follows. First, the least square method, but not limited to, is applied to reconstruct the gradients for the macroscopic values in the meshes i and i+1. Then, local coordinate system x Oy ′ ′ is constructed at the middle point of side e; the x′ and y′ directions in the local coordinate system are the same as those of vectors n and t, respectively. The relationship of the macroscopic values and their gradients in the two coordinate systems x Oy ′ ′ and xOy are listed as follows: Note that the gradient of other quantities can be computed similarly to Equation (9). The algorithm presented in Reference [23] is implemented in the local coordinate system x Oy ′ ′ to obtain the local fluxes ′ = through the side e. Finally, the fluxes in the local coordinate system x Oy ′ ′ must turn into the fluxes in the reference coordinate system xOy by rotation and translation, as follows: In Figure 1, W i−1 , · · · , W i+2 denote the average macroscopic quantities for the i-1 th , . . . , i+2 th triangle grid, respectively, and e is one of the sides for the i th grid, with unit outer normal vector n = n x , n y T and unit tangent vector t = −n y , n x T . A crucial step in the application of the GKS on a triangle mesh is to compute the fluxes through the grid sides. Taking side e as an example, the computational method for the fluxes through the side e is as follows. First, the least square method, but not limited to, is applied to reconstruct the gradients for the macroscopic values in the meshes i and i+1. Then, local coordinate system x Oy is constructed at the middle point of side e; the x and y directions in the local coordinate system are the same as those of vectors n and t, respectively. The relationship of the macroscopic values and their gradients in the two coordinate systems x Oy and xOy are listed as follows: Note that the gradient of other quantities can be computed similarly to Equation (9). The algorithm presented in Reference [23] is implemented in the local coordinate system x Oy to obtain the local fluxes F = F ρ , F ρ u , F ρ v , F ρ ε T through the side e. Finally, the fluxes in the local Appl. Sci. 2020, 10, 1694 5 of 18 coordinate system x Oy must turn into the fluxes in the reference coordinate system xOy by rotation and translation, as follows: After all fluxes through the sides of the mesh are computed, Equation (5) can be employed to update the averaged macroscopic quantities. This is a crucial procedure for GKS dealing with the variable cross-section problems.

Physical Model and Boundary Conditions
The physical model, shown in Figure 2, is a piston-driving exponential-shaped resonator which is closed rigidly at the small end (left end, x = 0). The vibration of the piston generates the oscillatory flow field and induces standing-waves at the lowest acoustic mode in the resonator. We assume that the cross-sectional area s is a function of x. That is s(x) = s 0 exp(2mx/L), x ∈ [0, L]. Here, m is the shape parameter of the resonator, L is the length, and s 0 is the area of the resonator at the closed end. The velocity of the vibrating piston is given as u(t) = lω sin(ωt) = u 0 sin(θ), where ω is the angular frequency of the piston's vibration.
Appl. Sci. 2020, 10, x 5 of 17 After all fluxes through the sides of the mesh are computed, Equation (5) can be employed to update the averaged macroscopic quantities. This is a crucial procedure for GKS dealing with the variable cross-section problems.

Physical Model and Boundary Conditions
The physical model, shown in Figure 2, is a piston-driving exponential-shaped resonator which is closed rigidly at the small end (left end, x = 0). The vibration of the piston generates the oscillatory flow field and induces standing-waves at the lowest acoustic mode in the resonator. We assume that the cross-sectional area s is a function of x. That is Here, m is the shape parameter of the resonator, L is the length, and s0 is the area of the resonator at the closed end. The velocity of the vibrating piston is given as ( ) where ω is the angular frequency of the piston's vibration. An earlier study [29] describes eight different exponential resonators with m = 0, 0.6, 2.0, 2.1, 2.2, 2.3, 2.4, and 3.0; the authors revealed that the m = 2.2 exponential tube exhibits better sound-pressure amplification ability and suppression effects on shock waves, acoustic streaming, and vortices. Accordingly, in this study, the m = 2.2 exponential resonator was investigated under different acoustic vibration conditions.
The initial state of the computation is unperturbed air, i.e., initial density here, Ω is the resonance angular frequency of the exponential tube and is a function of m [21].
No-slip and adiabatic boundary conditions are applied. To compute the flux through the wall of b, which relates the triangle mesh i, as shown in Figure 3, a ghost triangle mesh g is constructed first.
be the average density, pressure, and velocity over mesh i, respectively, the corresponding values in the ghost cell are computed as follows: ; ; where subscript g indicates the ghost cell. Note that the above equation should be modified at driving end, as follows: An earlier study [29] describes eight different exponential resonators with m = 0, 0.6, 2.0, 2.1, 2.2, 2.3, 2.4, and 3.0; the authors revealed that the m = 2.2 exponential tube exhibits better sound-pressure amplification ability and suppression effects on shock waves, acoustic streaming, and vortices. Accordingly, in this study, the m = 2.2 exponential resonator was investigated under different acoustic vibration conditions.
The initial state of the computation is unperturbed air, i.e., initial density ρ 0 = 1.226 kg/m 3 , pressure P 0 = 101,300 Pa, kinetic viscosity coefficient µ = 1.86 × 10 −5 kg/(m.s), gas constant R = 287 J/(kg.K), ratio of specific heat γ = 1.4, Prandtl number P r = 0.71, and acoustic speed c 0 = 340.1 m/s. Other parameters are set as follows: the radius at the closed end is r 0 = 0.004 m, the length of resonator is L = 0.2 m, and the angular frequency of the piston vibration is ω = Ω = c 0 √ m 2 + π 2 /L; here, Ω is the resonance angular frequency of the exponential tube and is a function of m [21].
No-slip and adiabatic boundary conditions are applied. To compute the flux through the wall of b, which relates the triangle mesh i, as shown in Figure 3, a ghost triangle mesh g is constructed first. Letting ρ i , p i , and U i = (u i , v i ) be the average density, pressure, and velocity over mesh i, respectively, the corresponding values in the ghost cell are computed as follows: where subscript g indicates the ghost cell. Note that the above equation should be modified at driving end, as follows:

Results and Discussion
Three hundred cycles of wave motion are enough to achieve a steady state for Cases A-C. However, even five hundred cycles past, the pressure and velocity at the same moment of different cycles are different for Cases D and E. Therefore, all curves and figures shown in the following for Cases A-C are obtained at a steady periodic state, but at the 500th acoustic cycle for Cases D and E.
All results presented in the following sections were verified to be grid independent, by comparing them with different mesh sizes.

Validation of the Model
To verify the present GKS model and the relevant code, the numerical results from this work were compared with values predicted by using the Galerkin method. Figure 4 shows the curves of the pressure waves during two acoustic periods, calculated from these two methods. The pressure, P, in Figure 4 is P0 + p, where p is the oscillated pressure. The non-dimensional parameter 2 t π Ω expresses the acoustic cycle time under steady state, which is the same as that in Figure 2 of Reference [19]. The solid line denotes the results from the GKS model with shape parameter m = 2.0, and the circular symbol represents the results obtained with the Galerkin method, from Reference [19], with flare constant α = 4.0. Here, a shape parameter, m = 2.0, is corresponding to a flare constant, α = 4.0, and the non-dimensional velocity amplitude of the piston is U* = 5 × 10 −4 . The physical model and parameters settings are the same as in the literature [19], for convenient comparison. A slight discrepancy can be observed between the two waveforms, but overall, they seem to be in good agreement. The possible cause for the discrepancy may come from the different governing equations. The results from the Galerkin method are based on a one-dimensional wave equation, while the present results from the GKS are two-dimensional.

Results and Discussion
Three hundred cycles of wave motion are enough to achieve a steady state for Cases A-C. However, even five hundred cycles past, the pressure and velocity at the same moment of different cycles are different for Cases D and E. Therefore, all curves and figures shown in the following for Cases A-C are obtained at a steady periodic state, but at the 500th acoustic cycle for Cases D and E.
All results presented in the following sections were verified to be grid independent, by comparing them with different mesh sizes.

Validation of the Model
To verify the present GKS model and the relevant code, the numerical results from this work were compared with values predicted by using the Galerkin method. Figure 4 shows the curves of the pressure waves during two acoustic periods, calculated from these two methods. The pressure, P, in Figure 4 is P 0 + p, where p is the oscillated pressure. The non-dimensional parameter tΩ/2π expresses the acoustic cycle time under steady state, which is the same as that in Figure 2 of Reference [19]. The solid line denotes the results from the GKS model with shape parameter m = 2.0, and the circular symbol represents the results obtained with the Galerkin method, from Reference [19], with flare constant α = 4.0. Here, a shape parameter, m = 2.0, is corresponding to a flare constant, α = 4.0, and the non-dimensional velocity amplitude of the piston is U* = 5 × 10 −4 . The physical model and parameters settings are the same as in the literature [19], for convenient comparison. A slight discrepancy can be observed between the two waveforms, but overall, they seem to be in good agreement. The possible cause for the discrepancy may come from the different governing equations. The results from the Galerkin method are based on a one-dimensional wave equation, while the present results from the GKS are two-dimensional.

Results and Discussion
Three hundred cycles of wave motion are enough to achieve a steady state for Cases A-C. However, even five hundred cycles past, the pressure and velocity at the same moment of different cycles are different for Cases D and E. Therefore, all curves and figures shown in the following for Cases A-C are obtained at a steady periodic state, but at the 500th acoustic cycle for Cases D and E.
All results presented in the following sections were verified to be grid independent, by comparing them with different mesh sizes.

Validation of the Model
To verify the present GKS model and the relevant code, the numerical results from this work were compared with values predicted by using the Galerkin method. Figure 4 shows the curves of the pressure waves during two acoustic periods, calculated from these two methods. The pressure, P, in Figure 4 is P0 + p, where p is the oscillated pressure. The non-dimensional parameter 2 t π Ω expresses the acoustic cycle time under steady state, which is the same as that in Figure 2 of Reference [19]. The solid line denotes the results from the GKS model with shape parameter m = 2.0, and the circular symbol represents the results obtained with the Galerkin method, from Reference [19], with flare constant α = 4.0. Here, a shape parameter, m = 2.0, is corresponding to a flare constant, α = 4.0, and the non-dimensional velocity amplitude of the piston is U* = 5 × 10 −4 . The physical model and parameters settings are the same as in the literature [19], for convenient comparison. A slight discrepancy can be observed between the two waveforms, but overall, they seem to be in good agreement. The possible cause for the discrepancy may come from the different governing equations. The results from the Galerkin method are based on a one-dimensional wave equation, while the present results from the GKS are two-dimensional.   Figure 5 shows the pressure fluctuation for eight m exponential ducts (m = 0, 0.6, 2.0, 2.1, 2.2, 2.3, 2.4, and 3.0) near position (x = 0, r(x) = 0), with piston driving amplitude l = 100 µm. Figure 5 shows that the pressure waveform and amplitude are strongly influenced by the shape parameter m. The value m producing highest pressure amplitude can also be estimated from Figure 5. A periodic shock wave appears in the exponential ducts of m = 0 (s L /s 0 = 1, cylindrical tube) and up to m = 0.6 (s L /s 0 = 3). A shock-free and smooth sine wave can be obtained in the exponential duct with m = 2.0 to 3.0. A shock-free and high-amplitude resonance with nearly equal amplitude can be observed with m = 2.1 and 2.2. In other words, the amplitude of the pressure wave increases with the increase of m when m ≤ 2.2 (s L /s 0 ≤ 81), but decreases with the increase of m when m ≥ 2.3 (s L /s 0 ≥ 100). Hossain [21] also reported that shock waves appear in exponential ducts up to s L /s 0 = 2, and the pressure amplitude increases with increasing area contraction ratio when s L /s 0 ≤ 100.

Effect of Shape Parameter
Appl. Sci. 2020, 10, x 7 of 17 Figure 5 shows the pressure fluctuation for eight m exponential ducts (m = 0, 0.6, 2.0, 2.1, 2.2, 2.3, 2.4, and 3.0) near position (x = 0, r(x) = 0), with piston driving amplitude l = 100 μm. Figure 5 shows that the pressure waveform and amplitude are strongly influenced by the shape parameter m. The value m producing highest pressure amplitude can also be estimated from Figure 5. A periodic shock wave appears in the exponential ducts of m = 0 (sL/s0 = 1, cylindrical tube) and up to m = 0.6 (sL/s0 = 3). A shock-free and smooth sine wave can be obtained in the exponential duct with m = 2.0 to 3.0. A shock-free and high-amplitude resonance with nearly equal amplitude can be observed with m = 2.1 and 2.2. In other words, the amplitude of the pressure wave increases with the increase of m when m ≤ 2.2 (sL/s0 ≤ 81), but decreases with the increase of m when m ≥ 2.3 (sL/s0 ≥ 100). Hossain [21] also reported that shock waves appear in exponential ducts up to sL/s0 = 2, and the pressure amplitude increases with increasing area contraction ratio when sL/s0 ≤ 100. As we know, the non-linear effect in propagating shock waves suppresses the growth of the pressure amplitude by the so-called acoustic saturation phenomenon. Therefore, in order to obtain high-intensity sound, the shock wave must be suppressed. The results shown in Figure 5 indicate that some exponential tubes are more suitable for this task. These results also show that, for the exponential tube, the optimal shape parameter was 2.2, with which the highest value of the pressure amplitude can be produced. An earlier study [29] also revealed that the m = 2.2 exponential tube has a better energy-gathering ability and shock-suppression effect.

Effect of Driving Strength
In this section, we study the effect of the driving strength on the non-linear features for the optimal m exponential duct. Five calculation cases with different driving strengths are involved, ranging from the quasi-linear to highly non-linear regions. The values of l are 100, 200, 300, 400, and 500 μm, which correspond to Cases A-E, respectively. The sound and flow fields are also analyzed numerically for the m = 2.2 exponential-shaped resonator, and an optimal driving strength is suggested. Figure 6 shows the evolution of the non-linear pressure and velocity distributions along the symmetry axis (r(x) = 0) of the resonator, at eight different times ( 0, θ = 4, π 2 4, π 3 4, π 4 4, π 5 4, π 6 4, π and 7 4 π ) during one cycle under various driving strengths. One can see from Figure   6 that the pressure amplitudes and their gradients are larger at the closed end of the resonator than As we know, the non-linear effect in propagating shock waves suppresses the growth of the pressure amplitude by the so-called acoustic saturation phenomenon. Therefore, in order to obtain high-intensity sound, the shock wave must be suppressed. The results shown in Figure 5 indicate that some exponential tubes are more suitable for this task. These results also show that, for the exponential tube, the optimal shape parameter was 2.2, with which the highest value of the pressure amplitude can be produced. An earlier study [29] also revealed that the m = 2.2 exponential tube has a better energy-gathering ability and shock-suppression effect.

Effect of Driving Strength
In this section, we study the effect of the driving strength on the non-linear features for the optimal m exponential duct. Five calculation cases with different driving strengths are involved, ranging from the quasi-linear to highly non-linear regions. The values of l are 100, 200, 300, 400, and 500 µm, which correspond to Cases A-E, respectively. The sound and flow fields are also analyzed numerically for the m = 2.2 exponential-shaped resonator, and an optimal driving strength is suggested. Figure 6 shows the evolution of the non-linear pressure and velocity distributions along the symmetry axis (r(x) = 0) of the resonator, at eight different times (θ = 0, π/4, 2π/4, 3π/4, 4π/4, 5π/4, 6π/4, and 7π/4) during one cycle under various driving strengths. One can see from Figure 6 that the pressure amplitudes and their gradients are larger at the closed end of the resonator than at the piston end, indicating that the contraction in the cross-sectional area can increase the pressure magnitude to satisfy the need of high-intensity acoustics [21]. Meanwhile, the profiles of the spatial distribution of pressure and velocity waves are dramatically changed with increasing driving strength. For example, the waveforms change from smooth waves at all instants in the quasi-linear case (Figure 6a) to non-smooth and distorted forms at some instants in Case B and Case C (Figure 6b,c, respectively); non-smooth and distorted waves are also observed at all instants in Figure 6d,e. However, no shock wave appears in the duct, even at high driving strengths, which is very different from the case in cylindrical duct [27].

Sound Fields
Appl. Sci. 2020, 10, x 8 of 17 at the piston end, indicating that the contraction in the cross-sectional area can increase the pressure magnitude to satisfy the need of high-intensity acoustics [21]. Meanwhile, the profiles of the spatial distribution of pressure and velocity waves are dramatically changed with increasing driving strength. For example, the waveforms change from smooth waves at all instants in the quasi-linear case (Figure 6a) to non-smooth and distorted forms at some instants in Case B and Case C ( Figure  6b,c, respectively); non-smooth and distorted waves are also observed at all instants in Figure 6d,e. However, no shock wave appears in the duct, even at high driving strengths, which is very different from the case in cylindrical duct [27].  Figure 7 shows the variation of temperature near the position (x = 0, r(x) = 0) over four standing wave cycles from Cases A-E for the m = 2.2 exponential tube. Several changes are shown in Figure 7 for temperature waveform and its amplitude and phase as the driving strength increases. From Cases A-C, little changes in the waveforms and amplitudes are observed, and the phases shift slowly. However, drastic changes in waveforms and amplitudes occur from Cases C and D and, in particular, Case E. Narrower peaks and broader troughs can be seen in Case D and Case E; distortion waves first appear at the troughs in Cases D and E that then affect the peaks in Case E. These phenomena imply a rapid growth of energy dissipation for Cases D and E. However, no "discontinuity wave" or "shock wave" or their analogs can be observed in Figure 7, even at large driving amplitudes.   Figure 7 for temperature waveform and its amplitude and phase as the driving strength increases. From Cases A-C, little changes in the waveforms and amplitudes are observed, and the phases shift slowly. However, drastic changes in waveforms and amplitudes occur from Cases C and D and, in particular, Case E. Narrower peaks and broader troughs can be seen in Case D and Case E; distortion waves first appear at the troughs in Cases D and E that then affect the peaks in Case E. These phenomena imply a rapid growth of energy dissipation for Cases D and E. However, no "discontinuity wave" or "shock wave" or their analogs can be observed in Figure 7, even at large driving amplitudes. = are more complicated than those at the other two moments in one acoustic cycle; a vortex pair is formed on either side of the tube where it is enlarged and contracted. It appears that the vortices generate and vanish periodically. Vortices are accompanied by a larger temperature gradient and a smaller velocity gradient at 0 θ = and θ π = , but a nearly uniform temperature field and greater velocity fluctuation, as well as smoother streamline, can been observed at 2 θ π = and θ = 3 2 π . The streaming that was superposed on the main flow developed owing to the generation, evolution, and shedding of vortices, which is referred to as the third type of direct current (DC) flow in Reference [12]. The DC flow is driven by the asymmetry of the fluid flow and temperature and transmitted by the vortices, and it accompanies energy storing and releasing in a cycle.
The flow regime for Case C shown in Figure 9 is different from that for Case A in Figure 8. For example, the flow fields in Figure 9 have a greater velocity gradient and smoother streamlines at 0 θ = and θ π = , and more varied temperature fields at 2 θ π = and 3 2 θ π = , which indicates a greater minor loss and energy dissipation; the vortex only appears at the large end of the resonator at certain instants ( 0 θ = ) and is accompanied by the biggest temperature variation in one cycle, because the pressure reaches its maximum at this moment in one cycle. The vortex shedding corresponds to an increase in the difference of pressure amplitude between the accelerating and decelerating phases, so that the cycle symmetry of the flow field is broken. The asymmetrical and non-smooth pressure and velocity fluctuation (see Figure 6) in the resonator is one of the causes of the scale and shape of the vortices, and the vortex propagation speed contributes to the position of the vortices, so the streaming flow and the process of energy storage and release in Figure 9 appear different from those in Figure 8.  Figures 8 and 9 show the transient flow fields in terms of streamline at typical moments (θ = 0, π/2, π, and 3π/2.) during one acoustic cycle in the m = 2.2 exponential tube for Cases A and C, respectively. The color contours indicate the distribution of the temperature (K), and the reference velocity is the u velocity (m/s). Considering the symmetry, only the top half part of the flow fields is depicted here. The evolution processes of the flow and streaming for different driving strengths are very different, which is the most important and interesting point. The flow direction and phase shown in these figures change with increasing driving amplitude, because the fundamental modal frequency of the exponential resonator increases with increasing driving strength [10].

Transient Flow Fields
The figures in Figure 8 indicate that the flow fields for Case A at θ = 0 and θ = π are more complicated than those at the other two moments in one acoustic cycle; a vortex pair is formed on either side of the tube where it is enlarged and contracted. It appears that the vortices generate and vanish periodically. Vortices are accompanied by a larger temperature gradient and a smaller velocity gradient at θ = 0 and θ = π, but a nearly uniform temperature field and greater velocity fluctuation, as well as smoother streamline, can been observed at θ = π/2 and θ = 3π/2. The streaming that was superposed on the main flow developed owing to the generation, evolution, and shedding of vortices, which is referred to as the third type of direct current (DC) flow in Reference [12]. The DC flow is driven by the asymmetry of the fluid flow and temperature and transmitted by the vortices, and it accompanies energy storing and releasing in a cycle.
The flow regime for Case C shown in Figure 9 is different from that for Case A in Figure 8. For example, the flow fields in Figure 9 have a greater velocity gradient and smoother streamlines at θ = 0 and θ = π, and more varied temperature fields at θ = π/2 and θ = 3π/2, which indicates a greater minor loss and energy dissipation; the vortex only appears at the large end of the resonator at certain instants (θ = 0) and is accompanied by the biggest temperature variation in one cycle, because the pressure reaches its maximum at this moment in one cycle. The vortex shedding corresponds to an increase in the difference of pressure amplitude between the accelerating and decelerating phases, so that the cycle symmetry of the flow field is broken. The asymmetrical and non-smooth pressure and velocity fluctuation (see Figure 6) in the resonator is one of the causes of the scale and shape of the vortices, and the vortex propagation speed contributes to the position of the vortices, so the streaming flow and the process of energy storage and release in Figure 9 appear different from those in Figure 8.  Figures 8 and 9 indicate that the DC flow is very sensitive to the intensity of the flow field because the pressure differential, wave amplitude, and the corresponding vortex propagation speed are different. Here, a detailed analysis is given for Case E, with large driving amplitudes. The instantaneous velocity distributions at eight different instants (θ = 0, π/4, π/2, 3π/4, π, 5π/4, 3π/2, and 7π/4.) during one cycle, along with the streamlines, are shown in Figure 10; the whole flow fields are depicted in the figure for better flow visualization.
The flow field and streaming structure in Figure 10 are very different from those in Figures 8 and 9. A very strong pressure differential occurred between positive and negative half periods that broke the vortex, distorted the flow direction, and generated turbulence. Thus, no vortices can be observed in Figure 10, and the process of energy storage and release is also unclear. Meanwhile, the changes in the velocity and temperature from one moment to another are gentler than those in the previous cases, implying that the bigger the harmonic distortion is (see Figures 6 and 7), the more homogeneous the distributions of the flow and temperature fields inside the resonator are [32]. Moreover, chaos and irregularity of the streaming overlaid on the main flows can speed up the uniformity distribution of the flow field in space and time. The turbulent motion also causes the time-averaged pressure, velocity, and temperature to fluctuate irregularly [27], which will consume the energy and prevent its storage and release.
Appl. Sci. 2020, 10, x 12 of 17 θ π = The color contours indicate the temperature distributions (K); the reference velocity is the u velocity (m/s). Figures 8 and 9 indicate that the DC flow is very sensitive to the intensity of the flow field because the pressure differential, wave amplitude, and the corresponding vortex propagation speed are different. Here, a detailed analysis is given for Case E, with large driving amplitudes. The instantaneous velocity distributions at eight different instants ( 0 θ = , 4, π 2, π 3 4, π , π 5 4, π 3 2, π and 7 4. π ) during one cycle, along with the streamlines, are shown in Figure 10; the whole flow fields are depicted in the figure for better flow visualization. The flow field and streaming structure in Figure 10 are very different from those in Figures 8  and 9. A very strong pressure differential occurred between positive and negative half periods that broke the vortex, distorted the flow direction, and generated turbulence. Thus, no vortices can be observed in Figure 10, and the process of energy storage and release is also unclear. Meanwhile, the changes in the velocity and temperature from one moment to another are gentler than those in the previous cases, implying that the bigger the harmonic distortion is (see Figures 6 and 7), the more homogeneous the distributions of the flow and temperature fields inside the resonator are [32]. Moreover, chaos and irregularity of the streaming overlaid on the main flows can speed up the uniformity distribution of the flow field in space and time. The turbulent motion also causes the time-averaged pressure, velocity, and temperature to fluctuate irregularly [27], which will consume the energy and prevent its storage and release. The numerical simulation of velocity field with turbulent motion for high-amplitude cases is also reported in References [27,33] for cylindrical resonators, and the critical Reynolds numbers for the transition to turbulence reported are 374 and 400, respectively. In this work, however, the Reynolds number for the m = 2.2 exponential resonator in Case E is 177, which is smaller than that for cylindrical resonators. This may be because the streaming in the exponential resonator is hydrodynamically asymmetric.

Mean Flow Fields
The streaming velocity can be calculated based on the average mass transport velocities [27]: where, ust and vst are the axial and radial components of the streaming velocity, respectively, and the symbol ⋅ means time averaging. The secondary flow (mean flow), which superposes on the resonant oscillation and forms the fluid motion, can be described by the streaming velocities ust and vst.
The mean flow fields for the m = 2.2 exponential resonator under five cases are shown in Figure  11, where the color contours indicate the distribution of the axial streaming velocity ust (m/s). Only the upper half of the flow fields is depicted because the mean flow fields are all nearly axisymmetric. We can find from Figure 11 that, along with increasing driving strength, the "vortex fission" phenomenon occurs near the small end of the resonator, but little changes can be found for the vortex structure near the driving end from Cases A-D. However, the big vortex ring near the large The numerical simulation of velocity field with turbulent motion for high-amplitude cases is also reported in References [27,33] for cylindrical resonators, and the critical Reynolds numbers for the transition to turbulence reported are 374 and 400, respectively. In this work, however, the Reynolds number for the m = 2.2 exponential resonator in Case E is 177, which is smaller than that for cylindrical resonators. This may be because the streaming in the exponential resonator is hydrodynamically asymmetric.

Mean Flow Fields
The streaming velocity can be calculated based on the average mass transport velocities [27]: where, u st and v st are the axial and radial components of the streaming velocity, respectively, and the symbol · means time averaging. The secondary flow (mean flow), which superposes on the resonant oscillation and forms the fluid motion, can be described by the streaming velocities u st and v st . The mean flow fields for the m = 2.2 exponential resonator under five cases are shown in Figure 11, where the color contours indicate the distribution of the axial streaming velocity u st (m/s). Only the upper half of the flow fields is depicted because the mean flow fields are all nearly axisymmetric. We can find from Figure 11 that, along with increasing driving strength, the "vortex fission" phenomenon occurs near the small end of the resonator, but little changes can be found for the vortex structure near the driving end from Cases A-D. However, the big vortex ring near the large end also split into two vortices for Case E, so that the streaming rolls in Figure 11e spread to the middle region from the two ends and divide into secondary vortices of various scales. end also split into two vortices for Case E, so that the streaming rolls in Figure 11e spread to the middle region from the two ends and divide into secondary vortices of various scales.
In brief, Figure 11 indicates that vortices only appear at the two ends of the resonator in Figure  11a-d, but irregular and unsteady vortices of various scales occur throughout the resonator in Figure 11e. An earlier study [32] -argued that the refrigeration performance can be enhanced by a so-called bi-cellular structure, which looks like that shown in Figure 11a-d. He et al. [22] also observed a better performance of the tube when the flow was less uniformly distributed. Therefore, to avoid detrimental effects of secondary streaming, an optimal driving strength must be found. All information indicates that this optimal driving strength is l ≤ 400 μm.  Table 1 lists the increments of ΔP, |umax|, and ΔT and their growth ratios with each increment of l = 100 μm. Here,  In brief, Figure 11 indicates that vortices only appear at the two ends of the resonator in Figure 11a-d, but irregular and unsteady vortices of various scales occur throughout the resonator in Figure 11e. An earlier study [32]argued that the refrigeration performance can be enhanced by a so-called bi-cellular structure, which looks like that shown in Figure 11a-d. He et al. [22] also observed a better performance of the tube when the flow was less uniformly distributed. Therefore, to avoid detrimental effects of secondary streaming, an optimal driving strength must be found. All information indicates that this optimal driving strength is l ≤ 400 µm. Table 1 lists the increments of ∆P, |u max |, and ∆T and their growth ratios with each increment of l = 100 µm. Here, ∆P = P max − P min and ∆T = T max − T min are the peak-to-peak pressure and temperature amplitude, respectively, and |u max | denotes the absolute value of the maximum axial flow velocities. Although the values of ∆P, ∆T, and |u max | all increase with increasing driving strength, their growth ratios do not always increase (some values are negative). From Table 1, it can be easily found that there is a sudden decrease for the increments of ∆P and |u max | but a sudden increase for the increments of ∆T under the condition of Case D to Case E, although they keep increasing from Cases A to Case D. Therefore, ∆P and |u max | have maximal increments when l increases from 300 to 400 µm and minimal increments when l increases from 400 to 500 µm. Moreover, the biggest growth ratios for ∆T occur when ∆P and |u max | have the smallest growth ratios, implying a more rapid increase of energy dissipation under the condition change from Case D to Case E than other condition changes.

Optimal Driving Strength
This suggests that the optimal exponential resonator has the biggest ratio of output work to input work when the driving displacement amplitude is between 300 and 400 µm. Furthermore, if the driving amplitude is larger than 400 µm, the sound pressure amplification ability will become significantly weak, owing to the acceleration of energy dissipation.

Conclusions
In this study, the flow and streaming motions in an exponential-type resonator were studied, using the gas-kinetic BGK model. The geometry and driving strength were found to be the main factors affecting the streaming patterns and the performance of the system. The observations are as follows: (1) An optimal shape parameter, i.e., m = 2.2, for the exponential-shaped resonator, can be used to obtain the largest oscillatory pressure amplitude. Therefore, some numerical experiments are carried out for an m = 2.2 exponential-shaped resonator under five different driving strengths, which lead the acoustic fields to change from quasi-linear to highly non-linear. (2) Vortex formation and evolution are closely linked to the waveform of pressure, velocity, and temperature. The third type of DC streaming is found under quasi-linear and moderate non-linear conditions. Driven by the asymmetry of fluid flow and temperature, the DC flow is formed as a result of the generation, evolution, and shedding of vortices, and it accompanies a periodic energy storage and release. (3) The flow and streaming patterns are found to be very sensitive to the driving strength. With increasing driving strength, a very large pressure and temperature asymmetry will break the vortex and disorder the flow direction and waveforms. In particular, a turbulent flow with irregular and unsteady direction and secondary flow may occur if the driving amplitude increases beyond a critical value. (4) The largest growth ratio for peak-to-peak pressure amplitude and axial flow velocity occurs when the driving strength varies from 300 to 400 µm, while the smallest occurs when the driving strength varies from 400 to 500 µm. Therefore, to obtain a larger output power and avoid