A Dynamic Model for the Study and Simulation of the Pantograph–Rigid Catenary Interaction with an Overlapping Span

In this paper, the authors present a mathematical and engineering model to optimally calculate the dynamic equation on the pantograph–catenary interaction when considering a rigid catenary with an overlapping span. The model starts from well-known methods adapted to the special features of rigid catenary. As a result, an algorithm for the integration of a dynamic equation based on explicit methods is provided. Moreover, from this algorithm, a reliable, efficient, and user-friendly software tool called RICATI is developed in order to approach the model to railway-based companies. The results show the usefulness of an application. such as RICATI, to check the behavior of the configuration initially established for a catenary, allowing solutions to be obtained for the problems encountered when simulating the passage of the pantograph (or pantographs), not only for the overlapping span but also for the entire catenary. That encourages us to continue future works.


Introduction
Rail transportation is fundamentally important in modern societies. As a result, the development of railway technologies that allow the efficient circulation of people and goods faces challenges related to a variety of technical problems. Addressing these challenges is the object of study in research centers and universities. Throughout the history of the railway, different types of traction have been proposed for the tractor units. Traction by electric power is currently the most widely used system, in which the power is mainly supplied via an overhead contact line or catenary, and the pantograph is the mechanism of the tractor unit used to capture this energy.
Two types of catenaries can be considered. Flexible catenary, or elastic catenary, is an extensional component (bending is almost negligible but has to be considered). On the other hand, a rigid catenary is a bending component (axial extension is negligible). The study of the pantograph-catenary dynamic interaction to achieve optimum operating conditions for railways is therefore a fundamental technological problem. For this purpose, the contact force must be as uniform as possible. Contact losses or take-offs must be avoided, for which it is necessary to establish optimal assembly parameters. Thus, there is significant research interest in the development of a mathematical model that allows realistic simulations, and that is also computationally efficient.
The problem of pantograph-catenary dynamic interaction is highly important when the pantograph interacts with the elastic catenary, which has been the subject of a wide range of publications (see references [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18]). This number of studies on flexible catenary is due to the wide use of this type of catenary in high-speed railways that require higher performance. However, this issue also has special relevance when the interaction occurs with a rigid catenary, as is the case in tunnels and underground roads; however, few scientific studies have been conducted for this problem. In the case of this power system, the element that transmits the electrical energy is not a wire but a rigid rail. Due to its weight, this rigid rail cannot be maintained parallel to the track using tension or by suspending it from another wire with an amount of deflection; rather, it is necessary to increase the number of supports on which it must be suspended to decrease the distance between them. For example, a typical span of a rigid catenary has a length of about 10-12 m, whereas the span for the flexible catenary is around 50-60 m.
It is well known that, in the cases of rigid catenary, the speed of the units is restricted to maximum values of around 100-120 km/h. A small number of studies have been published for this type of line (see references [19][20][21]). Thus, it is an important goal of scientific research to raise this maximum speed limit.
In addition, to achieve greater rigidity of the catenary, which may be a cause of a significant fluctuation of the contact force, the sudden variation of this force must be added when the pantograph circulates through the spans of transition between two series of spans. These variations are more frequent because the length of the spans is notably less than in the case of flexible catenaries. No publications in the scientific literature address this problem in a satisfactory and detailed manner. In the case of overlapping spans, it is also necessary to configure the beams with a special slope-shaped geometry that ensures that the flow of the pantograph is as smooth as possible between the sequences of spans.
In this study, we addressed the problem of dynamic pantograph-rigid catenary interaction with an overlapping span. This problem has significant practical importance and, as previously indicated, has not been closely studied in the scientific literature. For this purpose, we present a methodology based on the combination of well-known mathematical models. Based on this model, a computer application for the study and simulation of the dynamic interaction of the pantograph with a rigid catenary is developed. This model is applied to two catenary series with an overlapping span and several pantographs in circulation. In the first section of this paper, the mathematical foundations of the model are explained and the algorithm of the developed computer program is presented. In the subsequent section, the computational aspects and results of the simulations are discussed.
The main novelty of this work lies in the application of well-known methods to study the dynamic interaction between the pantograph and rigid catenary, and the establishment of a development methodology. Therefore, the methods used in this paper are not novel. However, the application of all of them together to solve the problem of pantographrigid catenary interaction is a novelty. Much of the work found in the literature has been developed for flexible catenary. However, in this paper, existing methods and algorithms have been adapted in order to solve the case of pantograph-rigid catenary interaction. It should be noted that the mechanical characteristics of a rigid catenary are different from those of a flexible catenary. As a result, a specific study is necessary for this type of catenary because it is not feasible to consider a flexible catenary with extreme rigidity In addition, a software tool is introduced that addresses the problem of pantographrigid catenary dynamic interaction practically and efficiently from the perspective of both the accuracy of the results and the response time.
To sum up, these results for rigid catenary, as well as the underlying application, constitute the novelty of this method.

Dynamic Equations of the Catenary Pantograph System
The pantograph-catenary system is initially composed of two subsystems interacting with each other under restricted conditions. The dynamic equations for a time instant t n are shown in Equation (1), where M is the mass matrix, assumed to be constant; C n is the damping matrix; K n is the stiffness matrix; ϕ n is the matrix of the constraint conditions; R n is the vector of external loads on the system; q n is the vector of coordinates that define the position of the system; and λ n is the vector corresponding to the restriction forces. All of the latter elements are assumed to be variables, either in totality or in some of their terms for each instant t n . This results in a non-linear system of differential equations of second order whose integration in time will yield the evolution of the coordinates q n and contact forces λ n . The subscript n represents the value of variables at the instant t n .
When considering two series of spans, as shown in Figure 1, it is necessary to partition the system coordinates and the different elements of the equation. In the case of two pantographs, as represented in Figure 1, this results in the following structure: where q c1 , q c2 , q p1 , and q p2 , represent, respectively, the coordinates of the two series of spans (subscript c for catenary) and the two pantographs (subscript p). A similar notation may be used for the stiffness matrix and the other terms of Equation (1). The mass and damping matrices will have a structure similar to that of the stiffness matrix according to Equation (2). The vector λ n of the restraining forces or contact forces is divided into terms corresponding to the internal forces between each of the pantographs and catenaries of the two series of spans. The subscript n represents the value of variables at an instant t n .

