A method for the analysis of masonry arches

E. Ricci &E. Sacco

Università degli Studi di Cassino e del Lazio Meridionale, Cassino, Italy
M.D. Piccioni
Politecnico di Bari, Bari, Italy

ABSTRACT: In this paper, a numerical method for assessing the safety of masonry arches is presented. The
catenary and the thrust curve of existing arch differential equations are revised. The variational formulation of
the differential equations of the catenary and of the thrust curve is presented with its finite element approxi-
mation. The optimization problem, solving the nonlinear constrained minimization which allows determining
the thrust curve in equilibrium with a given load for a prescribed geometry, corresponding to the minimum
value of the thrust force is discussed. Numerical examples are presented in order to assess the ability of the
presented numerical procedure in solving the problem of the determination of the thrust curve, comparing the
analytical, when available, with the numerical solutions. The values of the minimum and maximum allowable
thrust forces are also provided.

1 INTRODUCTION case, the design loading conditions are assigned and

the shape of the arch can be derived to safety carry
1.1 Masonry arches the prescribed loads. In the second case, both the
Masonry buildings represent most of the histori- arch geometry and the loading conditions are given
cal and monumental constructions in the World. The and the problem consists in evaluating the degree of
evaluation of the safety state of this cultural heritage safety of the existing structure.
represents an actual and important structural prob-
lem.
The masonry material is characterized by a high 1.2 Analysis methods for arched masonry structures
compressive strength respect to the tensile one. This The history of the development of the methodol-
peculiarity turns out in the difficulty of building ogies for the analysis of arched masonry structures is
wide span roofing systems, which has been over- long and full of interesting ideas that promoted new
come by ancients through the use of arched and researches (Benvenuto, 1991; Heyman, 1998; Kur-
vaulted structures. The curved profile and the con- rer, 2008; Como, 2016).
trast between elements have the aim to make arches In the second half of seventeenth century, Robert
and vaults to be subjected mainly to compression. Hooke (1675) approached the problem of the ma-
The arch represents one of the most fascinating thematical determination of the thrust curve to eva-
structure to be studied. The correct understanding of luate the stability of the arch. Then, David Gregory
his static behavior let researchers to approach the (1697) wrote about the catenary and gave a mathe-
study of the vaults mechanical behavior; in fact, matical formulation of this curve applied to the arch
vaults and can be modeled as a series of arches analysis, stating that arches are stable when a cate-
linked together along one direction, as in the case of nary is found to be placed within their thickness.
barrel vaults, or more, as for the cloister vaults. This Then, during the eighteenth century in France, a dif-
points out the importance of studying the static of ferent view of the problem was presented by La
masonry arches. Hire, Couplet and Coulomb, who studied the arch as
It is central to distinguish two different problems, a structural element composed of wedges whose
when approaching the modelling of arches: the de- equilibrium can be influenced by friction (e.g. Cou-
sign of new structures and the valuation of the safety lomb 1763). Studies on the safety of the arch were
of the old ones for a possible restoration. In the first also proposed adopting a kinematic approach, i.e.
describing the possible collapse mechanisms (Baldi, O’Dwyer (1999) used the so-called Force Net-
1621). It was shown that an arch collapses when a work Approach to find the thrust surfaces which fit
certain mechanism activates and leads to the devel- in the vault thickness, by dividing the structure in
opment of four hinges at least, since the collapse of parts whose mass is considered to be lumped in
an arch is mostly due to the relative rotation among nodes linked together. It is assumed that the nodes
blocks instead of sliding, which rarely happens can move only vertically on axes inside the thick-
(Smoljanovic et al. 2013). ness vault while the equilibrium of the structure has
One of the most important contributions to the to be assured. It seems to be a very interesting me-
study of arched masonry structures was given by J. thod and its performances were demonstrated for
Heymann (1982). He stated that if it is possible to two-dimensional structures like domes.
find a thrust curve for an arch in equilibrium with In his PhD dissertation, Block proposed the
the applied loads, including the deadweight and this Thrust Network Analysis (TNA) methodology. The
line is completely inside the thickness of the arch, TNA, based on projective geometry, duality theory
then the arch is stable. This result is based on three and optimization techniques, is able to find the poss-
hypotheses: ible funicular solution for a vault subjected to vertic-
- sliding between the blocks is not allowed; al load. It turned out to be quite intuitive and easy to
- the tensile strength of masonry is null; use.
Fraternali (2010) formulated a network approach
- the compressive strength of masonry tends to
to identify the geometry of the thrust surface and the
infinite. associated stress field by using a predictor-corrector
These criteria can be applied to arches whichever procedure. The method allows the user to estimate
shaped and, on the basis of the static theorem of the the regions of the structure exposed to possible dam-
limit analysis, frees from the necessity of finding the age.
real thrust curve among all the admissible ones. In Limiting the interest at the analysis of masonry
fact, it is sufficient to find just one thrust curve res- arches, it could be remarked that several and some-
pecting this criteria to determine the arch stability. times complex geometries can be found for existing
With the spread of the Finite Element (FE) me- arches. As consequence, the analytical solution of
thod, the limit analysis was left apart in favor of a the differential equation governing the profile of the
series of models based on the FE analysis, which let thrust curve could not exist and the use of a numeri-
to take into account the constitutive law of material cal procedure can become necessary for the evalua-
and its linear and nonlinear behavior in a very de- tion of the state of safety of the arch. In this frame-
tailed way. work, Sacco (2015) proposed a simple procedure
Nowadays, a large use of nonlinear FE analysis is able to derive the thrust curve by enforcing the equi-
made for the evaluation of the state of stress and for librium equation of a discrete number of points of
the identification of the damage evolution for maso- the arch.
nry structures. This approach, able to lead to detailed Indeed, numerical methods for arches based on
information on the safety state of the construction, static techniques could become artificial when ex-
requires the use of material models characterized of- tended to vaulted structures. Thus, in order to derive
ten by many parameters, which could be not so sim- a numerical method that can be extended from
ple to determine. Moreover, strain and damage loca- arches to vaults, a clear mathematical formulation,
lization problems arise when damage models with such as the Finite Element method, is preferred.
softening are used in the computations. Due to this In the following, because of the complexity in the
drawback of the FE approach, many researchers modelling of masonry vaults, it has been considered
faced the evaluation of the state of safety of vaulted more appropriate to start with the analysis of a one-
masonry structures restudying the old techniques dimensional problem. In fact, a new stress based me-
based on the fundamental assumption that masonry thodology is proposed in this paper concerning the
is not capable of withstanding tensile stresses. analysis of masonry arches, with the aim of extend-
Moreover, the importance of determining an equili- ing these approaches to the study of vaulted struc-
brated solution, neglecting the kinematic compatibil- tures in the next future.
ity, has been rediscovered. The catenary problem is presented and the case of
In the last years, several methods based on the de- the thrust curve for an existing arch is discussed. In
finition of an equilibrated thrust network have been both cases, the differential governing equations are
proposed especially for vaults. Most of them are given and the analytical solutions are explicitly re-
mainly based on the idea proposed in the work by ported, when available. Then, the variational formu-
Schek (1974) concerning the Force Density Method; lations of the two differential equations are given
it considers a pin-jointed network structure whose with its finite element approximations. In the first
position of the network nodes is defined enforcing case, i.e. for the catenary problem, the nonlinear FE
the equilibrium condition between internal and ex- approach is discussed and an iterative procedure for
ternal forces. determination of the solution is presented. In the
second case, the FE problem becomes linear and the The analytical solution of the catenary equation is
problem of the determination of the thrust curves well known:
corresponding to the minimum and maximum value
of the thrust force is addressed. H  f 
y cosh  x  c1   c2 , (8)
The obtained tool appears easy to handle since it f H 
is function of simple parameters such as the geome-
tric position of the pinned nodes. Examples of the where c1 and c2 are constants integration that are de-
application of this method are also provided. termined accounting for suitable conditions.

2 VARIATIONAL FORMULATION OF THE

THRUST CURVE
2.1 The catenary equation
Let an arch be considered. The reference system
is fixed so that the x axis is located at the arch im-
post and the y axis corresponds to the symmetry
axis. In the following the introduced Cartesian coor-
dinate system is adopted; polar coordinates are not
used, as the objective of this work is to solve arch
problems characterized by complex geometry and
not specifically circular arches.
The arch is supposed to be subjected to a distri-
buted vertical load, f.
With reference to Figure 1, the equilibrium of the
arch along the horizontal and vertical directions
leads to: Figure 1. Equilibrium of the axial forces along the thrust curve.

d ( N cos  )
0, (1) When an existing arch, characterized by a given
ds
geometrical shape, is considered, the thrust curve
d ( N sin  ) and the mid-curve of the arch, where the distributed
 f . (2)
ds load is applied, are different curves, as schematically
illustrated in Figure 2.
From the Equation (1) it can be noticed that the In this case, the differential equation of the thrust
component of normal axial force along the x axis is curve (7) becomes:
constant. Therefore, it can be set equal to H, which
assumes the meaning of the thrust of the arch: Hy' ' f 1   ~
y '  0 ,
2
(9)
N cos   H , (3)
where y is the thrust unknown function, while ~ y is
so that the axial force can be written as: the mid-curve of the arch. If the arch is circular, the
H analytical solution of the differential Equation (9)
N . (4) exists and is represented in the form:
cos 
f R x 
Substituting Equation (4) in the equilibrium Equa- y  x arcsin  R2  x2 
tion (2), it gives: H  R2  x2  , (10)
d ( N tan  )  c1 x  c2
 f. (5)
ds where c1 and c2 are constants integration that are de-
Taking into account that: termined accounting for suitable conditions and R is
the radius of the arch.
tan    y ' , ds  dx 1   y '  ,
2
(6) Note that in the case, characterized by circular shape
of the arch, the analytical solution of the Equation
Equation (2) becomes: (9) can be determined. In general, if the function y is
a polynomial or describes an elliptical arch, the ana-
Hy '' f 1   y '   0 , lytical solution cannot be found. Moreover, even if
2
(7)
the function y can be integrated, the analytical solu-
representing the differential equation of the catenary. tion is not always the easiest way to solve the prob-
lem. The following example can clarify this concept.
 where L is the total span length of the arch. The
formulation (11) represents the basis for the devel-
opment of the FE approach for the determination of
approximated solution of the catenary curve. The
domain is discretized by dividing the span of the
arch in ne two-node elements with length l(e). Each
node has one degree of freedom, representing the
coordinate of the thrust curve.
In the typical e-th element, defined by the two
nodes located at x1(e) and x2(e), the y coordinate de-
scribing the catenary is approximated as:
y  y1( e )1  y2( e ) 2 , (12)
(e) (e)
where y1 and y2 are the values of the unknown
Figure 2. Equilibrium of the axial forces when the midline of the function y at x1(e)and x2(e); the approximation func-
arch differs from the thrust curve. tions are set as:
x2  x x  x1
1  , 2  , (13)
l (e) l (e)
Let a polycentric arch be considered. The mid-curve
is described by a function y which is circular and it with l(e) = x2(e)-x1(e) length of the element.
is characterized by different centers and different ra- Substituting the approximation (12) into the Equa-
dii. Graphically a curve like this can be easily repro- tion (11), the following system of equation is ob-
duced, once the centers and the radii are known, but tained in residual form for each finite element:
analytically it is piecewise defined. It means that the R (e)  H K (e) y (e)  F(e)  P(e) , (14)
solution is different for each part of the arch in
which the centers and the radii are different. Moreo- with
ver, when a survey has to be done, the distance of
the points of the arch from a reference system is 1  1 1
K (e)  ,
measured in terms of Cartesian coordinates xy, l ( e )  1 1
which is clearly the easiest way. For these reasons it x2 2  
seems more convenient to use a method that can be F ( e )   f 1   y1( e )1 ' y2( e )2 '  1  dx ,
generally and easily applied and whose validity does x1 2 
not depend on the shape of the considered arch. This (15)
 P (e) 
is the case of the numerical methods which let the P ( e )   1( e )  ,
problem of the determination of the thrust curve be  P2 
solved for whichever shape and loading conditions.  y1( e ) 
In the present article, only the case of arches sub- y   (e)  ,
(e)

jected to vertical distributed loads is considered.  y2 

where P ( e ) is the vector accounting for the interaction
2.2 Variational formulation among adjacent elements and R(e) is the residual
term that has to be equal to {0 0}T in solution. Eq-
Though the analytical solution of the catenary
uation (14) represents a nonlinear algebraic system
problem is available, the determination of the thrust
of equations that can be solved developing an itera-
curve for more complex arch geometries or loading
tive procedure based, for instance, on the Newton
conditions could become difficult and numerical me-
method. In particular, the series development of Eq-
thods can be developed.
uation (14) gives:
The variational formulation of the problem of the
catenary is derived by multiplying all the terms of R ( e )
,k

Equation (7) by a variation δy of the unknown func- R (e)

R ( e ), k
 (e)
y
y (e)
 y ( e ), k 
(16)
tion y defining the thrust curve, and integrating
along the span of the arch:  R ( e ), k  K t ( e )  y ( e )  y ( e ), k   0 ,

 
L
 Hy '' f 1   y '   y dx
2
0 where the subscript k denotes the iteration number.
0
L L
(11) The 2  2 tangent matrix K t( e ) is obtained as:
  H  y '  y ' dx   f 1   y '   y dx  H  y '  y 0
2 L
,k
F ( e )
0 0
K (e)
t HK (e)
 (e) . (17)
