Numerical investigation on an array of Helmholtz resonators for the reduction of micro-pressure waves in modern and future high-speed rail tunnel systems

Previous research has proposed that an array of Helmholtz resonators may be an effective method for suppressing the propagation of pressure and sound waves, generated by a high-speed train entering and moving in a tunnel. The array can be used to counteract environmental noise from tunnel portals and also the emergence of a shock wave in the tunnel. The implementation of an array of Helmholtz resonators in current and future high-speed train-tunnel systems is studied. Wave propagation in the tunnel is modelled using a quasi-one-dimensional formulation, accounting for non-linear effects, wall friction and the diffusivity of sound. A multi-objective genetic algorithm is then used to optimise the design of the array, subject to the geometric constraints of a demonstrative tunnel system and the incident wavefront in order to attenuate the propagation of pressure waves. It is shown that an array of Helmholtz resonators can be an effective countermeasure for various tunnel lengths. In addition, the array can be designed to function effectively over a wide operating envelope, ensuring it will still function effectively as train speeds increase into the future. & 2017 The Authors. Published by Elsevier Ltd. All rights reserved.


Introduction
Issues related to tunnel aerodynamics and acoustics are expected to become increasingly important with the trend for increasing speeds, combined with a greater use of tunnels in future high-speed rail projects as a result of environmental noise, land disputes and challenging terrain. Tunnels act as a wave guide for pressure waves generated by trains entering and crossing each other. The severity of the pressure disturbances increases with the speed of the train and the train-totunnel area ratio. Issues associated with this include the emission of micro-pressure waves (MPWs) from tunnel portals and pressure loads on the train and tunnel infrastructure. It is also theorised that shock waves could form far ahead of the train in long tunnels, due to non-linear steepening of the pressure wave [1]. This paper investigates the use of an optimised array of Helmholtz resonators embedded in redundant tunnel space as a countermeasure to these issues.
The layout of the paper is as follows: (1) A methodology for discretising and solving the system of equations used to represent the propagation of pressure waves in a tunnel with Helmholtz resonators is set-out and validated against experimental data; (2) A suitable optimisation procedure will be applied to a demonstrative tunnel system, which allows the array to be optimised for a variety of operating conditions, tunnel geometries and constraints; (3) A sensitivity study of the optimised array to different initial conditions is then performed in order to assess its practicality; (4) An optimal configuration will then be compared against inviscid 3-D predictions, in order to test the assumptions of the quasi-1-D model.

Overview of physical phenomena
The region around the train can be divided into the near and far field. The near field is dominated by complex 3-D aerodynamics. Multiple sources of aero-acoustic noise over a large frequency range can be identified, which are related to the geometry of the train (e.g. pantograph and inter-car spaces). High frequencies are rapidly attenuated further away from the train owing to dissipation, resulting in near one-dimensional propagation in the far field. The strongest pressure waves tend to be as a result of train entry into the tunnel. The pressure rise of the entry compression wave can be estimated under steady inviscid flow using [2]:  far-field space is the Mach number of the train, and a 0 and ρ 0 are the speed of sound and density of air at ambient conditions, respectively. The value of p 1 may deviate slightly depending on the streamlining of the nose and the associated stagnation pressure losses [3]. The shape of the entry compression wave is governed by 3-D effects and is strongly influenced by the train nose and tunnel entrance profile. However, the characteristic frequency of the entry wave can be approximated by υ λ = a / 0 , where λ is the characteristic wavelength. The amplitude of the signal is greatest at characteristic frequency, which is estimated to be in the order of 1-10 Hz [1]. Such low frequencies are less affected by dissipation and it is posited that non-linearity over the course of propagation could result in the emergence of shock waves. However, one of the most common problems associated with compression waves is the emission of MPWs from tunnel exits. The amplitude of MPWs is proportional to the rate of change of pressure of the wavefront in the low frequency and far-field approximations [4]: where p 2 is the acoustic pressure reaching the exit of the tunnel and r o is the distance from the exit to the observer. This equation assumes that the MPW radiates into semi-infinite space. MPW countermeasures, such as tunnel entrance hoods and train nose shape optimisation, function by reducing the gradient of the compression wave at the point of generation. This is adequate for short tunnels, where the accumulated nonlinear steepening of the pressure wavefront is likely to be minimal. However, they may be insufficient for preventing unacceptable MPWs for long tunnels with higher operating speeds in the future. For example, the rate of steepening increases with the amplitude of the initial wavefront. Alternatively, ballasted track has been observed to be an effective method for reducing the severity of MPWs by counteracting the wave steepening process in long tunnels through dispersion and friction. It has been proposed that the behaviour of ballast is analogous to an array of Helmholtz resonators [5,6]. However, the use of slab track in modern tunnel infrastructure is often preferred. It requires less maintenance and can reduce overall tunnel constructions costs (by around 30% [7]), while it avoids damage to the underside of the train and railheads due to ballast being thrown up at high speeds [8]. It has been proposed that an array of Helmholtz resonators could be connected to a tunnel [1], combining the benefits of slab track with the dissipative properties of ballast.

Physical parameters
The system studied consists of a tunnel connected via a neck to an array of identical equispaced Helmholtz resonators (see Fig. 1). The geometric parameters of the resonator array are the hydraulic diameter d; length l; and height h. The radius, r is expressed as half of the hydraulic diameter. Derived parameters are volume V and cross-section area S. Subscripts are used to denote the cavity, c; neck, n; tunnel, t; and the train, z. The calculated physical parameters of the model are the dissipation in the boundary layer ; the density at equilibrium; the linear speed of sound; the characteristic angular frequency ω; the diffusivity of sound ν d ; the linear acoustic pressure in the tunnel p t ; and the Mach number, which are defined by:  where  is the specific gas constant, T 0 is the ambient temperature, γ is the ratio of specific heats, Pr is the Prandtl number, μ μ / d is the ratio of shear and bulk viscosities, ν is the kinematic viscosity and u is the flow velocity in the tunnel.

Near-field model
The modelling approach follows the work of Sugimoto [1] and Lombard et al. [9] and is valid assuming: plane wave propagation under the low frequency approximation (φ < 1); weak acoustic non-linearity (ϵ ≪ 1); negligible interaction between neighbouring resonators (κ ≪ 1); averaged geometric properties, i.e. a continuum distribution of resonators (σ ≪ 1); lumped element approximation for the resonator response (ϖ ≪ 1). These non-dimensional parameters are defined by: where υ c is the first non-planar cut-on frequency for the tunnel section and χ is the cavity spacing. The cut-on frequency for a circular duct [10] is used for simplicity, as it was found to be within ∼10 % of the value predicted from a 2-D eigenvalue analysis for the tunnel section in Fig. 1. For ϵ it is assumed that The assumption of negligible interaction between the resonators (κ ≪ 1) allows the size of the equation system to be reduced by disregarding reflected waves propagating in the upstream direction (i.e. ≪ − + p p ), which are not of interest for the purpose of this work. The assumptions ϵ ≪ 1 and κ ≪ 1 are expected to be satisfied for most operational scenarios. For example, ϵ ≈ 0.14 for the extreme case of = u 600 z km/h and β = 0.3. The conservative approximation κ ≈ ≪ S S / 1 c t applies for the geometry in Fig. 1, assuming approximately equispaced and identical cavities which are placed end-to-end, and that the available cross-section space for the cavities will be relatively small (see Appendix A). The continuum distribution approximation simplifies the discrete response of each resonator to a spatially averaged distribution along the length of the array. As a result, the reflection and interaction of waves by each resonator is not modelled. This is acceptable for the purpose of this work, provided the degree of reflection is small (i.e. κ ≪ 1). The lumped element approximation parameter ϖ has been introduced as the characteristic wavelength may be comparable to the resonator dimensions, in particular the cavity length, l c . This assumption simplifies the spatio-temporal response of each resonator to a time dependent ODE, which ignores the finite propagation time of the pressure waves in the cavity. A sufficiently large value of ϖ would indicate a phase error between the assumed and actual response of the resonator. The value of the parameters in Eq. (4) for which the assumptions no longer hold largely depends on the required accuracy of the simulation, which could be found through comparison with experiments or 3-D models. Further consideration of this is provided in Section 6.
The system consists of two coupled fractional order differential equations: a non-linear PDE describing the forward propagation of acoustic waves in a tube and a non-linear ordinary differential equation (ODE) describing oscillations within a Helmholtz resonator [1]: where p c is the pressure in the resonator cavity. Non-linear advection is accounted through coefficients A and B. Unsteady viscous losses due to a laminar boundary layer are incorporated through coefficient C. It is noted that this does not incorporate the effect of steady wall friction due to surface roughness (see [6]), which plays a role in the attenuation of pressure waves in real slab track tunnels. Therefore, Eq. (5a) will under-predict the attenuation of compression waves in long tunnels with rough walls. Coefficient D accounts for the dissipation due to the diffusivity of sound. The non-linear resonator response accounts for the unsteady viscous boundary layer in the neck of the resonator and adiabatic processes in the cavity through F and H, respectively. Jet losses at the resonator neck ends are incorporated through coefficient I. Coupling between the two equations can be removed by setting E¼0, which implies no resonators are attached. The coefficients are defined by: where * r t is the reduced tunnel radius, which accounts for the effect of a continuous distribution of resonators on the boundary layer: where r c is the radius of the cavity. Similarly, the neck-tunnel correction length is defined for a flanged end by [13]: The Riemann-Liouville (RL) integral with fractional order α = 1/2 is given by: where derivatives of orders 1/2 and 3/2 are found by differentiating Eq. (12) once and twice respectively, with respect to t. The non-local nature of the fractional derivative/integral increases the computational cost of the problem as the entire pressure-time history must be integrated.

Far-field model
This study is only interested in the local distortion of the initial pressure wave over the course of the array. Therefore, a change in the reference frame to one moving with the ambient speed of sound (called the far field from here-on) can be used to avoid the large domain size of the near-field model over large distances. The resulting system is written in ( )   , coordinates, where  and  are non-dimensional space and time variables, respectively: The non-dimensional variables f and g, which are O(1) are also used:   where δ t and δ n are the ratios of the boundary-layer thickness to the tunnel and neck radius, respectively. The linear resonator response is recovered when = = U W 0 in Eq. (15b) and the viscous neck length corrections are discounted. A dimensional analysis of the dissipative terms β δ , t , and δ n has been studied in the linear regime [14,1], where the effect of β was found to be negligible. The system is mainly governed by the ratio of geometric dispersion K and the tuning parameter Ω.

Numerical treatment
The treatment of the far-field equation system in Eq. (15) is considered. We are interested in efficient treatment of the fractional derivatives and application of accurate shock capturing schemes, whilst minimising computational cost in order to make the scheme suitable for optimisation. MATLAB [15] is used for the 1-D numerical treatment and optimisation procedure. However, alternative codes with support for differential algebraic equation (DAE) solvers and derivative free optimisation can also be used.

Fractional derivatives
A diffusive representation of the fractional derivatives has been used, as in Ref. [9]. This is valid for Caputo-type fractional differential operators of the order α > 0, α ∉  and assuming [16]. Therefore, the diffusive representation cannot be applied for a theoretical discontinuity, which cannot physically occur. This presents a more numerically efficient approach, compared to integrating the entire time-history using the RL integral. The diffusive representations for the derivatives of order 1/2 and 3/2 can be approximated using a quadrature scheme with N l integration points: where η l and θ l are the quadrature weights and abscissa, respectively. The diffusive variables ψ and ξ satisfy the first order ODEs: The values of η l and θ l can be found using a suitable quadrature scheme for improper integrals, such as Gauss-Laguerre, or Gauss-Jacobi after a change in integration limits [9]. Alternatively, Richoux et al. [17] use an objective function based on the linear dispersion relation of Eq. (5a) with a non-linear optimisation procedure and found it to be more accurate for the same N l , compared to the aforementioned quadrature schemes. We adopt the latter approach using the fmincon optimisation function with sequential quadratic programming. Initial values for η l and θ l are provided using modified Gauss-Jacobi quadrature, while the constraints θ ≥ 0 l , η ≥ 0 l and θ θ ≤ l max are also imposed to ensure positivity and prevent divergence, respectively.

First order system
The method of lines provides a high degree of flexibility in terms of temporal and spatial discretisation approaches. This requires Eq. (15) to be written as a system of first-order partial differential algebraic equations. Expanding the second order derivative of g 2 in Eq. (15b) using the product rule and substituting-in Eq. (17)-(20) gives after some rearrangement: The advective flux can be treated using an appropriate shock capturing scheme. For example, Total Variation Diminishing (TVD) schemes [18] are effective for preventing the emergence of non-physical oscillations, at the cost of degenerating to first order accuracy near smooth extrema. Alternatively, Weighted Essentially Non-Oscillatory (WENO) schemes [19] can achieve uniformly high order accuracy in smooth regions (including at smooth extrema) by relaxing the TVD requirement, whilst maintaining essentially non-oscillatory shock transition. The fifth order accurate WENO-Z þ scheme [20] is used, which is less dissipative and more accurate in smooth regions than the original WENO scheme [19]. A comparison of TVD and WENO schemes is provided in Appendix B for a steep wavefront. Other derivatives are discretised using 5th order backwards differences. At least one upstream stencil point was found to be necessary to maintain stability, particularly for large Δ . However, this now requires an upstream boundary condition (see Section 4.5).

Space marching
The first order system requires a numerical solver which can handle semi-explicit DAEs of index 1. The MATLAB solvers ode15s (semi-explicit index-1 DAE) and ode15i (fully implicit DAE) [21] are used for this paper. These are based on the variable-step variable-order (1-5) backwards differentiation formulas, which are well suited for stiff problems. ode15s was generally found to exhibit faster convergence, while ode15i was less likely to fail in the presence of stiffness. The Jacobian sparsity pattern is provided in order to accelerate the solvers. This is beneficial when using the diffusive representation of the fractional derivatives, which significantly increases the number of variables. An analytical sparsity pattern is difficult to implement for this system and instead is estimated prior to running the solver based on the numerical calculation of the Jacobian. This approach is valid as long as the sparsity pattern does not change with the space step. = − N 8 10 l was found to provide a good compromise between the desired accuracy of the fractional derivatives, versus solver speed. Relative and absolute integration tolerances of − 1e 3and − 1e 6, respectively were found to provide adequate spatial convergence over a range of cases.

Generation of initial conditions
The initial pressure waveform is applied at ( ) =  f f 0, 0 after performing a change to the far-field variables using Eqs. (13) and (14). Initial conditions for the remaining variables are found by solving the system of ODEs representing the resonator response (Eqs. (21b)-(21e)) by time-marching, subject to ( ) =  f f 0, 0 and zero-valued starting points at =  0 for all other variables. These values are used as a good initial guess for the generation of consistent initial conditions for implicit and semi-explicit DAE solvers.

Boundary conditions
A disadvantage of the far-field regime is the complexity of implementing physically representative boundary conditions (BCs). Fortunately, for the purpose of optimisation this work is only concerned with the distortion of the initial wavefront and not in the effect of wave reflections from the tunnel ends. All variables are null-valued at =  0, assuming undisturbed air at ambient conditions downstream of the wavefront. The solution is not affected by the upstream boundary at =   end , as only downstream wave propagation is modelled. However, a "dummy" BC is still required here due to the finite difference stencil. The use of backwards differencing for dissipative terms helps to localise boundary effects at the upstream boundary. An example of the effect of fixed-value (FV) and zero-gradient (ZG) upstream BCs is shown in Fig. 2, for an arbitrary array configuration. The FV BC takes the value at the boundary for the initial condition (e.g.
. Boundary effects can be minimised by placing the boundary sufficiently far from the wavefront. Zero-gradient boundary conditions were chosen for the upstream boundary for all variables.

Experimental validation
The numerical model has been validated against experimental data provided by Olivier Richoux et al. [17]. The test rig consists of a wave guide connected to a lattice of equi-distant and identical resonators. The measured input signal at x¼ 0 (see Fig. 3a) is sampled at υ = 20 kHz s using a spline fit, in order to be used as the initial condition, f 0 to Eq. (21). Predictions from the far-field model are then time-shifted to account for the distance between the input signal and measurement position. The time shift is estimated as: where Δx i is the space step size and N x is the number of space steps. The results are normalised by the magnitude of the initial condition p s . Fig. 3b shows that the initial disturbance evolves towards a shock without the presence of the resonator array, due to the process of non-linear steepening. Viscothermic losses are insufficient to prevent this from occurring. In contrast, Figs. 3c and 3d show that the wavefront has evolved into a soliton with the presence of the array, due to the cancellation of dispersive effects and non-linear steepening. The presence of solitons is confirmed by Sugimoto [1] in the inviscid regime with Ω ≫ 1 in Eq. (15), which gives the Korteweg-de Vries equation. A shock wave can no longer propagate in this case. The parameters in Eq. (4) for this system are: φ = 0.08, ϵ = 0.32, κ = 0.5, σ = 0.09 and ϖ = 0.07, with υ ≈ 300 Hz from the power spectrum of the initial disturbance. Good agreement with experiments is obtained even with relatively large values for κ and ϵ, which are larger than what is expected for high-speed rail operations in tunnels.

Application to a demonstrative system
In this section the developed model is applied to a demonstrative high-speed train-tunnel system in order to provide some indication about the effectiveness of such an array, subject to different operating regimes.

Array specification
It is envisaged that redundant space in the tunnel could be used to house the resonators. For example, hollow prefabricated slab sections could be slotted into the tunnel at certain intervals, or over a stretch of the total tunnel length. The cavities could be positioned underneath walkways, or under the track given sufficient space. The demonstrative tunnel cross-section used for this analysis is shown in Fig. 1, based loosely on the type of single-bore tunnel sections planned for stretches of the High Speed Two (HS2) route between London and Birmingham [22]. The optimisation variables selected for this study are shown in Table 1, where the tilde accent is used to denote a nondimensional quantity. The minimum of the cavity height h c and length l c is set as the maximum possible neck diameter. The volume of the cavity is then = − V S l S l c cc nn , where S n is scaled to account for the number of neck holes N n for each cavity. The total length of the resonator array is then the sum of the cavity lengths. No additional spacing between the cavities is used for the optimisation procedure in this paper, i.e. χ = l c .

Initial condition
The initial condition for this system represents a typical wavefront generated by a high-speed train entering a tunnel. This qualitatively resembles the form shown in Fig. 4, which is based on experimental measurements [23,24]. The pressure where p 1 is the difference between the asymptotic pressures at = ± ∞ x and L is a length used to characterise the maximum steepness of the wavefront (not to be confused with the characteristic wavelength, λ). Considering linear wave propagation, the substitution ∼ − x a t 0 is made in order to define Eq. (24) on time. Assuming negligible distortion of the wavefront between the tunnel entrance and the start of the array, p 1 can be approximated using Eq. (1) and L can be approximated by [26]: where ϕ is an empirical coefficient typically in the range ϕ < < 0.75 1.25, which accounts for the nose and tunnel entrance shape. The value of ϕ can be much larger for modern trains with long noses and entrance hoods. In this study a conservative approximation of˜=̇∼ L L d / 1 t is used, which is shorter than typically seen in tunnels with long train noses. Therefore, additional numerical experiments should be performed for larger values of L, which are representative of modern trains. It is noted that Eq. (24) is equivalent to the expression used in Ref. [24] with (∂ ∂ ) ≈ ( ) p t p a L / / max 1 0 and similar to the form originally presented by Yamamoto [4], with ϕ ∼ 1 and ≪ M 1 z . To enable Eq. (24) to be used as an initial condition its range is restricted to − < < b x b and it is scaled so that ( − ) = p b p 1 and ( ) = p b 0 (where =ḃ d 10 /3 t ), in-line with previous work by Aoki et al. [25]. Since this signal is comprised of many frequencies the characteristic wavelength is approximated by λ ≈ ( ) L b mean , 2 . The grid resolution is given by υ Δ < t 1/4 L , where υ = a L / L 0 , in order to sufficiently resolve the steepest part of the wavefront.

Optimisation strategy
The design of a resonator array is a complex non-linear optimisation problem, where the derivative of the objective function cannot be derived analytically. Therefore, a derivative free and non-linear multi-objective algorithm which can handle mixed-integer constraints is required. The multi-objective Genetic Algorithm (MOGA) based on the non-dominated Sorting Genetic Algorithm-II (NSGA-II) [27] has been chosen from the MATLAB Optimisation Toolbox. This allows a globally Pareto optimal set to be obtained for two or more objective functions. The Pareto optimal represents the best known solution with respect to all objectives and can not be improved in one objective without worsening another. If our goal is to reduce the amplitude of MPWs generated at the tunnel exit, then according to Eq. (2) we should aim to minimise the maximum gradient of pressure waves. Therefore, the objective function is chosen as the ratio of the maximum absolute rate of change of pressure across the wavefront, with and without the presence of the array: For example, a 50% reduction in the maximum pressure wave gradient ( = Obj 0.5) approximately corresponds to a 6 dB reduction in SPL of the MPW. The second objective function has been chosen as the total array length normalised by the diameter of the tunnel bore,˜=l l d / A A t . For example, a shorter array may be less expensive to install and maintain. Upper and lower bounds have been prescribed for the input variables, based on realistic tunnel geometry constraints (see Table 1). The upper bound onl c is limited by the continuum distribution and lumped element approximations.
The physical and geometric parameters used for the remainder of the paper are shown in Table 2. A sensitivity analysis of the optimisation variables in Table 1 on the objective function has been carried out using the Global Sensitivity Analysis Toolbox (GSAT) [28] in MATLAB. This uses the Fourier Amplitude Sensitivity Test (FAST) to categorise the sensitivity of the objective function by its variance. As shown in Table 3, the neck diameter, number of resonator cavities and the cavity length have the largest sensitivity indexes and so serve as the best parameters to tweak the performance of the array.

Results of genetic algorithm optimisation
The geometric configuration of the array has been optimised for two cases, as shown in Table 4. The first case uses an incident wavefront representative of demonstrative typical operating conditions. However, it is likely that a real array will be exposed to a range of incident wavefronts. For example, train operating speeds may increase after the array is constructed. Alternatively, new rolling stock may be introduced with different nose profiles, which would affect the wavefront shape. Therefore, the array has also been optimised for a range of incident wavefronts, representing an extreme envelope of operating conditions. In this case, the objective function is the average of Eq. (26), evaluated for several initial wavefronts based on various combinations ofL and u z . A margin of ∓100 km/h on u z and 750% onL has been chosen to demonstrate this. Fig. 5 shows the Pareto optimal sets of array configurations after at least 500 generations of the MOGA, for the typical regime and extreme envelope. Each point on the pareto front represents an optimal array configuration with respect to the two objectives. The corresponding resonator dimensions for some of the Pareto points are provided in Table 5. The Pareto fronts diverge after˜≈ l 7 A , where the arrays optimised under the extreme operating envelope can no longer match the performance of the typical regime. This behaviour is expected as the extreme envelope provides a compromise between absolute performance and greater robustness to different incident wavefronts. This is also seen in Fig. 6, where there is greater variability in the pressure gradients for the longer Pareto array configurations under the extreme operating envelope, compared to the typical regime. The wavy behaviour seen without the array in Fig. 6c is associated with the numerical error in the calculation of steep gradients. This behaviour is reduced by refining the grid.
Examples of waveforms corresponding to the smallest objective values from the Pareto fronts are shown in Fig. 7 under typical operating conditions and Fig. 8 under the extreme envelope. The resultant waveforms are compared both with and without the array, at a distance corresponding to the end of the array. It is noted that the non-dimensional time and space scales differ in Fig. 8, due to the different values of L and u z . Non-linear steepening is clearly evident in the waveforms without the array present. In contrast, the resonator array strongly disperses the wavefront and introduces oscillatory Table 2 Physical and constant geometric parameters used for the case study, with reference to Fig. 1.  Table 3 First order global sensitivity coefficients (%) for the optimisation variables and bounds in Table 1  behaviour in its wake, corresponding to the natural frequency of the resonators. However, the amplitude of the oscillations is much less than that of the initial pressure rise. The sensitivity of the longest Pareto optimal array configurations (˜≈ l 24 A ) in Fig. 5 to the shape of the incident pressure wavefront has been tested for the typical regime and extreme envelope. The objective function has been evaluated over a range of train speeds (acoustic pressures) and characteristic wavefront lengths for the two optimal array configurations, as shown in Fig. 9. Both of the optimal array configurations exhibit robust performance over the range of speeds and characteristic wavelengths considered. However, the performance of the array optimised under the typical operating regime is more sensitive to the incident wavefront shape, judging by the greater range in objective function values. The array optimised under the extreme operating envelope exhibits a slightly higher mean objective value (0.48 vs. 0.47) over the range of velocities and wavefront lengths. However, it is more robust, with a smaller standard deviation from the mean value (4% vs. 13%). Greater robustness for the extreme operating envelope may be obtained at the cost of increased computation time by including more optimisation points (˜) L u , z , for example at ( ) 0.5, 300 , ( ) 1.5, 500 and ( ) 1.0, 400 .

3D modelling
The far-field model has been validated against experimental data in Section 4.6. However, it is unfeasible to validate it against a real high speed train-tunnel system. Instead, the 1-D predictions have been benchmarked against 3-D CFD in order to provide some insight into the validity of the underlying assumptions. OpenFOAM 3.0.x [29] has been used to mesh and solve the 3-D models.

Modelling overview
The propagation of pressure waves in the tunnel is modelled using the semi-implicit solver pisoCentralFoam [30], which blends the compressible pressure based PISO algorithm [31] and density based central-upwind scheme of Kurganov & Tadmor [32], depending on the local Mach number and CFL condition. The hybrid scheme can operate efficiently over a wide Table 5 Variation of optimised resonator dimensions (in metres) with non-dimensional array length for the typical operating regime (a) and extreme operating envelope (b), based on the Pareto front in Fig. 5.  range of Mach numbers, compared to using either scheme alone. The TVD condition is enforced using the van Leer slope limiter. Crank-Nicolson and unbounded 2nd order backwards temporal schemes were found to introduce non-physical pressure oscillations into the solution. Instead, a blended (1st -2nd order) Euler-Crank-Nicolson scheme is used.   Table 4.
Non-physical pressure oscillations were not apparent with a blending factor of less than 0.4 (0 ¼pure Euler). The models are solved without turbulence and viscosity, as this is expected to have a negligible effect on the solution. For example, an acoustic pressure drop of less than 0.5 % is predicted due to friction over the length of the array (˜< l 7 A ) by using the Colebrook-White equation [33], for a concrete roughness height of 1.0 mm.

Model set-up
The 3-D models have been created in the SolidWorks CAD software based on the geometry in Fig. 1 and meshed using the snappyHexMesh tool in OpenFOAM. All internal walls have a nominal thickness of 0.05 m. The maximum cell size is given by Δ = x L /10 max , which provides at least 10 cells over the steepest part of the initial condition (see Fig. 4). Additional mesh refinement is prescribed around the resonator necks and tunnel walls (see Fig. 10), with a minimum cell size of Δx /128 max . A mesh convergence test for Case 4 (see Table 6) indicated that this provided sufficient spatial accuracy.
The scaled form of the initial condition (see Section 5.2) is defined as an initial field in the model for pressure ( + p p 0 ), in a region just upstream of the array (see Fig. 11). The corresponding velocity field is given by the linear approximation ρ ( ) = ( ) u x p x a / 0 0 , assuming no disturbance has propagated upstream. The time-step size is given by:    was selected in order to encompass the natural frequency of the resonators and characteristic frequency of the initial condition.
Non-reflective boundaries conditions (NRBCs) are used both up and downstream from the array for all fields, using the waveTransmissive patch type. This assumes that there are no reflections from the far-field into the working section, i.e. an infinitely long and smooth tunnel outside of the boundaries. A direct comparison with the 1-D model cannot be made, as the boundaries are prescribed on time (see Section 4.5). However, the 1-D far-field model is analogous to spatial boundaries which are moving at approximately the same acoustic speed as the wavefront, and so the wavefront will never reach the downstream boundary. Zero-gradient and slip boundary conditions are assigned for the pressure and velocity at the walls, respectively. The linear equation system is solved using the iterative diagonal incomplete-LU preconditioned Bi-conjugate gradient stabilised solver (PBiCGSTAB) for all fields, except density which is solved explicitly. Mesh skewness and non-orthogonality are corrected for by evaluating gradients using the method of least-squares and applying one non-orthogonal corrector on the pressure equation [31,34]. The models, with between 1-4 million cells each were solved in parallel on an 8-core Intel Xeon processor. Table 6 lists the array configurations used for the comparative study. The resonator geometry and incident wavefront shapes have been chosen in order to test the 1-D modelling assumptions when applied to a demonstrative tunnel system, and are not based on the optimised geometry in Section 5.4. In particular, the lumped element (ϖ) and continuum distribution (s) approximations are tested for Cases 1 & 2 by varying the cavity length and inter-cavity spacing. Cases 3 & 4 are geometrically identical and compare the effect of the steepness of the initial wavefront through parameter L. The approximations ϵ ≤ ≪ 0.04 1 and κ ≤ ≪ 0.04 1 are satisfied for all cases. They can be tested as part of future work by increasing the wavefront amplitude and the cross-sectional area of the cavities, respectively. The tunnel geometry and physical parameters are shown in Table 2.

Comparison of 1-D and 3-D models
Power Spectral Density plots of the initial pressure waveforms used in the 3-D analysis (see Fig. 12) indicate that a significant proportion of the signal power lies below both the approximate characteristic frequency and the first non-planar cut-on frequency. Therefore, the continuum distribution (s), lumped element (ϖ), and non-planar cut-on frequency (φ) parameters (calculated using Eq. (4)) are conservative approximations for the incident pressure waveform used in this paper.
The predicted waveforms at the end of the array are compared in Fig. 13, accounting for the 1-D viscous and inviscid predictions. The 1-D inviscid predictions are found by setting ν to a very small number (e.g. ν = − 1e 50). The 1-D waveforms are first changed to the near-field variables using Eq. (13) and (14) and then time-shifted using Eq. (23). The 3-D waveforms  are obtained from the area weighted average pressure across the downstream boundary at each time step (see Fig. 11). The difference between the maximum and minimum acoustic pressures across the downstream boundary at a given time step generally does not exceed ∼15%, indicating that the pressure wave is largely planar.  There is negligible difference between the 1-D viscous and inviscid waveforms. Good agreement is obtained between the 1-D and 3-D predictions for Cases 3 and 4. The amplitude of the upstream oscillations in Case 4 are comparable to the initial pressure rise. In this case, the use of a variable cavity length or neck radius, and hence a variable natural frequency could be used to break-up these oscillations into several smaller peaks. σ ≈ 0.7 for Cases 1 and 2, which would indicate that the continuum approximation is no longer valid. However, good agreement between the 1-D and 3-D predictions is still obtained for Case 2. On the other hand, poor agreement is obtained for Case 1, which also violates the lumped element approximation. This can be clearly observed in the 3-D model for Case 1 (see Fig. 14). The comparable cavity and characteristic wavefront length gives rise to a non-negligible wave propagation time in the cavity, and hence an apparent delay in the resonator response, as seen in the pressure waveform in Fig. 14e. This appears to result in a shift between the 1-D and 3-D predictions for the upstream oscillation in Fig. 13a.
These results indicate that the accuracy of the 1-D model is more sensitive to the lumped element approximation, compared to the assumption of a continuum distribution of resonators. However, the effect of violating the continuum distribution approximation may become more pronounced by increasing the cross-sectional area of the cavities and hence the value of κ. Removing the lumped element approximation in the 1-D model would be computationally costly. A simpler approach is to limit the maximum cavity length to some fraction of the characteristic wavelength. These results demonstrate that the 1-D model can be used as a reliable tool for rapidly predicting the effect of an array of Helmholtz resonators with complex geometries on a pressure wave in a full-scale tunnel system.

Conclusions
The use of an array of Helmholtz resonators embedded in redundant tunnel space was numerically investigated in this paper as a countermeasure for the emission of micro-pressure waves into the environment, due to the entry compression wave generated by a high-speed train. A multi-objective Genetic Algorithm was used with a non-linear 1-D model of the array to generate a Pareto-optimal set, which allows the length of the array to be traded-off against its effectiveness. The effectiveness is measured by the reduction in the maximum rate of change of pressure of the resultant wavefront, compared to an equivalent distance without the array. The geometric constraints used for the analysis are approximately representative of the type of single-track double-bore tunnels planned for the HS2 link between London and Birmingham. The Pareto front indicates that a ∼100 m long array could provide up to a 50 % reduction in the maximum rate of change of pressure for a steep incident wavefront representative of a short nosed train entering a tunnel with no entrance hood at 400 km/h. This approximately corresponds to a 6 dB reduction in the sound pressure level of the resultant micro-pressure wave. The array can avoid the emergence of shock waves, but it is still insufficient in this case for suppressing the emission of unacceptable micro-pressure waves. Therefore, the combination of an array of Helmholtz resonators and existing countermeasures should be considered.
A sensitivity study indicated that the arrays could be optimised to work effectively over a wide range of incident pressure waves, future-proofing the array for increases in train speeds. Using the modelling framework developed in this paper, further consideration could be given to: the effectiveness of an array for a greater range of incident wavefront lengths; how an array compares with other countermeasures (such as an entrance hood) in order to counteract the emergence of shock waves and meet MPW requirements, particularly for long tunnels.
A justification is provided for the use of WENO instead of TVD schemes for the treatment of the advective flux (see Section 4.2). Fig. B.15 compares predictions for a steep incident wavefront (i.e.˜= L 0.1 so that ( ) > p t d /d 100 max kPa/s) using the 5th order accurate WENO Zþ scheme [20] and the 2nd order TVD scheme of Kurganov & Tadmor [32] with the van Leer flux limiter. In this case, the WENO scheme is able to resolve the steep gradients more accurately than the TVD scheme across a range of sample frequencies.