1. Introduction
The assessment of the equilibrium of masonry structure is a topic of high relevance, given the diffuse presence of masonry buildings in Italy and Europe, often with high historical value.
Over the centuries, different methodologies have been developed for the assessment of the equilibrium of masonry constructions and to quantify their degree of safety, starting from practical rules [
1], until the definition of graphical methods applied to two dimensional structures (such as portals, bridges, and flying buttresses [
2,
3]) and to axis-symmetric vaults under monoaxial stress regimes, such as masonry domes [
4]. Graphical methods have been extended to consider also biaxial stress states, again under the hypothesis of axial symmetry of the structure [
5].
For masonry material, the starting hypotheses have been clearly formulated by Heyman in [
6,
7], and assume that the material is infinitely resistant in compression and that sliding between rigid blocks does not occur inside the material. Under these assumptions, the theorems of Limit Analysis hold, and an equilibrium solution is defined as admissible if it is balanced with the given loads and is everywhere of pure compression. In particular, based on the Safe Theorem, the impossibility of collapse of a masonry structure, under the given loads, is ensured if at least one admissible equilibrium solution can be found.
A wide number of methods have been developed for the assessment of masonry vaults within the Heyman’s theory, both considering the kinematic and the Static Theorem. Concerning the Kinematic theorem of Limit Analysis, the proposed methods are based on energy criteria minimization. Indeed, the minimum of the total potential energy can be used to look at both stable and non-stable mechanisms [
8,
9]. Moreover, it can be also used to couple internal stress states with the corresponding mechanisms as in [
10]. Other methods, such as the Thrust Network Analysis (TNA), are based on the application of the Static Theorem. TNA, introduced in [
11] and developed in [
12,
13,
14,
15,
16], consists in the discretization of a masonry vault through a finite number of compressed bars, in equilibrium with the load of the vault, and located in-between the intrados and the extrados [
17].
Among different continuum approach methods, here we focus on the Membrane Equilibrium Analysis (MEA), first proposed in [
18], and further developed in [
19]. For the assessment of the equilibrium of a vault, MEA requires to find at least one membrane surface comprised in-between the geometrical boundary of the vault, balancing the external loads through a negative semi-definite membrane stress field.
In MEA, the equilibrium is reduced to a scalar second order partial differential equation with non-constant coefficients. If the shape of the membrane is assigned, the coefficients are the components of the curvature of the membrane, and the unknown of the problem is the stress potential of the membrane stress projected onto a reference plane (Pucher stress). On the contrary, if the stress potential is assigned, the coefficients of the equation are the components of the Hessian of the stress potential and the unknown is the scalar function defining the elevation of the surface with respect to the reference plane. In both cases, the equilibrium equation may change its behaviour (hyperbolic-parabolic-elliptic) depending on the assigned shape or on the assigned stress potential. This has a strong implication on the numerical solution of the equilibrium equation, that can be solved as a boundary value problem when the equation is elliptic, or as an evolutionary problem when the equilibrium equation is hyperbolic at least in a part of the domain.
MEA has been adopted mainly for the assessment of the equilibrium of masonry constructions [
20,
21,
22,
23,
24], but it is a general tool for computing the stress field on curved membrane surfaces. For example, it has been used for the evaluation of the stress field of soft tissues, such as the human cornea [
25].
An approach similar to MEA, based on the search of the Airy stress function, has been proposed in [
26], where NURBS are used for the discretization, exploiting their high-order continuity.
In [
27] a numerical iterative procedure for the assessment of the equilibrium of masonry vaults, based on LSM (a FEM-like approximation of the governing equation, see [
28]), has been introduced. Starting from a tentative membrane shape, included in the volume of the vault, the method looks for an equilibrated concave stress potential, ensuring the negative semi-definiteness of the corresponding stress. This procedure, however, is not easily applied to non-concave shapes, such as cross vaults given the non-ellipticity of the corresponding equilibrium equation. Recently, in [
29] a numerical procedure based on standard FEM, for the design of purely compressed shells, has been put forward. With this method, the starting point is a class of concave stress potentials, and the equilibrium equation is purely elliptic.
In the present paper, focusing on the assessment of the equilibrium of existing masonry vaults, we start from a family of concave stress potentials in order to ensure the admissibility of the stress field. This implies the ellipticity of the equilibrium equation: Therefore, we look for an optimal membrane shape comprised between the geometrical boundary of the vault, as required by the Safe Theorem of Limit Analysis. We approach this problem through an optimization algorithm, minimizing an error function first with respect to the parameters of the stress potential, and then with respect to the boundary conditions of the the membrane shape.
Optimization is a recurring tool in different engineering applications. In the present work, we find an optimal ideal membrane shape, which best satisfies the requirements of the Safe Theorem of Limit Analysis, which result in some mathematical constraints. We point out that this is not the case of a shape optimization, since the shape of the analyzed vaults is fixed, and we just optimize an ideal membrane. In other engineering application, however, the shape of structures is optimized under the constraints of material saving, or compliance minimization. This is the case of the Structural Topology Optimization [
30], used for optimally distribute the material to minimize the compliance o structures under given loads [
31,
32,
33]. Optimization is also used in geotechnical applications [
34], optimization of geodesic domes [
35], of soil-steel structures [
36], and lightweight structure [
37]
3. Outline of the Method
The Membrane Equilibrium Analysis (MEA) applied to masonry vaults, requires to find a membrane surface located in between the intrados and the extrados of the vault and, according to Heyman’s assumptions, sustaining the loads through compressive stresses. The membrane is a determined structure, and this last constraint is generally not satisfied on a given membrane shape for given loads. In such a case, the membrane shape has to be changed iteratively from an initial guess shape, until one that satisfies the unilateral restrictions on the stress is found.
3.1. Unknowns and Data in the Governing Equilibrium Equation
On adopting a numerical discretization method based on a variational approximation of the governing Equation (
6), such as the Finite Element Method, or the Lumped Stress Method defined below, the differential equation must be elliptic. Therefore, if the unknown is the stress and the given form of the membrane is not everywhere convex or concave, these methods of approximation fail and cannot be applied. A possible way out, already proposed in [
29], consists in switching the roles of unknowns and data, by assuming a parametric family of concave stress potentials as given, and considering the shape as the unknown. Then, the equation is elliptic for any choice of the parameters defining a specific member of the given family of stress potentials, and the corresponding shape can be obtained by prescribing convenient boundary data for the shape. Besides, for any such concave stress potential, there exists a small repertoire of purely compressed equilibrium shapes, that one can also explore, depending on the possibly different boundary data, that can be imposed on the shape. In the present paper, the equilibrium solution is tackled by constructing an LSM approximation of the fundamental second order differential equation of equilibrium (
6). This computer routine, which returns a shape associated to any given stress function, can be manipulated by changing, following a minimization algorithm, the parameters controlling the stress potential and the boundary data for the shape.
3.2. Numerical Discretization and Algorithm
The transverse equilibrium Equation (
6) is solved numerically by discretizing the domain into a finite number of triangular elements, and adopting a Lumped Stress Method [
28] approximation of the Hessian of both the stress potential
F and of the membrane shape
f. With the LSM approximation, the two functions are described by simplicial surfaces, defined by the nodal values
and
, with
,
N being the total number of nodes of the given triangulation.
The vectors collecting the nodal values of the shape function f and of the stress potential F are denoted and . is divided into two vectors and , of length and , distinguishing between internal and boundary nodes, respectively.
The search of a statically admissible stress regime for the vault is carried out through an algorithm divided in several steps. Given the geometry of the vault, defined by its intrados and extrados, we start on assigning an appropriate functional form for the stress potential, depending on a limited number
n of parameters, such that
A first step is the parametric solution of the discretized version of Equation (
6), with suitable boundary conditions, in terms of the membrane shape defined by
f; therefore,
.
The second step consists in the optimization of the parameters , performed by minimizing an error function I defined as the integral of the quadratic approximation of the characteristic function of the vault domain.
Given a point
, the characteristic function indicates whether the corresponding values of the rise of the membrane is comprised between the volume of the vault, or not. In particular
Its quadratic approximation reads
For a generic point
of the domain, given
and
, the approximated characteristic function is represented in
Figure 2.
We denote the subregion of the entire planform obtained by removing from the domain a boundary strip of size c, defined as a fraction of the diameter L of (in the applications we set ). The rationale behind this procedure is linked to the nature of the governing equation: due to ellipticity, the solution is only weakly dependent on the boundary data, that is, altering the boundary values has a detectable effect only in a narrow band close to the boundary. So, we can split the optimization into an interior and an exterior part, postponing the optimization of the boundary data.
The error function
I is then defined as follows:
and
being the parts of the domain
, for which the membrane surface is below the intrados, or above the extrados, respectively.
In the third step, the values of
, found in the previous step, are fixed and the boundary values of the discretized membrane shape are considered as variables. Then, the discretized version of Equation (
6) is solved parametrically with respect to these boundary values. In the vector notation introduced above, the transverse equilibrium equation reads
where the summation convention on repeated indices is assumed on the lower-case indices
, whilst the upper-case index
is not summed, and
,
, and
are coefficient matrices computed via the LSM scheme, representing the discrete counterparts of the second derivatives
,
, and
.
Equation (
15) can be rewritten in a compact form by factoring out the unknown term
, as
where
Expression (
16) can be also rewritten distinguishing between boundary conditions and domain equations, in the form
where
and
are submatrices of
corresponding to rows and columns associated to internal and boundary nodes. Finally,
is the vector of the boundary parameters. From Equation (
18),
can be evaluated as
On introducing (
19) into a discrete version
of the objective function
I, we get
We then optimize with respect to .
Notice that it is possible to compute explicitly the gradient of
through the chain rule, in the form
that is
The membrane shape obtained at the end of this minimization process can be compared to the actual geometry of the vault. If the final value of the objective function
is zero, the membrane is completely comprised between the intrados and the extrados of the considered vault, and therefore the solution is statically admissible and proves the stability of the vault under the given loads, according to the Safe Theorem of Limit Analysis. In such a case, the degree of safeness of the structure can be assessed by reducing the slenderness of the vault and computing a geometrical safety factor, as proposed by Heyman [
6].
4. Application
4.1. Cloister Vault
Consider the case of a cloister vault under its self-weight, whose geometry is described by the volume enclosed in-between two surfaces, the intrados and the extrados, defined canonically inside the domain
as follows:
where
,
,
, and
are the subregions of
depicted in
Figure 3.
A perspective representation of the midsurface of this vault is represented in
Figure 4.
The values adopted in Equations (
23) and (
24) are
m,
m, and
m. The values have been taken by the geometry of the vault analyzed in [
21], and located in the interior of Palazzo Caracciolo di Avellino in Naples. From [
21] we also assume
kN/m
as distributed load per unit projected area.
By following the methodology described in
Section 3, the transverse equilibrium Equation (
6) can be solved for the geometry
, by assuming a special parametric form for the stress potential.
For the case at hand, inspired by the analysis described in [
21], we identify a central zone
of the domain
, in which the projected stress is a plane pressure.
is delimited by a curve
defined as
where
is a parameter,
is set to
and
P to 40. In
Figure 5 we show the contour
obtained for two representative values of
, namely
and
.
Inside
, the Pucher stress is assumed as isotropic, that is
The corresponding stress potential is therefore easily computed as
where
is a parameter defining the intensity of the stress potential (i.e., of the stress) that will be subsequently optimized. In the external zone
the stress is assumed to be quasi-uniaxial in the radial direction. This means that, on defining
as the direction of the radii emanating from the geometrical center of the domain and directed toward the external boundary, the stress in
can be written as
where
is related to
and denotes the stress intensity along the radius, and
is a small stress which has the function to give a small extra-curvature to the potential in the direction
orthogonal to
, in order to preserve the ellipticity of Equation (
6) when solved for the shape
.
The stress potential associated to the stress distribution (
28) in the external zone is obtained by imposing a linear growth in the radial direction, enforcing the continuity of the potential
F and of its slope across the curve
:
where
s is a parameter indicating the distance from
along the radius.
In
Figure 6, we show the the stress potential for two representative values of
. Adopting these two shapes as the stress potential, and taking as boundary condition for the membrane geometry
the value of the rise of the middle surface between the intrados and the extrados of the vault, Equation (
6) can be solved using LSM, as described in
Section 3.2. By constructing a structured triangulation based on a regular node distribution over the planform
, the membrane shapes corresponding to the stress potentials of
Figure 6, are represented in
Figure 7.
Figure 7 gives an idea on how the choice of the parameter
(and thus, the extension of the central zone characterized by a biaxial stress state) influences the final shape of the membrane geometry. In the second case (
), the membrane geometry resembles the shape of a cloister vault more than in the case of
. It is also evident that the height of the two shapes (in the direction of the
-axis) are not comparable with the height of the cloister vault defined from Equations (
23) and (
24), which is
m.
The minimization of the error function (
20) with respect to the parameters
and
is performed as follows: first,
is given a finite number of values in the range
, i.e.,
. Then, for each value of
, we search the
that minimizes
. In
Figure 8 we show the considered values of the
on the horizontal axis, and on the vertical axes we show the corresponding optimal value of
.
The lowest value of
is reached when
, which correspond to
. Therefore the parameters of the stress potential are fixed to
and
. The best-fitting membrane geometry is then represented in
Figure 9.
In
Figure 10 we also show two sections of the membrane and of the intrados and extrados, the first one taken at
, and the second one taken along the diagonal of the domain, that is for
.
We notice from
Figure 10b that some points of the membrane fall outside the volume of the vault; this is confirmed also by
Figure 11, where we show the color map of the function
defined as
for which positive values indicate points in which the membrane shape is above the extrados of the vault, negative values indicate when the membrane is below the intrados, and vanishing values indicate nodes that are comprised within the volume of the vault.
We then perform the last step of the optimization, with respect to the boundary values of the membrane shape, obtaining a new membrane geometry whose sections for
and
are shown in
Figure 12. In this case, we notice that the points that were falling outside the volume of the vault before this optimization step, are now comprised between the intrados and the extrados.
The whole algorithm has been repeated assuming the same loads, but reduced thickness, in order to assess the geometrical safety factor. The higher thickness reduction for which the proposed algorithm still produces admissible solutions, in the sense of Limit Analysis, is of the original thickness, corresponding to a geometrical safety factor of .
4.2. Cross Vault of the Anterior Portico of the Church of St. Pietro in Vineis in Anagni
In this section we consider the case of the anterior portico of the church of San Pietro in Vineis in Anagni. The portico consists of three similar cross vaults: the geometrical description of each vault, referring for notation to
Figure 13 and
Figure 14, is given by:
The values given to the geometrical parameters introduced in Equations (
31) and (
32), are reported in
Table 1.
For this kind of vault a concave membrane surface comprised between the intrados and the extrados does not exist; as a consequence, using Equation (
6) with a known shape and unknown stress potential would lead to a non-elliptic problem, that cannot be solved using a FEM-like procedure. The inverse problem, consisting in the search of the shape corresponding to a given concave stress potential, is elliptic instead, and the FEM-like procedure proposed in [
21] can be applied.
It is a well-known result (see [
27]) that the stress potential corresponding to a standard cross vault has the shape of a cloister vault. The case at hand, however, is not a standard cross vault since there is a curvature not only orthogonally to each median of the planform, but also in the median direction. We expect, therefore, a biaxial stress field, with components both along the direction of the medians, and in its orthogonal direction.
We propose therefore the following stress potential
For sake of example, in
Figure 15 we present three different stress potentials obtained for
,
, and
and
.
We assume the same definition of the error introduced in Equation (
20) and evaluate it for different values of
, and then minimizing
with respect to
. The values of
considered in this case are
The corresponding plot of the error
versus
is reported in
Figure 16 in a logarithmic scale along the error axis.
The minimum occurs at
, and corresponds to an error of
. The corresponding stress potential is shown in
Figure 17.
The corresponding shape of the membrane, obtained through the numerical solution of Equation (
6), and imposing as boundary values for
the values of the midsurface (between the intrados and the extrados), is shown in
Figure 18.
To evaluate the position of the membrane surface after this first optimization phase, we show in
Figure 19 the colormap of
, defined like in Equation (
30), applied to this case study. The presence in
Figure 19 of points for which the error function is non zero indicates that the membrane found after the first optimization step is not fully comprised between the intrados and the extrados. Therefore we need a further optimization step to achieve this goal.
The second step of the minimization process is then performed with respect to the boundary values of the membrane shape. The value of the error function , obtained as a result of the second optimization phase, is zero everywhere, meaning that the membrane fits in-between the bounding surfaces of the vault.
For sake of example, we show in
Figure 20 a graph of a section of the intrados, the extrados and the membrane shape along the diagonal line
, that, after the first minimization phase, exhibits points lying outside the vault geometry. As expected, there are no points of the membrane section falling outside the geometry of the vault.
The influence of changing the boundary values
is made evident by the color map of the difference
defined as
and shown in
Figure 21, where
f is the membrane geometry at the end of the first minimization step, and
is the geometry after the final minimization with respect to the boundary conditions.
From
Figure 21 we see that the main effects of the modification of the boundary conditions is that the rise change of the membrane, represented by Equation (
35), after the second optimization phase, are visible mainly in the boundary part of the domain, as expected from the discussion about elliptic partial differential equations. Moreover, we see that
is negative in the lowest-right and upper-left corners, meaning that the membrane rise is diminished after this optimization, fixing the error shown in
Figure 19.
In order to quantify the Geometrical Safety Factor for this vault, the algorithm has been repeated progressively reducing the thickness. The maximum possible reduction, given the assumed shape for the stress potential, is , corresponding to a GSF .
4.3. Discussion
Concerning the algorithm adopted in both examples, we point out that it is an automatic procedure that is easily adaptable to the specific case. The only different step between the proposed algorithms, beyond the geometries of the vault, is the parametric description of the stress potential, which is specific to each case study, and consists, from an operative point of view, in replacing a computer routine to another one.
Concerning the Geometrical Safety Factor, computed at the end of each application, we remark thaat it is an intrinsic characteristic of each structure, and does not depend, in general, on the adopted membrane shape or its associated equilibrated stress. With computational models, such as MEA, it is only possible to try to find the best value, that, in turn, is associated to specific choices on the stress potential (i.e., its shape, or its parametric description), or on the discretization, or on the specific algorithm.
In our examples, we made two precise choices for the shape of the stress potential, based on experience and on literature examples. A richer description of the stress potential could have improved the obtained value of the Geometrical Safety Factor, at the price of a more complicate and slow optimization algorithm, especially for the first optimization pahse, which is characterized, for each evaluation of the objective function, by a matrix inversion, whose computational cost, in turn, depends in the adopted discretization. Therefore, the proposed values for the GSF have to be considered as conservative values, and could be improved by using more refined optimization algorithms, not yet developed in the framework of the Membrane Equilibrium Analysis.
5. Conclusions
In this work we present a method for the assessment of the equilibrium of masonry vaults based on the search for a suitable membrane shape in equilibrium with a given compressive stress field. The membrane shape is required to be entirely comprised between the intrados and the extrados of a vault, according to the Safe Theorem of Limit Analysis.
The proposed method aims at finding a concave stress potential in equilibrium with the applied loads on a membrane geometry entirely included in the volume of the membrane. To match this objective, the algorithm starts from a parametric description of the stress potential and follows two subsequent optimization steps: in the first one, the parameters of the stress potential are optimized, through the solution of the transverse equilibrium equation of a membrane, assuming fixed boundary conditions for the membrane shape; in the second step, the boundary values of the membrane shape are selected, by minimizing an error function defined as a quadratic approximation of the characteristic function over the planform of the membrane.
The main advantage of this approach is the possibility of solving, in each optimization phase, an elliptic problem, since the stress potential is assumed to be known and concave. This strategy permits the solution of problems for which a concave shape, comprised in the volume of the vault, does not exist, such as in the case of cross vaults. For such vaults, in fact, the search of a stress field starting from a possible membrane shape (in-between the intrados and the extrados of the vault) is described by a non-elliptic second order partial differential equation, that cannot be solved by the simple inversion of the matrix arising from its discretization.
The algorithm is applied to two benchmarks: a cloister vault, with the dimensions of a vault in Palazzo Caracciolo di Avellino, in Naples, and the cross vault of the anterior portico of San Pietro in Vineis, in Anagni (Italy). For both problems, given a suitable parametric description for the stress potential, a membrane shape is found, fulfilling the requirements of Limit Analysis with acceptable geometrical safety factors, demonstrating the validity of the proposed methodology in assessing the stability of vaults.
Future perspectives for this work will include the development of an algorithm for which the stress potential and the membrane shape are optimized simultaneously, under the constraints of the equilibrium and of the Safe Theorem of Limit Analysis.