y
Assembling Equations (16) for all the elements of  H min  max H ej
the discretization, the full linear system of equations H min  H  H max with  . (25)
 H max  min H j
i
at the iteration k+1 results:
R  R ,k  K t  y  y ,k   0 , (18) In this way, a possible range of value for the thrust
H is found.
with evident meaning of the symbols. The problem of the determination of the thrust
The finite element formulation of the differential curve corresponding to the minimum or maximum
Equation (9) leads to a linear algebraic problem. For thrust force can be found by minimizing or max-
the typical e-th element, Equation (9) becomes: imizing the value of H under the conditions (22) and
(25), taking into account positions (23). In particu-
H K ( e ) y ( e )  F ( e )  P ( e )  0 , (19) lar, in the case of minimum thrust force, the follow-
ing optimization problem is recovered:
where
x2
  
min H 2 / χ  K 1 F ,
F ( e )   f 1   y '   1  dx .
2 H
(20)
x1 2  j
H 0, (26)
The equation governing the problem of the whole y ej
structural system obtained assembling the element j
Equations (19), takes the form:   H  0 with j  1, 2,..., nd  .
y ij
H K y  F , (21)
with evident meaning of the symbols.
The linear system of equations (21) admits a 3 APPLICATIONS
unique solution once the value of the thrust H is as-
signed and suitable boundary conditions are pre- In order to assess the proposed method for the de-
scribed. As consequence, for any assigned value of termination of the catenary and of the thrust curve
the thrust it is possible to compute a solution, i.e. a for an existing arch, some numerical applications
shape for the thrust curve. have been developed.
Setting χ  H y , equation (21) gives: Let a round arch be considered. The internal and
external radii are Rin and Rex respectively, so the
χ  K 1 F . (22) midline radius is Rm = (Rin+Rex)/2. The arch is sub-
jected to a vertical distributed load f. All the geome-
Note that the solution of this problem is unique in trical and loading data are listed in Table 1.
terms of the new introduced variables j. The following two boundary conditions have
Let an arch characterized by a given thickness be been set to solve Equation (21):
considered. In this case, the profiles of the extrados
and of the intrados are defined by the functions ye(x) x   Rm  y0
(27)
and yi(x), respectively. The values of these functions x  Rm  y0
in correspondence of the finite element nodes allow
to introduce the two vectors ye and yi, respectively. so the first and the last node of the discretization are
An interesting problem concerns the determina- fixed; in such a way, the extremes of the thrust curve
tion of the thrust curve corresponding the minimum coincide with the midpoint at the impost of the arch.
possible value of the thrust force. To this end, each
component of the vector χ is divided by the corres-
ponding component of ye and yi, obtaining two vec- Table 1. Geometrical and loading data.
tors grouping the values of possible thrust forces He Radius Span Load
and Hi, respectively: [mm] [mm] N/mm
Extrados 1700 3400
j j
H  e
j , H  i
j with j  1,.., nd , (23) Intrados 1300 2600
y ej y ij Midline 1500 3000
f 2
with nd the total number of node of the FE discreti-
zation. In particular, it can be noticed that:
y ej  y ij  H ej  H ij with j  1,.., nd , (24) In Figure 3, the thrust curve obtained solving the
problem governed by Equation (21) with the boun-
so that the thrust H can assume a value included be- dary conditions (27) is plotted for a value of the
tween the maximum of Hje and the minimum of Hji, thrust force H = 1856 N. In the same figure, the ca-
i.e. it is: tenary curve (dashed line) is also reported for the
same value of the thrust force. The prescribed value
of the thrust force ensures that both the thrust and
the catenary curves are completely inside the arch.
The differences between the two curves clearly ap-
pears.

Figure 5. Thrust curve corresponding to the maximum thrust

force.

The method is also applied to the analysis of a Tudor

arch characterized by the geometrical and load con-
Figure 3. Line of thrust for a round arch with set load (continuous
line). First attempt. Comparison with the catenary curve (dashed line). ditions listed in Table 2:

Table 2. Tudor arch. Geometrical and loading data.

