3.1 Base model

The baseline wind tunnel model follows the work of Fung et al. [Reference Fung, Settles and Ray1]. We describe the blow-down wind tunnel with two control volumes, one for the reservoir and one for the plenum chamber (see Fig. 1). We then use basic conservation laws to describe the system. Using conservation of mass and energy at both control volumes, we obtain four differential equations:

(1) \begin{align}\frac{{d{M_r}}}{{dt}} = - {\dot m_v}\end{align}

(2) \begin{align}\frac{{{M_p}}}{{dt}} = {\dot m_v} - {\dot m_n}\end{align}

(3) \begin{align}\frac{{d{E_r}}}{{dt}} = - {\dot H_v}\end{align}

(4) \begin{align}\frac{{d{E_p}}}{{dt}} = {\dot H_v} - {\dot H_n}\end{align}

Equations (1) and (2), which describe mass conservation in the reservoir and plenum respectively, assume that all gas flows from the reservoir into the plenum and then through the nozzle, which is typical for most facilities. If some gas is lost or injected from a separate source (which is not the case for the present facility), the equations must be adjusted accordingly. Equations (3) and (4) describe energy conservation and assume the process is adiabatic, meaning there is no heat transfer between the environment and the gas at either the reservoir or the plenum.

To solve this system of equations, we first need to express the unknowns in terms of thermodynamic variables. The mass in both tanks is represented by the perfect gas law, giving us:

(5) \begin{align}M = V\rho = \frac{{VP}}{{RT}},\end{align}

where $M$ is mass, $V$ is volume, $P$ is pressure, $R$ is the specific gas constant, $\rho $ is density and $T$ is temperature. Similarly, the internal energy in both volumes can be expressed as:

(6) \begin{align}E = Me = \frac{{VP{C_v}}}{R},\end{align}

where $e$ is the specific energy, and ${C_v}$ is the specific heat at constant volume.

The mass flows in Equations (1) and (2) also need to be expressed in terms of thermodynamic variables. In this work, we assume that the nozzle is always choked. While this assumption does not hold at the start and end of a test, it is reasonable given that these intervals are relatively short and that the nozzle can achieve choking at comparatively low pressures, a scenario commonly encountered in the literature. With this assumption:

(7) \begin{align}{\dot m_n} = {A_n}\frac{{{P_p}}}{{\sqrt {{T_p}} }}\sqrt {\frac{\gamma }{R}} {\left( {\frac{2}{{\gamma + 1}}} \right)^{\frac{{\gamma + 1}}{{2\left( {\gamma - 1} \right)}}}}\end{align}

A choked nozzle is a common approach in literature [Reference Hwang and Hsu2], however, different methods are employed for modeling the flow through the control valve. In our facility, where the reservoir pressure is significantly higher than the target plenum pressure, the control valve also operates in a choked condition during normal operations. We can then express the mass flow through the valve as a function of the total pressure and temperature at the valve inlet. The important distinction is that, while the nozzle throat area is fixed (at least in our facility), the control valve throat area varies continuously. For our control valve, we determined the throat area as a function of the percentage valve opening using tables provided by the manufacturer.

Substituting these expressions into Equations (1)–(4) we obtain a system of four differential equations with four unknowns; plenum total pressure and temperature, and reservoir total pressure and temperature. The equations can then be solved through a numerical solver.

However, the system has an additional unknown through the valve area (or equivalently, the valve opening percentage). In general, this can be resolved in one of two ways: either by introducing additional equations to compute the valve area as part of a control strategy, as is commonly done in the literature to regulate plenum pressure (e.g. via a PID controller), or by prescribing the valve area directly as an input to the model. Because the objective of the present work is to develop and validate a new model, rather than a new controller, the latter approach is adopted. Accordingly, the valve opening is specified as a function of time, using measurements recorded during each experimental run.

In deriving this system, we have made a series of assumptions. We assume that the flow is adiabatic throughout the system and that there are no pressure losses, except in the control valve. Additionally, we have neglected the dynamics of the piping that connects the reservoir to the plenum. This leads us to conclude that the mass exiting the reservoir is identical to the mass entering the valve. Therefore, in Equation (1), we used ${\dot m_v}$ rather than ${\dot m_r}$ .