Dynamic Equations of the Catenary Pantograph System
The pantograph-catenary system is initially composed of two subsystems interacting with each other under restricted conditions. The dynamic equations for a time instant tn are shown in Equation (1), where M is the mass matrix, assumed to be constant; Cn is the damping matrix; Kn is the stiffness matrix; φn is the matrix of the constraint conditions; Rn is the vector of external loads on the system; qn is the vector of coordinates that define the position of the system; and λn is the vector corresponding to the restriction forces. All of the latter elements are assumed to be variables, either in totality or in some of their terms for each instant tn. This results in a non-linear system of differential equations of second order whose integration in time will yield the evolution of the coordinates qn and contact forces λn. The subscript n represents the value of variables at the instant tn.
When considering two series of spans, as shown in Figure 1, it is necessary to partition the system coordinates and the different elements of the equation. In the case of two pantographs, as represented in Figure 1, this results in the following structure: where qc1, qc2, qp1, and qp2, represent, respectively, the coordinates of the two series of spans (subscript c for catenary) and the two pantographs (subscript p). A similar notation may be used for the stiffness matrix and the other terms of Equation (1). The mass and damping matrices will have a structure similar to that of the stiffness matrix according to Equation (2). The vector λn of the restraining forces or contact forces is divided into terms corresponding to the internal forces between each of the pantographs and catenaries of the two series of spans. The subscript n represents the value of variables at an instant tn.

Beam Model for Rigid Catenary
In this work, the rigid catenary is modeled as a flexible beam. To derive a model for the catenary, the finite element method is used, as explained in references [22,23]. As the flexible beam is a simpler structure, the stiffness matrix is a band matrix with fewer elements which has a positive impact, not only on the computational cost in terms of calculation, but also in terms of storage, as it is treated as a sparse matrix.
The total length of the catenary is divided into a series of segments that behave as beams with flexural rigidity, according to the Euler-Bernoulli equation. In the discretization,

Beam Model for Rigid Catenary
In this work, the rigid catenary is modeled as a flexible beam. To derive a model for the catenary, the finite element method is used, as explained in references [22,23]. As the flexible beam is a simpler structure, the stiffness matrix is a band matrix with fewer elements which has a positive impact, not only on the computational cost in terms of calculation, but also in terms of storage, as it is treated as a sparse matrix.
The total length of the catenary is divided into a series of segments that behave as beams with flexural rigidity, according to the Euler-Bernoulli equation. In the discretization, each node presents two coordinates corresponding to the displacement and rotation, resulting in matrices of 4 × 4 for each element, according to Figure 2.
where yi and θi represent the displacement and rotation in the node i, respectively; L represents the length of the element; m represents its mass; p is the weight of the beam per unit of length; E is the elastic modulus of the beam material; and I is the inertia moment of the section. These matrices are repeated for the different elements of the two catenaries. It can be seen that the different terms of the mass and stiffness matrices, and the vector of external loads, are constant in this case. For the catenary damping matrix, Rayleigh damping is assumed according to references [22,23], wherein the said matrix is a linear combination of the stiffness and mass matrix, per the equation: where α and β are two constants that are determined from two assumed values of the critical damping for two different frequencies of interest. In our case, for 1 Hz, a critical damping percentage of 0.5% is assumed, and for 15 Hz, the assumed percentage is 1%. These values are obtained experimentally by trial and error. The experiment consists of the following. A video recording is made on two points along the line and the movement of the two points after passing the pantograph is recorded. Then, two points on the line coinciding with the real points are considered in the simulation program, and for two frequencies of interest (1 Hz and 15 Hz in this case) various damping percentages are A punctual mass matrix is assumed for the masses. Thus, the components of the coordinates and matrices for the beam element are: where y i and θ i represent the displacement and rotation in the node i, respectively; L represents the length of the element; m represents its mass; p is the weight of the beam per unit of length; E is the elastic modulus of the beam material; and I is the inertia moment of the section. These matrices are repeated for the different elements of the two catenaries. It can be seen that the different terms of the mass and stiffness matrices, and the vector of external loads, are constant in this case. For the catenary damping matrix, Rayleigh damping is assumed according to references [22,23], wherein the said matrix is a linear combination of the stiffness and mass matrix, per the equation: where α and β are two constants that are determined from two assumed values of the critical damping for two different frequencies of interest. In our case, for 1 Hz, a critical damping percentage of 0.5% is assumed, and for 15 Hz, the assumed percentage is 1%. These values are obtained experimentally by trial and error. The experiment consists of the following. A video recording is made on two points along the line and the movement of the two points after passing the pantograph is recorded. Then, two points on the line coinciding with the real points are considered in the simulation program, and for two frequencies of interest (1 Hz and 15 Hz in this case) various damping percentages are assumed. The results of the simulation are compared with those obtained in the recording; it is then verified that the percentages given in the article are the ones that most closely match the real results. Figure 3 shows the details and geometric characteristics of a section of Furrer rigid beam, which is used in Spanish railways. The figure shows a normal section and a section with a joint flange. It can be seen that the normal beam is composed of two elements: the body of the beam itself, which is made of aluminum; and the contact wire, made of copper, through which the electric energy flows [20,21]. Obviously, Figure 3 consists of two elements: an aluminum beam and a copper wire. As we have two materials, we should have two moments of inertia and two elastic moduli. However, according to [24], it is possible to obtain an equivalent single beam representing the composite beam. This approach is considered in this paper in order to simplify the model. Figure 3 shows the details and geometric characteristics of a section of Furrer rigid beam, which is used in Spanish railways. The figure shows a normal section and a section with a joint flange. It can be seen that the normal beam is composed of two elements: the body of the beam itself, which is made of aluminum; and the contact wire, made of copper, through which the electric energy flows [20,21]. Obviously, Figure 3 consists of two elements: an aluminum beam and a copper wire. As we have two materials, we should have two moments of inertia and two elastic moduli. However, according to [24], it is possible to obtain an equivalent single beam representing the composite beam. This approach is considered in this paper in order to simplify the model. As in the treatment by finite elements, according to Equation (3), we assume that the beam is composed of a single material. It is necessary to determine the equivalent values for the section of the area, the position of the center of gravity, the moment of inertia, and the density, as if this section is formed by a single material. Thus, assuming that the entire section is exclusively aluminum, we first define the relationship between the elastic modulus or Young's modulus, where the subscript 1 corresponds to copper, and subscript 2 to aluminum. This notation is also adopted for the remainder of the properties, and thus we have:

Properties of the Rigid Catenary Section and the Supports
In Equation (5), m describes the relationship between the moduli. To simplify the notation, the subscripts are excluded.
Equation (6) shows the area of the equivalent section Aeq; the new position of the center of gravity of section G with respect to the position of the center of gravity of the aluminum section G2 yeq; the moment of inertia of the equivalent section with respect to the x-axis that passes through G Ixeq; and the density of the equivalent section: As in the treatment by finite elements, according to Equation (3), we assume that the beam is composed of a single material. It is necessary to determine the equivalent values for the section of the area, the position of the center of gravity, the moment of inertia, and the density, as if this section is formed by a single material. Thus, assuming that the entire section is exclusively aluminum, we first define the relationship between the elastic modulus or Young's modulus, where the subscript 1 corresponds to copper, and subscript 2 to aluminum. This notation is also adopted for the remainder of the properties, and thus we have: In Equation (5), m describes the relationship between the moduli. To simplify the notation, the subscripts are excluded.
Equation (6) shows the area of the equivalent section A eq ; the new position of the center of gravity of section G with respect to the position of the center of gravity of the aluminum section G 2 y eq ; the moment of inertia of the equivalent section with respect to the x-axis that passes through G I xeq ; and the density of the equivalent section: The details of the section of the rigid catenary section are shown in Figure 3, and the values for the Furrer beam are taken from reference [20,21]. After applying Equations (5) and (6) for the materials and sections of the beam, respectively, moments of inertia for the normal section of I xeq = 4,228,994 mm 4 and the flanged section of I xeq = 5,421,787 mm 4 are obtained. More details can be found in reference [20]. In addition, the spans of the beam are delimited by the supports, the detail of which is shown in Figure 3. The effect of the supports on the dynamic equations is assumed to be equivalent to that of a linear spring with stiffness k s = 1 × 10 5 N/m.

