NeTrainSim: a network-level simulator for modeling freight train longitudinal motion and energy consumption

Abstract

Although train modeling research is vast, most available simulation tools are confined to city- or trip-scale analysis, primarily offering micro-level simulations of network segments. This paper addresses this void by developing the NeTrainSim simulator for heavy long-haul freight trains on a network of multiple intersecting tracks. The main objective of this simulator is to enable a comprehensive analysis of energy consumption and the associated carbon footprint for the entire train system. Four case studies were conducted to demonstrate the simulator’s performance. The first case study validates the model by comparing NeTrainSim output to empirical trajectory data. The results demonstrate that the simulated trajectory is precise enough to estimate the train energy consumption and carbon dioxide emissions. The second application demonstrates the train-following model considering six trains following each other. The results showcase the model ability to maintain safe-following distances between successive trains. The next study highlights the simulator’s ability to resolve train conflicts for different scenarios. Finally, the suitability of the NeTrainSim for modeling realistic railroad networks is verified through the modeling of the entire US network and comparing alternative powertrains on the fleet energy consumption.

Introduction

The transportation sector is the largest consumer of total energy accounting for 26% of the US energy use in 2020 [1]. Railroad transportation results in a 75% cost reduction compared to other ground transportation modes and transports roughly 40% of long-distance freight volume [2]. In 2021, freight trains consumed a staggering 3,082 million gallons of diesel, leading to an estimated environmental impact of approximately 35.2 million tons of carbon dioxide (CO₂) emissions [3, 4].

The energy consumption of trains is influenced by various factors, including logistical, technical, and operational factors. Logistical factors are related to the trainload and network characteristics. Technical factors include vehicle physical characteristics such as the fuel type and aerodynamic parameters. Finally, operational factors include speed and driving dynamics [5].

The purpose of this paper is to describe an open-source simulator for heavy long-haul freight trains on a network where they interact with each other while producing valid instantaneous energy consumption estimates allowing for the estimation of the carbon footprint. The simulator, named NeTrainSim (Network Train Simulator), is built specifically for energy consumption and greenhouse gas emissions prediction of trains considering the main logistical, technical, and operational factors impacting them. Each car or locomotive in the train is considered a point mass positioned at the vehicle’s center of gravity with only a longitudinal degree of freedom while ignoring lateral and vertical dynamics. The reason behind limiting the degree of freedom to the longitudinal component is justified by an expectation of a significant reduction in the computational simulation time with minimum impacts on energy consumption estimates.

The resistance forces, consisting of the aerodynamics, rolling, curve, and grade resistance forces, corresponding to each locomotive and car in the train are modeled at their specific location on the track. Furthermore, the simulator takes into account additional inputs for the mathematical representation of the train motion such as the network structure, track characteristics, and train parameters. NeTrainSim allows locomotives to be distributed along the train length in three locations: the head, inside, and the end of the train. The distribution of cars with custom loads can easily be modified. NeTrainSim also includes a graphical user interface to facilitate the user’s experience.

With regard to the paper layout, initially, the train motion model is presented. Subsequently, a description of the NeTrainSim simulator is provided. Lastly, example case studies of four routes are presented to illustrate the simulator’s capabilities.

Literature review

While train simulation research is widespread, most existing simulators are unable to model instantaneous train movements at scale. Specifically, multiple train simulators typically ignore the instantaneous motion of the train to achieve scalability. Alternatively, detailed train simulators are developed to simulate longitudinal train dynamics (LTD) considering the motion of the train as a whole and/or any relative motion between vehicles in the direction of the train movement [6]. Simulators found in the literature are of two types: whole-trip simulators and sectional or short-trip simulators. As summarized by Ref. [7], whole-trip LTD simulators [8, 9] replicate one fixed-configuration train running on a fixed route. Whole-trip simulators such as those developed by Refs. [10, 11] are focused on calculating the in-train forces and their patterns with the vehicle connection system and draft gear behavior taken into consideration on a single track. Similarly, Qi et al. [12] provided a positioner model that optimizes the speed of the train to protect wagons from damage. According to Ref. [7], a drawback of these simulators is related to their lengthy simulation time due to the complexity of the involved models and computing strategies. The complexity of these models comes from the numerical solvers of differential equations—such as Runge–Kutta [13, 14], Park method [15], and others—that have been incorporated into these simulators.

A short-trip simulator has the same limitations as whole-trip simulators. Yet, they run relatively fast compared to their counterparts on account of their simplified dynamics models. While whole-trip simulators provide a more detailed assessment, short-trip simulators provide a microanalysis of a single train vehicle or the train as a whole [16]. Other researchers proposed discrete mathematical models for the simulation of specific train systems. Varazhun et al. [15] provided a predictive model for couplers’ forces in train cars due to electrodynamic braking. Li et al. [17] developed a different type of simulation; their model optimized the train trajectory, number of vehicles, and hauled weights based on the track profile. Finally, Wei and Lin [18] developed a simulator to predict pressure values in the air brake system.

Another simulator type is the one that was developed by the Federal Railroad Administration (FRA). The FRA-sponsored simulator is the Train Energy and Dynamics Simulator (TEDS). TEDS was developed for multiple purposes including conducting safety and risk evaluations, energy consumption studies, incident investigations, train operation studies, and ride quality evaluations. TEDS simulates the behavior of the train along the centerline of an ideal track with one degree of freedom (longitudinal motion) while discarding the vertical and lateral motion [19]. Despite the robustness of the simulator, it only simulates one train on a single track.

The FRA sponsored another simulator named ATTIF to perform accident investigation, train configuration evaluation, and assist in the training of train operators. The ATTIF simulator uses simplified nonlinear dynamics of railroad vehicles that allow for maintaining a fair degree of accuracy and a relatively short simulation time [20]. According to Ref. [21], ATTIF integrated a detailed multi-body dynamics coupler system model starting in 2012. In addition, the Train Dynamics and Energy Analyzer/train Simulator (TDEAS) was developed by the Chinese State Key Laboratory of Traction Power to perform detailed whole-trip longitudinal train dynamics and energy analyzes [9]. Wu et al. [21] summarized other simulators of friction draft gear modeling. Nevertheless, none of these simulators consider train energy consumption with respect to train forces and terrain topology at scale considering train interactions.

Lastly, Cipek et al. [22] convert and simulate a conventional 103-ton and 1.6-MW heavy-haul diesel-electric locomotive to a battery hybrid equivalent and derive fuel consumption and related greenhouse gas emissions models. The results of this research are an accurate representation of train fuel and energy consumption. However, as concluded, the model cannot be generalized but could be considered as a basis for later studies.

In the context of greenhouse gas (GHG) emissions, Kirschstein and Meisel [23] proposed a framework for developing GHG-emission models specifically tailored to rail freight transportation. The framework takes into account factors such as load characteristics, distance traveled, fuel consumption, and energy efficiency. It also considers the impact of different rail infrastructure conditions on emissions. It is noteworthy to mention that the authors’ research primarily focuses on the German context, suggesting that other countries should develop their own tailored factors to account for specific regional trains and infrastructure characteristics. In a study conducted by Graver and Frey [24], dynamometers were employed to measure the GHG emissions of trains. The authors showed the results for three locomotives. However, they did not develop a mathematical model for greenhouse gas emissions.

In summary, despite the presence of various comprehensive train simulation packages in the market, a tool capable of assessing energy consumption and carbon footprint generated by a network of trains remains conspicuously absent. These simulation tools, which are mostly limited to city- or trip-scale research, provide mostly micro-level simulations of network segments, resulting in a lack of full analysis and making intelligent decisions on a broader, system-wide scale difficult. This necessitates not only the simulation of train movements and operations but also an intricate understanding of the fuel efficiencies of different train models, the varying emission factors of different fuels, and the energy consumption patterns of the trains under diverse operating conditions. This narrow emphasis, along with the lack of a tool to estimate the carbon emissions produced by the locomotives inside the modeled network, highlights a critical gap in train modeling technology at the moment. This gap impedes a more robust understanding of the environmental impact of railway transportation systems and emphasizes the need for a more comprehensive simulation tool capable of modeling the complete network and supporting system-level analysis. The paper, and its resultant simulator NeTrainSim, present an attempt toward simulating an entire train system for energy consumption and carbon footprint emissions analysis, thus striving to fill this current market void.

