An Enhanced Beam Model for the Analysis of Masonry Walls

Background: Some typologies of masonry constructions (e.g. towers or walls with openings) can be reasonably studied through simple beam or frame-like models. For these structures, shear mechanisms often play an important role inducing failure and collapse. Objective: The paper presents an enriched beam model for studying the in-plane response of masonry walls. Initially formulated for masonry columns, towers and masonry slender structures in general, the model is now modified in order to also capture the shear failure mechanisms, in addition to the flexural ones. Methods: Starting with a one-dimensional no-tension model, a strength domain in the plane of the axial and tangential stress of the beam has been added, which has been defined by limiting both the stress shear component with respect to any possible direction and the main compressive stress. Results: The model, implemented in the FEM computational code MADY, allows for short computational times in studying the response of single panels as well as walls with openings. Conclusion: Comparisons with some experimental results from literature and some numerical results from more refined 2D models show the effectiveness and accuracy of the model’s predictions in terms of global and local response.


INTRODUCTION
The issue of modelling masonry structures is still an open one in structural engineering. This is due to the complexity of the material's mechanical behaviour together with the great variety of structural types and constituents, that is mortar and bricks, and their textures. Meanwhile, seismic assessment, repair and strengthening of ancient masonry buildings is a priority, as such structures represent a large part of the cultural heritage of the Mediterranean basin. In particular, masonry's lateral load capacity is crucial in seismic areas, where old masonry buildings have often exhibited high vulnerability to earthquakes. literature [1,2]. A thorough discussion concerning the features and benefits of using such models in comparison with others can be found in a few studies [3 -5].
A large number of models following the equivalent-frame method have been proposed in the literature over the last decade [6 -10]. Most of these make use of macro-elements to represent piers and spandrels. Often, they are developed assuming the classical hypotheses of a beam and make use of a limited number of nodes. Furthermore, they generally rely on a phenomenological description of the constitutive behavior of the panel, especially for shear resistance behavior. The Interested readers may refer to another study [11] for a comprehensive state-of-the-art review.
This paper presents an advanced one-dimensional continuum model for non-linear static analysis of masonry shear walls.
In the previous research by the authors, a numerical model was defined, which performed non-linear static and dynamic analyses on masonry columns, arches, towers, and slender structures in general. The model relies on a costitutive equation formulated in terms of generalized strain and stress for beams with various cross-sections, under the hypotheses that the material is non-linear elastic, unable to withstand tension with limited compressive strength [12,14]. Furthermore, to describe the cyclical behavior and seismic response of masonry elements more realistically, a damage process under compression has been introduced [15]. The model has been implemented in the FEM computational code MADY and applied to the study of several slender masonry structures, which typically exhibit flexural behavior and consequent failure mechanism [16,17].
Indeed, in the past one-dimensional no-tension models have been used extensively to study problems of the stability of masonry columns [18 -20], the static and dynamic response of unreinforced and reinforced arches and vaults [21 -23], as well as the out-of-plane behavior of masonry panels and walls. The same does not hold for the study of the in-plane behavior of masonry panels under lateral actions. In this case, as is widely known, panels are likely to undergo various distinct failure mechanisms, or a combination of some of them: for flexural behaviour (rocking/toe-crushing failure), as well as for shear with diagonal cracking or sliding along bed-joints [3,24], which cannot be modelled using classical no-tension beam models. The proposed constitutive model has been enhanced to address shear behaviour and account for failure with diagonal cracking in masonry panels. To this aim, the beam is formulated according to the Timoshenko theory to account for shear deformations; then, a constitutive law is defined for the shear force, in addition to that proposed for the axial one and bending moment.
Essentially, the constitutive model is characterized by a strength domain in the plane of the axial and tangential stresses of the beam, which has been defined by limiting the stress shear component with respect to any possible plane of the panel; namely, not only with respect to the one orthogonal to the beam axis. Furthermore, in analogy to the proposal for the shear stress, the principal compressive stress, rather than the axial stress component of the beam, has been limited to achieving maximum strength.
In the following, the proposed model is first presented, and then some numerical results are provided and compared to some experimental data available in the literature.
Then, some numerical analyses are conducted in order to highlight the potential of the model to approximate the 2D solution. To this end, some results obtained through the proposed models are compared with those obtained via the commercial code ANSYS [25] and via the formulated bidimensional masonry-like model [26].
Lasly, preliminary results on an idealized wall with a centered opening are provided as an example for practical applications.