Pantograph Model
The pantograph is an articulated system whose purpose is to capture the electrical energy circulating through the contact wire and transmit it to the traction unit. To take into account its effect on dynamic equations, it is usually modeled as a set of masses, springs, and shock absorbers, on which different forces can act: pushing force, aerodynamic force, friction, etc. The manufacturer usually provides the values of these parameters. Figure 4 shows a drawing of a pantograph interacting with two rigid catenaries and some models of point masses. In the proposed simulation, a pantograph of a single head mass is assumed, whose properties correspond to those in the European standard EN50318 (see reference [25]), as shown in Table 1.
normal section of Ixeq = 4,228,994 mm 4 and the flanged section of Ixeq = 5,421,787 mm 4 are obtained. More details can be found in reference [20]. In addition, the spans of the beam are delimited by the supports, the detail of which is shown in Figure 3. The effect of the supports on the dynamic equations is assumed to be equivalent to that of a linear spring with stiffness ks = 1 × 10 5 N/m.

Pantograph Model
The pantograph is an articulated system whose purpose is to capture the electrical energy circulating through the contact wire and transmit it to the traction unit. To take into account its effect on dynamic equations, it is usually modeled as a set of masses, springs, and shock absorbers, on which different forces can act: pushing force, aerodynamic force, friction, etc. The manufacturer usually provides the values of these parameters. Figure 4 shows a drawing of a pantograph interacting with two rigid catenaries and some models of point masses. In the proposed simulation, a pantograph of a single head mass is assumed, whose properties correspond to those in the European standard EN50318 (see reference [25]), as shown in Table 1.     Thrust force applied to the base: F emp = 120 N According to the cited EN50318 standard, to establish the contact model with the catenary, it is necessary to add a contact element of zero mass with a spring to the head mass of the collector m 2 . Furthermore, because the pantograph can interact with two catenaries in our case, we assumed two contact elements with their corresponding stiffening springs, initially k 3 = k 4 = k cont , although these values can be canceled out when take-off occurs, as explained below. In addition, it is assumed that each mass has a coordinate corresponding to its vertical displacement, resulting in four coordinates, including the contact elements, for the standard model. In the case in which the pantograph has two head masses, two contact elements must be added to each mass with the corresponding springs, according to Figure 4. Based on the above, the matrices of mass, rigidity, and the vector of external forces corresponding to the point mass model of the norm are: The damping matrix has a structure similar to that of the stiffness matrix. It should also be noted that the terms of the stiffness matrix are not constant because, when take-off occurs, the contact disappears and the stiffness of the contact element k 3 or k 4 is canceled out. Furthermore, in the pantograph path, there are intervals in which the pantograph interacts only with a single catenary. In contrast, when the pantograph crosses the common transition zone, it interacts with the two catenaries, which also modify the structure of the stiffness matrix, as explained in later sections. In addition, in the data provided by the manufacturers, non-linearities can appear due to variations in the stiffness and damping, or the presence of friction forces. The mass matrix is always constant. Notice that the mass matrix is rank-deficient, meaning it is not possible to calculate its inverse. In this case, only the inverse of non-zero diagonal elements is calculated. More explanations about the calculation of the inverse of mass matrix will be given in Section 2.7.

Contact Pantograph/Catenary Model
The systems of pantograph and catenary are dependent but interact with each other under constraint conditions. In particular, the pantograph contact element slides on the contact wire of the catenary. Thus, a relationship between the positions is established with the pantograph contact element, represented by the contact element, and the position of the nodes of the catenary corresponding to the discretization using finite elements. In addition, this relationship varies depending on the time as the pantograph advances. The relationship between the coordinates of the catenary and the contact element of the pantograph corresponds to the term ϕ n that appears in the dynamic Equation (1). We suppose the distribution of the contact using a law that facilitates the sliding of the pantograph contact element along with the nodes of the catenary without any abruptness, as represented in Figure 5. This situation is also equivalent to assuming that the contact force is not punctual, but it is distributed over the surface of the pantograph contact element, according to the supposed distribution law (see Figure 5). For convenience, an Oxy axis system is assumed which moves with the pantograph, centered on the surface of the pantograph contact element. The variable x represents the position of the different points of the catenary along the horizontal axis, in contact with the plate. If L is the length of the friction surface, y3 is the coordinate corresponding to the vertical position of the pantograph contact element, and yc(x) is the vertical position of the different points of the catenary, then y = f(x), the function of contact distribution, will be fulfilled as follows: with the condition: Equations (8) and (9)  For convenience, an O xy axis system is assumed which moves with the pantograph, centered on the surface of the pantograph contact element. The variable x represents the position of the different points of the catenary along the horizontal axis, in contact with the plate. If L is the length of the friction surface, y 3 is the coordinate corresponding to the vertical position of the pantograph contact element, and y c (x) is the vertical position of the different points of the catenary, then y = f(x), the function of contact distribution, will be fulfilled as follows: with the condition: Equations (8) and (9) allow the position of the contact element of the pantograph to be obtained, when the configuration of the catenary is known at time t n , as a weighted average  (10) This results in a similar equation for (y 4 ) n . For f(x), we assumed a polynomial function defined as follows: The vertical position of any point of the catenary on the rubbing surface can be established from the position of the discretization nodes and the interpolation functions according to the finite element method, finally yielding the following relationship:

Formulation of Restriction Conditions
Equation (12) should be repeated for all system restrictions, corresponding to all pantograph-catenary contacts. The number of restrictions will depend on the number of possible contacts. Thus, for the pantograph with one head mass, two possible contacts exist because the pantograph can interact with up to two catenaries when crossing the common contact zone.
If the pantograph has two head masses, there will be four possible contacts, and if two pantographs are considered in the simulation, the restriction conditions must be repeated for the second pantograph. Thus, Equation (12) is a matrix equation whose number of rows corresponds to the number of restriction conditions and whose number of columns corresponds to the number of degrees of freedom.
To simplify the model, we assume a single pantograph with a single head mass. Because there can be up to two contacts, the matrix of constraint conditions expressed for a moment t n has the form: where the row vector ϕ 1n represents the restriction conditions of contact element 1 on catenary 1, and ϕ 2n represents the restriction conditions of contact element 2 on catenary 2, according to Figure 5. In the formulation of Equation (13), the following cases are considered, as shown in Figure 6. 1.
The pantograph circulates through catenary 1 without reaching the common transition zone of the two catenaries, defined from the overlapping spans. In this case, interaction only exists with catenary 1 through contact element 1, and contact element 1 is activated and contact element 2 is deactivated.

