LOAD CARRYING CAPACITY OF MASONRY ARCH RAILWAY BRIDGES AT THE SERVICEABILITY LIMIT STATE

. This article focuses on masonry arch railway bridges. The main topic is the algorithm, which was developed for the purpose of a non-linear arch bridge analysis. The sensitivity analysis of the input parameters of the calculation was carried out. Subsequently, the study was carried out for the comparison of four calculation methods. For this comparison, a set of arch bridges were chosen. Two methods used commercial software, and the other two methods were developed. One of the developed algorithms uses a simple linear calculation. The second developed algorithm uses an iterative approach to handle material and geometrical non-linearity. The first commercial software method uses 2D non-linear models, the second method is a limit collapse analysis using the software LimitState:RING. These three methods were developed to handle the SLS (serviceability limit state) criteria, and their results were compared to the result of LCC at the ULS (ultimate LS) using LimitState:RING.


Introduction
Masonry arch bridges are some of the oldest types of bridges in the world. Within the last century, these structures stopped being constructed. Therefore, the arch bridges, which are now still used as railway bridges, are from fifty to one hundred and fifty years old. For this reason, the tensile strength of the mortar should be considered to be at a zero level. On that account, the geometry of a thrust line of an arch is impacted. The linear calculation cannot handle the crack opening, which causes significant changes to the final thrust line.
The developed algorithm (hereinafter referred to as MVo) has two modes. In the first mode, just one linear calculation is completed, and the forces are verified directly after the first step. In the second one, the arch verification is done after the convergence of the steps of the non-linear analysis. The developed algorithm uses a beam model. This has two advantages: the algorithm is simplified and there is a lesser amount of computation time. The non-linear calculation means that the load carrying capacity, hereinafter referred to as LCC, is calculated iteratively for all load steps of the moveable load, which is time consuming even when using the beam model (for this article, around a hundred arch bridges were analysed using a usual computer, and the total computation time was around four hundred hours). The disadvantage is that the model parameters, such as the geometry of nodes and cross-section properties, must be completely renewed in every step of the calculation, therefore creating the algorithm of such a calculation was time-consuming for the programmer. The modelling of arch bridges using non-linear beam elements was proposed by [1] and used, for example, in [2].
The LCC is assessed according to [3] as a multiple of load model 71, which can pass the bridge safely (under conditions described in Section 5). This multiple is called Z LM 71 .
The masonry can be modelled in several ways. For example, both the atomic model and detailed 3D model are simulating the masonry elements and mortar separately. The simplest method is to homogenize, in a suitable manner, the masonry in one material. This approach is referred to in [4], the main code of masonry structures, which is, according to [5], called "yield line theory" and is based on laboratory tests conducted on brick masonry walls subjected to lateral loads, showing that failure takes place along a well-defined pattern of lines. Other homogenization techniques are described in the paper [6], which addresses many different homogenization methods available in the literature, showing their advantages and disadvantages. The strategy employed in [7] is based on the well-established first-order homogenization schemes, e.g. [8].
The material non-linear behaviour of masonry is impacted mainly by the low or zero tensile strength of the masonry [9]. The tensile strength can be usually neglected, but, for example, in the article [10], the tensile strength is used and the increase of bearing capacity is shown. At the SLS (serviceability limit state), the crushing of masonry is not allowed, at the ULS, the crushing can occur and in RING (LimitState:RING -see [11] software), even plastic hinges and the sliding of blocks of the arch can occur, until the moveable mechanism of structure, and therefore collapse of the arch occurs. The geometric non-linearity is handled only in the created algorithm, MVo. The effect of geometric non-linearity cannot be Δ Figure 1. The model used in the first step of the calculation.
neglected just for the case of arches of a large span and low ratio sagitta over the span of the arch. In this article, a sensitivity analysis of key input parameters is going to be performed using the MVo code. Afterwards, the results of four calculation methods will be compared: • MVo code -modelling using beam elements, ▷ Linear (SLS analysis), ▷ Non-linear (SLS analysis), • Scia -planar elements were used (SLS analysis), • RING (ULS analysis).

