{'fraction_non_alphanumeric': 0.026066298097204645, 'fraction_numerical': 0.010324385532538467, 'mean_word_length': 4.698380566801619, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Tracking multi-directional intersecting cracks in numerical modelling of masonry shear walls under cyclic loading

In-plane cyclic loading of masonry walls induces a complex failure pattern composed of multiple diagonal shear cracks, as well as flexural cracks. The realistic modelling of such induced localized cracking necessitates the use of costly direct numerical simulations with detailed information on both the properties and geometry of masonry components. On the contrary, computationally efficient macro-models using standard smeared-crack approaches often result in a poor representation of fracture in the simulated material, not properly localized and biased by the finite element mesh orientation. This work proposes a possible remedy to these drawbacks of macro-models through the use of a crack-tracking algorithm. The macro-modelling approach results in an affordable computational cost, while the tracking algorithm aids the mesh-bias independent and localized representation of cracking. A novel methodology is presented that allows the simulation of intersecting and multi-directional cracks using tracking algorithms. This development extends the use of localized crack approaches using tracking algorithms to a wider field of applications exhibiting multiple, arbitrary and interacting cracking. The paper presents also a novel formulation including into an orthotropic damage model the description of irreversible deformations under shear loading. The proposed approach is calibrated through the comparison with an experimental test on a masonry shear wall against in-plane cyclic loading.


Introduction
Damage in masonry appears in the form of localized cracks. The propagation of these cracks within the structure determines the development of the collapse mechanism and consequently the capacity of the structure under the given loading conditions. This structural behaviour has been taken into account for the development of both classical and advanced analysis tools for the structural assessment of masonry constructions [1][2][3][4].
Numerical modelling of masonry structures has undergone important developments in the last decades. The inventory of numerical methods today includes simplified techniques such as limit analysis software [5][6][7][8], equivalent frame methods [9][10][11] and macro-modelling approaches [12][13][14] to sophisticated ones as micro-modelling [15,16] and multi-scale finite element models [17][18][19]. Simplified numerical methods, like the macro-modelling ones, consider masonry as a continuum material with average properties. This simplification allows for the numerical analysis of large scale structures at an affordable computational cost with however limited accuracy in the damage representation. Reciprocally, the explicit consideration of the masonry micro-structure in micro-modelling techniques results in a more realistic representation of cracking, at the cost of increased computation time and memory requirements that limits up to date their application to small structural elements rather than whole structures.
An affordable alternative to the simulation of localized cracking in large scale structures was proposed by Cervera and co-workers [20][21][22]. The method bases on the enhancement of the smeared crack approach with the use of a crack-tracking algorithm. The outcome of this approach is the realistic representation of evolving localized cracks within the analysed structure for a computational cost comparable to that of the classical smeared crack approach [23,24].
Despite the above recent applications on masonry structures, tracking algorithms have a long history in the simulation of localized failure in quasi-brittle materials. They are used for the enhancement of problems involving the progression of both weak, or smeared (e.g. [22,25,26]) and strong, or discrete, discontinuities (e.g. [27][28][29][30]). Independently of the chosen approach, their application aims to circumvent the problem of the mesh-dependency that the irreducible standard finite element formulation presents. An additional feature is the realistic representation of cracking, since the localization of damage occurs within a narrow band of the mesh. Notwithstanding the aforementioned advantages, tracking algorithms are commonly designed with important assumptions regarding the initiation, propagation and completion of the simulated cracks. Such assumptions have limited the use of these approaches to the simulation of propagating cracks under monotonic loading. This paper presents a methodology to simulate efficiently the localized cracking in masonry structures induced by in-plane cyclic loading using a continuum damage mechanics model enhanced with a crack-tracking algorithm. Cyclic loading of masonry structures may result in complex damage patterns characterized by multiple and intersecting cracks. To appropriately model this, a novel numerical methodology is proposed, on the basis of the tracking algorithm developed in previous works [22,31,32], that allows the simulation of intersecting and multidirectional cracking, extending the use of tracking algorithms to a large field of applications. At the constitutive level, the paper presents a novel formulation in the context of a continuum damage model for the description of irreversible deformations under shear loading.
After this introduction, the paper is organized in the following way. Section 2 presents the used constitutive damage model. Following that, the novel tracking algorithm is detailed in Sect. 3. The proposed approach is calibrated in Sect. 4 through the simulation of an experimentally tested masonry wall under cyclic shear loading. The same section addresses the issue of mesh objectivity of the tracking algorithm. The paper ends with the concluding Sect. 5.

Constitutive model
The material non-linearity is represented in the numerical simulations through the use of the constitutive model presented in [33]. It is a strain-based continuum damage model with two scalar damage variables for distinguishing between tensile and compressive damage. The decomposition of the effective stress tensor into positive and negative parts, as well as the use of the two damage variables, results in a unilateral effect under tension-compression cyclic loading. The model can account for irreversible deformations through a simple and efficient algorithm, within the context of continuum damage mechanics. The above properties, as well as the limited number of involved material parameters and the simple explicit formulation, make it a suitable model for the seismic analysis of large concrete and masonry structures [34][35][36][37].
For the sake of completeness, this section presents the basic ingredients of the aforementioned tension-compression damage model. The shear behaviour of the model is carefully considered through the proper definition of the damage surfaces and the evolution of the irreversible strains. The modelling of the irreversible strains is achieved by improving previous formulations mainly intended for concrete structures in [33,34,38]. The section ends with a discussion over the effect of the proposed alterations to the evolution of the irreversible deformations under shear loading.

Constitutive law
The Cauchy stress tensor (of rank two) r is computed in the tension-compression damage model through the constitutive relation where d þ and d À are scalar variables denoting the compressive and tensile damage affecting the positive r þ and negative parts r À of the effective stress tensor respectively. The latter are calculated through the distinguishing stress split according to [38] In the above equations, r j is the principal stress corresponding to the eigenvector p j of the effective stress tensor r and the symbols hÁi are the Macaulay The model uses two scalar variables s AE , the equivalent stresses, to distinguish between loading, unloading and reloading. The adopted values of s þ and s À determine the shape of the d þ and d À damage surfaces. In this work we have adopted the values proposed in [19,39] In the above, f þ and f À denote the tensile and compressive uniaxial strengths respectively and f À b the biaxial compressive strength.
I 1 is the first invariant of the effective stress tensor and J 2 the second invariant of the deviatoric effective stress tensor. Finally, r max and r min stand for the maximum and minimum principal effective stresses respectively. The use of the Heaviside step function (H½x ¼ 1 if x ! 0 and H½x ¼ 0 if x\0) activates the s þ for tensile and shear stress states (1st, 2nd and 4th quadrants in Fig. 1), and the s À for compressive and shear shear stress states (2nd, 3rd and 4th quadrants in Fig. 1).
The j 1 variable controls the shape of the failure surface in the shear quadrants (i.e. 2nd and 4th quadrants in Fig. 1). Its value varies between 0 (i.e. the Drucker-Prager criterion) and 1 (i.e. the criterion proposed in [40]). Note that the use of the Macaulay brackets for r max in Eq. (5) makes the function s À independent of j 1 for pure compression. This, in addition to the fact that j 1 does not appear in Eq. (4), make the tensile and compressive response in the 1st and 3rd quadrants of Fig. 1 independent of j 1 . As shown by Petracca et al. [39,41], the different choice in the definition of the equivalent stresses from those proposed in [35,38,40] results in a better control of the model behaviour under shear monotonic loading. Fig. 1 Effect of the s À function to the shape of the d À failure surface in the shear quadrants for the plane stress case As it will be shown later, the choice of such surfaces in the current model results also in the consideration of the irreversible deformations under shear loading.
Following the above, the damage criteria can be expressed as where the internal variables r þ and r À represent the current damage thresholds and the respective expansion of the d þ and d À damage surfaces. Their initial values are defined to be equal to the uniaxial tensile strength r þ 0 ¼ f þ and the uniaxial compressive strength r À 0 ¼ f À . After damage is triggered, their values are equal to the maximum attained values by the equivalent stresses with n representing the number of the current (load / displacement) increment. The effect of the irreversible deformations is accounted using the following definition of the effective stress tensor where C 0 denotes the (rank four) elastic constitutive tensor and i the irreversible part of the second order total strain tensor and e the elastic one. Contrary to the total strain tensor, which is an external variable of the model, the irreversible strain tensor is an internal variable. The following expression is proposed in this work for the rate of the irreversible strains or in incremental form as In the above equations, n and n þ 1 subscripts refer to subsequent time steps. The symbol ''D'' denotes the increment of a specific quantity from step n to n?1. The variables r and r À ð rÞ represent, respectively, the effective stress and the stress threshold in the absence of increment of irreversible strains at time step n þ 1 rð rÞ ¼ max r À n ; s À ð rÞ Note that the quantity r À ð rÞ À r À n is always nonnegative due to Eq. (14), an essential thermodynamic requirement for ensuring non-negative dissipation. Finally, the material parameter b, with values between 0 and 1, defines the rate of the irreversible deformation by scaling the current increment of the equivalent stresses, and consequently the stress thresholds, through the following equation ¼ ð1 À bÞ ðr À ð rÞ À r À n Þ ð 15bÞ It is convenient to consider the effect of the irreversible strains through the definition of the effective stresses, since the Cauchy stresses are computed using the values of the effective stresses r AE and the two damage variables d AE . This can be performed through the update of the effective stresses using Eqs. (10) and (11) in a similar way to [33] r nþ1 ¼k r ð16Þ The computation of the effective stresses at each increment of the analysis is presented in the flowchart of Table 1.
Following the above, the evolution of the damage variables are defined, which in this work follow the exponential relationship proposed in [42].
Tension and compression evolution laws consider respectively the positive G þ f and negative G À f fracture energies, as well as the characteristic finite element width l dis through the corresponding discrete softening parameter H AE d ensuring mesh-size independent energy dissipation according to the crack-band theory [43]. For the case of tension this is [44] For the case of compressive damage, the discrete softening parameter H À d is defined as (see Appendix) The above definition is consistent with the crack bandwidth approach yielding objective results for different values of b by considering the contributions to the dissipated energy of the evolution of the irreversible strains and the compressive damage. In the above, the material characteristic length l AE mat for tension and compression is The use of the discrete softening parameter of Eq. (20) results in an objective representation of the dissipated energy Gf À . This is illustrated in Fig. 2 for the case of a cyclic uniaxial compressive loading using the proposed H À d of Eq. (20) instead of the one adopted in [45], which is the negative counterpart of Eq. (19). The use of the new softening modulus (Fig. 2a) results in a consistent energy dissipation according to the input value of Gf À independently of the choice of b.
The effect of b can now be visualised using Eqs. (11) and (15) for the two limit cases of b ! 0; 1 The first case, with b ¼ 0, recovers the standard continuum damage mechanics model without any evolution of irreversible strains. Contrariwise, for b ! 1 the effective stresses are updated such that there is no evolution of the equivalent stress s À and correspondingly r À . For this case, there is no compressive damage evolution and the total incremental strains at this step will be irreversible. A value of b within 0 and 1, gives a proportional evolution of both irreversible strains and compressive damage. The above effect is illustrated in Fig. 2a.

Irreversible strains under shear loading
As implied by their definition in Eq. (11), irreversible strains i are influenced by the selection of the compressive damage criterion and thus of the (iii) Calculate s À ð rÞ using Eq. (5) and r À ð rÞ using Eq. (14).
Effect of the choice of the discrete softening parameter for compression to the softening behaviour for different values of b: H À d a using Eq. (20), b using the formulation in [45] compressive equivalent stress s À . This means that the evolution of the irreversible strains is possible only when compressive damage is increasing. The effect of this choice is investigated here locally under pure shear. The assumed material parameters for this demonstrative example are presented in Table 2.
For comparison purposes, four cases have been considered with different definitions of s À . In the first one, referred hereafter as Case-1, the s À is calculated using only the negative effective stress tensor ( J 2 À , I 1 À , r À max , r À min ) as The above definition using the negative part of the effective stress tensor, results in an ''open'' d À damage surface in the shear quadrants, as shown in Fig. 1 for the plane-stress case. This means that in this local test, compressive damage under pure shear is triggered only when the minimum principal stress reaches the uniaxial compressive strength of the material. For the other three cases, the equivalent stress of Eq. (5) is used that is a function of the total effective stress tensor ( r) with j 1 =0.0; 0.5; 1.0 for Case-2, Case-3 and Case-4, respectively. The shapes of the four damage surfaces are presented qualitatively in Fig. 1. Each case has been subjected to a shear deformation composed of two stages: (a) loading up to a shear strain of c xy ¼ 0:5 Â 10 À3 and unloading to zero, (b) loading up to a shear strain of c xy ¼ À2:5 Â 10 À3 . Figure 3 presents the response under shear for the four different cases. The letters appearing in the graphs refer to the following conditions occurring during the two loading stages: Point ''A'' triggering of tensile damage, Point ''B'' triggering of compressive damage, Point ''C' maximum shear deformation for the first stage of loading, Point ''D'' Irreversible deformation for s xy ¼ 0, Point ''E'' activation of the tensile damage during the second stage of loading, Point ''F'' activation of the compressive damage during the second stage of loading.
The part O-A in all graphs is the same and it represents the linear behaviour until reaching the d þ damage surface. After this point, the response in the four models depends on the specific assumption of the negative equivalent stress s À , i.e. to the specific shape of the d À failure surface in the the second and fourth quadrants of Fig. 1. For Case-1 and Case-2, the compressive failure criterion U À is not reached for the given maximum shear strain of c xy ¼ 0:5 Â 10 À3 . This results in the increase of the shear capacity during the loading path A-C (see Figs. 3a, b). An additional consequence of this is that there is no evolution of irreversible strains since the d À failure surface is not reached (paths C-D in Figs. 3a, b). On the contrary, evolution of irreversible strains occurs for both Case-3 and Case-4, since there is evolution of the compressive damage d À within the loading path A-C (Figs. 3c, d).
In Case-4, the overlapping of the tensile d þ and compressive d À failure surfaces results in the triggering of the tensile and compressive damage at the same time (B C) and in the earlier evolution of the irreversible strains. For this reason, the unloading path (C-D, Figs. 3c, d) is less sharp than the one in Case-3, resulting in larger irreversible strains.
The damage occurring during the first stage of a cyclic static loading analysis (O-D path) has an important effect during the reverse loading at the second stage. The stiffness during the reverse loading depends on the d þ and d À values attained during the first stage of the cyclic loading analysis. Damage increases once again when the effective stresses s AE overcome again the attained thresholds r AE (Point E for d þ and Point F for d À ). For Case-3 and Case-4, this happens at the same increment for both s AE (E F in Figs. 3c, d). The capacity of Case 2 keeps increasing up to the satisfaction of the failure criterion U À (Point F, Fig. 3b). However, this does not happen for Case-1, resulting in an increasing shear capacity.
The above analyses show that the description of the unilateral effect under shear cyclic loading, and thus the definition of the unloading stiffness, requires the careful calibration of the equivalent stresses s AE as well as the magnitude of the irreversible strains through the b parameter.

A tracking algorithm for intersecting and multidirectional cracks
The simulation of localized cracks within the macromodelling approach is achieved through the use of a local crack-tracking algorithm. The aim of the algorithm is to identify the elements crossed by propagating cracks at each time/load increment. The crack is then simulated through the use of a nonlinear stressstrain relationship for the elements within the crack path, while the ones outside maintain a linear elastic response. Tracking algorithms are used in different numerical approaches developed for the simulation of localized cracks in quasi-brittle materials (see [46]). Despite the diversity of the applied fields, their function is always the same, i.e. the localization of strains in a limited part of the structure, where cracks are expected to propagate. The methodology used here is based on the algorithm presented in [20][21][22]32]. The local nature of this crack-tracking technique enables the identification of the crack trajectories from information already available at each step of a finite element analysis, i.e. the stress state of each finite element. In this way, computational cost, in terms of time and memory resources, is kept to values comparable to the classical smeared crack approach (see [32]). This strategy aims to a Shear strain γ xy × 1000 Fig. 3 Response under shear reverse loading for different assumptions about the compressive damage criterion. a Case -1 (s À ¼ s À ð r À Þ), b case -2 (s À ¼ s À ð rÞ and j 1 ¼ 0:0), (c) case -3 (s À ¼ s À ð rÞ and j 1 ¼ 0:5), d case -4 (s À ¼ s À ð rÞ and j 1 ¼ 1:0). Point A: Triggering of tensile damage, Point B: Triggering of the compressive damage, Point C: Maximum shear deformation for the first stage of the loading, Point D: Irreversible deformation for s xy ¼ 0, Point E: Activation of the tensile damage, Point F: Activation of the compressive damage good balance between accuracy and efficiency, which is of primary importance in problems involving multiple and intersecting cracks, as is the case of cyclic loads in masonry structures.
This section presents the novel developments making possible the simulation of intersecting and multi-directional cracks using a local crack-tracking algorithm. The new algorithm is a further development of the one presented in [32], allowing the modelling of arbitrary cracking, starting from the boundary or the interior of the analysed structure and propagating towards one or two directions. The latter enhancement is fundamental to simulate the response of structures experiencing multiple and diverse cracking due to flexural, tensile or shear stress states.
The following of this section focuses on the details of the novel algorithm. Firstly, a general introduction on the tasks of the tracking algorithm is presented. This is followed by the up to date assumptions regarding the modelling of crack initiation under monotonic loading. Modelling of crack propagation and intersection is presented in Sect. 3.3. Finally, the section closes with the novel approach to simulate multi-directional cracking.