MODEL SPECIFICATIONS
In order to introduce a constitutive equation also for shear behaviour, in addition to axial/flexural one, shear deformations are accounted for according to the Timoshenko theory (see Section 2.1).
Then, a nonlinear-elastic constitutive model is assumed for describing masonry's response. Specifically, a strength domain in the plane of the stress components of a 2D beam (σ z -τ zy ) has been defined, as explained in detail in Section 6.
Given the kinematic constrains and the material constitutive laws, the patterns of the axial stress σ z and the tangential stress τ zy over the beam's cross-section can thus be determined for any strain pattern, as shown by the examples provided in Section 8. Finally, by simply integrating these stresses, a constitutive relation which associates the generalized stresses axial force N, shear V and bending moment M-to the generalized strains -ε, κ and γ-has been obtained.

Kinematic Hypothesis
According to the Timoshenko beam theory, cross-sections remain plane after bending but not necessarily normal to the deformed axis of the beam. Hence, the displacement field is identified through the longitudinal and transverse displacements of the beam's axis, u, υ and the rotation of the crosssection φ, which is independent of υ.
The total slope of the deflected beam's axis is due to both bending and shear, and is given by: where z denotes the abscissa along the beam axis and γ is the transverse shear strain. Overall, the generalized strains are given by:

(2)
where ε is the extensional strain and κ the curvature of the beam axis.
Lastly, it is worth noting that, in view of the hypothesis of plane sections, each longitudinal fiber undergoes the axial strain: where y is the distance between the fiber and the geometrical centre of the cross-section.

Constitutive Models and Definition of the Strength Domain σ z -τ zy
With reference to the coordinate system shown in Fig. (1), let's assume that the stress state at any point of any beam's cross section is completely described by the components σ z (axial stress) and τ zy (tangential stress). The formulation of the model stems from the idea of defining a strength domain σ z -τ zy while keeping in mind a plane state of stress, where σ y is supposed to be nil (σ y = 0). In particular, the beam's strength domain σ z -τ zy has been defined by requiring some conditions on the stress to be satisfied. As widely used in the framework of no-tension beam models [27], the material does not withstand traction in the zdirection, i.e. σ z is assumed to be non-positive. Moreover, with reference to a plane state of stress, the new condition assumed herein -referred to in the following as Condition (1)-requires that, at any point of the panel, the shear stress component τ n on any plane -defined by its normal direction n ( Fig. 1)-must respect the Mohr-Coulomb criterion, that is: where σ n denotes the normal stress component parallel to n, and τ o > 0 and m < 0 are two constant material parameters, usually known as cohesion and internal friction, respectively.
With reference to the plane (σ, τ), the foregoing condition corresponds to requiring that any state of stress in the panel, which can be represented graphically through the Mohr circle, will be admissible if it is internal to the cone with vertex at point (-τ o /m,0) and identified by the half-lines τ = ±(τ o +mσ). Hence, by requiring that the circle be tangent to these two halflines ( Fig. 2), it is a simple matter to verify that Eq. (4) can also be written as where σ and σ + denote the principal stresses.
The second condition i.e. Condition (2) consists of assuming that the minimum negative normal stress (compression) on any plane is greater than the minimum admissible compressive stress τ, namely Referring to the reference system (o,y,z) depicted in Fig. (1), where the z-direction corresponds to the beam's longitudinal axis, let's assume that the plane stress state in the panel is

(7)
where σ y is assumed to be negligible (quite an acceptable hypothesis especially for modest stress states [3]).
Consequently, the principal stresses turn out to be (8) and the domains D 1 and D 2 that respectively guarantee fulfillment of Conditions (1) and (2) can be defined in the plane (σ z , τ zy ).
Specifically, by substituting Eq. (8)  with center at point (2mτ zy , 0) and half minor and major axis respectively equal to τ o and (Fig. 3). The ellipse intersects the reference σ z -axis at the points (10) and the τ zy -axis at the points

(11)
As per the condition given on σ z (i.e. σ z ≤0), the domain of interest is limited to the portion in the half-plane of the negative σ z (hatched area in Fig. 3).

