Next Article in Journal
Review of Renewable Energy-Based Charging Infrastructure for Electric Vehicles
Next Article in Special Issue
Brunelleschi’s Dome: A New Estimate of the Thrust and Stresses in the Underlying Piers
Previous Article in Journal
Geometrical Degrees of Freedom for Cellular Structures Generation: A New Classification Paradigm
Previous Article in Special Issue
Graphical and Analytical Quantitative Comparison in the Domes Assessment: The Case of San Francesco di Paola
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

From Stress to Shape: Equilibrium of Cloister and Cross Vaults

by
Andrea Montanino
1,*,†,
Carlo Olivieri
2,†,
Giulio Zuccaro
1,† and
Maurizio Angelillo
2,†
1
Department of Structures for Engineering and Architecture, University of Naples “Federico II”, Via Claudio 21, 80143 Napoli, Italy
2
Department of Civil Engineering, DICIV, University of Salerno, Via Giovanni Paolo II 132, 84084 Fisciano, Italy
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Appl. Sci. 2021, 11(9), 3846; https://doi.org/10.3390/app11093846
Submission received: 20 March 2021 / Revised: 21 April 2021 / Accepted: 21 April 2021 / Published: 23 April 2021
(This article belongs to the Special Issue Theory and Modelling of Historic Masonry Architecture)

Abstract

:
The assessment of the equilibrium and the safety of masonry vaults is of high relevance for the conservation and restoration of historical heritage. In the literature many approaches have been proposed for this tasks, starting from the 17th century. In this work we focus on the Membrane Equilibrium Analysis, developed under the Heyman’s theory of Limit Analysis. Within this theory, the equilibrium of a vault is assessed if it is possible to find at least one membrane surface, between the volume of the vaults, being in equilibrium under the given loads through a purely compressive stress field. The equilibrium of membranes is described by a second order partial differential equation, which is definitely elliptic only when a negative semidefinite stress is assigned, and the shape is the unknown of the problem. The proposed algorithm aims at finding membrane shapes, entirely comprised between the geometry of the vault, in equilibrium with admissible stress fields, through the minimization of an error function with respect to shape parameters of the stress potential, and then, with respect to the boundary values of the membrane shape. The application to two test cases shows the viability of this tool for the assessment of the equilibrium of existing masonry vaults.

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]

2. Mathematical Preliminaries

In this section some mathematical details concerning the formulation of the equilibrium equation of a membrane are recalled, the notation is introduced, and the unilateral constraints on the material are introduced.

2.1. Geometry and Equilibrium

