398 Int. J. Masonry Research and Innovation, Vol. 1, No.

4, 2016

Horizontal reaction components of pointed vaults

Gábor Lengyel and Katalin Bagi*

Department of Structural Mechanics,
Budapest University of Technology and Economics,
1111 Budapest, Műegyetem rkp. 3. K.mf.63., Hungary
Email: [email protected]
*Corresponding author

Abstract: The magnitude of horizontal reactions under self-weight are

analysed for pointed barrel vaults with different eccentricities. The magnitude
is expressed as the function of the eccentricity (deviation of the pointed
generator curve from the circle). Vaults with fixed supports and with increasing
outwards support displacements are considered. An analytical model based on
continuous Timoshenko beam theory is suggested for barrel vaults and then
checked with computer-simulated experiments. The results show that with
increasing support separation, the horizontal reaction tends to be at a constant
value when reaching the statically determinate state. After this, the horizontal
reaction varies only because of the gradual modification of the geometry due to
large support displacements.

Keywords: 3DEC; geometry; masonry; pointed arc; support displacement.

Reference to this paper should be made as follows: Lengyel, G. and Bagi, K.
(2016) ‘Horizontal reaction components of pointed vaults’, Int. J. Masonry
Research and Innovation, Vol. 1, No. 4, pp.398–420.

Biographical notes: G. Lengyel is a PhD student at TU Budapest. His field of

study is the mechanics of masonry vaults.
K. Bagi is a full professor in Structural Mechanics at TU Budapest. She
specialises in the discrete element modelling of granular assemblies and
masonry structures.
This paper is a revised and expanded version of a paper entitled ‘Analysis of
the horizontal reaction of pointed barrel and cross vaults’, presented at the
8th International Congress of the Croatian Society of Mechanics, 29 Sept–
02 Oct 2015, Opatija, Croatia.

1 Introduction

1.1 The analysed problem

Arches and vaults based on pointed generator curves are essential components of gothic
architecture. Pointed structural shapes have preliminaries in Islamic architecture, and the
idea may have been brought to the Christian Europe from Moorish Spain or from the near
East by the crusaders. The beauty of pointed shapes surely significantly contributed to the
increasing popularity of this idea in the cathedral art of the Christian world from the 12th
century, but without clear mechanical advantages such a rapid and widespread

Copyright © 2016 Inderscience Enterprises Ltd.

 Horizontal reaction components of pointed vaults 399

assimilation of this style would not have been possible. Pointed structures are known to
express smaller lateral thrusts to their supports than circular vaults, which means that
arches and vaults could be built higher and more slender, larger distances could be
covered with shapes pointing towards the sky, and more light could be allowed to enter
through the large windows than ever in the constructions of the previous centuries.
Pointed curves (which will serve below as generators for the vault shapes to be
analysed) can most easily be produced with the help of circular arcs in the way shown in
Figure 1: two arcs, having their radii R larger than half of the span L, and their centres O
being located eccentrically (i.e. not coinciding with the middle of the span), intersect in a
symmetrical manner. The larger is the radius of the arcs, the larger is the eccentricity (i.e.
the deviation of the arc centre from the middle point of the span, quantified with the help
of the scalar α) and the higher is the resulting pointed curve. Note that R = (0.5 + α) L,
which means that α = 0 corresponds to the semicircular shape and α = 0.5 belongs to the
pointed arch based on the equilateral triangle. As long as the arc centre remains on the
horizontal base line, the tangent is vertical at the two supports, but the generator curves
become higher with increasing eccentricity.

Figure 1 The analysed pointed curves and the definition of the eccentricity parameter α (see
online version for colours)

This study focusses on the magnitude of lateral thrusts expressed on the supporting
structural members (walls and pillars) by pointed barrel vaults.
In the present study analytical investigations and discrete element simulations are
performed, to assess how the magnitude of eccentricity influences the magnitude of
horizontal thrust (normalised with the vertical reaction component, which is equal to the
weight of the half of the structure). The analysed barrel vault is shown in Figure 2. An
analytical model based on the Timoshenko beam theory is developed, and the results are
compared with DEM-simulated experimental results.

Figure 2 The analysed structure

400 G. Lengyel and K. Bagi

The geometry of a vault is never perfect in reality. In addition to inevitable imperfections

of the construction procedure, a primary reason is that displacements of the supports
(particularly those caused by the forces expressed by the vault itself) lead to gradually
increasing modifications of the vault shape and position. The dominant support
displacement is usually the outwards translation of the underlying structural elements. To
analyse this effect, during the computations, the supports of the simulated vaults were
gradually moved outwards in a quasi-static manner, and the horizontal reactions were
recorded throughout the process. An important objective of the investigations was to see
how the reactions varied before and after reaching the statically determinate state.

1.2 The applied methods

The present paper compares analytical results on the horizontal reactions to the outcomes
of numerical simulations.
The two most important classes of analytical approaches to study the forces in
masonry arches are limit state analysis and deformable continuum models.
Being inspired by the classical works of Heyman (1966, 1969), which were based on
the idea of Kooharian (1952), several researchers applied limit state analysis to develop
different techniques for the analysis of masonry arches. Following this approach,
Romano and Ochsendorf (2010) focussed on pointed masonry arches and analysed the
effect of the deviation from circular shape, regarding minimum thickness, collapse loads,
and the admissible domain of horizontal thrust. The authors determined the statically
admissible range of the thrust through graphic statics, according to the static theorem of
limit state analysis. Their results will be compared below to the outcomes of the present
investigations.
Alexakis and Makris (2013) analysed the minimum thickness of elliptical arches, with
the help of Limit State Analysis. A most important finding was that the minimum
thickness depends on the stereotomy.
The investigations of De Rosa and Galizia (2007) were also based on the static
theorem. Arches of different eccentricities (from semicircular to pointed arch based on
the equilateral triangle) were analysed under self-weight plus live load, and the highest
possible live load multipliers were searched for. Their method took into consideration the
frictional limit and the compression strength of the voussoirs (‘yield domains’ described
the admissible relations between the different components of the contact forces and
moment). Solving the optimisation problem to maximise the load multiplier while
satisfying the equilibrium conditions as well as the restrictions for contact forces and
moments, it was found that the equilateral arch showed the best resistance to live loads. It
was also pointed out that sliding occurred only in case of small frictional resistance at the
contacts (below 17−33°, depending on the arch shape).
The kinematic theorem of Limit State Analysis was applied by Coccia et al. (2015)
for circular arches. The investigations pointed out that the arch thickness significantly
influenced the support displacements and the horisontal reactions belonging to the
collapse state.
The application of deformable continuum models for masonry is perhaps less
accepted in general, though the widely applied MEXE method, for example, is based on
the elastic beam analysis of Pippard (1948). It was emphasised by Huerta (2001) citing
the studies of Heyman that the elastic solution (and particularly the internal force system)
is rather sensitive to small disturbances of the geometry, for example, to support
 Horizontal reaction components of pointed vaults 401