(12)
with axis of symmetry coincident with the σ z -axis, vertex at point (σ o ,0) and intersecting the τ zy -axis at points (0,±σ o ). Also in this case the domain is limited to the portion belonging to the half-plane of negative σ z (Fig. 4).
has been denoted with ; in this case, fulfillment of Condition (1) also ensures Condition (2). Conversely, if = , then turns out to be the one depicted in Fig. (5b). However, the first situation is of greater interest in practical applications, as it is more likely to occur for the usual values of parameters σ o , τ o , and m. Hence, hereafter we always refer to such case, i.e.Ɗ ≡ Ɗ 1 .
It is worth emphasizing that the strength domain Ɗ differs substantially from the one obtained by applying the Mohr-Coulomb criterion solely to the shear stress τ zy acting on the panel's horizontal planes, i.e. on the beam's cross-section. In other words, the proposed model differs greatly from that accounting only for horizontal bed joint sliding. Fig. (6) shows a comparison.
According to the strength domain assumed and, as already stated, referring to the case Ɗ ≡ Ɗ 1 for the sake of simplicity, the constitutive equation can be written as follows: (13) with and E the Young's modulus; moreover, let so that Eq. (9) can be written as , and denoting

(15)
with G the shear modulus, the constitutive law for the tangential stress turns out to be

BEHAVIOR OF THE BEAM'S CROSS SECTION
Given the kinematic assumptions set forth in Section 2.1 and the assumed constitutive law, the normal stress σ z has a piecewise linear diagram over any rectangular cross section and can be determined. Various possible σ z patterns can occur, each of which corresponds to a different region in the plane of the generalized strains (ε, κ).
This constitutive equation has been presented in details in [16] for the simpler case in which the effects of shear deformation can be neglected and is not provided herein for the sake of brevity.  Fig. (7). σ z and τ zy patterns over a non-partialized beam cross-section.
As regards the τ zy patterns, some clarification is in order. It is firstly worth noting that, once the values of the parameters τ o and m are fixed, the strength domain given by Eq. (9) is a function of σ z alone; hence, in view of Eqs. (13) and (3), turns out to be a function of two generalized strain only, i.e.
. Furthermore, it is a simple matter to verifiy that, for a fixed pair of values of ε and κ, , depends on the y-coordinate with the same law of dependence on σ z . Thus, the shear stress boundary (y) is still an ellipse over the beam cross-section, with center C and half major axis a equal to: The intersections of the ellipse with the y-axis may or may not fall within the cross-section, depending on the values of ε, κ, in addition to those of the mechanical parameters E, m and τ o , as shown by the following examples.
Lastly, from Eq. (16) and assuming for the sake of simplicity that the shear strain of each fiber is equal to an average value given by the generalized strain, i.e. γ zy = γ, the shear stress τ zy turn out to be a function of ε, κ and γ, i.e. τ zy = τ zy (ε, κ, γ). Overall, many situations must be considered in order to define the patterns of the axial and tangential component σ z and τ zy , depending on the values of ε, κ and γ. Some of those deemed more significant are shown by way of example in the following.
With reference to a rectangular cross-section of dimensions b and 2h, let us suppose that κ > 0 and the section is not partialized i.e. the ε z pattern is the one shown in Fig. (7), where the neutral axis does not cross the section.
In view of Eq. (17), the ordinate of the vertices of ellipse (y) must always be on the opposite sides with respect to y n . Hence, in this case, with y n > h, two qualitatively different circumstances can occurr depending on the values of ε and κ: both vertices of fall outside the section, or just one of them, as shown in Fig. (7) respectively in rows A and B, together with the corresponding σ z pattern. Moreover, by simply comparing γ to of Eq. (15), any possible τ zy pattern can also be obtained, as shown in the Fig. (7).
Essentially, different situations, leading to different τ zy pattern, can occurr. With reference to the situation shown in row A of Fig. (7), it may turn out that γ≤ over the whole cross-section, i.e. for any y ϵ [-h,h], so that τ zy (y)= Gγ; or we can have γ> over the whole cross-section, and in this case τ zy (y) = (y); lastly, γ> may hold solely over a part of the section, with the consequent τ zy pattern shown in Fig.  (7). It is worth just noting that, in case of different positions of the ellipse with respect to the cross-section i.e. for various different values of of the strains (ε,κ)-, the latter situation may provide further possible τ zy pattern, not shown for the sake of brevity.
Similar situations can occur if the section is partialized, as shown in Fig. (8). Note that some situations, such as the entirely elastic case, are not consistent with the position of the ellipse and, thus, can not occurr.
Finally, once the σ z and τ zy patterns are known, for any value of (ε, κ, γ), the generalized stress -N, V and M are determined as