Referring to Figure 1 for notation, the geometry of a generic membrane surface S contained inside a typical shell structure can be described in Monge form in terms of the position vector x of the points of S as follows:
x = x 1 e 1 + x 2 e 2 + f e 3
where { e 1 , e 2 , e 3 } is the orthonormal triad coherent with a given Cartesian frame { O ; x 1 , x 2 , x 3 } , the couple of coordinates ( x 1 , x 2 ) , which are the parameters of the curvilinear description (1), belongs to a region Ω of the plane { O ; x 1 , x 2 } called planform of S, and f = f ( x 1 , x 2 ) is the rise of the membrane with respect to the x 1 x 2 plane. We assume f C 0 ( Ω ) —that is to say that the considered surface S is continuous but not necessarily smooth.
An efficient way to describe the membrane equilibrium of a membrane under a surface load q per unit area, is the formulation proposed by Pucher in [38]. Pucher analysis is based on the introduction of the projected stress S , in terms of which two of the three equilibrium equations for the membrane can be made independent of the membrane shape.
The membrane stresses are tangent in each point to the membrane surface S. Therefore it is convenient to introduce a curvilinear basis on the membrane, which are computed as
a α = x x α , α = 1 , 2
In addition, it is useful to introduce the contravariant basis a α . α = 1 , 2 defined such that a α · a β = δ α β , being δ α β the Kronecker delta.
Therefore we can decompose the membrane generalized stress tensor as
T = T α β a α a β ,
where T α β are the contravariant components of the generalized stress. The projected (Pucher) stress components S α β (see Figure 1) are easily defined in terms of T α β as S α β = J T α β , J being the ratio between the surface area and the projected area. In terms of Pucher stress components, the vector equilibrium equation, projected into the non-orthonormal vector basis e 1 , e 2 , n = a 3 , n being the unit normal to S, becomes:
S 11 , 1 + S 12 , 2 + p 1 = 0 , S 21 , 1 + S 22 , 2 + p 2 = 0 , S α β f , α β p γ f , γ + p 3 = 0 .
where { p 1 , p 2 , p 3 } are the Cartesian components in the given reference of p = J q , that is the load per unit projected area. By using the projected stress, the first two equilibrium equations do not depend on the membrane shape and are the same as for the plane stress problem. In the case of pure vertical loading, say p = p e 3 , the first two equations can be solved introducing an Airy stress potential F ( x 1 , x 2 ) —that here we only assume to be continuous—defined such that
S 11 = F , 22 , S 22 = F , 11 , S 12 = S 12 = F , 12 .
The third equation of (4) expresses the balance of the vertical component of the force p 3 = p · e 3 with the scalar product between the Pucher stress tensor and the Hessian of the function f, that defines the curvature of the surface. On introducing (5) into (4), we obtiain, as transversal equilibrium equation
F , 11 f , 22 + F , 22 f , 11 2 F , 12 f , 12 + p 3 = 0
Equation (6) is a second order partial differential equation, that can be either elliptic, parabolic or hyperbolic, depending on its (usually non-constant) coefficients. If the problem is hyperbolic or parabolic, the characteristic lines are real and the problem must be formulated as an evolutionary problems by assigning two scalar conditions on a non-characteristic line.
If Equation (6) is elliptic, the problem can be formulated as a boundary value problem and the boundary conditions for this equation can be either of the Dirichlet type, i.e.,
F | Ω = p
or Neumann type, i.e.,
F n | Ω = q
where p and q are known function on the boundary Ω , and n is the outward unit vector normal to the boundary.
The advantage of boundary value problems is that the solution of the corresponding discretized partial differential equation can be obtained by a simple matrix inversion.

2.2. Unilateral Constrains and Singular Stresses