The analytical solutions of the catenary, Equation Rise Span Load
(8), and of the thrust curve, Equation (10), are also [mm] [mm] N/mm
reported in the same Figure 3, resulting completely Extrados 1208 3300
superimposed to the results obtained by the two nu- Intrados 988 2700
merical procedures. Midline 1098 3000
Now, the interest is devoted to find the solution f 1
of the thrust curve corresponding to the minimum
and the maximum value of the thrust force.
To this end, the optimization problem (26) is ap- The Tudor arch is a polycentric arch characte-
proached and solved. The thrust curve corresponding rized by a ratio rise/span quite low. It means that the
to the minimum value of thrust force is obtained for range of value the parameter H can assume to de-
H = 1518.76 N and it is plotted in Figure 4. termine a line of thrust fitting the width of the arch is
Analogously, following a similar procedure, the very small. As for the circular arch, the first and the
thrust curve corresponding to the maximum value of last node are bounded so that they coincide with the
thrust force can be determined. The solution is plot- extremes of the mid-curve. The minimum and the
ted in Figure 5, corresponding to H = 1851.83 N. maximum value that H can assume are equal to 1030
N and 1100 N, respectively. In Figure 6 the analyzed
arch and the thrust line obtained for the minimum
thrust is plotted, while in Figure 7 the thrust line cor-
responding to the maximum value of H is plotted.
It can be noticed that the interval of values that H
can assume is small as a consequence of the particu-
lar shape of the arch, which has a very small rise if
compared to the span.

Figure 4. Thrust curve corresponding to the minimum thrust force.

 cation, or when it is subjected to tricky loading con-
ditions.
Moreover, the proposed procedure represents a
starting point for the development of a numerical
formulation of two-dimensional vaulted structures.

REFERENCES
Benvenuto, E. 1991. An Introduction to the History of Struc-
tural Mechanics. Part I, Statics and Resistance of Solids;
Part II, Vaulted Structures and Elastic Systems. Springer -
Verlag, New York.
Figure 6. Thrust curve for a Tudor arch corresponding to the min- Block, P. 2009. Thrust Network Analysis: exploring three-
imum thrust force. dimensional equilibrium. PhD Thesis Massachusetts Insti-
tute of Technology.
Como, M. 2016.Statics of Historic Masonry Constructions.
Springer International Publishing Switzerland.
Coulomb, C.A. 1773. Essai sur une application des règles des
Maximis et Minimis àquelques Problèmes de statique rela-
tivs àl’Arquitecture. Mémoires de mathématique et de phy-
sique présentés àl’Academie Royal des Sciences per divers
savans7: 343–382.
Fraternali, F. 2010. A thrust network approach to the equili-
brium problem of unreinforced masonry vaults via polyhe-
dral stress functions. Mechanics research Communica-
tions37(2): 198-204.
Gregory, D. 1697. Catenaria. Philosophical Transactions of the
Royal Society 19(231): 637–652.
Heyman, J. 1982. The Masonry Arch. Chichester: Ellis Hor-
wood Ltd.
Figure 7. Thrust curve for a Tudor arch corresponding to the max-
imum thrust force.
Heyman, J. 1998. Structural Analysis: A Historical Approach.
Cambridge University Press, Cambridge.
Kurrer, K.E. 2008. The History of the Theory of Structures:
From Arch Analysis to Computational Mechanics. Ernst
&Sohn, A Wiley Company, Berlin.
O'Dwyer, D. 1999. Funicular analysis of masonry vaults. Com-
4 CONCLUSIONS puters & Structures 73(1-5): 187-197.
Sacco, E. 2015. Some Aspects on the Statics of Masonry
The variational formulation and its finite element Arches. In Aita D., O. Pedemonte & K. Williams (eds),
Masonry Structures: Between Mechanics and Architecture:
discretization of the differential equation governing 265-290. Springer International Publishing Switzerland.
the classical problem of the catenary and of the Schek, H.J. 1974. The force density method for form finding
thrust curve for an existing arch are presented. A and computation of general networks. Computer Methods
numerical algorithm for solving the nonlinear alge- in Applied Mechanics and Engineering 3(1): 115-134.
braic problem of the catenary is developed. Smoljanovic, H.,N. Živaljić, Z. Nikolić2013. Overview of the
methods for the modelling of historical masonry structures.
It is found that, once suitable boundary condition Gradevinar 65(7): 603-618.
are prescribed, the thrust curve is uniquely defined
when the value of the thrust force is assigned. It is
demonstrated that for round arches, the numerical
solutions are in perfect agreement with the analytical
ones, both for the catenary and for the thrust curve.
An optimization problem is formulated which al-
lows determining the thrust curve corresponding to
the minimum value of the thrust force; analogously,
the thrust curve for the maximum value of H is also
determined.
The proposed numerical procedure appears sim-
ple and reliable for the determination of the admissi-
ble and equilibrated stress states ensuring the stabili-
ty of the arch. The illustrated applications represent
the assessment of the presented numerical proce-
dure. Of course, its use is justified when the arch is
characterized by complex geometries, as in the case
of the Tudor arch presented as last numerical appli-