NUMERICAL RESULTS
The proposed model has been implemented into the finite element code MADY [12]. To obtain a locking-free Timoshenko beam element, only one Gauss point has been used for the numerical integration [28]. For the rest, the numerical procedure used is standard and has been described in detail elsewhere [12,15,16,29], with reference to other types of "beam" elements.
To highlight the ability of the proposed beam model to capture the in-plane behaviour of masonry panels, some comparisons have been conducted with the numerical results obtained by means of two different 2D finite element models.
The first bi-dimensional model used has been developed by the authors themselves and is presented in detail in [26]. The model is an extension of the classical masonry-like or notension model [30,31]: besides limiting the compressive strength, the model accounts, for a bound to the shear stress component on each plane, that depends on the acting normal stress component. With reference to shear, the idea of the material constitutive behaviour is, therefore, the same used for developing the proposed beam model. Naturally, the two models differ essentially in the kinematic constraints and the consequent stress state given by Eq. 7. Moreover, they differ with regard to the hypothesis on the cohesive force acting on the cracked parts of the structure, being nil in the beam model but provided for the 2D model. The 2D model has been implemented into the code MADY and has been used in the following analyses with four-node elements in the plane stress framework.
The second model used a 2D model available within the FE commercial code ANSYS. Specifically, 2D geometries under the assumption of plane stress have been discretized by using PLANE182 finite elements, while the Mohr-Coulomb model with cut-off has been selected to represent the masonry's constitutive behavior [25].
It is also worth mentioning that the constitutive law of the two 2D models used are similar but not equivalent [32].

Comparisons for Single Panels
In this section, some comparisons with experimental results for a single panel are provided. Specifically, referring to the experimental research conducted at the Joint Research Centre in Ispra in 1995 [10, 24, 33 -36], two geometries have been considered: A squat panel with 1.35 m height, and a slender panel with 2.0 m height, both of them having a rectangular cross-section 1 m in width and 0.25 m in thickness. As in the experimental tests, the panels are perfectly clamped at their base, with a further restraint to rotation at the top. Moreover, they have first been subjected to a constant vertical load P equal to 150kN, which entails a normal compressive stress σ z equal to -0.6MPa. The self-weight has also been considered, assuming for the mass density a value ρ = 1800 kg/m 3 [24]. Then, all the analyses have been conducted by applying a monotonic increasing horizontal displacement at the top of the panels.
As per to some indications given in the literature [10,24,36], the main mechanical characteristics assumed are as follows: E = 1800MPa, G = 820MPa (i.e.ν = 0.1), σ o = -6.2MPa, and for shear, the values τ o = 0.72MPa and m = -0.4, have been chosen. It is worth noting that such values are reasonable, and fall within the range of values assumed in the aforementioned research papers to calibrate other numerical models. Naturally, the issue of assuming appropriate values for these material parameters is an open one, which requires comparisons with a large amount of experimental data, a task which is beyond the scope of the present paper.
For the numerical simulations, both panels have been discretized into 80 beam elements. The number of finite elements has been chosen to obtain the most accurate response provided by the model, also in view of the comparisons of the stress state given by the more refined 2D models. Decreasing the element number -down to very few elements, has shown to cause slight variations in the results -which are in any case conservative; conversely increasing the number of finite elements has not lead to any significant variation.
The global response of the analysed panels has been presented through pushover curves, i.e. shear force vs top displacement curves. Fig. (9a) shows the results obtained for the squat panel. The analogous diagrams for the slender panel are shown in Fig. (9b).
As shown in Figs. (9a and b), the responses obtained via the beam model are consistent with those from the 2D-models. Unsurprisingly, the response obtained via the beam model exhibits a slightly larger global stiffness with respect to that predicted through the 2D models. Nevertheless, the predicted total lateral strength of the panels are in very good agreement. It is quite interesting to note that no significant variation in the agreement of results emerges between the two structural cases considered, that is by varying the slenderness of the panels.
The accuracy of the beam model's predictions has been also verified from a local perspective. With reference to the squat case, Figs. (10a and b) provide a comparison with results obtained via the 2D models, respectively in terms of minimum principal stress σ min -previously denoted by σ --and maximum shear stress τ max attained throughout the panel. Such comparison is carried out at the loading step where a displacement of 5 mm is attained. As shown in Fig. (10), the beam model yields quite a different result in terms of stress state, despite the similar global response depicted in Fig. (9). The degree of stresses achieved at the top and the base of the panel provided by the beam model is acceptable, if compared with those predicted by the 2D models. However, the predicted stress state throughout the panel is more widespread if compared to that given by the two 2D models, which is confined to a limited area along the diagonal of the panel. Such a distribution cannot be captured by the beam model due to the embedded kinematic contraints.
Like the global responses, the local results for the slender panel are also qualitatively similar to the analogous ones shown in Fig. (10) for the squat case, and hence have not been presented here for the sake of brevity.
In the following, some further comparisons have been carried out by varying the boundary conditions, as they have proven to affect the panels response significantly.
Specifically, for each geometry, the following boundary conditions -widely assumed to represent the actual constraints in real buildings -have also been considered: the panels are perfectly constrained at the base and [i] totally free at the top, a condition that in the following will be referred to as C, i.e. cantilever; the panels, still clamped at their base, have further [ii] restraints on the axial displacement and rotation at the top nodes after application of the vertical load; such cases will be referred to as DF, which stands for double fixed.
Figs. (11a and b) show the pushover curves for the squat and the slender cantilevers (C panels), respectively. The trends of the results are similar to those already shown in Figs. (9a  and b), thus confirming the ability of the simple model to provide an acceptable picture of the global response. Fig. (12) shows the results obtained for the C slender panel in terms of stress state. In this case, contrary to what happens with the panel when top rotation is constrained, the proposed model appears to be more accurate in describing the state of stress in the panel body. Indeed, the stress state prediction is quite satisfactory if compared to that given by the 2D models.
Results for the C squat panel are substantially similar and thus omitted here.

