This paper developed a detailed fluid dynamics model and a parallel computing scheme for air brake systems on long freight trains. The model consists of subsystem models for pipes, locomotive brake valves, and wagon brake valves. A new efficient hose connection boundary condition that considers pressure loss across the connection was developed. Simulations with 150 sets of wagon brake systems were conducted and validated against experimental data; the simulated results and measured results reached an agreement with the maximum difference of 15%; all important air brake system features were well simulated. Computing time was compared for simulations with and without parallel computing. The computing time for the conventional sequential computing scheme was about 6.7 times slower than real-time. Parallel computing using four computing cores decreased the computing time by 70%. Real-time simulations were achieved by parallel computing using eight computer cores.
Brake is a critical aspect of railway train operations. Railway brake systems  are mainly of the air brake type which relies on pressurized air to push brake shoes or pads against wheel treads or brake disks. This paper focuses on conventional automatic (fail-safe) air brake systems on long freight trains; electronically controlled pneumatic (ECP) brake systems  and other alternative brake systems  are not discussed. Comparatively, the automatic air brake systems on long freight trains pose greater challenges for modeling and simulations in the sense of complicated fluid dynamics modeling and intensive computing loads. The modeling approach presented in this paper can also be used to model passenger train air brake systems; it can also be modified to simulate ECP brake systems.
Numerical modeling and simulations are indispensable approaches to study and improve air brake systems. A review regarding air brake modeling prior to the 1990s can be found in Ref. . Another more recent review  regarding longitudinal train dynamics has classified air brake models into empirical models, fluid dynamics models, and empirical–fluid dynamics models. Empirical models [5–8] use mathematical equations to fit measured characteristics of air brake systems. These empirical models do not necessarily follow fluid dynamics principles behind the measured characteristics. In comparison, fluid dynamics models [9–15] describe air flows according to fluid dynamics principles. The flows described include flows in brake pipes and those between various components of brake valves. Empirical–fluid dynamics models [16,17] blend the fluid dynamics models and empirical models by retaining fluid dynamics pipe models but modeling brake valves using empirical functions.
Fluid dynamics air brake models are clearly the most advanced among all three types of models. However, they do have a well-known but not widely discussed disadvantage: low computing speed. The authors' experience with both empirical and fluid dynamics models indicates that, with the same step-size, the computing speed of an empirical model can be 100 times faster than a fluid dynamics model. It is difficult to achieve real-time simulations using fluid dynamics models for long heavy haul trains (e.g., 241 vehicles) due to the large number of equations to solve. The computational efficiency issue of fluid dynamics models is also an important reason why empirical models are still commonly used [6–8]. Empirical–fluid dynamics models [16,17] can evidently improve the computing speed but they also inherit some limitations of the empirical models. For example, the empirical valve models only cover the specific measured valve types and working conditions used to develop such models.
This paper develops a detailed fluid dynamics air brake model for long heavy haul trains. The model includes subsystem models for pipes, locomotive brake valves, and wagon brake valves. Experimental data are presented for validation. The novelty of this study is also contributed by the parallel computing scheme that significantly improves the computing speed of the fluid dynamics model. Computing time is compared for simulations with and without parallel computing. A new hose connection boundary condition is developed to consider air pressure loss while crossing the hose connection; a pipe tee boundary condition is improved to also consider air pressure loss.
Automatic Air Brake
The automatic air brake studied in this paper is illustrated as Fig. 1. It consists of one (or more for distributed power trains) locomotive brake system and a number of wagon brake systems (one system per wagon). All systems are connected by hoses which allow intervehicle movements. Brake application and release signals are transmitted via the pressurized air in the brake pipe.
Locomotive Brake System.
The locomotive brake system is shown as Fig. 2. Note that this figure shows only the train automatic brake parts; locomotive independent brakes that are specifically used for locomotives are not discussed in this paper. In this system, the air compressor detects the pressure of the main reservoir and automatically starts or stops the feed of pressurized air to the main reservoir as required to maintain pressure. The regulating valve controls the openings of and according to the position of the driver's handle and the pressure in the equalizing reservoir. The equalizing reservoir has a small volume (15 ) and acts as a pilot for the pressure change of the brake pipe. The relay valve detects the pressure difference between the equalizing reservoir and the brake pipe so as to change the openings of and . When the equalizing reservoir pressure is higher than the pipe pressure, is opened and the main reservoir charges the brake pipe until an equilibrium is reached between the equalizing reservoir and the brake pipe. When the pressure difference is otherwise, is opened and the brake pipe discharges to atmosphere until an equilibrium is reached.
Wagon Brake System.
The wagon brake system utilizes the well-known principle of a “triple valve” : brakes are applied when brake pipe pressure is lower than auxiliary reservoir pressure, while brakes are released when the pressure difference is otherwise. These basic actions are achieved via the sliding valve and graduating valve as shown in Fig. 3. The service piston divides the volume of a distributor into an upper chamber and a lower chamber which are, respectively, connected to brake pipe and auxiliary reservoir. The service piston changes the positions of the graduating valve and the sliding valve according to different pressure situations in the upper and lower chambers. There are various ports, orifices, and connecting grooves in the graduating valve, sliding valve, and sliding valve seat. Therefore, different positions of valves can form different air passages.
The studied wagon brake system (see Fig. 4) has the following features: quick service (two stages), braking, lapping, quick release, release, regulated recharge, and emergency braking (two stages). Note that in Fig. 4, , , and are uncontrolled openings (orifices or chokes); all other openings are controlled by various valves.
Quick service (stage 1 or preliminary): The graduating valve and the sliding valve uncover port . The quick service bulb has a small volume (0.6 L) and is opened to the atmosphere via a choke . Therefore, a small amount of pressurized air can be vented from the brake pipe to atmosphere via the upper chamber of the distributor. The quick service feature quickly decreases the brake pipe pressure so as to accelerate the braking process in the wagon brake system; it also accelerates the brake propagation rate in the brake pipe. Note that the pipe pressure reduction during this stage is small (about 10 kPa).
Quick service (stage 2): The graduating valve and the sliding valve close port . The second stage of quick service is achieved via the quick service limiting valve (or in-shot valve) which detects the pressure of the brake cylinder. When the cylinder pressure is lower than a certain value (60 kPa), the valve keeps port open and allows pressurized air flows into the brake cylinder from the brake pipe via the upper chamber. The second stage of quick service has the same functions as the first stage; in addition, it helps increase the cylinder pressure (in-shot) to help with extending the push rod.
Braking: The braking process occurs simultaneously with the second stage of quick service. During braking, the graduating valve and the sliding valve uncover which allows pressurized air flows from the auxiliary reservoir to the brake cylinder via the lower chamber.
Lapping: When the pressure difference between the upper chamber and the lower chamber is not enough to hold the graduating valve in its up position, the graduating valve will drop and close port which stops the pressure increase in the brake cylinder. Note that the sliding valve is not moving during lapping; therefore, brake release is not activated.
Quick release: When brake release is activated, the pressurized air in the brake cylinder will be exhausted to the atmosphere. When the quick release valve has detected a certain rate of out flow from the brake cylinder, the quick release valve can open the check valve which allows pressurized air that is stored in the quick release reservoir to be charged to the brake pipe via the upper chamber. The quick release feature accelerates the release process in the wagon brake system; it also increases the propagation rate of brake release in the brake pipe.
Release: The graduating valve and the sliding valve close port and uncover port . The pressurized air in the brake cylinder is exhausted to the atmosphere. Note that the release cannot be lapped, which means graduated release is not an option for this wagon brake system. During the brake release, ports , , and are also opened; therefore, the auxiliary reservoir and the quick release reservoir can be charged to operating pressure (500 kPa or 600 kPa) and be prepared for future actions. Port can also be opened if there was any pressure loss in the emergency bulb during the brake application.
Regulated recharge: During a release process, when the upper chamber pressure is significantly higher than that of the lower chamber, the cage will be pushed by the sliding valve and the return spring will be deflected. The sliding valve will then decrease the opening size of port so as to decrease the recharging rate of the auxiliary reservoir. This feature slows the recharging of auxiliary reservoirs toward the front of the train during brake release processes, which can help to achieve a uniform recharging rate and brake release along the train.
Emergency braking (stage 1): When the emergency valve detects significant pressure difference between the emergency bulb and upper chamber, the emergency valve can open the check valve which exhausts the pressurized air in the brake pipe to the atmosphere via the upper chamber. During the emergency braking, auxiliary reservoir, lower chamber, and brake cylinder have the same air passages as those during the braking process that has been described previously. The first stage of emergency braking is very quick and drastic; it can increase brake cylinder pressure from 0 to 160 kPa in about 2 s. The fast increase of pressure is limited to about 160 kPa then a slower response in stage 2 initiates.
Emergency braking (stage 2): As the first stage of emergency braking is quick and drastic, it is not desirable from the train dynamics point of view to keep applying brakes at this rate. For this reason, the wagon brake system has a dual-stage emergency braking feature which decreases the opening size of (indirectly) when the cylinder pressure has reached 160 kPa. This feature is designed and tuned to make sure that the train has sufficient braking capability and acceptable train dynamics performance during emergency braking.
Pipes that were modeled (with the consideration of air flows) for each wagon include two sections of brake pipe and one section of branch pipe (see Fig. 4). The modeling implies conservation laws and various boundary conditions. Hose connections were modeled as boundary conditions. Pipes that connect the quick release reservoir, auxiliary reservoir, and brake cylinder were simplified as additional volumes .
The pipe model is based on a gas dynamics model developed by Benson et al.  for the purpose of internal combustion engine simulations. The model is illustrated as Fig. 5 while the governing equations are expressed as 
where is the reference pressure. Numerical steps to determine the Riemann variables are given in the Appendix.
As can be seen from Figs. 2 and 4, four types of boundary conditions are needed for the pipe modeling: (1) closed end boundary which can be used to simulate the end of the brake pipe; (2) partially open boundary which can be used to simulate the interface between the branch pipe and upper chamber or between the brake pipe (choked) and atmosphere; (3) hose connection; and (4) pipe tee. The closed end and partially open boundaries developed in Ref.  are used for this air brake model; they are not discussed in this paper. This section develops a new hose connection boundary with the consideration of pressure loss and improves the pipe tee boundary condition to also consider pressure loss.
As reviewed previously, the resistance caused by hose connections can be emulated by using additional friction  or decreasing the pipe cross section  at the positions of hose connections. Both methods are effective but are not suitable for this study as they treat the brake pipe of the whole train as one single piece. To be able to formulate the model to conduct parallel computing, explicit boundaries are needed to separate individual vehicles. Therefore, hose connections are better to be modeled as boundary conditions for this study. Benson et al.  have developed a boundary condition for devices with adiabatic pressure loss which can be used to simulate hose connections. However, this boundary condition was developed for nonhomentropic flow and an iterative process must be used to numerically solve the boundary condition. For fluid dynamics air brake models that have low computing speed, a more efficient boundary condition is desirable. Under this situation, this study develops a new hose connection boundary condition that does not need an iterative solving process.
As shown in Fig. 6, if pressure loss is not considered while the air flow crosses the connection, Riemann variables will have the same values at the two connected ends, i.e., and , where is the Riemann variable that enters the boundary while is the Riemann variable that is reflected back from the boundary. As the air brake model has assumed an isothermal process and that there is no entropy change in the boundary condition, both and remain unchanged during the process. Remembering that the air pressure is determined by the sum of and (Eq. (8)), the pressure changes can therefore be simulated by changing the values of . It is understandable that the resistance of the hose connection will result in pressure build up upstream of the connection (reflecting pressure waves) and pressure loss downstream as shown in Fig. 6. In other words, with the consideration of pressure loss, the upstream pressure will be higher than the case without pressure loss while the downstream pressure will be lower than the case without pressure loss. To quantify the pressure loss, the definition of “pressure loss coefficient” that has been used by Benson et al.  is referred. Benson et al.  related the pressure loss to the density and velocity of the flow. In this study, the following expression was developed:
where is the variation for and is the pressure loss coefficient. As the pressure loss is small for each hose connection and the volumes of the upstream pipe and downstream pipe are the same, it was assumed that the pressure increase upstream and the pressure loss downstream are of the same value as Eq. (9).
where and indicate the sequence of pipe section; is the cross-sectional area of the pipe section; is the sum of all three cross-sectional areas. By integrating Eq. (9) into Eq. (10), a pipe tee boundary condition that considers pressure loss can be developed. The same as for the pipe connection boundary condition, in the pipe tee boundary condition, remain the same during the process due to the homentropic assumption while is modified according to Eq. (9) and its corresponding flow direction.
Brake Valve Models
Modeling of brake valves implies modeling of moving parts (valves and pistons) and air flows between various volumes. This section selects the two most complicated modeling cases to exemplify the method of brake valve modeling: modeling the service portion of the wagon brake system including the service piston, graduating valve, and sliding valve (Fig. 3) and modeling air flows between the auxiliary reservoir and brake cylinder.
where is the physical limitation for the position of the piston, is the step-size, and is the limit for displacement change rate.
where is the physical limit for the lowest position of the sliding valve and is the limit for displacement change rate. Having determined the positions of the sliding valve and the graduating valve, the opening size of the ports can be determined according to the design of the air brake system.
The case of an auxiliary reservoir charging a brake cylinder (Fig. 7) is used to exemplify modeling of air flows between various volumes in the brake system. Both cylinder and auxiliary reservoir are connected to the distributor via pipes. In this model, the reservoir pipe and the cylinder pipe are treated as additional volumes and added into the volumes of the auxiliary reservoir and brake cylinder, respectively. The distributor controls the size of the opening that connects the reservoir and the cylinder. Air flows between the two volumes can be modeled as 
where is the additional volume converted from the brake cylinder pipe.
The air brake model is validated against measured data from the test facility in CRRC Brake Science & Technology Co., Ltd., Meishan, China. The test facility, shown in Fig. 8, has the capability of testing up to 150 sets of connected wagon air brake systems at a time and is able to conduct tests at operating pressures of both 500 kPa and 600 kPa. All presented pressures in this paper indicate gauge pressures while absolute pressures have to be used in the modeling process. Brake pipe pressures and brake cylinder pressures under different brake situations were measured for both 500 kPa and 600 kPa operating pressures. This section presents the measured and simulated results for the 600 kPa operating pressure case which is the norm in Chinese heavy haul railways. Key parameters of the brake system are listed in Table 1 and the results are presented in Fig. 9. It can be seen that the model has considered 150 sets of wagon air brake systems and the total brake pipe length is over 2 km.
Figure 9 presents the brake pipe pressures and brake cylinder pressures at three wagon positions: the first wagon (001), the last wagon (150), and the one in the middle (075). Pressures that have the initial value of 600 kPa are brake pipe pressures while those that have the initial value of 0 kPa are brake cylinder pressures. Three typical brake situations are presented: the minimum service braking with 50 kPa pressure reduction in the brake pipe, the full service braking with 170 kPa pressure reduction in the brake pipe, and the emergency braking with a full 600 kPa pressure reduction in the brake pipe. Brake release is only measured and simulated for the minimum service braking case as the minimum braking is the common case used to correct train speeds during train operations. The other two cases of brake applications will surely stop the train and the brake release is of less interest when the train is stopped. During the simulations, the actuating time for brake application and release is determined from the brake pipe pressures of wagon 001 of the measured data. The step-size of 0.1 ms (0.0001 s) was used for numerical integrations. The mesh sizes for brake pipes and branch pipes were 0.96 m and 0.75 m, respectively.
Figures 9(a) and 9(c) have been marked up with a series of air brake system features that have been discussed in Sec. 2. The quick service feature is manifested in both brake pipe pressures and brake cylinder pressures. In brake pipe pressures, the quick service is manifested by the quick pressure decreases at the beginning of brake application; the magnitude of pressure reductions during the quick service is 30–40 kPa. In brake cylinder pressures, the quick service is reflected by the quick increases of cylinder pressures at the beginning of brake application (see Fig. 9(c)); the magnitude of pressure increases is about 60 kPa. The quick service feature cannot be seen in the emergency case as the emergency case already has very quick pressure changes.
The lapping feature is manifested by the stagnation of pressure increases in the brake cylinder during brake application. The quick release feature is manifested by the abrupt pressure increases in the brake pipe at the beginning of brake release. At the end of brake release, the push rods will be retracted. During this process, as brake cylinder volumes are continuously decreasing, brake cylinder pressures decrease at an evidently slower rate as shown in Fig. 9(a). Figure 9(e) clearly shows the dual-stage emergency braking feature. At the first stage of the emergency braking, brake cylinder pressures rapidly increase to about 160 kPa. And then, at the second stage, cylinder pressures gradually increase to about 450 kPa. Comparing the measured results and simulated results, it can be seen that all discussed features have been well simulated.
Table 2 compares some key values for the brake applications between the measured results and the simulated results. Due to the ambiguity of the ends of pressure ascending and descending processes, the ascending time of cylinder pressure for the emergency case is determined at 420 kPa; the pipe pressure descending time for the minimum braking and full braking cases is determined at 555 kPa and 435 kPa, respectively. From an engineering perspective, the simulations have reached a good overall agreement with the measurements. Three relatively evident differences are noted: (1) the simulated maximum brake cylinder pressure in the minimum braking case is about 10% lower than the measured result; (2) the simulated pipe pressure descending time in the full braking case is about 15% shorter than measured results at wagon 075 and wagon 150; and (3) the simulated brake application propagation rate is about 13% faster than the measured results in the minimum and full braking cases.
The following factors could affect the simulation differences: (1) pressure loss modeling in boundary conditions (Eq. (9)), (2) pipe wall friction modeling (a constant coefficient in this paper), (3) temperature (homentropic) modeling, and (4) brake pipe leakage. In this paper, pressure loss and pipe wall friction were modeled empirically as for all other air brake system studies. It is understandable that pressure loss coefficient and pipe wall friction can affect brake application propagation rate and the descending process of pipe pressures. Air temperature variation, which was assumed constant in this paper, can also affect the brake application propagation rate and magnitudes of air pressures. Pipe leakage is not considered in this paper, which can slightly decrease the brake pipe pressures and further increase brake cylinder pressures during brake applications.
Table 3 compares the measured results and simulated results for brake release (Figs. 9(a) and 9(b)). This comparison focuses on three key points of the brake release process. Point 1 is the time when brake release is activated; point 2 is the time when the push rod starts retracting; point 3 is the time when brake cylinder pressure drops to zero. From the time difference between the brake release activations at wagon 001 and wagon 150, the brake release propagation rate can be calculated. The measured brake release propagation rate is 162.52 m/s while the simulated rate is 163.31 m/s. Also from the results in Table 3, it can be seen that the simulated results and the measured results have reached a good agreement.
Parallel Computing Scheme
It is well-known that fluid dynamics air brake models are computationally expensive. For longitudinal train dynamics simulations [20,21], which is a common application area of air brake modeling, fluid dynamics air brake models can be the most computationally intensive part. The long computing time is required by the strong nonlinearity of fluid dynamics and some numerically stiff components in the system. Nonlinear fluid dynamics requires small mesh sizes and small computing step-sizes. In some early applications, one mesh section per vehicle [11,22] and step-sizes larger than 1 s  were used. From today's point of view, these are very rough simulations. The authors' experience indicates that a mesh size of at least 1 m per section is recommended to simulate the details of air brake system behavior. This mesh size is also recommended by Pugi et al. . In terms of numerically stiff components in air brake systems, these are components that have small volumes, such as the upper chamber in the distributor and a low pressure brake cylinder. Figure 10 shows the simulated brake cylinder pressures in the emergency braking case using different step-sizes. During the simulations, the mesh sizes for brake pipes and branch pipes were 0.96 m and 0.75 m, respectively. The 20 ms step-size is about the maximum step-size for a stable solution of the pipe model. It can be seen that the details of push rod extensions are simulated only when the step-size reaches 0.1 ms. With the step-size of 0.1 ms, the simulation time for a 30 s brake application is 202.75 s; the simulation time is about 6.75 times slower than real-time.
Parallel computing is an important and effective way to improve computational efficiency . In the area of railway train dynamics, parallel computing has been successfully used to improve computational efficiency for railway dynamics simulations [20,25,26] and optimizations [27,28]. In this paper, a parallel computing scheme as shown in Fig. 11 was developed in which the wagon brake systems are divided into a number of groups. This figure gives an example of four wagon brake system groups (groups 0–3). The computing of each group is processed by an independent computer core; all groups are simulated in parallel. Different cores are synchronized by time. One core (core 0) has also been used to compute the locomotive brake and boundary conditions. Among the computer cores, core 0 is called the master core while all other cores are called slave cores. In each time-step, each slave core will project two characteristics, i.e., Riemann variables () into the boundaries that are at the two brake pipe ends of the wagon brake group. All entering Riemann variables () are collected by the master core and then determine the corresponding reflecting Riemann variables (). Those reflecting Riemann variables are then sent to corresponding slave cores to proceed with the computing.
In this paper, all simulation codes were written in the fortran language and the Central Queensland University, high performance computing system  was used; the message passing interface  technique was used to enable the parallel computing. The same simulations as shown in Fig. 9 were conducted using the traditional sequential computing scheme (one core computing all brake systems) and two parallel computing schemes: (1) four cores with four wagon groups and (2) eight cores with eight wagon groups. In the first parallel scheme, core 0 was assigned with 30 wagons while all other cores were assigned with 40 wagons each. In the second parallel scheme, core 0 was assigned with ten wagons while all other cores were assigned with 20 wagons each. These arrangements were designed to balance the computing loads among all computing cores as discussed in Ref. . The computing time is listed in Table 4. It can be seen that the computing time of the sequential scheme is about 6.7 times slower than real-time. Using four cores and four wagon groups can decrease the computing time by 72% and real-time simulations can be achieved by using eight cores and eight wagon groups.
A detailed fluid dynamics air brake system was developed, which includes subsystem models for pipes, wagon brake valves, and locomotive brake valves. The brake pipe model is based on conservation of mass and conservation of momentum; it was solved using the method of characteristics. A new efficient hose connection boundary condition that considers pressure loss was developed. The brake valve models have considered movements of valves and pistons as well as air flows between various volumes. A brake cylinder model with variable volume was also developed.
The model was validated using measured data from a test facility with 150 sets of connected wagon brake systems. Simulated results and measured results reached an agreement with the maximum difference of 15%. It is expected that these differences could be further reduced with further adjustment of parameters as studies continue. All important features of the air brake system have been simulated in detail.
Parallel computing schemes were developed to improve the computing speed. Computing time was compared for simulations with and without parallel computing. The computing time for the conventional sequential computing scheme was about 6.7 times slower than real-time. Parallel computing with four computer cores decreased the computing time by 70%. Real-time simulations were achieved by parallel computing using eight computer cores.
The authors gratefully acknowledge the financial support from the State Key Laboratory of Traction Power, Open Projects (TPL1504) and the National Natural Science Foundation of China (Grant No. 51575458). The authors are also grateful to Dr. Shulei Sun (China Academy of Engineering Physics) for providing experimental data. Acknowledgment is also made to Professor Wei Wei (Dalian Jiaotong University) for providing advice during the development of the air brake model. The editing contribution of Mr. Tim McSweeney (Adjunct Research Fellow, Centre for Railway Engineering) is gratefully acknowledged.
where and for calculating for , and and for calculating for . Figure 12 is provided as an example to demonstrate the process of numerical integration. Note that this figure only shows the case of subsonic flow as supersonic flow does not occur in brake pipes according to the authors' experience. Three main nodes are considered: , , and ; the objective is to determine and for . As all Riemann variables are already known for the three main nodes at the time point of , Riemann variables at and can be interpolated if and are known. As the values of Riemann variables are also the slopes of and in Fig. 12, therefore