2.
The pantograph circulates through catenary 2, having exceeded the common transition zone. In this case, there will only be interaction with catenary 2 through contact element 2; thus, contact element 2 is activated and contact element 1 is deactivated.

3.
The pantograph circulates in the common transition zone between the two catenaries. In this case, the pantograph interacts with the two catenaries, and the two contact elements are activated.
In the first case, the vector ϕ 1n is calculated according to Equation (8), whereas the position of contact element 1, (y 3 ) n , is calculated according to Equation (10). Because contact element 2 is deactivated, the restriction condition ϕ 2n is formulated by requiring that the position of contact element 2, (y 4 ) n , and the position of the head mass, (y 2 ) n , at time t n , be coincident; that is: The second case is similar to the first, but with alternating conditions: ϕ 1n is obtained by imposing the condition that the position of contact element 1, (y 3 ) n , and the position of the head mass, (y 2 ) n , at instant t n are coincident: Matrix ϕ 2n and the position of contact element 2, (y 4 ) n , are calculated with equations similar to Equations (8) and (10). Finally, it must be taken into account that, regardless of the catenary with which the interaction occurs, the contact force can be canceled out because of a take-off or can appear suddenly after a take-off when a new coupling occurs. These circumstances must be considered when integrating Equation (1), as explained in the following section. moment tn has the form: where the row vector φ1n represents the restriction conditions of contact element 1 on catenary 1, and φ2n represents the restriction conditions of contact element 2 on catenary 2, according to Figure 5. In the formulation of Equation (13), the following cases are considered, as shown in Figure 6.
1. The pantograph circulates through catenary 1 without reaching the common transition zone of the two catenaries, defined from the overlapping spans. In this case, interaction only exists with catenary 1 through contact element 1, and contact element 1 is activated and contact element 2 is deactivated. 2. The pantograph circulates through catenary 2, having exceeded the common transition zone. In this case, there will only be interaction with catenary 2 through contact element 2; thus, contact element 2 is activated and contact element 1 is deactivated. 3. The pantograph circulates in the common transition zone between the two catenaries.
In this case, the pantograph interacts with the two catenaries, and the two contact elements are activated. In the first case, the vector φ1n is calculated according to Equation (8), whereas the position of contact element 1, (y3)n, is calculated according to Equation (10). Because contact element 2 is deactivated, the restriction condition φ2n is formulated by requiring that the position of contact element 2, (y4)n, and the position of the head mass, (y2)n, at time tn, be coincident; that is:

Integration of Dynamic Equations
To proceed with the integration of the differential Equation (1), we applied the explicit method of the central differences. Explicit methods are more suitable than implicit methods for dealing with non-linearities, in the particular case of the pantograph-catenary interaction [23], because they correctly take into account the variation of the stiffness matrix. Non-linearities occur due to debonding and coupling. Due to the obvious stability problems of explicit methods, the integration step must be carefully adjusted, which is undertaken in this development. The critical integration step that conditions the stability can be estimated as the element length divided by the wave propagation speed, as shown in Equations (13.10)-(16) (p. 401) of [23]. Regarding the element length, in this work all elements have the same length. If this is not the case, the shortest element would have to be taken for the calculation of the step. With respect to the propagation speed, it depends on the mechanical characteristics of the beam section. However, the authors have deduced the critical integration step by trial and error so that, when the integration step is above the critical time, the algorithm degenerates. Otherwise, the algorithm gives reasonable results. Therefore, different increments of t have tried to determine the critical value below when the algorithm is stable, and the results are feasible.
Thus, for this proposal, the speed and acceleration in the coordinates are initially approximated for a moment of time t n , and an interval ∆t, in accordance with: . q n = 1 2∆t (q n+1 − q n−1 ) ..
Substituting these equations into Equation (1) yields the following expression: C n q n+1 = R n − K n q n − ∅ t n λ n + 1 ∆t 2 M(2q n − q n−1 ) + 1 2∆t C n q n−1 (17) Equation (17) allows the coordinates q n+1 at time t n+1 to be obtained from the elements of the equation and the coordinates at the previous instant t n . However, for this calculation, it is necessary to solve a system of linear equations because, although the mass matrix M is diagonal, the same does not occur with the damping matrix C n . After all, we assumed a Rayleigh damping model for the catenary according to Equation (4). Although it is possible to assume a diagonal matrix for C n , we preferred to use the Rayleigh model because we think it is a more realistic approach. The same problem exists for the pantographdamping matrix.
One of the main advantages of the explicit integration methods is that the integration variables can be obtained directly without having to solve a system of equations, resulting in a more computationally efficient integration. However, for this purpose, it is necessary to slightly modify the dynamic equation by correcting the velocity vector in the coordinates by half a step, which results in: Equation (18) involves the calculation of the inverse of a rank-deficient mass matrix. However, the algorithm actually works by firstly calculating the position of the nodes, for which the values of the mass matrix are non-zero following the method of central differences. This method appears in [23]. Thus, to calculate the inverse of a diagonal matrix of non-zero elements, it is only necessary to calculate the inverse of the particular element m i , that is, 1/m i .
Once these positions are calculated, the contact forces are calculated for which the mass matrix is no longer necessary following Hooke's law where the nodes, which are the ends of the contact spring, and the stiffness of the spring are known.
In Equation (18), these terms appear together although they can be calculated, as mentioned above, in a decoupled manner.
According to reference [23], this approach introduces a small error that can be ignored in structural systems with low damping, as seen in our case. Next, the velocity and acceleration in the coordinates are approximated by the following equations: ..
Substituting the previous equations into Equation (18) and solving for the vector q n+1 yields: In this case, it is possible to calculate the vector q n+1 directly without the need to solve a system of equations because the mass matrix is diagonal. However, some variables remain to be determined: the position of the contact elements (y 3 ) n+1 and (y 4 ) n+1 , and contact forces (λ 1 ) n+1 and (λ 2 ) n+1 . The mass of the contact elements is null and Equation (20) cannot be applied; however, the positions of the catenary nodes in the pantograph environment are known at t n+1 . Thus, when the contact elements are activated, we use Equation (10) particularized in t n+1 for (y 3 ) n+1 and a similar equation for (y 4 ) n+1 . When the contact element is deactivated, its position coincides with that of the head mass; in this case, we will use Equation (14) or (15).
For the contact forces in t n+1 , these forces are calculated according to the deformation in the corresponding contact spring, according to the expressions: As can be seen in the previous equations, an additional correction term is added, given by the variables v 1 and v 2 , to the theoretical deformation of the spring given by the difference in the positions of the contact element and the head mass. This term takes into account the effect of the slope in the sectioning span as a consequence of the previous configuration of the beam. The values of functions v 1 and v 2 are discussed in the following section.
Because the contact force is a compression force, each integration step in Equation (21) must be solved by checking that the forces are negative. In this case, the value of the force is correct and the corresponding stiffness of the contact spring k 3 or k 4 is the theoretical value k cont of the stiffness matrix in Equation (7), and this value is maintained for the next integration step. However, if a value of the force is positive, it means that the spring works in traction, and therefore a take-off, has taken place. In this case, the contact force becomes zero and the corresponding rigidity is annulled for the next step. This is equivalent to eliminating the contact catenary-pantograph. Based on the above, the algorithm for detecting take-off at instant t n+1 is: According to Equation (22), it takes into account the possibility of contact and take-off by considering the impact between the catenary and the pantograph. Thus, for catenary 1, if the contact force is negative or equal to 0 ((λ 1 ) n+1 ≤ 0), then there is contact, and the stiffness matrix is configured according to this condition for the next step. If the contact force is positive, no contact is considered, the contact force is 0, and the stiffness matrix is set for the next step with that condition.
For the second catenary, it is the same procedure ((λ 2 ) n+1 ≤ 0). It is evident that, when the pantograph does not circulate through the common transition area, some contact elements are deactivated, and the corresponding force will be annulled. In addition, when circulating through the common transition area, the two contact elements are activated, and take-offs and couplings can be detected in all cases. To initialize the integration algorithm, it is necessary to provide a value for the coordinates and their derivatives, and for the contact forces in the zero instant. It is assumed that, at the initial instant, the velocities and accelerations in the coordinates are zero: Substituting these values into Equation (18) results in a linear system that can be resolved for q 0 and λ 0 , as follows: Equations (19)- (22), which are based on the method of explicit integration of central differences, make it possible to obtain the system variables directly by simply treating the non-linearities caused by take-offs and the effect on the overlapping, thus adapting the different variable elements of the dynamic equation at each integration step.