A Wall with an Opening
Simplified models are especially useful for analyzing more complex structures than simple panels, for which sophisticated models may involve very high computational costs. Hence, in order to highlight the possible use of the proposed model in practical applications, some preliminary results are presented for an idealized masonry wall, provided for illustrative purposes only. In this framework, the wall is represented as an idealized frame, in which the piers and spandrels are connected through rigid offsets.
More in detail, the analysed structure is a wall with a centered opening (Fig. 14). The wall's main geometrical characteristics are as follows: height = 3. In addition to the self-weight, a distributed vertical loadwith total resultant of 220kN-is firstly applied at the top of the wall; then, it is subjected to an increasing horizontal displacement, keeping rotation at the top costrained.
A sketch of the idealized wall used in the finite element analysis is depicted in Fig. (14). The dimensions of the rigid nodes have been defined according to the relevant literature [2].
Using the proposed beam model, each pier has been modeled through 40 finite elements, while the spandrel has been discretized into 20 elements. Each node is instead represented by 6 (vertical) plus 6 (horizontal) rigid beam elements.
(a) (b) Fig. (14). The analyzed wall: geometry and beam model idealization. The results have been compared with those obtained with the 2D models used in the foregoing, i.e. the masonry-like model implemented in MADY and the one available in ANSYS, assuming the same mechanical properties, load and constraint conditions. Fig. (15), which shows the total base shear vs top displacement curves obtained via the three models, highlights a consistent trend with some slight difference in the total lateral capacity.
Obviously, although some promising results have been obtained, a number of issues still remain to be settled in order to enable using the proposed model for the analyses of real masonry structures, especially with regards to defining the behavior and geometries of the spandrels and rigid offsets.

CONCLUSION
The need to preserve an enormous heritage of masonry buildings, which as demonstrated by recent seismic events are extremely vulnerable to horizontal actions, calls for simplified models and low-computation numerical methods to perform preliminary and/or wide-area investigations.
It is to this end that a refined beam model has been proposed. Starting with a no-tension material approach, the model has been enhanced by endowing it with the ability to describe shear behavior and predict the collapse mechanisms typically exhibited by relatively squat panels loaded in their plane.
Some comparisons with the results obtained through two different 2D models (which also account for a limit to shear stress) have been conducted to evaluate the accuracy of the results predicted by the proposed model. Despite its simplicity, the beam model is quite capable of capturing the panels global response. This is true regardless of the slender-ness and constraint conditions. The model can also provide some broad indications on the local stress distribution in the panels. Also from a local perspective, the predictive capacity of the model is practically independent of the slenderness of the walls. Conversely, the accuracy of these results is affected by the panel's constraint conditions, being more acceptable for panels which are free at the top.
Lastly, some preliminary results have also been carried out for an idealized wall with an opening, in order to show the possible application of the proposed model to assess the seismic capacity of a masonry wall and facade through nonlinear static analyses.

CONSENT FOR PUBLICATION
Not applicable.