Out-of-plane local mechanism analysis with finite element method**

The stability check of masonry structures is a debated problem in Italy that poses serious problems for its extensive use. Indeed, the danger of out of plane collapse of masonry walls, which is one of the more challenging to evaluate, is traditionally addressed not using finite element models (FEM). The power of FEM is not properly used and some simplified method are preferred. In this paper the use of the thrust surface is suggested. This concept allows to to evaluate the eccentricity of the membrane stresses using the FEM method. For this purpose a sophisticated, layered, finite element with a no-tension material is used. To model a no-tension material we used the smeared crack method as it is not mesh-dependent and it is well known since the early ’80 in an ASCE Report [1]. The described element has been implemented by the author in the program Nòlian by Softing.


Introduction
The aim of this research is not theoretical but based on the possibility of applying the computational mechanics' methods to the daily work of the designer. The research started with this query regarding the analysis of masonry structures: is it necessary to use complex and not general methods such as the kinematic analysis for out-plane mechanisms? We will show that the thrust line method, applied through Finite Element analysis using layered elements with no tension material can be a practical solution. This article provides some hints in order to suggest a practical solution to a complex and common problem.
The aim of this research is, therefore, to give to structural engineers a complete tool for the analysis of masonry structures.
Beside the stress analysis of the in-plane stresses, there is, in fact, the problem of the out-of-plane mechanism, which is normally solved with kinematic virtual work approach [8]. This implies that there is not a unified method to deal with these two issues.
Indeed, the out-of-plane problem too could be addressed with FEM, and we will try to show that this is possible. The thrust-surface concept is classic, but the original way illustrated in this paper is to address the problem of determining the thrust-surface by a FEM approach through a layered, no-tension finite element. This is an original approach.

The finite element employed for this study
A so called "degenerate shell element" that is a solid element obtained by extrusion of a plane element (Figures 1 and 2) is used. In the degenerate shell element the three dimensional field equations are used and the term "degen-  erate" is used for the degeneration of a continuum to a surface structure . This is to preserve the computational advantages of the single-layer Mindlin-Reissener finite element computation. This approach was first proposed by Ramm et al. [9] and implemented in the last two decades. Moreover this element is a layered one ( Figure 3) so that if a no-tension material is assigned to each layer, the position of the membrane stress resultant is correctly obtained (Figure 4). The finite element was developed for the following reasons: we needed a shell element with high rate of convergence, that could model curved surfaces with varying thickness, i.e. a general element for shell structures that could be used in our "commercial" software, Nòlian. Moreover we needed an element suitable to model a concrete structure so a layered model was the best choice. Layers of reinforcement could be modelled, FRP could be inserted  too and so on, and this with a great easiness. For these reasons we decided for a degenerate shell element. Taking account of the degenerate shell assumption, the displacement field is described by the five degrees of freedom of a normal to the node: three translations of the middle node and two rotations. This model allows, moreover, the layer model for the element.
From a software point of view, we have a so called "material" for each Gauss point of every layer so that the history of displacement and stress is stored into the informations of the "material". A hierarchical formulation has been adopted to implement a 9-nodes Serendipity and Heterosis element. The shape function of the eight boundary nodes are Serendipity shape functions and the hierarchic shape function for the ninth node is a bubble function. A reduced integration is possible and a selective integration can be used.
It should be noted that an appropriate integration through the thickness has to be carried out in order to obtain a reasonable accuracy if the behaviour of the material is non-linear and the constitutive law has to be satisfied through the thickness. This is very important. This problem is solved with the layered schema. The stress component in each layer is considered constant in the thickness of the layer so the stress distribution in the shell is modelled by a piecewise constant approximation. The layer thickness is expressed in curvilinear coordinates so that the layer thickness varies as the shell thickness varies. The stress resultants are obtained by integrating the corresponding stress component in each layer.
The element, over the material non linearity, has a geometrical non-linear capability too, obtained by adding a geometrical stiffness matrix, As said, the material is assigned to the layer so that many materials can be used without modifying the element formulation, the material we used in this work has the aforementioned characteristic.
As transverse shears are taken into account, a triaxial yield criterion must be employed. The criterion is formulated in terms of the first two stress invariants J 1 e J 2 and only two parameters are involved.
The two parameters may be assigned by the user and, in this case, we assume α = 0.355σ 0 and β = 1.355σ 0 from experimental results [7]. Figure 5 illustrates the two dimensional stress dominion. The tensile stress is assumed to be linear elastic until the tensile stress limit is reached, cracks are assumed to be perpendicular to the direction of maximum principal tensile stress. After cracking is occurred, the elastic-  ity modulus and the Poisson's ratio are both reduced to zero. A smeared representation of cracked material is assumed which implies that cracks are not discrete, but distributed across a region and this avoids a "mesh dependent" behaviour [6]. A gradual release of the stress component normal to the cracked plane is adopted in order to avoid "numerical shock" in the stiffness matrix [ Figure 6].

Thrust line and FEM
In the layered element the element stresses are obtained in each Gauss point by a simple integration of the stress contribution of each layer. This method is accurate enough for stress analysis. The exact direction and position of the resultant of the membrane stress, i.e. the thrust line, needs a more accurate procedure. This is due to the fact that each layer is independent and generally the directions of the main stresses are not parallel in each layer. For this reason, we preferred to evaluate the main stress' direction in each layer and find, by projection, the direction of the resultant of the stresses. Then, it is possible to obtain the position of the resultant, i.e., the depth of it in the thickness of the element. It is straight forward that it is possible that the two main stresses are both in compression. In this case we obtain a thrust net and not only a thrust line. In our implementation, and in the images in this paper, we represented the greater value of the two values, if any. As you can see, we do not use the concept of element flexural stress as it is not a layer value but the result of an integration. Shortly, we could assume the main directions of the stresses in the exterior layer, then project membrane and flexural stresses in that direction and obtain the eccentricity of the thrust line by a simple ratio. But, as we said, we discovered that this simple method in many situation is not accurate.

Application to the arc
The problem of the line of thrust has been addressed with FEM for arcs. Shortly we will show that this is possible if a finite element is enough sophisticate to have plastic behaviour with limit to tensile stress and rotating or smeared crack.
We take the data for this example from [11] where the classical theory is applied to determine the minimum thickness of a semi-circular arch subjected to gravity load only. The data for this example are the following: span = 17 m, minimum thickness = 0.90 m, specific weight = 18,000 N/m 3 . The target value for the position of the hinge is 0.6 rad. The value 1.0 is reached (red color). Figure 8 shows the results of an analysis with elastic material. Figure 9 shows instead the same arch but with a material that has non-tension behavior. As you see, the no-tension solution achieves the correct value for the angle

The dome as a benchmark
As a benchmark for the proposed method, we use a hemispherical dome, this as the study of the dome is well documented in literature. A problem arises, using a smeared rotating crack no-tension material because the hoop tensile strength makes the matrix of the element singular. To avoid singularity, we drugged the element matrix for the degree of freedom related to the hoop direction. The problem of cracking in the base of the dome is well known and, in effect vertical cracking transforms the dome from a twodimensional to a one-dimensional structure. Considering that the vertical cracks are typical, the masonry dome is a shell only from the architectural point of view, while it is a  (Figure 11, main stress directions in the dome, in red tension stress, Figure 12, an historical picture of a vertical crack in the dome of the Pantheon in Rome, 1925, during restoration). Nevertheless, we assume the dome as a good benchmark.
As target values we assume the research of the minimal thickness [8] that suggests t/R = 4.32%.
With this thickness, we reach a thrust surface in good accordance with the target values (Figure 13 a section of the dome and the thrust surface and Figure 14 from [5]). In fact the maximum eccentricity of the thrust is about unitary and the shape of the thrust line is what is expected. In Figure  15 a wedge of the dome shows the meaning of the above assertion about a dome as a system of arches.

This method applied to a prestigious architectural case
An interesting study on the Basento bridge [6], designed by Sergio Musmeci and realized in 1971-76,the proposed method was used for a thrust-line analysis to show that in this structure the flexural action is very small. Why this interest for the amount of flexural stress in the bridge? Musmeci was strongly interested in natural shapes. He was worried by the use of arbitrary shapes in architecture and what could ensure that a shape is natural, fluent and not an arbitrary act of human presumption? The natural shapes use the least material possible so the minimum weight structure has a good assumption to be not gratuitous. Certainly a minimum weight structure has no relevant flexural actions.
With the study of the eccentricity of the thrust surface in the bridge, it has been possible to evaluate the truthfulness [ Figure 16] of this attempt. The maximum eccentricity shows, infact, that the membrane resultant is far from the central line of the membrane for less than 1 mm! (analysis of the bridge by the author with Nòlian by Softing).

The corner analysis: an example of a general application
As already stated, the analysis of domes is not the main purpose of our research as our intention is more practical than theoretical and it is not so common for a professional to deal with domes. So, we are now exposing a problem of local mechanism that is among the most tedious to deal with kinematic methods. We will see that good results can be obtained, giving in addition a general and complete image of the behaviour of a masonry structure both for in-plane and out-plane actions.
The classic approach to the out-of-plane mechanism, is based on equilibrium of a part of the structure due to resisting and acting momentum, this means that it is generally accepted the theory of thrust line which is based on the same concept. Therefore, the shear stress limit is not considered and the out-plane mechanism is due to eccentricity of the forces only. The method that we propose is to find, for each Gauss point of the element, the resultant of the compressive stress for each layer. This procedure is accurate and is very similar to the method of obtaining the eccentricity of the resultant through the relation: e = M / N. The stresses considered are the main compressive stresses in the Gauss point for each layer.
Of course, the case of a simple wall is straightforward and elastic or plastic analysis lead to the same results, this happens for isostatic models.
For this reason and for the sake of brevity, we choose the case of the corner of a building that is a more complex problem and that will show that the usual kinematic model (Figure 17, [7]) is perhaps not satisfactory in general cases. Moreover, using a nonlinear analysis with a no-tension finite element, the loss of convergence does not necessary means that a limit state is achieved, but we can monitor the level of damage of each element so that we can detect if a limit state is achieved.
It is useful to point out that, to define the material properties, we do not use a micro-model. The characteristics of the material are instead obtained by an homogenized Figure 18: The trust line eccentricity in a masonry corner continuous: bricks and mortar are homogenized through the Eshelby tensor. Figure 18 shows the eccentricity of the thrust line in color. The value reached is 0.6181 that means a collapse multiplier of 1 − e = 0.31. With the kinematic model (Figure 16, Milani et al.) the value, computed with the geometric shape obtained by the FEM model, is 0.35.

Conclusion
This method shows the analyst the real behaviour of the masonry structure even for very complex shapes. And, what is important for the analyst, a single method covers all the aspects of the behaviour of the structure.
Naturally, it is clear that in a FEM model all the additional items can be modelled as well as contact, friction between walls, rods and so on, and this is impossible in simplified models. It is also obvious that this method can be used in dynamic non-linear analysis simply monitoring and recording at each step the value of eccentricity that exceeds a predetermined value. In this way we can also obtain, in a simple and rational way, also the so called "dynamic kinematic multiplier".

Funding information:
The authors state no funding involved.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.

Conflict of interest:
The authors state no conflict of interest.