displacements. In addition, if linearly elastic material behaviour is applied in the model,

significant deviations can be expected from the real structural behaviour where the
masonry can hardly resist tension. In spite of these doubts, a few authors took the
challenge and applied continuous beam theories for masonry arches. The papers of
Audenaert, Peremans and De Wilde (2004) and Audenaert, Peremans and Reniers (2007)
are excellent examples for this approach: the authors developed a linearly elastic,
perfectly plastic, no-tension model in the frame of a classical beam bending approach.
The discrete built-up of the masonry and the continuum-mechanical modelling
approach are not necessarily in unresolvable contradiction. Finite element models can
also be applied to study masonry arch behaviour, and the discrete nature of the masonry
can be represented in different ways. Individual voussoirs and the mortar layers between
them can be modelled, for instance, with 3D finite elements having different material
properties for the mortar and for the blocks where contact elements represent the mortar-
block interface. Hejazi and Jafari (2010) applied FEM this way to compare different
bricklaying techniques and analysed the stress distributions in the mortar and in the
blocks. Alternatively, contact elements with suitable constitutive behaviour can also
represent the joints. When such a detailed modelling is not necessary, homogenised
continua can be defined, either as 1D beam models, or 2D shell or 3D body elements.
The seeming contradiction between limit state analysis with rigid blocks and micro-
level continuum models where the blocks and the joints have finite resistance to different
stress components is partly resolved by the approach proposed by Foce and Sinopoli
(2001) and Foce and Aita (2003). The graphostatic method originally introduced by
Durand-Claye (1867) searches for a range of admissible contact force systems that satisfy
both the equilibrium conditions and the restrictions for compressive, tensional, and
frictional stresses in the joints. Using the method, Aita et al. (2004) investigated the
collapse modes (rotational and sliding) of pointed arches under self-weight, for different
values of thickness and friction coefficients. The mechanical state (i.e. elastic, cracked for
tension, crushed for compression) of points of arches having different circular, pointed,
or horizontally elongated ellipsoidal shapes were analysed in Aita, Barsotti and Bennati
(2009) and their load bearing capacity was compared. The higher load bearing of pointed
arches was clearly shown. The arch-and-pier system was analysed recently in Aita,
Barsotti and Bennatti (2015).
Continuous beam models are usually thought to be unreliable for masonry arches. To
query this belief, in the present paper a continuous beam model (inspired by Oh, Lee and
Lee, 1999) will be developed, and the range of validity of its predictions will be checked
with DEM-simulated experiments.
The computer simulations in the present study were done with 3DEC, a commercial
discrete element code. The discrete element method (DEM) considers the structure to be
a collection of separate blocks, ‘discrete elements’, each of which is able to move and - in
most software - to deform independently of each other. The blocks may come into
contact with each other; hence, distributed forces can be transmitted from one block to
another, causing stresses and deformations in the blocks. They can slide along each other,
be separated, or the contacts may be partially cracked. Because of its ability to follow the
characteristic failure modes and collapse histories of masonry structures in detail, DEM
can serve as a tool to check theoretical results with the help of computer-simulated
experiments. A characteristic application of DEM can be found in Rizzi, Rusconi and
Cocchetti (2014): a discrete element code was used for the comparison of the theoretical
402 G. Lengyel and K. Bagi

predictions for circular arches by Milankovitch (1907), Heyman (1969), and Cochetti,
Colasante and Rizzi (2011).
In that version of the 3DEC code that is applied in the present study the discrete
elements are made deformable by being subdivided into uniform-strain tetrahedra whose
nodal translations are the basic unknowns of the analysis. The discrete elements in 3DEC
may have any polyhedral shape. The displacements of the nodes of the tetrahedra are to
be calculated from Newton’s law of force and acceleration. In the equations of motion the
mass of a gridpoint is defined as the volume of the Voronoi-cell around that gridpoint
multiplied with the material density of the element. Different forces acting on the
Voronoi-cell like self-weight, drag force, distributed forces from neighbouring Voronoi
cells inside the same discrete element, and contact forces expressed by the neighbouring
discrete elements are reduced to the gridpoint.
The joints between the elements may have, in principle, two different roles in DEM,
depending on the problem to be modelled. The first option is that the joints represent
some kind of mortar layer having a finite thickness in reality. In this case, the material
parameters of the joints express the deformability of these layers and their failure criteria.
Therefore, a normal and a tangential stiffness (i.e. resistance to relative translation) have
to be prescribed along with fracture criteria (normal and shear strength) and perhaps a
friction coefficient in order to define the conditions when the joints fail.
The other option (the one used in the present study) is most suitable if there is no
material layer modelled between the blocks in the real system. In this case the aim is to
simulate dry contacts or perhaps to neglect the weak mortar in an old masonry structure.
The friction coefficient of the joints still has a real physical meaning, expressing the
sliding criterion of the two contacting bricks or stone blocks along each other. A normal
and a shear stiffness express the resistance to relative displacements of contacting
material points in the contacts: the increments of the forces transmitted from one block to
the other are the product of the increment of the relative normal or tangential translation
of the two material points and of the normal or shear stiffness. These contact stiffnesses
are, in this case, artificial numerical penalty parameters only.
The equations of motion are solved (the displacement increments of the nodes during
consecutive finite Δt time intervals are determined) using a central difference scheme for
numerical time integration. From the displacement increments done during a time step
the strain increments, changes of stresses and forces, and so on are then calculated: the
simulation of the state-changing is achieved by step-by-step calculation of the
incremental motions of the gridpoints.
Technical details of the applied 3DEC code are described in several papers, books
and manuals, for example, in Lemos (2016), Lengyel and Bagi (2015), or Szakaly,
Hortobagyi and Bagi (2016).
The paper is organised as follows. Section 2 focusses on deriving theoretical
predictions for the horizontal reaction of barrel vaults, treating them as continuous, two-
dimensional structures (so a unit slice of the vault is considered), assuming that the
analysed slice is in plane stress state. The derivations are then checked with the help of
3DEC simulations in Section 3, where the plane stress and the plane strain behaviour are
also compared. Finally, in Section 4 the most important conclusions are drawn.
 Horizontal reaction components of pointed vaults 403

2 Theoretical analysis of the horizontal reaction of pointed barrel vaults

2.1 Basic assumptions and notations

The middle surface of the analysed pointed barrel vault is cylindrical, produced by two
intersecting circular arcs that serve as the generator curve of the middle surface
(Figure 3a). Consider a unit thickness slice of the barrel vault with thickness t and
material density ρ; this slice is modelled as a Timoshenko beam with curved axis
(Henrych, 1981; Borg, 1959). The calculations are based on the following assumptions:
• The material of the arch is linearly elastic, and it resists arbitrary tension as well as
compression and shear. (This assumption contradicts real masonry behaviour so its
consequences will be analysed with the help of comparing the theoretical predictions
to discrete element simulation results.)
• The translations and rotations are small. (Analysed support displacements will not
exceed 0.03% of the span, which corresponds to 3 mm for a 10 m span.)
• The arch is represented with its axis defined as the curve in the half of the thickness.
The basis of the derivations will be the behaviour of an infinitesimally short arch
element shown in Figure 3b, whose length ds is measured along the axis. (The mass
centre of the arch element is located slightly outside the axis, but this deviation will
be neglected.)
• According to the Timoshenko beam theory, the material points that form planar
cross-sections being originally perpendicular to the axis of the beam will remain
planar, though not necessarily perpendicular to the deformed axis.

Figure 3 (a) The geometrical parameters of the analysed arch; (b) the arch element and the forces
acting on it (see online version for colours)

The position of the beginning cross-section of the arch element is described by the angle
φ. A local coordinate system (ξ, η) is assigned to its centroid (Figure 3b). This local
coordinate system rotates while proceeding along the axis in such a way that the ξ and η
axes always point in the tangential and radial directions respectively, so the static and
404 G. Lengyel and K. Bagi

kinematical variables will be written in terms of components according to directions

varying along the axis.

2.2 Statical equations

Every arch element is in equilibrium under its weight and the internal forces acting at its
beginning and ending cross-sections. Figure 3b shows the positive directions of N(φ),
V(φ), and M(φ), the normal force (negative for compression), the shear force, and the
bending moment at φ, respectively The equilibrium of the arch element can be expressed
by the following three equations:

∑ F ξ : − N + ( N + dN ) cos(dφ ) + (V + dV )sin(dφ ) − ρ gtds cos(φ ) = 0

i

∑ Fη : − V + (V + dV ) cos(dφ ) + ( N + dN ) sin(dφ ) − ρ gtds sin(φ ) = 0

i (1)
∑ M ζ : + M − ( M + dM ) + NR − ( N + dN ) R + ρ gtds cos(φ ) R = 0
i
O

Since the angle dφ is small (cos (dφ) ≈ 1 and sin (dφ) ≈ dφ), after some rearrangements
the three equations can be written as
dN
+ V = R ρ gt cos(φ )
dφ
dV
−N + = R ρ gt sin(φ ) (2)
dφ
dN 1 dM
− − = − R ρ gt cos(φ )
dφ R dφ

These equations are, in principle, identical to those in Audenaert, Peremans and De Wilde
(2004) and Audenaert, Peremans and Reniers (2007).

2.3 Kinematical equations

Three displacement variables will be used to characterise the kinematics of a cross-
section at φ (Figure 3b). The translations u(φ) and v(φ) denote the tangential and radial
translation components of the centroid of the cross-section, while θ(φ) expresses the
rotation of the material region forming the cross-section that is perpendicular to the axis
in the initial state.
In case of the Timoshenko theory, the rotation θ(φ) of the planar material region
forming the cross-section is not equal to the rotation of the axis, ψ(φ). Instead, θ(φ) is the
combined result of two effects:
1 As the planar set of material points forming the cross-section rotates while the arch
element has zero shear deformation, the rotation of the material region is equal to
ψ(φ).
2 The shear deformation γ(φ) of the arch element (without any rotation of the axis)
leads to the rotation of the material region equal to −γ(φ).
The rotation of the material region forming the cross-section is θ(φ) = ψ(φ) − γ(φ). The
rotation of the axis, that is, ψ(φ) = θ(φ) + γ(φ), can alternatively be expressed in terms of
 Horizontal reaction components of pointed vaults 405

the translation components of the centroids at φ and at φ + dφ , taking into account that
the angle dφ is small and that ds = R ⋅ dφ (Figure 4):
(v + dv) ⋅ cos(dφ ) − (u + du ) ⋅ sin(dφ ) − v dv u ⋅ dφ
ψ= = − (3)
ds ds ds
or equivalently:
1 ⎛ dv ⎞
ψ = ⎜ −u⎟ (4)
R ⎝ dφ ⎠

Figure 4 Calculation of the rotation of the axis in a curved element: (a) the translations of the
axis; (b) the rotation components of the cross-section (see online version for colours)

From the kinematic variables above, the deformations of the arch element, that is, the
elongation (ε), the shear deformation (γ), and the curvature (κ), have to be derived with
the help of the derivatives with respect to the length parameter of the arch element, ds.
These deformations of the arch element will be expressed now in terms of the relative
displacements of the cross-sections at φ and at φ + dφ.
First, the elongation of the axis is
(u + du ) ⋅ cos(dφ ) − (v + dv) ⋅ sin(dφ ) − u du v ⋅ dφ
ε= = + (5a)
ds ds ds
or
1 ⎛ du ⎞
ε= ⎜ + v⎟ (5b)
R ⎝ dφ ⎠
Second, the shear deformation of the arch element, γ(φ) = ψ(φ) − θ(φ), can also be
expressed in terms of the basic kinematical variables:
1 ⎛ dv ⎞
γ =ψ −θ = ⎜ − u − Rθ ⎟ (6)
R ⎝ dφ ⎠
Finally, the curvature of the cross-sections at φ and at φ+dφ is expressed as follows:
406 G. Lengyel and K. Bagi

1 dθ
κ =− (7)
R dφ

The constitutive equations will relate these three characteristics to the internal forces.

2.4 Constitutive equations

According to Borg and Gennaro (1959), in a curved beam with small thickness the
internal forces N(φ), V(φ), and M(φ) can be related to the deformations ε (φ), γ (φ), and
κ (φ) in the following way:
EI
N = EA ⋅ ε + κ ; V = kGA ⋅ γ ; M = EI ⋅ κ . (8)
R
Here E and G denote the Young's modulus and the shear modulus of the material of the
arch. The Timoshenko shear coefficient k is chosen to 5/6 because of the arch having a
rectangular cross-section. The geometrical characteristics I and A (inertia and area of the
cross-section) are assumed to be constant along the arch, which means that the possibility
of partial cracking or complete hinging of the cross-sections is excluded from the
analysis. (Note that the Bresse equations applied in Audenaert, Peremans and Reniers
(2007) do not contain the effect of shear deformations, but on the other hand, in those
papers a sophisticated elasto-plastic no-tension model is used for studying the behaviour
of the cross-sections.)
In a statically indeterminate structure the results may, in general, strongly be
influenced by the actual material properties. In every calculation throughout the paper the
Young’s modulus E = 14.8 GPa, Poisson’s ratio ν = 0.20, and material density
ρ = 2348 kg/m3 (travertine limestone) were applied. In addition, a sensitivity test will also
be presented below in Section 2.6 to see how the material properties influence the statical
behaviour.

2.5 Boundary conditions

The mechanical model of the half of the pointed arch is shown in Figure 5, and the
boundary conditions are listed in Table 1 (dy and Rz denote horizontal outwards
translation and vertical internal force component, respectively).

