Limit analysis of masonry walls using discontinuity layout optimization and homogenization

A numerical limit analysis model for masonry walls subject to in‐plane loading is posed as a discontinuity layout optimization (DLO) problem, with the masonry conveniently modeled using a smeared continuum (“macromodeling”) approach and a homogenized yield surface. Unlike finite element limit analysis, DLO is formulated entirely in terms of discontinuities and can produce accurate solutions for problems involving singularities naturally, without the need for mesh refinement. In the homogenized model presented, masonry joints are reduced to interfaces, with sliding governed by an associative friction flow rule and blocks are assumed to be infinitely resistant. The model takes account of the interlock ratio of the masonry blocks, their aspect ratio and the cohesion and coefficient of friction of interfaces in both the vertical and horizontal directions. Results from the proposed model are compared with those from the literature, showing that complex failure mechanisms can be identified and that safe estimates of load carrying capacity can be obtained. Finally, to demonstrate the utility of the proposed modeling approach, it is applied to more complex problems involving interactions with other elements, such as voussoir arches and weak underlying soil layers.

encountered at failure, with tailored meshes, adaptive mesh refinement or other techniques required to overcome this. Discontinuity layout optimization (DLO) provides an alternative to upper bound finite element limit analysis formulations and has the benefit of allowing highly discontinuous modes of response to be modeled without difficulty. In the case of masonry, when used in conjunction with a homogenized yield surface, the DLO modeling approach described herein involves adopting a homogenized representation of the composite masonry material, but with the inherently discontinuous modes of response that occur at failure successfully captured. DLO has previously been used to model a range of in-plane, [12][13][14] out-of-plane, 15 and 3D 16 limit analysis problems; here the focus will be on in-plane problems, though the methodology developed is potentially transferable to other cases.
For the masonry wall problems considered here, joints are for simplicity reduced to interfaces and modeled using a Mohr-Coulomb failure criterion, with a prescribed angle of friction and low or zero cohesion, together with a no tension rocking criterion. Masonry units are assumed to be rigid. The aforementioned assumptions have been made widely by previous workers in the context of rigid block limit analysis, considering either associative 17 or non-associative 18,19 friction models at interfaces. Of the non-associative solutions previously reported in the literature (e.g., , it has generally been found that as the number of masonry units increases, the difference between associative and non-associative solutions reduces. When using homogenization theory the representative volume element (RVE) employed is assumed to be very small in comparison to the dimensions of the wall, such that a given modeled wall can be assumed to comprise an infinite number of infinitely small masonry units. As this can be expected to reduce the influence of the flow rule, and in the interests of computational efficiency, an associative flow rule has therefore been adopted in this article. Various homogenization formulations for masonry have been proposed; 23 here the approach proposed by de Buhan and de Felice 1 is adopted, but extended to enable the interlock ratio, , to be varied, and to allow the Mohr-Coulomb material properties to be changed (specifically, the cohesion, c, and the coefficient of friction, , are allowed to be different in the vertical and horizontal directions). By combining homogenization with DLO the aim of the present work is to provide a general and systematic means of identifying arbitrary in-plane failure mechanisms in masonry structures, overcoming the mesh sensitivity associated with finite element limit analysis-based approaches (e.g., the rigid finite element model proposed by Munro and da Fonseca, 9 previously applied to various out-of-plane masonry wall problems, 24,25 is known to provide nonconservative results 26 ), and the limited search space associated with meta heuristic optimization formulations that involve adjustment of the geometries of predefined macro-block failure mechanism topologies, recently applied to the analysis of various masonry constructions. 27,28 Unlike finite element mesh subdivision, DLO is not sensitive to mesh orientation as it considers an array of overlapping discontinuity lines oriented in different directions. However, as with all other methods, computational time is sensitive to model refinement. The proposed formulation also appears to be more generally applicable than a recently proposed linear programming (LP) formulation that allows a range of failure plane orientations to be considered. 29 The article is structured as follows: in Section 2, an extended homogenized macroscopic strength criterion for masonry is derived; in Section 3, the proposed homogenized DLO formulation for masonry structures is presented; in Section 4, results from the DLO model developed are compared with analytical solutions and with discretized solutions obtained by Orduña and Lourenço, 22 Ferris and Tin-Loi, 21 and Gilbert et al, 20 with the homogenized DLO model then applied to a structure also involving voussoir arches and to a settled facade wall problem previously studied by Pepe et al. 30 and Tiberti et al.. 31 Finally, conclusions are drawn in Section 5.

Homogenization theory for limit analysis problems
The fundamentals of the homogenization method for periodic media have been described by Suquet 32 and de Buhan. 33 For masonry structures, assumed to comprise infinitely strong blocks with joints reduced to interfaces, a popular homogenized model was developed by de Buhan and de Felice, 1 and more recently Cecchi et al. 2 developed this theory for in-plane and out-of-plane behavior based on Reissner homogenized plates.
Here the starting point is the model developed by de Buhan and de Felice, 1 but with additional parameters then included. A unit cell is identified as a parallelogram spanning four blocks within a wall with wall aspect ratio r w = w∕h, as shown in Figure 1. Local coordinate axes e 1 and e 2 are indicated and the blocks in the unit cell are denoted A 1 , A 2 , In the kinematic approach the potential failure mechanisms that can occur within this unit cell parallelogram are clearly of interest. The properties of the mortar are the cohesion, c, and angle of friction, , defined c v and v for vertical joints and c h and h for horizontal joints. Other parameters are the masonry unit (block) aspect ratio, r b = w b ∕h b , and the interlock ratio, (0 < < 1). Note that for convenience the terms "strain" and "displacement" are used herein as shorthand for "strain rate" and "displacement rate", respectively. Homogenization theory is based on oscillating functions that prescribe that the macroscopic stress and strain tensors must be the averages of the corresponding microscopic quantities. The theory is based on the assumption that the size of the unit cell, or RVE, is very small compared to the overall dimensions of the bodies, in these case walls, under consideration.
According to the yield design homogenization method, 34 the macroscopic criterion of a periodic composite material is obtained from solving an auxiliary yield design problem associated with the representative unit cell  of the composite material. Denoting G hom to be the macroscopic strength domain and F (.) to be the associated yield strength function, their definition reads: where f (.) is the yield strength function at any point of . is equal to the volume averaged value of over the unit cell: Taking advantage of the virtual work principle, a dualization of the equilibrium equations may be performed, which allows the macroscopic strength domain to be redefined as follows; for a given macroscopic strain tensor D, G hom may be defined as: where U is a displacement jump across a discontinuity surface S disc , d is a microscopic strain tensor and u is a displacement field, kinematically admissible with the macroscopic strain tensor, of the form: where u per (x) is a periodic fluctuation. As denoted in de Buhan and de Felice, 1 the functional ⟨ ⟩, the average of the support functions, is usually denoted maximum resisting work.

Analytical kinematic relations for unit cell
In periodic homogenization theory adjacent cells must fit together in their common deformed state. This means that at every vertex of the periodic cell the strains must be compatible and the stress vectors must be equal and opposite. In other words, for an admissible strain field the displacement fields must be strain periodic. In the kinematic approach described here the masonry blocks possess infinite strength, that is, are rigid, with mortar joints satisfying the Mohr-Coulomb criteria and obeying the following displacement field u(x 1 , x 2 ): where D is a symmetric macroscopic strain matrix, is a rotation matrix, v(x 1 , x 2 ) is a periodic displacement field, and x contains the Cartesian coordinates, where: The periodic displacement must be equal on any two boundary points related by periodicity. This implies that: where 1 and 2 correspond to the components of the periodicity patch . The relative displacement can be obtained from Equation (6): According to Figure 1, the periodicity patch can be expressed as For the four vertices of the cell the position vectors, A 1 being arbitrary chosen as reference, are given by: where To suppress the two global translations, reference point A 1 has been arbitrarily chosen not to displace. The displacement field is given by: The failure mechanisms for unit values of strains D and rigid rotation Ω are shown in Figure 2.

Upper bound estimate of macroscopic strength
The macroscopic strength criterion depends on the normal n and tangential t components of the displacement jump at each joint, j. The following notation is introduced for the displacement jump vector U: Within the unit cell one can identify five different joints between the blocks, with zero strain rate within each block as shown in Figure 1B with zero strain rate within each block (d = 0). According to Equation (4), the maximum resisting work ⟨ ⟩ is defined by: where l j is the length of the joint. Adopting the Mohr-Coulomb failure criterion (or, strictly speaking, simply the "Coulomb failure criterion," given the presence of a joint of fixed orientation), the failure criterion for each joint was defined by Suquet 32 by the support function (U, n) where (c j , j ) are the Mohr-Coulomb strength parameters of joint j. Different material characteristics for horizontal and vertical joints are denoted by suffixes h and v, respectively, such that the cohesion is, respectively, c h and c v and friction h = tan h and v = tan v . Table 1 provides details of the expressions involved in the calculation of the support function (U, n) and the kinematic admissibility from Equation (12).
The maximum resisting work equation (14) evaluates to:  Table 1. Using the expressions in Equation (15), it can be shown that the upper bound of the macroscopic strength criterion is defined by:

Homogenized yield surfaces for masonry walls
The homogenized yield surfaces for masonry walls, taking into account interlock ratio and block aspect ratio, can be derived from Equation (17). The homogenized masonry yield surfaces for = 0.7 are shown in Figure 3, where Case A is applicable for values The equations of the four planes defining Case A and six planes defining Case B are derived in Appendix A (see Equations (A10)-(A13) and (A10), (A14)-(A17), respectively).

Equivalent homogenized strains in DLO
The stages in the DLO procedure are outlined diagrammatically in Figure 4. With DLO a solid body is discretized using nodes interconnected via potential discontinuities, with compatibility at nodes enforced explicitly in the kinematic formulation. Although many discontinuities will cross over one another, it can be shown that compatibility at crossover points is enforced implicitly. 12 Also, an adaptive solution procedure can potentially be used to reduce the number of potential discontinuities that need to be represented in the problem constraint matrix at any given time, 12 though this was not employed in the present study.
Steps in the DLO analysis of a body subjected to a horizontal (seismic) load: (A) Body discretized with nodes; (B) potential discontinuities interconnecting the nodes; (C) discontinuities identified at failure when using the proposed homogenized failure model

F I G U R E 5 Displacement jump and rotation at a discontinuity
Generally in DLO the discontinuities may be straight lines or more complex shaped lines. With reference to Figure 5, the angle between the discontinuity line and e 1 is then denoted by , where and are the x-and y-axis direction cosines, respectively, for the discontinuity line. The components of the tangent and the normal to the discontinuity line are expressed as: The displacement discontinuity may also be expressed in terms of the components along the discontinuity line, denoted by U t , and the normal to the discontinuity line denoted by U n , such that: The flow rule must be respected at the two extremal points of each discontinuity line. Hence at one extremal point U t = s and U n = −n + l∕2, and at the other extremal point U t = s and U n = −n − l∕2, where s is the tangential relative strain, n is the normal relative strain and is the angular rotation of the discontinuity. These terms are collected in the relative displacement matrix d where n is positive for compressive normal strains. These can be expressed as The equivalent strain rate tensor D along a discontinuity line may be expressed as in Hassen et al.: 35 Thus the relationships between the components of D and the components of d depend on the orientation of the discontinuity . The relationship at the two extreme points may be written as: From Equations (19)- (21), we can derive the transformation matrix C: Hence,

Optimization problem incorporating macroscopic strength criterion
A range of strength criteria can be handled using DLO; for example, see Smith and Gilbert 12 for details of how the standard Mohr-Coulomb criterion can be implemented for simple plane plasticity problems. However, for the in-plane masonry problems of interest here, from the macroscopic strength criterion (17) each inequality may be rewritten in terms of six extra positive parameters, or "plastic multipliers", at each end of the discontinuity: These equations can be combined to eliminate Ω. The components of macroscopic strain rate D are then expressed as functions of six plastic multipliers, related by two extra equations: A combination of the plastic multipliers is used to express D 11 and D 22 such that the friction coefficients h and v disappear in the formulation of the support function. The upper bound of the macroscopic support function may then be rewritten in terms of six plastic multipliers: Since the flow rule must be respected at the two extremal points of each discontinuity line, six plastic multipliers are introduced at each end, denoted by p 1 and p 2 . Combining Equation (22) and the three first conditions of Equation (26), the kinematic form of the in-plane problem may be expressed such that the unknowns d are related to the plastic multipliers via the following expression: where The homogenized masonry flow rule in Equation (26) also involves a relationship between the plastic multipliers. This is taken into account by the following expression: together with the g matrix of work coefficients from Equation (27): The optimization problem can then be expressed as follows: where f T L are the live loads, f T D are the dead loads, p is a plastic multiplier matrix, B is a compatibility matrix, and is a load or "adequacy" factor. Finally, p and g are vectors of plastic multipliers and their corresponding work equation coefficients, Q and M are matrices of plastic multiplier terms and C is a suitable transformation matrix.
Alternatively the optimization problem can be posed in the form of the original DLO formulation 12,13,15 (see Appendix B for details): where N is a flow matrix. In this latter formulation it is more convenient to add other flow rules and potentially to implement the adaptive solution process proposed by Gilbert et al.. 15 This formulation has therefore been incorporated in a research version of the LimitState:GEO 36 DLO analysis software, now widely used in industry. This software has been used to obtain solutions to the numerical example problems described in the next section.

De Buhan and de Felice wall
The first example comprises a simple rectangular wall with aspect ratio r w , comprising blocks with two different aspect ratios (r b = 2.286 and r b = 4.571). Adequacy factors have been computed using the proposed homogenized DLO model and compared with both experimental results and with analytical solutions based on a mechanism involving a fracture line radiating at an angle from the toe of wall. The first analytical model considered employs the homogenization theory proposed by de Buhan and de Felice, 1 giving a solution of: and where: The second analytical model considered is one proposed by Heyman, 37 giving the solution of: anal. = min( sliding , rocking ) with Heyman's solutions are based on zero tensile strength elastic theory, with the block aspect ratio r b and material properties c and not considered explicitly. Instead only the wall aspect ratio r w is considered, with the assumption being that not all the mass of the wall will be mobilized in sections of the wall that go into tension, thereby reducing the restoring moment.
Thus analytical, DLO and experimental results are presented in Figure 6, taking the coefficient of friction = 0.6 and with the adequacy factor plotted against the reciprocal of wall aspect ratio (1∕r w ) for the two block aspect ratios considered, that is, r b = 2.286 and r b = 4.571, using a nodal resolution comprising approx. 1000 nodes. These plots are similar in form to those presented by de Buhan and de Felice, 1 with the experimental results again obtained from tests carried out at the Ingegneria Strutturale e Geotecnica Department of the University of Rome on model walls tested to failure by tilting the supporting plane of the structure. 1 To obtain the homogenized DLO results the reported material properties were used, with the joint cohesion c = 0. Figure 7 shows a sample homogenized DLO failure mechanism, together with analytical and rigid block (associative friction) failure mechanisms for comparison (wall aspect ratio r w = 1.0 and block aspect ratio r b = 2.286).
It is evident from Figure 6 that the analytical adequacy factors are higher than those obtained using the DLO homogenization model. This is partly because the analytical equations derived by Heyman, 37 and from homogenization theory of de Buhan and de Felice, 1 assume only a single discontinuity line, whereas the proposed DLO homogenization model allows a number of discontinuity lines to be involved in the failure mechanism. This leads to lower associated adequacy factors , due to a smaller mass of wall being mobilized to provide a restoring moment. Also, since Heyman's analytical equation does not take into account block aspect ratio, it provides poor predictions for the block aspect ratio r b = 2.286 case.

Orduna and Lourenco wall
In the second example, the rectangular wall originally studied by Orduña and Lourenço 22 is considered, taking joint cohesion c = 0 and coefficient of friction = 0.75. First, it is of interest to explore how the adequacy factor obtained using the homogenized DLO model is influenced by the number of nodes used to discretize the wall; Figure 8 shows the convergence characteristics for a wall of aspect ratio r w = 1, containing blocks of aspect ratio r b = 3.0 and with an interlock ratio = 0.5. Sample DLO failure mechanisms are shown in Figure 9, together with discretized block representations of the DLO solutions. The images shown in Figure 9D-F were generated by mapping the centroid of the blocks to the DLO rigid block regions and applying the DLO rigid block rotation and translation to the mapped blocks; interpenetrations between blocks may occur when the blocks lie close to adjacent DLO rigid blocks.
It is also of interest to compare results obtained from the proposed homogenized DLO model with results obtained from associative and non-associative friction rigid block models; results obtained using the methods described by   Gilbert et al. 20 are presented in Table 2 and in Figure 10. Note that Orduña and Lourenço 22 presented in their paper a failure mechanism for a problem comprising 48 courses, with the corresponding non-associative load factor being 0.539 20 (corresponding non-associative and associative load factors of 0.53886 and 0.545034, respectively, were reported by Gilbert et al. 20 ). The results show that the homogenized DLO model produces the lowest adequacy factors while the associative friction rigid block model produces the highest adequacy factors. However, as the number of blocks in the wall is increased, differences in the computed adequacy factors diminish, with the rigid block solutions tending toward the computed converged homogenized DLO model solution ( ∞ = 0.4869). This is to be expected, as the homogenized model is effectively modeling the scenario where the wall includes an infinite number of blocks. Next, as considered by Orduña and Lourenço, 22 the influence of various problem input parameters on the computed adequacy factors is considered. Thus block aspect ratio r b , interlock ratio , and wall aspect ratio r w were varied; see Figure 11, which for comparative purposes also includes analytical homogenized and rigid block (associative friction) solutions, with approx. 800 rigid blocks used in the latter case. As expected, the adequacy factor increases with increasing wall aspect ratio r w , and also increasing block aspect ratio r b . The highest valid adequacy factors are obtained when the factor versus interlock ratio, taking block aspect ratio r b = 3. Wall aspect ratios r w as indicated. anal., analytical homogenized solution; block, rigid block solution; DLO, homogenized DLO solution, Equation (29) interlock ratio = 0.5, since for interlock ratios < 0.5 a reversal in the loading direction would lead to reduced adequacy factors (e.g., for an interlock ratio = 0.3, when the load is reversed the solution will correspond to the = 1 − 0.3 = 0.7 value provided).

Ferris and Tin-Loi walls
Next, the simple rigid block walls previously studied by Ferris and Tin-Loi 21 are considered. Unlike the previous examples considered, some of these walls include openings and were variously modeled on the assumption of associative and non-associative joints in the original study; 21 the walls were later reanalyzed using various discrete analysis models by Gilbert et al. 13,20 Here the wall geometries provided in Gilbert et al. 20 were used, taking the block aspect ratio r b = 3, interlock ratio = 0.5, cohesion c = 0 and coefficient of friction = 0.75 in the present homogenized DLO model. Note that areas immediately above openings needed to be modeled explicitly using rigid blocks in the homogenized DLO model in order to avoid localized failure of the masonry. Note that a rigid block can readily be modeled in DLO by including in the formulation potential discontinuities lying on the perimeter of the block, but no potential discontinuities that cross the region occupied by the block. Previously obtained results are compared to homogenized DLO results in Table 3. For wall problems comprising a small number of blocks and contacts it is evident that the computed homogenized DLO solutions are quite conservative (with the adequacy factor being up to 32.5% lower than the corresponding rigid block (non-associative) adequacy factor); however, the degree of conservatism can be observed to diminish for walls containing a higher number of blocks and contacts (disparity down to 12.2% in the case of the largest problem). A sample solution is presented in Figure 12

Masonry buttress wall
The stability of an old masonry buttress wall under earthquake loading is now considered, with results from rigid block and homogenized DLO models again compared. The wall comprises masonry arches atop buttresses, loosely based on the facade of the Odeon of Herodes Atticus in Athens, Greece, 38 as shown in Figure 13. As in the previous example, in the homogenized DLO model the masonry, in this case arches, immediately above the openings were modeled using rigid blocks. Dry mortar with a cohesion c = 0 and a coefficient of friction = 0.6 was assumed in the analyses.  In the rigid block (associative friction) model the computed adequacy factor = 0.2449 whereas in the homogenized DLO model the adequacy factor = 0.2321, a difference of less than 6 percent. It is also evident from Figure 13 that the failure mechanisms are similar.

Building facade subject to support movement
The next example involves a three-storey wall forming the facade of a model-scale building with six openings, as previously studied by Pepe et al. 30 and Tiberti et al.. 31 The facade is 1.5 m wide and 1.8 m high, and contains 0.2 m wide and 0.3 m high openings. The rightmost section of the wall is subject to a vertical foundation settlement, analyzed via both rigid block and homogenized DLO models, as shown in Figure 14. The material properties were as employed by Tiberti et al., 31 except that for sake of simplicity the blocks were here assumed to be non-tensile resistant and infinitely resistant in compression (in Tiberti et al. 31 the masonry was assumed to have a small cohesive and tensile strength of 0.01 MPa and a compressive strength of f c = 90 MPa; considering the small scale of the wall considered in comparison to the crushing strength, the latter can be expected to have a negligible influence on the solutions obtained). Also, for the homogenized DLO models large blocks representing lintels were positioned above openings, as indicated in Figure 14, to avoid localized failure, with the interface between these and the surrounding masonry modeled with a Mohr-Coulomb material with cohesion c = 0 and coefficient of friction = 0.4; these properties were also assumed for all interfaces between units in the surrounding masonry. All wall elements were assumed to have a unit weight of 12kN/m 3 . The runs were carried out using a specially Two rigid block models were developed, one with "large" blocks (0.1 m long × 0.05 m high × 0.05 m thick) and another with "small" blocks (0.05 m long × 0.025 m high × 0.05 m thick). Thus for the corresponding homogenized DLO model a block aspect ratio r b = 2 and interlock ratio = 0.5 were used to represent the main masonry elements. The support settlement was imposed by replacing the rightmost fixed support with a weightless rigid block constrained to translate vertically and subjected to an upward vertical dead load force of f Ds = 20 kN/m 2 , and a downward variable load f Ls , where f Ls = −f Ds . The resulting foundation reaction force was then calculated from f s = (1 − )f Ds , as outlined in Portioli and Cascini. 39 Computed foundation reaction forces (expressed as kN per metre wall thickness) and identified crack patterns are presented in Figure 14 and in Table 4 for the three models. Results obtained using various numbers of DLO nodes are also included in the table.
As expected, as the size of the constituent blocks is reduced, initially via the rigid block models, and then via the homogenized DLO model, which at the limit can be considered to represent an assemblage of infinitely small blocks, a progressively larger support reaction force was found to be required to support the rightmost section of the wall. (In this case, the 1000 node and 250 node DLO solutions approximately correspond to the small and large unit rigid block model solutions, respectively, with by inference the 2000 node DLO solution corresponding to a rigid block model comprising very small units.) In Table 5, rigid block and homogenized DLO solutions are compared with solutions presented by Tiberti et al. 31 The objective function (OF) defined in Tiberti et al. 31 is the work done by the limiting foundation reaction force for a unit wall thickness, assuming a vertical displacement of 1 mm. In the case of the homogeneous genetic algorithm (GA) models described by Tiberti et al., 31 the OF value increases as the GA population size is increased. In the case of these models the identified failure mechanisms necessarily comprise only a small number of large rigid blocks, a function of the initial coarse triangulated mesh used as a starting point in the proposed adaptive GA-based solution strategy; this is required to keep computational costs manageable. In comparison, the failure mechanisms obtained from the homogenized DLO models include regions of closely spaced fracture lines, which can be identified due to the far wider search space available. This also means that larger OF values are obtained using DLO, at much lower computational cost than when using the GA-based scheme (e.g., using 250 DLO nodes an OF value of 0.00286 was obtained in 1.1 s; a comparable OF value, of 0.00267, obtained using a GA population of 200, required 4736 s to solve, i.e., >4000 times longer, with differences in the CPUs used [Intel i7-3635QM here vs Intel i7-5500U in the previous study] likely accounting for only a very small part of this difference). In future there is also scope to use the adaptive solution procedure employed by Smith and Gilbert 12 to obtain DLO solutions significantly more quickly.

4.6
Building facade subject to soil induced settlement In the final example a soil structure interaction example is considered. The building facade considered in the previous section is now modeled with more realistic dimensions, by scaling up the building dimensions by a factor of six, and also selecting more realistic material properties. A spread foundation is also introduced, composed of the same type of masonry From Figure 15, it is evident that the predicted overall modes of response of the three building facades are similar, though somewhat different to those shown in Figure 14, where only vertical movement of one part of the support was allowed. The example clearly demonstrates the capability of the homogenized DLO model developed herein to be applied to more complex scenarios, in this case involving soil-structure interaction.

Commentary
In many of the examples, it has been demonstrated that, as the number of blocks in a wall increases, results tend toward those obtainable from the proposed homogenized DLO model. This was found to be the case even when a non-associative friction model was employed in the comparative rigid block models-notwithstanding that all the homogenized DLO models considered herein employed an associative friction model. Thus the use of homogenization theory in conjunction with the kinematic limit analysis approach described permits engineers to rapidly analyze in-plane masonry wall problems without need to discretize the wall into a finite number of blocks and mortar joints. As expected, the adequacy factor was found to decrease with increasing wall aspect ratio r w and decreasing block aspect ratio r b . Also, the highest adequacy factors were obtained with an interlock ratio = 0.5 (i.e., regular running bond). Although for sake of simplicity in the present study, constituent blocks have been assumed to be infinitely resistant, a finite masonry strength criterion can potentially be added to the proposed homogenized DLO model, permitting a more diverse range of failure mechanisms to be identified, including, for example, shear cracking in blocks. 40 The approach can potentially be further extended to accommodate material crushing and torsion interaction 22,[41][42][43] and to model out-of-plane problems such as those considered in the context of homogenization by various workers. [44][45][46][47][48][49]

CONCLUSIONS
A new numerical limit analysis model suitable for modeling masonry walls subject to in-plane loading has been presented. This uses DLO in conjunction with a homogenized yield surface to provide a smeared continuum ("macromodeling") representation of the constituent masonry. The homogenized model presented takes account of the interlock ratio of the masonry blocks with, for sake of simplicity, masonry joints reduced to interfaces and blocks assumed to be infinitely resistant.
Since DLO is formulated entirely in terms of discontinuities it can produce accurate solutions for problems involving singularities. However, unlike rigid block representations of masonry, it is not necessary to specify the locations of potential failure surfaces in advance. Thus when DLO is used in conjunction with a homogenized yield surface, the size of the constituent blocks effectively reduces to zero. Practically speaking, this obviates the need to include a large number of blocks in a numerical model in order to safely predict the response a real-world wall. This has been demonstrated in the example problems considered in the present contribution, which include walls with and without openings, subject to horizontal loadings and support settlements. This was found to be the case even when a non-associative friction model was employed in comparative rigid block models, notwithstanding that an associative friction model was used in the homogenized DLO model described herein.

DATA AVAILABILITY STATEMENT
The research version of the LimitState:GEO software used to solve the problems described herein can be made available for research use by contacting the corresponding author. The associated data input files can also be supplied upon reasonable request.

R E F E R E N C E S APPENDIX A. HOMOGENIZED MASONRY YIELD SURFACES
The yield surfaces for masonry, taking into account interlock ratio and block aspect ratio, can be determined from Equation (17). The homogenized masonry yield surfaces are shown in Figure 3 for = 0.7. The set of macroscopic stress fields is a polyhedral cone and is formed by all the convex combinations of the normal to the strain planes. The homogenized strain surfaces can be derived from the following inequalities: From: From: From: The intersection between planes yields two different polyhedral surfaces having four or six planes depending on the values of r b , h , v and . The following two intersections vectors are the same for both planes. Vector 1 is obtained from (A2) and (A3): Vector 2 is obtained from (A4) and (A9): For the four plane yield surface we have the two additional equations as follows. Vector 3 is obtained from (A4) and (A8) and for h − r b (1 − ) < 0: Vector 4 is obtained from (A6) and (A9) and for h − r b < 0: For the six plane yield surface, we have four additional equations as follows. Vector 3 is obtained from (A7) and (A9) and for h − r b (1 − ) > 0: Vector 4 is obtained from (A4) and (A5) and for h − r b > 0 > 0: Vector 5 is obtained from (A3) and (A7) and for h − r b (1 − ) > 0: Vector 6 is obtained from (A2) and (A5) and for h − r b > 0:

APPENDIX B. ALTERNATIVE DLO FORMULATION
In this alternative formulation equation (28) where When = 0 or = 0 a division by zero would occur. Hence special solutions are required as follows.
For the case when = 0, Equation (28) is reformulated by performing one Gaussian elimination step: