On the Modelling of Asymptotic Wavefronts in Long Ducts with Chambers

Abstract

A novel method of determining the possible shapes of pressure wavefronts in ducts after they have travelled sufficient distances to evolve to asymptotic states is introduced. Although it is possible in principle to achieve the same outcome by simulating complete flow histories from the time of the creation of the wavefronts, this can be impracticable. It is especially unsuitable to use such methods when extremely small grid lengths are needed to represent the final outcome adequately. The new method does not simulate the propagation phase at all. Instead, it explores what final end states are possible, but gives no information about the initiating disturbance or the wavefront evolution towards the assessed asymptotic state. Accordingly, the two methods do not overlap, but instead are complementary to each other. A typical case in which the new capability has high potential is described and used to illustrate the purpose and use of the methodology. However, the primary focus is on the presentation and assessment of the method, not on any particular phenomenon. It is shown that the required computational resources are far smaller than those needed for conventional unsteady flow simulations of propagating wavefronts. The potential numerical limitations of the method are highlighted and, with one exception, are shown to be either of no consequence or easily reduced to acceptable levels. Special attention is paid to the one exception because it cannot be proven to be unimportant and, indeed, it would be unsafe to use it in general analyses of wave propagation. However, strong evidence is presented of its acceptability for the study of asymptotic wavefronts.

1. Introduction

In many, perhaps most, studies of unsteady fluid flows in a pipe or duct, it is sufficient to treat the whole flow as moving almost as a solid body. This is also true for individual pipes in networks of connected pipes. However, when studying the consequences of disturbances induced suddenly—e.g., rapid valve closures—account needs to be taken of fluid compressibility, which enables rapid accelerations to occur without inducing excessive pressure changes (Fox 1977 [1], Wylie and Streeter 1993 [2], Chaudhry 2014 [3]). In such cases—commonly known as ‘water hammer’, even when the fluid is not water—the initiating event first influences only the fluid closest to it. This then influences the adjacent fluid before that influences its adjacent fluid, and so on. That is, the region of flow that has been influenced by the event continually increases in length. The interface between this region and the region of flow that has not yet been affected typically propagates along a pipe at approximately the speed of sound. In this paper, this interface is at the leading part of a well-defined wavefront causing a finite increase in pressure. The trailing part of such wavefronts is usually less clearly defined. Loosely, however, it is where the region of rapid change ceases and a subsequent region of much slower change begins.

Although pressure wavefronts typically propagate at approximately the speed of sound, they usually change in shape as they do so. This is because the local speed of propagation depends upon the local flow state, which necessarily varies between the leading tip and the trailing edge. There are two primary reasons why this dependency exists. Attention is commonly drawn first to the fact that the speed of sound is itself a variable. It depends not only on the intrinsic properties of the fluid, but also on the local pressure and/or temperature. This dependence is more pronounced in gases than in liquids, but it exists for both. In practical flows, the speed of sound increases with increasing pressure; so, if the initiating event causes increasing pressure, then, in the absence of any other influence, the higher-pressure regions of the wavefront will tend to travel faster than its toe. Given sufficient time, they will catch up to the toe, and the wavefront will then be a single, discontinuous step, namely a shock. In contrast, if the initiating event causes decreasing pressure, the wavefront will lengthen indefinitely.

A second—and usually more important—reason for the differing speeds of propagation is that waves move at the speed of sound relative to the fluid, not relative to a fixed reference such as a pipe. Now, the initiating event necessarily causes changes in velocity, and these must be allowed for when estimating speeds of propagation. Again, therefore, the consequence will be a tendency for the wavefront to either lengthen or shorten, potentially into a shock.

Wavefront shortening and lengthening is rarely of importance in the case of liquid flows, except in special cases such as cavitating flows. However, it can be of critical importance in gas flows. Figure 1 shows pressure histories measured at successive locations along a railway tunnel in Germany. The graphs are offset in time to enable the shapes of the individual histories to be shown clearly. There are strong secondary effects before and after the main part of the wavefront, but the shortening effect is demonstrated clearly.

The challenge for engineers is to ensure that wavefronts do not become so steep that subsequent reflections at tunnel portals will cause unacceptably large audible phenomena in the surrounding environment. To do this, the first task is to assess whether the secondary effects will constrain the steepening process sufficiently during propagation. This depends strongly on the type of track. Historically, all track was ballasted, but slab track construction is increasing chosen because of its less demanding maintenance requirements. Reports of analyses and experimental measurements of wavefront propagation in these different types of tunnel include Mashimo et al. (1997) [4], Fukuda et al. (2006) [5], Nakamura et al. (2011) [6], Wang et al. (2018) [7], Fukuda et al. (2020) [8] and Iyer et al. (2021) [9]. When such investigations do not confirm that the resulting steepness will be acceptable, the most common design response is a modification of the tunnel entrance region to reduce the initial steepness of train-entry wavefronts. Alternatively, a reduction in the train-entry speed can have the same effect. Another possibility, so far explored only theoretically, involves utilising unused space beneath tunnel walkways to create air chambers that mimic a key effect of ballast tracks (Liu et al. 2021 [10]). The methodology developed below has been devised to address a difficulty encountered when assessing the final shapes of such wavefronts after propagating long distances in such a tunnel. The method is potentially of use for studying other asymptotic wavefronts, but the only specific application discussed below is the particular case of air chambers.

1.1. Practical Flows

In real flows of finite amplitude, the characteristics of an evolving wavefront depend on the environment in which it is propagating, as well as upon the initiating event and the inertial behaviour discussed so far. An almost universal complication arises as a consequence of fluid viscosity and of the so-called ‘no-slip’ condition that causes fluid velocities on a solid surface to remain the same as the velocity of the surface itself. This causes frictional resistance to flow that influences the pressure–flowrate dependency and hence changes the evolution of wavefronts. This phenomenon is quite complex because its consequences can depend upon rates of change of flow, as well as upon absolute values of these changes. Broadly speaking, however, long regions of slowly varying flow are influenced primarily by low-frequency behaviour of friction, whereas short wavefronts are influenced primarily by high-frequency behaviour.

Other phenomena that influence some flows include (i) a non-uniform pipe area, (ii) heat transfers with pipe walls and (iii) the leakage of mass through the pipe walls (whether inflow or outflow). The author expects the generic methodology described herein to be used for a range of purposes, and this is reflected in the presentation given below. However, in the particular example presented in Section 5 to illustrate the use of the methodology, all influences of friction, non-uniform area and heat transfer are neglected. Attention is focussed on compressive wavefronts and on a particular form of leakage that, in principle, has the potential to exactly counterbalance the inertial-shortening process illustrated in Figure 1 and thereby prevent the adverse consequences associated with very short wavefronts.

The chosen geometrical configuration is a duct that has wall holes through which air exchange with adjacent air chambers is possible (see Figure 2). The chambers are isolated from one another; so, although they contribute to the effective bulk compressibility of air in the duct, they do not provide a path for axial flows along the duct. This configuration is being investigated to assess its potential for countering wavefront steepening in railway tunnels, thereby reducing the likelihood of sudden pressure disturbances inside a tunnel causing unacceptable disturbances outside it. In a nutshell, the challenge for tunnel designers is to prevent wavefronts becoming steeper than a prescribed upper limiting value.

1.2. Structure of the Paper

The remainder of this paper begins with a qualitative discussion of some characteristics of wavefronts, with special attention being paid to a condition that might be approached asymptotically in long ducts. Then, two sections are devoted to describing the numerical method used herein to investigate this behaviour. Unusually, instead of prescribing a cause and predicting its effect, the method assesses (the characteristics of) a possible effect without making any assumption about its possible cause. This enables the asymptotic behaviour to be studied in fine detail with high accuracy, without the complications that usually arise in studies of wavefront propagation in long ducts. The method is then used to demonstrate a particular characteristic of asymptotic wavefronts, and the paper closes with a short section giving a brief summary and highlighting the key conclusions.

2. Wavefront Shapes

In most water hammer-like applications, wavefronts are quite long, and this is recognised implicitly in the terminology used to define the phenomenon. For example, in the classical reservoir pipe valve case, any valve closure that causes a wavefront length to be shorter than two pipe lengths is termed ‘rapid’. Then, neglecting the influences of other causes of change, the maximum pressure reached at the valve will be almost independent of the detailed sequence of the closure. In contrast, if the valve closure time exceeds that required for the tip of the initial disturbance to travel all the way along the pipe and back to the valve before the closure is complete, the maximum pressure will be smaller. The amount of the reduction will depend on the closure sequence and hence upon the shape of the wavefront—e.g., ‘steep near the tip and less steep afterwards’ or ‘shallow near the tip and steeper afterwards’, as illustrated in Figure 3.

An important feature of this illustration is that the shape of the wavefront exerts its influence only when the wavefront reaches a boundary. Another example of this is described by Tijsseling (2021) [11] in relation to the reflection of a wavefront at an open end of a pipe connected to a reservoir. During such reflections, most of the internal energy of the fluid is ‘converted’ to kinetic energy, and the reflected wavefront travels upstream to induce this change. However, a small amount of energy is released into the reservoir during the reflection process, causing waves to radiate outwards, and, in this context, the shape of the wavefront is important. The amplitudes of the radiated waves depend upon the variation in steepness along the wavefront. That is, they depend upon the rates of change of pressure and velocity. Counter-intuitively, the amplitude of the incident wavefront is not important in its own right. Instead, what matters is the rate at which the change occurs.

Waves radiated into water reservoirs during wavefront reflections are of little practical significance except, perhaps, in highly specialist activities such as sonic surveillance. However, the same behaviour can be hugely important in air. An obvious instance is the design of musical instruments that exploit the phenomenon to delightful effect. Unfortunately, this contrasts sharply with the railway tunnel example cited above. That process is essentially the same, but the consequence can be far from delightful. The most common initial cause of unacceptable wavefronts occurs when the nose of a train enters a tunnel at high speed. Initially, the length of the wavefront might be in the order of, say, ten tunnel diameters, but it can shorten greatly as it propagates. If the tunnel is sufficiently long, the wavefront can be very steep when it reaches the exit portal of the tunnel. In extreme cases, a loud bang can be emitted from the tunnel. It is neither musical nor pleasant, so it has to be prevented. A YouTube link [12] to an example of this phenomenon is given in the references at the end of the paper.

Various counter-measures exist to reduce the amplitudes of the emitted disturbances. Most are designed to reduce the abruptness of the process of nose-entry to the tunnel. These include, for example, reducing the train speed (unpopular with operators), tapering the entrance region, allowing some air to escape to the atmosphere through holes in the walls of tunnel extensions (Winslow and Howe 2005 [13], Zhang et al., 2018) [14] or extending train noses (Kikuchi et al. 2011 [15], Miyachi et al. 2022 [16]). Perhaps unexpectedly, these methods are more reliable and more effective than similar measures at the tunnel exit, although the latter have received attention—even including active measures, as well as passive ones analogous to those at tunnel entrances (Nishimura et al., 1994 [17], Aoki et al., 1999, [18], Wang et al., 2015 [19]). A third approach is to implement measures along the tunnel to prevent the unacceptable steepness of wavefronts arriving at the exit portal (Hieke et al., 2009 [20], Liu et al., 2021) [10]. The last of these investigates the influence of air chambers by studying pressure wave propagation over large distances. However, the time steps used in its numerical calculations are considered to be too coarse for high confidence to be placed in the detailed shapes of the leading parts of the resulting wavefronts, even though useful conclusions could be drawn.

Asymptotic Wavefronts

As indicated above, inertial forces tend to cause compressive wavefronts to shorten, and hence steepen, as they propagate, whereas other phenomena tend to have the opposite effect. For some time after its initiation, the shape of a wavefront will depend primarily on the initiating event itself—e.g., on rates of change of flow during a valve closure. Gradually, however, the shape will depend increasingly on the environment in which the wavefront is propagating and less and less on its initial shape. In sufficiently long ducts, an approximate balance will gradually develop between the steepening and flattening influences, and, thereafter, the wavefront will propagate in an almost unchanged form. Thus, for example, all three wavefronts in Figure 3 would tend to evolve towards the same eventual shape. An important corollary of this fact is that it is not possible to infer the initial shapes from a knowledge of the asymptotic states. However, the characteristics of the asymptotic conditions are of interest in their own right, independently of how they were created and irrespective of particular practical applications.

In principle, it is possible to study asymptotic wavefronts by simulating the whole process of propagation until an asymptotic condition is reached. Indeed, this is exactly how the present author has studied them over a period spanning several decades. However, that method has both practical and technical disadvantages. The most obvious practical one arises from the need to use small grid lengths to represent conditions along steep wavefronts in appropriate detail. Unfortunately, this same grid size must be used along the whole length of the duct. Furthermore, the corresponding integration time steps must also be small. Another big disadvantage is that current understanding of pressure wave damping during propagation over long distances inside ducts is weak. Furthermore, there is a risk of the accumulation of numerical errors at the wavefronts themselves. At best, these difficulties reduce the level of confidence that can realistically be justified in predicting the detailed shapes of wavefronts when they reach the exit.

Vardy and Brown (2000) [21] reduced some of the disadvantages by simulating only within a limited band of (x,t) space moving with the wavefront. This greatly reduced CPU requirements, but it also broke the link between the wavefront and the upstream flows, thereby partially negating the apparent benefit. More important, in all of the solutions, the leading region of the simulated wavefronts developed a ramp-like shape that the authors were unable to explain with confidence. They hypothesised that it was probably a numerical consequence of interpolation, but were unable to rule out the possibility of it being a true physical phenomenon.

In the following sections, problems associated with simulating wavefront propagation per se are avoided completely. Instead of asking what final state will result from prescribed initial conditions, the analysis focusses exclusively on possible end results. The aim is to see what final shapes are possible for wavefronts that have travelled sufficiently far along ducts/tunnels of a prescribed geometry. As indicated above, the answer to this question does not depend upon the process that initiated the wavefront—although this will influence the length of duct required for the asymptotic state to be approached closely.

This approach might seem like ‘putting the cart before the horse’, but it is potentially valuable because it enables highly accurate predictions of the range of possible asymptotic states. At the very least, it will enable the outcomes of more conventional analyses to be assessed with greater confidence than hitherto. Importantly, the method of achieving this benefit requires only trivial CPU resources, so large numbers of cases could be studied rapidly.

For completeness, it is emphasised that the analysis cannot be extended backwards in time to identify what initial conditions could have led to the eventual outcome. The primary reason for this is that the method enforces asymptotic behaviour everywhere, regardless of how far back in time (or how far back in distance from the wavefront tip), the analysis is extended. Conversely, of course, analyses that begin with the initiating event cannot stipulate the characteristics of the resulting wavefronts a priori.

3. Numerical Grid

This numerical simulation is a particular form of the Method of Characteristics (MoC), which is widely used for the simulation of waves in pipelines. It is especially popular in studies of low Mach number flows in which strong changes occur over short timescales. It has a well-known inability to satisfy conservation laws exactly, but, for the purposes of studying details of highly localised wavefronts, that is far less is less important than the spurious oscillations that can arise when using the control volume or finite element methods. Furthermore, the most important aspect of the following development is not the MoC equations themselves, but the method of specifying the numerical grid. Accordingly, this is discussed first, and the presentation of the relevant equations is deferred to the following Section.

Figure 4a depicts the profile of an isolated, asymptotic wavefront at three instants as it travels at a steady velocity U_w (NB: all points on the wavefront are travelling at the same speed; otherwise, asymptotic conditions would not exist). At the time indicated by the line A-A, the wavefront has almost passed location C, but its tip has not yet reached location D. At the second instant, the whole wavefront is between these locations, and nominally steady conditions exist at both of them. Ambient conditions still prevail at D, whereas, at C, the pressure and velocity have changed in response to the wavefront. At the third instant, depicted by the line B-B, the wavefront has begun to disturb conditions at D. Importantly, because the wavefront is prescribed as being asymptotic, all three curves in the figure are identical. That is, the spatial profiles at successive instants differ only by a uniform distance. Likewise, the temporal histories at different locations differ only by a uniform time interval. Moreover, these statements are true regardless of whether the intervals are small or large. The chosen numerical grid takes advantage of this unique characteristic behaviour of the asymptotic state.

Figure 4b shows a fixed rectangular grid. It is designed to ensure that the ratio ∆x/∆t is exactly equal to the speed of the wavefront, U_w. As a consequence, valid solutions at every grid point will necessarily be identical except for offsets in time and space. The methodology described in the following section enforces this behaviour. It is worth emphasising that this grid structure is possible even though the wavefronts under investigation are of finite amplitude. However, it would not be possible if the wavefront were not in an asymptotic state or if other disturbances were travelling in the opposite direction. Either of these complications would lead to different physical behaviours at different locations, just as they do in more general simulations of pressure surges in pipelines. The only other instance in which a fixed rectangular grid can exactly match the physical propagation of waves is with vanishingly small-amplitude disturbances. In that case, there is no need to stipulate either asymptotic behaviour or an absence of waves travelling in the opposite direction. However, when these restrictions are removed, it is no longer the case that identical solutions exist at every location in space and time.

The use of a numerical grid structure satisfying ∆x/∆t = U_w has two important benefits and one important disbenefit. These are discussed fully below (Section 4.2), but they can be summarised as follows:

Benefit-1: It enables the numerical methodology to ensure that all points along the whole of the finite-amplitude wavefront propagate at an identical, constant, physically-correct speed. A loose analogy can be drawn with the widely accepted practice of forcing all points on general wavefronts in low Mach number, transient flows to propagate at exactly the same speed (usually chosen as the ambient speed of sound). However, that is an approximation to the true physics. It is made for simplicity and is justified only when the resulting errors are acceptable.

Benefit-2: It enables the speed of propagation of the wavefront to be chosen a priori by the software user. Indeed, it enforces that action. Use is made of this below (Section 5.3) when illustrating a practical application of the methodology. It is worth noting here that the chosen speed must be below the speed of sound. Propagation at a greater speed would imply that the wavefront is an abrupt discontinuity, namely a normal shock, whereas the analysis used herein is valid only in the absence of discontinuities.

Disbenefit: The downside of choosing the numerical grid to satisfy ∆x/∆t = U_w is that this makes the use of extrapolation inevitable—because the gradients of the MoC characteristic lines are greater than the speed of sound. This is a serious matter; numerical methods that involve extrapolation are usually inherently unstable. Accordingly, repeated attention is drawn to this concern below, including an assessment in which the potential errors are magnified tenfold.

In the particular software used to obtain the outcomes presented below, the numerical solution is obtained at successive times at a single spatial location, with use being made of the asymptotic property to infer values at adjacent locations and times. Hence, the numerical grid needs to include only this location and grid points very close to it. Nevertheless, the analysis has to include every time step in which the wavefront has an influence at the chosen location, and the number of such increments could be quite large when a fine grid pattern is used. However, in the author’s software, the grid exists only at time steps close to the current solution instant. This is made possible by nesting the whole solution at the end of each time step. As a result, the complete solution is obtained using only 4 × 4 arrays. It is emphasised, however, that this choice has been made solely for programming convenience. It has no influence whatsoever on the resulting solution, so little purpose would be served by expanding on the relevant programming details herein.

4. Numerical Method

The cross-sectional-averaged continuity and momentum equations for 1D unsteady flow in the presence of mass addition/subtraction (through the duct walls) and skin friction (on the wall surfaces) can be combined and expressed in characteristic form as

$$1\mathrm{D}\mathrm{MoC}:{\displaystyle \frac{dp}{dt}}\pm \rho c{\displaystyle \frac{dU}{dt}}=-c^2{\displaystyle \frac{{\dot{m}}_y}{A}}\mp c\left\{-{\displaystyle \frac{{\dot{m}}_y}{A}}\left(U_{yx}-U\right)-{\displaystyle \frac{{\tau }_w\mathcal{l}}{A}}\right\}$$

(1)

in which all symbols are defined in the nomenclature.

In general, the wall shear stress τ_w can have both steady-flow and unsteady-flow components, but, for simplicity, it is treated as zero herein because the principal purpose is to introduce a new methodology, not to study particular practical applications comprehensively. However, since τ_w = 0, the lateral mass flowrate ${\dot{m}}_y$ cannot also be zero. If it were, the right-hand side of Equation (1) would be zero and nothing would exist to counter the inevitable long-term consequence of inertial steepening, namely the evolution of the wavefront into a normal shock. At the other extreme, if the holes in the duct wall led directly to the atmosphere, the cumulative mass flowrate would be unconstrained and, in effect, the wavefront would evolve in a dispersive manner. Indeed, assuming ambient conditions prevail indefinitely on the other side of the wall, the amplitude of the wavefront would decrease indefinitely with an increasing distance from its source. Instead, therefore, the wall holes lead to chambers of finite volume. Flows into these cause their internal pressure to increase, imposing increased resistance to further flows into them. After the wavefront has passed, the pressures in the tunnel and chamber are equal.

To preserve generality, at least provisionally, the mass flowrate per unit length of wall is represented generically by

$$\mathrm{Mass}\mathrm{flowrate}:{\dot{m}}_y={\dot{m}}_y\left\{a_y,V_{CH},{\rho }_A,{\rho }_{CH},p_A,p_{CH},K_y\right\}$$

(2)

where a_y is the area of holes per unit length of the duct, V_CH is the chamber volume per unit length and K_y is a resistance coefficient for flows through the wall. Details of the mechanism and its analysis are given by Liu et al. (2021) [10] and, using a different representation of the flow, by Vardy et al. (2022) [22]. In both cases, the importance of the resistance characteristics of the wall openings is highlighted. NB: If there were no resistance, the chambers would behave like Helmholtz resonators.

4.1. Illustration

Figure 5 shows indicative pressure histories at an arbitrary location during the passage of an asymptotic wavefront along a duct with air chambers. This example is discussed in greater detail in Section 5, which is focussed on the application of the methodology in the context of railway tunnels, but the qualitative behaviour is of general relevance in its own right. The example is included here for reference purposes when discussing the numerical solution process. It is seen that the pressure in the chambers lags behind that in the open duct. Indeed, a lag is inevitable whenever any wavefront propagates along such a duct, regardless of whether asymptotic conditions have developed. The magnitude of the lag depends upon the time required for air to flow into the chambers in response to the evolving pressure. The resulting wavefront can be asymptotic only if the total rate of mass flow into the chambers during the wavefront passage is exactly equal to the difference between the longitudinal mass flowrates upstream and downstream of the wavefront. If the lateral rate is insufficient, the wavefront will shorten. If it is too great, the wavefront will lengthen. The stated balance needs to exist relative to axes moving with the wavefront itself, not relative to axes that are stationary relative to the tunnel.

4.2. Details of the Numerical Integration Process

The focus now moves to the numerical integration process. This is the heart of the solution methodology, so it is described in detail. The methodology has many features in common with conventional MoC solutions using a fixed-grid approach. In principle, it is a time-stepping process in which the known solution up to some instant t_old is used, together with boundary conditions, to deduce the corresponding solution at a later instant t_new. In wider use, either implicit or explicit schemes may be used. With the former, the solution is obtained simultaneously at every grid point at t_new. With explicit schemes, it is obtained independently at each grid point. The present scheme is nominally explicit. Furthermore, since the whole solution needs to be obtained at only one spatial location, only one solution point exists at any instant t_new. In Figure 6, the point ‘A’ is the one at which a solution is to be obtained. The remaining points are as follows:

• M,N: Grid points at the most recent solution instant. The values of the flow parameters at these location are known from solutions obtained in the previous two time steps;
• B,C: Identical locations along the wavefront as A, but at different times. The flow parameters at these points are not known, but it is known that they are identical to the values at A (because, as stipulated in Section 3, ∆x/∆t = U_w);
• L,R: Locations on characteristic wave paths projected backwards in time from A;
• L2: The same spatial location as L, but at the previous time step;
• L′,L2′: The same locations along the wavefront as L and L2, but one time step later. The flow parameters at L′ & L2′ are identical to those at L and L2.

The shaded zone in the figure includes all points in (x,t) space in the current time interval on which the physical solution at A depends. Notice, in particular, locations outside the range x_L < x < x_R at the instant t = t_M are too remote for waves originating at them to reach x_A in one time step. This is a crucially important complication to which reference is made frequently below, even though its practical influence on the presented solutions is small.

Solution Sequence

The iterative solution sequence leading to values of the flow parameters at point A can be summarised as follows:

• Make an initial estimate of the pressure at A. Designate it p_AG;
• Infer the corresponding value of c_A, assuming isentropic change between A and the tip of the wavefront;
• Estimate the location of R and interpolate between M and N to obtain implied values of p_R, U_R and ${\dot{m}}_{y,R}$;
• Solve Equation (1) along RA and Equation (2) along MA simultaneously. Hence, deduce values of U_AR, ${\dot{m}}_{y,A}$ and p_CH,A consistent with the estimated p_AG. NB: The velocity at A is given the suffix AR to distinguish it from a second estimate of U_A obtained independently below;
• It is now possible to make an improved estimate of the location of R implied by the estimated value of p_AG. Hence, repeats steps 2 and 3 until an acceptable convergence is achieved. Then proceed to the following steps;
• Using a suitable equation for the gradient of an MoC wave path, derive a relationship between x_L and the local gradient of (U + c)_L. In the author’s software, the gradient of the wave path is evaluated as ½[(U + c)_L + (U + c)_A}, with U_A equal to the current estimate for U_AR obtained in step 4 above;
• Obtain a second relationship between x_L and the local gradient of (U + c)_L by extrapolation from M and B, noting that all values at B are the same as at A;
• Solve these two equations to give estimates of x_L and (U + c)_L. Use this solution to improve the estimated gradient of the MoC wave path, and repeat the process until convergence is achieved. Then extrapolate for the value of ${\dot{m}}_{y,L}$.
  NB: The use of extrapolation in this step will justifiably ring warning bells in the minds of readers with appropriate experience. A detailed discussion of it is given later;
• It is necessary to know U_L and c_L independently, whereas the preceding step has yielded only their sum. Solve this together with a standard relationship for a simple wave-like disturbance between L and B. Then use an isentropic relationship to infer p_L from c_L;
• Interpolate between B and M to obtain implied values of p_L2′, u_L2′ and ${\dot{m}}_{y,L2’}$ and note that identical values are applicable at L2;
• Solve Equation (2) between L2 and L to obtain the implied value of p_CH,L;
• Solve Equation (1) along RA to obtain a second estimate of the velocity at A, namely U_AL;
• At this stage, except for the dual value of U_A, everything is consistent with the chosen pressure p_AG. Accordingly, any difference between the two values of U_A will exist only because of an error in the estimated value of p_A. Therefore, a unique solution would be found if steps 1–12 were repeated with the true value of P_A. In practice, this value is found by repeating the whole process with successively improved estimates of p_AG. Convergence is achieved fairly rapidly because it is known a priori that increases in p_AG will increase U_AR but will decrease U_AL, and vice versa. Hence the sign of (U_AR − U_AL) determines whether the estimated p_AG needs to be increased or decreased;
• After completing step 13, the final solution at A is known, and the whole sequence is repeated at successive grid times until the whole of the wavefront has passed the chosen location.

It should be noted that some of the above steps have been left somewhat open in that they could be undertaken in any of a number of ways. However, provided that suitable choices are made, accurate solutions will be achievable by simply reducing the grid size until further reductions have a negligible influence. Unlike many computational studies, computer resources are of almost no consequence for the present purposes. As an example, the CPU time for each of the examples given in Section 5 was in the order of one second on a simple laptop computer.

4.3. Numerical Limitations

Before illustrating the use of the methodology, five possible (or apparent) limitations are raised and methods of minimising or eliminating potential adverse consequences are described.

4.3.1. Asymptotic Wave Speed U_w

Perhaps the most obvious apparent limitation is that the whole process begins by stipulating that the grid structure must satisfy ∆x/∆t = U_w exactly. This appears to conflict with the inevitable fact that designers will never know a priori what wave speed will exist in any particular application. However, this is misleading, because the whole basis of the analysis is to choose an outcome and then to explore how it could arise, not vice versa. As a consequence, we are completely free to choose whatever value of the wave speed we wish and to use the analysis to determine the remaining properties of wavefronts that could cause it. Naturally, it makes sense to focus on ranges of values that are expected to be physically meaningful, but that is a practical matter, not a formal requirement.

4.3.2. Wavefront Amplitude

Likewise, the overall amplitude of the wavefront cannot be known a priori, but again, it is misleading to regard this as a limitation. Unlike the wavefront speed, we are not free to choose it a priori, but there is no need to do so. Indeed, it is properly regarded as an output of the methodology, not as an input into it. This is illustrated by the role of the parameter p_AG in the detailed procedure presented above for the solution in a typical timestep. From this perspective, the final pressure after the whole wavefront has passed is simply the end result of a sequence in which the pressures at all intermediate times have been determined.

4.3.3. Interpolation

The third potential limitation, namely the need for interpolation to determine the location of, and parameter values at, the point R at the base of the MoC characteristic line RA in Figure 6, is real, not merely apparent. First, since the parameter values need to be inferred from values at the adjacent grid points M and N, they cannot always be exactly correct. Second, even if exactly correct values could be deduced, the location of point R itself would be subject to error because it has to be inferred from data at points A and R alone, even though, physically, it depends on local values everywhere along the true path between R and A. This issue is present in any MoC analysis based on numerical discretisation. Its influence cannot be eliminated completely, but it can be reduced greatly by using sufficiently small values of ∆x and ∆t. It is usually assumed to be very small if a series of simulations with successively smaller integration steps converge to a common solution. This is the approach adopted in the illustrative solutions presented below.

4.3.4. Extrapolation

Although errors caused by interpolation can reasonably be expected to reduce with reducing the grid size, this is not necessarily also true for errors resulting from extrapolation. Therefore, special attention needs be paid to the use of extrapolation to locate point L and to infer parameter values at it. In more general MoC analyses, it is usual to use one of two ways to avoid this problem. One includes regarding c as a constant and disregarding convective terms in the analytical equations, thereby yielding wave paths satisfying dx/dt = c instead of dx/dt = U ± c. However, that would be a nonsense in the present context because it would prevent the numerical simulation of the inertial steepening that the chambers are used to counteract.. The other common way of avoiding extrapolation involves using greater time steps than those illustrated in Figure 6, even though this simultaneously increases the influence of interpolation at points such as R. However, this approach is not possible herein because the whole solution process is applicable only for grids that satisfy ∆x/∆t = U_w. Instead, therefore, a more direct method of assessing the influence of extrapolation is used below.

Specifically, solutions obtained using simple linear extrapolation are compared with alternative solutions in which the changes implied by linear extrapolation are deliberately magnified. This is a demanding test in its own right, but the risk of serious error is further reduced by using extrapolation only for two parameters that can be expected to vary only slowly, namely ${\dot{m}}_y$ and the sum (U + c). The assessment process is illustrated in Section 5.2 below, where it is seen that differences between solutions obtained with and without magnification by factors of ±10 are small in comparison with the physical effect under study. It is not claimed that the test demonstrates that extrapolation will be safe in other contexts, but it is a strong indicator that its influence for the present purposes is weak, indeed probably negligible. It is certainly tiny in comparison with the physically important differences illustrated in Section 5.3.

4.3.5. Triggering the Integration

The final limitation considered in this section is more academic than practical, but it nevertheless merits attention. It relates to the definition of the leading tip of the wavefront. So far, the discussion has been based on an implicit assumption that the whole of the wavefront can propagate at a uniform speed. However, this cannot be exactly true, especially for the examples considered herein because, in all cases, the assumed wave speed is smaller than the speed of sound. In reality, therefore, there must be a leading portion between the tip of the simulated wavefront and a more distant location that will have been reached by an infinitesimal disturbance travelling at the speed of sound. Herein, this issue is circumvented by defining the location of the leading tip of the asymptotic wavefront to be at a position where the pressure has already increased by a small finite value. One of the examples in Section 5.2 compares outcomes obtained using different values of this increase.

5. Numerical Examples

The use of the methodology is now illustrated. First, features of an asymptotic wavefront corresponding to particular geometrical parameters are described in detail. Then, the sensitivity of the numerical processes to key factors is assessed and is shown to be acceptably small. Having established this, the method is used to investigate the physical relationship between the asymptotic wavefront and the amplitude of the upstream pressure.

5.1. Base Case Asymptote

Figure 7 shows predicted histories of the pressure, velocity and lateral mass flowrate during the passage of an asymptotic wavefront along a tunnel with air chambers. The origin of the time axis is at the leading tip of the wavefront, so these graphs are valid at all axial locations. Also, since the whole wavefront is propagating at a uniform speed, U_w, it would be equally acceptable to interpret the graphs as spatial profiles (of a wavefront travelling from right to left). The only change required for this is a scaling of the horizontal axis by U_w.

The cross-sectional area of the tunnel is 100 m², and the volume of the uniformly distributed chambers is equal to 3% of the corresponding volume of the tunnel. The rate of flow through connections between the tunnel and the chambers is prescribed to vary as the square of the flowrate. The wavefront speed is stipulated to be 2 m/s slower than the ambient speed of sound, so U_w = 338.4 m/s. Therefore, since the duration of the period when the rate of change of pressure becomes negligible (Figure 7c) is approximately 0.036 s, the length of the wavefront itself may be regarded as approximately 12 m. This is only slightly greater than the tunnel diameter, and hence the use of a one-dimensional methodology is open to question. However, this will depend upon the geometry of the flowpaths between the tunnel and the chambers. Although shown as discrete holes in Figure 2 above, it is expected that continuous slots will be preferable for practical reasons. Also, it is important to bear in mind that the asymptotic condition is the cumulative consequence of large numbers of chambers. Each individual chamber has only a small influence.

In this particular example, the pressure rises along the wavefront to a maximum of approximately 1.5 kPa, and the same maximum necessarily applies in the chambers too. Inevitably, the increase in the chambers lags behind that in the tunnel—because flow into the chambers can exist only whilst the pressure inside them is smaller than that in the tunnel. By inspection, in this instance, the duration of the period required for the pressures to almost equalise is slightly more than 10% of the duration of the overall pressure change caused by the wavefront itself. The graphs of the axial velocity (inside the tunnel) and the lateral mass flowrate per unit length into the chambers are included primarily for completeness. However, it is noteworthy that the maximum rate of flow into the chambers occurs long after the instant corresponding to the maximum rate of change of pressure in the tunnel.

The most obvious practical significance of the chambers is that they facilitate the existence of the asymptotic behaviour. That is, they prevent the wavefront steepening into a shock. Even so, the maximum steepness of the wavefront in this example, namely about 76 kPa/s, is significantly greater than would be acceptable for a railway tunnel portal close to buildings. It follows that chambers with the particular properties chosen here would not be suitable in a practical application, so other combinations of chamber volume and connection properties would need to be assessed. However, since the simulations require such small CPU resources, large numbers of such simulations can be undertaken in a short time. This is a huge practical benefit, because the alternative of exploring options by means of simulations involving propagation along the whole tunnel requires CPU resources that are many thousands of times greater (sometimes millions of times). For the avoidance of doubt, it is not claimed that the need for any such simulations is completely eliminated. Some will be necessary for a range of other reasons. Nevertheless, the required number of them will be greatly reduced.

5.2. Numerical Accuracy

Figure 8 illustrates the sensitivity of the base case solution to various numerical factors. This is done by comparing rates of change of pressure, not the pressure itself. There are two reasons for this. The more important of these for the theoretical modelling purposes of this paper is that the rates of change are more sensitive than the absolute values, so errors can be seen more easily. However, it is also true that the physical application that triggered the research depends far more strongly on rates of change of pressure than on the pressure itself.

Each of the four boxes in the figure shows three graphs, and, in each box, the continuous line denotes a base case that is common to all four boxes. The broken lines show corresponding solutions with greater and smaller values of the parameter under consideration. The upper row shows the influence of (a) the grid size and (b) the convergence tolerance used to control the inbuilt iterative processes. By inspection, numerical errors due to these causes are so small that it is difficult to distinguish between the curves.

The same is true for only two of the curves in Figure 8c, which assesses the sensitivity to the extrapolation process to which attention has been drawn above. In the base case, a simple linear extrapolation of the relevant flow parameter from values at points M and B in Figure 6 is used to estimate the required value at point L. In the other two cases, the prescribed difference between the values at B and L is 10 times greater than the difference corresponding to linear extrapolation. Using $\phi$ to denote any extrapolated parameter, the values used at the extrapolated point L in Figure 6 above satisfy

$$\mathrm{Extrapolation}:{\phi }_L-{\phi }_B=\alpha \left[\left({\phi }_B-{\phi }_M\right)\left\{\left(x_L-x_B\right)/\left(x_B-x_M\right)\right\}\right]$$

(3)

in which $\alpha$ is a numerical amplifier. For simple linear extrapolation, $\alpha =1,$ and for the cases used to assess the consequences of extrapolation, $\alpha =+10$ and $\alpha =-10$. The author considers these differences to be implausibly large and expects the base case to be the closest of the three to the true numerical solution of the equations.