Structure of the tracking algorithm
The crack-tracking algorithm presented here is called at the beginning of each numerical increment before the calculation of the stress-state and has two tasks: (a) to identify the origin of new cracks and (b) to recognize the crack path and the finite elements pertaining to it. The output of the crack-tracking subroutine is a flag system used to distinguish elements within a crack-path from those outside. In this way, the two groups of elements follow different constitutive behaviour. Here, the continuum damage model presented in Sect. 2 is used for the elements within a crack, while for the rest a linear stress-strain relationship is maintained.

Crack initiation I-monotonic loading
New cracks originate from elements that satisfy the failure condition as defined by Eq. (8). No restrictions are posed to the location of the crack origins in the used mesh as presented in [32]. For cracks initiating at corner or at internal elements, the centroid of the triangular element is selected as the crack origin. For the case of a crack starting at an element with one side on the boundary, the crack origin is located at the midpoint of the boundary edge. A minimum distance is defined between existing and new cracks, referred hereafter as exclusion radius r excl [22]. If two finite elements satisfy the failure criterion at the same time/ load increment and the distance between them is less than the exclusion radius, the one with the highest value of the tensile equivalent stress s þ [see Eq. (4)] is selected as a crack origin. The value of the exclusion radius can be defined according to the structural problem or the fabric of the material analysed as shown in [22,23,32]. Apart from this functionality, the definition of an exclusion radius averts the initiation of secondary spurious cracking surrounding the crack-tip [27,47], a pathology that may spoil the correct damage localization.

Crack propagation and crack intersection
The propagation direction is defined as orthogonal to the direction of the maximum principal effective stresses. This direction is calculated for each element at the crack front (crack tip element) using a stress tensor constructed from the average values of the stresses at the nodes of each element. Starting from the crack entry point, and using a vector with the aforementioned direction (V e in Fig. 4), the next potential element of the crack is identified and the crack entry coordinates are saved. Figure 4a presents an example with the propagation of an internal crack, while in Fig. 4b shows a detail of the crack propagation and the crack entry coordinates on the wake of a crack tip. If the crack origin is an element lying on the boundary of the structure, the crack propagates in a single direction. This is the case of the flexural cracks at the corner of openings or at the base of piers. On the contrary, shear cracks can initiate from the interior of a masonry wall or a spandrel. To allow this, when the tracking algorithm recognizes that the crack origin is an element at the interior of the structure, the propagation of the crack is monitored towards two opposite orientations, as shown in Fig. 4a. In this way, both shear and flexural cracks can be adequately modelled, as in the masonry frame wall presented in Fig. 5.
In structures experiencing multiple cracking, it is unavoidable that cracks will meet other cracks during their propagation. The assumption used in [22,32] that a crack is complete when meeting another crack at its propagation front is reconsidered in this work in order to provide a realist simulation of intersecting cracks in cyclic in-plane loading of shear walls. In the proposed algorithm, a crack continues unimpeded its evolution following the same procedure as if the crack-tip was not a part of another crack. The crack direction is identified as orthogonal to the direction of the maximum principal effective stresses of each finite element and the next elements of the crack are sought. An example of such case is presented in the following section.
In order to avoid sudden changes in the crack direction due to the poor estimation of the stresses arising from the use of standard finite elements, a correction of the crack propagation can be opportunely applied through the use of a maximum curvature criterion [20][21][22]. This criterion compares the crack propagation direction of the crack-tip element V e with the average crack direction V c (Fig. 4a). The latter is obtained as the average sum of the propagation directions of the elements in the crack path lying within a distance lower than r neigh from the crack-tip. If the propagation direction diverges with an angle a larger than a max , then the next element of the crack is sought using the average crack direction V c instead of V e . This correction is illustrated in Fig. 4. The value of the r neigh defines the length of the crack behind the crack tip that will be considered in the computation of the vector V c (see Fig. 4a). The maximum curvature a Tensile damage contour, and b maximum principal strains at the end of the analysis angle a max determines the margins of the allowed crack rotation. The effect of these numerical parameters, as well as the criteria for their determination, are detailed in [22,32].