Model novelty

This paper makes the following contributions to the current state-of-the-art in train modeling:

• It develops and introduces the first model for simulating train-following behavior that accounts for the external forces acting on the train while modeling the train as a sequence of point masses (a point mass for each locomotive and car that constitute the train) and capturing the significant latencies associated with the braking of long freight trains.

• NeTrainSim is the first open-source simulator that focuses on energy consumption considering six different powertrain technologies: diesel, biodiesel, their hybrid variants, electric, and hydrogen fuel cells.

• NeTrainSim offers advanced network modeling capabilities, allowing for the simulation of entire rail networks, including country-scale simulations.

• It provides a micro-simulation model for train motion at the network scale. The simulator supports a conflict resolution mechanism in the simulation network of many intersecting lines/tracks utilizing a FIFO (first-in, first-out) strategy.

• Given that we model the train motion second by second, the proposed simulator achieves scalability with minimum fidelity sacrifice in either modeling the train motion or in calculating the trains’ energy consumption and CO₂ emissions.

• NeTrainSim is developed in the C ++ programming language, which makes it more stable, faster, and able to handle large networks.

Mathematical model

The proposed train motion model is developed based on the 1992 Canadian National variation for resistance forces cited in [25], and refers to the models proposed in [26, 27] for a tractive force and train-following model. Table 1 shows the model variable definitions.

Table 1. Model variables Definition

Variable	Definition
${\tilde{a}}_n\left(t\right)$	Smoothed acceleration of train $n$ at instant $t$ (m/s2)
$A_{c,l}$	Frontal area of car $c$ or locomotive $l$ (m2)
$a_n\left(t\right)$	Acceleration of train $n$ at instant $t$ (m/s2)
${\text{CF}}_{\text{Fuel Type}}$	Conversion factor from energy consumption kWh to fuel quantity by fuel type
$C_{c,l}$	Track curvature of car $c$ or locomotive $l$ (°)
$c+l$	Number of cars and locomotives in the subject train
${\text{FCS}}_n\left(t\right)$	Fuel cell status of train n in time t (%)
${\text{FD}}_{\text{Fuel type}|n}\left(t\right)$	Fuel depletion of fuel type for train n in time t
$F_n\left(t\right)$	Tractive force of train $n$ at instant $t$ (N)
${\text{TGCD}}_n\left(t\right)$	Total Grid Consumption/Delivery of train n at time t (kWh)
$G_{c,l}$	Track gradient of car $c$ or locomotive $l$ (%)
$K_{c,l}$	Canadian National streamlining coefficient of car $c$ or locomotive $l$
$N_{\text{max}}$	Number of notches in the given locomotive
$P_l^{\text{max}}$	Maximum engine power of locomotive $l$ (kW)
$R_n\left(t\right)$	Resistive force of train $n$ at instant $t$ (N)
$T_n$	The time it takes to activate the brakes of the train plus the operator perception reaction time (s)
$m_{c,l}$	Total mass of car $c$ or locomotive $l$ (kg)
$m_c$	Total mass of car $c$ (kg)
$m_{\text{c}}^{\text{a}}$	Mass on single axle of car $c$ (kg)
$m_l$	Total mass of locomotive $l$ (kg)
$m_l^a$	Mass on a single axle of locomotive $l$ (kg)
${\text{SOC}}_n\left(t\right)$	The battery state of charge of train n at time t (%)
$s_n\left(t\right)$	Spacing from the rear bumper of train $n$ to the rear bumper of train $n-1$ and is computed as $x_{n-1}\left(t\right)-x_n\left(t\right)$ (m)
$s_n^{\text{j}}$	Train spacing at jam density (m). Equal to the length of train $n$ plus a buffer (taken to be 2m)
$t_1,t_2,t_3$	Calibration parameters for the throttle input level
$u_{\text{d}}\left(t\right)$	Desired speed or max speed a train can go by at instant t (m/s)
$u_{\text{f}}$	Track free-flow velocity (km/h)
$u_{\text{m}}\left(t\right)$	Train speed at maximum throttle at instant t (m/s)
$u_n\left(t\right)$	Speed of train $n$ at instant $t$ (m/s)
$x_n\left(t\right)$	Position of the back of train $n$ relative to the start of the trip (m)
${\lambda }^{\ast }$	Throttle level that equates resistance forces at instant t ($0\le \lambda \le 1$)
${\lambda }_n\left(t\right)$	Throttle level of train $n$ at instant $t$ ($0\le \lambda \le 1$)
$\mathrm{\Delta }t$	The solution time step (s)
$G\left(t\right)$	Grade of track at instant $t$ (%)
$N$	Notch number
$g$	Gravitational acceleration (9.8066 m/s2)
$m$	Train total mass $m=\sum _{c,l}{m}_{c,l}$ (sum of locomotive and car masses) (kg)
$\eta$	Mechanical efficiency of the transmission and gear
$\mu$	Coefficient of friction between the wheel and the track

Traction force model

The motion model proposed in Ref. [27] is based on a prior one developed by Ref. [26]. In both models, the throttle position is assumed to be hyperbolically proportional to vehicle speed. The throttle level increases up to a maximum with increasing speed and decreases when the speed approaches the desired speed [27]. Equation (1) demonstrates the proposed hyperbolic throttle function.

1

$${\lambda }_n\left(t\right)=\left\{\begin{array}{cc}\frac{\frac{u_n\left(t\right)}{u_{\text{d}}\left(t\right)}}{t_1+\frac{t_2}{1-\frac{u_n\left(t\right)}{u_{\text{d}}\left(t\right)}}+t_3\frac{u_n\left(t\right)}{u_{\text{d}}\left(t\right)}},0\le u_n\left(t\right)\le u_{\text{m}}\left(t\right)\hfill & \\ max\left(\frac{\frac{u_n\left(t\right)}{u_{\text{d}}\left(t\right)}}{t_1+\frac{t_2}{1-\frac{u_n\left(t\right)}{u_{\text{d}}\left(t\right)}}+t_3\frac{u_n\left(t\right)}{u_{\text{d}}\left(t\right)}},{\lambda }^{\ast }\right),u_{\text{m}}\left(t\right)<u_n\left(t\right)\le u_{\text{d}}\left(t\right)\hfill & \end{array}\right),$$where variables $t_1,t_2\text{and}{t}_3$ are calibrated parameters that were originally introduced in Ref. [26]. These parameters are calculated based on the fact that the full capability of the vehicle motor power is never used and only around 60% of the vehicle capacity is used. However, in trains, this does not apply. Accordingly, these parameters are calibrated to reflect the full usage of the train power as proposed in the model [27], of which parameters t₁, t₂, and t₃ equal 0.190, 0.152, and 0.050, respectively. However, these values were obtained for passenger trains. Freight trains use more aggressive throttle levels. Thus, these values were adjusted to be 0.001, 0.050, and 0.030, respectively. It should be noted that the user can alter these default values as needed.

${\lambda }^{\ast }$ is the minimum throttle that is allowed to overcome the resistance forces when the desired speed is reached. Equation (2) is used to calculate the value of ${\lambda }^{\ast }$ that allows the vehicle to remain at the desired speed.

2

$$\frac{min\left(\frac{1000\eta {\lambda }_n\left(t\right)P_l^{max}}{u_n\left(t\right)},\mu m_{l}g\right)-R_n\left(t,u_{\text{d}}\right)}{m}=0.$$

Unlike motor vehicles controlled by continuous throttle behavior functions, trains are controlled by discrete throttle notches, which results in incremental changes in throttle with running speed [27]. Accordingly, the continuous function in Eq. (1) is not to be directly applied to the train throttle forces without discretizing it first. To discretize Eq. (1), ${\left(\frac{N}{N_{\text{max}}}\right)}^2$ is used in the approach described in Ref. [6]. The rest of the procedure is addressed in Ref. [27].