Figure 5 Boundaries of the half arch (see online version for colours)
 Horizontal reaction components of pointed vaults 407

Table 1 Boundary conditions of the half arch

Boundary u V θ dy Ry
Γ0b – d 0 – –
Γsym. 0 – 0 0 0

Because of the symmetry of the problem, at the top of the arch the resultant of N(φmax)
and V(φmax) should be horizontal; a statical boundary condition is received this way. Five
kinematic boundary conditions can also be specified. The bottom left cross-section of the
arch remains horizontal and does not move vertically; its horizontal translation is zero or
a prescribed value d. The top cross-section neither rotate nor translate horizontally,
because of the symmetry of the problem again. (Note that these conditions exclude the
occurrence of cracks at the support or the crown in the masonry.) The six boundary
conditions can be written in terms of the statical and kinematical variables in the
following way:

u φ =0 = 0, u φ =φmax sin(φmax ) − v φ = φmax cos(φmax ) = 0

v φ =0 = d, V φ =φmax sin(φmax ) + N φ = φmax cos(φmax ) = 0 (9)
θ φ =0 = 0, θ φ =φmax =0

2.6 Solution and the results

After substituting the kinematical equations (5b), (6), and (7) into the constitutive
equation (8), three first-order differential equations are received, while the statical
equation (2) also constitutes three first-order differential equations. A system of six first-
order differential equations is received, where the unknowns are the displacements,
rotations, and internal forces, Equation (10):
⎡ N '⎤ ⎡ 0 −1 0 0 0 0 ⎤ ⎡ N ⎤ ⎡ ρ gRt cos(φ ) ⎤
⎢V ' ⎥ ⎢ 1 0 0 0 0 0 ⎥⎥ ⎢⎢ V ⎥⎥ ⎢⎢ ρ gRt sin(φ ) ⎥⎥
⎢ ⎥ ⎢
⎢ M '⎥ ⎢ 0 R 0 0 0 0 ⎥ ⎢M ⎥ ⎢ 0 ⎥
⎢ ⎥=⎢ ⎥⎢ ⎥+⎢ ⎥
⎢ u ' ⎥ ⎢ NR / EA 0 M / EA 0 −1 0 ⎥ ⎢ u ⎥ ⎢ 0 ⎥ (10)
⎢ v' ⎥ ⎢ 0 R / kGA 0 1 0 R⎥ ⎢ v ⎥ ⎢ 0 ⎥
⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥
⎣⎢ θ ' ⎦⎥ ⎣⎢14444444
{
0 0 − R / EI 0 0 0 ⎦⎥ ⎣⎢ θ ⎦⎥ ⎣⎢
4244444444 3 { 144244
0
3⎦
⎥
x' (φ ) A x (φ ) b (φ )

This system can be written in the form x '(φ ) = Ax(φ ) + b(φ ). The complete solution of
this system consists of a homogeneous and a particular part: x(φ ) = xh (φ ) + xp (φ ).
The homogeneous system has six linearly independent solutions (x i (φ ), i = 1K 6).
A 6 × 6 matrix can be defined that consists of these six solution vectors:
X(φ ) = [ x1 (φ ) x 2 (φ ) x 3 (φ ) x 4 (φ ) x 5 (φ ) x6 (φ ) ] (11)

(Note that bold capitals denote matrices, so that X is a matrix whose columns are the six
independent xi vectors.) The solution is the linear combination of the xi(φ) vectors, so
408 G. Lengyel and K. Bagi

xh (φ ) = X(φ )c . This system can be solved according to Arnold (1973) and Tufekci and
Arpaci (1998). The eigenvalues of the matrix A are as follows:
λ1,2 = −i; λ3,4 = 0; λ5,6 = i (i is the imaginary unit), so three duplicities in the
eigenvalues are received. The unknown matrix X(φ) can be written in the following form:
X(φ ) = e − iφ H 1 ( A) + φ e − iφ H 2 ( A) + H 3 ( A) + φ H 4 ( A) + eiφ H 5 ( A) + φ eiφ H 6 ( A) (12)

where H1, H2, … H6 are the Hermitian matrix polynomials of matrix A and can be
determined in the following way. Introduce the auxiliary complex variable, z. Every Hj
should satisfy H j ( z ) = α j 5 z + α j 4 z + α j 3 z + α j 2 z + α j1 z + α j 0 = z α j , where
5 4 3 2 T

T T
z = ⎡⎣ z 5 z4 z3 z 1⎤⎦ , and α j = ⎡⎣α j 5 α j 4 α j 3 α j 2 α j1 α j 0 ⎤⎦ . The
z2
coefficients αj0, αj1, … αj5 belonging to Hj are searched for, separately for every j. These
coefficients can be calculated from the linear equation system, Eq. (13), where the actual
values z = −i, z = 0, and z = i are the eigenvalues of matrix A:

⎡ z T |z =− i ⎤
⎢ ⎥
⎢ d zT | ⎥ ⎡ −i 1 i −1 −i 1⎤
⎢ dz z =− i
⎥ ⎢ 5 4i −3 −2i 1 0 ⎥⎥
⎢ zT | ⎥ ⎢
⎢ z =0 ⎥ ⎢0 0 0 0 0 1⎥
⎢ d T ⎥α j = ⎢ ⎥α j = e j (13)
⎢ z |z = 0 ⎥ ⎢0 0 0 0 1 0⎥
⎢ dz T ⎥ ⎢i 1 −i −1 i 1⎥
⎢ z |z = i ⎥ ⎢ ⎥
⎢ d ⎥ ⎣⎢ 5 −4i −3 2i 1 0 ⎦⎥
⎢ z T |z = i ⎥
⎣ dz ⎦

On the right side, the j-th unit vector can be seen for H1 , e1 = [1 0 0 0 0 0] ,
H 2 , e2 = [ 0 1 0 0 0 0] , and so on.
To derive the particular part, the method of undetermined coefficients can be applied,
) )
where xp (φ ) = X(φ )c (φ ) . The vector c (φ ) is a function of position along the arch, ϕ. By
substituting this expression into x '(φ ) = Ax(φ ) + b(φ ), the following relation is received:
) ) )
X '(φ )c (φ ) + X(φ )c '(φ ) = AX(φ )c (φ ) + b(φ ) (14)
) )
But matrix X(φ ) solves the homogeneous equation, so X '(φ )c (φ ) = AX(φ )c (φ ).
)
Consequently, c '(φ ) = X −1 (φ )b(φ ), which yields the particular solution:

xp (φ ) = X(φ ) ∫ X −1(φ )b(φ ) dφ (15)

Finally, the complete solution is the sum of the homogeneous and the particular
part: x(φ ) = xh (φ ) + xp (φ )
As mentioned above the homogeneous solution is the linear combination of the
solution vectors. The unknown vector c (remember that xh (φ ) = X(φ )c ) has six scalar
components, which have to be determined with the help of the six boundary conditions,
see Eq. (9), as follows:
 Horizontal reaction components of pointed vaults 409