Crack initiation II-multi-directional cracks
The criteria for crack initiation, presented in Sect. 3.2, define that new cracks can start from a certain distance from existing ones, which is determined by the exclusion radius r excl . This procedure is robust for structural cases with monotonic loading, where the crack directions do not show sudden changes during the loading history (see previous works [22,23,31,32]). However, cyclic loading, such in earthquakes or wind loads, may induce multi-directional cracking to the structure, referring to cracks initiating in the same region of existing cracks but with a different direction. To illustrate this, consider the simplified case of a shear wall (Fig. 6) subjected to cyclic shear loading at its top (restraining the rotation). The top displacement towards the left produces a shear crack, starting from the middle of the wall and propagating towards its two ends (Fig. 6a). The reversal of the loading cycle finds the wall at its initial configuration with the crack closed (Fig. 6b). The loading in the opposite direction results in the opening of a new crack that propagates diagonally towards the other two corners (Fig. 6c). If the first crack is completely closed, the structural condition before the opening of the two cracks will very be similar resulting to the possible origin of the second crack very close to the origin of the first one.
The accurate simulation of such cases needs a different approach regarding crack initiation from the one presented in Sect. 3.2. For the specific case of the shear wall in Fig. 6, this approach must allow the initiation of the second crack at the finite element that satisfies the failure criterion after the reversal of the loading cycle, independently of its distance to the closed crack.
The modelling of cases involving the initiation of cracks due to different stress-states at a region of existing cracks is possible through the categorization of cracks into ''active'' and ''inactive''. Prior to the definition of these two categories, two new versors are introduced: • Crack directionv cr : a versor denoting the average direction of a consolidated crack. The crack direction of a crack is calculated using the crack entry coordinates at each element within the total length of the crack. • Crack closure directionv cl : a versor denoting the average direction orthogonal to the maximum principal effective stress of all the elements within the crack. This versor shows the average crack propagation direction of all the elements within a crack at each time/load increment of the analysis.
These two versors are updated at the end of each time/ load increment for each of the consolidated cracks. A crack i is defined as ''active'' if the angle between its crack direction versorv i cr and its crack closure direction versorv i cl is lower than a predefined limit angle a lim . The crack status is alternatively expressed through the use of the following function H i ½jv i cr Áv i cl j À cosða lim Þ where (Á) is the dot product, the subscript i ¼ 1; N cr is used to denote the versor values for each crack of the total existing cracks N cr and H½ is the Heaviside step function. Following the above, an element will be the crack origin when all of the following criteria are satisfied: (i) The failure condition, as defined by the constitutive model (Eq. 8), is reached. (ii) There is no ''active'' crack within a distance defined by the exclusion radius r excl .
This means that new cracks can initiate at locations where the current stress-state has resulted in the closing of the pre-existing cracks. To illustrate this, Fig. 7