NeTrainSim uses Eq. (1) to calculate the throttle level based on the desired speed. The desired speed is a variable set to the maximum speed the locomotive can theoretically achieve. This variable value within the simulation framework is user-specified, however, based on a comprehensive review of relevant literature, a default value of 120 km/h has been established.

Figure 1a shows the discretized throttle level based on eight train notches and the desired speed. The resulting discretization is aligned with the train number of notches in Table 2. The throttle level is then used to calculate the train tractive forces. Figure 1b shows the tractive forces at different notches. It is noticeable that this curve does not align with the curve in Ref. [28]. This is because Fig. 1b is a multiplication of power and the curve [28] is a multiplication of acceleration.

Fig. 1 [Images not available. See PDF.]

Engine characteristics: a throttle/notch level; b tractive force

Table 2. Trains characteristics used in scenarios I and II

Train characteristics	Scenario I	Scenario II
Track length (km)	162	322
Stopping stations at (km)	40, 42, 88, 150	40, 42, 88, 150
Transmission efficiency	0.98	0.82
Max locomotive power (kW)	3262.0	2445.9
Number of locomotives	3	11
Number of of axles per locomotive	6	6
Coefficient of friction	0.4	0.4
First locomotive $K_{\text{l}}$ value	24	24
Other locomotive $K_{\text{l}}$ value	5.5	5.5
Car $K_{\text{c}}$ value	5	5
Locomotive frontal area (m2)	14.8645	14.8645
Car frontal area (m2)	12.0774	11.1484
Number of cars	71	139
Number of car axials	4	4
Locomotive length (m)	22.3	23.0
Car length (m)	29.0	20.7
Locomotive weight (ton)	198	190
Car weight (ton)	44	100
Grade (%)	0–2.4	0–2.0
Curvature (%)	0	0–5

The net tractive force after overcoming the resistance forces is available to accelerate the train forward. The resistance forces change instantaneously on the track for each locomotive/car based on their attributes and location on the track. When the tractive forces are equal to the resistance forces, the train is unable to accelerate and travels at a constant speed. Alternatively, when the resistance force is higher than the tractive force, the train decelerates.

Resistance forces model

The modified Davis Equation coefficients have been updated to reflect modern trains as demonstrated in Eq. (3):

3

$$R_r=1.5+\frac{18N}{m_{c,l}}+0.03u_n\left(t\right)+\frac{K_{c,l}{A}_{c,l}{u}_n^2}{10,000m_{c,l}}.$$

Longitudinal motion model

This section is dedicated to proposing a model for train-following dynamics, which encompasses the integration of external force vectors, including resistive and throttle forces. Furthermore, the model accounts for the operator perceptual latency and the temporal lag inherent in the train’s braking system actuation. The resultant formulation yields a temporal profile of the train acceleration as a function of the inter-train spacing and the relative velocities of both the following and leading trains at discrete time intervals.

The tractive force on each locomotive is computed using Eq. (4). The model includes the basic tractive forces and the maximum force that can be sustained between the locomotive wheels and the track and includes a throttle function, as proposed by Ref. [26]. The throttle function is discretized as described previously. The max train acceleration $a_n^{\text{max}}\left(t\right)$ is computed using Eq. (5) as the difference between the total tractive force $F_n\left(t\right)$ and the total resistance forces $R_n\left(t\right)$ relative to the total mass $m_n$. The total tractive force is computed as the summation of the tractive forces on all the locomotives ($l$) using Eq. (4). The throttle input used in Eq. (4) is assumed to be the same for all locomotives.

4

$$F_{t|n}\left(t\right)=\sum _{l}min\left(\frac{1000{\eta }_n{\lambda }_n\left(t\right)P_l^{max}}{u_n\left(t\right)},\mu m_{l}g\right),$$

5

$$a_n^{\text{max}}\left(t\right)=\frac{F_n\left(t\right)-R_n\left(t\right)}{m_n}.$$

Models that calculate traction, braking forces, and resistive forces are attuned to the specific characteristics of each vehicle. These models are adaptable, permitting adjustments to accommodate various vehicle designs, reflecting the performance of different models under varied operational conditions.

Here $a_n^{\text{max}}\left(t\right)$ is train n maximum acceleration in m/s² at instant t, $F_{t|n}\left(t\right)\text{and}{R}_n\left(t\right)$ are the train tractive and resistance forces in Newtons at instant t, $m_n$ is train total mass in kg, ${\eta }_n$ is train mechanical efficiency of the transmission and gear $\left(0\le \eta \le 1\right)$, $P_l^{\text{max}}$ is the maximum locomotive engine power of train n in kW, $u_n\left(t\right)$ is velocity of train n at instant t in m/s, $\mu$ is coefficient of friction between the wheel and the track, $m_{\text{l}}$ is total weight of locomotives in kg, and $g$ is gravitational acceleration (9.8066 m/s²).

The resistance forces are computed using the Canadian National variation of the Davis equation for both the locomotives and rail cars as in Eq. (3). The gradient resistance force is added. Curve resistance is converted to an equivalent grade resistance by assuming that the unit resistance of a 1° curve is the same as the resistance of a 0.04% grade [29].

Hence the final resistance for each locomotive ($l$) or car ($c$) is$$1.5+\frac{18N}{W}+0.03u_n+\frac{K_{c,l}{A}_{c,l}{u}_n^2}{10,000W}+20\left(G_{c,l}\left(t\right)+0.04\left|C_{c,l}\left(t\right)\right|\right).$$

Given that the Davis equation generates the resistance force in lbs, the unit conversion (4.4482) is necessary to convert from units of lbs to Newtons. Equation (6) is the result of this conversion:

6

$$R_{\text{r}}=\frac{4.44822\times 1.10231}{1000}\sum _{c,l}{m}_{c,l}\left(1.5+\frac{16329.34}{m_{c,l}^a}+0.0671u_n\left(t\right)+\frac{48862.37A_{c,l}{K}_{c,l}{u}_n{\left(t\right)}^2}{m_{c,l}}+20\left(G_{c,l}\left(t\right)+0.04\left|C_{c,l}\left(t\right)\right|\right)\right).$$

The modeling of train deceleration considers a constant deceleration $d_{\text{des}}$, which is user-specified but typically set at 0.2 m/s². We use a simple linear train-following model to compute the safe spacing between trains at steady-state conditions, $s$ using Eq. (7):

7

$$s_n\left(t\right)=s_n^{\text{j}}+T_n{u}_n\left(t\right),$$$\text{where}{s}_n^{\text{j}}$ is the spacing when stopped, which is taken to be the length of the train n plus a buffer of 2 m; $T_n$ is the time it takes to activate the brakes plus the operator perception reaction time, and $u_n\left(t\right)$ is the train velocity. $T_n$ is estimated by

8

$$T_n=\frac{L_{\text{c}}^{max}}{u_{\text{s}}}+t_{\text{pr}},$$where $L_{\text{c}}^{\text{max}}$ is the longest distance the brake signal needs to travel from the controlling locomotive to the last car in the batch of cars that are controlled by that set of locomotives. The brake signal is assumed to travel at the speed of sound ($u_{\text{s}}$) taken to be 343 m/s. $t_{\text{pr}}$ is the operator perception reaction time (taken to be 4.5 s in this paper as an average of what was found in Ref. [30] but can be user-specified). Using Eq. (7), the terms are re-arranged to estimate the train following speed the next time step based on current spacing, as demonstrated in Eq. (9):

9

$${\tilde{u}}_n\left(t+\mathrm{\Delta }t\right)=\text{min}\left(\frac{s_n\left(t\right)-s_n^{\text{j}}}{T_n},u_{\text{f}}\right),$$$\text{where}{u}_{\text{f}}$ is the free-flow velocity of the track ahead of the train. The time-to-collision ($\text{TTC}$) is computed assuming the train continues at its current speed, as shown in Eq. (10).

10

$$\text{TTC}=min\left(\frac{s_n\left(t\right)-s_n^{\text{j}}}{max\left(u_n\left(t\right)-u_{n-1}\left(t\right),0.0001\right)},{\text{TTC}}_{max}\right).$$