Functions for Correcting the Contact Position Due to the Effect of Overlapping Span
Although the assembly of a rigid catenary canton is undertaken using a series of spans in which the beams are assumed to be straight, in the case of overlapping spans, it may be advisable to previously deform the beam by giving it a sloped shape, and thus obtain the smoothest possible transition between cantons. This makes it difficult to simulate because the finite element model used is based on the hypothesis of straight beams. To overcome this disadvantage, it is assumed that the rigidity of the beam is not altered by the previous deformation. Thus, the integration of the differential Equation (18) is carried out as if the catenary is configured as a perfectly straight beam, both in the normal spans and in the section span. Then, to take into account the previous deformation of the beam and its influence on the contact force, the vertical position of the contact element calculated from Equation (10) is supplemented by the correction functions v 1 and v 2 . These functions take into account the variation of the position of these elements to then calculate the contact forces according to Equations (21) and (22). These functions depend on the geometry of the deformation and their effect is only taken into account when the pantograph is circulating through the section slopes. In the simulation proposed in the work, a linear slope is also assumed due to its simplicity; however, this methodology can also be adapted to any beam configuration with slight modifications.
In the case of a linear slope, the variable x defines the horizontal displacement of the pantograph from the beginning of the route, x 1 is the position where the ascending slope begins in the canton of outgoing catenary 1, and x 2 is the starting position of the descending slope in the canton of incoming catenary 2, with x 1 > x 2 , as shown in Figure 7. Let L 1 be the length of catenary slope 1, L 2 be the length of catenary slope 2, and h 1 and h 2 be the respective heights at the end of each slope. This is compiled for v 1 as: and for v 2 as: Appl. Sci. 2021, 11, x FOR PEER REVIEW 13 of 27 the length of catenary slope 1, L2 be the length of catenary slope 2, and h1 and h2 be the respective heights at the end of each slope. This is compiled for v1 as: and for v2 as:

Algorithm for the Integration of Dynamic Equation
According to the description in the previous sections, the algorithm for the integration of the dynamic equation in the modified form of Equation (18) is as follows: 1. Input data for catenary and pantograph. Regarding the catenary: the number of spans of the cantons, lengths of spans, moments of inertia of the elements of the Furrer beam section (body of the beam and the contact wire), characteristics of the materials, the geometry of the sections, characteristics of the flanges, etc. The characteristics of the equivalent section of the beam are then determined according to Equations (5) and (6). For pantographs, the following are entered: number and type of pantographs, masses, spring stiffnesses, dampers, etc. 2. The setting of the different terms of the dynamic equation at the initial instant: the stiffness matrix, damping matrix, mass matrix, and external load vector, according to Equations (2)-(4) and (7). The matrix of restriction conditions at the initial position of

Algorithm for the Integration of Dynamic Equation
According to the description in the previous sections, the algorithm for the integration of the dynamic equation in the modified form of Equation (18) is as follows: 1.
Input data for catenary and pantograph. Regarding the catenary: the number of spans of the cantons, lengths of spans, moments of inertia of the elements of the Furrer beam section (body of the beam and the contact wire), characteristics of the materials, the geometry of the sections, characteristics of the flanges, etc. The characteristics of the equivalent section of the beam are then determined according to Equations (5) and (6). For pantographs, the following are entered: number and type of pantographs, masses, spring stiffnesses, dampers, etc.

2.
The setting of the different terms of the dynamic equation at the initial instant: the stiffness matrix, damping matrix, mass matrix, and external load vector, according to Equations (2)-(4) and (7). The matrix of restriction conditions at the initial position of the pantograph(s) is also established according to Equations (8)-(15).

3.
The initial value of the integration variables is set to t 0 = 0: coordinates of the nodes q 0 and their derivatives, and the initial contact force λ 0 according to Equations (23) and (24). 4.
The integration cycle begins at the time of the differential Equation (18) in t n+1 . The cycle varies for different time values with n = 0,1,2 . . . , adding in each cycle the value of ∆t to t n : 1.
Calculation of the coordinates corresponding to the positions of the catenary nodes and the pantograph masses q n+1 , according to Equation (20).

2.
Calculation of the position of the contact elements of the pantograph(s): (y 3 ) n+1 and (y 4 ) n+1 according to Equation (10). If the contact element is deactivated, its position coincides with that of the pantograph head mass; in this case, Equation (14) or (15) shall be used.

5.
For the calculation of the forces, the positions of the contact elements, and the restriction conditions, the three cases of activation of the contact element are accounted for, according to the position of the pantograph along the line, as explained in Section 2.6. 6.
The remainder of the terms of the dynamic Equation (18) is also set in t n+1 : C n+1 , K n+1 , and R n+1 , and, in particular, the pantograph stiffness matrix must be modified if there is a take-off or a coupling according to Equation (22). 7.
Calculation of the velocity vector in the coordinates for the next integration cycle in t n+1/2 by the mid-pass correction of the velocity vector, using Equation (19).