The restrictions that the generalized stress T is negative semi-definite and does not work for the corresponding strain E , that is positive semi-definite, define Rigid No-Tension (RNT) materials in the sense of Heyman. In formulae:
T S y m , E S y m + , T · E = 0 .
The first application of Pucher’s transformation for NT masonry vaults can be found in [18], where it is shown that, due to the NT constraint, both the surface generalized stress and the matrix of the projected stresses must be negative semi-definite. In terms of the stress function F, this condition can be written as
F , 11 + F , 22 0 , F , 11 F , 22 F , 12 2 0 .
The condition of semi-definiteness of the stress correspond to concavity of the stress potential.

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 F ^ i and f ^ i , with i = 1 , 2 , , N , 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 f ^ and F ^ . f ^ is divided into two vectors f ^ I and f ^ B , of length N I and N B , 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
F = F ( x 1 , x 2 ; a 1 , a 2 , , a n ) .
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, f = f ( x 1 , x 2 ; a 1 , a 2 , , a n ) .
The second step consists in the optimization of the parameters a 1 , a 2 , , a n , 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 x Ω , 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
C ( f ( x ) ) = 0 if f i n t r ( x ) < f ( x ) < f e x t r ( x ) + if f ( x ) < f i n t r ( x ) or f ( x ) > f e x t r ( x ) .
Its quadratic approximation reads
C ( f ( x ) ) = 0 f i n t r ( x ) < f ( x ) < f e x t r ( x ) f ( x ) f i n t r ( x ) 2 f ( x ) < f i n t r ( x ) f ( x ) f e x t r ( x ) 2 f ( x ) > f e x t r ( x ) .
For a generic point x of the domain, given f i n t r ( x ) and f e x t r ( x ) , the approximated characteristic function is represented in Figure 2.
We denote Ω c 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 c = L / 5 ). 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:
I = Ω a ( f f i n t r ) 2 d Ω + Ω b ( f f e x t r ) 2 d Ω
Ω a and Ω b being the parts of the domain Ω c , for which the membrane surface is below the intrados, or above the extrados, respectively.
In the third step, the values of a 1 , a 2 , , a n , 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
( K x x , H j f ^ j ) ( K y y , H h F ^ h ) + ( K y y , H j f ^ j ) K x x , H h F ^ h ) 2 ( K x y , H j f ^ j ) ( K x y , H h F ^ h ) + p ^ H = 0
where the summation convention on repeated indices is assumed on the lower-case indices ( j , h = 1 , 2 , , N ) , whilst the upper-case index H ( H = 1 , 2 , , N I ) is not summed, and K x x , K y y , and K x y are coefficient matrices computed via the LSM scheme, representing the discrete counterparts of the second derivatives 2 / x 1 2 , 2 / x 2 2 , and 2 / x 1 x 2 .
Equation (15) can be rewritten in a compact form by factoring out the unknown term f ^ , as
K i j f ^ j + p ^ i = 0 ,
where
K i j = ( K y y , H h F ^ h ) K x x , H j + ( K x x , H h F ^ h ) K y y , H j 2 ( K x y , H h F ^ h ) K x y , H j .
Expression (16) can be also rewritten distinguishing between boundary conditions and domain equations, in the form
K I I f ^ I + K I B f ^ B + p ^ I = 0 , f ^ B = f ¯ B ,
where K I I and K I B are submatrices of K corresponding to rows and columns associated to internal and boundary nodes. Finally, f ¯ B is the vector of the boundary parameters. From Equation (18), f ¯ I can be evaluated as
f ^ I = K I I 1 ( p ^ I K I B f ¯ B ) .
On introducing (19) into a discrete version I 2 of the objective function I, we get
I 2 = i = 1 N a ( f i f i n t r , i ) 2 + i = 1 N b ( f i f e x t r , i ) 2 .
We then optimize I 2 ( f ^ I ) with respect to f ¯ B .
Notice that it is possible to compute explicitly the gradient of I 2 ( f I ( f B ) ) through the chain rule, in the form
I 2 f B , j = I f I , h f I , h f B , j
that is
I 2 f B , j = 2 f i , h ( K I I K I B ) h j f i < f i , i n t r or f i > f i , e x t r 0 otherwise .
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 I 2 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 Ω = { ( x 1 , x 2 ) L , L × L , L } as follows:
f i n t r ( x 1 , x 2 ) = h x 1 L 2 ( x 1 , x 2 ) ( Ω 1 Ω 3 ) h x 2 L 2 ( x 1 , x 2 ) ( Ω 2 Ω 4 )
f e x t r ( x 1 , x 2 ) = h x 1 L + t 2 + t ( x 1 , x 2 ) ( Ω 1 Ω 3 ) h x 2 L + t 2 + t ( x 1 , x 2 ) ( Ω 2 Ω 4 )
where Ω 1 , Ω 2 , Ω 3 , and Ω 4 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 L = 2.5 m, h = 2.13 m, and t = 0.23 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 p = 8 kN/m 2 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 f ( x 1 , x 2 ) , 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
r c ( θ ) = r 0 2 ( 1 α ( c o s ( 2 θ ) ) P )
where r 0 is a parameter, α is set to 0.05 and P to 40. In Figure 5 we show the contour r c obtained for two representative values of r 0 , namely r 0 = 0.5 and r 0 = 1 .
Inside Ω * , the Pucher stress is assumed as isotropic, that is
S = p 0 ( e 1 e 1 + e 2 e 2 ) .
The corresponding stress potential is therefore easily computed as
F = 1 2 p 0 ( x 1 2 + x 2 2 ) ,
where p 0 is a parameter defining the intensity of the stress potential (i.e., of the stress) that will be subsequently optimized. In the external zone Ω e x t the stress is assumed to be quasi-uniaxial in the radial direction. This means that, on defining r ^ as the direction of the radii emanating from the geometrical center of the domain and directed toward the external boundary, the stress in Ω e x t can be written as
S = σ ( s ) r ^ r ^ + ε σ 0 t ^ t ^
where σ ( s ) is related to p 0 and denotes the stress intensity along the radius, and ε σ 0 is a small stress which has the function to give a small extra-curvature to the potential in the direction t ^ orthogonal to r ^ , in order to preserve the ellipticity of Equation (6) when solved for the shape f ( x 1 , x 2 ) .
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 Γ :
F e x t = F i n t | Γ · r ^ s + F | Γ ε σ 0 2 s 2 ,
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 r 0 . Adopting these two shapes as the stress potential, and taking as boundary condition for the membrane geometry f ( x 1 , x 2 ) 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 r 0 (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 ( r 0 = 1 ), the membrane geometry resembles the shape of a cloister vault more than in the case of r 0 = 0.5 . It is also evident that the height of the two shapes (in the direction of the x 3 -axis) are not comparable with the height of the cloister vault defined from Equations (23) and (24), which is 2.13 m.
The minimization of the error function (20) with respect to the parameters p 0 and r 0 is performed as follows: first, r 0 is given a finite number of values in the range [ 0 , L ] , i.e., r 0 , k = 1.1 + k 0.01 , k = 0 , , 20 . Then, for each value of r 0 , k , we search the p 0 , k that minimizes I k = I k ( r 0 , k , p 0 , k ) . In Figure 8 we show the considered values of the r 0 , k on the horizontal axis, and on the vertical axes we show the corresponding optimal value of I k .
The lowest value of I k is reached when r k = 1.14 , which correspond to p 0 , k = 3.45 × 10 3 . Therefore the parameters of the stress potential are fixed to r 0 = 1.14 and p 0 = 3.45 × 10 3 . 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 x 2 = 0 , and the second one taken along the diagonal of the domain, that is for x 2 = x 1 .
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 g ( x 1 , x 2 ) defined as
g ( x 1 , x 2 ) = f f i n t r f < f i n t r 0 f i n t r < f < f e x t r f f e x t r f > f e x t r
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 x 2 = 0 and x 2 = x 1 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 38 % of the original thickness, corresponding to a geometrical safety factor of 1.61 .

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:
f i n t r ( x 1 , x 2 ) = 4 h i , 1 b i 2 b i 2 4 x 2 2 + 4 h i h i , 1 a i 2 a i 2 4 x 1 2 ( x 1 , x 2 ) ( Ω 1 Ω 3 ) 4 h i , 2 b i 2 b i 2 4 x 2 2 + 4 h i h i , 2 a i 2 a i 2 4 x 1 2 ( x 1 , x 2 ) ( Ω 2 Ω 4 )
f e x t r ( x 1 , x 2 ) = 4 h e , 1 b e 2 b e 2 4 x 2 2 + 4 h e h e , 1 a e 2 a e 2 4 x 1 2 ( x 1 , x 2 ) ( Ω 1 Ω 3 ) 4 h e , 2 b e 2 b e 2 4 x 2 2 + 4 h e h e , 2 a e 2 a e 2 4 x 1 2 ( x 1 , x 2 ) ( Ω 2 Ω 4 )
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 F ( x 1 , x 2 )
F ( x 1 , x 2 ) = α 0 x 1 2 a i 2 4 + β x 2 2 b i 2 4 ( x 1 , x 2 ) ( Ω 1 Ω 3 ) α 0 x 2 2 b i 2 4 + β x 1 2 a i 2 4 ( x 1 , x 2 ) ( Ω 2 Ω 4 )
For sake of example, in Figure 15 we present three different stress potentials obtained for β = 0 , β = 0.3 , and β = 0.7 and α 0 = 1 .
We assume the same definition of the error introduced in Equation (20) and evaluate it for different values of β , and then minimizing I k with respect to α 0 . The values of β considered in this case are
β k = 0.3 + k 0.02 , k = 0 , . . . , 10
The corresponding plot of the error I k versus β is reported in Figure 16 in a logarithmic scale along the error axis.
The minimum occurs at β = 0.4 , and corresponds to an error of 3.94 × 10 2 . 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 f ( x 1 , x 2 ) 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 g ( x 1 , x 2 ) , 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 I = I ( f ¯ B ) , 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 x 2 = b i a i x 1 , 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 f ¯ B is made evident by the color map of the difference h ( x 1 , x 2 ) defined as
h ( x 1 , x 2 ) = f o p t ( x 1 , x 2 ) f ( x 1 , x 2 )
and shown in Figure 21, where f is the membrane geometry at the end of the first minimization step, and f o p t 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 h ( x 1 , x 2 ) 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 36 % , corresponding to a GSF = 1.56 .

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.

Author Contributions

Conceptualization, A.M., M.A.; methodology, A.M., M.A.; software, A.M.; validation, A.M., G.Z., M.A.; formal analysis, A.M.; data curation, C.O.; writing—original draft preparation, A.M., C.O.; writing—review and editing, G.Z., M.A.; supervision, G.Z., M.A.; funding acquisition, M.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research has been partially supported by the “AIM: Attraction and International Mobility”, PON R&I 2014–2020 Campania, n. 1849854–3 “Smart Secure and Inclusive Communities”. CUP: E66C19000230005.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Derand, F. L’architecture des Voutes; Sébastien Cramoisy: París, France, 1643. [Google Scholar]
  2. Huerta, S. The use of simple models in the teaching of the essentials of masonry arch behaviour. In Theory and Practice of Construction: Knowledge, Means, and Models; Mochi, G., Ed.; Ravenna, Italy, 2005; pp. 747–761. Available online: https://core.ac.uk/download/pdf/148651963.pdf (accessed on 20 March 2021).
  3. Block, P.; Ochsendorf, J. Lower-bound analysis of masonry vaults. In Structural Analysis of Historic Construction; Taylor & Francis Group: London, UK, 2008; pp. 593–600. [Google Scholar]
  4. Poleni, G. Memorie Istoriche Della Gran Cvpola del Tempio Vaticano, E De’Danni di Essa, E De’Ristoramenti Loro, Divise in Libri Cinqve; Nella Stamperia del Seminario: Padova, Italy, 1748. [Google Scholar]
  5. Lau, W.W. Equilibrium Analysis of Masonry Domes. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA, USA, 2006. [Google Scholar]
  6. Heyman, J. The stone skeleton. Int. J. Solids Struct. 1966, 2, 249–279. [Google Scholar] [CrossRef]
  7. Heyman, J. The Stone Skeleton: Structural Engineering of Masonry Architecture; Cambridge University Press: Cambridge, UK, 1997. [Google Scholar]
  8. Iannuzzo, A. Energy based fracture identification in masonry structures: The case study of the church of “Pietà dei Turchini”. J. Mech. Mater. Struct. 2019, 14, 683–702. [Google Scholar] [CrossRef]
  9. Iannuzzo, A.; Olivieri, C.; Fortunato, A. Displacement capacity of masonry structures under horizontal actions via PRD method. J. Mech. Mater. Struct. 2019, 14, 703–718. [Google Scholar] [CrossRef]
  10. Iannuzzo, A.; Van Mele, T.; Block, P. Piecewise rigid displacement (PRD) method: A limit analysis-based approach to detect mechanisms and internal forces through two dual energy criteria. Mech. Res. Commun. 2020, 107, 103557. [Google Scholar] [CrossRef]
  11. O’Dwyer, D. Funicular analysis of masonry vaults. Comput. Struct. 1999, 73, 187–197. [Google Scholar] [CrossRef]
  12. Block, P.; Ochsendorf, J. Thrust network analysis: A new methodology for three-dimensional equilibrium. J. Int. Assoc. Shell Spat. Struct. 2007, 48, 167–173. [Google Scholar]
  13. Block, P. Thrust Network Analysis: Exploring Three-Dimensional Equilibrium. Ph.D. Thesis, Massachusetts Institute of Technology, Cambridge, MA, USA, 2009. [Google Scholar]
  14. Block, P.; Lachauer, L.; Rippmann, M. Thrust network analysis. Shell Structures for Architecture: Form Finding and Optimization; Routledge: London, UK; New York, NY, USA, 2014; Volume 71. [Google Scholar]
  15. De Goes, F.; Alliez, P.; Owhadi, H.; Desbrun, M. On the equilibrium of simplicial masonry structures. ACM Trans. Graph. 2013, 32, 1–10. [Google Scholar] [CrossRef]
  16. Marmo, F.; Rosati, L. Reformulation and extension of the thrust network analysis. Comput. Struct. 2017, 182, 104–118. [Google Scholar] [CrossRef]
  17. Maia Avelino, R.; Iannuzzo, A.; Van Mele, T.; Block, P. Parametric Stability Analysis of Groin Vaults. Appl. Sci. 2021, 11, 3560. [Google Scholar] [CrossRef]
  18. Angelillo, M.; Fortunato, A. Equilibrium of masonry vaults. In Novel Approaches in Civil Engineering; Springer: Berlin/Heisenberg, Germany, 2004; pp. 105–111. [Google Scholar]
  19. Angelillo, M.; Babilio, E.; Fortunato, A. Singular stress fields for masonry-like vaults. Contin. Mech. Thermodyn. 2013, 25, 423–441. [Google Scholar] [CrossRef]
  20. Gesualdo, A.; Cennamo, C.; Fortunato, A.; Frunzio, G.; Monaco, M.; Angelillo, M. Equilibrium formulation of masonry helical stairs. Meccanica 2017, 52, 1963–1974. [Google Scholar] [CrossRef]
  21. Gesualdo, A.; Brandonisio, G.; De Luca, A.; Iannuzzo, A.; Montanino, A.; Olivieri, C. Limit analysis of cloister vaults: The case study of Palazzo Caracciolo di Avellino. J. Mech. Mater. Struct. 2019, 14, 739–750. [Google Scholar] [CrossRef]
  22. Cennamo, C.; Cusano, C.; Angelillo, M. A limit analysis approach for masonry domes: The basilica of San Francesco di Paola in Naples. Int. J. Mason. Res. Innov. 2019, 4, 227–242. [Google Scholar] [CrossRef]
  23. Cusano, C.; Montanino, A.; Zuccaro, G.; Cennamo, C. Considerations about the static response of masonry domes: A comparison between limit analysis and finite element method. Int. J. Mason. Res. Innov. 2021, in press. [Google Scholar]
  24. Cusano, C.; Montanino, A.; Olivieri, C.; Paris, V.; Cennamo, C. Graphical and Analytical Quantitative Comparison in the Domes Assessment: The Case of San Francesco di Paola. Appl. Sci. 2021, 11, 3622. [Google Scholar] [CrossRef]
  25. Angelillo, M.; Montanino, A.; Pandolfi, A. On the Connection Between Geometry and Statically Determined Membrane Stresses in the Human Cornea. J. Biomech. Eng. 2020, 142, 051006. [Google Scholar] [CrossRef] [PubMed]
  26. Miki, M.; Igarashi, T.; Block, P. Parametric self-supporting surfaces via direct computation of airy stress functions. ACM Trans. Graph. 2015, 34, 1–12. [Google Scholar] [CrossRef]
  27. De Chiara, E.; Cennamo, C.; Gesualdo, A.; Montanino, A.; Olivieri, C.; Fortunato, A. Automatic generation of statically admissible stress fields in masonry vaults. J. Mech. Mater. Struct. 2019, 14, 719–737. [Google Scholar] [CrossRef]
  28. Fraternali, F.; Angelillo, M.; Fortunato, A. A lumped stress method for plane elastic problems and the discrete-continuum approximation. Int. J. Solids Struct. 2002, 39, 6211–6240. [Google Scholar] [CrossRef]
  29. Olivieri, C.; Angelillo, M.; Gesualdo, A.; Iannuzzo, A.; Fortunato, A. Parametric design of purely compressed shells. Mech. Mater. 2021, 155, 103782. [Google Scholar] [CrossRef]
  30. Bendsøe, M.P. Optimal shape design as a material distribution problem. Struct. Optim. 1989, 1, 193–202. [Google Scholar] [CrossRef]
  31. Sigmund, O. A 99 line topology optimization code written in Matlab. Struct. Multidiscip. Optim. 2001, 21, 120–127. [Google Scholar] [CrossRef]
  32. Montanino, A.; Alaimo, G.; Lanzarone, E. A gradient-based optimization method with functional principal component analysis for efficient structural topology optimization. Struct. Multidiscip. Optim. 2021, 1–12. [Google Scholar] [CrossRef]
  33. Bruggi, M. On an alternative approach to stress constraints relaxation in topology optimization. Struct. Multidiscip. Optim. 2008, 36, 125–141. [Google Scholar] [CrossRef]
  34. Nguyen, T.; Ghabraie, K.; Tran-Cong, T. Simultaneous pattern and size optimisation of rock bolts for underground excavations. Comput. Geotech. 2015, 66, 264–277. [Google Scholar] [CrossRef]
  35. Pilarska, D. Two subdivision methods based on the regular octahedron for single-and double-layer spherical geodesic domes. Int. J. Space Struct. 2020, 35, 160–173. [Google Scholar] [CrossRef]
  36. Sobótka, M.; Łydżba, D. Shape optimization of soil-steel structure by simulated annealing. Procedia Eng. 2014, 91, 304–309. [Google Scholar] [CrossRef] [Green Version]
  37. Chen, Y.C.; Huang, B.K.; You, Z.T.; Chan, C.Y.; Huang, T.M. Optimization of lightweight structure and supporting bipod flexure for a space mirror. Appl. Opt. 2016, 55, 10382–10391. [Google Scholar] [CrossRef]
  38. Pucher, A. Uber der spannungzustand in gekrummten flachen. Beton Eisen 1934, 33, 298–304. [Google Scholar]
Figure 1. Geometry of a generic membrane, with its planform, the memebrane stress T and the Pucher (projected) stress S .
Figure 1. Geometry of a generic membrane, with its planform, the memebrane stress T and the Pucher (projected) stress S .
Applsci 11 03846 g001
Figure 2. Approximation of the characteristic function for a generic point x .
Figure 2. Approximation of the characteristic function for a generic point x .
Applsci 11 03846 g002
Figure 3. Partition of Ω into four portions.
Figure 3. Partition of Ω into four portions.
Applsci 11 03846 g003
Figure 4. Three-dimensional view of the midsurface of the vault described by Equations (23) and (24).
Figure 4. Three-dimensional view of the midsurface of the vault described by Equations (23) and (24).
Applsci 11 03846 g004
Figure 5. Representation of the curve Γ for r 0 = 0.5 and r 0 = 1 .
Figure 5. Representation of the curve Γ for r 0 = 0.5 and r 0 = 1 .
Applsci 11 03846 g005
Figure 6. Representation of the stress potential assuming (a) r 0 = 0.5 and p 0 = 1 ; (b) r 0 = 1 and p 0 = 1 .
Figure 6. Representation of the stress potential assuming (a) r 0 = 0.5 and p 0 = 1 ; (b) r 0 = 1 and p 0 = 1 .
Applsci 11 03846 g006
Figure 7. Membrane surfaces in equilibrium with the stress potentialof Figure 6: (a) r 0 = 0.5 and p 0 = 1 ; (b) r 0 = 1 and p 0 = 1 .
Figure 7. Membrane surfaces in equilibrium with the stress potentialof Figure 6: (a) r 0 = 0.5 and p 0 = 1 ; (b) r 0 = 1 and p 0 = 1 .
Applsci 11 03846 g007
Figure 8. Plot of I k as function of r k .
Figure 8. Plot of I k as function of r k .
Applsci 11 03846 g008
Figure 9. Best fitting membrane shape after the first optimization phase, for the problem of the cloister vault.
Figure 9. Best fitting membrane shape after the first optimization phase, for the problem of the cloister vault.
Applsci 11 03846 g009
Figure 10. Sections of the vault and of the membrane surface represented in Figure 9 along x 2 = 0 (a) and along the diagonal x 2 = x 1 (b), after optimization with respect to the shape of the stress potential.
Figure 10. Sections of the vault and of the membrane surface represented in Figure 9 along x 2 = 0 (a) and along the diagonal x 2 = x 1 (b), after optimization with respect to the shape of the stress potential.
Applsci 11 03846 g010
Figure 11. Colormap of g with respect to x 1 and x 2 . Yellow points represent regions of the domain for which the membrane is above the extrados. Blue points on the diagonal of the domain represent regions of the domain for which the membrane shape is below the intrados.
Figure 11. Colormap of g with respect to x 1 and x 2 . Yellow points represent regions of the domain for which the membrane is above the extrados. Blue points on the diagonal of the domain represent regions of the domain for which the membrane shape is below the intrados.
Applsci 11 03846 g011
Figure 12. Sections of the vault and of the membrane surface along x 2 = 0 (a) and along the diagonal x 2 = x 1 (b), after optimization with respect to the membrane boundary values.
Figure 12. Sections of the vault and of the membrane surface along x 2 = 0 (a) and along the diagonal x 2 = x 1 (b), after optimization with respect to the membrane boundary values.
Applsci 11 03846 g012
Figure 13. Parition of the domain Ω .
Figure 13. Parition of the domain Ω .
Applsci 11 03846 g013
Figure 14. Perspective of the cross vault with main geometrical parameters.
Figure 14. Perspective of the cross vault with main geometrical parameters.
Applsci 11 03846 g014
Figure 15. Stress potential for (a) β = 0 , (b) β = 0.3 , and (c) β = 0.7 .
Figure 15. Stress potential for (a) β = 0 , (b) β = 0.3 , and (c) β = 0.7 .
Applsci 11 03846 g015
Figure 16. Plot of I k as function of β .
Figure 16. Plot of I k as function of β .
Applsci 11 03846 g016
Figure 17. Optimized stress potential.
Figure 17. Optimized stress potential.
Applsci 11 03846 g017
Figure 18. Membrane surface in equilibrium with the stress potential of Figure 17.
Figure 18. Membrane surface in equilibrium with the stress potential of Figure 17.
Applsci 11 03846 g018
Figure 19. Colormap of g ( x 1 , x 2 ) for the case of the cross vault, after the first optimization phase. Yellow points represent regions of the domain for which the membrane is above the extrados. Blue points on the diagonal of the domain represent regions of the domain for which the membrane shape is below the intrados.
Figure 19. Colormap of g ( x 1 , x 2 ) for the case of the cross vault, after the first optimization phase. Yellow points represent regions of the domain for which the membrane is above the extrados. Blue points on the diagonal of the domain represent regions of the domain for which the membrane shape is below the intrados.
Applsci 11 03846 g019
Figure 20. Section of the membrane and of the limit surfaces of the cross vault along the diagonal x 2 = b i / a i x .
Figure 20. Section of the membrane and of the limit surfaces of the cross vault along the diagonal x 2 = b i / a i x .
Applsci 11 03846 g020
Figure 21. Effects of the change of boundary conditions on the shape of the membrane.
Figure 21. Effects of the change of boundary conditions on the shape of the membrane.
Applsci 11 03846 g021
Table 1. Intrados and extrados data referring to notation in Figure 14.
Table 1. Intrados and extrados data referring to notation in Figure 14.
a i [m] b i [m] h i [m] h i , 1 [m] h i , 2 [m]
4.684.592.161.842.05
a e [m] b e [m] h e [m] h e , 1 [m] h e , 2 [m]
4.884.792.362.042.25
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Montanino, A.; Olivieri, C.; Zuccaro, G.; Angelillo, M. From Stress to Shape: Equilibrium of Cloister and Cross Vaults. Appl. Sci. 2021, 11, 3846. https://doi.org/10.3390/app11093846

AMA Style

Montanino A, Olivieri C, Zuccaro G, Angelillo M. From Stress to Shape: Equilibrium of Cloister and Cross Vaults. Applied Sciences. 2021; 11(9):3846. https://doi.org/10.3390/app11093846

Chicago/Turabian Style

Montanino, Andrea, Carlo Olivieri, Giulio Zuccaro, and Maurizio Angelillo. 2021. "From Stress to Shape: Equilibrium of Cloister and Cross Vaults" Applied Sciences 11, no. 9: 3846. https://doi.org/10.3390/app11093846

APA Style

Montanino, A., Olivieri, C., Zuccaro, G., & Angelillo, M. (2021). From Stress to Shape: Equilibrium of Cloister and Cross Vaults. Applied Sciences, 11(9), 3846. https://doi.org/10.3390/app11093846

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop