A principle of virtual power-based beam model reveals discontinuities in elastic support as potential sources of stress peaks in tramway rails

This paper addresses the question to which extent the stiffness of an elastic foundation in general, and of a discontinuity occurring in the foundation in particular, affects the stresses to which a beam resting on such a foundation is exposed to. This is particularly relevant for tramway rails, where foundation discontinuities may be generated in the course of maintenance works. In the presented work, the underlying mathematical framework was derived based on the principle of virtual power, using a suitable virtual velocity field and a foundation-related traction force vector as input quantities. The resulting mathematical expressions for the virtual powers performed by the internal and external forces were approximated based on Finite Element discretizations, and eventually solved correspondingly, in the format of sequential Finite Element analyses. Numerical studies, performed on a tramway rail-shaped beam, have confirmed that foundation discontinuities indeed induce substantial increases in the stress tensor components, when moving from regions of a high-stiffness foundation to regions of a low-stiffness foundation; as it may occur if specific sections of the tramway network get renewed.


Introduction
Tramway networks usually comprise different variants of track superstructures, depending on the particular requirements concerning noise and vibration mitigation; see, e.g., [1] for a respective overview which is particularly relevant for the Viennese tramway network. While those variants differ from each other in terms of the types and thicknesses of the layers building up together the track superstructure, all of them have in common an elastic layer on which the tramway rails rest. Together with the subjacent structure, this elastic layer constitutes an elastic foundation supporting the tramway rail. Unarguably, the compliance of this layered foundation influences the mechanical behavior of the rails decisively. As a matter of fact, maintenance and repair works performed on tramway rail tracks sometimes involve the replacement of whole rail sections. Such activities are likely to cause discontinuities in the elastic foundation, given that the involved materials degrade and, even more, may have partly changed entirely over the past decades. Investigating the potential effects such discontinuities may have on the elastic deformations and stresses tramway rails experience in response to practically relevant mechanical loading was the main objective of the work presented in this paper.
Understanding and predicting the effects of elastic foundations on thereon resting structures, in terms of displacements and stresses arising in response to mechanical loading, has been a long-pursued goal of the theoretical and applied mechanics community. Many consider Winkler as founding father of the field; he assumed that the deflection of an elastically founded structure is linearly related to the forces acting upon this structure via a foundation modulus [2]. In mathematical terms, the classical Winkler foundation reads as p = k · w, where p and w denote the pressure and the vertical displacement distributed across the foundation surface, and k denotes the material-specific foundation modulus [3]. Later, aiming for a more realistic consideration of the interactions between mechanically loaded structures and elastic foundations, the Winkler foundation was extended by additional terms, resulting in much-noticed works published by some of the most famous mechanicians active in the twentieth century, including Biot [4], Filonenko-Borodič [5], Hetényi [6], Vlasov [7], Reissner [8], or Pasternak [9]; see also the article by Kerr [3] for a comprehensive review on the subject. More recent approaches focused on refined model assumptions and consideration of additional characteristics and loading modes, including, e.g., improved consideration of boundary conditions [10], viscoelastic foundations [11], dynamic effects of moving loads [12], and torsional vibration and buckling [13][14][15]. Due to the ever-growing relevance of computational mechanics from the 1970s onwards, numerical solution strategies have found their way into the field of elastic foundations [16,17], allowing for solution of the governing mathematical framework while investigating (more) complicated structural settings. As concerns rails, several studies have been published dealing with Vignole rails [18][19][20][21][22][23][24][25], whereas models related to tramway operation appear to be limited to investigations concerning the ground-borne noise and vibrations [26,27]. Nevertheless, despite the considerable degree of sophistication that has been reached by now, it seems that the mechanism of restrained torsion, which affects tramway rails in a potentially substantial manner, has been disregarded so far.
In this work, we took a new path towards consideration of elastically supported beams in general, and of tramway rails in particular, tying in with the approach recently published by the present paper's authors [28], where virtual kinematics were introduced reflecting stretching, bending, as well as primary and secondary warping deformations of a beam, and through rigorous application of the principle of virtual power [29] respective equilibrium equations and natural boundary conditions could be derived, including higher-order torsion-related force quantities. Eventually, this mathematical framework was utilized for computing displacement quantities along the longitudinal axis of the beam, and corresponding stress distributions in selected cross sections, both by means of the Finite Element method, see also [30] for further details. In order to extend this approach, elastic foundation-related traction force vectors are introduced, leading to respective extensions in the model-governing equations, see Sect. 2. Their numerical solution, based on sequentially performed oneand two-dimensional Finite Element analyses, is presented in Sect. 3. Based on numerical examples, presented in Sect. 4, the potentials and features of the new modeling approach presented in this paper are discussed in Sect. 5, containing also concluding remarks.

A brief introduction to the principle of virtual power
The principle of virtual power (PVP) [29,[31][32][33][34] was chosen to be the fundamental theoretical vehicle of the subsequently presented work. In general format, the PVP reads as where P int is the virtual power performed by internal forces, and P ext is the virtual power performed by external forces. Considering a three-dimensional (3D) continuum, P int and P ext are defined as follows: and In Eqs. (2) and (3), x denotes the location vector of points throughout the continuum with volume V and surface S,v(x) denotes an arbitrary continuous virtual velocity field,d(x) = ∇ Sv (x) denotes the virtual Eulerian strain rate tensor, f(x) and T(x, n) denote the volume and surface forces, and σ (x) denotes the Cauchy stress tensor. Thereby, the PVP is only confined in terms of the assumption thatv(x) is related to geometrically compatible virtual displacements viav(x) =u(x) = ∂û(x)/∂t, t being the time variable.

Model representation of tramway rails
First, we stipulate that, on the one hand, the length of the beam under investigation (which is considered to be straight and of prismatic shape) is much larger than its cross-sectional dimensions, and that on the other hand, cross sections which are plain and oriented perpendicular to the beam axis in the (undeformed) reference configuration, remain, also after the application of mechanical loads, plain and oriented perpendicular to the (then deformed) beam axis-or, in other words, that the classical assumptions of a Bernoulli beam are fulfilled [35][36][37]. In order to apply the PVP, as introduced in Sect. 2.1, a virtual velocity fieldv(x), which represents this beam reasonably well under the expected loading conditions, needs to be defined. This field is defined through componentsv x (x),v y (x), andv z (x), with respect to an orthonormal base frame spanned by vectors e x , e y , and e z , hence x = x e x + y e y + z e z , whereby x is the coordinate in direction of the beam axis, and coordinates y and z are measured from the geometrical center (abbreviated as GC) of the cross section; hencê v(x) =v x (x) e x +v y (x) e y +v z (x) e z . Figure 1 shows a schematic illustration of such a beam.
In line with Sapountzakis [38], who presented a displacement field representative for torsion-compliant beams, we consider the virtual velocity field introduced in [28 andv wherev GC represent the virtual velocities of the geometrical center of the cross section (at the longitudinal coordinate x), whileω x (x) stands for the angular virtual velocity with respect to the twist axis, which is the axis parallel to the beam axis passing through the center of twist (abbreviated as CT), being offset from the geometrical center by y CT and z CT ; determination of y CT and z CT has been demonstrated in [28]. Furthermore, ψ I (y, z) denotes the primary warping function describing the contribution of non-restrained, Saint-Venant-type torsion, whereas ψ II (y, z) denotes the secondary warping function, describing the effect of Fig. 1 External forces acting on a prismatic beam undergoing stretching, bending, shear, and torsional deformations, namely 3D continuum-specific traction vectors acting onto the lateral surface, T(x, n = n C ) = i T i (x, n = n C ) · e i , traction vectors acting onto the cross-sectional surfaces, T(x, n = e x ) = i T i (x, n = e x ) · e i and T(x, n = −e x ) = i T i (x, n = −e x ) · e i , volume force vectors acting throughout the entire volume, f(x) = i f i (x) · e i , as well as elastic foundation-related traction force vectors acting onto the bottom surface, T EF (x, n = n C bottom ) = T EF z (x, n = n C bottom ) · e z . Note that the depicted beam is, for simplicity, of rectangular cross-sectional shape (implying GC = CT), while the model is not limited to any particular kind of cross-sectional shape warping torsion. The related components of the virtual Eulerian strain rate tensor follow straightforwardly from the virtual velocity field defined by Eqs. (4)-(6), reading as follows: andd Finally, the internal and external forces taken into account in Eqs. (2) and (3) must be defined. If the simplifying assumptions related to a Bernoulli beam are sufficiently fulfilled, the stress tensor components σ x x (x), σ xy (x), and σ xz (x) are considerably larger than the remaining components, i.e. σ yy (x), σ zz (x), and σ yz (x), implying In Eq. (10), symbol ⊗ denotes the dyadic product. The volume force vector and the general traction force vector, both occurring in Eq. (3), are defined as follows: and ∀ x on S: T(x, n) = T x (x, n) · e x + T y (x, n) · e y + T z (x, n) · e z , (12) where f i (x) and T i (x, n), i = x, y, z, denote the components of f(x) and T(x, n). Importantly, the definition of the traction force vector includes the unit normal vector at location x on the surface. For example, n = n C indicates the lateral surface of the beam, n = −e x its starting cross section S 0 at coordinate x = x 0 , and n = e x its end cross section S 1 at coordinate x = x 1 , see Fig. 1.
Furthermore, as key novelty of this contribution, the effect of an elastic foundation is considered as well. To that end, one additional, elastic foundation-related traction force vector is taken into account, In Eq. (13), S bottom denotes the bottom surface of the beam on which it rests, k EF denotes the foundation modulus, u GC z (x) denotes the actual displacement of the geometrical center, and φ x (x) denotes the actual rotation angle with respect to the axis that is parallel to the beam axis and passes through the center of twist. Importantly, depending on how the beam is mounted to the subgrade, an elastic foundation may represent a tensional or tensionless interaction between the structure and the subgrade. The latter case is standardly considered through setting the foundation modulus to zero at all locations where the applied mechanical loading leads to a lift-off of the beam [21,39,40]. Owing to the fact that the beam studied in this paper may undergo not only flexural but also torsional deformations, a different approach is pursued in order to take into account a tensionless elastic foundation. Namely, only the sections of the bottom surface of the studied beam (which is actually in contact with the subjacent elastic foundation) are considered for calculation of T EF z (x), on which the following condition for the overall displacement in z-direction is fulfilled: More details on the implementation of this condition are provided in Sect. 2.3 of this paper.

Virtual powers performed by internal and external forces
Considering the virtual Eulerian strain rate tensor according to Eqs. (7)-(9), as well as the stress tensor according to Eq. (10) in Eq. (2) yields The terms on the right-hand side of Eq. (14) include a number of stress resultants provoked by the initially defined virtual beam kinematics, namely the axial force, N (x), where A is the cross-sectional area; the bending moments with respect to the y-axis, M y (x), and with respect to the z-axis, M z (x), with j ∈ {y, z}, as well as with β j = z if j = y and β j = −y if j = z; the classical torsional moment M T (x), and the warping-related fourth-order moment M ψ Considering, on the other hand, the virtual velocity field according to Eqs. (4)-(6), the volume and traction force vectors according to Eqs. (11) and (12), and the elastic foundation-related traction force according to

Eq. (13) in Eq. (3) yields
and with y 1 (x) and y 2 (x) as coordinates indicating the tensionless section of the beam's bottom surface at longitudinal coordinate x, see Fig. 2. Equation (22) includes a number of external forces, which are defined as distributed forces along the beam axis, or as forces acting onto the starting and end cross section of the beam, namely the distributed axial load per length n(x), with cross-sectional contour C; the axial forces N (x i ) acting onto the starting and end cross section of the beam, with i ∈ {0, 1}, as well as with β i = −1 if i = 0 and β i = 1 if i = 1, implying, when taking into account the traction-reaction law T(x, n) = −T(x, −n), Possible cases for the integration boundaries relating to calculation of T EF z (x), comprising a no lift-off, b partial lift-off, and c complete lift-off of the bottom surface, leading to corresponding choices of y 1 and y 2 , compare Eq. (23); the coordinates y 1 and y 2 , as well as L bottom , which is the contour of the bottom surface S bottom along the y-direction, are functions of the beam axis x. Note that the depicted cross section represents profile Ri60R1 [41] the distributed transverse (shear) loads per length in y-and z-direction, s y (x) and s z (x), with j ∈ {y, z}; the shear loads in y-and z-direction, S y (x i ) and S z (x i ), acting onto the starting and end cross section of the beam, with i ∈ {0, 1} and j ∈ {y, z}; the e y -and e z -oriented distributed bending moments per length, m y (x) and m z (x), with j ∈ {y, z}, whereas β j = z if j = y and β j = −y if j = z; the e y -and e z -oriented bending loads acting onto the starting and end cross section of the beam, M y ( with i ∈ {0, 1} and j ∈ {y, z}; the distributed torsional moment per length m T (x), the "classical" torsional moments M T (x i ) acting onto the starting and end cross section of the beam, the bimomental loads M ψ II (x i ) acting onto the starting and end cross section of the beam, with i ∈ {0, 1}; the distributed fourth-order warping moment per length, m and fourth-order moment loads M ψ IV (x i ) acting onto the starting and end cross section of the beam, with i ∈ {0, 1}.

Equilibrium equations and natural boundary conditions
Next, both Eqs. (14) and (22) are considered, and partial integration is applied to each integral term until the virtual velocity components, namelyv GC , do not occur anymore in differential form within the respective integral term. Inserting then the such adapted formulations of P int and P ext into Eq. (1) yields a set of equilibrium equations, as well as a set of natural boundary conditions. The equilibrium equations read as and whereas the natural boundary conditions read as and with i ∈ {0, 1}.

Consideration of actual beam deformations and linear elastic constitutive behavior
The nonzero components of the linearized strain tensor ε(x) corresponding to the virtual Eulerian strain rate components defined in Eqs. (7)-(9) read as and being hence defined through the components of the actual displacement vector related to the geometrical center, u GC x (x), u GC y (x), and u GC z (x), as well as through φ x (x), which is (as defined previously, in Sect. 2.2) the actual rotation angle around the center of twist. Inserting Eqs. (45)-(47) into Hooke's law specified for the stress state defined in Eq. (10) yields and where E is the Young's modulus of the material the beam is composed of, and G is its shear modulus. Furthermore, σ bend xy and σ bend xz are the shear stresses due to bending. The latter are not directly linked to beam axis displacements or cross-sectional rotations, but result from a boundary value problem based on the equilibrium equations for a classical 3D continuum specified for the aforementioned shear stresses and the bending-induced normal stresses. Details on how to compute σ bend xy and σ bend xz can be found in Sect. 3 of this paper and in [28].

Introduction of cross-sectional rigidities
Insertion of the virtual velocity field given in Eqs. (4)-(6), the virtual Eulerian strain rate field given in Eqs. (7)-(9), the stress tensor given in Eq. (10), and the constitutive relations given in Eqs. (48)-(50) into Eq. (2) allows for formulating P int in terms of so-called cross-sectional rigidities, see [28] for all related details. Proceeding as explained above, and considering d 4ω x /dx 4 ≈ 0, see [28], yields with axial rigidity E A, bending rigidities E A yy , E A zz , and E A yz , torsional rigidity G I TI , axial warping rigidity E A ψψI , and transverse warping rigidity G A R . Furthermore, these rigidities contain area moments which are yet to be defined. In particular, A yy , A zz , and A yz are the second-order area moments, with j, k ∈ {y, z} ; I TI is the second-order area moment of primary torsion, A ψψI is the second-order area moment of primary warping, and A R is the second-order area moment of secondary warping,

Numerical solution
The mathematical framework elaborated in Sect. 2 is solved numerically, in the form of a procedure comprising four computation steps. This strategy has been elaborated for the first time and in minute detail in [28]. In the following, a concise summary of the key aspects related to this numerical solution strategy is presented, with emphasis on the novelties of this paper.

Computation of primary and secondary warping functions
For computation of the primary and secondary warping function, as well as the coordinates of the center of twist, the shear stress components given in Eqs. (49) and (50) are split into primary shear stresses associated with ψ I (referred to as σ I xy and σ I xz ) and secondary shear stresses associated with ψ II (referred to as σ II xy and σ II xz ). Then, as shown in [28], the following boundary value problem related to the primary shear stresses can be derived: ∀ y, z on A: where ψ GC I (y, z) is the primary warping function quantified with respect to the geometrical center; note that ψ GC I (y, z) is related to ψ I (y, z) via [28] ψ GC I (y, z) = ψ I (y, z) − y CT · z + z CT · y + Furthermore, requiring traction-free lateral surfaces, a boundary condition can be derived, complementing Eq. (56), reading as ∀ y, z on C: n y · ∂ψ GC I (y, z) with n y and n z as the components of the normal vector of the lateral surface. The boundary value problem defined by Eqs. (56) and (58) can be transformed into the corresponding weak form, where (y, z) is a test function belonging to the Sobolev space [42]. Analogously, a boundary value problem related to the secondary shear stresses can be derived, defined by ∀ y, z on A: and ∀ y, z on C: C ∂ψ II (y, z) ∂ y · n y + ∂ψ II (y, z) ∂z · n z ds = 0.
The corresponding weak form reads as Both Eqs. (59) and (62) are solved based on a two-dimensional (2D) Finite Element (FE) analysis [43,44], using for that purpose isoparametric four-noded quadrilateral finite elements with bilinear shape functions [30,45], giving, eventually, access to ψ I (y, z) and ψ II (y, z), respectively, and to warping-related cross-sectional quantities such as the center of twist.

Computation of displacements and rotation angle along the beam axis
In order to actually compute the three displacement components, u GC x (x), u GC y (x), and u GC z (x), as well as the beam rotation angle, φ x (x), a one-dimensional (1D) FE analysis [43,44] is employed. To that end, the beam is divided into n E 1D finite elements E m , m ∈ [1, n E ]. Along each of the finite elements, the distributions of the virtual velocities are approximated through shape functions and nodal vectors. Before presenting the respective equations, we point out that for discretizing transversal virtual velocities (and displacements), cubic shape functions are used, based on so-called Hermite polynomials [46,47]. However, such shape functions are not compatible with a nonzero deviation moment A yz . In order to nevertheless apply them, the transversal virtual velocities (and displacements), as well as all related quantities occurring in Eq. (51) must be expressed in a base system with the transversal base vectors coinciding with the principal axes η and ζ of the cross section, see also A. Then, the element-specific virtual velocities read aŝ (1) , 0, 0, 0, 0, 0, 0,v GC,m x, (2) , 0, 0, 0, 0, 0, 0  (1) dx , 0, 0, 0, 0, 0,v GC,m ζ, (2) , − dv GC,m ζ, (2) dx , 0, 0 andω (1) dx , 0, 0, 0, 0, 0,ω m x, (2) , with indices (1) and (2) denoting the two nodes of finite element E m , and ξ denoting the local coordinate (measured with respect to the center of the element), which relates to the global coordinate x via the Jacobian determinant |J| and the element length l m , Furthermore, N lin , N cub posD , and N cub neg denote the vectors composed of the aforementioned linear and cubic shape functions used for approximation of the virtual velocities within one finite element, see [28] for a detailed description. Importantly, the actual displacement components and the actual beam rotation angle are discretized in full analogy to Eqs. (63)-(66). Inserting Eqs. (63)-(66), the analogously defined u GC and φ x (x), as well as the shape functions according to [28] into Eq. (51), formulated in terms of principal axes η and ζ , summarizing over all finite elements n E , and assembling all virtual and actual element-specific quantities into corresponding global ones (e.g., the element-specific virtual velocity vectorv GC,m x is assembled into one single virtual velocity vectorv GC x ) yields the following FE approximation of P int : with the (matrix) assembly operator A. Furthermore, K m K m E A ψψI and K m G A R denote the element-specific stiffness matrices, see [28] for details.
Similarly, the virtual power performed by the external forces, see Eq. (22), can be expressed in terms of a corresponding FE approximation, reading as In Eq. (69), the elastic foundation-related element matrices read as and and with η 1 (x) and η 2 (x) as the principal axes-related counterparts of y 1 (x) and y 2 (x), as introduced in Sect. 2.3 and Fig. 2. Furthermore, the element-specific load vectors occurring in Eq. (69) are defined as Inserting Eqs. (68) and (69) Equation (84) can be solved for the sought-after actual displacement quantities. For realizing a tensionless elastic foundation, the following iterative approach was implemented: Step I: Assume, as initial guess, that the whole bottom surface of the beam undergoes a deformation in zdirection being larger than zero, relating to case (a) depicted in Fig. 2. Hence, ∀x : , implying corresponding values for η 2 (x) and η 1 (x). Consider y 1 (x) and y 2 (x) as reference distributions y ref and φ x (x), and the respective derivatives, considering for that purpose Eq. (84), as well as y ref 1 (x) and y ref 2 (x).
Step III: Compute u z (x) and, based on the considerations shown in Fig. 2 are computed, see [28] for details.

Computation of cross-sectional stress distributions
Inserting the forces obtained as described above into Eq. (48) allows to define the bending-related share of Additionally, we adopted the ansatz for σ xy (x) bend and σ xz (x) bend proposed by Kourtis et al. [48], and where ν is the Poisson's ratio, and a y (x) and a z (x) are defined as and Equations (85)-(87), together with Eqs. (88) and (89), are inserted into a suitable equilibrium condition, and into a boundary condition reflecting traction-free lateral surfaces, both of which can be expressed in terms of a corresponding weak form, see [28] for details. The latter is again solved by means of a 2D FE analysis, using four-noded quadrilateral elements and bilinear shape functions, see [30] for details.

Model input data
We considered a double-clamped beam of length l = 50 m, representing a tramway rail of cross section Ri60R1 [41], see Fig. 2, for which steel-typical elastic constants were chosen, namely a Young's modulus of E = 2.1 × 10 5 N/mm 2 , and a Poisson's ratio of ν = 0.28. As for the elastic foundation, several scenarios were studied, ranging from a very stiff elastic foundation to a comparably compliant elastic foundation supporting the whole beam, and including also a discontinuity in the elastic foundation in the middle of the beam, representing a step-wise transition from a stiff to a compliant elastic foundation. Thereby, the upper limit of the foundation modulus was chosen according to the respective suggestions made in [49], k max EF = 0.4 N/mm 3 . The compliant (or soft) elastic foundation, in turn, was assumed to exhibit a considerably lower foundation modulus, namely k min EF = 0.1 N/mm 3 . Furthermore, an intermediate elastic foundation was also considered, exhibiting a foundation modulus of k med EF = 0.25 N/mm 3 .

Loading conditions
As for the applied loading, the loads to which the studied beams were exposed were chosen such that they reasonably represent the loads transferred from tramway wheels onto rails. In the Viennese tramway network, the maximum occurring wheel load amounts to P z = 59.69 kN [50]. When assuming that the steels making up the rail and the wheel feature the same elastic constants, a simplified form of the Hertzian contact theory [51,52] gives access to the mean contact pressure acting on a small rectangular area, In Eq. (90), r wheel denotes the radius of the wheel in contact with the rail, and 2 · b cont denotes the width of the rail/wheel contact area. Furthermore, typical trams used in the Viennese tramway network exhibit a wheel radius of r wheel = 345 mm [50] and the rail head radius of rail profile Ri60R1 [41,[52][53][54] implies b cont = 5 mm. Evaluating Eq. (90) accordingly yields q mean = 622.15 N/mm 2 , which corresponds to a contact area length of l cont = P z / (q mean 2 b cont ) ≈ 10 mm. The respective distributed transverse load s z (x) is obtained through multiplying the contact pressure q mean with the width of the contact area, s z (x) = 6221.5 N/mm, which was applied over the length l cont . Importantly, application of P z , and consequently of s z (x) occurs as shown in Fig. 3, implying an offset of 16.04 mm between the line of action of those forces and the center of twist. Hence, P z induces a torsional moment M T , which is taken into account in terms of a distributed torsional moment m T (x) = 99.79 kNmm/mm, applied between −l cont /2 and l cont /2 (both with respect to the point of load application). It is emphasized that the chosen beam length of l = 50 m ensured that the bearing conditions did not have any relevant effects on the resulting forces and stress distributions.

Load case I: point load on continuous elastic foundation
One wheel load, as defined in Sect. 4.2, was applied at x load = l/2, whereby the beam was resting on a continuous elastic foundation; with the corresponding stiffness varying as defined in Sect. 4.1, see Fig. 3.
As expected, a more compliant elastic foundation leads to higher displacements u y (x) and u z (x), see, e.g., the two plots in the top row of Fig. 4. A similar behavior is also observed for the bending moment M y (x) and the classical torsional moment M T (x) minus the primary warping-related torsional moment M ψ T (x), see the plots in the second row of Fig. 4. The distributions along the beam of the warping-related bimoment M ψ II (x) and trimoment M ψ III (x) are depicted in the third row of Fig. 4. It is particularly interesting and instructive to compute the cross-sectional distributions of the stress tensor components σ x x , σ xy , and σ xz . Notably, such distributions were shown in previous works [28,30]. In this work, for the sake of conciseness, we only show the minima and maxima along the beam axis, see the plots in the fourth, fifth, and sixth row of Fig. 4. In order to allow for better comparability of the stress states along the beam axis, the von Mises stress was computed, through The respective mean and maximum values in the cross sections along the beam axis are shown in the last row of Fig. 4, confirming that a lower stiffness of the elastic foundation implies higher stresses in the rail.

Load case II: point load on discontinuous elastic foundation
Again one wheel load was considered, acting onto a beam which was resting on a discontinuous elastic foundation, with a k EF (x < l/2) = k min EF and k EF (x > l/2) = k max EF . Furthermore, three locations of load application were considered, namely x load = 0.45l, 0.5l, 0.55l, see Fig. 3. The respective resulting distributions of displacements, bending moments, torsional moments, and stresses are shown in Fig. 5.
On the one hand, the qualitative trends described already in Sect. 4.3 are confirmed; hence, the stresses related to load application on a softer elastic foundation are higher than the stresses related to load application on a stiffer elastic foundation. On the other hand, it turns out that all computed quantities are clearly asymmetric, owing to the fact that they are influenced both by regions of softer and stiffer elastic foundation. This asymmetry turns out to be-not surprisingly-most pronounced for x load = 0.5l.

Summary and conclusions
In this paper, we have presented a new approach allowing for calculating how an elastic foundation affects the force and stress distributions of beams resting thereon. The derived mathematical framework, taking the effects of restrained torsion rigorously into account, was numerically evaluated by means of a sequential, two-step Finite Element scheme, solving the equations first on the cross-sectional level in order to obtain the crosssectional values, then in longitudinal direction of the beam, and finally again over the cross section at selected locations along the beam axis. As compared to conventional, full three-dimensional Finite Element analyses, this new approach is significantly more efficient, see also [28] for a quantitative comparison. Hence, the present paper adds a potentially valuable aspect to the growing body of scientific literature on beam theory-related Finite Element analysis; see, e.g., [55][56][57][58].
The question raised in Sect. 1 of this paper was whether discontinuities in the elastic foundation may be potential triggers of mechanical failures occurring in tramway rails. The results presented in Sect. 4.4 have not shown any particular stress peaks caused by such discontinuities themselves. However, an almost step-wise increase of stresses upon moving from a rail section resting on a stiff elastic foundation to a rail section resting on a soft elastic foundation was observed. The plots in the last row of Fig. 5 show that the maximum von Mises stress increases (for the studied jump in the foundation modulus) by more than 50 %, while the mean von Mises stress increases by more than 30 %. Hence, this study supports the hypothesis that maintenance works to which tramway networks are exposed to over the years may indeed promote the occurrence of stress peaks leading in further consequence to rail failures; especially when taking into account that materials used for elastic foundations are nowadays stiffer than in former times, while wheel loads have become higher.  ψ III , as well as the minima and maxima of stress tensor components σ xx , σ xy , and σ xz , and the maximum and mean values of the von Mises stress σ vM ; the abscissa of all diagrams indicates the coordinate in beam direction, given in meters which is used to replace in Eq. (98) dζ by a term containing dη instead, allows to eventually derive P ext,EF = − where 1 η(x), 2 η(x), and 3 η(x) have been defined already in Eq. (79). Note that the integral terms occurring in Eq. (100) can be found in discretized format in Eq. (69).