The six boundary conditions set requirements for the five functions u(φ), v(φ), θ(φ),
N(φ), and V(φ). These are components of x(φ) so they can be written as:
6
u (φ ) = ∑ ci uhi (φ ) + u p (φ )
i =1
6
v(φ ) = ∑ ci vhi (φ ) + v p (φ )
i =1
6
ϑ (φ ) = ∑ ciϑhi (φ ) + ϑ p (φ )
i =1
6
N (φ ) = ∑ ci N hi (φ ) + N p (φ )
i =1
6
V (φ ) = ∑ ciVhi (φ ) + V p (φ )
i =1

(Note that the scalar function M(φ) is not listed here, since no boundary condition is
specified for the bending moment.). Since the x1(φ), x2(φ), … x6(φ) components of the
homogeneous part and the xp(φ) particular part are already determined, the only
unknowns on the right side are the c1, c2, … c6 coefficients.
The six boundary conditions can be collected into a linear equation system
M⋅c +b = 0 where:

⎡ uh1 (0) uh2 (0) uh3 (0) uh4 (0) uh5 (0) uh6 (0)
⎤
⎢ vh1 (0) vh2 (0) vh3 (0) vh4 (0) vh5 (0) vh6 (0)
⎥
⎢ u (φ ) sin φ −
θh1 (0) θh2 (0) θh3 (0) θh4 (0) θh5 (0) θh6 (0)
⎥
M=⎢ (
h1 max max
⎢ vh1(φmax ) cos φmax
)( vh2 (φmax ) cos φmax
)(
uh2 (φmax ) sin φmax −
vh3 (φmax ) cos φmax
)(
uh3 (φmax ) sin φmax −
)(
uh4 (φmax ) sin φmax −
vh4 (φmax ) cos φmax
)(
uh5 (φmax ) sin φmax −
vh5 (φmax ) cos φmax
)
uh6 (φmax ) sin φmax −
vh6 (φmax ) cos φmax
⎥
⎥
(
⎢ Vh1(φmax ) sin φmax +
⎢ Nh1 (φmax)cos φmax
)( Nh2 (φmax )cos φmax
)(
Vh2 (φmax ) sin φmax +
)(
Vh3 (φmax ) sin φmax +
Nh3 (φmax )cos φmax
)(
Vh4 (φmax ) sin φmax +
Nh4 (φmax)cos φmax
)(
Vh5 (φmax ) sin φmax +
Nh5 (φmax)cos φmax
)
Vh6 (φmax) sin φmax +
Nh6 (φmax )cos φmax ⎥
⎥
⎣ θh1(φmax ) θh2 (φmax ) θh3 (φmax ) θh4 (φmax ) θh5 (φmax ) θh6 (φmax ) ⎦

and
⎡ u p (0) ⎤
⎢ v p (0) ⎥
⎢ ⎥
⎢ θ p (0) ⎥
b=⎢ ⎥
⎢ u p (φmax ) sin φmax − v p (φmax ) cos φmax ⎥
⎢V p (φmax ) sin φmax + N p (φmax )cos φmax ⎥
⎢ ⎥
⎢⎣ θ p (φmax ) ⎥⎦

Now the unknown c can easily be determined from c = − M−1 b. The results will be
visualised and illustrated in the rest of Section 2.6.
Now the results on two statical characteristics will be shown here. The first quantity
(denoted by δ, named the normalised horizontal reaction) is the magnitude of the
horizontal force component (Figure 6), divided with the magnitude of the vertical
component (i.e. with the weight of the half arch, shown in Figure 5):
410 G. Lengyel and K. Bagi

RH
δ := .
RV

Figure 6 Reaction components acting on the bottom boundary (see online version for colours)

This quantity was determined for different eccentricities, thicknesses, and for different
magnitudes of the horizontal support displacement d. The second quantity (e, the
eccentricity of the compression force) can be different in every cross-section:
(− M / N )
e := ⋅100%.
(t / 6)

This quantity expresses whether the whole cross-section is under compression: for
e < −100% or for e > 100% the resultant of the normal stresses is outside the kernel of the
cross-section, which means that tension occurs at a part of the cross-section according to
the applied theoretical model and in the case of a no-tension material cracks should
appear.
Regarding the normalised horizontal reaction, first a sensitivity study was performed
to clarify whether the material properties of the vaults make any influence on δ.
The reactions were determined for circular arches having the same geometry with
t = 0.08 L, using different Poisson’s ratios (0–0.5) and Young’s moduli (ranging over
several orders of magnitude). The colour scale in Figure 7 indicates the corresponding
values of δ. Obviously, the normalized horizontal reaction was practically independent
of ν and E. Indeed, though the analysed arch model is statically indeterminate, since the
material parameters were always constant along the arch in every case and the geometry
was also the same, it was found that the equilibrium was maintained with the same
direction and position of the reactions. After this test, all further calculations were based
on properties of travertine limestone given in Section 2.4.

Figure 7 Normalised horizontal reaction for different material parameters (see online version for
colours)
 Horizontal reaction components of pointed vaults 411

Figure 8 shows how the normalised horizontal reaction depends on the shape of the arch.
The calculations were done for the arch thickness t = 0.08 L. Proceeding leftwards along
the horizontal axis the arch shape becomes more and more strongly pointed. According to
the Timoshenko beam model, the horizontal reaction significantly decreases with
increasing pointedness α: it drops by about 45% when reaching the arch shape based on
the equilateral triangle (α = 0.5), and by 55% in the case of α = 1.0, if compared to the
horizontal reaction of the circular arch (α = 0).

Figure 8 Normalised horizontal reaction as a function of the eccentricity of the vault shape, for
self-weight with fixed support (d = 0), for t = 0.08 L (see online version for colours)

Figure 9 is also specific to fixed supports. Arches having different thicknesses were
analysed under self-weight. Similar to Figure 8, it was found that pointed arches have
smaller normalised horizontal reactions, and the increasing thickness of the arch makes
the direction of the resultant reaction more vertical.

Figure 9 Normalised horizontal reaction: (a) as a function of the eccentricity of the arch, for
different thicknesses; (b) as a function of the thickness of the arch, for different
eccentricity parameters (see online version for colours)

The next diagram, Figure 10, focusses on the effect of outwards support displacements of
arches with the same thickness, t = 0.08 L. The analytical derivation predicts that the
increase of the support translation causes the horizontal reaction to decrease
monotonically. Figure 10 shows a linear dependence of the normalised horizontal force
on the magnitude of the displacement d, for every analysed arch shapes. For the circular
arch (α = 0) the force decreases to zero and then turns to outwards direction already at
around d/L = 0.017%, which corresponds to 1.7 mm for a 10 m span. Indeed, Figure 10
confirms that the elastic solution for the forces is very sensitive to the support
displacements.
412 G. Lengyel and K. Bagi

Figure 10 Normalised horizontal reaction for different support displacements (see online version
for colours)

Larger α causes a slower decrease of the horizontal reaction with the support
displacement, but the tendency is the same. Note that a negative δ means that since the
horizontal reaction is directed outwards, significant tensioned zones exist in the vault.
This is in contrast to real masonry behaviour, so the tension-resisting theoretical model
cannot suitably reflect the horizontal reactions for circular or only slightly pointed shapes
with even very small support displacements. The theory loses its validity for smaller
displacement values in case of circular or slightly pointed arches than in case of strongly
pointed shapes. The issue will be analysed now in more detail with the help of the
eccentricity of the compression forces.
The diagrams in Figure 11 illustrate how the e eccentricity of the compression force
varies along the length of the arch submitted to its self-weight, without support
displacements. (Note that going rightwards along the horizontal axis means proceeding
upwards along the axis of the arch.) Solid lines belong to thin arches (t = 0.08 L); dotted
lines correspond to thick arches (t = 0.20 L). According to the theoretical model, tension
occurs within those cross-sections where the value of e falls outside the interval (−100%,
+100%). The diagrams show that strongly pointed arches and thick arches are, in general,
more favourable: the eccentricities are smaller at most cross-sections. The neighbourhood
of the crown carries large bending moments in all shapes, particularly in the case of small
thickness. From among the three shapes (circular, pointed based on equilateral triangle,
pointed with α = 1.0) the arch based on an equilateral triangle is the most advantageous
from the point of view of the eccentricity of the normal forces.

Figure 11 Eccentricity of compression, for self-weight (see online version for colours)
 Horizontal reaction components of pointed vaults 413

Interesting to note that in the case of thick vaults the reactions act inside the kernel
(e is smaller than 100%) for both circular and pointed shapes.
The four curves in Figure 12 belong to four different magnitudes of support
displacement (red curve for fixed support). All diagrams were calculated for an arch with
t/L = 0.08 (thin arch) and α =0.0 (circular arch). The continuous beam model predicts that
the eccentricity becomes very large around the crown already for the slightest
displacement: d/L = 0.03% means 3 mm translation for a 10 m span, which is in the range
of construction errors.

Figure 12 Eccentricity of compression, for increasing support displacement (see online version for
colours)

At first sight these results seem to suggest that the theoretical model is reliable only for
fixed supports and thick structures. However, it will turn out from the discrete element
simulation results in Section 3 that though the support displacements indeed make the
continuous model unreliable, the arch thickness does not have an appreciable influence
on the validity of the theoretical prediction of the reactions.

3 Discrete element analysis of the horizontal reaction of pointed barrel

vaults

3.1 Geometrical and material parameters

3DEC models of barrel vaults of different α eccentricities were prepared to have the
same reference surface with a span L = 10 m. The shell thickness was 60, 80, or 100 cm
in the different models (t = 0.06, 0.08, and 0.10 L). The longitudinal length of the
analysed barrel slice was equal to the quarter of the thickness. The contacts between the
bricks were planar and perpendicular to the (y, z) plane (Figure 13) and to the generator
curve.
414 G. Lengyel and K. Bagi

Figure 13 3DEC model of barrel vaults: (a) the blocks of the analysed slice; (b) nodes on both
faces (front and back) are supported against translation in the x direction, to produce
plane strain state; (c) nodes only on the front face are supported in the x direction, to
allow plane stress state (see online version for colours)

The trapezoid-shaped discrete elements were made deformable by subdividing them into
tetrahedral uniform-strain finite elements. The finite element mesh density was
approximately the same everywhere, with the maximal edge size of the tetrahedra equal
to one quarter of the thickness of the masonry shell.
The blocks were linearly elastic, isotropic in the 3DEC model, with infinite resistance
to tension, compression, and shear. Material parameters corresponding to travertine
limestone were applied again, with those data given in Section 2.4. The joints between
the blocks were considered to be dry so that there was no resistance against tension in the
contacts, and the friction angle was set to 30 or 40 in the different tests. Based on the
analysis in Simon and Bagi (2015), the normal and tangential stiffness of the joints was
set to 103 GPa/m.

3.2 Boundaries and the applied loads

Table 2 summarises the boundary conditions applied in the 3DEC simulations. While the
analysed arch was in a plane stress state in the theoretical model, in the DEM simulations
plane stress and plane strain states were both simulated, in order to see the difference.
Considering a long barrel, plane stress state may be a good approximation towards the
ends of the structure if ends of the barrel are not sufficiently supported in the longitudinal
direction. Plane strain state may characterise the middle region of the barrel where the
neighbourhood of the analysed slice does not allow for longitudinal translations and,
consequently, longitudinal strains are equal to zero. The plane strain state was achieved
in the simulations by fixing all the nodes against longitudinal translations along both
faces of the slice (see the dark dots in Figure 13a). To get a plane stress state, only the
front face was supported (Figure 13b). Regarding the in-plane behaviour of the analysed
arch (Figure 5), the vaults were imagined to be supported by two longitudinal walls under
the barrel whose contacts with the blocks of the barrel resisted no tension but arbitrarily
high frictional force (the degree of friction was set to 90°). This was achieved by
allowing arbitrary upwards (but not downwards) vertical displacements and prescribed
outwards displacements in the bottom nodes of the lowest block of the arch (these nodes
formed the Γ0b boundary). Consequently, hinging and separation was possible between
the barrel and the imagined supporting wall, but sliding could not occur. The Γsym
 Horizontal reaction components of pointed vaults 415

boundary (neighbourhood of the crown) was defined in the way shown in Figure 13a: the
top block represented the half of the crown stone, and the nodes located on the axis of
symmetry were fixed against y-directional translation.

Table 2 Boundary conditions in the 3DEC simulation of barrel vaults

Boundary v w
Γ0b d *
Γsym. 0 –

3.3 Simulation results for fixed supports, self-weight only

Several arch geometries were analysed under self-weight with no displacements of the
supports. After equilibrating the vaults, the reaction components and the internal forces
could be recorded and then analysed. In general, a surprisingly good agreement was
found between the analytical solutions and the discrete element simulations, in spite of
the fact that significant regions of the discrete systems were partially cracked.
Figure 14 compares the analytical prediction given by the continuous Timoshenko
beam to the discrete element simulation results for the plane stress state and also for the
plane strain state. The diagram shows the normalised horizontal reaction as a function of
the geometrical eccentricity parameter α, for 40° friction angle and arch thickness
t = 0.08 L (a relatively thin arch). Though it was clearly shown in Section 2.6 (Figure 11)
that significant regions carry tension according to the continuous model, the comparison
with DEM brought an unexpected finding: the discrete element results practically
perfectly coincide with the theoretical predictions. Consequently, the presence of the
cracks in the discrete model does not cause a considerable difference in the distribution
of internal forces from those internal forces that were predicted by the continuous model.
This statement is checked in Figure 15 with the help of the eccentricity parameter of the
compression forces. The value of e is shown along the length of the arch, for α = 0
(circular shape) and α = 0.5 (pointed arch based on the equilateral triangle). Solid lines
represent the discrete element simulations and dotted lines denote the analytical
predictions. Indeed, the results on the position of the compression force show a very close
coincidence between the cracked discrete systems and the continuous beams, even at
those regions where the eccentricity of the force is large.

Figure 14 Comparison of the analytical and DEM results on the normalised horizontal reaction
(t = 0.08 L) (see online version for colours)
416 G. Lengyel and K. Bagi

Figure 15 Eccentricity of the compression along the arch: analytical and DEM results (t = 0.08 L;
plane stress state) (see online version for colours)

This surprisingly good agreement remains valid for other arch thicknesses and friction
coefficients. Very low frictional resistance (30°) and different thicknesses (t = 0.06, 0.08,
and 0.10 L) were tested, for plane stress and for plane strain as well. The results on the
horizontal reaction were practically identical (Figure 16).

Figure 16 Normalised horizontal reaction for different arch thicknesses and frictional resistance in
the contacts: (a) plane stress state; (b) plane strain state (see online version for colours)

These normalised horizontal reactions can be compared to those minimal and maximal
admissible values found by Romano and Ochsendorf (2010) with the help of limit state
analysis. For arches with t/L = 0.15 the horizontal thrust was found by them to vary
between the minimal value 0.24 … 0.14 times the half weight, depending on the shape of
the arch (the first value corresponds to the semicircular arch, α = 0, according to the
notations of the present paper, while the second value corresponds to the arch based on
the equilateral triangle, α = 0.5). The maximal value in Romano and Ochsendorf (2010)
was 0.58 … 0.3 times the half weight, depending again on the eccentricity of the arch
shape. Figure 16 shows that the results of the discrete element simulations and the
analytical solution given in the present study vary from 0.4 to 0.18, so they are in the
interior of the admissible domain found by Romano and Ochsendorf (2010).
 Horizontal reaction components of pointed vaults 417

3.4 Simulation results for support displacements

The arches already equilibrated for self-weight were then submitted to slow support
displacements. A quasi-static support displacement procedure was produced by applying
prescribed velocity histories on the Γ0b boundary. The displacement was produced in
1 mm steps. In each step the initial velocity was set to 0.0054 m/s and then gradually
decreased until reaching the 1 mm displacement, after which the system was stopped and
balanced (until the equilibrium error dropped below 10−5). This loading process was the
result of a rather lengthy experimenting with different values of velocities and
displacement step lengths, and these parameters gave the shortest computation times
while keeping the structure always close to equilibrium throughout the whole process (the
ratio of largest unbalanced force magnitude vs. average applied force magnitude always
remained under 1%) The reactions and internal forces were determined in the equilibrium
state.
A characteristic result is shown in Figure 17: with increasing displacements, the
horizontal force was decreasing first and then at a definite displacement value (about
0.04% of the span, depending on the exact vault shape) it reached a steady magnitude
that remained constant for further small displacements. (Obviously, the initially statically
indeterminate arch became determinate at this stage of the displacement process because
of crack opening producing internal hinges.) The magnitude of this constant value of the
normalised horizontal reaction will be termed the residual horizontal reaction. (Note that
this aspect of the behaviour is already impossible to reflect with the help of the applied
continuous analytical model: as previously illustrated by Figure 12, the compression
eccentricities become huge for much smaller displacements.)

Figure 17 Normalised horizontal reaction for support displacements in the DEM model
(t = 0.08 L; friction angle 40°): (a) range of small displacements; (b) range of large
displacements (see online version for colours)

The residual reaction is shown in Figure 18, for plane strain and plane stress state,
respectively. The curves in different colours belong to different arch thicknesses and
friction coefficients. The thick black line represents the normalised horizontal reaction
predicted by the analytical solution for self-weight, for t = 0.08 L.
418 G. Lengyel and K. Bagi

Figure 18 The magnitude of normalised residual horizontal reaction for different arch geometries:
(a) plane strain state; (b) plane stress state. Black line represents the analytical
prediction for plane stress state, for fixed supports, t = 0.08 L (see online version for
colours)

The residual reaction is, in general, slightly smaller in the case of plane strain state (i.e.
the middle region of a long aisle vault) than in plane stress state. The friction coefficient
does not considerably affect the residual reaction, but the thickness of the arch is
important: arches with smaller thickness express a larger normalized horizontal reaction
on the supporting structural members below.
The residual thrust values can also be compared to those found by Romano and
Ochsendorf (2010) with the help of limit state analysis. As already mentioned above, they
found that for arches with t/L = 0.15 the horizontal thrust should be between the minimal
value 0.24 … 0.14 times the half weight, and the maximal value 0.58 … 0.3 times the
half weight, depending on the eccentricity parameter of the arch shape. (Remember that
the first value always belongs to the semicircular arch, α = 0, while the other value
corresponds to the arch based on the equilateral triangle, α = 0.5) That thickness,
t/L = 0.15, for which detailed comparable data were available, was bigger than even the
largest thickness for which discrete element simulations were done in the present study,
t/L = 0.10. The DEM simulations gave that for the thinnest arch, t = 0.06 L, the horizontal
thrust was 0.4 … 0.23, very close to the analytical Timoshenko prediction. With
increasing thickness, the normalised horizontal residual thrusts decreased and approached
towards the minimal normalised horizontal thrust predicted by the limit state analysis: for
t = 0.10 L DEM gave 0.3 … 0.17.
To summarise, the analytical solution gives an upper boundary to the values of δres:
note that the smallest analysed thickness, t = 0.06 L, is close to Heyman’s minimal wall
thickness for the circular arch. So it can be concluded that the continuum-based
analytical solution gives a safe (but not economic) estimation of the horizontal reactions,
not only for fixed supports but also for small support displacements. The simulations
could have been continued until failure in order to get an insight into the variation of the
horizontal force, but the huge computational times (i.e. months) did not allow for such an
analysis.

4 Conclusion

Vaults with pointed arches, in general, express smaller horizontal thrusts on their
supports than vaults with circular or only slightly pointed arches.
 Horizontal reaction components of pointed vaults 419

The horizontal reaction components of barrel vaults with fixed supports can reliably
be predicted with the help of the proposed linearly elastic continuous beam model, in
spite of the fact that the model does not take into account the possibility of cracking. This
statement holds only for the case of no displacements at the supports.
In the case of increasing outwards support displacements the barrel vaults reach a
state when they abruptly become statically determinate, and after this the reactions
remain practically constant for further small increase of the displacements (a residual
horizontal reaction magnitude can be assigned to the barrel vaults). Then the increasing
displacements modify the geometry of the vault so significantly that the horizontal
reactions increase again, until failure. In the range of small displacements the continuous
model gives a good approximation of the residual reaction for thin structures, and a safe
but not economic approximation for vaults with large thickness.

Acknowledgements

The authors would like to thank Dr. Fabian Dedecker of ITASCA Consulting Inc. for his
support as well as the ITASCA Education Partnership Program for providing a copy of
3DEC. Financial support of the OTKA 100770 project is gratefully acknowledged.

References
Aita, D., Barsott, R., Bennati, S. and Foce, F. (2004) ‘The statics of pointed masonry arches
between “limit” and “elastic” analysis’, in Roca, P., Molins, C. (Eds.): Arch Bridges IV,
CIMNE, Barcelona, pp.353–362.
Aita, D., Barsotti, R. and Bennati, S. (2009) ‘Load bearing capacity of circular, pointed and
elliptical masonry arches’, Procs. XIX Congresso AIMETA, 14–17 Sept 2009, Ancona, Italy.
Aita, D., Barsott, R. and Bennati, S. (2015) ‘Some explicit solutions for nonlinear elastic depressed
masonry arches loaded to collapse’, Procs. XXII Congresso AIMETA, 14–17 September 2015,
Genova, pp.354–359.
Alexakis, H. and Makris, N. (2013) ‘Minimum thickness of elliptical masonry arches’, Acta
Mechanica, Vol. 224, pp.2977–2991.
Arnold, V.I. (1973) Ordinary Differential Equations, MIT Press, Cambridge.
Audenaert, A., Peremans, H. and De Wilde, W.P. (2004) ‘Static determination of the internal forces
and displacements in arch bridges’, TMS (The Masonry Society) Journal, Vol. 22, No. 1,
pp.97–109.
Audenaert, A., Peremans, H. and Reniers, G. (2007) ‘An analytical model to determine the ultimate
load on masonry arch bridges’, Journal of Engineering Mathematics, Vol. 59, pp.323–336.
Borg, S. and Gennaro, J. (1959) Advanced Structural Analysis, Van Nostrand, New Jersey.
De Rosa, E. and Galizia, F. (2007) Evaluation of safety of pointed masonry arches through the
Static Theorem of Limit Analysis. Procs. Arch’07, 5th International Conference on Arch
Bridges, Madeira, Portugal, pp.659–668.
Coccia, S., Di Carlo, F. and Rinaldi, Z. (2015) ‘Collapse displacements for a mechanism of
spreading-induced supports in a masonry arch’, International Journal of Advanced Structural
Engeering Vol. 7, pp.307–320.
Cochetti, G., Colasante, G. and Rizzi, E. (2011) ‘On the analysis of minimum thickness in circular
masonry arches’, Applied Mechanics Reviews, Vol. 64, No. 5, #051002 (27 pages).
Durand-Claye, A. (1867) Note sur la vérification de la stabilité des voutes en maconnerie et sur
l’emploi des courbes de pression. Annales des Ponts et Chaussées, Vol. 13, pp.63–93.
420 G. Lengyel and K. Bagi

Foce, F. and Sinopoli, A. (2001) ‘Stability and strength of materials for static analysis of masonry
arches: Durand-Claye’s method’, in Abdunur, C. (Ed.): Arch 01, 3rd International Conference
on Arch Bridges, Presses de l’École Nationale des Ponts et Chaussées, 2001 Paris,
pp.437–443.
Foce, F. and Aita, D. (2003) ‘Betwen limit and elastic analysis. A critical re-examination of
Durand-Claye’s method’, in Huerta, S. (Ed.): Procs. First International Congress on
Construction History, Inst. Juan de Herrera, 20–24 January 2003 Madrid, Vol. II, pp.895–905.
Hejazi, M. and Jafari, F. (2010) ‘Structural effects of brick arrangement and span length on mid-
pointed arches’, Advanced Materials Research, Vol. 133–134, pp.411–416.
Henrych, J. (1981) The Dynamics of Arches and Frames, Elsevier, Amsterdam.
Heyman, J. (1966) ‘The stone skeleton’, International Journal of Solids and Structures, Vol. 2,
pp.249–279.
Heyman, J. (1969) ‘The safety of masonry arches’, International Journal of Mechanical Sciences,
Vol. 11, No. 4, pp.363–385.
Huerta, H. (2001) ‘Mechanics of masonry vaults: The equilibrium approach’, inLourenco, P.B.,
Roca, P. (Eds.): Historical Constructions, Guimaraes, pp.47–69.
Kooharian, A. (1952) ‘Limit analysis of voussoirs (segmental) and concrete arches’, Journal of the
American Concrete Institute, Vol. 24, No. 4, pp.317–328.
Lemos, J.V. (2016) ‘The basis for masonry analysis with UDEC and 3DEC’, in Sarhosis, V.,
Bagi, K., Lemos, J.V. and Milani, G. (Eds.): Computational Modeling of Masonry Structures
Using the Discrete Element Method, IGI Global, Hershey, USA.
Lengyel, G. and Bagi, K. (2015) ‘Numerical analysis of the role of ribs in masonry cross vaults’,
Computers and Structures, Vol. 158, pp.42–60.
Milankovitch, M. (1907) ‘Theorie der Druckkurven’, Zeitschrift für Mathematik und Physik,
Vol. 55, pp.1–27.
Oh, S.J, Lee, B.K. and Lee, I.W. (1999) ‘Natural frequencies of non-circular arches with rotatory
inertia and shear deformation’, Journal of Sound and Vibration, Vol. 219, No. 1, pp.23–33.
Pippard, A.J.S. (1948) ‘The approximate estimation of safe loads on masonry bridges’, Civil
Engineer in War Vol. 1, pp.365–372.
Rizzi, E., Rusconi, F. and Cocchetti, G. (2014) ‘Analytical and numerical DDA analysis on the
collapse mode of circular masonry arches’, Engineering Structures, Vol. 60, pp.241–257.
Romano, A. and Ochsendorf, J. (2010) ‘The mechanics of gothic masonry arches’, International
Journal of Architectural Heritage, Vol. 4, No. 1, pp.59–82.
Simon, J. and Bagi, K. (2016) ‘DEM analysis of the minimum thickness of oval masonry domes’,
International Journal of Architectural Heritage, Vol. 10, No. 4, pp.457–475.
Szakaly, F., Hortobagyi, Zs. and Bagi, K. (2016) ‘Discrete element analysis of the shear resistance
of planar walls with different bond patterns’, Open Construction and Building Technology
Journal, Vol. 10, pp.190–202.
Tufekci, E. and Arpaci, A. (1998) ‘Exact solution of in-plane vibrations of circular arches with
account taken of axial extension, transverse shear and rotary inertia effects’, Journal of Sound
and Vibration, Vol. 209, No. 5, pp.845–856.