(b) (a) (c)
Active crack Inactive crack Fig. 6 Example of cracking under the shear loading of a wall: a Opening of one crack during the loading towards the left, b closing of the crack when returning to the initial configuration and c opening of a second crack, diagonal to the first during the reverse loading illustrates how the algorithm models the two shear cracks at the middle of the wall for the example of Fig. 6. Figure 7a shows the stress-state within the crack at the moment of the opening of the crack during the movement towards the left. The two versors of the crack (v cr andv cl , shown in one of the two crack-tips of the crack) coincide, as the crack direction has been defined according to the crack-propagation direction of each element within it. The crack is active and can continue its propagation. Figure 7b presents the stress state for a step during the reversal of the cycle. The wall is now deforming towards the right and the existing crack is closed. This can be seen by the direction of the principal effective stresses which have rotated comparing to when the crack was open. The crack closure (v cr ) and crack direction versors (v cl ) form an angle which is greater than the limit angle a lim . The existing crack is therefore identified as inactive, and new cracks are allowed to open. This happens at a later instance, shown in Fig. 7c, when a crack initiates very close to the existing crack and propagates thereafter. The use of a limit angle as a way to identify the closure of existing cracks and the possibility of new ones to open resembles the concept used to model multiple cracks at constitutive level in the smeared crack models by de Borst and Nauta [48]. A similar concept was used also by Song & Belytschko aiming to model crack branching for propagating cracks [47] under dynamic fracture. Likewise, the limit angle introduced in this work (a lim ) is a numerical parameter aiming to deactivate the use of the exclusion radius when the strain state of the crack has exceeded the margins set by Eq. (27) and thus to allow new cracks to emerge close to existing ones. maintained constant during the whole experiment, aiming to represent gravitational loading acting on the top of the masonry wall during a seismic event. The seismic action was reproduced by applying an horizontal displacement at the top of the wall. Proper boundary conditions were applied to avoid its top and bottom rotation. Figure 8a shows a sketch of the experimental set-up adopted during the execution of the cyclical shear tests. The experimentally obtained shear load against displacement graph is presented in Fig. 8b. The wall presents the typical behaviour for shear dominating failure. This is characterised by an overall drop in the capacity with increasing load cycles and appreciable energy dissipation throughout the whole loading history, as represented by the area of the hysteretic cycles of the load-displacement curve. The quasi-brittle character of the shear failure is due to the decrease of strength produced by the formation of two systems of crossing diagonal cracks in the wall [49,50].
In this work, the experiment is reproduced by using the continuum damage model and the crack-tracking algorithm presented in Sects. 2 and 3. The wall has been discretized using 3-noded constant-strain triangles (2836 nodes). The material parameters are presented in Table 3. The compressive strength, the density and the Young's modulus are equal to the values obtained through experimental tests on stack bond prisms reported in [50,51]. The rest of the mechanical properties of the materials, and the parameters of the constitutive model have been chosen in order to calibrate correctly the shear response of the material, on the basis on the discussion presented in Sect. 2.2. To this aim, a value of j 1 ¼ 0:15 has been chosen which is similar to the value used in other available studies on masonry shear walls [19,39,41]. Regarding the crack-tracking parameters, the exclusion radius has been set equal to r excl ¼ 0:20m, which is a value between the length and the width of the masonry units. The limit angle used for the activationdeactivation of the cracks under the cyclic loading is a lim ¼ 40 . This value means that three additional cracks with an orientation varying by 40 can initiate next to an existing crack. To visualize this, and considering that the exclusion radius has been defined equivalent to the units size, the above value implies that a total of four cracks with different orientations can initiate from the same unit. For the investigated case of a shear wall under cyclic loading, this value seems reasonable as it can allow the formation of diagonal cracking due to shear failure. The parameters for the maximum curvature criterion are a max ¼ 25 , r neigh ¼ 0:60 m. These values of the parameters showed to be adequate to avoid sudden changes in the crack direction as reported in previous works of the authors [22].
The wall has been analysed under plane stress conditions and the simulation has been performed in two stages. The first one includes the application of the self-weight and the vertical compressive load at the top of the wall. In the second stage, a cyclic horizontal displacement is applied at the top of the wall with increasing magnitude, as in the experiment. Vertical displacement at the top is restrained, precluding the Fig. 8 a Experimental set-up and b experimentally obtained force against displacement graph of the tested masonry wall. Pictures adapted from [49] rotation of the wall. The numerical solution is performed in an incremental manner. At each time/ load increment, the nonlinear equations are solved with the use of a modified Newton-Raphson (using the secant stiffness matrix). Convergence is achieved when the ratio between the norm of the iterative residual forces and the norm of the total external forces is lower than 10 À2 . Calculations are performed using an enhanced version of the finite element software COMET [52], while pre-and post-processing with GiD [53], both developed in CIMNE, Barcelona. Figure 9 presents the contours of the tensile d þ and compressive d À damage indices, as well as that of the maximum principal strains at the end of the analysis. The wall is dominated by diagonal shear cracking, which is in agreement with the reported experimental damage pattern [49,50]. There is evidence of flexural cracking at the lower and bottom parts of the wall, as well as at the two sides. Despite the opened cracks during the loading history, the main degradation of the wall is concentrated at the two diagonal cracks that have opened in the middle of the wall, as shown by the contour of the maximum principal strains. The use of the crack-tracking algorithm results in a very clear representation of the localized cracks, aiding the identification of the occurring collapse mechanism at the wall. Figure 10a presents the damage evolution and Fig.  10b the state of the cracks for the opposite top displacements of the first cycle. For demonstration purposes, the contour value 1 (red) is used for the active cracks and -1 (white) for the inactive ones, while the elements outside of the crack maintain value equal to 0 (orange) (Fig. 10b). It is visible that all cracks are active for the maximum displacement of the first cycle with a loading direction towards the right. Two flexural cracks at the two corners of the wall consist the first occurring damage. After those, a shear crack initiated at the middle of the wall (left column of Fig. 10). The results with the same top displacement but opposite direction are presented in the right column of Fig. 10. All the previously opened cracks have now closed and are correctly identified as inactive. New horizontal cracks appear at the two opposite corners, as well as a diagonal shear crack at  The proposed approach is capable to simulate both the correct initiation of new cracks close to the existing ones, as well as their intersection. Figure 11 presents the damage pattern during the loading cycle with a maximum displacement at the top equal to d ¼ AE 3.0 mm. After the initial damage shown in Fig. 10, additional shear cracks have opened, with the deformation, however, localizing at the initial diagonal crack at the middle of the wall (see center column of Fig. 11). Some flexural cracks initiate at the two lateral sides of the wall and propagate horizontally towards the interior. Note that the cracks opened during an opposite loading direction have been correctly identified as closed (right column in Fig. 11).

Numerical results
The response of the simulated wall in terms of horizontal load against horizontal displacement is presented in Fig. 12. The numerical simulation captures the shear failure of the wall, which is characterized by an overall softening branch and stiffness degradation after exceeding the maximum capacity. The first cycles of the loading are characterized by an almost linear behaviour, due to the small amount of damage as shown in Fig. 10. A significant drop in the capacity occurs after the complete degradation of the diagonal cracks at the middle of the wall, which occurs during the loading cycles with the maximum displacement of 3.0 mm, as in the experiment. After that, the increase of the displacement results in the drop of the stiffness due to the degradation of the additional shear cracks, yielding the final damage pattern and strain state of Fig. 9. The used constitutive model is able to represent satisfactorily the cyclic damaging behaviour of the wall. The slight underestimation of the magnitude of the experimental envelope of the hysteretic loops is due to the specific nature of the used damage model, with the association of the evolution of the irreversible strains to only that of the compressive damage variable d À . In fact, as can be seen in Fig. 9c, compressive damage d À appears only at the two diagonal cracks at the middle of the wall, making them the only source of the irreversible deformations.

Mesh-dependency study
The tracking algorithm has shown to improve the mesh objectivity in structural analyses of small [22,31,32] and large scale structures [54] with monotonically propagating cracks. In this section, the shear wall is analysed with an additional mesh aiming to study the mesh dependency of the presented approach for the more demanding case of cyclic loading. The first column of Fig. 13 shows the different discretization used in the mesh dependency analyses. The first mesh, denoted hereafter as Mesh-A, is the mesh used for the aforementioned results with 5470 elements. The second mesh, Mesh-B, presents a constant structured pattern in the greatest part of the wall with equilateral triangles making angles of 30 , 90 and 150 with the horizontal axis and has a lower number of elements (3704) compared to Mesh-A. The different columns of Fig. 13 present the damage pattern at the end of the different cycles during the loading history. Diagonal shear cracking dominates the main part of the wall in both cases, while flexural cracks appear at the boundaries. The damage pattern is very similar for the greatest loading history of the wall, with some differences occurring in the last cycles, during which the coarser mesh (Mesh-B) presents some additional shear cracking. Figure 14 illustrates the force against the applied displacement graphs for the two meshes overlapped. The two analyses bare important similarities. The predicted capacity of the wall shows a variation of less than 3% and is attained at the second loading cycle for a displacement of 3.0 mm. Following that, the structure presents a softening response with dropping capacity and energy dissipation due to the further opening of the cracks. As presented in Fig. 13, after the first two loading cycles most of the cracks have already formed, and the diagonal cracks at the middle of the wall as well as the boundary cracks at the top and bottom ends of the wall have completely opened.
Upon this point, the numerical responses reflect the small differences regarding the positions of the cracks and the appearance of more cracks at the coarser mesh (Mesh-B), slightly alternating the residual response of the wall.

Conclusions
This paper has presented the enhancement of the macro-modelling finite element method with a cracktracking algorithm for the simulation of the behaviour of masonry structures under in-plane cyclic loading. A novel methodology has been proposed to simulate intersecting and multi-directional cracks that occur in masonry under cyclic actions with a local cracktracking technique. This has been achieved by setting suitable crack opening/closure criteria in the numerical algorithm. Cracking is represented at a local level using a continuum damage model with distinct representation of tensile and compressive damage. The constitutive model has been formulated to include the description of the evolution of irreversible deformations in a simplified way.
The enhanced macro-model has been used to simulate an experimentally tested shear masonry wall under cyclic quasi-static loading. The tracking algorithm has shown to be capable to simulate shear and flexural cracks propagating in multiple directions and intersecting among them. The results show that the used constitutive model can represent correctly the gradual decrease of the structural capacity during the cyclic loading process. The proposed approach can therefore properly account for complex damage typologies when cracks appear and propagate under shear cyclic condition.
The paper includes a discussion at constitutive level regarding the evolution of the irreversible strains under shear loading. A simple and explicit algorithmic formulation, together with a proper definition of the compressive damage surface, allow to consider the effect of irreversible strains both under shear and compression stress states.