The desired acceleration, at some time into the future using the spacing at time $t$ and incorporating it in the range policy presented in Eq. (9), is computed twice. First assuming the speed is achieved over a time interval $\text{TTC}$ (Eq. (11)) and the second is assumed to occur over a time interval $T_n$ (Eq. (12)).

11

$$a_{n,1-1}\left(t\right)=max\left(\frac{{\tilde{u}}_n\left(t+\mathrm{\Delta }t\right)-u_n\left(t\right)}{\text{TTC}},-\mu g\right),$$

12

$$a_{n,1-2}\left(t\right)=min\left(\frac{{\tilde{u}}_n\left(t+\mathrm{\Delta }t\right)-u_n\left(t\right)}{T_n},a_n^{max}\left(t\right)\right).$$

We then compute the train acceleration as a weighted combination of the two accelerations, where the term ${\beta }_1$ is computed using Eq. (14). The coefficient ${\beta }_1$ is a binary variable that is equal to zero when the acceleration is negative and equals one when the acceleration is either zero or positive. The first acceleration term is used for the train’s negative accelerations (decelerations) while the second term is used for positive accelerations.

13

$$a_{n,1-3}\left(t\right)=\left(1-{\beta }_1\right)a_{n,1-1}\left(t\right)+{\beta }_1{a}_{n,1-2}\left(t\right),$$

14

$${\beta }_1=\frac{a_{n,1-1}\left(t\right)+\left|a_{n,1-1}\left(t\right)\right|}{2\times max\left(\left|a_{n,1-1}\left(t\right)\right|,0.0001\right)}.$$

An alternate train acceleration is computed by taking the Lagrangian derivative (a vehicle-based derivative) of Eq. (9), as formulated in Eq. (15):

15

$$a_{n,1-4}\left(t\right)=max\left(min\left(\frac{u_{n-1}\left(t\right)-u_n\left(t\right)}{T_n},a_n^{max}\left(t\right)\right),-\mu g\right).$$

We then compute the train acceleration as a weighted combination of these two accelerations, where the term ${\beta }_2$ varies in the range $\left[0,1\right]$.. The first acceleration term ($a_{n,1-3}\left(t\right)$) ensures that the train spacing between it and the train ahead complies with the range policy presented in Eq. (9). The second acceleration term ($a_{n,1-4}\left(t\right)$) ensures that the train adjusts its speed to the speed of the train directly ahead of it.

16

$$a_{n,1}\left(t\right)={\beta }_2{a}_{n,1-3}\left(t\right)+\left(1-{\beta }_2\right)a_{n,1-4}\left(t\right).$$

The complete train longitudinal motion model is a modification of the Fadhloun-Rakha car-following model [26] that is formulated in Eq. (17). The first term computes the train acceleration when the speed of the train ahead of it is greater than or equal to its speed while the second term computes the train acceleration while approaching a slower-moving train. It ensures that the train attempts to decelerate at the desired deceleration level ($d_{\text{des}}$). Again the $\gamma$ term is a binary variable that is either 0 or 1. It equals zero if the speed of the train ahead is greater than or equal to the subject train speed and is 1 if the speed of the train ahead is less than the subject train speed.

17

$$a_n\left(t\right)=\left(1-\gamma \right)a_{n,1}\left(t\right)-\gamma a_{n,2}\left(t\right).$$

The other parameters are computed as

18

$$a_{n,2}\left(t\right)=min\left(\frac{{\left(u_n{\left(t\right)}^2-u_{n-1}{\left(t\right)}^2\right)}^2}{4{\left(max\left(s_n\left(t\right)-s_n^{\text{j}}-T_n{u}_n\left(t\right),0.0001\right)\right)}^2{d}_{\text{des}}},\mu g\right).$$

The $d_{\text{des}}$ in this equation is the desired deceleration level, which is assumed to be constant or could be retrieved from a brake model. In this research, we used a constant value of $0.2\text{m}/{\text{s}}^2$ and is user-specifiable.

19

$$\gamma =\frac{u_n\left(t\right)-u_{n-1}\left(t\right)+\sqrt{{\left(u_n\left(t\right)-u_{n-1}\left(t\right)\right)}^2}}{2\times max\left(\left|u_n\left(t\right)-u_{n-1}\left(t\right)\right|,0.0001,.\right)}.$$

When the train spacing is greater than $s_{\text{lad}}$, the movement of the train is assumed to be free of train interaction and is achieved using the vehicle dynamics alone, i.e., $a_n^{\text{max}}\left(t\right)$.

20

$$s_{\text{lad}}=s_{n-1}^{\text{j}}+T_n{u}_{\text{f}}+\frac{u_{\text{f}}^2}{2\times d_{\text{des}}}.$$

The acceleration is then constrained by the maximum jerk allowed ($j_{\text{max}}$), as shown in Eq. (21). This ensures that the train movement is smooth.

21

$$\left|a_n\left(t\right)\right|=min\left(\left|a_n\left(t\right)\right|,\left|a_n\left(t-\mathrm{\Delta }t\right)+j_{max}\mathrm{\Delta }t\right|\right).$$

The smoothed acceleration ${\tilde{a}}_n\left(t\right)$ is then computed using an exponential smoother, as demonstrated in Eq. (22). Here $\alpha$ is the exponential smoother. A smoothing factor value of 1.0 provides no smoothing and lower values provide more smoothing.

22

$${\tilde{a}}_n\left(t\right)=\alpha \times a_n\left(t\right)+\left(1-\alpha \right)\times {\tilde{a}}_n\left(t\right).$$

The train speed is then computed using the first-order Euler approximation, as formulated in Eq. (23):

23

$$u_n\left(t+\mathrm{\Delta }t\right)=max\left(min\left(u\left(t\right)+\tilde{a}\left(t\right)\times \mathrm{\Delta }t,u_{\text{f}}\right),0\right).$$

Energy and carbon emission models

As described in our previous work [31, 32], the quantification of CO₂ emissions generated during combustion depends on the molar quantities of each reactant. For fossil fuel combustion (C₁₂H₂₃) as an example, the balanced equation below (Eq. (24)) indicates that 2683.067 g of fully oxidized carbon (44.01 g/mol) are emitted by burning 1L of diesel. This number is obtained by multiplying the number of CO₂ moles in 1L—following Eq. (24)—which is 60.97 mol by the molar weight of CO₂. However, due to the incomplete oxidation of carbon in the locomotive engine, a small portion of the emissions exists as carbon monoxide (CO) and hydrocarbons (HC). As indicated by Ref. [24], it is observed that the resulting CO₂ emissions (with a molar mass of 44.01 g/mol) constitute an average of 95% of the fully oxidized CO emissions [24].

24

$$4{\text{C}}_{12}{\text{H}}_{13}+71{\text{O}}_2\to 48{\text{CO}}_2+46{\text{H}}_2\text{O}.$$

Figure 2 presents the dynamic changes in the ratio of carbon dioxide (CO₂) emissions compared to the total oxidation of carbon, as well as the notch number in diesel fuel. Specifically, when 1L of fossil diesel fuel undergoes complete combustion, it releases approximately 2559.5 g of CO₂ (53.04 mol/L). On the other hand, if biodiesel (C₁₆H₃₆O₂) with the same combustion efficiency is used, it generates approximately 2226.7 g of CO₂. This relationship is described by Eq. (25), which provides a mathematical model for the correlation between the proportion of CO₂ emissions and the position of the notch. However, since the model only explains 34.3% of the error, which is relatively small, we use a constant of 95% (average value as used earlier).

25

$$F_{{\text{CO}}_2}\left(t\right)=0.963+0.066·\text{Notch}-0.105·{\text{Notch}}^2.$$

Fig. 2 [Images not available. See PDF.]

The relative proportion of carbon dioxide (CO₂) content resulting from partial oxidation of carbon to complete oxidation of carbon, as a function of the notch position [24]

Train delay and number of stops estimation

The train delay is computed at each time step for each locomotive/car by comparing its travel time to its travel time if driven at the free-flow speed of the corresponding track it is on, as shown in Eq. (26). This is based on previous work done in the traffic domain and validated against empirical data [33]. The delay for a specific train/trip is then computed as the summation of the delays across all the time steps that constitute the trip. The total network delay is then computed as the summation of the total delay of all trains simulated.

26

$$d_n\left(t\right)=\frac{{\sum }_{c,l}\left(1-\frac{u_n\left(t\right)}{u_{f|c,l}}\right)\times \mathrm{\Delta }t}{c+l}.$$

Similarly, the number of stops the train incurs is computed based on work done in the traffic flow domain [34], as demonstrated in Eq. (27). This equation captures all partial stops incurred by the train each time step and then is summed up across all time steps to compute the number of stops experienced by the train. This is then summed up across all the trains to compute the total number of stops incurred across the network of trains.

27

$$s_n\left(t\right)=\left\{\begin{array}{cc}\frac{\sum _{c,l}\left(\frac{u_n\left(t-\mathrm{\Delta }t\right)-u_n\left(t\right)}{u_{f|c,l}}\right)}{\left|c+l\right|}\hfill & u_n\left(t-\mathrm{\Delta }t\right)>u_n\left(t\right)\\ 0\hfill & \text{otherwise}\end{array}\right).$$

Simulator description

Simulator logic

NeTrainSim simulates multiple trains on a given network. The network is defined as a graph connecting nodes with links of multiple intersecting lines. The simulator also allows the implementation of signals at specified nodes. The simulator is a motion-based longitudinal simulator in that it allows the modeling of both tractive and resistive forces acting on the trains, while also incorporating train-following models that govern their behavior when in the vicinity of each other. The simulator is a time-driven algorithm that calculates the movements of the different trains at each time step of the simulation. Once completed, a summary file is generated containing information pertaining to the train travel time, traveled distance, consumed energy, and fuel consumption.

Further, the simulator is designed to incorporate and model a mixed technology consists. In such a case, the energy consumption is assumed to be evenly distributed at the wheels of all active locomotives, taking into account their operational status. However, there may be instances where a locomotive power source is no longer able to provide energy, resulting in the locomotive becoming inoperative. In such cases, the energy consumption calculations do not consider the energy utilization of that locomotive. The simulation platform continuously updates its calculations based on the operational status of each locomotive, ensuring that the energy consumption estimations remain accurate and relevant.

This simulation technique enables the modeling of complicated mixed technology situations and gives significant insights into energy use patterns, assisting in energy consumption optimization. Furthermore, by constantly updating the computations based on real-time operating conditions, the NeTrainSim guarantees that its estimates are accurate, giving trustworthy data for decision-making. It should be noted that the model then computes the energy at the battery/tank, depending on the technology using the efficiency curves that were described earlier.

The simulator (Fig. 3) is divided into modules and each module handles a set of tasks. The network module handles the network calculations and defines the network structure. The train dynamics module defines the train characteristics, their paths, and their movement dynamics. The energy module handles energy consumption calculations with different energy sources. Lastly, the simulator module is the central component where all calculations are synchronized, and actual train movements are simulated.

Fig. 3 [Images not available. See PDF.]

Simulator schema

The network links are assumed to be linear (only in length calculations). This is because the simulator uses vectors to calculate the train coordinates and reduce the calculation time. Vectorization requires the links to be linear instead of curves since curves are composed of millions of approximate vectors.

Line segments are treated as a piecewise sequence of links. Each link has a constant grade, curvature, and speed limit. When these links are short relative to the train length, the train spans many links. Therefore, every car or locomotive has its specific grade, curvature, and maximum allowed speed. The train is not allowed to exceed the maximum speed of any of the train-spanned links. Furthermore, trains must reduce their speed before entering a link to ensure that they maintain a free-flow speed (speed limit) less than the train’s current speed.

NeTrainSim (Fig. 4) starts by setting the locations where the speed must be zero (e.g., for crew changes). If no stops are specified, which is the default, a stop is specified at the end of the route. The main driving point of the simulator is checking whether all trains reach their destination. The simulation ends when all trains have reached their destinations. A summary file is then generated along with an optional trajectory file. If at least one train does not reach its destination, the simulator first determines the trains for which the trip is still in progress and then runs the calculations specifically for those trains until they all reach their destinations.

Fig. 4 [Images not available. See PDF.]

Simulation flow chart

At each time step (Fig. 5), the simulator retrieves the grade, curvature, and free-flow speed for every unit in every train.

Fig. 5 [Images not available. See PDF.]

Time step simulation calculations

The simulator uses this information to calculate resistance forces. In addition, the simulator sets the maximum speed that each train can go based on the maximum allowable speed for all units in a given train. Simultaneously, the simulator calculates how far it is to the next stopping station, reduced-speed point, or train ahead. This ensures that the train reduces its speed appropriately without colliding with other trains.

All the gathered information is passed to the train dynamics module to calculate how much acceleration or deceleration is required. This speed is used to compute the incremental distance traveled during the current time step. This distance is added to the train’s cumulative traveled distance. Lastly, the energy consumption of the train is calculated based on the train characteristics as stated in Ref. [35]. All trains are advanced in the same manner.

The last step toward making NeTrainSim a comprehensive network simulator relates to the implementation of an effective conflict management strategy as conflict zones are quite common for trains. The adopted strategy consists of two main stages. First, any potential conflict zones are identified for each train prior to the start of the simulation. That is achieved by pre-defining the route that each train will take. Based on those routes, two possible scenarios could arise, namely: if the links defining the route of a specific train and the links defining the routes of the remaining trains are exclusive, the conflict management strategy will be deactivated for that specific train as the possibility of conflicts has been eliminated beforehand. Otherwise, it will be active for the entire trip of the train; thus, ensuring that conflicts with other trains are addressed appropriately as they arise. This two-tier activation strategy ensures a significantly faster simulation time for networks with no or few potential conflicts.

Generally speaking, conflict zones for trains, that can be encountered in the simulator, fall into one of three main categories as shown in Fig. 6:

• The first category occurs when two trains are approaching each other and have to share a single two-way track. In this case, the first train to reach either of the entry nodes of the conflict zone has priority to use the shared track while the other train is forced to stop until the track is cleared. If several trains are queued in both directions, then the trains in the same direction as the first crossing train will be cleared first. From a network perspective, this strategy is more effective in terms of minimizing the total average delay of the trains than a FIFO strategy. While the used conflict management strategy is basic in that it only accounts for the instantaneous positions of the trains in the decision-making process, the research team will complement it later on with an Eco-cruise algorithm and enhancements to the priority of train movements. The purpose of such an algorithm is to ensure that the trains clear the conflict zones in the most energy-efficient manner.

• The second category is quite similar to the first one with the main difference that there exist several links that the trains can utilize to cross the conflict zone. For instance, Fig. 6-case 2 illustrates a scenario in which the conflict zone could be cleared by traveling on one of three links depending on the direction of travel. The three links consist of two one-way tracks (one in each direction) as well as a bidirectional track that can serve trains traveling in either direction. The conflict management strategy presented earlier for the first category remains valid here with the main addition that the use of the two one-way tracks is prioritized over the use of the bidirectional link.

• It is noteworthy to mention, at this level, that the first and second categories concern scenarios in which any links present in the conflict zone connect its entry node to its exit node which can only occur when the train length is smaller than the link length. In other words, trains need to cross a single link to clear it. Hence, there is a need for a third category that addresses the scenarios in which the conflict zones span over several successive links. This last category, which could be cast as a generalization of the first two categories, can happen when a train spans over several links as shown in Fig. 6(case 3). Because of that, the adopted conflict management strategy includes a module for the correct identification of the extremities of any conflict zones in the simulated network.

Fig. 6 [Images not available. See PDF.]

Conflict zone resolution strategy

Simulator scalability

The simulator uses discrete time steps to update the speeds, accelerations, energy consumption levels, and all the dependent statistics of the different trains. The smaller the step size is, the more accurate the results are. Using a time step of three seconds reduces the total simulation time by ~ 60% (relative to the 1-s time-step case) while the resultant statistics are largely unchanged. Increasing the time step beyond three seconds creates more significant changes in the predictions.

