Two FORTRAN codes were developed based on the semi-implicit and nearly-implicit methods in order to analyze their performances in solving transient compressible flow. Pressure transient involves the strong coupling of pressure, density, velocity, and temperature, and it is used to testify the performances of the two methods. The transient condition is a step decrease of outlet pressure of a pipe from 2.7 MPa to 2.6 MPa. Other parameters, such as boundary conditions, time step, and space step, are shown in Table 1. The calculated pressure and temperature responses by the two methods are shown in Figures 2, 3, respectively.

It can be seen from Figure 2A that the semi-implicit method is capable of capturing the fast pressure responses with little numerical diffusion. However, it has poor stability in capturing the slow temperature responses, as seen in Figure 2B that temperature diverged at 0.02 ~ 0.04 s. It can be seen from Figure 3B that the nearly-implicit method has good stability in capturing the slow temperature responses. However, it cannot capture the fast pressure responses very accurately, as seen in Figure 3A that the pressure oscillations are suppressed and vanish very quickly.

One possible reason for the great numerical diffusion of the nearly-implicit method is that the time step Δτ = 10⁻⁵ s is too large to capture the fast pressure responses accurately. The variation of time step from 10⁻⁶ s to 10⁻⁸ s was studied, and some of the results are shown in Figure 4. It was found that if time step Δτ ≤ 10⁻⁷ s, the numerical diffusion of pressure can be very small. However, there are two weaknesses: (1) spurious oscillations exist at the pressure step change points; and (2) the time step is too small, which means a great computational effort.

As shown in Figures 2~4, neither the semi-implicit method nor nearly-implicit method provides a stable and accurate solution in solving transient compressible flow.

For the semi-implicit method, velocity is implicit while density is explicit in continuity equation. In energy equation, only velocity is implicit, the other parameters, such as density, pressure, and temperature, are all explicit. In momentum equation, the pressure and frictional terms are implicit, while the convective term is explicit. So the expressions of density, temperature, velocity, and pressure can be derived directly. In calculation, firstly, pressure, density and temperature are solved explicitly from the derived expressions, and then velocity is solved based on the solved pressure.

In the nearly-implicit method, the expression of pressure is derived from continuity and energy equations in the same way as in the semi-implicit method. Different from the semi-implicit method, the convective term in the momentum equation is treated implicitly. In calculation, velocity is solved first, then pressure is solved subsequently based on the solved velocity, and then density and temperature are solved implicitly by updating the solved pressure, velocity and density in continuity and energy equations.

There are mainly two differences between the semi-implicit and nearly-implicit methods: (1) the explicit or implicit treatment of the convective term in the momentum equation, which determines the different expressions of velocity, and (2) the explicit or implicit calculation of density and temperature. According to the two differences and the results in Figures 2–4, it is assumed that (1) in order to capture a detailed fast pressure transient response, the convective term in momentum equation should be treated explicitly, and (2) in order to get a stable solution for the slow temperature transients, coefficients in the continuity, and energy equations should be updated immediately using the previously solved parameters, such as pressure, velocity and density.

To testify this hypothesis, the convective term in momentum equation is treated implicitly as in the nearly-implicit method, and density and temperature are solved explicitly as in the semi-implicit method. If this hypothesis was correct, it is predicted that neither an accurate pressure response nor a stable temperature solution could be obtained. The results are shown in Figure 5, in which the time step is Δτ = 10⁻⁵ s. It can be seen that for the fast pressure response, the artificial diffusion is too big, while for the slow temperature response, the solution diverged. The results are the same as predicted, so the hypothesis is reasonable.

According to the hypothesis discussed above, a new HSI method is proposed. In the HSI method, the convective term in momentum equation is treated explicitly, and density and temperature are solved implicitly by updating coefficients in the continuity and energy equations using the previously solved parameters. The calculation strategy of the HSI method is in the following steps:

Under the same outlet pressure transient in section Analysis of Semi-Implicit and Nearly-Implicit Methods in Compressible Flow, and using parameter settings in Table 1, the pressure and temperature responses of the HSI method are shown in Figure 6. It can be seen from Figures 6A,B that the fast transient responses caused by pressure oscillations are captured and the slow temperature responses are stable and converged. As a result, the newly proposed HSI method has the advantages of both the semi-implicit and nearly-implicit methods but none of their weaknesses.