It is important to note that in our facility, the control valve is positioned right before the plenum, with the majority of the piping running from the reservoir to the valve – a typical configuration. Lastly, we assume that the working gas, which is typically air, behaves as a perfect gas.

Figure 2 shows a comparison of this base model with experimental data from a supersonic tunnel test run. Figure 2(a) compares the pressures predicted by the model with experimental measurements for both reservoir and plenum. The model’s plenum pressure compares well with the experimental measurements, with an average error of 5% for this particular test. Since most literature works were interested in the plenum pressure, it is then clear why the base model from Fung et al. [Reference Fung, Settles and Ray1] and the similar successive ones were deemed sufficient. However, Fig. 2(b) shows a different trend for temperatures. The adiabatic assumption leads to temperatures significantly lower than reality with a 34% error towards the end of the test. Additionally, the base model fails the capture the initial temperature drop observed in experimental data, which is due to Joule–Thomson effect. The error in temperature is reflected in the reservoir pressure in Fig. 2(a). Since mass flow out of the reservoir scales with $1/\sqrt T $ (Equation (7)), the model over-predicts the mass flow exiting the reservoir and thus its pressure decreases faster than reality. An interesting feature in Fig. 2(b) is the spike in plenum temperature at the beginning of a test. The sharp temperature increase is seen in both simulation and experimental measurements and it is due to the initial compression of the air [Reference Leavy19]. Lower values in the experimental case can be attributed to the thermal mass of the plenum (including flow spreader and screens).

Figure 2.

Comparison between base model and experimental data. (a) Pressures, (b) Temperatures.

The results in Fig. 2 show that the experimental temperature is not correctly captured by the base model. Two major differences are visible in the actual measured data. First, the process is diabatic since temperature remains significantly higher in both the valve inlet and the plenum, and second a sharp drop in plenum temperature is visible at the beginning of a test. These observations motivate the inclusion of non-ideal effects, which are introduced and discussed in the following section.

3.2 Heat addition

The primary driver of the diabatic process is the thermal mass of the tanks and pipes that make up the system. As the working gas expands and cools, it removes heat from the (mainly) steel. This is true in all parts of the wind tunnel, the reservoir, the piping, the plenum, nozzle, etc. However it is more severe in some parts than others. In our facility, heat addition is assumed to occur upstream of the control valve, primarily in the long piping connecting the reservoir to the valve, while heat transfer from the reservoir tanks is implicitly captured by the energy balance. This assumption is motivated by the forced-convection nature of the heat transfer process: although the flow remains well within the subsonic regime, gas velocities are higher in the piping than in the reservoir or plenum due to the smaller diameter.

Three observations corroborate this hypothesis. First, Fig. 2(b) shows there is little difference (maximum 7%) between trends in the temperature at the control valve inlet and in that in the plenum. This indicates that the heat transfer in the plenum is small and can be considered adiabatic. Second, observations during multiple experiments showed that the pipes were significantly colder than the reservoir or plenum at the end of a test, pointing to a much more significant heat removal. Third, a predictive version of the model, presented later, relies exclusively on the thermal mass of the piping and does not use experimental temperature data; despite this, it shows good agreement with measured results, further supporting the assumption that the dominant heat transfer occurs in the piping.

To model the heat addition we first rely on experimental data. The temperature at the valve inlet ${T_{vi}}$ is recorded during each test. We can then define ${\rm{\Delta }}T = {T_{vi}} - {T_{va}}$ , where ${T_{va}}$ is the valve inlet temperature under adiabatic conditions. ${\rm{\Delta }}T$ then is the temperature difference between adiabatic and diabatic expansions. If we assume a quasi-steady process, we can then express the instantaneous thermal power added in the pipes as:

(8) \begin{align}\dot q\left( t \right) = {\dot m_v}{\rm{\Delta }}h\left( t \right)\end{align}

where the specific enthalpy, for a perfect gas, is calculated as ${\rm{\Delta }}h\left( t \right) = {C_p}{\rm{\Delta }}T\left( t \right)$ . The issue with this approach is that each test has a different function of $\dot q\left( t \right)$ , thus not allowing the model to predict future experiments. This is shown by the different curves in Fig. 3(a). This figure shows the temperature difference for 12 wind tunnel tests. These tests were performed at different plenum pressure settings (thus with different mass flow rates). The large scattering of temperature differences at time equal 0 is due to the fact that while most of the piping is indoor (and thus at an almost constant temperature throughout the year), the reservoir is outdoor and so its temperature changes seasonally and daily based on sun exposure.