Results
As a result of the previous study, a software tool called RICATI was developed. In this section, the results obtained from this software and the simulation are discussed.

Computational Aspects
To facilitate the use of the algorithms, the RICATI application was created, which incorporates the algorithms with a user-friendly interface for the introduction of values from a database. The tool also presents graphic results that allow conclusions to be obtained without having to directly interpret a table with hundreds of thousands of numerical data points [26].
In order to achieve an efficient software tool, the enhancements carried out in this work are related to the way of storing the matrices involved in the resolution of the dynamic equation, as well as the use of suitable methods for solving the associated static problem.
In particular, the sparsity and symmetry of the stiffness matrix are exploited, improving the efficiency of the implementation, dramatically reducing the memory storage requirements, and making use of iterative methods based on Krylov projection methods. The execution time spent for solving the static problem is also dramatically reduced. For this purpose, the SPARSKIT library by Yousef Saad is used [27].
Therefore, a high-performance implementation has to take into account the features of current architectures, for example, cache memory. These features are particularly important when rebuilding the traditional algorithms to a block-oriented implementation. Blockoriented algorithms reduce drastically the data flow between main memory and secondary memory enhancing the performance of the final implementation. These good features are obtained by using standard libraries in the implementations. In particular, for this paper, several computational kernel from BLAS library [28,29], especially BLAS 3, are identified and used.
RICATI executes the algorithms based on the data of two cantons with the number of spans indicated for each canton. As a result, an interactive graphic interface is provided that allows the behavior of the catenary to be observed before the passage of the pantograph.
RICATI is structured into four basic components, as shown in Figure 8: 1. Interface: The interface allows all of the data associated with the assembly to be entered. The set of parameters is stored in a database so that they can be retrieved, modified, or copied to obtain different calculations or perform simulations. The interface is structured in blocks according to the kind of data that they contain: a. General parameters: Generic data values are introduced regarding integration, speed, and data common to the whole catenary. b. Values for the catenary damping frequency: The critical damping percentage is established for two reference frequencies. c. Section data: Geometry of the overlapping span. d. Cross-section properties: Mechanical properties of the profile, wire, and flanges. e. Pantographs: Data of the pantograph to be used in the interaction. Up to four pantographs can be used. The pantographs are coded in such a way that they can be used in different calculations without the need to re-introduce all of their characteristics. f.
Canton spans: Number of spans in each canton and for each span, length and height of the support above the nominal value, and distance from the flange to the left support.
2. Interaction Algorithms: The dynamic library containing the algorithms for calculating the pantograph-catenary interaction. As a result, among other values, a table is generated indicating at each instant, according to the defined integration time step, and for each pantograph, the vertical position of two catenary points preset in the calculation, the height of the contact point, and the effort in each canton and plate. 3. Graph management: As a result, an interactive graphical interface is displayed in which different graphs can be obtained by independently selecting the abscissa (distance, time) and ordinate (efforts, elevation) axes, and allowing different graphs to overlap, in addition to the magnification of the area of the graphs marked by the user, as shown in Figure 9.

1.
Interface: The interface allows all of the data associated with the assembly to be entered. The set of parameters is stored in a database so that they can be retrieved, modified, or copied to obtain different calculations or perform simulations. The interface is structured in blocks according to the kind of data that they contain: a. General parameters: Generic data values are introduced regarding integration, speed, and data common to the whole catenary. b.
Values for the catenary damping frequency: The critical damping percentage is established for two reference frequencies. c.
Section data: Geometry of the overlapping span. d.
Cross-section properties: Mechanical properties of the profile, wire, and flanges. e.
Pantographs: Data of the pantograph to be used in the interaction. Up to four pantographs can be used. The pantographs are coded in such a way that they can be used in different calculations without the need to re-introduce all of their characteristics. f.
Canton spans: Number of spans in each canton and for each span, length and height of the support above the nominal value, and distance from the flange to the left support.

2.
Interaction Algorithms: The dynamic library containing the algorithms for calculating the pantograph-catenary interaction. As a result, among other values, a table is generated indicating at each instant, according to the defined integration time step, and for each pantograph, the vertical position of two catenary points preset in the calculation, the height of the contact point, and the effort in each canton and plate. 3.
Graph management: As a result, an interactive graphical interface is displayed in which different graphs can be obtained by independently selecting the abscissa (distance, time) and ordinate (efforts, elevation) axes, and allowing different graphs to overlap, in addition to the magnification of the area of the graphs marked by the user, as shown in Figure 9.

4.
Simulation: This allows for repeated execution of the calculation by modifying, at an established interval and increment, parameters such as the length of the overlapping span, the height of the slope, and the distance from the beginning of the slope to the end of the slope. The application repeatedly executes the algorithm by varying these data, thus allowing for a table with the maximum, minimum, deviation, number of take-offs, total distance of take-offs, and total time of take-offs in each independent canton or considering the sum of the efforts in the common zone of the two cantons.

Simulation:
This allows for repeated execution of the calculation by modifying, at an established interval and increment, parameters such as the length of the overlapping span, the height of the slope, and the distance from the beginning of the slope to the end of the slope. The application repeatedly executes the algorithm by varying these data, thus allowing for a table with the maximum, minimum, deviation, number of take-offs, total distance of take-offs, and total time of take-offs in each independent canton or considering the sum of the efforts in the common zone of the two cantons.

Simulation Results
The objective is to check the influence of the geometric characteristics of the slopes in the common area of the catenary before the passage of the pantograph. For the execution of the simulations, a set of two cantons with the following characteristics is defined as the basis for the calculations: • Two cantons of 20 spans of 10 m length in the normal spans. A pantograph of a simple plate and two stages defined according to EN50318 whose data are in Table 1 of Section 2.4.

•
As in standard EN50318 for the flexible catenary, no damping is assumed for the catenary.

•
The interface for the remainder of the values defining the catenary.
Using the RICATI application, we carried out different simulations in which the behavior of the catenary at the passage of the pantograph is evaluated according to values such as the mean, deviation of efforts, and take-offs in the common area.
A good result is considered to include a minimum deviation, no take-offs, and an average of the forces close to the pushing force of the pantograph (approx. 120 N).

Simulation Results
The objective is to check the influence of the geometric characteristics of the slopes in the common area of the catenary before the passage of the pantograph. For the execution of the simulations, a set of two cantons with the following characteristics is defined as the basis for the calculations: A pantograph of a simple plate and two stages defined according to EN50318 whose data are in Table 1 of Section 2.4. • As in standard EN50318 for the flexible catenary, no damping is assumed for the catenary.

•
The interface for the remainder of the values defining the catenary.
Using the RICATI application, we carried out different simulations in which the behavior of the catenary at the passage of the pantograph is evaluated according to values such as the mean, deviation of efforts, and take-offs in the common area.
A good result is considered to include a minimum deviation, no take-offs, and an average of the forces close to the pushing force of the pantograph (approx. 120 N).
To check the most suitable configuration, four simulations are carried out: • An overlapping span of 1 m, slope length 1 m, and variation in the height of the slope at its end.

Simulation 1
We consider an overlapping span of 1 m in both cantons beginning the slope of the slope from the beginning of the span. In each calculation, the height of the end of the slope varies from 0 m (fully horizontal) to 0.2 m with an increase of 0.005 m.
The analysis of the results is carried out in the common area and, as a result, we obtain the following graphs that represent the average values of the pantograph effort on the contact wire in Newtons, the standard deviation, and the average distance of take-offs. In this case, there are no take-offs in the common area of the catenary, so this graph is omitted.
The abscissa axis of the graphs represents the height of the slope at its end in meters, and the ordinate axis represents the effort (mean, maximum, and deviation) in Newtons, as obtained in the integration calculated with the given elevation.
By analyzing the graphs (Figure 10), it is observed that the deviation of the efforts in the common zone is greater because we increase the slope. Analogously, the average effort and the maximum values increase with a behavior similar to the variation of the deviation.
The effort of the pantograph passing through the common area has a better behavior in the case of zero elevation from the horizontal slope and, therefore, with a 0 slope before its installation. This may seem counterintuitive when considering that, due to the weight of the beam, it should have a certain negative slope and therefore be an obstacle for the passage of the pantograph. However, in the static configuration of the catenary assembly, the beam has a certain positive elevation due to the weight of the beam in the adjoining span, which would justify the results of this simulation.
We then observe that a rise added to the beam on the slope could eventually be counterproductive. However, this raised the question about the influence of the length of the sectioning span on these results. Thus, we proposed the second simulation.

Simulation 2
We consider a sectioning span of 1 m in both cantons, while maintaining the elevation of the sloping end at 0.1 m. In each calculation, the length of the slope length is varied from 0.5 to 1 m with an increase of 0.05 m.
The analysis of the results is carried out in the common area and, as a result, we obtain the graphs shown in Figure 11, which represent the average, maximum, and minimum values of the pantograph effort on the contact wire in Newtons, the standard deviation, and the average distance of take-offs. In this case, there are no take-offs in the common area of the catenary, so this graph is omitted.
The abscissa axis of the graphs represents the length of the slope in meters and the ordinate axis represents the effort (mean, maximum, minimum, and deviation) in Newtons obtained in the integration calculated with the given elevation.
Analysis of the graphs ( Figure 11) shows that the deviation of the stresses in the common area is adequate in cases in which the length of the slope is less than or equal to 0.45 m. The deviation of the stresses in the common area is adequate in cases in which the

Simulation 2
We consider a sectioning span of 1 m in both cantons, while maintaining the elevation of the sloping end at 0.1 m. In each calculation, the length of the slope length is varied from 0.5 to 1 m with an increase of 0.05 m.
The analysis of the results is carried out in the common area and, as a result, we obtain the graphs shown in Figure 11, which represent the average, maximum, and minimum values of the pantograph effort on the contact wire in Newtons, the standard deviation, and the average distance of take-offs. In this case, there are no take-offs in the common area of the catenary, so this graph is omitted.
The abscissa axis of the graphs represents the length of the slope in meters and the ordinate axis represents the effort (mean, maximum, minimum, and deviation) in Newtons obtained in the integration calculated with the given elevation.
Analysis of the graphs (Figure 11) shows that the deviation of the stresses in the common area is adequate in cases in which the length of the slope is less than or equal to 0.45 m. The deviation of the stresses in the common area is adequate in cases in which the length of the slope is less than or equal to 0.45 m. From this length, we see how the maximum and the deviation increase as the length of the slope increases until it coincides with the whole section span. Similarly, the minimum stresses decrease as the maximum increases from 0.45 m in slope length, confirming the increase in deviation.

Simulation 3
We consider an overlapping span of 0.5 to 5 m in both cantons with an increase of 0.1 m in each calculation. In all cases, a slope elevation of 0 m is considered, i.e., we use a horizontal beam without slope for overlapping span to obtain the optimum length of the beam to establish a suitable elevation. We continue with the same configuration of normal spans of 10 m.
As in the previous cases, the analysis of the data obtained by the calculations is carried out in the common area, and, as a result, we obtain the graphs shown in Figure 12, which represent the average and maximum values of the pantograph stress on the contact wires in Newtons, the standard deviation, and the total distance of take-offs. In this case, The effort at the passage of the pantograph through the common area has a better behavior in the case of small lengths of the slope (considering an elevation of 0.1 m from the end). This result is congruent with simulation 1 because the case of 0 m of the slope is equivalent to a horizontal overlapping span in which only the elevation obtained in the initial static configuration is applied. We then observe that the length of the overlapping span can significantly influence the results obtained for a small elevation of the beam and that it could be sufficient with the elevation obtained in the initial static configuration. In this case, the length of the overlapping span profile influences the elevation of the beam and, therefore, the results obtained.
To test this hypothesis, the following two simulations start from the 0 m elevation of the slope and the variation of the overlapping span length. In simulation 3, we consider a configuration equal to that used previously (i.e., normal spans of 10 m), and in simulation 4 we change the length of the normal spans to 5 m to check the influence of the length of the normal spans on the result.

Simulation 3
We consider an overlapping span of 0.5 to 5 m in both cantons with an increase of 0.1 m in each calculation. In all cases, a slope elevation of 0 m is considered, i.e., we use a horizontal beam without slope for overlapping span to obtain the optimum length of the beam to establish a suitable elevation. We continue with the same configuration of normal spans of 10 m.
As in the previous cases, the analysis of the data obtained by the calculations is carried out in the common area, and, as a result, we obtain the graphs shown in Figure 12, which represent the average and maximum values of the pantograph stress on the contact wires in Newtons, the standard deviation, and the total distance of take-offs. In this case, and as expected, take-offs may occur in the common area of the catenary; thus, we show the graph indicating the total length of take-offs produced in this area.
The abscissa axis of the graphs represents the length of the overlapping span in meters and the ordinate axis represents the effort (mean, maximum, minimum, and deviation) in Newtons obtained in the integration calculated with the given elevation.
Analysis of the graphs ( Figure 12) shows that the deviation of the stresses in the common zone is adequate in cases in which the span length is below 3.5-3.7 m, with slightly worse behavior at values below 1.5 m. Above 3.7 m, the behavior worsens significantly, and take-offs appear over longer distances. This behavior is explained by the absence of the section elevation effect compensated by the "excessive" length of the overlapping span, which causes it to lose the elevation and fall below the nominal height on the track. Thus, the pantograph "collides" with this beam with stresses of more than 2000 N.
In the configuration analyzed, the optimum distance is around 3.5 m. From this length onwards, it is necessary to elevate the slope to achieve optimum values in the common area of the catenary. However, if we consider the results separated by the canton represented in Figure 13, we observe that, from 3.1 m, take-offs occur in the wire of canton 1 and, from 3.4 m, in the wire of canton 2. To ensure that the optimum distance is below these values, it is necessary to maintain contact of the plate with the two wires simultaneously in the common area.
It is expected that the optimum distance of the overlapping span length is influenced by the length of the adjacent normal spans. We can check this situation by repeating this simulation with normal spans of 5 m instead of 10 m, as shown in the following simulation.

