Calculation of Rail Internal Impedance by Using Finite Elements Methods and Complex Magnetic Permeability

The internal parameter of UNI 60 rail is calculated by using finite elements methods. Steel’s characterizations by its normal magnetization curve and by complex magnetic permeability are here considered and included into the proposed FEM models. Rail’s resistance and internal inductance in function of current and frequency are calculated using both FEM and analytical models. The results obtained at the frequency of 50 Hz are compared with few measurements available, and then they are extended to other frequencies.


Introduction
The railway electric system is usually modeled by a Multiconductor Transmission Line (MTL) [1]. The characterization of this model by the definition of the resistance, inductance, capacitance, and conductance per unit length (p.u.l.) matrices requires that the electric and magnetic behavior of all the materials must be linear and isotropic. In the frequency domain, the MTL model is characterized by the impedance and admittance matrices defined starting from the previous ones. All the elements of these matrices are defined by the spatial position of conductors, their crosssection, and the magnetic and electric properties of the whole system's materials. Only for MTL systems with isotropic and linear materials and circular cross-section conductor widely separated, the elements of these matrices can be analytically computed [2].
Rails are one of the most important elements in the railway system; they are both guide of the railway vehicle and conductors of the electric system. In fact, rails are the common conductors between traction and signalling circuits and the interface between these circuits and the ground.
Their shape and the steel's properties are chosen taking only mechanical requirements into account; electric and magnetic characteristics of the rails are a consequence of these primary choices. As a result, rail cross-section is far from the circular shape and the steel's magnetic behaviour is nonlinear and hysteretic. The electric and magnetic characteristics are also neither complete nor certain; generally, they are derived from empirical formulas and they are based on the little data available in the literature. The computation of the tracks' representative parameters is, at the same time, a fundamental and an extremely critical task in the calculation of the total system's parameters. It is important to note that the simulation of the railway electric system requires a rail large signal model.
The main problems in the calculation of the parameters associated with rails are related to the rail's cross-sectional geometry and to the magnetic behaviour of the steel.
As a matter of fact, the complexity of the cross-section does not allow an analytic solution of the magnetic field and current density inside the rail, even in presence of material with linear magnetic and electric behaviour. The accurate calculation of the rail's internal parameters must be carried out using numerical techniques; the rail's representation with an equivalent conductor with a circular cross-section must be considered the simplest, but the heaviest, approximation.
Furthermore, the nonlinear behaviour and hysteretic phenomena in the steel make the internal parameters of the rail dependent on both the frequency and the current. In order to include the rails into an MTL large signal model it International Journal of Vehicular Technology is necessary to introduce a proper linearization in the steel's magnetic behaviour.
Some papers try to address the problems underlines. In particular, in [3] FEM simulation results and measurements on a BS110A English rail type are reported, both for large and small signal. In [3], the rail steel hysteresis cycles are reported and the complex magnetic permeability is introduced, but it is not applied in the FEM proposed model. In [4] the measured values of the internal parameters of a UNI 60 rail for large signal at 50 Hz are reported and in [5] these values are recalled and collected together with other values obtained by many authors on different rail types and at different power frequencies. In [1,6] the small signal measurements on a UIC 60 rail are presented, in particular the behaviour of a UIC 60 rail, biased with a large AC 50 Hz current and a DC current, respectively, is investigated in the audiofrequency range.
In this paper, the methods proposed in the literature for the representation of the steel's magnetic permeability are introduced and successively applied to the calculation of the internal parameters for large signal simulation, with reference to a UNI 60 rail, by using FEM. In particular, two different FEM rail's models are proposed. These models are based on two different kinds of linearization of the magnetic material behaviour. The first one takes into account the normal magnetization curve of the material and adopt the linearization method based on the sinusoid equivalent concept. The second rail's FEM model includes both the effects of hysteresis and saturation by using the linearization based on complex magnetic permeability function.
The obtained results are also compared with experimentally measured data [4,5] and other methods introduced in the literature.

Calculation of the Magnetic Permeability
The internal impedance of a conductor is related to the distribution of the magnetic field and current density on the conductor's cross-section; these distributions depend both on the shape of the cross-section and on its electric and magnetic characteristics.
For a typical geometry with a semi-infinite slab of isotropic and linear conductive material in which sinusoidal current flows, the modulus of the magnetic field and current density fall down exponentially with the increase in the distance y from the surface of the conductor where The terms σ and μ represent the conductivity and magnetic permeability of the material, respectively, and f is the current frequency. The term α represents the attenuation constant and β corresponds to the phase constant. At a distance from the surface of the slab equal to the skin depth δ,  the magnitude of the fields is reduced by a factor of 1/e with respect to the magnitude of the same fields on the surface of the slab.
In a ferromagnetic material, both saturation and hysteresis are observed in the magnetic behaviour. These phenomena have been measured in a sample of rail steel [3,7,8]. Under these conditions, the skin depth as defined in (2) loses its significance. Although saturation and hysteresis phenomena modify the fields distribution into the conductor, they do not affect the qualitative aspects of the skin effect.
In the next part of this chapter, some techniques for deriving an approximated value of the magnetic permeability, starting from the hysteresis loop and/or the normal magnetization curve, will be introduced. These techniques allow the deduction of a permeability value that takes into account the quantitative aspects related to presence of saturation and/or hysteresis in the material's magnetic behaviour.

Evaluation of the Saturation Effects.
The effect of saturation on the magnetic field and current density distribution can be evaluated with reference to the typical geometry (semi-infinite slab) previously described. In this type of analysis the material magnetic behaviour is represented with an appropriate function that relates the induction B to the magnetic field H.
Steel is both nonlinear and hysteretic; in order to take into account only the saturation effects the normal magnetization curve is introduced into this model. Numerical results are here obtained taking into account this curve measured on a rail's steel sample in [9].
The distribution of magnetic field in the semi-infinite slab conductor previously introduced can be computed by solving the following differential equation: where International Journal of Vehicular Technology The numerical solution of (3) can be computed using the finite difference method (FDM). The associate difference equation is where The terms p and h are the time and space integration steps, respectively. The stability of this method is guaranteed for values of u less than 0.5. Saturation modifies the attenuation of the magnetic field inside the conductor with respect to the linear case and it introduces distortion in the waveform. The reduction of the magnetic field's maximum amplitude with the increase in distance from the surface is, however, very close to the linear case exponential decay. Figure 1 shows the magnetic field distribution into the semi-infinite slab calculated for a time period, when the initial transient has been extinguished.
The rapidity with which the amplitude of the magnetic field falls down depends both on the frequency and the maximum value of the magnetic field on the surface of the conductor. The effective skin depth can be defined as the distance from the surface where the magnetic field maximum value is reduced by a factor of 1/e with respect to the maximum value of the same field on the surface of the slab. From relation (2), it is immediate to define the function of effective permeability as It is important to underline that effective permeability does not depend on the frequency; it depends only on the magnetic field strength on the surface of the slab. This is true only if the non-linear BH function that represents the magnetic characteristic does not change with frequency. Figure 2 shows the effective permeability and the effective skin depth evaluated at 50 Hz as a function of surface magnetic field strength.

Evaluation of the Hysteresis Effects.
It is known that the hysteresis phenomenon causes time delay between the induction and the magnetic field, and that the loop area represents the energy converted in heat in a cycle. The shape of the hysteresis loop for a linear hysteretic material is an ellipse, in fact induction and magnetic field pulse at the same frequency and the angle between them can be represented, in the phasor domain, considering a complex value of magnetic permeability The component of the induction in phase with the magnetic field is associated with the average energy stored in magnetic field; therefore, the in-quadrature component is related to hysteresis power loss.
Referring to the semi-infinite slab previously described, the skin depth takes the following complex value For small arguments of the complex skin depth, the attenuation constant depends on an equivalent skin depth defined as The phase constant is given by The angle of complex permeability defined in (8) is always negative, so (10) and (11) show that the attenuation constant increases in the presence of hysteresis with respect to the case of material with real magnetic permeability equal to the modulus of complex permeability and the phase constant decrease. If the hysteresis losses are small, the angle of complex permeability is also small. In this case, it is possible to neglect the reduction of the skin depth introduced by hysteresis.

Saturation with Hysteresis.
In the case of the simultaneous presence of saturation and hysteresis, the behaviour of the material can be approximated using an opportune value of complex magnetic permeability evaluated as follows [3]. The angle of complex permeability can be determined from the energy analysis of the fields B and H. For a sinusoidal magnetic field, the induction is given by The lost power density due to hysteresis is given by Using a typical formulation of the analysis of the circuits under nonsinusoidal conditions, the apparent power density per unit volume is given by [10] S The angle of complex permeability is defined as The angle of the complex permeability depends on the maximum value of the magnetic field, which is related to the shape of the hysteresis cycle. The angle of complex permeability can be considered independent of the frequency until the shape of the hysteresis cycle does not change with frequency. The case of small hysteresis losses corresponds with little values of the complex magnetic permeability angle.
The evaluation of the modulus of complex magnetic permeability depends on the distribution of the fields into the material. In the condition of small hysteresis losses the attenuation is reasonably similar to the only saturated case described in Section 2.1, so that the modulus of complex permeability can be taken equal to the value of the effective magnetic permeability. This is a practical approximation; the numerical solution of the magnetic field in the presence of small hysteresis losses [3] confirms that the increment of attenuation is very little and than should be neglected.
The components of complex magnetic permeability defined in (7) and (15) are related to the magnetic field maximum value on the conductor's surface. In order to relate complex magnetic permeability to the current that flows into a conductor it is necessary to take into account a conductor with a cross-section that will produce an uniform magnetic field on its surface, like the semi-infinite slab surface. The simplest choice is the circular conductor, and Ampere's law for a cylindrical conductor must be taken into account: where r represents the radius of the conductor. The function of complex magnetic permeability becomes For conductors with cross-section different from the circular shape it is possible to apply relation (16) using proper values for r, as described in the next chapter. It is however an approximation because the circular equivalent conductor neglects the effect of corners that are related to the nonuniformity of the magnetic field on the conductor's surface.

Circular Equivalent Conductor Model
The substitution of the rail with a circular equivalent conductor avoids the problem of calculating the fields into the cross-sectional surface, that is related to the complexity of the geometry, and it requires the setting of a proper rail equivalent radius. The internal impedance of a circular crosssection conductor is defined by the following relation: In order to take into account both the effects of saturation and hysteresis losses into this model, the complex magnetic permeability is considered; its introduction in (18) make the internal parameters dependent on the current. Equation (18) modifies as The complex magnetic permeability function is related to the current by relation (16), in which the parameter r represents the equivalent radius chosen as follow.
The equivalent radius can be defined by imposing the equality, at a given frequency, between the actual value of one of the two internal impedance components and the same value calculated by the equivalent circular model. Generally, DC and AC high frequency resistance is the actual values used to define, respectively, the DC and AC high frequency equivalent radius. The calibration of this model is very simple and it do not require measurement or simulations on the actual conductor, but the regulation can be done only at one frequency and it is not possible to regulate both the resistance and internal inductance simultaneously. As a consequence, it is impossible to fit the real frequency dependant internal impedance components in a wide frequency range and model's results are accurate only for one component and in a small frequency range.
In order to achieve a more precise analytical model, some authors propose a 4 parameters tubular equivalent conductor [11,12]. This model fit the actual frequency-dependent parameters better than the circular equivalent conductor and International Journal of Vehicular Technology 5 its accuracy is extended in a wider frequency range. However, the calibration of this model requires information about the real parameters that are available only from measurements and/or simulations [12].
For DC calculation, two different values for the equivalent radius can be taken into account. The first one is the radius that gives the same rail area and it produces the same DC rail resistance For a UIC 60 rail, equivalent area radius is equal to 49.4 mm.
The second one is based on the geometric mean radius. The geometric mean distance (GMD) and geometric mean radius (GMR) of an irregular conductor are needed for inductance calculation under steady-state condition [13]. The use of GMR as rail equivalent radius will produce the same DC rail self-inductance. For an irregular shape like the rail cross-section, GMR is difficult to calculate because the two double integral [14]; for the UIC 60 rail this value is numerically evaluated and it is equal to 62.9 mm. Equivalent area radius and GMR methods are developed under DC conditions, but they can also be applied for the very low frequency AC calculation if the current density distribution over the cross-section is almost uniform. This is the classic assumption adopted in the calculus of the electric power lines parameters at power frequency.
As frequency goes up they the nonuniform distribution of the current density caused by the skin effect becomes relevant and these methods cannot be applied. Due to the presence of ferromagnetic material and the dimension of the rail, the upper limit of the very low frequency band is about a few Hertz. It is important to note that even at 1 Hz, a 100 A current produces a straightforward non-uniformity in the current density.
In AC high frequency the skin effect forces the current in the rail periphery. In this case a equivalent perimeter radius must be taken into account: For a UIC 60 rail, equivalent perimeter radius is equal to 107.7 mm. The lower limit of the high frequency band is represented by the minimum frequency that will produce a distribution of current that is almost peripheral and in which local effects at the corners can be neglected. This requires that the skin depth will be, at least, less than 1 mm, that corresponds to a frequency of about 1 kHz. In the transition range, included between the low and high frequency range boundaries, an intermediate equivalent radius must be taken into account and the area and perimeter equivalent radius represent the extreme values for the rail equivalent radius.
It is important to note that transition range include power frequency and the lower part of the harmonic spectra produced by traction converters. The definition of this frequency boundaries is however difficult, and considerations about the conductor's shape and dimensions, and the skin  depth will be helpful for this task. Finally, these limits depend on the current due to the complex magnetic permeability function. However it is possible to consider about three decades between these two limits.

FEM Models
The calculation of the internal parameters of a rail can be carried out also using finite element softwares. Such tools resolve the magneto-electric problem in a numerical way starting from the geometry of the problem, the characteristics of the materials, the sources, and the boundary conditions. Thanks to the uniform cross-section of the rail, the model can be performed in 2D on the cross-section surface, reducing consequently the simulations computational time without lack of precision. The model has been implemented using Ansoft Maxwell version 10. A rail type UNI 60 is represented in the vacuum in order to neglect any kind of proximity effects; this choice matches the conditions under which the circular equivalent conductor model is derived. Moreover, the measurements conditions in [3,7] are very close to the simulated case. The simulations are performed using "Eddy Current Field Solver." This solver assumes that all currents are sinusoids oscillating at the same frequency, and they produce a time-varying magnetic field in the plane perpendicular to the conductors in which currents flow. The magnetic field induces eddy currents in the source conductor and in any other parallel conductors. Eddy Current Field Solver uses FEM to compute the magnetic vector potential and the electric scalar potential in the phasor domain using the following relations: In a 2D problem, induction B lies in the xy plane and the electric field E has only z-component, so that the magnetic potential in (22) has only z-component and the electric potential is constant over the cross-section. Fields are computed starting from the magnetic vector potential numerically obtained.
The simulations are performed imposing the current on the rail's surface and taking into account that the empty space extends infinitely out, placing "balloon" condition on the boundary of the model.
The average internal energy and the average losses per unit length are calculated by integrating the average values of the specific correspondent quantities over the surface of the conductor obtained by FEM simulations. Under sinusoidal conditions, these quantities are given by The internal impedance components per unit length are determined from a port-process operation that consider the magnetic stored energy (23) and power losses (24). Under sinusoidal conditions, the resistance and the inductance per unit length are given by The FEM overcomes one of the main problems in the calculation of the parameters: the one associated to the rail's cross-sectional geometry. The problem linked to the magnetic behaviour of the steel can be faced developing two different FEM models based on two different kinds of linearization of the magnetic material behaviour.
The first one takes into account the first magnetization curve of the material and adopts the linearization method based on the sinusoid equivalent concept. Eddy current nonlinear solver calculates the effective permeability of nonlinear materials as the ratio of the effective induction B e and effective magnetic field H e of the original BH curve: These quantities are defined starting from energy considerations. In fact, assuming a sinusoidal magnetic field with amplitude H m , B e is defined as the amplitude of the sinusoidal induction that gives an average energy density equal to the one obtained considering the real induction, determined by the normal magnetization diagram of the material: The effective induction and effective magnetic field can be simple derived from (28) and (29) as follow: It is important to underline that the values of B m and H m are automatically computed during simulation by the nonlinear eddy current solver, and they are related to the electromagnetic quantities values in the problem. In this way it is possible to calculate internal parameters taking into account only the effects of saturation in an approximated way. Figures 4 and 5 show the resistance and internal inductance calculated by using FEM nonlinear solver. The second rail's FEM model includes both the effects of hysteresis and saturation by using the linearization based on complex magnetic permeability function.
In fact, as in this software it is not possible to define the hysteresis loop in the material properties, the magnetic losses due to a time varying field and the effects of saturation can be included in the model only using a proper value of complex magnetic permeability for each value of current.
The complex magnetic permeability function is calculated by (17) starting from the current flowing into the rail. This assumption implies however an approximation since relation (17) is derived starting from the hypothesis of circular conductor cross-section, that implies the uniformity of the magnetic field on the surface conductor. This condition, as it is noticeable in Figure 6, is far from the actual rail working conditions. The use of the function of complex magnetic permeability represents a very efficient method to calculate the rail parameters using the finite element software and taking into account, in approximated way, both saturation and hysteresis   losses. Figures 8 and 9 show the resistance and internal inductance calculated by using FEM complex permeability. Figures 6 and 7 show some examples of the solution of the magnetoelectric problem that is necessary to deduce the internal parameter using relation from (23) to (26). It is important to underline that both FEM models require the solution of the magneto-electric problem. For this reason, this kind of simulation needs more computation time than the equivalent circular conductor model. Finally, it is important to note that computation time of FEM models increases with the requested precision and with frequency because the number of triangles in the region near the rail's perimeter increases with the skin depth reduction.

Discussion of Results
The graphs that show the dependence of internal parameters on the effective current and the frequency obtained from the two FEM proposed models are quite similar in shape, but the difference in values is remarkable.
The comparison between models' results and measured internal parameters of a UNI 60 rail type can be carried out only at the frequency of 50 Hz [4,5] as reported in Figures  10 and 11, where both measurements and their interpolation are plotted.
The resistance calculated with the nonlinear FEM model is less than the resistance calculated with the linear FEM model with ferromagnetic material characterized by complex permeability; the internal inductance calculated with the nonlinear FEM model is bigger than the internal inductance calculated with the FEM model based on complex magnetic permeability. Such differences are associated with the hysteresis losses that are neglected in the nonlinear FEM model based on the normal magnetization curve. In fact the effects of hysteresis are to increase the AC resistance and decrease the internal inductance by reducing the average energy stored in the material [3].  Figure 7: Current density inside a rail and flux lines computed with a current of 100 A rms and 1000 A rms for the typical railway power frequencies. Resistance and internal inductance increases with current. The magnetic permeability increase with current from 0 A to about 800 A, and only above this current saturation will be noticeable.
In order to achieve an expression between internal parameters and the current at the frequency of 50 Hz, an interpolation of measured data has been computed. The resistance and internal inductance have been interpolated using third-and fourth-order polynomials, respectively [4],   It is important to underline that the measurement, of hysteresis cycles that characterize UNI 60 rail's steel and its resistivity are not available. These measurements are reported only for an English rail type BS110A, that is very similar to the UNI 60 rail. The measured hysteresis cycles in [3,8] [4]. In this work a mean value of normal magnetization curve and hysteresis loops in [3,8] is taken into account and a resistivity of 2.2510 −7 Ω · m· is considered. For all these reasons, a complete validation of the models is subordinated to the experimental measurements on a steel sample from UNI 60 rail. The FEM model based on the normal magnetization curve is more sensitive to the accuracy of the parameters than the model based on the complex magnetic permeability. Small changes in the normal magnetization curve, especially in its initial values, will produce large changes in the internal parameter computed using this model. Small changes in the shape of hysteresis loops will produce small changes in the complex magnetic permeability function.
Although the characterization of the rail steel is not available with reasonable precision, the results obtained with the FEM models are very close to the measured values. The model that gives better results, both in resistance and internal inductance, is the FEM one based on the complex magnetic permeability. The little difference between the internal parameters curves calculated using this model and the measured data interpolation one can be associated to the uncertainty both in the parameters used in the models and in the measured data. Internal inductance (mH/km)

Rail current (A rms )
Measured [4], [5] Interpolated [4], [5] FEM nonlinear FEM compl. perm. Figure 11: Internal inductance versus current at 50 Hz. In Table 1 a quantitative analysis of the differences between measured and calculated internal parameters is reported. These differences are calculated as follows: where R calc and L calc represent the values of the internal parameters calculated with the FEM models; R mis and L mis are the interpolated values calculated from (26) and (27).

Conclusion
Two methods based on the FEM software for the calculation of the internal parameters of rail with reference to the rail type UNI 60 have been introduced. The obtained results have been compared to data available in the literature. This comparison shows that the accurate evaluation of the internal parameter of a rail can be performed taking into account the rail geometry and hysteresis effect using FEM software.
It is important to emphasize that the results are related to the precision with which the magnetic and electric characteristics of the rail steel are known [9], with particular focus on its magnetization characteristic, especially in the initial low values of the magnetic field. The FEM model based on the linearization of the normal magnetization curve is highly influenced by the shape of the steel's B-H characteristic. Small variations of this curve, especially in the very low field region, imply high variations of the resistance and internal inductance.
The results here obtained depend on the shape of the hysteresis cycle. In order to achieve good accuracy, it is necessary to know exactly the hysteresis cycles' family for various values of excitation and frequencies.