Scalability trials, encompassing 1 to 200 train increments with 1 to 10-s time steps, demonstrated linear performance O(n) escalation to 50 trains, transitioning to quadratic O(n²) beyond that threshold. The 200-train scenario was executed in 5.87 h on an Intel© Core i7-8750H, 32 GB RAM laptop.

Case studies

Four scenarios are presented here (I, II, III, and IV), with train characteristics shown in Table 2. The first attempts to validate NeTrainSim against empirical train trajectory data. That is achieved by simulating the actual network on which the empirical trajectory was collected and setting the train characteristics similar to those of the train that completed the trip. The validation of the simulator can be, thus, achieved by comparing the simulated trajectory to the empirical field-based observations. The second scenario primarily demonstrates the performance of the train-following model implemented in the simulator in terms of regulating the longitudinal motion of the trains when following one another. In that regard, the second scenario involves six trains following each other. The second through sixth trains are shorter and lighter than the first, by 50% to ensure that they can catch up with the lead train.

Scenarios I and II are one-way tracks of lengths 162 and 322 km, respectively, with 4 intermediate stopping stations. These stations force the train to stop completely and then move again. The stops are distributed as shown in Table 2. The trains start and end their trips with a speed equal to zero. The tracks consist of 207 and 156 one-way links for scenarios I and II, respectively, with lengths varying between 0.3 and 9 km. Different grade, curvature, and maximum speed combinations are assigned to every link along the track (as shown in Figs. 8 and 12). In our case studies we assumed the coefficient of friction to be constant with a value of 0.4.

Next, we follow up with scenario III, which groups a total of six sample cases presented using a similar network configuration. Each of the presented examples aims to highlight a specific feature of the simulator in relation to its conflict resolution strategy and/or train-following algorithm. A simple network configuration with short links and low free-flow speeds is used to achieve that objective. A diagram of the considered network is presented in Fig. 7. It consists of four origin–destination zones (A, B, C, D), two intermediate nodes (E, F), and five links connecting the different nodes as shown with corresponding lengths of 10 km for link EF and 14.14 km for the others. It is noteworthy to mention, here, that the number of links connecting node E to node F can be varied between one and three links depending on the investigated case.

Fig. 7 [Images not available. See PDF.]

Diagram of the network for scenario III

Finally, the free-flow speed of all the links in the network is set to $10\text{m}/\text{s}$ which makes it easier to visualize how the conflict resolution strategy is working. The train characteristics used in this Scenario are the same as the ones used in scenario II.

Finally, scenario IV compares different powertrain technologies on the entire US train freight network. To facilitate an accurate comparison, the complete USA Freight Network was comprehensively simulated. The adopted methodology and results are described in detail in our previous work [36].

Scenario I

Figure 8 shows the train speed profile for the first scenario. The traveled distance is shown on the x-axis and the speed in m/s on the y-axis. The dotted line presents the field measurements. Note that at high-grade values (dashed line), the speed drops as the train decelerates as a result of the significant increase in the resistance forces. For instance, at a distance of 35 km (Fig. 8b), the speed drops from 25 to 22 m/s due to a grade of 2%.

Fig. 8 [Images not available. See PDF.]

Speed profile of the train in scenario I

Field data for a freight train on this alignment was available to the authors from a rail operator in North America. The speed profile from the field data is provided for comparison with the model’s results. The acceleration/deceleration decisions of the model logic are somewhat different from that of the train operator. In zone Fig. 8a, both our simulator and the driver showed similar behavior in accelerating; however, the operator is found to be less aggressive than our model prediction as described by Eq. (1). Figure 8c expressed a similar behavior except the driver is slightly more aggressive than in Fig. 8a. In Fig. 8b, both the driver and the simulator decreased their speed to accommodate the stopping station. Nevertheless, the simulator is found to be more forceful on the brakes than the driver. In Fig. 8d, the driver reduced his speed unlike what the simulator did.

Statistics results show that the simulated trajectory resulted in 10 min and 44 s of delay and 8.0 of stops for the leading train.

Figure 9 plots the simulated and empirical acceleration profiles as a function of the traveled distance. The figure shows that the operator continuously changes the train speed and thus there is more noise in the empirical data. However, there are spots like those highlighted in Fig. 9 where the acceleration decisions are similar.

Fig. 9 [Images not available. See PDF.]

Acceleration profile of the train in scenario I

The acceleration provided in Fig. 9 is the actual acceleration used to change the train speed. Another interesting acceleration profile refers to the one resulting from the application of the brakes to reduce the speed in order to not exceed the free-flow speed on a downgrade. This is another type of acceleration, referred to as virtual acceleration, which is used to regenerate energy, along with the observed deceleration, as indicated in Ref. [37]. Figure 10 shows this virtual acceleration profile. Figure 10a is an instance of applying the brakes at a downgrade section while maintaining the train speed.

Fig. 10 [Images not available. See PDF.]

Virtual acceleration profile of the train in scenario I

Figure 11 shows the rate of energy consumption (EC) of the train (further described in Ref. [37]). As shown, this rate is following the speed profile. The energy consumption rate is the highest when the train is accelerating from a speed of zero and is lowest when the train is decelerating. When the EC is below zero, this indicates the train is regenerating energy and storing it in its batteries. The model predicts the total energy consumed to be 10.15 MWh while the field data show an energy consumption of 10.58 MWh, which corresponds to a 4.5% difference. Furthermore, when the aforementioned train, maintaining identical configuration and weights, operates on diesel fuel, there is an approximate 150% rise in energy consumption, corresponding to 26.72 MWh or 2,685.3 L of diesel fuel. This increase in energy usage is associated with emissions of approximately 6,872.9 kg of carbon dioxide (CO₂) and 343.2 kg of carbon monoxide (CO) and hydrocarbons (HC).

Fig. 11 [Images not available. See PDF.]

Energy consumption profile of the train in scenario I

Scenario II

Similarly to scenario I, we start by presenting information pertaining to the leading train speed profile in Fig. 12 as a function of the distance traveled (solid line). The dashed lines present the speed limit (also known as the free-flow speed) on each of the links along the train path. The profile shows sections with high grades and decelerating speeds. For example, at a distance of approximately 14 km (Fig. 12a), the speed drops from 17 to 6 m/s due to an uphill grade of 2%.

Fig. 12 [Images not available. See PDF.]

Speed profile of the leading train (train 1) in scenario II

There are also sudden drops in speed like the one at approximately 64 km (Fig. 12b) because of stops at stations. Also, the train does not exceed the speed limit of any link. Moreover, the train reduces its speed before leaving a high-speed link and approaching a low-speed link as indicated at a distance of 90 km in Fig. 12c. Statistics results show the simulated trajectory resulted in 1 h and 20 min of delay and 11.8 of stops for the leading train.

The leading train acceleration profile is shown in Fig. 13. The maximum acceleration comes right after a complete stop, and it is relatively high when the train is increasing its speed. When the deceleration is constant, the jerk is near zero, due to the smoothing function in Eq. (22), which constrains the train ability to reach the maximum deceleration level as soon as the brakes are applied.

Fig. 13 [Images not available. See PDF.]

Acceleration profile of the leading train (train 1) in scenario II

Figure 14 shows the rate of energy consumption for the leading train in scenario II. The total energy consumed by the leading train is found to be 83.5 MWh. For the following trains, it is around 50.6 MWh. The same configuration of the leading train would consume 382 MWh with approximately 38,391 diesel liters, which is approximately equivalent to 98,261.8 kg of CO₂ and 4,913 kg of CO and HC.

Fig. 14 [Images not available. See PDF.]

Energy consumption profile of the leading train (train 1) in scenario II

Figure 15 shows the time–space diagrams of the different trains. The slope of each train trajectory allows the determination of the instantaneous speed at a particular time. The first train (train 1) moves freely. Its speed is limited only by the maximum speed and its propulsive and braking capabilities. The speed reductions are due to these limitations. The following trains are constrained by the lead train. The initial headways are 1000 s to allow each train to traverse a significant distance before being impeded by the train ahead of it. At 14 and 170 km, the resistance forces are large due to high grades; and train 1 slows down. Trains 2–6 slow at that location and follow each other at the minimum allowable headway. At 100 km, train 1 reaccelerates and the headways increase again until 200 km.

Fig. 15 [Images not available. See PDF.]

Train time–space diagram in scenario II

Figure 16 shows the headways between trains. Since train 1 does not follow another train, its headway is not shown. As can be seen, the headway trends are similar for trains 2–6 but displaced in time. At the beginning of each train trajectory, it travels at the maximum allowed speed until the headway is less than $s_{\text{lad}}$ as calculated by Eq. (20). At that point, the train reduces its speed to follow the train ahead. At approximately 600 min, train 1 reduces its speed because of the significant grade. The rest of the trains bunch up behind it as a consequence of their faster speeds. After this, train 1 re-accelerates and the following trains move freely. Once each train reaches its destination, the plot of its headway to the train ahead ceases to be plotted.

Fig. 16 [Images not available. See PDF.]

Train headway to leading train in scenario II

Scenario III

As mentioned earlier, scenario III consists of six sub-scenarios that aim to showcase the functioning of the three categories of conflicts, for which a management strategy is implemented in the simulator, using the network diagram presented in Fig. 7. For each of the simulated cases, the results consisting of a screenshot of the simulation, the traveled distances, and the speed profiles of the trains are presented in Fig. 17 through Fig. 22. The first four cases, relating to the first category of conflicts in which two trains are approaching each other and have to share a single two-way link (between node E and node F), are as follows:

• Scenario III.1.1: Two trains T₁ and T₂ traveling from zones A and C to zones D and B respectively are introduced in the network one after the other. The results, presented in Fig. 17, demonstrate the simulator conflict management strategy is working as expected. Given that train T₁ started its trip before train T₂. It was given priority to cross the conflict area EF. As confirmed by the presented distance and speed profiles, Train T₂ was forced to stop at node F before the shared track until it was cleared. At that point, it was allowed to continue proceeding toward its destination.

• Scenario III.1.2: This scenario is similar to the previous scenario with the only difference that train T₂ was introduced in the network first. Figure 18 confirms that the results of this simulation are consistent with the expected outcome that train T₁ is the one to stop at the entrance of the shared track.

• Scenario III.1.3: This scenario is also similar to scenario III.1.1. However, the trains were introduced in the network at precalculated times in such a way that train T₁ clears the shared link while train T₂ is reducing its speed but before completely stopping. The speed profile of Fig. 19 confirms that was the case as train T₂ reduced its speed from the free-flow-speed of 10 m/s to around 3.5 m/s before accelerating again when train T₁ cleared link EF.

• Scenario III.1.4: Same as scenario III.1.3 except that train T₁ is traveling from zone D to zone A and train T₂ is traveling from zone B to zone C. Again, the results were as expected as shown in Fig. 20.

Fig. 17 [Images not available. See PDF.]

Scenario III.1.1 results: a time–distance diagram; b time–speed diagram

Fig. 18 [Images not available. See PDF.]

Scenario III.1.2 results: a time–distance diagram; b time–speed diagram

Fig. 19 [Images not available. See PDF.]

Scenario III.1.3 results: a time–distance diagram; b time–speed diagram

Fig. 20 [Images not available. See PDF.]

Scenario III.1.4 results: a time–distance diagram; b time–speed diagram

Next, we present scenario III.2. This scenario is set up similarly to scenario III.1.1 except that two one-way links (one in each direction) connect nodes E and F. As expected, the two trains were able to cross segment EF simultaneously without any conflicts. That is confirmed by the overlap of the distance and speed profiles presented in Fig. 21. The observed overlap is quite normal given the symmetry of the network and the fact that the two trains started their trips at the same time. Finally, scenario III.3 illustrates the case where the conflict zone stretches over several successive links.

Fig. 21 [Images not available. See PDF.]

Scenario III.2 results: a time–distance diagram; b time–speed diagram

Figure 22 confirms that the simulator handled this type of conflict successfully given that train T₂ stopped at what was determined to be the entrance of the conflict area.

Fig. 22 [Images not available. See PDF.]

Scenario III.3 results: a time–distance diagram; b time–speed diagram

Scenario IV

The final scenario entails using the NeTrainSim simulator to model the entire US freight train network to quantify the fuel/energy consumption for different train technologies.

This scenario uses this scalability in demonstrating the effectiveness of the simulation framework in assessing diverse powertrain technologies on the energy consumption at the tank, as summarized in Table 3. The simulation’s accuracy is established by an initial comparison with actual energy consumption data of diesel trains from the Association of American Railroads [4], which documents an expenditure of 3.118 billion gallons of diesel fuel in 2021. These empirical figures corroborate the simulated values presented in Table 3, thereby affirming the reliability of the simulation with an error of only 0.9%. The table further extends the energy consumption analysis to encompass a spectrum of powertrain alternatives, maintaining consistent train configurations for an equitable comparison.

Table 3. US simulated energy consumption comparison using NeTrainSim

Type	Total consumption per year	Eq. energy per year (kWh)	Relative to diesel (%)	CO2 emissions ×106 ton/year
Diesel	3,089,551,518 Gal	116,370,395,749	100.0	7.88
Diesel hybrid	2,619,659,798 Gal	98,671,553,356	84.8	6.68
Biodiesel	3,323,887,151 Gal	116,433,053,998	100.1	7.37
Biodiesel hybrid	2,790,214,948 Gal	97,738,952,319	84.0	6.19
Battery electric	53,733,886,668 kWh	53,733,886,668	46.2	–
Hyrogen FC	2,947,102,078 kg	98,226,912,261	84.4	–
Catenary electric	49,020,429,884 kWh	49,020,429,884	42.1	–

In the first column, Table 3 lists the annual total energy consumption for the entirety of the US freight train network in units specific to each powertrain technology. The second column standardizes these figures into kilowatt-hours (kWh), to enable a straightforward comparative analysis. The third column quantifies the proportional energy usage of each alternative relative to the conventional diesel benchmark. Finally, the fourth column provides the annual CO₂ emissions for each powertrain technology, measured in million tons, thus offering insights into the environmental implications of each technology. For an in-depth exposition of these results, readers are directed to consult Ref. [36].

Conclusions

In this paper, we present NeTrainSim, a unique open-source multi-train simulator for fuel, energy, and greenhouse gas emissions prediction of diesel, biodiesel, their hybrid variants, hydrogen, and electric freight trains. NeTrainSim is a motion-based simulator that uses train-following strategies adapted from traffic flow theory in combination with train motion modeling to control the longitudinal motion behavior of the different trains. That is complemented by a conflict management strategy that allows the simulator to deal with the main types of conflicts encountered in railroad networks. Given that freight trains are very long, the model decomposes the train into its constituent locomotives and cars while computing the resistance forces on the train. The simulator outputs consist of different metrics such as the instantaneous accelerations, speeds, positions, fuel/energy consumption levels, CO₂ emissions, delays, and stops of all the trains in the simulated network. The tool is demonstrated to produce results similar to empirical data demonstrating that the model produces defensible results for travel times, speed trends, and energy consumption of the trains. As of now, the simulator does not capture the so-called energy consumption “at the well” or at the ultimate energy source. Nevertheless, the simulator can compare energy consumption characteristics for different energy technologies considering onboard energy use. The tool is demonstrated to be scalable with computational times in the O(n²), where n is the number of trains. In future studies, the simulator would benefit from incorporating a more advanced brake model and considering a variable coefficient of friction along the track.

Acknowledgements

The information, data, or work presented herein was funded in part by the Advanced Research Projects Agency-Energy (ARPA-E), U.S. Department of Energy, under award number DE-AR0001471. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

Data availability

The NeTrainSim source code can be accessed from https://github.com/VTTI CSM/NeTrainSim.

References

1. Nalley, S; LaRose, A. Annual energy outlook 2022; 2022; Washington, U.S Energy Information Administration:

2. Association of American Railroads (2021) Freight railroads & climate change. https://www.aar.org/wp-content/uploads/2021/02/AAR-Climate-Change-Report.pdf, Accessed 18 May 2022

3. United States Environmental Protection Agency (EPA) (2023) Fast facts: U.S. transportation sector greenhouse gas emissions 1990–2021. https://www.epa.gov/system/files/documents/2023-06/420f23016.pdf. Accessed 17 May 2022