Figure 3.

(a) Temperature difference vs time, and (b) total energy vs mass.

An alternative approach is to express the integral enthalpy (Equation (9)) as a function of the integral mass (Equation (10)) through the control valve. The integral mass is independent of the specific testing parameters and as such is more generally applicable. Similarly, we hypothesised that the energy extracted from the piping should be a function of the mass of gas processed. Figure 3(b) shows the integral enthalpy against the integral mass for the same 12 tests from Fig. 3(a). In this case the data collapse on a quadratic curve. Since this curve is not function of the specific testing conditions, we can now use it in our model to predict heat addition.

(9) \begin{align}{H_I}\left( t \right) = \mathop \smallint \nolimits_0^t {\dot m_v}{\rm{\Delta }}hdt\end{align}

(10) \begin{align}{M_{vI}}\left( t \right) = \mathop \smallint \nolimits_0^t {\dot m_v}dt\end{align}

To do this, we use the data in Fig. 3(b) to extract a best fit polynomial curve. For this model we selected to express the enthalpy as ${H_I} = BM_{vI}^2$ , where B is a coefficient. The y-intercept is zero by definition. The linear term was also imposed as zero (it will shortly be explained why). The coefficient B is then calculated from a best fit of the existing data. For the values in Fig. 3(b), the resulting curve has a ${R^2} = 0.9648$ , showing a good approximation.

(11) \begin{align}\frac{{d{H_I}}}{{d{M_{vI}}}} = {\rm{\Delta }}h\end{align}

(12) \begin{align}{\rm{\Delta }}T = \frac{{2B{M_{vI}}}}{{{C_p}}}\end{align}

In order to use this expression, we can take the derivative of the total enthalpy with respect to the mass, as shown in Equation (11). Doing the same with the polynomial expression, we can finally express the change in temperature as a linear function of the integral mass as shown in Equation (12). The last equation shows why the linear term in the polynomial approximation was set to zero, so that at 0 mass, there is no heat addition. Note that, while Equation (12) can now be used to predict the effects of heat addition, the coefficient $B$ is facility-specific and requires experimental data to be calculated.

Finally, by adding Equation (12) to the model, we can now calculate the correct valve inlet temperature (since the mass flowrate through the valve is already calculated, it is easy to obtain the integral mass). Figure 4 shows the comparison between the model with heat addition and the same experimental data from Fig. 2. The pressure plot (Fig. 4(a)) shows little difference from the base model, with only the reservoir pressure prediction improving. However, the temperature plot (Fig. 4(b)) now shows a very good match between the simulated and the measured valve inlet temperatures. The plenum temperature prediction also improved, albeit it is now over-predicted. The reason (and the solution) for this is in the following section.

Figure 4.

Comparison between base model with heat addition and experimental data. (a) Pressures, (b) Temperatures.

The experimental cases used to determine the heat-addition coefficient $B$ span the full operational total-pressure range of the facility, as illustrated by the 12 tunnel runs presented in this study. The good agreement observed across these conditions indicates that the proposed heat-addition model remains applicable over a wide range of operating pressures. It is noted, however, that the coefficient $B$ is facility-specific, as it represents the combined heat-transfer characteristics of a given tunnel configuration, including geometry, materials and convective environment; nevertheless, the methodology used to model heat addition and to determine $B$ from experimental data is general and may be applied to other blow-down facilities, and in the present case a limited number of experimental runs (12) was sufficient to determine it. Finally, while the present facility is characterised by a large thermal mass upstream of the control valve, the same modeling approach may be extended to facilities with significant piping downstream of the valve by introducing an additional energy balance describing the temperature evolution between the valve and the plenum, along with a second heat-addition coefficient (e.g. ${B_2}$ ), without altering the underlying modeling framework.

One issue remains with the presented solution for heat addition, it requires experimental measurements to obtain the coefficient B, as such it can’t be used in the design phase of a new facility. However, a physical explanation of the process can be provided. The model assumes that the heat is removed from the pipes connecting the reservoir to the control valve. Since the length of these pipes is known (and could be estimated during the design phase), as well as their mass per unit length, we can calculate their thermal mass. We assume that the pipes do not exchange heat with the surrounding, such that heat is exchanged only between pipes and internal gas. This is a reasonable assumption since the forced convection inside the pipe is considerably larger than the natural convection between the pipes and environment. We can also assume that, since the convective heat transfer is large and the wall thickness of the pipes is not large, the pipes are at the same temperature of the gas. We can then express the valve inlet temperature as:

(13) \begin{align}\frac{{d{T_{vi}}}}{{dt}} = \frac{{{{\dot m}_v}{C_p}\left( {{T_r} - {T_{vi}}} \right)}}{{{M_{pipes}}{C_{pipes}}}}\end{align}

where ${M_{pipes}}$ and ${C_{pipes}}$ denote the mass and specific heat of the pipes, respectively. Results from implementing Equation (13) into the base model are shown in Fig. 5. The pressure plots (Fig. 5(a)) show minimal variation relative to experimental pressure traces. The valve inlet temperature prediction in Fig. 5(b) initially aligns with experimental data but gradually deviates, reaching a difference of approximately 10 K by the end of the run. These results were obtained without relying on any experimental data. Consequently, both approaches have distinct merits: Equation (13) is suited for the design phase, whereas Equation (12) yields more accurate predictions once experimental data are available. For the remainder of the manuscript Equation (12) will be used. It is important to note that, while in our facility the thermal mass is provided by the long piping system, other facilities may incorporate thermal mass directly within the reservoir or in a dedicated thermal mass tank [8]. In such cases, Equation (13) remains applicable, with ${M_{pipes}}$ and ${C_{pipes}}$ replaced by the corresponding properties of the respective thermal masses.

Figure 5.

Comparison between base model with predicted heat addition and experimental data. (a) Pressures, (b) Temperatures.

3.3 Joule–Thomson effect

The Joule–Thomson effect is a real gas effect where total temperature is not conserved in an isenthalpic process. The consequence is that when a gas is expanded, we can observe a change in temperature. Although air is typically considered a perfect gas, large pressure changes, such as those occurring when throttling through a control valve, can make the Joule–Thomson effect significant. For temperatures and pressure ratios typically encountered in supersonic wind tunnels, the Joule–Thomson effect decreases the total temperature after the control valve. Joule–Thomson is a well-documented phenomenon in wind tunnel design literature [Reference Pope and Goin8–Reference Lukasiewicz11].

The Joule–Thomson effect is more significant for larger pressure ratios. Thus, for blow-down wind tunnels, it is most severe at the beginning of a test, when the reservoir pressure is at its maximum, and decreases as the reservoir pressure drops. This trend is opposite to that of gas cooling due to expansion in the reservoir, which becomes more pronounced as the reservoir pressure approaches its minimum. Consequently, while the Joule–Thomson effect causes an initial temperature drop across the valve, it mitigates plenum temperature drift toward the end of a test and helps maintain a more constant temperature overall.

Figure 2(b) shows the presence of the Joule–Thomson effect in our facility. We can observe a drop in temperature between the valve inlet and the plenum temperatures. This drop is highest at the beginning of the test and slowly decreases as the reservoir pressure drops. On the other hand, the present numerical model predicts that valve inlet and plenum temperature are approximately the same during the whole test, which will now be addressed.

To model the Joule–Thomson effect and real gas effects in general, we replace the calculation of energy and enthalpy in the model. Up to now, under the assumption of a calorically perfect gas, energy and enthalpy were calculated using constant specific heats. Instead, we now use real gas tables from the work of Lemmon et al. [Reference Lemmon, Jacobsen, Penoncello and Friend20], for which enthalpy (or energy) is expressed as $h = f\left( {P,T} \right)$ . Importantly, since Equations (1)–(4) are expressed in terms of energy and enthalpy, the original model remains valid; only the calculation of energy and enthalpy has changed.

3.3.1 Heat addition for real gases