The used method -MVo
The algorithm which was created using the MAT-LAB software can handle circular arch bridges or arch bridges given by the set of points representing arch axis and input of the desired degree of a polynomial, which will be used for finding a smooth geometry of the axis of the arch by the least square method. It is assumed that such a set of points can be obtained, for example, by geodetic measurements. In this article, the results of circular arch bridges will be presented. See the example of the Legion Bridge in the literature [12], for which the geometry was fitted by the least square method of the polynomial of 8th degree. The main property of the model is that it uses beam elements. From 128 to 256 elements are used. The axis of the original arch is considered for the first step of the calculation. Due to crack opening, the geometry changes. According to these geometry changes, the cross-section properties (area and moment of inertia) are changed correspondingly. The new geometry of every node of every step is "guessed" as a geometry of a thrust line of the previous step.
The springs representing the backfill behaviour act linearly for compression in the soil and the stiffness in the tension is considered to be zero. The fact that the stress in the soil is either tension or compression is assessed by the deformation from the live load. The spring acts only in a horizontal direction. The stiffness of the spring is calculated from E def :  Figure 3. View of whole model.
Where K spring is the stiffness of the spring in the given node, E def is the deformation modulus of the soil, ∆Z is one-half of the horizontal projection of distance between two adjacent nodes to the given node, B is the width of the vault, ∆L is the length of the substituted soil (the distance to each node), see Figure 3.

Assumptions of the calculation:
• The supports at the arch springing are considered infinitely stiff, no deformation is allowed, i.e., no uneven settlements of supports are considered.
• The Bernoulli-Navier hypothesis for beam elements is considered.
• Constant soil characteristics are considered for all points of the backfill.
• There is no inflection point in the geometry of the arch in the case of the polynomial shape.
• The new position of nodes must lie on a normal line of the curve of the original axis geometry.
• In the study of the sensitivity of input parameters, the arch is symmetric, the longitudinal slope of alignment of a railway is considered to be zero. The thickness of the arch is constant.
See the loading conditions in Figures 4 to 7. The structure is loaded by self-weight, the weight of the backfill, ballast, earth pressure, and live loads.
The live load dispersion (load distribution) is done in the same way as in LimitState:RING. The Boussinesq distribution is used, see the [11]. The load is distributed at angle ϕ Ball in the ballast and ϕ Back in the backfill. This distribution is displayed in Figures 6 The developed algorithm can be used for the calculation of LCC at the SLS limits or the ULS limits. In this article, only the SLS is handled. The collapse load can be calculated in other software. In this article, RING is used, see [11] and [13]. The used equations are described in [14], and [15]. The geometric non-linearity is handled by the second order analysis, the material non-linearity is handled by the crack opening and iterative changes of the crosssection. The smooth geometry is needed for all steps of a non-linear calculation. This is a typical property of the beam model. Even a small aberration from the smooth geometry causes the unreal values of internal forces -especially bending moment and shear force. In the algorithm, during tens of steps, the small aberration always occurred, which led to unreal results. There are two ways of handling this problem. The first is to use many (at least thousands) elements, the second way of handling this is to smooth all the new coordinates of nodes. The first way is very   time-demanding. Therefore, smoothing was used. On the grounds of that, in circular arch bridges, the functions of bending moments are functions of sine and cosine functions, Fourier curve fitting was used. This curve fitting for the geometry of the new calculation step leads to the fastest calculation convergence. During the first few steps, the fitted data have some errors (as can be seen in the figure). The error is equal to the chosen precision ε, when the calculation converges. eT h4F it is e T h before fitting, f itted is after the curve fitting. In Figure 8, the fitting is done for an example of the third step of vault bridge calculation. In the next steps, the error is usually close to ε and the curves look identical. The used parameters in the study using MVo code:

The control method taken from article [15] -Scia
The second group of used models is the group using 2D planar elements. The masonry blocks are represented by planar elements with a linear behaviour. The joints between the blocks are modelled by a set of beam elements that are "compression-only". The soil is modelled by elements with Young modulus equal to E def , which allows us to also model the live load dispersion and all other effects of backfill, except for passive earth pressure. For the calculation of the SLS, these deformations are small, and passive earth pressure is not activated. As mentioned before, this analysis was done using the Scia commercial software [16].
For the details of this model, see [12], [15]. Scia parameters used in this article: E def , E masonry are chosen the same as in MVo code, ν M ason = 0.2, ν backf ill = 0.333; ν -Poisson's ratio.  mechanism. The resulting collapse mechanisms can be divided into the following cases: collapse by the crushing of masonry; collapse by the opening of cracks; collapse by shearing between the blocks. The method was described in [17] and [18], the modes of the vault collapse are also described in [19]. Figure 10 depicts the M-N diagram -the combinations of the acceptable moment and axial capacities of a structural member at the ULS.

The RING method
The program RING gives two options for the analysis. The compressive force in the vault is transferred by the joints: (i) Through an infinitely thin strip of stone at the edge of the vault (external when collapsing towards the inside of the arch, internal when deflecting outside the arch) if infinite strength of the masonry is assumed.
(ii) Through a rectangular strip that represents the stress at which the masonry is crushed.
For the case of infinitely stiff blocks, the curve of the M-N diagram is linear. The black dot is the point of rotation of blocks when creating a moveable mechanism.

The method of verification
The first main indicator of how the cross-section of height H is loaded is the position of the thrust line. The thrust line is a locus of points through which the resultant force goes. For each point of the vault, it can be (for the case of the beam model) calculated as e T h = M/N , where M is the bending moment and N is the normal force acting on the cross-section -from the load in a given combination. The higher bending moment acting in the cross-section, the higher is e T h and the lesser is the compressed area. For fulfilling the criteria of the SLS, e T h must be lower than H/3 according to [20], for fulfilling the criteria of the ULS, e T h must be lower than H/2. Eccentricity larger than H/2 means that the thrust line lies out of the crosssection, which leads to the collapse of the structure (if the tensile strength is neglected). This is illustrated by Figure 11: The method of obtaining the e T h depends on the chosen model. For the case of 2D planar elements, the verification can be done by checking that the height of the compressed area is at maximum 0.5H. For the case of the linear beam model, e T h is obtained in a single step. For the case of the non-linear beam model, the e T h varies for each step of the calculation and only the steps that fulfill the criterion of convergence (and equilibrium conditions as well as limit strain conditions) can be taken as a final value of e T h . The second main indicator of the state of the vault is normal stress that is acting on the cross-section from the given load. In the SLS, the stress should be less or equal to 0.45f k . This criterion assures that, during the usual loading conditions, no point of the structures will crush. The crushing is permitted in the case of the ULS, a state of collapse of the structure.

Deformation modulus of the soil E def
The chosen scope of E def corresponds to the variation from the clay-sand to compacted gravel of an ideal grain size. It can be seen that the parameter is most sensitive for the case of arch bridges with a high ratio of v/L. The lower the ratio v/L, the lower is the sensitivity to E def .

Coefficient of friction µ
The sensitivity analysis of the coefficient of friction between the blocks of masonry elements is done just for circular bridges. As it was shown in the article [13], the shear resistance depends mainly on the shape of the vault. The circular bridges usually have enough shear resistance. This fact is proven also by Figure 13, which shows the results of a calculation in the SLS. The ratio V /N -shear force over normal force -which should be less or equal to 0.4 due to [4] and which should be less or equal to 0.6 due to experimental data are for all the investigated arch bridges not higher than 0.225. Hence the LCC is not impacted by shear strength. However, for the lower values of the coefficient of friction, the sensitivity to change of this parameter is high.

The thickness of the arch
The thickness of the arch is the most sensitive parameter. That is why there should be put an effort in getting this parameter during the diagnostic survey. The thickness of the arch is sensitive both due to the maximal stress and the eccentricity of the load. The higher the thickness, the higher the range where the thrust line can occur. In the non-linear calculation, the crack opening, and therefore the possibility of finding the ideal geometry raises with increased thickness.

Characteristic strength of masonry
The characteristic strength of masonry is a parameter with similar sensitivity to the change in E def . In the SLS, there is always some limit for which increasing the strength does not increase the LCC, because the maximal eccentricity is a decisive criterion. This can be seen from Figure 15, but mainly from the study and comparison of the three mentioned methods for handling the SLS criteria in Section 7.

The results of the study of the set of arch bridges -A comparison of the four forementioned methods
The LCC of bridge spans L = 2.5, 5 and 10 m was calculated, the characteristic strength of masonry was (1 for special cases) 2, 4, 6, 8, and 10 MPa, sagitta v was always considered to be L/2 and L/4. The depth of backfill p was always considered to be 0.5, 1, 1.5 m.
Owing to the fact that Scia results were considered for a control reason, the LCC of the medial thickness (of total 3) was not calculated. The difference between the MVo code and the Scia model is that in Scia, the behaviour of the backfill is linear -in the horizontal direction as well as the vertical direction. Talking about the arches, where the height of the backfill is similar or even larger than the span of the arch, the results are significantly impacted by the soil behaviour. The problem becomes more influenced by the impact of soil and its modelling rather than the behaviour of the masonry arch. See the result of RING modelling in Figure 18, MVo non-linear code in Figure 19, the results of MVo linear code in Figure 20, the results of control Scia calculation in Figure 21 and a comparison of the RING method with the MVo code in Figure 22. For the purpose of clarity, the comparison of four forementioned method was added in two separate plots: 16 and 17.

Result discussion
It was confirmed that the results of the two models, which should have similar results, are, indeed, similar. The biggest difference is for the smallest arch of the span of 2.5 m. For this geometry, the calculations are impacted greatly by the behaviour of the backfill. The linear backfill behaviour helps the arch in both vertical and horizontal directions. For the arch with a high effect of backfill, special modelling should be done using special geotechnical software. It should be noted that the results of the control modelthe Scia model -are imprecise due to the reading of the graphical results. Especially the verification of the maximal eccentricity is sensitive to interpretation of the height of the compressed area. Seeking the LCC iteratively is time-demanding when using the Scia software. Conversely, the MVo code is created for calculating the LCC and results are obtained simply, and straightforwardly. The verification is done with precise numbers. c The linear calculation must obviously give a lower LCC. In the non-linear calculation, the axis of the arch is "optimised", the geometry is updated to decrease the bending moments in the next step. The limits of changing the geometry are given by the arch extrados, intrados, and the SLS criteria. The resulting bending moments depend mainly on the crack opening. If the height of the cross-section decreases and is close to the limit -one-half of the cross-section height, then the moment of inertia decreases eight times. The bending moment in the next step decreases up to eight times. The optimization of the geometry is done for every load step, every position The same calculation using the MVo code can take the engineer about 5 minutes. In the MVo code, all the input parameters can be changed easily, even the thickness of the arch or other geometry inputs. The iterative process of finding the multiple Z LM 71 (including moving load) is implemented and the engineer gets the result directly. Finally, the advantage of the MVo code is that the geometric non-linearity of the structure and the non-linear behaviour of the backfill are included.
The authors tried to find the dependence of LCC using RING on LCC using the MVo code. Such an obvious dependence was not found. The curves of the LCC ratio Z RIN G /Z M V o are similar in shape, but the tendency to grow or decrease differs for every analysed arch. For most of the cases, the ratio is the lowest for the highest depths of the backfill. For the majority of the cases, the LCC is higher at the ULS, usually up to two times higher. A comparison of the sensitivity of input parameters due to the SLS (MVo) and the ULS (RING -see in the [15]) give very similar results when comparing the sensitivity to the change of E def in the SLS and sensitivity to the change of angle of the internal friction in the ULS. The same applies to the coefficient of friction between the blocks of masonry, ratio H/L, and strength of masonry. Individual curves give a slightly different behaviour, but in general, the accordance of the compared plots was very good. In the RING result plots, there is no obvious breakpoint between the growing and constant segment, which can be seen in the SLS methods. The reason is that in the RING analysis, the criteria are different: low strength arches collapse by the crushing of the whole cross-section. The middle range of strengths is impacted by the crushing, and with growing strength, the crushed area decreases. When the strength is higher than usual masonry elements can have or is close to infinite, the criterion of the maximal eccentricity of the load is decisive and the constant segment of the plot occurs.

Conclusions
The comparison of sensitivity to change of the input parameters due to the SLS and the ULS gives very similar results. This comparison was done for key parameters, which were proven to be the most sensitive in the article [13]. The comparison of analysis using Scia software with 2D planar elements to the MVo results showed a good agreement except for small spans of the arch, where the behaviour of backfill affects the final LCC significantly. For this reason, a further study of the backfill behaviour with an advanced methods designed specifically for the soil behaviuor should be conducted. The comparison of analysis using the MVo non-linear method to the RING results showed no particular dependence. The LCC from the RING method is mostly higher than the LCC from the MVo method. The RING results are up to two times higher. The linear calculation gives very conservative results, LCC is usually two to three times lower than from the non-linear model. This is because of the fact that the non-linear model can change its geometry to minimize the bending moments and the only boundary conditions for the geometry of the arch are the extrados and intrados of the assessed arch. The effect of crack opening, which reduces the bending moment at the critical points of the arch, is also increasing the final LCC.

List of symbols
ZLM71 is the multiple of load model 71, which can pass the railway bridge safely, ZRING is resulting ZLM71 from the software Limit-State:RING, ZMV o is resulting ZLM71 from the MVo code, Kspring is the horizontal stiffness of the spring in the given node, E def is the deformation modulus of the soil, Emasonry is Young's modulus of masonry, ∆Z is one-half of the horizontal projection of distance between two adjacent nodes to the given node, B is the width of the vault, H is the cross-section height, L is the span length (intrados of the vault), p is the depth of the backfill at the top of the arch, v is the sagitta (rise) of the arch, ∆L is the length of substituted soil (the distance to each node), see Figure 3, ϕ Ball is the angle of load distribution in the ballast, ϕ Back is the angle of load distribution in the backfill, ϕ is the angle of internal friction of the soil, K0 is the coefficient of earth pressure at rest, q is the vertical nodal force from the live load, ε is the chosen maximal error of resulting displacement, reaching this error terminates the iterative calculation, e T h = M/N is the eccentricity of the load in a given combination. Maximal acceptable e T h is displayed in Figure 11, M is the bending moment acting on the cross-sectionfrom the load in a given combination, N is the normal force acting on the cross-section -from the load in a given combination, V is the shear force acting on the cross-section -from the load in a given combination, eT h4F it is the function of the eccentricity of the load before the curve fitting, f itted is the function of eccentricity of the load after the curve fitting, γ F ill is the specific weight of the backfill, γStone is the specific weight of the masonry, γ Bal is the specific weight of the ballast, f k is the characteristic masonry strength, µ = 0.6 is the friction coefficient (between the masonry blocks), h Ballast is the height of the rail ballast (under the sleeper), L Sleeper is the length of the sleeper (in the transverse direction of the bridge), νMason is the Poisson's ratio of the masonry, ν backf ill is the Poisson's ratio of the backfill.