4. Association of American Railroads. Railroad facts; 2022; New York, Watertown:

5. García-Álvarez, A; Perez-Martinez, P; González-Franco, I. Energy consumption and carbon dioxide emissions in rail and road freight transport in Spain: a case study of car carriers and bulk petrochemicals. J Intell Transp Syst; 2012; 17, 3 pp. 233-244. [DOI: https://dx.doi.org/10.1080/15472450.2012.719456]

6. Iwnicki, S; Spiryagin, M; Cole, C et al. Handbook of railway vehicle dynamics; 2019; 2 Boca Raton, CRC Press: [DOI: https://dx.doi.org/10.1201/9780429469398]

7. Wu Q (2017) Optimisations of draft gear designs for heavy haul trains. Dissertation, Central Queensland University

8. Cole C (1999) Longitudinal train dynamics: characteristics, modeling, simulation and neural network prediction for Central Queensland coal trains. Dissertation, Central Queensland University

9. Wu, Q; Luo, S; Cole, C. Longitudinal dynamics and energy analysis for heavy haul trains. J Mod Transp; 2014; 22, 3 pp. 127-136. [DOI: https://dx.doi.org/10.1007/s40534-014-0055-x]

10. Shabana, AA; Aboubakr, AK; Ding, L. Use of the non-inertial coordinates in the analysis of train longitudinal forces. J Comput Nonlinear Dyn; 2012; 7, 1 011001. [DOI: https://dx.doi.org/10.1115/1.4004122]

11. Kovalev, R; Sakalo, A; Yazykov, V et al. Simulation of longitudinal dynamics of a freight train operating through a car dumper. Veh Syst Dyn; 2016; 54, 6 pp. 707-722. [DOI: https://dx.doi.org/10.1080/00423114.2016.1153115]

12. Qi, Z; Huang, Z; Kong, X. Simulation of longitudinal dynamics of long freight trains in positioning operations. Veh Syst Dyn; 2012; 50, 9 pp. 1409-1433. [DOI: https://dx.doi.org/10.1080/00423114.2012.661063]

13. Chen, C; Han, M; Han, Y. A numerical model for railroad freight car-to-car end impact. Discret Dyn Nat Soc; 2012; 2012, 927592. [DOI: https://dx.doi.org/10.1155/2012/927592]

14. Jin, X; Luo, Y. The mathematic description of features of the friction type draft gears. Roll Stock; 2011; 49, pp. 1-4. (in Chinese)

15. Varazhun, I; Shimanovsky, A; Zavarotny, A. Determination of longitudinal forces in the cars automatic couplers at train electrodynamic braking. Procedia Eng; 2016; 134, pp. 415-421. [DOI: https://dx.doi.org/10.1016/j.proeng.2016.01.032]

16. Evans, J; Berg, M. Challenges in simulation of rail vehicle dynamics. Veh Syst Dyn; 2009; 47, 8 pp. 1023-1048. [DOI: https://dx.doi.org/10.1080/00423110903071674]

17. Li, W; Jiang, S; Jin, M. Multi-Objective optimization and weight selection method for heavy haul trains trajectory. IEEE Access; 2022; 10, pp. 41152-41163. [DOI: https://dx.doi.org/10.1109/ACCESS.2022.3167155]

18. Wei, W; Lin, Y. Simulation of a freight train brake system with 120 valves. Proc Inst Mech F J Rail Rapid Transit; 2009; 223, 1 pp. 85-92. [DOI: https://dx.doi.org/10.1243/09544097JRRT119]

19. Andersen DR, Booth GF, Vithani AR et al (2012) Train energy and dynamics simulator (teds): a state-of-the-art longitudinal train dynamics simulator. In: ASME 2012 Rail Transportation Division Fall Technical Conference, Omaha, pp 57–63

20. Sanborn GG, Heineman JR, Shabana AA (2007) A low computational cost nonlinear formulation for multibody railroad vehicle systems. In: ASME 2007 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, Las Vegas, pp 1847–1856

21. Wu, Q; Cole, C; Luo, S; Spiryagin, M. A review of dynamics modeling of friction draft gear. Veh Syst Dyn; 2014; 52, 6 pp. 733-758. [DOI: https://dx.doi.org/10.1080/00423114.2014.894199]

22. Cipek, M; Pavković, D; Kljaić, Z et al. Assessment of battery-hybrid diesel-electric locomotive fuel savings and emission reduction potentials based on a realistic mountainous rail route. Energy; 2019; 173, pp. 1154-1171. [DOI: https://dx.doi.org/10.1016/j.energy.2019.02.144]

23. Kirschstein, T; Meisel, F. GHG-emission models for assessing the eco-friendliness of road and rail freight transports. Transp Res Part B Methodol; 2015; 73, pp. 13-33. [DOI: https://dx.doi.org/10.1016/j.trb.2014.12.004]

24. Graver, B; Frey, H. Comparison of locomotive emissions measured during dynamometer versus rail yard engine load tests. Transp Res Rec; 2013; 2341, 1 pp. 23-33. [DOI: https://dx.doi.org/10.3141/2341-03]

25. American railway engineering and maintenance-of-way association (2021) AREMA manual for railway engineering

26. Fadhloun, K; Rakha, H. A novel vehicle dynamics and human behavior car-following model: model development and preliminary testing. Int J Transp Sci Technol; 2020; 9, 1 pp. 14-28. [DOI: https://dx.doi.org/10.1016/j.ijtst.2019.05.004]

27. Wang, J; Rakha, HA. Longitudinal train dynamics model for a rail transit simulation system. Transp Res Part C Emerg Technol; 2018; 86, pp. 111-123. [DOI: https://dx.doi.org/10.1016/j.trc.2017.10.011]

28. Spiryagin M, Wu Q, Cole C et al (2016) Advanced studies on locomotive dynamics behavior utilizing co-simulation between multibody and train dynamics packages. In: CORE 2016, Maintaining the Momentum, Conference on Railway Excellence, Melbourne

29. Hay, WW. Railroad engineering; 1991; New York, Wiley:

30. Brandenburger, N; Jipp, M. Effects of expertise for automatic train operations. Cogn Technol Work; 2017; 19, pp. 699-709. [DOI: https://dx.doi.org/10.1007/s10111-017-0434-2]

31. Aredah A, Fadhloun K, Rakha HA et al (2022) NeTrainSim: A longitudinal freight train dynamics simulator for electric energy consumption. In: Transportation Research Board Annual Meeting 2022, Washington DC

32. Ahn, K; Aredah, A; Rakha, H et al. Simple diesel train fuel consumption model for real-time train applications. Energies; 2023; 16, 8 3555. [DOI: https://dx.doi.org/10.3390/en16083555]

33. Dion, F; Rakha, H; Kang, YS. Comparison of delay estimates at under-saturated and over-saturated pre-timed signalized intersections. Transp Res Part B Methodol; 2004; 38, 2 pp. 99-122. [DOI: https://dx.doi.org/10.1016/S0191-2615(03)00003-1]

34. Rakha, H; Kang, YS; Dion, F. Estimating vehicle stops at undersaturated and oversaturated fixed-time signalized intersections. Transp Res Rec; 2001; 1776, 1 pp. 128-137. [DOI: https://dx.doi.org/10.3141/1776-17]

35. Wang, J; Ghanem, A; Rakha, H et al. A rail transit simulation system for multi-modal energy-efficient routing applications. Int J Sustain Transp; 2021; 15, 3 pp. 187-202. [DOI: https://dx.doi.org/10.1080/15568318.2020.1718809]

36. Aredah, A; Du, J; Hegazi, M et al. Comparative analysis of alternative powertrain technologies in freight trains: a numerical examination toward sustainable rail transport. Appl Energy; 2024; 356, 122411. [DOI: https://dx.doi.org/10.1016/j.apenergy.2023.122411]

37. Aredah A, Fadhloun K, Rakha H et al (2022) NeTrainSim: a network freight train simulator for estimating energy/fuel consumption. Preprints 2022080518. https://doi.org/10.20944/preprints202208.0518.v1. Accessed 17 May 2022