Micromodels for the in-plane failure analysis of masonry walls: Limit Analysis, FEM and FEM/DEM approaches

In the last decades the modeling of masonry structures has become an argument particularly appealing for many researchers and a large variety of numerical techniques have been formulated with the aim to produce practical applications in civil engineering, with special reference to the preservation and restoration of cultural heritage. Nevertheless, the question appears today still far from being resolved in a general way. The characteristics of fragility, heterogeneity and anisotropy of masonry, as well as the extreme variety of the building/construction rules strongly compromise the possibility of a unified description of its mechanical behavior. In this work a comparison of different models and techniques for the assessment of the mechanical behavior of two-dimensional block masonry walls subjected to the static action of in-plane loads is presented. Different approaches and numerical models are considered: a Limit Analysis approach (LA), a FEM/DEM procedure and a non-linear heterogeneous Finite Element analysis (FE). Here a standard Limit Analysis is adopted via a homemade procedure based on Linear Mathematical Programming, considering friction at interfaces. Analyses are performed referring to benchmark examples from literature.


INTRODUCTION
asonry is a composite (heterogeneous) material obtained by assembling natural or artificial blocks by means of mortar layers or dry joints and it is one of the more common structural materials adopted for centuries for ordinary or monumental constructions. The investigation of its mechanical behavior plays a fundamental role in view of the protection and conservation of architectures of historical and archaeological interest. However, to deal with the structural response of historical masonry structures is a complex task. In the last decades a large variety of numerical models and approaches have been proposed in literature, but no one can be applied in a general manner regardless the constructive typology. The selection of the most appropriate modelling strategy is indeed strictly related with the nature of the object to analyze. Depending on the adopted model for analysis it is possible introduce three distinct categories: micro-mechanical, macromechanical models and multiscale models. The choice of a micro-mechanical model involves a distinct representation of masonry constituents (units, mortar and unit/mortar interface) whose properties are obtained from experimental tests on small masonry specimens. A micro-mechanical model is suitable for a very detailed response [1][2][3][4][5], but this approach has a limitation represented by the great computational effort due to the high number of degrees of freedom connected to each unit and joint in case of real masonry structures, characterized by considerable number of units. Macro-mechanical models describe masonry as homogeneous continuum, use phenomenological constitutive laws for constituents, including also some inner variables for damage, and friction coefficients. The material parameters are derived by means of experimental tests on small masonry specimens or directly on the single constituents. Macro-mechanical models are characterized by high computational efficiency, as they do not provide an accurate description of the internal structure of masonry material [6][7][8].
Multiscale continuum models represent a very promising approach for the analysis of masonry structures since they can accurately retain memory of the mechanical and geometrical properties of the material (micro-structure) together with the capability to contain the computational effort compared to a fully micro-mechanical model [9][10][11][12]. These models are often derived by considering only two material scales: a micro-scale where, after deducing the mechanical properties of the components, preferably through experimental tests, a material representative volume element (RVE) is defined and a macro-scale continuum model, is obtained by performing a homogenization procedure in most cases based on the solution of notable boundary conditions problems for the RVE. Among these approaches, both concurrent and semiconcurrent multiscale models have been proposed in the literature for masonry-like materials, the first ones referring to a strong coupling between the micro-and macro-levels [13][14][15][16], the latter ones, often referred to as computational homogenization models, characterized by an only weak (although two-way) coupling between them [17][18][19][20]. Most of the existing approaches are devoted to the mechanical behavior of periodic (i.e. regular) masonries, for which a suitably defined unit cell plays the role of RVE, but there exist also different homogenization techniques for both linear and nonlinear analyses of random microstructures, already applied or directly applicable to irregular masonry structures [21][22][23]. Other multiscale strategies have been proposed that exploit different homogenization techniques based on the so-called Cauchy rule, and its, generalizations [24] that allowed the derivation of both classical and generalized continua able to properly represent scale effects, that in masonry materials are prominent [25][26][27][28][29]. In this work the attention is mainly focused on the category of micromodels particularly focusing on Limit Analysis, which represents a very effective tool to estimate the collapse load and collapse mechanism for one-leaf masonry structures [30][31][32][33][34][35]. In particular the model here presented considers an associative flow-rule connected to a dilatant behaviour of the joints. This hypothesis has been successfully adopted to study the behaviour of historical masonry structures [36][37][38] even if it may occur that the flow rule is quasi-associated under shear actions, depending by the level of pre-compression of masonry. It must be pointed out that a limitation of Limit Analysis (both associated and not associated) concerns ductility that is assumed infinite, and this hypothesis not always fits well the softening behaviour of masonry [39]. A validation of the proposed model is provided, via suitable comparisons with two models that finely describe the microstructure and already adopted to the evaluation of the in-plane failure behavior of masonry panels: (i) a discrete model based on a combined Finite/Discrete Element method (FEM/DEM) [40,41], and (ii) a heterogeneous nonlinear Finite Element model (FEM), derived from the works [13,15,16]. The three models are used to reproduce the analysis performed in [42] -here regarded as a benchmark -on several masonry panels with different height-to-width ratio and with or without openings. The comparison of the numerical M results shows the efficacy of Limit Analysis for a fast and reliable assessment of the in-plane failure analysis of masonry walls.

LIMIT ANALYSIS
ollowing and reinterpreting the model formulated in [30,31], masonry is described as a system of n rigid blocks and m joints unable to carry tension and resistant to sliding by friction. The blocks can translate and rotate about the edges of the contact blocks (hinging) as well as slide along the joints.
the orthonormal basis in the three-dimensional space. We consider the two blocks in Fig. 1.
Loads are applied to the centroid of each rigid block th i : static dead loads are collected in vector , , . For the whole structure it results The vector of the load over the whole system is . The vector   j  ε ε refers to the whole structure and corresponds in a virtual work sense to the vector of static variables σ .  where Eqn. (1.1) 1 represents the kinematic compatibility for the whole system of interfaces and blocks in which B is the compatibility matrix as in [42], but also in nonlinear static analysis by [43], Eqn. (1.1) 2 defines the equilibrium for the whole structure, Eqn. (1.1)3 describes the generalized yield domain of the system where N is the block-diagonal gradient matrix, Eqn. (1.1) 4 represents the flow rule which express the vector ε as a linear combination with non-negative coefficients λ, called inelastic multiplier and M is the block-diagonal matrix of the modes of failures. Eqn. (1.1) 5 is the complementarity condition which defines the plastic behavior of contact surface. Moreover, the collapse mechanism must be characterized by a non-negative work of the live loads, defined by Eqn. (1.1) 6 . Within the framework of the holonomic perfect plasticity, the same relations govern the problem of a non-standard rigid-plastic discrete materials. Resorting this formal analogy, the collapse load for a masonry structure, under the hypothesis of proportional load with the factor 0   , can be determined.
After some algebra the authors obtained the following non-linear and non-convex programming problem (NLNCP) with the unknowns  , 2 σ , λ and the bounds 0  λ , where 2 σ are the undetermined unknown of the system which represent the statically undetermined term of the generalised contact stress [30] In order to deal with the NLNCP, authors [31] developed a specific code, called ALMA (Analisi Limite Murature Attritive), which used a two-step procedure to solve the problem: in the first step a linear programming problem (LP), obtained by replacing friction with dilatancy, is solved; in the second step, the NLNCP solution is approached using, as initial guess for the unknowns of the problem, the solution of the first step. However, the problem of Limit Analysis of structures with frictional interfaces (non-standard LA) is numerically very difficult to be solved. The solution could not exist and when it exists it could be a local minimum instead of the global one [44]. On the other hand, due to the presence of non-associative flow rules, the Drucker stability postulate no longer holds and the solution in terms of contact actions and collapse load factor loses its uniqueness. Moreover, bi-dimensional or three-dimensional real structures, characterized by many degrees of freedom, increase the computational complexity of the problem. From Limit Analysis theory it is well known that if normality rule holds, i.e. the vector of inelastic strain results normal to the yield surface, the static and kinematic theorems of limit analysis could be formulated in a linear programming context, resulting in two dual problems, which lead to a unique solution. In particular, following the approach of [31], results of this work refer to a linear programming optimization problem related to the kinematic approach of Limit Analysis, which provides the collapse multiplier and the corresponding mechanism. Friction is considered in term of dilatancy. The kinematic problem is defined as with the bounds on the unknowns 0  λ . To overcome some computational limits of the original code ALMA, mainly related to the number of blocks and interfaces involved into analysis, a new version of the code, ALMA 2.0, was implemented using MATLAB for linear optimization and a Python interface for pre-and post-processing operations.

FEM/DEM ANALYSIS
he limit analysis model is compared with the models adopted in [40,41], here regarded as a benchmark. In the referred works, two micro-mechanical models have been proposed for the in-plane failure analysis of masonry walls: a discrete element method (DEM) and a combined finite/discrete element method (FEM/DEM). Both λ A N f models considered fall in the field of discrete or distinct element methods, which have been proved to be particularly suitable for the study of masonry structures [45]. DEM model is based on the original numerical method formulated by [46], and recently developed by [43]. The model is based on the assumption of rigid block and mortar joints modeled as zero-thickness elastic-plastic interfaces, adopting a Mohr-Coulomb yield criterion. Masonry is seen as a system of rigid blocks, whose interactions are represented by forces and moments depending on their relative displacements and rotations. FEM/DEM method is a combination of discrete elements, originally formulated by [47], and developed by [48], it consists in a discrete element method in which the individual elements are meshed into finite elements, adopting a triangular discretization of the domain with embedded crack elements that activate whenever the peak strength is reached. The method, initially developed in the field of geo-mechanics, has been adopted to study the behavior of historical masonry [49]. Differently from the DEM described above, blocks can be assumed to behave as rigid or elastic bodies. Mortar joints might be idealized as elastic or elastic-plastic zero-thickness Mohr-Coulomb interfaces. Blocks are modeled by means of finite elements while interfaces are modeled as discrete elements. In this work, discrete models are realized adopting the FEM/DEM method.

HETEROGENEOUS FEM ANALYSIS
he second comparison model is a heterogeneous finite element model, by which masonry is described as a system of n deformable blocks and m dry joints in which all the nonlinearities due to unilateral contact together with decohesion and friction phenomena are lumped. In detail, masonry blocks are made by bulk finite elements with linearly elastic and isotropic material properties, whereas dry joints are represented as zero-thickness damage-based cohesive interface elements placed in between, equipped with an intrinsic mixed-mode displacement-based traction-separation law (TSL). The variational formulation for the equilibrium problem at the microscopic scale is rather classic and, thus, not reported here for the sake of brevity. However, with the aim to find more theoretical and computational details, the reader is referred to the works [13,16,17], in which this formulation is embedded within a more general multiscale framework. In this paper the following mixed-mode TSL of the type   = T T  , proposed by [50], is chosen (see Fig. 2), incorporating unilateral contact with friction in an approximated manner, exploiting a penalty approach to enforce the highly nonlinear kinematical constraint in the normal direction to the joint surface: over the entire deformation history; finally,  is the friction function, defined as:  being the friction coefficient. It is worth noting that the adopted friction function is different than that appearing in [46], due to the presence of an exponent (in this work set as 0.5, in general less than 1) applied to the dimensionless mode-II separation. This novel formulation is able to account for the rapid increase of friction forces between the two sides of masonry joints during its decohesion, coherently with what suggested in [51]. The adopted TSL is versatile, meaning that it can be applied in this form for any masonry type (with dry or mortar joints), characterized by a hybrid cohesive/frictional mechanical behavior. However, for masonries characterized by a dominant frictional behavior, as in the case of dry joints analyzed in this paper, the cohesive part of the TSL (1.4) tends to vanish, so that the adoption of the friction function (1.7) allows friction forces to develop up to their maximum value from the beginning of numerical simulations, i.e. at low separations between adjacent units. Results are similar to those of the smaller panel previously analyzed, with a hinging behavior which exhibits a rotation of the upper right portion of the panel and three 'stair-stepped' cracks of the wall now obtained with all the three models. Also, in this case the FEM-DEM present a sliding of the blocks not observed into result of the other models but, by the adoption of a fictitious cohesion value for head joints, could allow to obtain a mechanism less influenced by sliding. Fig. 6 refers to results obtained for Example 3, characterized by the presence of an opening. In this case the results obtained with the three models are exactly the same, with a hinging behavior which implies the rotation of half panel around a hinge positioned about at the lower right corner of the panel. All the models compute a principal diagonal crack passing through the opening. Some slight differences could be pointed out for the mechanism corresponding to the FEM, as for example the position of the hinge which is located at ground level (for Limit Analysis and FEM/DEM it is positioned upon the first row of blocks) and the absence of movement of the portion of the panel to the left of the opening (the other two models detect a slight rotation of this macro-block around a hinge positioned exactly at the left lower corner of the window).
(a) (b) (c)  Fig. 9 refers to results obtained for Example 6, characterized by the presence of four openings. The collapse of the panel is due to a hinging behavior with the formation of several diagonal cracks which divide the structure into various macroblocks. About the position of the cracks, the models provide results slightly different, except for those that identify the macro-blocks positioned to the right of the panel. Another difference is observed for FEM that compute the left part of the panel as stable, while Limit Analysis and FEM/DEM reveal in that portion different cracks and rotating macro-blocks. The collapse multiplier obtained with Limit Analysis, FEM/DEM and FEM has been compared with the results related to the associative (LP) and non-associative (Mixed Complementarity Problem MCP and Mathematical Program with Equilibrium Constraints MPEC) problems solved by Ferris and Tin-Loi [42]. As expected, results obtained with ALMA2.0 are very close to the one obtained with the associative case.
Referring to results of Example 1, Fig. 10a, it could be noticed how ALMA 2.0 as well as FEM provide results in good agreement with benchmark ones, with negligible numerical differences. FEM/DEM instead returns a lower value of the collapse multiplier. The same is for Example 2, Fig. 10b, where it could be observed as the Limit Analysis model reaches the same value of LP problem resulting little higher respect to those of the non-linear problem. FEM on the contrary provide a lower value of , closer to MCP and MPEC problems. Also, in this case FEM/DEM returns a lower value of collapse multiplier. The response of ALMA 2.0 for Example 3, Fig. 10c, is always closer to the result of the LP problem but in this case, it is little higher. FEM confirms its capacity to get results in good agreement with those of MCP and MPEC problems, while FEM/DEM provides a slightly higher value of collapse multiplier.
The difference between models' responses is more evident referring to Figs. 10d, 10e and 10f, where results of Examples 4, 5 and 6 confirms, as expected, that the value of collapse multiplier provided by the solution of the linearized problem treated with ALMA 2.0 is very close to those corresponding to the associative problem. On the other hand, the capacity of FEM to get collapse multiplier values closer to the non-linear models becomes clearer as well as the results provided by FEM/DEM, even if with some negligible numerical difference.

CONCLUSIONS
he present work focuses the attention on the comparison between different modelling techniques. In particular, Limit Analysis (LA), FEM/DEM and FEM strategies have been analyzed and applied to the benchmark masonry panels introduced by Ferris and Tin-Loi [42]. The three models object of this study are able to reproduce correctly the collapse mechanism, taking into account the real texture of masonry walls, describing accurately each block's geometry and disposition, obtaining useful information about the collapse mechanism and the potential crack patterns that may develop. The comparison of the results obtained by means the different discrete models points out some advantages of using LA approaches to model the structural response of masonry panels with respect to the other techniques presented in this work. The LA approach requires, indeed, less computational effort with respect to FEM/DEM and FEM techniques, which could present a critical issue especially for structure with a large number of degrees of freedom. Moreover, another advantage of the LA approach concerns the limited number of mechanical parameters to be introduced as inputs of the numerical model-. Indeed, unlike FEM/DEM and FEM, which require more mechanical information, the only parameter needed using LA is the friction angle. Anyway, by the comparison of results in terms of collapse multiplier, it has been noticed that FEM, unlike the other models, is able to get values closer to those related to the non-associative response obtained by Ferris and Tin-Loi [42]. The next step of this research will be focused on the analysis of more complex geometries, including a linearized procedure to take into account pure shear according to Coulomb friction as well as crushing of the blocks. The assumption of infinite compressive strength, typical of unreinforced masonry structures, could be removed improving the code with the possibility of a crushing failure of masonry. This aspect is particularly interesting for structures reinforced with metallic ties where a concentration of stress could arise. According to several literature contributes [52] a procedure that iteratively modifies the yield surface, by adding a series of new constraints that consider the limited compressive strength of joints, could be introduced. Another further development, that has been recently implemented in the code, is the introduction of contact surfaces with cohesion in the formulation. The study about the influence of settlements on the global structural response is also an ongoing research.