Less qualitative arguments in support of approximately linear extrapolation are possible. One is based on the simple observation that the dominant process is inertial and that the lateral flows merely modulate it. This is important because the amplified extrapolations magnify the whole change, including the dominant inertial contribution. Another is that errors in the assumed values of the lateral mass flowrate have a self-correcting tendency. Overestimation in one time step will cause the pressure in the chambers to increase too much, and hence the driving force in the next time interval will be smaller than it would have been if the overestimation had not occurred. Notwithstanding these observations, it is not possible to provide a rigorous proof that the extrapolation is 100% safe. Instead, readers must decide for themselves whether the author’s confidence in the presented results is justified.

The final box in Figure 8 illustrates the influence of the numerical method of kick-starting the whole process. Ahead of the wavefront, ambient conditions prevail, and these would persist indefinitely unless some method of triggering a disturbance existed. In the software, this is done by arbitrarily deeming a small pressure change to have occurred and then tracking its consequences upstream. Of necessity, this initial step will not be physically meaningful, and so the chosen value must be acceptably small. Figure 8d compares predicted profiles based on three trigger values. By inspection, differences can be seen, even though the trigger values are in the order of 1000 times smaller than the final amplitude of the wavefront. Nevertheless, the differences are small, and, in particular, the variation in the predicted maximum rate of change of pressure is small. Attention is drawn to this because, as explained in Section 2 above, (∂p/∂t)_max is the parameter of greatest practical importance in the context of pressure disturbances emanating from railway tunnel portals.

5.3. Practical Use of the Methodology

Although the primary purpose of this paper is to introduce the new methodology in its own right, there is merit in presenting a simple example to illustrate its potential practical value. Figure 9 shows predicted histories of the pressure and its rates of change for three cases with identical values of all prescribed values except the asymptotic wave speed. The continuous blue line is the base case used above. For it, the prescribed wave speed is 2 m/s slower than the ambient speed of sound. The other two cases correspond to 1 m/s and 3 m/s deficits, respectively. By inspection, the smaller the prescribed speed difference, the greater the pressure that can be sustained indefinitely. For practical design purposes, however, it is more useful to rephrase this to reflect the parameters that designers and operators can influence directly. At the outset, they can choose a cross-sectional area of the tunnel to limit the maximum pressure amplitude. Later, adjustments can be made to the train speed for the same purpose. Then, graphs such as those in Figure 9 can be used to estimate the maximum rate of change of pressure at the tunnel exit portal. The exact speed of the wavefront is of no direct practical interest in this context. However, it is the input parameter used in the analysis to assess the influence of pressure on the properties of the asymptotic wavefronts.

For completeness, it is important to recognise that large differences may exist between the amplitude of the initial disturbance and the subsequent amplitude of the wavefront as it nears the tunnel exit. Many factors can influence changes during propagation along tunnels, including, for instance, skin friction, the track geometry and airshafts. To allow for these, designers will still need to undertake simulations of wave propagation along the whole tunnel. However, for the required purpose of estimating the resulting pressure amplitude, it will be sufficient for such analyses to use a coarse grid size—perhaps a few tunnel diameters. Such simulations necessitate vastly smaller CPU demands than simulations designed to give accurate predictions of wavefront shapes—and hence of rates of change of pressure. Instead, the method introduced herein can assess the likelihood of a wavefront with the expected final amplitude being sufficiently steep to causes an unacceptable disturbance outside the tunnel. This complementary use of the two types of analysis has the potential to reduce greatly the resources needed to obtain reliable predictions, especially in the case of long tunnels.

6. Summary and Conclusions

The possibility of wavefronts evolving towards an asymptotic state as they propagate along a duct has been demonstrated, and a method of predicting the range of possible detailed forms of such asymptotes has been developed. The method, which requires minimal CPU resources, is unusual in that it predicts possible outcomes without any need to know what has caused them. Indeed, since the method is tailored expressly to mimic asymptotic conditions, it is wholly unsuitable for studying wavefronts during preceding periods of evolution before asymptotic forms have developed. In this sense, the method complements conventional analyses in which prescribed boundary conditions define initial ‘causes’ and simulations are used to deduce the resulting ‘effects’. In such methods, the continual accumulation of numerical errors during propagation makes them poorly suited to modelling eventual asymptotic states.

The new analysis needs to simulate only short lengths of duct—e.g., a few duct diameters—and only the small times needed for wavefronts to travel such distances. It needs only trivial CPU resources, even when tiny grid sizes are used to achieve high accuracy—e.g., one-thousandth of the duct diameter. The use of such grid sizes in simulations along complete ducts would be impracticable.

The methodology is inherently simple, and it has been shown to be highly stable numerically. Importantly, this is true even for the method’s only obvious drawback, namely a need for two parameters to be determined by means of extrapolation, not interpolation. The consequences of this drawback have been assessed in detail and shown to be minimal, even when the influence of the extrapolations has been magnified by factors of ±10.

The use of this method has been illustrated for the particular case of a duct that has small, enclosed air chambers into which air can flow, thereby modifying the dynamics of airflow along the duct. This configuration has been studied previously for a particular practical purpose, but the simulations have been of uncertain accuracy because of the CPU limitations cited above. It has been found that one-to-one correlations exist between the pressure amplitude of an asymptotic wavefront and (i) its speed and (ii) its maximum steepness. Unique relations such as these do not exist in wavefronts that have not evolved to an asymptotic state.