The shock tube problem is a typical Riemann problem, and has an exact solution (Toro, 2009). In the shock tube problem, a closed tube with a length of 100 m was separated into two equal parts by a membrane film, as shown in Figure 7. The left part of the tube is helium with a pressure of 2.0 MPa and a temperature of 400 K, and the right part of the tube is helium with a pressure of 1.0 MPa and a temperature of 400 K. If the film suddenly breaks in the center, the helium in two parts mixes together, and shock wave, rarefaction wave, and discontinuity wave are produced in this process.

There is no exact closed-form solution to the Riemann problem for the Euler equations of shock-tube problem. However, it is possible to devise iterative schemes whereby the solution can be computed numerically to high accuracy. The space step is set as 0.05 m in both the iterative scheme and the HSI method. Comparison of the exact solution with that of the HSI method is shown in Figure 8. Three time steps of 10⁻⁴ s, 10⁻⁵ s, and 10⁻⁶ s are calculated using the HSI method. It was found that when Δτ = 10⁻⁴ s the solution is divergent, and when Δτ = 10⁻⁵ s and Δτ = 10⁻⁶ s, the results of the HSI method are almost identical, as seen in Figure 8. So the time step is set as 10⁻⁵ s in this case to get a convergent solution with shorter computational time in the following discussion.

It can be seen from Figure 8 that the HSI method can capture the pressure transient with good accuracy. The slope of each wave nearly equals to the exact solution and there are no spurious oscillations at pressure step change points.

In this section, the HSI method was compared with FLUENT simulations in which a two-dimensional laminar flow model was adopted. In the FLUENT simulations, constant outlet pressure, inlet mass flow, and inlet temperature were adopted for the boundary conditions of compressible flow. In order to test the performances of the HSI method under different transient scenarios, three typical transient conditions are studied and compared with FLUENT simulations, namely an outlet pressure transient, inlet mass flow transient, and inlet temperature transient. Among the three transients, the pressure transient is a fast transient with a very short characteristic time that involves significant changes of pressure, velocity, and temperature. The mass flow transient is also a fast transient but does not involve significant changes of pressure or temperature. The temperature transient is a slow transient with a long characteristic time, and causes mild changes of pressure, velocity, and temperature.

The maximum relative density difference is

In which ρ₁ is solved implicitly using the continuity equation and ρ₂ is solved by state equation. Δρ_max is the maximum relative density difference between ρ₁ and ρ₂. If Δρ_max ≤ Δρ_L (the lower limit of Δρ_max), then the errors under this time step are small enough to allow greater time step for faster calculation while still retaining high accuracy. On the other hand, if Δρ_max > Δρ_H (Δρ_H is the upper limit of Δρ_max), then the relative density error may be too large to converge, and the time step should decrease to ensure convergence. Based on the a lot of testing, the recommended values of Δρ_L and Δρ_H here are set as Δρ_L = 10⁻⁸ and Δρ_H = 10⁻³ for the HIS method for the transients of the HTRGTS.

As mentioned in section Maximum Relative Density Difference, a suitable value of Δρ_L saves computational time, and that of Δρ_H ensures convergence. However, the temperature response at x = 0.75 m, as shown in Figure 14, indicates that the solution is unstable under the outlet pressure transient condition.

To ensure stability, the CFL condition must be considered to control the time step below a certain value, so that the Courant number is lower than its maximum value (c_max). Under a fixed mass velocity of 38.90 kg/(m²·s), the relationship between c_max and the number of nodes under different pipe lengths are shown in Figure 15. It can be seen that (1) the value of c_max decreases with the increase of node number, and (2) c_max was not influenced by pipe length. If the node number is >300, the value of c_max ranges between 1.2 and 1.3. To ensure stability, the value of c_max is set as 1.1 for different node numbers. After the CFL condition was considered, the relationship between time step and time is shown in Figure 16A, and the temperature response at x=0.75 m is shown in Figure 16B. It can be seen from Figure 16A that time step first increased to a certain value and then stabilizes, and a stable solution of temperature is obtained, as seen in Figure 16B.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.