Abstract
This paper describes a numerical model based on discontinuous Galerkin method for thermoacoustic investigation. Numerical investigation was conducted to study the behaviour of thermoacoustic wave propagations induced by thermal effects in 2-dimensional enclosure. The compressible Navier-Stokes equations are used as the governing equations. The spatial domain was discretized by using unstructured discontinuous Galerkin method and the explicit fourth-order Runge-Kutta was used to integrate the temporal domain. The accuracy of the numerical results was assessed against other numerical methods and linear solutions. The comparisons with linear solutions show excellent agreement and comparisons with flux corrected transport (FCT) method show fair agreement.
1. Introduction
In addition to being commonly caused by mechanical vibration, the initiation of acoustic waves can also be triggered by sudden temperature change. When a compressible fluid experiences rapid heating, the fluid will expand. This expansion results in a local pressure disturbance which then propagates. Due to viscous and thermal effects, the amplitude of the thermoacoustic induced wave will eventually be dampened. The oscillating motion of the particles due to wave propagation also results in convective heat transfer, thus resulting in heat transfer enhancement [1]. To better understand the mechanism of this phenomenon, numerical studies have been conducted to investigate the thermoacoustic wave.
Previous numerical studies with regard to thermoacoustics have mostly adopted control volume and finite volume methods. An early study was conducted by Worlikar et al. using a finite difference low Mach-number model for simulation of unsteady adiabatic and stratified flow in a thermoacoustic stack [2, 3]. By using the computed results, they analyzed the configuration stacks and the nonlinear response of the flow due to different acoustic driving amplitudes and frequencies. Their results also showed that the computed results have good agreement with experimental results and linear theory. Hantschk and Vortmeyer [4] used the fluent finite volume model to simulate a Rijke tube. In this model, the thermoacoustic wave was generated with one heated wire screen. They showed that the obtained results are in good agreement with experiments. Farouk et al. [5, 6] studied the behaviour of thermoacoustic waves in a nitrogen-filled two-dimensional cavity. The thermoacoustic waves were generated by heating or cooling the vertical walls of the cavity. The wall temperatures were altered both impulsively and gradually. A fully correct transport (FCT) method was used to solve the full compressible Navier-Stokes equation. Flow patterns were shown to strongly correlate to the rapidity of the wall heating process. Lin [1] expanded the work of Farouk et al. [5] investigating acoustic wave induced convection and transport in gases under normal and microgravity conditions. He concluded that the acoustic streaming can enhance the heat transfer and the temperature along a stack and cooling effect can be predicted. Nijeholt et al. [7] used the CFX finite volume model to study a traveling-wave thermoacoustic engine. They concluded that CFD codes could be used in the future to predict and optimize thermoacoustic systems. Ke et al. [8] carried out 2-dimensional numerical simulations of a thermoacoustic refrigerator driven at large amplitude. Their computation was based on a pressure-correction (SIMPLE) algorithm for compressible flows. It was revealed that the optimal length of heat exchanger should be proportional to the amplitude of gas in it. In addition, the optimal heat exchanger length should be close to the peak-to-peak displacement amplitude of the gas medium. Lourier et al. [9] used a semi-implicit pressure-based fractional step method and characteristic boundary conditions (NSCBC). They showed that the acoustic CFL limitation can be removed by using implicit NSBC.
The discontinuous Galerkin method is a combination of the finite element method with the finite volume method. As such, this method can handle a complex spatial domain with a good level of accuracy. In 1972, Reed and Hill [10] introduced a DG method to solve the neutron transport equation. Bassi and Rebay developed the DG method for compressible gas dynamics in which the DG method was used to solve the nonlinear Euler equation [11] and then extended to solve the compressible Navier-Stokes equation [12]. Their study showed that the DG method had a high level of accuracy and was robust for supersonic flows. Liu et al. [13] developed a reconstructed DG method for a solution of compressible Navier-Stokes equations on 3D hybrid grids. Viscous and inviscid flux used formulations developed by Bassi and Rebay [11, 12]. The study was capable of increasing the accuracy of the DG method by one order higher than the existing DG methods for compressible turbulent flows with low Mach numbers and high Reynolds numbers.
However, studies on thermoacoustic waves using the discontinuous Galerkin (DG) method are very scarce in literature. Gineste simulated a one-dimensional thermoacoustic system by using the discontinuous Galerkin [DG] method with linearized Euler equations [14]. His results showed that the DG method is effective and quite flexible to be applied to thermoacoustic problems. Based on the study carried out by Gineste, we have extended the solution approach to a two-dimensional problem by using full compressible Navier-Stokes equation as the governing equations.
2. Mathematical Model and Numerical Scheme
The thermoacoustic wave propagation is modelled by the compressible Navier-Stokes equations. These equations can be expressed in a vector conservation form as follows:where is the vector of the conservative variables, are the nonlinear convective fluxes, and are the diffusive fluxes.where is the density of the fluid, and are the components of the momentum, is the pressure, is the total energy, is thermal conductivity, is dynamic viscosity, and is the stress tensor. The stress tensor is defined byIt is assumed that the fluid is ideal gas and the properties are constant.where is specific gas constant.
The diffusive fluxes contain second-order derivatives which are then converted to a system of first-order equations by adding auxiliary variables as conducted by Bassi and Rebay [12]. Thus, (1) can be written as [15]Auxiliary variables have 2 components, that is,By adding convective fluxes and diffusive fluxes into one partthen (5) becomes simpler:Because the discontinuous Galerkin methods are a combination of finite element methods and finite volume methods, the numerical steps are similar to those two methods. The 2-dimensional physical domain is composed of a number of triangular elements.The numerical solution in each element is approximated by nodal high order polynomial basis function .where are spatial coordinates, is the number of nodal points in each triangular element, and is the order of polynomial basis function.
The DG methods allow the solution to be noncontinuous across the element interfaces. The surface integral at the element interfaces, which is known as the numerical flux, is not neglected as it is in finite volume methods. The weak formulation can be obtained by applying the Galerkin integration by parts of the flux terms of the mathematical model as follows:The local Lax-Friedrich flux is used to approximate the inviscid fluxes and the central flux without penalty for diffusive fluxes [16]. Here denotes inner product, for example, .
The integration of the temporal domain was conducted using an explicit fourth-order Runge-Kutta [17]. We have implemented the numerical scheme into computer program based on the NUDG framework (https://github.com/tcew/nodal-dg) of which the details are described in [16].
3. Results and Discussion
3.1. Case 1: Point Heating Induced Wave
Case 1 investigates the effect of a temperature disturbance in initiating pressure wave propagation. The dimensions of domain are 5 m × 5 m (Figure 1) which consisted of 650 triangular elements. The gas medium was CO2 and the fields were initialized with
As can be seen, the initial temperature of the domain was set at 310 K. At , the temperature at the center of the domain was suddenly increased using special Gaussian distribution with a peak of 316.2 K. The sudden rise in temperature resulted in a pressure disturbance that propagated from the center of the domain.
Figures 2(a), 2(b), and 2(c) show contour plots for the pressure perturbations at ms, 3 ms, and 5 ms. The wave propagates radially, so the region at ms is larger than that at ms. As can be seen, due to the expanding front of the pressure wave, the amplitude of the wave reduces as the wave travels further from the center of the domain. The wave propagation velocity was revealed to be around 270 m/s which complies to the speed of sound for CO2 at atmospheric pressure and a temperature of 300 K.
(a) Pressure at ms
(b) Pressure at ms
(c) Pressure at ms
Figure 3 shows the evolution of the temperature along the mid-horizontal plane. At , the temperature at the center of the domain was raised to 316 K. Afterwards, the temperature reduces due to temperature dissipation. As can be seen, there is a slight observance of temperature disturbance propagation from the domain center.
To validate the results of the model, we have compared the simulation results to the results obtained using a linearized Euler finite difference time domain (FDTD) model [18]. We have compiled results using a second-order and fourth-order ( and ) discontinuous Galerkin method in solving the full compressible Navier-Stokes equation. Figure 4 shows the comparison of the results for pressure distribution along the mid-horizontal plane at ms, while Figure 5 shows the comparison for pressure history at the center of the domain. Since the pressure and temperature perturbation that occurs is relatively small, the wave propagation tends to be linear. This is confirmed by the match of the FDTD results and the DG model.
3.2. Case 2: Thermoacoustic Wave Propagation in a 2-Dimensional Enclosure
Having established the validity of the discontinuous Galerkin model, for the second case we proceed to simulate a thermoacoustically induced wave caused by the heating of a solid wall. Again, the gas medium used is CO2 at an initial pressure of 0.1 MPa and temperature of 310 K. Here we will compare a similar case by Lin [1] who used a finite volume method with flux corrected transport (FCT). The domain was set to have an area of 1 mm2 and consisted of 6082 triangular elements (Figure 6).
The upper and lower walls are insulated and a no slip boundary condition was used for all solid walls. The temperature of the right wall is kept constant at 310 K. For Case 2, the wall temperature at the left wall in which the wall temperature increases using a time dependent boundary condition is as follows:where is the time constant of the wall heating. First, we conducted a trial in which the left wall experiences a step change. Here we set the value of . Thus, at , the temperature increases suddenly from 310 K to 413 K.
Figures 7(a) and 7(b) show the pressure at s and 3.2 s for . Heating on the left wall results in the expansion of the fluid and a local increase of pressure at the left wall and eventually a pressure wave which propagates to the right. The sudden increase of temperature results in a sharp and well-defined pressure wave peak. The wave propagation can be seen clearly along with the attenuation of the pressure wave due to viscous effects.
(a)
(b)
When the pressure wave hits the right wall, it is reflected and its amplitude increases as can be seen in Figures 8(a) and 8(b). This occurs again and again as the wave hits the opposing walls multiple times until the wave is fully attenuated. From Figures 7(a), 7(b), 8(a), and 8(b), we can also observe traces circular fronts of weak secondary waves caused by the interaction of the main wave with the upper and lower walls.
(a)
(b)
In addition to the step temperature change on the left wall, we also conducted trials for gradual heating cases. We used , , and for gradual heating cases, where s.
Figures 9(a) and 9(b) show the pressure wave front at s and 3.2 s for . The pressure waves exhibit much lower amplitude than the pressure wave for (step change in temperature). However, the wave length is much longer.
(a)
(b)
Figures 10(a) and 10(b) show the pressure for after the wave is reflected on the right wall. Due to the smaller range of the pressure, it can be clearly seen that the minimum pressure increases. This was also present for albeit less apparent due to the larger peak of the pressure wave. The increase in minimum pressure is most likely due to the increase of chamber temperature caused by conduction from the heated wall.
(a)
(b)
Figure 11 compares the pressure history at the center of the domain for the values tested. As expected, the larger the heating time constant, the smaller the pressure wave peak. The increase of pressure due to the increase of the temperature chamber can be seen clearly here.
To examine the calculation convergence, results for polynomial orders , , and were compared. Figure 12 shows that the larger the value, the better the convergence of the solution.
Comparing the DG method to the finite volume method with FCT developed by Lin [1], we can see that the DG method has a much sharper wave front compared to FCT (Figure 13). This may be caused by the larger dissipative effect present in FCT, thus resulting in a lower wave peak due to this dissipation.
4. Conclusion
We have developed a model to investigate the initiation of thermoacoustic phenomena using the discontinuous Galerkin method. We have run the model to simulate two cases first of which was to validate the model with regard to the linearity when the perturbations are small. Comparisons with the numerical results of linearized Euler equations solved by the FDTD method have shown excellent agreement. For the second case, we have modelled the thermoacoustic wave initiated by the heating of a solid wall of a small chamber. The heating parameters were varied for a range of heating times from a step change in wall temperature to gradual heating. The results have shown that, for step change in wall temperature, the thermoacoustic effect induces a sharp wavefront with a large pressure amplitude, while gradual heating results in much lower pressure amplitudes. We have compared the results to a similar case solved by a finite volume FCT method. The comparison revealed that the DG method resulted in sharper wavefronts compared to the FCT method due to the presence of dissipative effects in the FCT method, thus demonstrating the better accuracy of the DG method. All in all, the DG method is a promising and suitable method for simulation of thermoacoustic phenomena. The higher order of the DG method compared to finite volume methods also allows a greater degree of accuracy.
Conflicts of Interest
The authors declare that they have no conflicts of interest.