In the previous section, we discussed heat addition effects for a calorically perfect gas, although real gas effects also influence this process. Equations (9) and (11) are written in terms of enthalpy, so they remain valid. Similarly to the previous section, we then modify the calculation of enthalpy. Instead of using ${\rm{\Delta }}h\left( t \right) = {C_p}{\rm{\Delta }}T\left( t \right)$ for a calorically perfect gas, we now calculate enthalpy from real gas tables [Reference Lemmon, Jacobsen, Penoncello and Friend20]. It is therefore more convenient to re-express Equation (12) in terms of enthalpy rather than temperature:

(14) \begin{align}{\rm{\Delta }}h = 2B{M_{vI}}\end{align}

Equation (14) provides the change in enthalpy due to heat addition. However, using the same methods previously discussed, the coefficient B must now be determined from experimental data using real gas tables to calculate enthalpy. The results obtained with real gas considerations collapse again onto a parabolic curve like in Fig. 3(b). However, they now yield an improved correlation coefficient of ${R^2} = 0.9903$ . Thus, incorporating real gas effects improves the accuracy of heat addition predictions.

Figure 6 shows the results for the model with real gas effects. Once again, pressure plots in Fig. 6(a) are not affected significantly. However, Fig. 6(b) now shows a much better agreement on the plenum temperature. The drop observed from the experiments is now replicated by the model. It is important to note that real gas effects do not require any previous experimental data to be modeled and thus can be implemented in the design phase.

Figure 6.

Comparison between model with heat addition and Joule–Thomson effect and experimental data. (a) Pressures, (b) Temperatures.

3.4 Pressure losses

The last effect we are going to add to the base model is not about temperature, but rather pressure. Figure 7(a) shows the measured total pressure difference between the reservoir and the valve inlet versus time. This figure shows that the pressure losses in the piping system can be significant. For high plenum pressure testing, which means larger mass flow rates (and thus velocities) through the pipes, pressure losses can exceed 1 MPa, which is on the same order of magnitude of the target plenum pressure. Thus, pressure losses should be included in the wind tunnel model.

Figure 7.

(a) Measured pressure losses vs time. Each colour represents a different experiment. (b) Pressure losses vs mass flow times velocity.

In the current model we have assumed that the pressure at the control valve inlet is equal to that in the reservoir ( ${P_{vi}} = {P_r}$ ), we are now going to derive a model that calculates the pressure at the valve based on losses in the pipes. The flow in the piping system can be considered quasi-steady, and due to the low speed involved, incompressible. We can then use Darcy-Weisbach equation to model the pressure losses in the system.

(15) \begin{align}{\rm{\Delta }}P = f\frac{L}{D}2\rho {v^2}\end{align}

By replacing density using continuity, and expressing in term of valve inlet pressure, we obtain:

(16) \begin{align}{\rm{\Delta }}P = {P_r} - {P_{vi}} = 2f\frac{L}{{DA}}\dot mv\end{align}

In Equation (16) $L$ , $D$ and $A$ are the length, diameter and cross-sectional area of the pipes. The friction factor $f$ is normally a function of the Reynolds number, however, since this is in the ${10^7}$ range, the friction factor is a constant [Reference Moody21]. Minor losses associated with fittings and geometric changes have the same functional dependence on velocity and are therefore implicitly captured within the pressure-loss formulation adopted here. We then expect a linear relation between the pressure losses and the product of mass flow and velocity. Figure 7(b) shows the measured pressure loss (from the same 12 tests of Fig. 7(a)) vs mass flow times velocity ( $\dot mu$ ) as well as a linear best fit. Although experimental measurements depart slightly from linearity at higher mass flowrates, the data fit well with the linear approximation with an ${R^2} = 0.9804$ . This shows that the pressure difference is due to frictional losses in the pipes. Equation (16) can then be added to the model to predict pressure losses in the pipes.

We can then calculate the pressure at the valve inlet from the reservoir pressure minus the pressure loss. This adds one algebraic equation to the model. Figure 8 shows the results from this model, which now includes head addition, Joule–Thomson and pressure losses. The temperature plots in Fig. 8(b) do not change significantly. However, Fig. 8(a) shows an improved agreement between all pressures. Note that the valve inlet pressure is now also plotted, it was omitted before since the model assumed it to be identical to the reservoir pressure.

Figure 8.

Comparison between model with heat addition, Joule–Thomson effect, and pressure losses and experimental data. (a) Pressures, (b) Temperatures.