Numerical Study of the Unsteady Aerodynamic Performance of Two Maglev Trains Passing Each Other in Open Air Using Different Turbulence Models

The strong change in the flow fields around two maglev trains (MTs) passing each other in open air may affect their manoeuvrability and passengers’ comfort. In this study, we evaluated the aerodynamic performance of two MTs passing each other via shear stress transport (SST) k–ω model and improved delayed detached eddy simulations based on the Spalart–Allmaras model (SA−IDDES) and the SST k–ω model (SST−IDDES). The accuracy of the numerical simulation method was verified using experimental data acquired from a moving model test. The results showed that the difference in the amplitude of the transient pressure obtained with the different turbulence models was less than 5%. The wake vortex structures on the intersection side were found to interact, and their intensity consequently decreased. The SST−IDDES model produced smaller-scale vortices than the SA−IDDES model, particularly in the near-wake region. There were large differences in the drag and lift forces obtained using the different turbulence models. Among them, the lift force of the tail car was more sensitive to the turbulence model, and its maximum value obtained with the SST−IDDES model was 11% larger than that obtained with the SA−IDDES model.


Introduction
It is common for trains travelling along a double-track railway line in open air to pass each other. As the trains pass each other, the flow fields around them are strongly disturbed, and a strong pressure wave is formed on the intersection side of the train surface, particularly when the head and tail cars pass. Furthermore, a lateral aerodynamic force is generated during this process, causing the trains to shake. This affects the safe operation of the trains, jeopardising the comfort of passengers [1][2][3].
Currently, single high-speed trains (HSTs) travel at speeds higher than 300 km/h (their relative speed exceeds 600 km/h), which exacerbates aerodynamic problems when two trains pass each other. For example, the magnitude and rate of change of the pressure pulse, side force, and overturning moment increase dramatically. Therefore, to ensure the safety of two trains passing each other, full-scale tests [4,5], moving model tests [6][7][8], and computational fluid dynamics (CFD) simulations [2,3,9] have been performed to investigate the aerodynamic performance of the trains. It was found that parameters including the train speed, line spacing, nose shape, and marshalling length affect the aerodynamic performance of the passing trains. Moreover, the operation surroundings of two trains passing each other under crosswinds are more severe than those without or with moderate ambient wind conditions. Hence, the aerodynamic performance of the two trains passing each other in crosswinds was investigated in detail [10,11].
passing each other under crosswinds are more severe than those without or with moderate ambient wind conditions. Hence, the aerodynamic performance of the two trains passing each other in crosswinds was investigated in detail [10,11].
The aforementioned studies mainly focused on the aerodynamic performance of two HSTs passing each other in open air. In recent years, maglev trains (MTs) have been recognised as safe, comfortable, low energy-consumption, and environmentally friendly ground vehicles [12,13]. Germany, Japan, the UK, South Korea, and China have developed MT technologies. Similar to HSTs, MTs also have a large length-to-width ratio body. However, there are significant differences between HSTs and MTs, as shown in Figure 1. On the one hand, the operating boundary conditions near the ground are different as the HST wheels are in contact with the rail, and the space at the bottom of the train is open. Meanwhile, MTs levitate without contact with the ground, and the space at the bottom of the train is relatively confined. On the other hand, the geometry of the body surface is different as the HST contains bogies, cavities, pantographs, etc., whereas the external surface of MTs is smooth. These differences may cause the flow fields around these trains to be notably different. In addition, the MT speed is higher than the HST speed. For example, the Shanghai commercial maglev line reaches 430 km/h. The effect of time-varying lift forces on the HST is almost negligible as two trains pass each other because of the large mass [14]. However, the variation of the lift force during the passage of MTs can impair the levitation stability. This issue has motivated several investigations of the aerodynamic performance of two MTs when passing each other. Huang et al. [1] used the unsteady RNG k-ε turbulence model to investigate the transient pressure and slipstream velocity induced by two high-speed MTs meeting at a speed of 430 km/h in open air. Chen [15] analysed the pressure pulse and aerodynamic force of two high-speed MTs passing each other using the SST k-ω model and the overset mesh method. Zhou et al. [16] utilised the two k-ε equation models to simulate the aerodynamic performance of two MTs passing each other and improved the aerodynamic performance of them by increasing the streamlined nose length and optimising the body shape.  Although the aerodynamic performance, including the transient pressure and the aerodynamic force, of two MTs passing each other in open air has been well investigated, relatively few studies have been conducted on the evolution mechanisms and basic laws of the airflow around the train, especially the unsteady wake flows. Wake flows when a single train runs in open air have also been investigated [17][18][19]. The surrounding flow of two trains passing each other in open air is different from that of a single train moving in Appl. Sci. 2021, 11,11894 3 of 20 open air. The present study aims to explore the flow around two MTs passing each other and to identify the vital flow mechanisms for their aerodynamic performance.
To the best of our knowledge, few studies have considered the effect of the turbulence model on the aerodynamic performance of two meeting MTs. In previous studies, the unsteady Reynolds-averaged Navier-Stokes (RANS) model with large time steps has been used, wherein the damping of the fine-scale turbulent fluctuation of the flow was neglected. Such RANS simulations result in poor predictions of separated flows [20,21]. In this study, we adopted a hybrid RANS/large eddy simulation (LES) approach to capture the smallscale vortex structures around the trains. This approach has been widely employed to capture wake vortex structures and is highly effective [12,18,19,22,23]. We selected the improved delayed detached eddy simulation (IDDES) model based on Spalart-Allmaras (SA), denoted as SA−IDDES, and shear stress transport (SST) k-ω, denoted as SST−IDDES. The SST k-ω model was used to compare the unsteady RANS results. The main objective of this study was to compare the differences in the aerodynamic characteristics of MTs simulated by the SST k-ω, SA−IDDES, and SST−IDDES methods and to reveal the flow mechanisms around two trains passing each other along a double track. The remainder of this paper is organised as follows: Section 2 introduces the methodology, including the geometry model, computational domain, numerical setup, mesh generation, and mesh sensitivity; Section 3 presents the numerical results and analysis, including experimental validation, transient pressure, velocity profiles, wake flows, and aerodynamic force; and Section 4 presents the conclusions.

Geometry Model
To ensure consistency with the scale of the moving model test, decrease the number of the computational cells, and meet the y + value demand of IDDES, two 1/16th scale MT models were used in this study. The marshalling length of MTs in operation is currently 2-3 cars to achieve rapid transit of the MT and reduce energy dissipation [22]. Therefore, consistent with the MT model by Luo et al. [12], a two-car grouped model per train, namely the head and tail cars, was used in this study, as presented in Figure 2. The train height H, defined as the characteristic length, was 0.24 m (3.84 m at full scale). The total length and width of the MT model were 8.36 H and 0.75 H, respectively. The line spacing is 1.04 H. Moreover, the clearance between the train and track was extremely small (0.005 H), which may have had a greater impact on the aerodynamic performance of the train, specifically in the region near the bottom of the train, as illustrated in Figure 2. Hence, the track in the numerical simulation could not be neglected.
of the airflow around the train, especially the unsteady wake flows. Wake flows when a single train runs in open air have also been investigated [17][18][19]. The surrounding flow of two trains passing each other in open air is different from that of a single train moving in open air. The present study aims to explore the flow around two MTs passing each other and to identify the vital flow mechanisms for their aerodynamic performance.
To the best of our knowledge, few studies have considered the effect of the turbulence model on the aerodynamic performance of two meeting MTs. In previous studies, the unsteady Reynolds-averaged Navier-Stokes (RANS) model with large time steps has been used, wherein the damping of the fine-scale turbulent fluctuation of the flow was neglected. Such RANS simulations result in poor predictions of separated flows [20,21]. In this study, we adopted a hybrid RANS/large eddy simulation (LES) approach to capture the small-scale vortex structures around the trains. This approach has been widely employed to capture wake vortex structures and is highly effective [12,18,19,22,23]. We selected the improved delayed detached eddy simulation (IDDES) model based on Spalart-Allmaras (SA), denoted as SA−IDDES, and shear stress transport (SST) k-ω, denoted as SST−IDDES. The SST k-ω model was used to compare the unsteady RANS results. The main objective of this study was to compare the differences in the aerodynamic characteristics of MTs simulated by the SST k-ω, SA−IDDES, and SST−IDDES methods and to reveal the flow mechanisms around two trains passing each other along a double track. The remainder of this paper is organised as follows: Section 2 introduces the methodology, including the geometry model, computational domain, numerical setup, mesh generation, and mesh sensitivity; Section 3 presents the numerical results and analysis, including experimental validation, transient pressure, velocity profiles, wake flows, and aerodynamic force; and Section 4 presents the conclusions.

Geometry Model
To ensure consistency with the scale of the moving model test, decrease the number of the computational cells, and meet the y + value demand of IDDES, two 1/16th scale MT models were used in this study. The marshalling length of MTs in operation is currently 2-3 cars to achieve rapid transit of the MT and reduce energy dissipation [22]. Therefore, consistent with the MT model by Luo et al. [12], a two-car grouped model per train, namely the head and tail cars, was used in this study, as presented in Figure 2. The train height H, defined as the characteristic length, was 0.24 m (3.84 m at full scale). The total length and width of the MT model were 8.36 H and 0.75 H, respectively. The line spacing is 1.04 H. Moreover, the clearance between the train and track was extremely small (0.005 H), which may have had a greater impact on the aerodynamic performance of the train, specifically in the region near the bottom of the train, as illustrated in Figure 2. Hence, the track in the numerical simulation could not be neglected.   Figure 3 shows the distribution of the monitoring points on the train surface and surrounding space. Eighteen points are arranged on the intersecting side of the train surface, as shown in Figure 3a. Among them, fifteen points are distributed in the longitudinal direction of the train at a distance of 0.094 m (0.39 H) from the bottom of the train, and four points are arranged in the vertical direction between the train cars. We arranged twelve stationary measurement points around the track at x = 0 (the intersection point of two head noses), as shown in Figure 3b. The transient pressure was obtained at the monitoring points on the train surface, and the slipstream velocity was obtained at the monitoring points in the surrounding space.  Figure 3 shows the distribution of the monitoring points on the train surface and surrounding space. Eighteen points are arranged on the intersecting side of the train surface, as shown in Figure 3a. Among them, fifteen points are distributed in the longitudinal direction of the train at a distance of 0.094 m (0.39 H) from the bottom of the train, and four points are arranged in the vertical direction between the train cars. We arranged twelve stationary measurement points around the track at x = 0 (the intersection point of two head noses), as shown in Figure 3b. The transient pressure was obtained at the monitoring points on the train surface, and the slipstream velocity was obtained at the monitoring points in the surrounding space.

Computational Domain and Boundary Conditions
As shown in Figure 4, the computational domain consists of three zones: the background domain, denoted as Zone 1, and the overset domains, denoted as Zone 2 and Zone 3. The length, width, and height of the background domain are 93.75 H, 26.04 H, and 13.02 H, respectively, yielding a small blockage ratio, which is defined as the area of the two trains projected onto a plane in the longitudinal direction divided by the cross-sectional area of the background domain, 0.37%, which meets the requirements of CEN (2013) [24]. A pressure outlet boundary condition was employed in the four sides and upper faces of the background domain, and the static pressure P∞ was set to zero. The ground and two rails were set as the stationary walls. The dimensions of the overset domain were 19.53 H (length) × 1.28 H (width) × 1.43 H (height), and all faces were overset boundary conditions. To ensure full development of the wake flow fields, the computational domain was extended by 30.7 H backward from the tail nose of the train. Zone 2 and Zone 3 were moved in opposite directions at 200 km/h. In Zone 2, Train A travelled in the positive x-direction, while in Zone 3, Train B travelled in the negative x-direction.
The relative speed between the trains was 400 km/h, resulting in a Mach number (Ma = U/c, U is relative speed of two trains, c is the local sound velocity, which is 340 m/s usually) larger than 0.3. According to Niu et al. [25], the compressibility of the flow had to be considered. To obtain a relatively steady flow field before the two MTs met, the initial distance between the head nose of the two trains was set to 15.63 H. Based on the train velocity and train height, the Reynolds number (Re) was 7.47 × 10 5 , which is larger than the minimum Re (0.25 × 10 5 ) recommended by CEN (2013) [24] for the scale-model test.

Computational Domain and Boundary Conditions
As shown in Figure 4, the computational domain consists of three zones: the background domain, denoted as Zone 1, and the overset domains, denoted as Zone 2 and Zone 3. The length, width, and height of the background domain are 93.75 H, 26.04 H, and 13.02 H, respectively, yielding a small blockage ratio, which is defined as the area of the two trains projected onto a plane in the longitudinal direction divided by the cross-sectional area of the background domain, 0.37%, which meets the requirements of CEN (2013) [24]. A pressure outlet boundary condition was employed in the four sides and upper faces of the background domain, and the static pressure P ∞ was set to zero. The ground and two rails were set as the stationary walls. The dimensions of the overset domain were 19.53 H (length) × 1.28 H (width) × 1.43 H (height), and all faces were overset boundary conditions. To ensure full development of the wake flow fields, the computational domain was extended by 30.7 H backward from the tail nose of the train. Zone 2 and Zone 3 were moved in opposite directions at 200 km/h. In Zone 2, Train A travelled in the positive x-direction, while in Zone 3, Train B travelled in the negative x-direction.
The relative speed between the trains was 400 km/h, resulting in a Mach number (Ma = U/c, U is relative speed of two trains, c is the local sound velocity, which is 340 m/s usually) larger than 0.3. According to Niu et al. [25], the compressibility of the flow had to be considered. To obtain a relatively steady flow field before the two MTs met, the initial distance between the head nose of the two trains was set to 15.63 H. Based on the train velocity and train height, the Reynolds number (Re) was 7.47 × 10 5 , which is larger than the minimum Re (0.25 × 10 5 ) recommended by CEN (2013) [24] for the scale-model test.

Numerical Method
Two numerical simulation methods were used to capture complex small-scale vortex structures around a train, i.e., the LES [26,27] and hybrid RANS/LES [12,19,22] methods. However, the former has relatively high requirements for grid quality and computing resources. In this study, we employed the IDDES method proposed by Shur et al. [28] as a hybrid RANS/LES approach to simulate the unsteady aerodynamic performance of two MTs passing each other. In the hybrid method, the RANS method is used to model attached flows in the near-wall regions, whereas the LES method is used to resolve largely separated flows in the remaining computational domain regions [28,29]. Furthermore, this method can effectively solve the existing problems of classical detached eddy simulations or delayed detached eddy simulation models that contain model stress dissipation, gridinduced separation, and logarithmic layer mismatch [28,30]. Consequently, IDDES has been widely used to simulate external flows around a train, particularly small-scale vortex structures [12,18,19,22,23]. To combine the delayed detached eddy simulation and wall modelling LES, the IDDES length scale is redefined as follows: where lRANS, lLES, and lhyb are the lengths of the hybrid, RANS, and LES models, respectively; where fdt and fB are the delay and empirical blending functions, respectively, and fe is an elevating function. We selected the IDDES method based on the SA and SST k-ω models to compare the difference in capturing flow fields based on the one-equation model (SA) and the two-equation model (SST k-ω). In the SA model, is the empirical constant and ψ is a low-Reynolds number correction. In the SST k-ω model, where k is the turbulent kinetic energy and ω is the turbulent dissipation rate. The sub-grid length scale Δ is defined as follows: where the empirical constant Cw is equal to 0.15. hmax is the maximum size of the cell in three directions, or wall distance, and hwn is the mesh step in the direction normal to the wall.
To facilitate the comparison, we also chose the unsteady RANS method to calculate the aerodynamic performance of two MTs passing each other. In this study, the two-equation model SST k-ω was used to solve the RANS equations, with the turbulent eddy viscosity μt calculated as follows:

Numerical Method
Two numerical simulation methods were used to capture complex small-scale vortex structures around a train, i.e., the LES [26,27] and hybrid RANS/LES [12,19,22] methods. However, the former has relatively high requirements for grid quality and computing resources. In this study, we employed the IDDES method proposed by Shur et al. [28] as a hybrid RANS/LES approach to simulate the unsteady aerodynamic performance of two MTs passing each other. In the hybrid method, the RANS method is used to model attached flows in the near-wall regions, whereas the LES method is used to resolve largely separated flows in the remaining computational domain regions [28,29]. Furthermore, this method can effectively solve the existing problems of classical detached eddy simulations or delayed detached eddy simulation models that contain model stress dissipation, gridinduced separation, and logarithmic layer mismatch [28,30]. Consequently, IDDES has been widely used to simulate external flows around a train, particularly small-scale vortex structures [12,18,19,22,23]. To combine the delayed detached eddy simulation and wall modelling LES, the IDDES length scale is redefined as follows: where l RANS , l LES , and l hyb are the lengths of the hybrid, RANS, and LES models, respec- where f dt and f B are the delay and empirical blending functions, respectively, and f e is an elevating function. We selected the IDDES method based on the SA and SST k-ω models to compare the difference in capturing flow fields based on the one-equation model (SA) and the two-equation model (SST k-ω). In the SA model, l RANS = d w , where d w is the distance to the wall, and l LES = C DES ψ∆, where C DES is the empirical constant and ψ is a low-Reynolds number correction. In the SST k-ω model, where k is the turbulent kinetic energy and ω is the turbulent dissipation rate. The sub-grid length scale ∆ is defined as follows: where the empirical constant C w is equal to 0.15. h max is the maximum size of the cell in three directions, or wall distance, and h wn is the mesh step in the direction normal to the wall.
To facilitate the comparison, we also chose the unsteady RANS method to calculate the aerodynamic performance of two MTs passing each other. In this study, the two-equation model SST k-ω was used to solve the RANS equations, with the turbulent eddy viscosity µ t calculated as follows: Appl. Sci. 2021, 11, 11894 6 of 20 Unsteady compressible Navier-Stokes equations were solved using a finite-volume solver STAR−CCM+. The time step ∆t was fixed at 5 × 10 −5 s, which was smaller than that reported previously [1,31]. After the initial development of the flow, the two MTs were operated with a total time step of 2070. For each time step, the maximum number of iterations was set as 30, and the residuals of all physical equations were below 1 × 10 −4 . Less than 1% of cells in the computational domain had a Courant-Friedrichs-Lewy (CFL) number greater than one. According to Krajnovic [32], a minor violation of CFL conditions has a minimal effect on the flow field. In this study, five cases (including three turbulence models, as well as the coarse and fine mesh of the SST-IDDES model) were performed at the National Supercomputing Centre in Wuxi, China.

Mesh Generation
In the numerical simulation, two trains passing each other and the relative motion between the train and other railway infrastructures (e.g., tunnels, wind barriers, and bridges) can be analysed using the dynamic mesh method [33,34], sliding mesh method [1,[35][36][37], and overset mesh method [12,21,38,39]. Nevertheless, the dynamic mesh method needs to regenerate the mesh at each time step, which may increase the calculation time. To capture small-scale flow structures, we often need to conduct local refinement around a train, particularly in the near wake. However, the number of grids required for the sliding mesh method is greater than that required for the overset mesh method with the same grid refinement method. Furthermore, the overset mesh method is suitable for complex and irregular moving models owing to its decomposition and overlapping techniques [39]. Consequently, the overset mesh method was applied in this study.
In the present research, the computational fluid dynamics commercial software STAR CCM+ was used to generate the mesh and for the calculations and post-processing. The computational domain contained two types of mesh: the near-wall of the train surface and rails were discretised with prism layers; the remaining computational domain utilised hexahedral grids. Three sets of grids, namely coarse, medium, and fine grids, were generated to verify the grid sensitivity, as presented in Figure 5a, and their detailed information is presented in Table 1. Here, the growth rate was 1.15 for the three meshes. As shown later in the paper, the medium grid was found to be sufficient for this study. Figure 5b shows a schematic of the medium mesh. To ensure a high resolution around the train and reduce the number of computational nodes, four types of mesh regions were constructed in the computation domain, namely wake refinement, fine, medium, and coarse areas. A similar refined method was also utilised in coarse and fine grids. The two MTs always met in the fine region. The local refinement region in the train wake of length and width of 9.11 H and 1.43 H, respectively, was used to capture small-scale vortices. Figure 6 shows the wall y+ of the train surface of the medium mesh. It can be observed that the y+ value of most regions on the train surface is less than 2 except for some transition regions, where it is still less than 5. Thus, it meets the requirements of the SST k-ω turbulence model.

Mesh Sensitivity
The pressure coefficient C P is defined by Equation (4), where ρ is the air density (1.225 kg/m 3 ), P is the static pressure, P ∞ is the reference pressure (equal to zero), and u tr is the train speed (55.56 m/s). The results of grids with different densities were computed using the IDDES method based on the SST k-ω model to validate the grid sensitivity. The time-history pressure coefficients of monitoring points No. 3 and No. 8 on the train surface for different density grids are presented in Figure 7. The peak-to-peak pressure coefficient (∆C P ) is shown in Figure 7a. For monitoring point No. 3, the difference in ∆C P between the coarse grid and fine grid was found to be 4.1%, and that between the medium grid and fine grid was 0.5%. For monitoring point No. 8, the difference in ∆C P between the coarse grid and fine grid was 3.1%, and that between the medium grid and fine grid was 0.6%. The differences in ∆C P between the medium and fine grids were found to be rather small.

Mesh Sensitivity
The pressure coefficient CP is defined by Equation (4), where ρ is the air density (1.225 kg/m 3 ), P is the static pressure, P∞ is the reference pressure (equal to zero), and utr is the train speed (55.56 m/s). The results of grids with different densities were computed using the IDDES method based on the SST k-ω model to validate the grid sensitivity. The timehistory pressure coefficients of monitoring points No. 3 and No. 8 on the train surface for different density grids are presented in Figure 7. The peak-to-peak pressure coefficient (ΔCP) is shown in Figure 7a. For monitoring point No. 3, the difference in ΔCP between the coarse grid and fine grid was found to be 4.1%, and that between the medium grid and fine grid was 0.5%. For monitoring point No. 8, the difference in ΔCP between the coarse grid and fine grid was 3.1%, and that between the medium grid and fine grid was 0.6%. The differences in ΔCP between the medium and fine grids were found to be rather small.    Figure 3). The nondimensional resultant slipstream velocity U was defined as follows: where u, v, and w are the streamwise, spanwise, and vertical velocity, respectively, and utr is the train velocity. The resultant velocity predicted by three meshes is consistent each other before the tail nose passing, while it shows a large deviation between the coarse and medium as well as fine meshes when wake flow passing.
The iso-surfaces of the second invariant of the velocity gradient, Q, is an effective approach to describe the three-dimensional instantaneous flow field structures around a train. In this study, we used it to investigate the basic law of the interaction of the wake vortex on the intersection side when two MTs pass each other. According to Hunt, et al. [40], the Q can be expressed as follows:   Figure 3). The non-dimensional resultant slipstream velocity U was defined as follows: where u, v, and w are the streamwise, spanwise, and vertical velocity, respectively, and u tr is the train velocity. The resultant velocity predicted by three meshes is consistent each other before the tail nose passing, while it shows a large deviation between the coarse and medium as well as fine meshes when wake flow passing. The iso-surfaces of the second invariant of the velocity gradient, Q, is an effective approach to describe the three-dimensional instantaneous flow field structures around a train. In this study, we used it to investigate the basic law of the interaction of the wake vortex on the intersection side when two MTs pass each other. According to Hunt, et al. [40], the Q can be expressed as follows: where Ω is the rate-of-rotation tensor, and S is the rate-strain-tensor. Figure 9 shows the instantaneous iso-surfaces of Q = 20,000 s −2 of three sets of meshes at the t = 0.08 s position, coloured according to the resultant velocity U. The V 2 − V 2 and V 1 − V 1 regions in the figure indicate that the fine grid can capture more small-scale vortices than the medium and coarse grids with the same Q value. In the near-wake region, the number of small-scale vortices captured by the coarse grid is significantly smaller than that captured by the fine grid, whereas the difference in the flow structures between the medium and fine grid was found to be rather small. Combing with the mesh independence study on the pressure, the resultant velocity and the Q criterion, the medium grid was found to have sufficient resolution for the purpose of the present study. where Ω is the rate-of-rotation tensor, and S is the rate-strain-tensor. Figure 9 shows the instantaneous iso-surfaces of Q = 20,000 s −2 of three sets of meshes at the t = 0.08 s position, coloured according to the resultant velocity U. The V2 − V2′ and V1 − V1′ regions in the figure indicate that the fine grid can capture more small-scale vortices than the medium and coarse grids with the same Q value. In the near-wake region, the number of small-scale vortices captured by the coarse grid is significantly smaller than that captured by the fine grid, whereas the difference in the flow structures between the medium and fine grid was found to be rather small. Combing with the mesh independence study on the pressure, the resultant velocity and the Q criterion, the medium grid was found to have sufficient resolution for the purpose of the present study.  Figure 3).   Figure 3). where Ω is the rate-of-rotation tensor, and S is the rate-strain-tensor. Figure 9 shows the instantaneous iso-surfaces of Q = 20,000 s −2 of three sets of meshes at the t = 0.08 s position, coloured according to the resultant velocity U. The V2 − V2′ and V1 − V1′ regions in the figure indicate that the fine grid can capture more small-scale vortices than the medium and coarse grids with the same Q value. In the near-wake region, the number of small-scale vortices captured by the coarse grid is significantly smaller than that captured by the fine grid, whereas the difference in the flow structures between the medium and fine grid was found to be rather small. Combing with the mesh independence study on the pressure, the resultant velocity and the Q criterion, the medium grid was found to have sufficient resolution for the purpose of the present study.  Figure 3).

Experimental Validation and Transient Pressure
A moving model test was used to obtain the experimental data. The real moving model test scene and the comparison between the test results and numerical simulation results are presented in Figure 10. This test was performed at the Central South University in China. Zhang et al. [41] had provided detailed information about the experiment. Two 1/16th-scale MTs with two-car marshalling for each train were adopted, and the trains met at an equal speed of 200 km/h for each train. The transient pressure on the train surface was measured using pressure sensors with a sampling frequency of 10 4 Hz.
was measured using pressure sensors with a sampling frequency of 10 4 Hz.
The pressure in the monitoring points on the train surface was computed using different turbulence models (SST, SA−IDDES, and SST−IDDES), and the results were compared to the moving model test results for validating the numerical simulations. Figure  10c, d show the time-history curves of pressure at monitoring point No. 5 on the head car surface and monitoring point No. 17 on the tail car surface. As shown, the change trend of the time-history pressure curves obtained by the three models is consistent with the test results, and only the local region of the curves exhibits a small error due to the wake flow passing through it. The values of the negative peak pressure coefficient ΔCPmin, the positive peak pressure coefficient ΔCPmax, and peak-to-peak pressure coefficient ΔCP simulated by the three turbulence models are slightly different from the experimental results. Among them, the differences in ΔCP at monitoring point No. 5 between SST, SA−IDDES, and SST−IDDES and the experimental data were 1.2%, 1.7%, and 0.7%, respectively, and the differences in ΔCP at monitoring point No. 17 between SST, SA−IDDES, and SST−ID-DES and the experimental data were 2.9%, 3.2%, and 3.6%, respectively. Thus, the simulated results are in good agreement with the experimental data, and the numerical simulation can be considered sufficiently reliable. The results of the simulation suggest that the turbulence model had minimal effect on the transient pressure. Therefore, the pressure change amplitude on the intersection side of the train surface was obtained only with the SST−IDDES model; the results are shown in Figure 11. For monitoring points along the longitudinal direction of the train surface, the absolute values of ΔCPmin, ΔCPmax, and ΔCP were higher in the streamlined head and tail nose regions than in other regions of the train surface, and the difference in the The pressure in the monitoring points on the train surface was computed using different turbulence models (SST, SA−IDDES, and SST−IDDES), and the results were compared to the moving model test results for validating the numerical simulations. Figure 10c, d show the time-history curves of pressure at monitoring point No. 5 on the head car surface and monitoring point No. 17 on the tail car surface. As shown, the change trend of the time-history pressure curves obtained by the three models is consistent with the test results, and only the local region of the curves exhibits a small error due to the wake flow passing through it. The values of the negative peak pressure coefficient ∆C Pmin , the positive peak pressure coefficient ∆C Pmax , and peak-to-peak pressure coefficient ∆C P simulated by the three turbulence models are slightly different from the experimental results. Among them, the differences in ∆C P at monitoring point No. 5 between SST, SA−IDDES, and SST−IDDES and the experimental data were 1.2%, 1.7%, and 0.7%, respectively, and the differences in ∆C P at monitoring point No. 17 between SST, SA−IDDES, and SST−IDDES and the experimental data were 2.9%, 3.2%, and 3.6%, respectively. Thus, the simulated results are in good agreement with the experimental data, and the numerical simulation can be considered sufficiently reliable.
The results of the simulation suggest that the turbulence model had minimal effect on the transient pressure. Therefore, the pressure change amplitude on the intersection side of the train surface was obtained only with the SST−IDDES model; the results are shown in Figure 11. For monitoring points along the longitudinal direction of the train surface, the absolute values of ∆C Pmin , ∆C Pmax , and ∆C P were higher in the streamlined head and tail nose regions than in other regions of the train surface, and the difference in the pressure variation amplitude in the middle region of the train surface was found to be small. The maximum value of ∆C P for these monitoring points was measured at point No. 1. For monitoring points along the vertical distribution in the middle of the train surface, the absolute values of ∆C Pmin and ∆C P first increased and then decreased, and ∆C Pmax increased with increasing height of the monitoring points. At these four measurement points, the maximum value of ∆C P was measured at point No. 10.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 11 of 19 pressure variation amplitude in the middle region of the train surface was found to be small. The maximum value of ΔCP for these monitoring points was measured at point No. 1. For monitoring points along the vertical distribution in the middle of the train surface, the absolute values of ΔCPmin and ΔCP first increased and then decreased, and ΔCPmax increased with increasing height of the monitoring points. At these four measurement points, the maximum value of ΔCP was measured at point No. 10.  Figure 12 shows time-history curves of U with different turbulence models at y = 0.375 m and z = 0.188 m (Figure 3). In the time-history curves of U simulated by different turbulence models, a peak appears when the head and tail cars pass by, and the first peak value was found to be larger than the second peak value, in agreement with the observations of Huang et al. [1]. However, the value of U is different when wake flows pass by the monitoring point. Wang et al. [42] reported similar observations.  Figure 3). Figure 12 also depicts that the main difference in the slipstream velocity obtained by different turbulence models is in the wake region. Hence, the velocity fields in the wake for different time periods are analysed. To facilitate the analysis, the three positions after the two MTs passed each other are defined and denoted as t1, t2, and t3, as shown in Figure  13. Figure 14 presents the flow fields of resultant velocity U at positions t2 and t3. As shown, the vortex structures of the intersection side interact, leading to a decrease in their intensity. Among them, the SA−IDDES and SST−IDDES methods can capture extra small-  Figure 12 shows time-history curves of U with different turbulence models at y = 0.375 m and z = 0.188 m (Figure 3). In the time-history curves of U simulated by different turbulence models, a peak appears when the head and tail cars pass by, and the first peak value was found to be larger than the second peak value, in agreement with the observations of Huang et al. [1]. However, the value of U is different when wake flows pass by the monitoring point. Wang et al. [42] reported similar observations.   Figure 12 shows time-history curves of U with different turbulence models at y = 0.375 m and z = 0.188 m (Figure 3). In the time-history curves of U simulated by different turbulence models, a peak appears when the head and tail cars pass by, and the first peak value was found to be larger than the second peak value, in agreement with the observations of Huang et al. [1]. However, the value of U is different when wake flows pass by the monitoring point. Wang et al. [42] reported similar observations.  Figure 3). Figure 12 also depicts that the main difference in the slipstream velocity obtained by different turbulence models is in the wake region. Hence, the velocity fields in the wake for different time periods are analysed. To facilitate the analysis, the three positions after the two MTs passed each other are defined and denoted as t1, t2, and t3, as shown in Figure  13. Figure 14 presents the flow fields of resultant velocity U at positions t2 and t3. As shown, the vortex structures of the intersection side interact, leading to a decrease in their intensity. Among them, the SA−IDDES and SST−IDDES methods can capture extra small-  Figure 3). Figure 12 also depicts that the main difference in the slipstream velocity obtained by different turbulence models is in the wake region. Hence, the velocity fields in the wake for different time periods are analysed. To facilitate the analysis, the three positions after the two MTs passed each other are defined and denoted as t 1 , t 2 , and t 3 , as shown in Figure 13. Figure 14 presents the flow fields of resultant velocity U at positions t 2 and t 3 . As shown, the vortex structures of the intersection side interact, leading to a decrease in their intensity. Among them, the SA−IDDES and SST−IDDES methods can capture extra small-scale vortices, whereas the SST k-ω cannot capture these small-scale vortices. For quantitative analysis, two horizontal sampling lines were extracted in the train wake at positions t 2 and t 3 , as shown in Figure 14. Figure 15 shows the U values at the four lines ( Figure 14). In contrast to the case of a single train, the U value on the intersection side is significantly smaller than that on the non-intersection side when the two MTs passing each other. Furthermore, the U value calculated with the SST k-ω method is significantly lower than that calculated with the SA−IDDES and SST−IDDES models. scale vortices, whereas the SST k-ω cannot capture these small-scale vortices. For quantitative analysis, two horizontal sampling lines were extracted in the train wake at positions t2 and t3, as shown in Figure 14. Figure 15 shows the U values at the four lines ( Figure 14). In contrast to the case of a single train, the U value on the intersection side is significantly smaller than that on the non-intersection side when the two MTs passing each other. Furthermore, the U value calculated with the SST k-ω method is significantly lower than that calculated with the SA−IDDES and SST−IDDES models.    scale vortices, whereas the SST k-ω cannot capture these small-scale vortices. For quantitative analysis, two horizontal sampling lines were extracted in the train wake at positions t2 and t3, as shown in Figure 14. Figure 15 shows the U values at the four lines ( Figure 14). In contrast to the case of a single train, the U value on the intersection side is significantly smaller than that on the non-intersection side when the two MTs passing each other. Furthermore, the U value calculated with the SST k-ω method is significantly lower than that calculated with the SA−IDDES and SST−IDDES models.   scale vortices, whereas the SST k-ω cannot capture these small-scale vortices. For quantitative analysis, two horizontal sampling lines were extracted in the train wake at positions t2 and t3, as shown in Figure 14. Figure 15 shows the U values at the four lines ( Figure 14). In contrast to the case of a single train, the U value on the intersection side is significantly smaller than that on the non-intersection side when the two MTs passing each other. Furthermore, the U value calculated with the SST k-ω method is significantly lower than that calculated with the SA−IDDES and SST−IDDES models.    Figure 16 shows the instantaneous iso-surface of the second invariant of the velocity gradient Q = 20,000 s −2 simulated with different turbulence models at positions t 1 and t 3 coloured using the resultant velocity U. Figure 16a presents the iso-surface of Q at position t 1 . A pair of twin counter-rotating longitudinal vortices (V 1 and V 2 or V 1 and V 2 ) were observed for all three cases. However, the iso-surface of Q calculated by different turbulence models was found to be significantly different. Among them, the SST k-ω model could not capture small-scale vortices. However, a large number of small-scale vortices were observed in the SA−IDDES and SST−IDDES models, and more small-scale vortices were captured with the SST−IDDES than with the SA−IDDES, particularly in the near-wake region. Figure 16b shows the iso-surface of Q at position t 3 . The V 1 −V 1 region was induced by the interaction of two wake vortices on the intersection side. With the SST k-ω model, the interaction between the two wake vortices in this region was strong and resulted in an 'entanglement' phenomenon. A similar phenomenon was also observed in the simulations using the SA−IDDES and SST−IDDES models. Furthermore, the flows in this region exhibited strong turbulence. For vortices V 2 and V 2 , the SA−IDDES and SST−IDDES models resulted in the formation and development of the vortex structures, and their intensity gradually decreased until they dissipated with the increasing distance from the tail car. However, because the RANS model dissipated the turbulent eddy viscosity, it could not resolve the small vortices Li, et al. [43]. Therefore, we did not observe small-scale vortices at distances from the tail car larger than 7.58 H in the SST k-ω model.

Wake Flows
Next, the streamlines projected on the vertical plane were used to study the local interaction of wake vortices as two MTs pass each other. Figure 17 shows the streamlines computed by different turbulence models in the plane X/H = 0 at t 2 and t 3 . The vertical planes were coloured by the dimensionless longitudinal velocity (U x = u/u tr ). The wake vortex of a single train travelling without crosswind is symmetrically distributed [17,18,22]. However, Figure 17 shows that the interaction between vortices V c1 and V c1 on the intersection side causes vortices V c1 and V c2 (or V c1 and V c2 ) to be asymmetric. This phenomenon is observed with three turbulence models. The two vortices V c1 and V c1 on the intersection side interact, causing their intensity to decrease and their width to be much smaller than that of the vortices on the non-intersection side. To provide a more quantitative comparison of the size of longitudinal vortices, the height of the vortex on the intersection side was denoted as h 0 , and the height and distance of the vortex core on the intersection side were denoted as h 1 and w, respectively, which is similar to the definition of Wang et al. [29], as shown in Figure 17b. Figure 17a shows the streamlines at position t 2 . It is evident that the turbulence model has minimal effect on wake vortices V c2 and V c2 on the non-intersection side. However, vortices V c1 and V c1 on the intersection side obtained by different turbulence models are significantly different. Figure 18 compares the h 0 , h, and w values for the different models at different moments. In Figure 18a, it can be observed that the differences in vortex core height h computed by different models are small, while height h 0 of vortex V c1 and distance w between the two vortex cores are significantly different. The height h obtained with the SST−IDDES model is 10% greater than that obtained with the SA−IDDES model, and the shortest and longest distances w of the two vortex cores are obtained with the SST−IDDES and SA−IDDES, respectively, with a difference of 10.3%. be rather similar. Nevertheless, the h0, h, and w values are significantly different, as trated in Figure 18b. The highest and lowest values of h0 were obtained with SST k-SA-IDDES, respectively, with a difference of 11.8%, and the highest and lowest val h were obtained with SST−IDDES and SA−IDDES, respectively, with a difference o Furthermore, the highest and lowest values of w were obtained with SST−IDDES an k-ω, respectively, with a difference of 30.9%.

Aerodynamic Forces
The drag-, side-, and lift-force coefficients were calculated as follows: where S is the cross-section area in the middle of the train; Cd, Cs, and Cl are the drag, side, and lift force coefficients, respectively; ρ is air density; and utr is train speed. Fd, Fs, and Fl are the drag, side, and lift forces, respectively. A MT running in the positive direction of x coordinates was used to analyse the aerodynamic force.
The results of the drag, side, and lift force coefficients obtained with different turbulence models are shown in Figure 19. Figure 19a,b show the time-history of the drag force coefficients of the head and tail cars, respectively. It can be observed from the figure that the maximum absolute value of the drag force coefficient of the head car is around the position where two head noses meet, and the maximum absolute value of the drag force coefficient of the tail car is around the position where two tail noses meet. The drag force coefficient obtained with the SA−IDDES turbulence model is remarkably different from that obtained with the SST k-ω and SST-IDDES turbulence models. For the head and tail cars, the SA−IDDES and SST−IDDES turbulence models exhibit the maximum differences

Aerodynamic Forces
The drag-, side-, and lift-force coefficients were calculated as follows: where S is the cross-section area in the middle of the train; Cd, Cs, and Cl are the drag, side, and lift force coefficients, respectively; ρ is air density; and utr is train speed. Fd, Fs, and Fl are the drag, side, and lift forces, respectively. A MT running in the positive direction of x coordinates was used to analyse the aerodynamic force.
The results of the drag, side, and lift force coefficients obtained with different turbulence models are shown in Figure 19. Figure 19a,b show the time-history of the drag force coefficients of the head and tail cars, respectively. It can be observed from the figure that the maximum absolute value of the drag force coefficient of the head car is around the position where two head noses meet, and the maximum absolute value of the drag force coefficient of the tail car is around the position where two tail noses meet. The drag force coefficient obtained with the SA−IDDES turbulence model is remarkably different from that obtained with the SST k-ω and SST-IDDES turbulence models. For the head and tail cars, the SA−IDDES and SST−IDDES turbulence models exhibit the maximum differences  Figure 17b presents the streamlines at position t 3 . As the distance from the tail car increases, the vortex size increases and the interaction between vortices V c1 and V c1 becomes stronger, leading to a greater change in their shape. The height and width of the vortices are found to be significantly larger at position t 3 compared to position t 2 . Vortices V c2 and V c2 on the non-intersection side with the different turbulence models are found to be rather similar. Nevertheless, the h 0 , h, and w values are significantly different, as illustrated in Figure 18b. The highest and lowest values of h 0 were obtained with SST k-ω and SA-IDDES, respectively, with a difference of 11.8%, and the highest and lowest values of h were obtained with SST−IDDES and SA−IDDES, respectively, with a difference of 17%. Furthermore, the highest and lowest values of w were obtained with SST−IDDES and SST k-ω, respectively, with a difference of 30.9%.

Aerodynamic Forces
The drag-, side-, and lift-force coefficients were calculated as follows: where S is the cross-section area in the middle of the train; C d , C s , and C l are the drag, side, and lift force coefficients, respectively; ρ is air density; and u tr is train speed. F d , F s , and F l are the drag, side, and lift forces, respectively. A MT running in the positive direction of x coordinates was used to analyse the aerodynamic force.
The results of the drag, side, and lift force coefficients obtained with different turbulence models are shown in Figure 19. Figure 19a,b show the time-history of the drag force coefficients of the head and tail cars, respectively. It can be observed from the figure that the maximum absolute value of the drag force coefficient of the head car is around the position where two head noses meet, and the maximum absolute value of the drag force coefficient of the tail car is around the position where two tail noses meet. The drag force coefficient obtained with the SA−IDDES turbulence model is remarkably different from that obtained with the SST k-ω and SST−IDDES turbulence models. For the head and tail cars, the SA−IDDES and SST−IDDES turbulence models exhibit the maximum differences in the peak values of the drag force coefficient (as depicted in the green ellipse), which are 9.3% and 4.5%, respectively. peak amplitude of approximately −0.64 is generated near a position where the two he noses meet. A positive peak amplitude of approximately 0.69 is generated near the po tion where the head and tail noses meet. For the tail car, a negative peak amplitude approximately −0.42 is generated near a position where the two tail noses meet, and positive peak amplitude of approximately 0.7 is generated near a position where the he and tail noses meet. Figure 19e,f show the time histories of the lift force coefficients of the head and t cars, respectively. As can be observed from Figure 19e, the turbulence model has minim effect on the lift force of the head car before the head nose of a train arrives at the midd of the other train, while it has a strong influence after that position. The positive and ne ative peak amplitudes obtained with the different turbulence models were approximate 0.14 and −0.14, respectively. As shown in Figure 19f, the lift force of the tail car is sign cantly different for the different turbulence models. For the maximum value of the force (green ellipse in Figure 19f), the value computed with the SST−IDDES model is 11 larger than that computed with the SA−IDDES model. The lift force of the tail car is close related to the wake vortex. There are large differences in the lift force of the tail car amo the different turbulence models because different turbulence models have different ab ties to capture wake flow. In addition, a minimum value appears when two tail nos meet. Overall, the variation trends of the drag, side, and lift forces computed by differe models are the same, but their amplitudes are different.  Figure 19c,d show the time histories of the side force coefficients of the head and tail car, respectively. The differences in the variation trend and the amplitude of the side force of the head and tail cars obtained by different turbulence models are considerably small, and the differences in the peak values are all lower than 6%. For the head car, a negative peak amplitude of approximately −0.64 is generated near a position where the two head noses meet. A positive peak amplitude of approximately 0.69 is generated near the position where the head and tail noses meet. For the tail car, a negative peak amplitude of approximately −0.42 is generated near a position where the two tail noses meet, and a positive peak amplitude of approximately 0.7 is generated near a position where the head and tail noses meet. Figure 19e,f show the time histories of the lift force coefficients of the head and tail cars, respectively. As can be observed from Figure 19e, the turbulence model has minimal effect on the lift force of the head car before the head nose of a train arrives at the middle of the other train, while it has a strong influence after that position. The positive and negative peak amplitudes obtained with the different turbulence models were approximately 0.14 and −0.14, respectively. As shown in Figure 19f, the lift force of the tail car is significantly different for the different turbulence models. For the maximum value of the lift force (green ellipse in Figure 19f), the value computed with the SST−IDDES model is 11% larger than that computed with the SA−IDDES model. The lift force of the tail car is closely related to the wake vortex. There are large differences in the lift force of the tail car among the different turbulence models because different turbulence models have different abilities to capture wake flow. In addition, a minimum value appears when two tail noses meet. Overall, the variation trends of the drag, side, and lift forces computed by different models are the same, but their amplitudes are different.

Conclusions
In this study, three turbulence models (SST k-ω, SA−IDDES, and SST−IDDES) were used to simulate the aerodynamic performance of two MTs passing by each other. The accuracy of the numerical simulation method was verified against the data obtained from the moving model test. The results of the simulations and experimental data showed good agreement. The main conclusions are as follows.
Only small differences were found in the transient pressure simulated by different turbulence models. Among them, the differences in the peak-to-peak values of the pressure were found to be less than 5%. When the head and tail cars passed each other, a peak value in the slipstream velocity appeared, and the first peak value was larger than the second. The choice of the turbulence model was found to have a significant effect on the slipstream. The resultant velocity values in the wake calculated with the SST k-ω model were significantly lower than those calculated with the SA−IDDES and SST−IDDES models.
The complex wake flows computed using the different turbulence models were significantly different. Among them, the SA−IDDES and SST−IDDES models could capture two main vortices, as well as small-scale vortices. Large amounts of small-scale vortices were obtained with the SST−IDDES model than with the SA−IDDES model, specifically in the near-wake region. The SST k-ω model failed to capture small-scale vortices. The streamlines projected on the plane were used to analyse the local region of the wake vortex. For the streamlines projected on the plane at the x = 0 position, the vortices on the intersection side of the trains were found to interact, resulting in their weakening. Different turbulence models resulted in a larger difference in the vortex structure on the intersection side but a lower difference in the vortex structure on the non-intersection side.
The drag force simulated with the SA−IDDES model was significantly different from that simulated with the SST k-ω and SST−IDDES models for both head and tail cars, particularly the peak value. For the head car, the largest discrepancy of the maximum absolute value of the drag force was between the SST-IDDES model and SA−IDDES model, with a difference of 9.3%. For the tail car, the difference between the maximum absolute peak value (SA−IDDES) and minimum absolute peak value (SST−IDDES) of the drag force was 4.5%. The turbulence model had little effect on the side force, where the variation trend is consistent and the difference in amplitude is small (less than 6%). The lift force was more sensitive to the choice of turbulence model, particularly that of the tail car. The largest discrepancy in the maximum lift force of the tail car was observed between the SST−IDDES and SA−IDDES models, with a difference of 11%.
The RANS model is more suitable in terms of the transient pressure, drag and side forces because their values predicted by different models have a small difference, while its computation time is merely about 0.64 times that of the IDDES model based on the k−ω. If we want to observe changes in the flow field around trains when two MTs passing each other, such as the fluid flow between the train and the track as well as the wake flow, the IDDES model based on k-ω is more recommendable because it can capture more small vortices. In future work, more turbulence models will be considered to simulate aerodynamic performance when two MTs passing each other, such as large eddy simulation (LES), embedded large eddy simulation (ELES) and partially−averaged Navier-Stokes (PANS).