Simulation 4
We consider a sectioning span of 0.5 to 5 m in both cantons with an increase of 0.1 m in each calculation. In all cases, a slope elevation of 0 m is considered, i.e., we use a horizontal slope in the sectioning to obtain the optimum length of the beam and, thus, to ensure a suitable configuration. In this case, the normal spans of both cantons are considered to have a length of 5 m to check the differences with the previous simulation. As in the previous cases, the analysis of the data obtained by the calculations is carried out in the common area and, as a result, we obtain the graphs shown in Figure 14, which represent the average, maximum, and minimum values of the pantograph effort on the contact wire in Newtons, the standard deviation, and the average distance of take-offs. In this case, and as expected, take-offs may occur in the common area of the catenary; thus, we show the graph indicating the total length of take-offs produced in this area.

Simulation 4
We consider a sectioning span of 0.5 to 5 m in both cantons with an increase of 0.1 m in each calculation. In all cases, a slope elevation of 0 m is considered, i.e., we use a horizontal slope in the sectioning to obtain the optimum length of the beam and, thus, to ensure a suitable configuration. In this case, the normal spans of both cantons are considered to have a length of 5 m to check the differences with the previous simulation.
As in the previous cases, the analysis of the data obtained by the calculations is carried out in the common area and, as a result, we obtain the graphs shown in Figure 14, which represent the average, maximum, and minimum values of the pantograph effort on the contact wire in Newtons, the standard deviation, and the average distance of takeoffs. In this case, and as expected, take-offs may occur in the common area of the catenary; thus, we show the graph indicating the total length of take-offs produced in this area.

Discussion
The inclusion of the slope in an overlapping span may be necessary from a critical length onwards to achieve optimum values of force in the common area. However, in general, for the assumed speed, it is not necessary to introduce a slope height to obtain good results. This is because, in the initial static configuration, the section already presents a certain inclination as a consequence of the action of gravity. However, the length of the section should not exceed a certain value from which the results of the contact force sig- The abscissa axis of the graphs represents the length of the overlapping span in meters and the ordinate axis represents the effort (mean, maximum, minimum, and deviation) in Newtons obtained in the integration calculated with the given elevation.
An analysis of the graphs ( Figure 14) shows that the deviation of the stresses is optimal for the span length of approximately 1.4 m or less, with slightly worse behavior at values below 1 m. From approximately 2 m, the behavior begins to worsen significantly and take-offs appear over greater distances.
The behavior is similar to that of simulation 3, but because the normal spans are shorter, the effect of lifting the beam in the section is dissipated before the critical point at 1.4 m, from which the lifting of the beam is compensated by starting to fall below the nominal height, which causes problems with the passage of the pantograph.

Discussion
The inclusion of the slope in an overlapping span may be necessary from a critical length onwards to achieve optimum values of force in the common area. However, in general, for the assumed speed, it is not necessary to introduce a slope height to obtain good results. This is because, in the initial static configuration, the section already presents a certain inclination as a consequence of the action of gravity. However, the length of the section should not exceed a certain value from which the results of the contact force significantly fluctuate.
As the authors state at the beginning of this paper, the main contribution of this paper consists in a model for the pantograph-rigid catenary interaction considering an overlapping span. This model is a particular model for this kind of catenary and not an elastic catenary where the stiffness is increased.
In order to build this model, the authors start from well-known methods that appear in Finite Element Method literature, such as [22,23], but particularized for the elements of a rigid catenary and its underlying problems. In fact, according to the literature, the model uses an explicit method as an integration method.
The RICATI tool allows the study of dynamic behavior in which new profile designs are considered and compared with existing designs. Therefore, RICATI allows multiple simulations to determine the optimal catenary structure for any configuration. For example, with the same characteristics as the previous simulations, it is possible to vary the speed of the unit to observe the behavior at different speeds. We observed that, at speeds of 80 and 120 km/h, the behavior of the examples is similar. However, at 200 km/h, Figure 15 shows variations in the behavior and indicates two optimal points in the length of the sectioning around 1.2 and 3.9 m.
This example shows the usefulness of an application, such as RICATI, to check the behavior of the initially established configuration for a catenary. The tool also allows solutions to be obtained for the problems encountered when simulating the passage of the pantograph (or pantographs), not only for the section, but also for the entire catenary.
The proposed methodology provides a tool to simulate the pantograph-catenary behavior in a rigid catenary considering the transition between two cantons of the catenary with overlapping spans. To advance the presented method, the following lines of work are proposed: • Study of the behavior for different configurations in the geometry of the sectioning span. • Development of the model in three dimensions. • Study of the behavior and simulation in the transition between elastic and rigid catenary.

•
Optimization of the algorithm, taking into account the data structures underlying the problem, in addition to the parallelization of the calculations, to obtain results in significantly shorter times. • Development of RICATI following the paradigm of Software-as-a-Service. This model of the application in the cloud would allow access by more companies, from any point at which there is connectivity, and at a lower cost. In addition, this solution would allow different tests to be carried out simultaneously by requesting resources from the provider of the cloud service, and significantly more quickly, to find better solutions.
In addition, as an important future work, the authors will address the problem of the irregularities in the beam following similar ideas of those described in [30][31][32].
Appl. Sci. 2021, 11, x FOR PEER REVIEW 27 of 29 simulations to determine the optimal catenary structure for any configuration. For example, with the same characteristics as the previous simulations, it is possible to vary the speed of the unit to observe the behavior at different speeds. We observed that, at speeds of 80 and 120 km/h, the behavior of the examples is similar. However, at 200 km/h, Figure  15 shows variations in the behavior and indicates two optimal points in the length of the sectioning around 1.2 and 3.9 m. This example shows the usefulness of an application, such as RICATI, to check the behavior of the initially established configuration for a catenary. The tool also allows solutions to be obtained for the problems encountered when simulating the passage of the pantograph (or pantographs), not only for the section, but also for the entire catenary.
(a) Total distance of take-offs (b) Standard deviation of pantograph effort Figure 15. Graphs at total distance take-offs (a) and standard deviation (b) when the train is circulating at a speed of 200 km/h.
The proposed methodology provides a tool to simulate the pantograph-catenary behavior in a rigid catenary considering the transition between two cantons of the catenary with overlapping spans. To advance the presented method, the following lines of work are proposed: