Acta Polytechnica https://doi.org/10.14311/AP.2022.62.0293 Acta Polytechnica 62(2):293–302, 2022 © 2022 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague A MINDLIN SHELL FINITE ELEMENT FOR STONE MASONRY BRIDGES WITH BACKFILL Petr Řeřicha Czech Technical University in Prague, Faculty of Civil Engineering, Thákurova 7, Prague, Czech Republic correspondence: petr.rericha@fsv.cvut.cz Abstract. Stone masonry bridges are difficult to analyse with commercial finite element (FE) packages for their specific heterogeneous composition. The stone arch is best modelled as a thick shell where there are predestined directions of tension failure, normal to the bed joints. A dedicated, very simple, Mindlin shell finite element is developed with five translational degrees of freedom per node. It features compatibility with linear isoparametric or constant strain elements for the backfill. Most bridges can be analysed with a sufficient accuracy assuming plain strain conditions. The element then simplifies to a Timoshenko beam element with three translational degrees of freedom per node. An application of the latter one to the bridge at Poniklá is presented. Keywords: Masonry arches, no tension joints, thick shell elements, octave script. 1. Introduction Masonry arch bridges still constitute a considerable part of the bridge stock of the transport infrastructure in Czechia. Many of them are valueable monuments of the technical expertise and craftsmenship of the past generations as well. Design standards for this sort of bridges do not exist worlwide, guides and manuals for load rating are provided instead, like [1] or [2]. They usually contain conservative approximate methods for a routine assessment and some guidance and recommendations for more complicated individual assessments based on advanced structural analysis methods. A structural analysis of masonry arch bridges is specific in that the ultimate limit state forces and stresses cannot be solved for assuming homogeneous isotropic linear elastic material. In contrast to that, design standards for reinforced concrete and steel structures admit such solutions. EN 1996-1-1 codifies in clause 6.1.1(2) that plane sections remain plane and the tensile strength of masonry perpendicular to bed joints is zero in the ultimate limit state. The standard is not compulsory for bridges but its conditions should be perceived as the minimum for them. Cracks in the bed joints of the masonry arches affect the stress distribution in a way that does not admit a linear solution. The simplest material model that can be accepted is the no-tension linear material where the no-tension condition applies to the normal stress in the planes parallel to the bed joints. This specific behaviour is difficult to simulate by material models available in general purpose program packages, [3]. Commer- cial packages do not include such material models. There were attempts to achieve the specific proper- ties of the masonry arches with a trivial model of a heterogeneous material – individual meshing of vous- soirs and joints by standard continuum finite elements with homogeneous material models. It is possible in academic works but not acceptable in design prac- tice for demands on the labour, software and input data on materials properties. This experience testi- fies that shell/beam elements using the rigid normals assumption (Timoshenko, Mindlin, Bernoulli-Navier and other assumptions) are indispensable for the solu- tion of the masonry arch bridges. Several dedicated codes have been written around 1990 for the analysis of masonry arches based on the Castigliano principle, CTAP [4], [5], rigid block as- sumption RING , [6], mechanism based ARCHIE, [7] and others. All of them assume plane strain condi- tions. These models simulate very well the observed behaviour of masonry arches in situ and in large scale model tests. They have been applied worldwide in practice despite their common drawback that they underestimate the interaction of the barrel, backfill and pavement. Two simple finite elements for thick shells are de- veloped based on the Timoshenko/Mindlin kinematic assumptions, one for plane stress/plane strain condi- tions, the other for 3D shells. In order to facilitate a seamless interaction with the backfill continuum, only translational degrees of freedom (DOFs) are em- ployed. The no-tension linear elastic constitutive equa- tions are integrated in a closed form for the normal stress in planes parallel to the bed joints. The shape functions of the elements are 2D and 3D variants of the same concept. The simpler 2D variant is presented first with an application to a bridge, a triangle thick shell element follows. 2. The TT element, Timoshenko beam with translational DOFs The displacement shape functions are sketched in Fig- ure 1 for the axial displacement ui,t of node i, for the 293 https://doi.org/10.14311/AP.2022.62.0293 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en Petr Řeřicha Acta Polytechnica u i,b u i,b u u i,t j,b X XX i,b j,b j,ti,t u j,t v vi Y Yi j l i,j d i, j vj i j X x y Figure 1. Shape functions, DOFs and conjugate nodal forces of the TT element lateral displacement vj of node j and for the axial displacement ui,b. Axial and lateral directions denote the element local x and y directions, indices t and b denote the top and bottom faces later also extrados and intrados of the arch. The local x axis goes from node i to node j. The deflection line is parabolic and symmetric. It’s single free parameter is determined so that the shear strain energy is minimum for any imposed node displacements. Analogous shape func- tions are valid for all four axial DOFs ui/j,t/b and the same applies to the two lateral DOFs vi/j . The cross-sections remain straight but not normal to the deflected beam axis. In an isoparametric beam element (linear variation of the cross-section rotations and shear deformation, axial strain and lateral de- flections along the beam axis), there would be no deflections for the ui/j,t DOFs. The element is known to lock for small depths of the beam owing to excessive work of the shear deformation. The present element features a constant shear deformation since the anti- symmetric (with respect to the element centre) part of the cross-section rotations is assigned to a constant curvature of the deflection line. The element DOFs are ordered in the {u} matrix, {u} = { {u}i {u}j } , {u}n =   un,t vn un,b   , n = i,j Linear normal and constant shear element strains are specified by top and bottom face normal strains εt, εb and shear strain γ, ordered in the matrix of internal strains ε = {εt,εb,γ}T . Shear strain γ and conjugate shear stress τ have opposite signs than usual in elasticity, γ = −(∂ux/∂y+∂uy/∂x)/2. The element geometric matrix [B] reads: [B] = 1 l [ −1 0 0 1 0 0 0 0 −1 0 0 1 l/2/d −1 −l/2/d l/2/d 1 −l/2/d ] Conjugate to node displacements {u} are the nodal forces {X}, conjugate to {ε} is the matrix of inter- nal forces {S}. The nodal forces in terms of nodal displacements read {X} = [B]T {S}({ε}) = [B]T {S}([B]{u}) (1) Constitutive equations for the cross-section are devel- oped separately for the bending, the first two members of {ε} and {S}int. ε ξ d d ζ t σ: >0 ε σ: εt <0 tz ( ) z ( )ζ b b σ: 0 ε 1 1 ξ 1 1 ε b >0<0 ξ d ξ ζ Figure 2. Stress/strain diagram, strain shape func- tions and normal stress diagrams of the TT element 2.1. Constitutive equations for bending No-tension linear elastic material is assumed, see the stress-strain diagram in Figure 2. The no-tension condition is applied to the normal stress in cross- sections and bed joints of the arch. Though simple, it makes possible the most frequent type of failure of the stone masonry arch bridges. The normal strain ε is determined by the top and bottom face strains εt, εb through the shape functions shown in the centre of Figure 2. ε = zt(ζ)εt + zb(ζ)εb Constitutive equations of the cross-sections are de- fined in terms of dimensionless functions st() and sb() of dimensionless arguments: st(εt,εb) = ∫ 1 0 zt(ζ)s(εtzt(ζ) + εb,xzb(ζ))dζ sb(εt,εb) = ∫ 1 0 zb(ζ)s(εtzt(ζ) + εbzb(ζ))dζ, (2) The cross-section forces are St/b = Ebdst/b = EAst/b (valid for both t and b subscripts, where E is the Young modulus, b the cross-section width, d the cross- section depth and A its area. Functions st/b() and their derivatives with respect to εt, εb have to be developed separately for four possible strain states: Algorithm start When εt > 0 and εb > 0 then all internal forces are 0, the cross-section is totally disintegrated. All stresses and internal forces vanish in the element. When εt > 0 and εb < 0, then (see the left strain diagram in Figure 2): ξ = εb εb − εt st = εbξ2/6, sb = (0.5 − ξ/6)ξεb ∂st ∂εt = ξ3/3, ∂st ∂εb = ( 1 − 2ξ εt εb ) ξ2 6 ∂sb ∂εt = (0.5−ξ/3)ξ2, ∂sb ∂εb = (0.5− ξ 6 )ξ−(0.5− ξ 3 )ξ2 εt εb When εt < 0 and εb > 0, then (see the right strain diagram in Figure 2): ξ = εt εt − εb 294 vol. 62 no. 2/2022 A Mindlin Shell Finite Element for Stone Masonry Bridges with Backfill i X =S +Vl/(2d)X =−S +Vl/(2d) X =−S −Vl/(2d)b t l d X =S −Vl/(2d) bj tjti bi V Y =−Vi j Y =Vj t b Figure 3. Element internal forces {S} and nodal forces {X} st = (0.5 − ξ 6 )ξεt, sb = εtξ2/6 ∂st ∂εt = (0.5− ξ 6 )ξ−(0.5− ξ 3 )ξ2 εb εt , ∂st ∂εb = (0.5−ξ/3)ξ2 ∂sb ∂εt = ( 1 − 2ξ εb εt ) ξ2 6 , ∂sb ∂εb = ξ3/3 When εt < 0 and εb < 0, then the cross-section is linear elastic and dimension- less generalised forces are st/b = (2εt/b + εb/t)/6 The cross-section stiffness matrix: ∂{sbend} ∂{εbend} = 1 6 [ 2 1 1 2 ] (3) Algorithm end The bending stiffness tangent matrix of the cross- section is denoted for a later reference [Dbend] = ∂{sbend} ∂{εbend} = Ebd [ ∂st ∂εt ∂st ∂εb ∂sb ∂εt ∂sb ∂εb ] (4) 2.2. Full cross-section stiffness matrix The shear force V is computed from the shear defor- mation γ in the Timoshenko beam. Linear elastic behaviour V = Gebdγ is assumed in the model, independent on the normal stress in the cross-sections. This assumption is mostly sufficient for cross-sections of the stone masonry arch bridges. The full cross-section stiffness matrix is then [Dcs] = [ [Dbend] {0} {0}T Gebd ] The matrix has to be adapted when more involved cross-section constitutive equations are necessary, for instance, shear failure of bed joints. The nodal forces {X} = {Xi,t,Yi,Xi,b,Xj,t,Yj,Xj,b}T of the element can be obtained in terms of the internal forces {S} either by expansion of [B]T {S} or by the equilib- rium conditions of the element in Figure 3: The di- mensionless cross-section stiffness matrix for bending [Dbend] = ∂{Sbend/∂{εbend} is not symmetric when the neutral axis lies inside the cross-section so that the x xe f bot tom fac e bottom face element e node i element f t t r rR 1 1 R B i X b X b X b X b X b X b f e f f e e Figure 4. Equilibrium of node i element and structure matrices will not be symmetric, too. The element stiffness matrix is [K]straight = Ebdl[B]T [Dcs][B] It is ready for application in straight beams, but for arches and other curved beams, it must be modified. 2.3. Arch nodes equilibrium The equilibrium of the nodal forces X⃗ibe = Xibex⃗e, X⃗ibf = Xibf x⃗f of the adjacent elements e and f at arch node i simplifies to Xibe + Xibf = 0 when the element’s unit length direction vectors x⃗l,e and x⃗l,f coincide, see the Figure 4. The node index i is omitted for brevity as long as the the development concerns just a single node i. Xbe and Xbf are the magnitudes of the respective forces with positive senses in ele- ment’s direction vectors. These forces are stored in the third components of the element nodal forces ma- trices and in case of parallel elements, the equilibrium equation at a node is {X}e + {X}f = {0}, where indices e and f indicate the left and right element, respectively. When the adjacent elements are not parallel, the resultants of the two forces add up to force R⃗. The force is decomposed into component R in the direction r and force T in the t direction. The component R can act anywhere on ray r since the arch is considered rigid in the transverse direction. It is thus added to forces acting at node i. Component T must vanish to keep the null moment with respect to joint i. T = (Xbex⃗e + Xbf x⃗f )·⃗t = = (Xbex⃗e + Xbf x⃗f )·(x⃗e + x⃗f )/|x⃗e + x⃗f | = 0, which implies Xbe + Xbf = 0, provided the two elements are not normal to each other. Force R is R = (Xbex⃗e − Xbf x⃗f )·r⃗ = (Xbex⃗e − Xbf x⃗f )·(x⃗e − x⃗f )/|x⃗e − x⃗f | = (Xbe − Xbf ) √ (1 − x⃗e·x⃗f )/2 and its vector R⃗ = Rr⃗ = 0.5(Xbe − Xbf )(x⃗e − x⃗f ) 295 Petr Řeřicha Acta Polytechnica The ray r of the force must connect points i and Bi, the intersection of the X⃗be and X⃗bf forces rays, otherwise node i would not be in equilibrium. Ray r is approximately radial to the arch. Vector g⃗ = 0.5(x⃗e − x⃗f ) (5) is an important property of node i. The matrix ex- pression for vector R⃗ is in terms of g⃗ R⃗ = {g⃗, −g⃗} { Xbe Xbf } Recall that Xbe denotes the third element of the nodal forces matrix of element e at node i and the analogue holds for Xbf . Conjugate displacements to Xbe and Xbf are ube and ubf , conjugate to R⃗ is the displace- ment vector u⃗ of node i. The principle of virtual work thus implies { ube ubf } = { g⃗ −g⃗ } ·u⃗ The same result can be obtained when the part of u⃗ in r⃗ direction, (u⃗·r⃗)r⃗, is projected on the x⃗e direction to obtain ube: ube = (u⃗·r⃗)r⃗·x⃗e = g⃗·u⃗. Projection on x⃗e yields the same expression but for a negative sign at g⃗. Expression for εb becomes εb =(ub,j + g⃗j ·u⃗j − ub,i + g⃗i·u⃗i)/li,j =(ub,j + {gj }T { ut,j vj } − ub,i + {gi}T { ut,i vi } )/li,j (6) where indices i and j denote the nodes of the element with the element local x axis from i to j. In the wake of it, the geometric matrix is modified to [B] = 1 l   −1 0 0 1 0 0gi,x gi,y −1 gj,x gj,y 1 λ −1 −λ λ 1 −λ   (7) with λ = l/d/2 and subscripts x and y indicating components of the vectors g⃗ in the global coordinate system. Different element lengths have to be used in the geo- metric Equation 7 and the volume integration implicit in 8. The top face nodes distance l is used in the for- mer case, the reduced length lred = lRax/Rt = lrratio is used in the latter case. Rax denotes the radius of the arch central axis, Rt is the radius of the top face curve. The element stiffness matrix for an arch element is [K]e = Ebdlred[B]T [Dcs][B] (8) -4 -2 0 2 4 -32 -30 -28 -26 -24 2587 2176 18311501 1182 930.61146 13401500 16411786 18891889 17731596 14151256 1104963.2840.8 739.7 732.1810873.5924.2963.4992.21012102410291030102710231018 1014 1013 1015 1021 1031 1043 32.79 33.3625.514.772.337 10.6723.2834.9345.5855.4464.5973.7681.7488.2593.2895.794.1388.34 122.6162.2170.4 173.1318.9525.3723.4734.6209.217011612211.5125.5 1667171753.04490.5 563.2447.9301.5 190.5135.7 Figure 5. Bridge at Ponikla, compression stresses [kPa] and depths of the cracks in the bed joints 3. Application to an arch bridge Exclusive translational DOFs facilitate the combina- tion of the TT element with simple 2D continuum elements. A dedicated matlab/octave code has been written to utilize this combination for the analysis of masonry arch bridges. The code features the sim- plest possible way how the interaction of the masonry arch, backfill and pavement can be assessed parallel with the principal pattern of failure of the masonry arch. The arch and pavement are modelled by the TT elements, the backfill by the classic constant strain triangle (CST) elements. Meshing of the backfill pro- vides the Persson and Strang mesh generator, [8]. The no-tension cross-section model defined in Section 2.1 supports the cracking of the arch bed joints, which is the dominant pattern of the masonry arch bridge failure. The code includes a simple, purely vector graphics output. It consists of 1150 octave command lines, not including the mesh generator. An application to the sandstone arch bridge at Ponikla in north Bohemia illustrates the code output in Figure 5. The span of the arch is 11.4 m. Two prin- cipal quantities are shown in each arch element, the maximum compression normal stress value inserted at the arch face where it occurs and the depth of the cracks in the bed joint. The load case is the design dead load combined with the tandem axle live loads with twice the design intensity level, i.e., 360 kN per wheel in lane 1 according to EN 1991-2. The TT element and the masonry arch FEM models based on it claim to be the simplest models available to assess the most frequent failure mode of the stone masonry arches and the interaction of the arch with backfill and pavement. 4. 3D Mindlin shell element TM The 2D models are often not adequate for the solution of masonry arches whose widths are mostly compa- rable to spans which makes the transverse variation of stresses and displacements nonuniform. A 3D ana- logue to the 2D TT element is, therefore, derived. The Mindlin assumptions, [9], and the simplest triangular facet shape are adopted, see Figure 6. The element local coordinate system features z axis normal to the facet plane. There is always a preference direction in the facet plane, the direction closest to the normals 296 vol. 62 no. 2/2022 A Mindlin Shell Finite Element for Stone Masonry Bridges with Backfill k x j d y z i u w t,x ut,y u u b,x b,y Figure 6. Geometry of the shell facet triangle ele- ment, the element local coordinate system and the DOFs at node i to the bed joints planes of the masonry arch. It is assumed forthwith that the arch is a cylindric segment shell and the bed joints normals are tangents to the directrix of the cylinder so that their directions lie in a single plane. The global x axis is assumed to lie in that plane, too. In terms of the bridge nomencla- ture, the global x axis is in the direction of the arch span. The element local x axis is chosen as the direc- tion closest to the global x within the element plane. The no-tension condition is applied to the σx (local x) normal stress component . Other options are possible and can be rather easily implemented in the code. It may be necessary in case of other than cylindric shells or more general failure modes and constitutive equations. Local y axis completes the element local right-handed cartesian coordinate system. All compo- nents are in this system up to Section 4.4, where two other systems are necessary. Five translational DOFs per node are shown in Figure 6 for node i. They are ordered in matrix {u} at each node, with node index omitted: {u}T = {ut,x,ut,y,w,ub,x,ub,y } Deflection w and in-plane displacement vector u⃗ are approximated by linear independent functions of natural triangular coordinates ξi. w = ∑ i wiξi, u⃗ = ∑ i u⃗iξi In literature, a common term for this kind of approxi- mation is isoparametric element, not to be confused with true isoparametric elements for 2D and 3D con- tinua. Derivatives of ξi are constants in the element area, ∂ξi ∂x = bi 2A , ∂ξi ∂y = ai 2A or ∂ξi ∂x⃗ = −n⃗i li 2A , where A is the element area (note that li/2A is the slope of the ξi surface upon the element area) and n⃗i = − 1 li { bi ai } . γ y γ x i j k a aa b b ij k i k y z w i φ i wjdw/dx d w /d y i j k A A A s n i i b j x w w Figure 7. Triangular coordinates and shear deforma- tion owing to deflection alone Chain differentiation yields then γw,x = − 1 2A ∑ i wiai, γw,y = 1 2A ∑ i wibi (9) The shear deformation γ⃗w owing to deflection alone is γw,x = − ∂w ∂y , γw,y = ∂w ∂x , (10) note that positive sense is indicated in the figure for each of the four rotations and γ⃗ is the rotation of the normal to the element plane embedded in the material before deformation towards the normal to the element plane after deformation with standard sign convention of its components. In particular, the motion (rotation) shown in the section view right to left in Figure 7 goes from the inclined dashed line to the horizontal and γx is thus negative. As in the TT element, the bending and shear defor- mation, are treated separately. 4.1. Bending and in-plane deformation Pseudocurvatures define the bending deformation. They are the derivatives of rotations of the normals to the element plane. The rotations of the normals in their turn depend on the in-plane displacements u⃗ of the top and bottom faces of the element. The element in-plane strains are specified in terms of the top and bottom faces (extrados and intrados) strain matrices {εt} and {εb} in analogy to the top and bottom uniaxial strains εt and εb in the Timoshenko 2D beam element. The element triangle geometry is defined in terms of distances ai = xk − xj , bi = yj − yk, where indices i, j, k are subjected to cyclic permutation and xi, yi specify node’s i positions in the local system, see Figure 7. The geometric equations 297 Petr Řeřicha Acta Polytechnica of the CST triangle applied to the top face read {εt} =   εt,x εt,y 2εt/b,xy   = 1 2A ∑ i   bi 00 ai ai bi   { ui,t,x ui,t,y } , (11) substitute b for t subscript to get the bottom face strain matrix. The in-plane strains vary linearly across the element depth d, {ε} = zt{εt} + zb{εb} (12) where zt = (d + z)/d = 1 + ζ and zb = −z/d = −ζ are linear shape functions of the element in terms of the element local 0 > z > −d coordinate with z = 0 at the top face of the element or of the relative dimensionless 0 > ζ > −1 coordinate. Plain stress conditions with σzz = 0 are assumed in the shell, then, the material stiffness matrix is [D] = E 1 − ν2   1 ν 01 0 1−ν 2   (13) for a linear elastic isotropic material. When cracked bed joints occur in planes normal to element’s x axis, the conditions become vague. Adjacent to the cracked joints an approximately uniaxial stress prevails. It is thus assumed forthwith that the material stiffness matrix in the cracked elements is [D({ε})] = E   s(εx) 0 01 0 1−ν 2   (14) where s() is the normalized bilinear stress/strain dia- gram in the left sketch of Figure 2. It is obvious that the element local x axis must approximately coincide with the normals to the bed joints of the arch barrel. The virtual work equation for cracked elements δw = ∑ i (δ{ui,t}T {gi,t} + δ{ui,b}T {gi,b}) = ∫ V δ{ε}T [D({ε})]{ε}dV = E ∫ V (δεxs(εx) + δεyεy + (1 − ν)δεxyεxy dV (15) delivers the expressions for the nodal forces {gi,t}, {gi,b} when strains are expanded in terms of the virtual nodal displacements using Equation 11. Functions st(εt,x,εb,x) and sb(), see 2 in Section 2.1, return dimensionles top and bottom cross-section forces conjugate to top and bottom strains (input parameters) of a unit depth cross-section with E = 1 and no-tension stress-strain diagram. With the aid of the functions, the integrals in the virtual work expression can be expressed in a closed form: δw = Ed 2 (δεt,xst + δεb,xsb + δεt,y (εt,y/3 + εb,y/6) + δεb,y (εt,y/6 + εb,y/3) + (δεt,xy (εt,xy/3 + εb,xy/6) + δεb,xy (εt,xy/6 + εb,xy/3))(1 − ν)) (16) Substitution by 11 and selection of individual nonzero virtual nodal displacements yields the expressions for the element nodal forces: gi,t,x = Ed2 (bist + ai(εt,xy/3 + εb,xy/6)(1 − ν)) gi,t,y = Ed2 (ai(εt,y/3 + εb,y/6)+ bi(εt,xy/3 + εb,xy/6)(1 − ν)) gi,b,x = Ed2 (bisb + ai(εt,xy/6 + εb,xy/3)(1 − ν)) gi,b,y = Ed2 (ai(εt,y/6 + εb,y/3)+ bi(εt,xy/6 + εb,xy/3)(1 − ν)) (17) Derivatives of the nodal forces gi,t/b,x/y with respect to the DOFs can be written in terms of the st/b functions, for instance ∂gi,t,x ∂uj,t,x = ∂gi,t,x ∂εt,x ∂εt,x ∂uj,t,x + ∂gi,t,x ∂εt,xy ∂εt,xy ∂uj,t,x = Ed 4A (bibj ∂st ∂εt,x + aiaj 1 − ν 2 1 3 ) Functions st, sb depend on x components of εt, εb only. The subscript x can thus be omitted in their derivatives. Formulas for the derivatives are in Sec- tion 2.1. When µ = 1−ν2 is introduced and common factor Ed4A is omitted for brevity, submatrices [Kbend,i,j ] of the element stiffness matrix associated with the in- plane DOFs {uplane}T = {ui,t,x,ui,t,y,ui,b,x,ui,b,y } are: i, i : [ bibj ∂st ∂εt + aiajµ/3 aibjµ/3 ajbiµ/3 bibjµ/3 + aiaj/3 ] i,j : [ bibj ∂st ∂εb + aiajµ/6 aibjµ/6 ajbiµ/6 bibjµ/6 + aiaj/6 ] j,i : [ bibj ∂sb ∂εb + aiajµ/6 aibjµ/6 ajbiµ/6 bibj/6µ + aiaj/6 ] j,j : [ bibj ∂sb ∂εb + aiajµ/3 aibjµ/3 ajbiµ/3 bibjµ/3 + aiaj/3 ] (18) For a linear element (uncracked bed joints) [Ki,j ]T = [Kj,i], which implies symmetry of the whole matrix [K]. For cracked bed joints ht,b ̸= hb,t, the symmetry is lost. 298 vol. 62 no. 2/2022 A Mindlin Shell Finite Element for Stone Masonry Bridges with Backfill 4.2. Shear owing to the rotations of normals The rotation of normals φ⃗ induces shear deformation, too. It would not induce any deflection in a consis- tent isoparametric element and then φ⃗ would simply be added to γ⃗ to obtain the total shear deformation. Such displacement pattern is known to imply locking in thin elements, like it does in the isoparametric vari- ant of the TT element. In the TM element, φ⃗ induces a quadratic field of deflection wφ, which must have zero values at all joints. The deflection should be compatible with the neighbouring elements along the boundaries. The antisymmetric part of the normal component of φ⃗ along an element’s side is assigned to deflection, the symmetric part to shear deformation, in analogue, to the TT element. Perfect compatibility of deflections along boundaries with neighbouring ele- ments is attained for elements lying in a plane. The implementation of the outlined displacement pattern is rather lengthy and it is skipped here. The resulting constant transverse shear owing to normals rotation is γu = ∑ i [ψi]({ui,t} − {ui,b}) [ψi] = ([ 0 −1/3 1/3 0 ] − 1 12A [ a2k − a 2 j bjaj − bkak −bkak + bjaj b2j − b 2 k ]) 1 d (19) The full transverse shear of the TM element is then {γ} = ∑ i ( [ψi], 1 2A [ −ai bi ] , −[ψi] ) {u}i (20) where indices i,j,k are subject to cyclic replacement. Formula 20 is necessary with thin shells where shear locking can occur. Masonry arches and other burried shells are mostly thick shells. Matrices [ψi] can be simplified by omitting the second summand in such applications. Experience with tests and applications testifies that keeping the second summand improves the element convergence. 4.3. TM flat slab element Simple linear constitutive equations are assumed for the transverse shear, Vx = Gγx, Vy = Gγy. (21) Vx is the standard shear force in element local y − z plane when looking along x axis and Vy is the standard shear force in element local x − z plane when looking along y axis. The contribution of the shear stiffness to the element stiffness matrix is written down in terms of the 5x5 submatrices associated with the two top vertex, one transverse and two bottom vertices translational DOFs of each node. Column matrices [φi] = [−ai,bi]T /2A and 2x2 symmetric matrices ψij = [ψi]T [ψj ] help to write down the shear contribution to the submatrices: [Kshear,i,j ] = GAd2   [ψij ] [ψi]T [φj ] −[ψij ][φi]T [ψj ] [φi]T [φj ] −[φi]T [ψj ] −[ψij ] −[ψTi [φj ] [ψij ]   (22) Complete 5×5 stiffness submatrices are the sums of the bending stiffness submatrices 18 extended by an empty third column and row inserted between the 2×2 subsubmatrices and submatrices 22: [Ki,j ] = [Kbend,i,j ]extended + [Kshear,i,j ] (23) The whole element stiffness matrix is the size 15×15 composition of submatrices [Ki,j ], but the submatrices are assembled in the global stiffness matrix so that the whole matrix is never set up. Cracking in the bed joints followed by development of virtual hinges in the arch have been considered the dominant failure pattern of masonry bridges since the early attempts [10], [11] to assess the load capacity. Other failure modes like shear sliding of the bed joints or transverse tension cracking of the arch require so- phisticated material models and properties which are almost impossible to obtain in the design practice. This stiffness matrix can be used for the solution of flat slabs loaded both in and out of plane but in shells, the bottom face DOFs ub,x, ub,y require a special treatment since they are not simply shared by the elements attached to a node. 4.4. TM shell element The DOFs at the bottom vertices of the elements connecting to a node do not lie in the same plane and are not independent of the DOFs at the top vertices. The equilibrium of the nodal forces at the top vertices is affected by the nodal forces at the bottom vertices. These two conjugate deficiencies must be removed. The normal to the shell surface is needed for a con- sistent formulation of the TM element. The normal is rigid in the frame of the Timoshenko-Mindlin shell theory, which implies that the displacements of all points of the normal share the same lateral displace- ment component. The exact definition of the normal direction n⃗ would be through the shell surface mathe- matical definition. For practical applications of the TM element, it is sufficient to define the unit normal n⃗i at a node i as the normalized sum of the normals of all elements connecting to the node, each element normal length proportional to the sine of the angle of the two adjacent sides of the respective element. Posi- tion vectors of nodes are denoted x⃗i and nodes of the connecting elements i, j, k, ordered counterclockwise when viewed from the top side of the shell. The sum is n⃗s,i = ∑ e n⃗e sin(αe) = 299 Petr Řeřicha Acta Polytechnica∑ e (x⃗j − x⃗i) × (x⃗k − x⃗i)/(|(x⃗j − x⃗i)||(x⃗k − x⃗i)|)e n⃗i = n⃗s,i/|n⃗s,i| Vector n⃗i’s positive direction is outwards from top surface of the shell. The vector also defines the tan- gential plane τi at node i, τi ⊥ n⃗i. The top vertices of the connecting elements share the displacement r⃗t,i of the top end of the rigid normal. The global components of vector r⃗t,i constitute the first three DOFs of the node i. The index of the node is omitted forthwith since a single node is considered. The compatibility of displacements at the bottom facets of the connecting elements demands that the projections on the τ plane of the bottom vertices displacements of all connecting elements share the same displacement vector r⃗b,τ ⊂ τ. A coordinate system is established in the τ plane such that axis 1 direction vector n⃗1 lies in the intersection of the global coordinates plane x−z and the τ plane, oriented as the global x axis. Axis 2 direction vector is n⃗2 = n⃗ × n⃗1 (vector product). The component expansion of vector r⃗b,τ is defined this way. The two components in the τ plane define the two complementary DOFs of the node. Note that they are not the displacement components in the global coordinate system x − y − z but in the local ‘tangential’ plane instead. The transformation matrix from global to these tangential vector components is [T ]τ =   {n1}T{n2}T {n}T   (24) The whole bottom vertices displacement vector is r⃗b = r⃗b,τ + n⃗(n⃗·r⃗t) = r⃗b,τ + n⃗ ⊗ n⃗r⃗t, but just the two components of r⃗b,τ ⊂ τ constitute the two complementary DOFs of the node. The out-of-τ- plane component depends on the basic three DOFs r⃗t of the node. Several coordinate systems (all cartesian) are employed. The global one, common for all, the τ system, common for a node and the local e systems of individual elements. The convention is adopted for component matrices that the matrix with components in the global system has no subscript τ, in the τ system, and subscripts e in the element systems. Furthermore, the transformation matrix [T ]e from the global to the element system has subscript e. The matrix form of the expression for the element DOFs is{ {rt}e {rb}e } = [ [T ]e 0 [T ]e{n}{n}T [T ]e[T ]Tτ ] { {rt} {rb}τ } (25) Recall that both the displacement component matrices {rb}e and {rb}τ have zero third components so that the third column and row of the transformation matrix in 25 can be omitted. The nodal forces conjugate to {rt}e and {rb}e are denoted {gt}e and {gb}e. Rigid normals to the shell surface imply that the equilibrium equations of the top and bottom joints of a node are not independent. Just a single equilibrium equation can be written in the direction n⃗ and it is assigned to the top end of the rigid normal – the top joint of the node. The two remaining equations at the bottom joint of the node include components acting in the τ plane. Component decompositions of the element nodal forces in the coordinate axes tripod n⃗1, n⃗2, n⃗ – the node local τ system, are added in the matrix equilibrium equation of the bottom joint of the node: ∑ e {gb}τ,e = ∑ e [T ]τ [T]Te {gb}e =   0 0 {n}{n}T {gb}e   The sum includes all elements connecting at the node. The third (tranverse) components are added to forces acting at the top vertex of the node to obtain the final nodal forces of the element:{ {gt} {gb}τ } = ∑ e [ [T ]Te {n}{n}T [T ]Te 0 [T ]τ [T ]Te ] { {gt}e {gb}e } (26) The transformation matrices of the nodal displace- ments and forces are transpose of each other, which testifies to their correctness. Note that the last scalar equations in 25 and 26 can be omitted and so can be the last columns of the tranformation matrices. The transformation matrices are then 5×5 in size. They are specific for each node of an element since the {n} and [T]τ matrices are different at each node. For an easy reference, they are denoted [Ti] for node i forthwith. The submatrices [Ki,j ] are transformed to the global coordinate system [Ki,j ]g = [Ti]T [Ki,j ][Tj ] (27) and assembled in the system matrix. 4.5. TM element code and test The TM element has been implemented in a dedi- cated Matlab/Octave code including a simple graphic output. Just the shell reference surface is drawn to keep the view readable. The maximum compression stresses in the bed joints and the relative depth of the cracks are inserted in each element. These values are sufficient to decide on the arch load capacity in the context of the present model. The graphic output uses vector graphics so that pictures can be zoomed in with a stable resolution. The mesh generator [8] is used for the shell surface meshing. The code is entirely self-contained, no input data file and pre- or postprocessing is necessary. It contains just about 660 octave command lines, not including the mesh generator. The pinched cylinder with rigid end diaphragms in Figure 8 is a popular benchmark to test shell ele- ments. The shell is thin, d/R = 0.01 so that the test is a severe one for a thick shell element. The bench- mark solution by double Fourier series with 80 terms in both directions was first presented in [12] based on 300 vol. 62 no. 2/2022 A Mindlin Shell Finite Element for Stone Masonry Bridges with Backfill 600 6 0 0 1 1 sym. 3 ri g id sy m . sym.x,y z Figure 8. Benchmark shell 0.2 0.4 0.6 0.8 1.0 0 16 32 488 elements along directrixr a ti o F E M d e fl e c ti o n /b e n c h m a rk Figure 9. Convergence of the FEM deflections, ratio to the benchmark deflection the work [13]. The Kirchhoff-Love kinematic assump- tions are used in the benchmark so that a definite deflection is obtained. Differences can be expected between the present and benchmark solutions in the vicinity of the pin force, in particular for fine meshes. The Mindlin assumptions imply infinite deflection be- neath the force in the continuum model, thus making the discretized models mesh dependent and inherently non-convergent. Owing to the three symmetry planes just shaded, 1/8 of the cylinder can be considered. Four mesh densities were considered with 8, 16, 32 and 48 elements along the directrix. The ratio of the loaded node deflection of the TM finite element solutions to the benchmark deflection 1.825 · 10−5 is shown in Figure 9. In spite of the correction of the transverse shear strain in Equation 20, the residual shear locking still persisted and affected the results in this thin shell, in particular, for low densiiy meshes. In order to reduce it further, the material shear stiffness in the transverse direction (element local x − z and y − z planes) was selectively lowered eight times. This has a negligible effect on the overall response of the FEM model since the contribution of the transverse shear compliance to it is small. At the denser mesh side of the Figure 9, the curve already tends to adopt to the infinite deflection of the continuum Mindlin shell and break the Kirchhoff-Love benchmark limit. The deformed mesh is shown in Figure 10. The ra- dial component of the deflection along the directrix from the benchmark solution [12] in Figure 11 com- pares well to the deformed mesh in Figure 10, note 50 100 150 200 250 300 -300 -250 -200 -150 -100 -50 0 0 100 200 300 Figure 10. Deformed mesh, 48 elements along the directrix Figure 11. Radial component of deflection along the directrix, benchmark by the solid line, [12] the shalow inward deflection in the lower part of the front directrix. The precision and convergence of the deflection beneath the pin force is, as expected, worse than for elements based on the Kirchhoff-Love assumptions. A solution of the Ponikla bridge arch loaded by characteristic arch selfweight and standard LM1 tan- dem axle with wheel forces of 250 kN is provided for the illustration of the code output in Figure 12. Dis- placements are 2000 times scaled up. In comparison to the 2D solution of the same bridge arch in Section 3, the interaction and selfweight of the fill and pavement are not accounted for, the tandem axle load and mesh density are slightly different, too. The tandem axle position in the transverse direction is the extreme eccentric one within the bounds defined by the EN 1991-2. The differences between the loaded/unloaded sides of the arch barrel are 1261/405 kPa in extreme normal stresses in the bed joints and 0.41/0.25 in relative crack depths. The backfill and pavement stiff- ness and selfweight reduce the differences in the real bridge, but the example testifies that the 2D models need corrections. The 3D bridge model analogue to the one in Section 3 is currently being worked on. The absence of rotational DOFs improves the con- vergence of the iterations. The no-tension, history independent material model admits a single load in- crement strategy. In the illustration example, the 301 Petr Řeřicha Acta Polytechnica 0 2 4 6 8 0 -10.5 -10 -9.5 2 4 6 8 10 -1758||0.34-1647||0.3 -963.8||0 -1222||0.36 -1493||0.36 -1368||0.36 -795.9||0 -1594||0.35 -862.4||0.35-1046||0.36 -593.7||0 -680.6||0 368.1||0 -1663||0.35-1727||0.35 -921.3||0 -899.1||0 -884.5||0 -832.6||0 -891.1||0 893.6||0 924||0 1036||0.21 1139||0.27637.6||0 740.3||0 678.6||0 726.5||0.16 693.5||0 993.4||0.2 952.8||0.26 952.8||0.26 875.1||0 1071||0.221139||0.27 1087||0.2 1141||0.35 1141||0.33 1186||0.33 -499.6||0 333.4||0 370.9||0 727.8||0.15 891.2||0.3 849.8||0.33 599||0.24847.9||0.33 889.3||0.31 594.3||0.23 605.2||0.099 608||0.27 495.5||0.27608.1||0.27 610||0.1 -622.6||0.33 -466.7||0.36 -346.7||0.4 -527.3||0.35 -331.1||0.087 -719.4||0.34 -459.8||0 -301.5||0 -301.1||0 349.1||0 342.1||0 159||0 336.4||0 -342||0.13 -218.1||0 162.3||0 189.8||0 329.3||0.096 324.7||0.09 495.6||0.27 459.9||0.2 405.9||0.25 329.8||0.17405.8||0.25 459.9||0.2 459.3||0.015 461.2||0.021 351.8||0.21 358.1||0.11 362||0 362.5||0 358.1||0.11 291.1||0.039 286.7||0 287||0 220.9||0 1261||0.41 1066||0.34 732||0.048 607||0 1068||0.34 773.3||0.096607.5||0 781.3||0.1 471.9||0 529.7||0 529.4||0 413||0 413.8||0 326.9||0 467||0 467.6||0 364.7||0 364.7||0 -426.6||0 -497.2||0 -557.2||0 -585||0 -625||0 -609||0-409.8||0 336.9||0 -409.8||0 337||0 -497.2||0 -290.2||0 -420.3||0 -420.3||0 -520.3||0 -585||0 -520.3||0 -573||0 -499||0 -399.5||0 -498.9||0 -462.6||0 -573.2||0 -528.8||0 -529.8||0 -462.6||0 -290.1||0 -272.1||0 -399.5||0 -367.1||0 -609||0 -560.4||0 -560.5||0 -492.6||0 -415.9||0 -481||0 -493.3||0 -480.3||0 -415.9||0 -412.4||0 -373.3||0 -416.4||0 -416||0 -411.6||0 -371.7||0 -391.1||0 327.4||0 264.6||0 -271.6||0 -252.6||0 -366.9||0 264.8||0 209.4||0 209.4||0 -252.2||0 -331.3||0 -157.2||0 -243.2||0 -331.3||0 -242.7||0 -237||0 -305.4||0 -308.9||0 -373||0 -283.4||0 -356.3||0 -344.6||0 -329.8||0 -369.7||0 -339.5||0 -325.6||0-312.6||0 -321.4||0 -316.1||0 -277.8||0 -283.5||0 322.6||0 301||0 277.8||0 -275.3||0 -279.3||0 -630.6||0 -566.2||0 -566.2||0-469.2||0 -470.3||0 -569.9||0-451.3||0 -452.3||0 -428.4||0 425.4||0 456.5||0 447.1||0 420.9||0 -395.3||0 -398.2||0 411.7||0 393.8||0 366.5||0 345.6||0 Figure 12. Bridge at Ponikla, 3D shell, compression stresses [kPa] and relative depths of the cracks in the bed joints ratio 0.005 of the RMS norms of the imbalance and load was reached in 3 iterations. 5. Conclusions The simplest finite elements for 2D and 3D thick shells have been developed for applications in masonry arch bridges limit load analysis. They feature seamless com- bination with 2D CST and 3D tetrahedron elements for the bridge’s backfill. Exclusively translational DOFs are used, which improves the convergence of iterations. The no-tension condition is applied to the normal stress in the bed joints planes. The output minimum normal stresses and cracks depths in these planes can be used to assess the ultimate limit loads of a bridge. A sample analysis of an arch of the bridge at Poniklá testifies that this material law and shell elements develop a characteristic failure mode of the stone masonry arch bridges, the gradual creation of the virtual hinges. Acknowledgements The author was supported by the grant NAKI DG20P02OVV001, ‘Tools for the preservation of historical values and functions of arch and vaulted bridges’ provided by the Ministry of Culture of the Czech Republic. References [1] The Highways Agency, London. Design manual for roads and bridges, vol. 3, Highway structures: inspection and maintenance, section 4, Assessment, Part 3, The assessment of highway bridges and structures, 2001. [2] Ministry of transport of the Czech Republic. Zatížitelnost zděných klenbových mostů - Load rating of masonry arch bridges, 2008. [3] T. E. Ford, C. E. Augarde, S. S. Tuxford. Modelling masonry arch bridges using commercial finite element software. In Proceeding of the 9th International Conference on Civil and Structural Engineering Computing. Civil-Comp Press, Stirling, 2001. [4] Mott-MacDonald, 20-26 Wellesley road, Croydon, CR9 2UL, UK. CTAP manual for the assessment of masonry arch bridges, 1990. [5] R. Bridle, T. Hughes. An energy method for arch bridge analysis. Proceedings of the Institution of Civil Engineers 89:375–385, 1990. https://doi.org/10.1680/iicep.1990.9397. [6] M. Gilbert. RING: a 2D rigid block analysis program for masonry arch bridges. In Proceedings of the 3rd International Arch Bridges Conference, pp. 459–464. 2001. [7] W. Harvey. Application of the mechanism analysis to masonry arches. Structural engineer 66:77–84, 1988. [8] P.-O. Persson, G. Strang. Simple mesh generator in MATLAB. SIAM Review 46(2):329–345. https://doi.org/10.1137/S0036144503429121. [9] R. D. Mindlin. Influence of rotatory inertia and shear on flexural motions of isotropic, elastic plates. ASME Journal of Applied Mechanics 18(1):31–38, 1951. https://doi.org/10.1115/1.4010217. [10] A. Pipard, R. Ashby. An experimental study of a voussoir arch. Proceedings of the Institution of Civil Engineers 10:383, 1938. [11] A. J. S. Pippard. The civil engineer in war, vol. 1, chap. The Approximate Estimation of Safe Loads on Masonry Bridges., pp. 365–372. 1948. https://doi.org/10.1680/ciwv1.45170.0021. [12] G. Lindberg, M. D. Olson, G. R. Cowper. New developments in the finite element analysis of shells, 1969. https://apps.dtic.mil/sti/pdfs/AD0707780.pdf. [13] W. Flugge. Stresses in shells. Springer Verlag, 1962. 302 https://doi.org/10.1680/iicep.1990.9397 https://doi.org/10.1137/S0036144503429121 https://doi.org/10.1115/1.4010217 https://doi.org/10.1680/ciwv1.45170.0021 https://apps.dtic.mil/sti/pdfs/AD0707780.pdf Acta Polytechnica 62(2):293–302, 2022 1 Introduction 2 The TT element, Timoshenko beam with translational DOFs 2.1 Constitutive equations for bending 2.2 Full cross-section stiffness matrix 2.3 Arch nodes equilibrium 3 Application to an arch bridge 4 3D Mindlin shell element TM 4.1 Bending and in-plane deformation 4.2 Shear owing to the rotations of normals 4.3 TM flat slab element 4.4 TM shell element 4.5 TM element code and test 5 Conclusions Acknowledgements References