Orthotropic multisurface model with damage for macromechanical analysis of masonry structures

European Journal of Mechanics / A Solids 102 (2023) 105077 Available online 13 July 2023 0997-7538/© 2023 The Authors. Published by Elsevier Masson SAS. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect European Journal of Mechanics / A Solids journal homepage: www.elsevier.com/locate/ejmsol Orthotropic multisurface model with damage for macromechanical analysis of masonry structures C. Gatta∗, D. Addessi Department of Structural and Geotechnical Engineering, Sapienza University of Rome, Via Eudossiana 18, 00184, Rome, Italy ARTICLE INFO Keywords: Masonry Damage Orthotropic response Macromechanical approach Finite element Nonlocal integral regularization ABSTRACT A novel macromechanical model with damage for the analysis of masonry structures in-plane loaded is presented. The model accounts for the directional mechanical properties typically characterizing response of masonry with regular texture. Indeed, the real heterogeneous material is modeled as a fictitious homogenized medium with orthotropic elastic constitutive behavior along the masonry natural axes, identified as the parallel and normal directions to bed joints orientation. The different strength characteristics along each material axis are taken into account by properly defining a damage matrix, which accounts for failure mechanisms due to axial tensile and compressive states, as well as shear. A suitable criterion is introduced, resulting in a damage limit surface geometrically defined in the space of the damage associated variables by the intersection of two ellipsoids and an hyperboloid. The model is implemented into a finite element procedure where the mesh-dependency numerical issue is avoided by adopting a nonlocal integral formulation. Validation examples, involving simple uni-axial and bi-axial tests, as well as more complex loading conditions, are provided to prove the model performances at both material and structural scale. 1. Introduction In the last decades many experimental and numerical studies were devoted to understanding and predicting the response of masonry structures, in view of their seismic assessment. In fact, masonry is the most ancient, but still widely used, construction building material. At the conventional microscopic scale, it is a composite material obtained by assembling blocks, with various nature and shape, by means of mortar layers or dry joints. Geometry, sizes, mechanical properties and arrangement of the constituent materials strongly affect the global structural response. Hence, the most natural and accurate modeling approach appears to be the so called micromechanical strategy, which separately describes each masonry component and, possibly, their interaction behavior (Oliveira and Lourenço,2004;Addessi and Sacco, 2014;Baraldi et al.,2018;Minga et al.,2018;Addessi et al.,2020a; Nocera et al.,2021). Accurate geometric and constitutive descriptions are obtained at the cost of computationally expensive numerical analyses, thus restricting the applicability of such approach to the study of small elements or structural details. Alternately, the scientific community proposed a large variety of continuum models which consider masonry as a homogenized medium where the constituents are no longer distinguishable. The homogenized material is usually modeled by resorting to the classical Cauchy continuum (Lourenço et al.,1997;Berto et al.,2002;Karapitta et al., ∗Corresponding author. E-mail address: [email protected] (C. Gatta). 2011;Pelà et al.,2013;Gatta et al.,2018;Addessi et al.,2021b; Malena et al.,2022), but also the micropolar Cosserat was successfully applied, especially to account for the effect of the characteristic microstructure length on the masonry macroscopic response (Addessi, 2014;Tuna and Trovalusci,2020). Anyway, the correct identification of the constitutive behavior of the homogenized material remains an open issue due to the uncertainty in the calibration of the evolutive laws of the inner variables governing the nonlinear mechanisms. To this end, direct approaches and homogenization and multiscale procedures were applied (Sacco et al.,2018;D’Altri et al.,2020). The latters deduce the material constitutive response of the homogeneous model adopted at the structural scale from the accurate analysis of a properly selected masonry representative volume element (RVE), accounting for the detailed description of components, geometry and arrangement. Weak coupling between the material and structural scale is established if a priori homogenization is performed (Anthoine,1995;Milani et al., 2006), whereas a stronger connection is obtained if step-by-step (Marfia and Sacco,2012;Petracca et al.,2016;Addessi et al.,2020b,2021a) or adaptive (Leonetti et al.,2018) multiscale procedures are adopted. Instead, the direct models calibrate the material properties and evolution laws of the inelastic variables through experimental data on masonry assemblages. https://doi.org/10.1016/j.euromechsol.2023.105077 Received 24 August 2022; Received in revised form 28 June 2023; Accepted 9 July 2023 European Journal of Mechanics / A Solids 102 (2023) 105077 2 C. Gatta and D. Addessi Continuum approaches based on phenomenological constitutive laws capable of describing the main features of the mechanical response without resorting to the nested step-by-step RVE-based homogenization procedures are also referred to as macromechanical models. In such context, several formulations were proposed, including damage models, plasticity models, coupled damage-plasticity models and smeared-crack models (Fusco et al.,2022). Despite most of these approximate the real masonry anisotropic response with the simplified hypothesis of isotropic behavior, these models give a fair compromise between accuracy and computational effort and were successfully applied to analyze large scale structures (Tiberti et al.,2016;Addessi et al.,2021b). However, when dealing with periodic well-organized masonry, the assumption of isotropic response might to be too simplistic, as regular masonry exhibits substantial discrepancy between properties observed in different material directions. The pioneer experimental campaign conducted by Page et al. on running bond panels (Page,1981,1983; Dhanasekar et al.,1985) clearly showed that mortar joints act as plane of weakness and their orientation with respect to the applied loads strongly affects the material strength. Moreover, the anisotropic response emerges already in the elastic range and, usually, reduces to an orthotropic-type. This was proved by the correlation determined experimentally by Cavaleri et al. (2014) between the ratios Young’s modulus-to-Poisson’s coefficient defined along head and bed joints directions, considered as masonry natural axes. Some attempts were made to include the effect of anisotropy in macro-models. For instance, Lourenço et al. (1997) proposed an orthotropic constitutive law fully based on the plasticity theory, which employs a Rankine-type and a Hill-type criterion to simulate tensile and compressive behavior, respectively. Berto et al. (2002) presented an orthotropic damage model for the analysis of brittle masonry subjected to in-plane loading, assuming the masonry natural axes as damage principal axes. They considered the equivalent effective stress measures as damage associated variables, by distinguishing the positive and negative values along the directions parallel and normal to the bed joints. Moreover, for tensile states, they introduced the dependence of the damage evolution parameters on the tensile specific fracture energy. Karapitta et al. (2011) adopted a smeared-crack constitutive model capable of discerning failure modes of unreinforced masonry due to tension normal and parallel to the bed joints, masonry crushing normal and parallel to the bed joints, and masonry shear under compressive vertical stress. Pelà et al. (2013) developed a damage model exploiting the concept of mapped tensors from the anisotropic field to an auxiliary isotropic workspace and, then, combined the model with the crack-tracking technique to reproduce the propagation of localized cracks (Pelà et al.,2014). Similarly, Bilko and Małyszko (2020) established masonry constitutive response in the framework of the elasto-plasticity theory by including a generalization of the Hoffman failure criterion in plane stress state. More recently, Tisserand et al. (2022) formulated an orthotropic thermodynamics-based model including damage, unilateral effect and internal sliding and friction. The reliability of the models described above strongly depends on the adopted failure criteria, whose definition is a hard task, given the complexity of masonry mechanical response. One of the first attempt to identify a proper failure surface for brick masonry under bi-axial stresses dates back to Dhanasekar et al. (1985), which defined the surface as the intersection of three elliptic cones in the space of stresses expressed in the natural axes. Stemming from this proposal, many other formulations were suggested. Berto et al. (2002) modeled the material damage space as a double pyramid with rectangular base in the equivalent effective stress space. Lourenço et al. (1997) composed the limiting surface with Rankin’s and Hill’s yield surfaces. Syrmakezis and Asteris (2001) described mathematically the surface by means of a third order polynomial, which provided satisfactory results in case of both compressive principal stresses. Further developments conducted to more complex failure criteria in order to overcome some limitations of the previous formulations. For instance, Lishak et al. (2012) proposed a very intricate surface shape composed of five parts, each corresponding to different masonry failure mode. Asteris and Plevris (2017) used the neural networks to approximate the limit surface in dimensionless form by obtaining an ‘onion’ shape. Bilko and Małyszko (2020) extended the approach presented in Lourenço et al. (1997) considering two orthotropic Hoffman-type failure criteria. Recently, Malena et al. (2022) modified the isotropic yield function proposed by Bigoni and Piccolroaz (2004) to take into account the masonry anisotropy. Relying on the above considerations, this work presents a novel macromechanical model with damage for the analysis of masonry structures in-plane loaded. The proposed constitutive law introduces an orthotropic description of the material elastic behavior. Then, the stiffness degradation due to cracking, crushing and shear is captured by properly defining a damage matrix, written in terms of independent scalar damage variables, whose evolution is ruled by equivalent strain measures. Moreover, a suitable damage criterion is introduced to account for the variation of the mechanical properties in the different material directions. The failure criterion results into a limit surface geometrically defined by the intersect of two ellipsoids and one hyperboloid in the space of the damage associated variables. The model is implemented in a finite element (FE) code, where the typical meshdependency issue of the numerical solution is overcome by adopting a nonlocal integral formulation. As observed by Bažant and Jirásek (2002), who presented a comprehensive survey of the nonlocal integral procedures applied in the field of plasticity and damage constitutive formulations, nonlocal models were mainly developed to: describe the nonlocal effects in presence of material heterogeneity; regularize the boundary value problem preventing ill posedness in presence of strain-softening and, then, obtain objective numerical solutions; capture size effects observed in experiments and in discrete simulations. Also, according to Bacigalupo and Gambarotta (2012), nonlocal constitutive models permit to include geometric and material length scales to account for the influence of block size and prevent pathological localizations related to the strain-softening nature of the constitutive equations adopted for brittle masonry. They analyzed running bond and English bond masonry by varying the stiffness ratio between brick and mortar and evaluated the characteristic lengths associated to the shear and extensional strains. These resulted as a fraction of the periodic cell size and characterized by different values along the direction parallel and normal to the bed mortar joints. In this work, the integral definition of the strain measures driving the damage variables evolution is introduced to regularize the problem in presence of strain-softening constitutive behavior and guarantee the mesh-independency of the FE results. The size of the region involved in the nonlocal procedure is determined by the nonlocal radius defining the Gaussian weighting function. The paper is organized as follows. Section 2describes the adopted constitutive relationship, focusing on the proposed damage criterion and evolution laws of the damage variables. Section 3provides some computational aspects related to the finite element formulation and regularization technique adopted. Section 4presents the numerical applications. First, the model capability of capturing masonry nonlinear response is evaluated at material level by performing simple uni-axial monotonic and cyclic tests. Then, exploration is moved towards more complex bi-axial loading conditions. Finally, the structural response of masonry walls is investigated comparing the obtained results, in terms of failure mechanisms and global load–displacement response curves, with those recovered from experimental investigations. Section 5concludes with some remarks. In the following the Voigt notation is adopted in the bi-dimensional (2D) framework, representing the second order tensors as 3-component vectors and the fourth order tensors as matrices. European Journal of Mechanics / A Solids 102 (2023) 105077 3 C. Gatta and D. Addessi Fig. 1. Material (𝑇 , 𝑁) and global (𝑋, 𝑌 ) axes of the orthotropic masonry material. Fig. 2. Schematic representation of failure modes associated to (a) 𝐷𝑇 𝑡, (b) 𝐷𝑁𝑡 , (c) 𝐷𝑇 𝑐 , (d) 𝐷𝑁𝑐 and (e) 𝐷𝑇 𝑁 . 2. Damage model 2.1. Damage-based constitutive law To account for masonry anisotropic macroscopic response, the real heterogeneous material is modeled as a fictitious 2D orthotropic medium under plane stress assumptions and the hypothesis of small displacements and strains. The material/intrinsic axes (𝑇 , 𝑁), parallel and normal to the direction of bed joints, are considered as axes of orthotropy. First, the constitutive law is defined in the reference system 𝑇−𝑁, then, it is expressed in the global coordinate system 𝑋−𝑌, by applying standard transformation rules (see Fig. 1). The relation between stresses, 𝜮𝑇 𝑁 , and strains, 𝐄𝑇 𝑁 , referred to the material system, and the corresponding global quantities, 𝜮𝑋𝑌 and 𝐄𝑋𝑌 , classically results as: 𝜮𝑋𝑌 =𝜱𝜮𝑇 𝑁 ,𝐄𝑋𝑌 =𝜳𝐄𝑇 𝑁 ,(1) where 𝜱and 𝜳are the rotation matrices expressed as: 𝜱=⎡⎢⎢⎣ 𝑚2𝑛2−2𝑚𝑛 𝑛2𝑚22𝑚𝑛 𝑚𝑛 −𝑚𝑛 𝑚2−𝑛2⎤⎥⎥⎦ ,𝜳=⎡⎢⎢⎣ 𝑚2𝑛2−𝑚𝑛 𝑛2𝑚2𝑚𝑛 2𝑚𝑛 −2𝑚𝑛 𝑚2−𝑛2⎤⎥⎥⎦ ,(2) being 𝑚= cos 𝜗and 𝑛= sin 𝜗, with the angle 𝜗measured counter clockwise from 𝑋to 𝑇-axis. The stress–strain relationship is defined as: 𝜮𝑇 𝑁 = (𝐈−𝐃)𝐂𝑇 𝑁 (𝐈−𝐃)𝑇𝐄𝑇 𝑁 = 𝐂𝑇 𝑁 𝐄𝑇 𝑁 ,(3) with the 2D strain vector 𝐄𝑇 𝑁 ={𝐸𝑇𝐸𝑁𝛤𝑇 𝑁 }𝑇collecting the axial elongations along 𝑇and 𝑁directions, 𝐸𝑇and 𝐸𝑁, and the shear strain 𝛤𝑇 𝑁 . The stress vector 𝜮𝑇 𝑁 ={𝛴𝑇𝛴𝑁𝛴𝑇 𝑁 }𝑇contains the work-conjugate quantities.  𝐂𝑇 𝑁 is the effective material stiffness matrix deduced by the energy equivalence principle of damage mechanics (Cordebois and Sidoroff,1982), which, unlike the strain equivalence concept, leads to symmetric stiffness matrix for any damage operator (Chow and Wang,1987;Ghrib and Tinawi,1995). According to formulas in (1), the constitutive effective matrix in the global system,  𝐂𝑋𝑌 , is computed as:  𝐂𝑋𝑌 =𝜱 𝐂𝑇 𝑁 𝜳−1 .(4) In Eq. (3) 𝐂𝑇 𝑁 is the orthotropic elastic constitutive matrix of the undamaged material, depending on the Young’s and shear moduli, 𝐸𝑇 𝑇 , 𝐸𝑁𝑁 and 𝐺𝑇 𝑁 , and Poisson’s ratios, 𝜈𝑇 𝑁 ,𝜈𝑁𝑇 ,𝐈is the 3 ×3 identity matrix and 𝐃is the 3 ×3 damage operator which is defined as follows: 𝐃=⎡⎢⎢⎣ 𝐷𝑇0 0 0𝐷𝑁0 0 0 𝐷𝑇 𝑁 ⎤⎥⎥⎦ ,(5) being 𝐷𝑇,𝐷𝑁and 𝐷𝑇 𝑁 three scalar damage variables. It is worth mentioning that principal values of damage are defined in most of anisotropic damage models and, consequently, damage affecting shear components results a proper combination of these quantities (Chow and Wang,1987;Yazdchi et al.,1996;Berto et al.,2002;Williams et al., 2003). Here, following the approach used by Matzenmiller et al. (1995) and Simon et al. (2017) for fiber-reinforced and layered composites, an independent damage variable in shear, 𝐷𝑇 𝑁 , is introduced. This assumption can be justified by the different damaged areas for normal and shear stresses and allows for higher model versatility. As concerns the damage variables 𝐷𝑇and 𝐷𝑁, these are defined on the basis of damage parameters accounting for tensile, 𝐷𝑖𝑡, and compressive, 𝐷𝑖𝑐 , (𝑖=𝑇 , 𝑁) strain states, as follows: 𝐷𝑇=𝛼𝑇𝐷𝑇 𝑡 + (1 − 𝛼𝑇)𝐷𝑇 𝑐 , 𝐷𝑁=𝛼𝑁𝐷𝑁𝑡 + (1 − 𝛼𝑁)𝐷𝑁𝑐 .(6) The weighting coefficients 𝛼𝑇and 𝛼𝑁, defined later on, are introduced to describe the unilateral stiffness recovery due to the re-closure of the tensile cracks under compressive states during cyclic loading histories. In other words, the material degradation is irreversible but its effect on the mechanical response can be activated or inactivated depending on the applied load (Pelà et al.,2013;Addessi,2014;Gatta et al.,2018). According to their physical meaning, all damage parameters, 𝐷𝑖𝑡, 𝐷𝑖𝑐 (𝑖=𝑇 , 𝑁) and 𝐷𝑇 𝑁 , can range between 0 and 1, representing European Journal of Mechanics / A Solids 102 (2023) 105077 4 C. Gatta and D. Addessi Fig. 3. Damage limit surface at the onset of the damaging mechanism: (a) 3D representation, (b) and (c) sections. the undamaged and completely degraded state, respectively. Moreover, the irreversible thermodynamic condition is enforced, such that  𝐷𝑖𝑡 ≥ 0, 𝐷𝑖𝑐 ≥0and  𝐷𝑇 𝑁 ≥0, together with the physical constraint 𝐷𝑖𝑡 ≥𝐷𝑖𝑐 (𝑖=𝑇 , 𝑁). Each damage variable is associated to a distinct failure mode, as schematically shown in Fig. 2 by the typical cracking patterns due to tensile and compressive states along each natural axis, and shear state. Accordingly, associated variables 𝑌𝑖(𝑖=𝑇 , 𝑁, 𝑇 𝑁) are introduced. These are equivalent strain measures ruling onset and evolution of the damage parameters, as clarified later in Section 2.2. The following expressions are assumed: 𝑌𝑇=𝐸𝑇+𝜈𝑁𝑇 𝐸𝑁, 𝑌𝑁=𝐸𝑁+𝜈𝑇 𝑁 𝐸𝑇, 𝑌𝑇 𝑁 =|𝛤𝑇 𝑁 |, (7) where 𝜈𝑁𝑇 = [(1− 𝐷𝑁)∕(1 −𝐷𝑇)]𝜈𝑁𝑇 and 𝜈𝑇 𝑁 = [(1− 𝐷𝑇)∕(1 − 𝐷𝑁)]𝜈𝑇 𝑁 are the degraded Poisson ratios, introduced to independently describe the axial damaging processes under uni-axial stress states along the material axes. From Eq. (7) it is clear that cracking and crushing failure modes, associated to 𝑌𝑇and 𝑌𝑁, depend on normal strains, whereas shear failure is solely controlled by shear deformation. On the basis of quantities in Eq. (7), the weighting coefficients 𝛼𝑇 and 𝛼𝑁in Eq. (6) are expressed as: 𝛼𝑖=𝐻(𝑌𝑖)with 𝑖= (𝑇 , 𝑁),(8) with 𝐻(∙) denoting the Heaviside function (i.e. 𝐻(∙) = 1 if ∙≥0, otherwise 𝐻(∙) = 0). It emerges that the model assumes no crack reclosure effect associated to 𝐷𝑇 𝑁 , as shear damage is caused mainly by transverse cracks which do not close under reversal shear loads. 2.2. Damage onset and evolution The definition of a proper failure criterion and evolution laws for the inelastic variables is a fundamental step to predict load bearing capacity of regular masonry. Inspired by the pioneer work of Dhanasekar et al. (1985), the adopted damage criterion accounts for the material state referring to the natural axes. Indeed, the proposed limit surface is geometrically defined by the intersect of two ellipsoids, 𝐹𝐷 1and 𝐹𝐷 3, and one hyperboloid, 𝐹𝐷 2, in the space of the damage associated variables. Few material parameters are needed to construct the surface, that is the initial uni-axial damage thresholds in the directions tangential, 𝑌𝑇 𝑡0and 𝑌𝑇 𝑐0, and normal, 𝑌𝑁𝑡0and 𝑌𝑁𝑐0, to the bed joints, European Journal of Mechanics / A Solids 102 (2023) 105077 5 C. Gatta and D. Addessi by distinguishing them to account for the non-symmetric behavior in tension and compression (as the subscripts ‘𝑡’ and ‘𝑐’ indicate), and the pure shear threshold 𝑌𝑠0.Fig. 3 shows the 3D representation of the damage limit surface in the positive 𝑌𝑇 𝑁 -semi space (Fig. 3(a)) and two sections corresponding to 𝑌𝑇 𝑁 = 0 and 𝑌𝑇=𝑌𝑁where the mentioned thresholds are indicated (Fig. 3(b,c)). Dashed black lines in Fig. 3(b,c) identify the three regions of the space in which each surface 𝐹𝐷 𝑖(𝑖= 1,2,3) defines the limit function. The resulting domain represents the damage limit surface: inner points correspond to material elastic states, points lying on the boundary indicate the onset of possible damaging mechanisms and require evolution of the surface so that the updated material states belong to the new surface boundary. The ellipsoids 𝐹𝐷 1and 𝐹𝐷 3, ruling states of bi-axial tension and compression coupled to shear, are expressed as: 𝐹𝐷 1=(𝑌𝑇 𝐴1)2 +(𝑌𝑁 𝐵1)2 +(𝑌𝑇 𝑁 𝐶1)2 − 1 ,(9) 𝐹𝐷 3=(𝑌𝛼𝐴 𝐴3)2 +(𝑌𝛼𝐵 𝐵3)2 +(𝑌𝑇 𝑁 𝐶3)2 − 1 ,(10) with: 𝑌𝛼𝐴 = cos 𝛼(𝑌𝑇−𝑂𝑇)+ sin 𝛼(𝑌𝑁−𝑂𝑁), 𝑌𝛼𝐵 = − sin 𝛼(𝑌𝑇−𝑂𝑇)+ cos 𝛼(𝑌𝑁−𝑂𝑁).(11) At the onset of the damaging process, 𝐴1=𝑌𝑇 𝑡0,𝐵1=𝑌𝑁𝑡0,𝐵3= √𝑌2 𝑇 𝑐0+𝑌2 𝑁𝑐0∕2 and 𝐴3=𝛽 𝐵3, with 𝛽the material parameter affecting the shape of the limit surface in compression (Fig. 3(b,c)). In Eq. (11), 𝑂𝑇and 𝑂𝑁are the coordinates of the central point 𝐎=(𝑂𝑇, 𝑂𝑁,0)of 𝐹𝐷 3and 𝛼is the angle defined on the basis of the uni-axial compressive thresholds, as illustrated in Fig. 3(b). Finally, quantities 𝐶1and 𝐶3in Eqs. (9) and (10) are determined so as to properly connect the two ellipsoids to 𝐹𝐷 2. This latter is defined as: 𝐹𝐷 2=𝐴2𝑌2 𝑇+𝐵2𝑌2 𝑁+𝐶2𝑌2 𝑇 𝑁 −𝐷2𝑌𝑇𝑌𝑁+𝐸2𝑌𝑇+𝐹2𝑌𝑁− 1 ,(12) with 𝐴2=1 𝑌𝑇 𝑡0𝑌𝑇 𝑐0 , 𝐵2=1 𝑌𝑁𝑡0𝑌𝑁𝑐0 , 𝐶2=1 𝑌2 𝑠0 , 𝐷2=(1 𝑌𝑇 𝑡0𝑌𝑇 𝑐0 +1 𝑌𝑁𝑡0𝑌𝑁𝑐0), 𝐸2=𝑌𝑇 𝑐0−𝑌𝑇 𝑡0 𝑌𝑇 𝑐0𝑌𝑇 𝑡0 , 𝐹2=𝑌𝑁𝑐0−𝑌𝑁𝑡0 𝑌𝑁𝑐0𝑌𝑁𝑡0 . (13) Region outside the damage domain represents no admissible material states. The limit surface evolution, satisfying admissibility conditions of the material states, follows a not homothetic transformation, as illustrated for some examples in Fig. 4. Here, for the sake of simplicity, null shear strains are considered and the surface evolution is plotted only in the 𝑌𝑇−𝑌𝑁plane. Extension to the 3D case is straightforward. In Fig. 4(a) the black square point denotes a case for which evolution of damage surface occurs. The point is characterized by 𝑌𝑇>0, 𝑌𝑁<0and 𝑌𝑇 𝑁 = 0 and, consequently, only the tangential tensile, 𝑌𝑇 𝑡0, and normal compressive, 𝑌𝑁𝑐0, thresholds evolve through the parameter 𝐴(with 𝐴≥1) properly determined, keeping constant 𝑌𝑇 𝑐0, 𝑌𝑁𝑡0and 𝑌𝑠0. The corresponding limit surface after evolution is shown in Fig. 4(b) with dashed lines. Obviously, various combination of strain states can occur, leading to different types of domain evolution. For instance, in Fig. 4(c) and (d) the cases of bi-axial tension and uni-axial compression along 𝑇-axis are illustrated. The solution of the nonlinear evolution problem of damage variables requires to distinguish between damage in tension and compression along each material axis on the basis of the sign of 𝑌𝑖(𝑖=𝑇 , 𝑁), and evaluate the current damage thresholds, 𝑌𝑇0,𝑌𝑁0and 𝑌𝑇 𝑁0, following the procedure illustrated in Fig. 4(a). The threshold quantities coincide with the coordinates of the point given by the intersection between the damage surface and the line connecting the axes origin and the material point exceeding the surface and, consequently, these can vary during the loading history also due to the evolution of the damage limit domain. Once evaluated the thresholds, the following evolutive rules are adopted:  𝐷𝑖𝑡 =𝑌𝑖−𝑌𝑖0 𝑌𝑖+𝑏𝑖𝑡 𝑌𝑖0 (𝑌𝑖≥0), 𝐷𝑖𝑐 =|𝑌𝑖|−|𝑌𝑖0| |𝑌𝑖|+𝑏𝑖𝑐 |𝑌𝑖0|(𝑌𝑖<0),  𝐷𝑇 𝑁 =𝑌𝑇 𝑁 −𝑌𝑇 𝑁0 𝑌𝑇 𝑁 +𝑏𝑇 𝑁 𝑌𝑇 𝑁0 , (14) with 𝑖=𝑇 , 𝑁 and 𝑏𝑇 𝑡,𝑏𝑁𝑡,𝑏𝑇 𝑐 ,𝑏𝑁 𝑐 ,𝑏𝑇 𝑁 material parameters governing the growth rate of the damaging processes: the higher the parameters, the slower the damage progression and greater the resistance and fracture energy density of the material. As an example, Fig. 5 shows the effect of the parameter 𝑏𝑇 𝑡 on both the evolution of the corresponding damage parameter with respect to the associated variable (Fig. 5(a)) and the uni-axial tensile stress–strain relationship (Fig. 5(b)). Detailed description of the time discretization algorithm developed for the damage evolution problem is provided in Section 3.2. To summarize, the proposed model requires to define, in addition to the elastic parameters, six parameters related to the damage criterion (𝑌𝑇 𝑡0, 𝑌𝑇 𝑐0, 𝑌𝑁𝑡0, 𝑌𝑁𝑐0, 𝑌𝑠0, 𝛽) and five parameters associated to the evolutive laws of the damage variables (𝑏𝑇 𝑡, 𝑏𝑇 𝑐 , 𝑏𝑁𝑡, 𝑏𝑁 𝑐 , 𝑏𝑇 𝑁 ). The resulting number of nonlinear parameters to be identify is comparable, and generally lower, to that required by other macromodels which take into account the anisotropic nonlinear response of masonry (Lourenço et al.,1997;Berto et al.,2002;Karapitta et al.,2011;Pelà et al.,2013). The flexibility of the model in reproducing various shapes of the constitutive stress–strain relationships could be further improved. In fact, although the peak strengths are strongly affected by the defined damage thresholds, a limit of the current formulation is that both the peak strengths and the fracture energy densities depend on the 𝑏𝑖 (𝑖=𝑇 𝑡, 𝑇 𝑐, 𝑁𝑡, 𝑁𝑐, 𝑇 𝑁) parameters. Hence, the evolutive laws of the damage variables could be modified to make the definition of strengths and fracture energies properties independent, but this would increase the number of parameters to be defined and the model complexity. 3. Computational aspects 3.1. Nonlocal FE formulation The model presented in the previous section was introduced in a 4-node isoparametric quadrilateral element, and implemented in a finite element procedure based on the classical displacement-based formulation. Each node is provided with two displacement degrees of freedom, and bi-linear interpolation functions are used for the two translation fields, 𝑈𝑋and 𝑈𝑌. As known, when modeling response of quasi-brittle materials characterized by strain-softening constitutive behavior, like masonry, the pathological mesh-dependence problem arises. Indeed, the strain may localize into narrow bands, whose width depends on the finite element size. To overcome this numerical issue several strategies were proposed based on the fracture energy concept (Govindjee et al.,1995;Petracca et al.,2016;Di Re et al.,2018), higher-order formulations (Addessi et al.,2002), Cosserat continuum (De Borst,1991;Addessi,2014) and nonlocal integral approach (Jirásek,1998;Marfia and Sacco,2012; Gatta et al.,2018). Among these, the last mentioned technique is adopted in this study, as it allows to obtain objective numerical results and relies on mechanical evidences. Giving up the principle of local action, it is assumed that the degrading process at each material point is influenced by the mechanical state of the points lying in a properly defined neighborhood. Hence, the integral definition of the damage associated variables in Eqs. (7) is introduced as:  𝑌𝑖(𝐗)=1 ∫𝛺𝜓(𝐗,𝐒)𝑑𝛺 (𝐒)∫𝛺 𝑌𝑖(𝐒)𝜓(𝐗,𝐒)𝑑𝛺 (𝐒) 𝑖=𝑇 , 𝑁, 𝑇 𝑁 , (15) European Journal of Mechanics / A Solids 102 (2023) 105077 6 C. Gatta and D. Addessi Fig. 4. Evolution of the damage surface considering different material states: (a,b) combined tension along 𝑇-axis and compression along 𝑁-axis, (c) bi-axial tension and (d) uni-axial compression. Solid and dashed lines refer to the initial and expanded damage surface, respectively. Fig. 5. Effect of the material parameter 𝑏𝑇 𝑡 on: (a) evolution of the damage parameter with respect to the corresponding associated variable and (b) stress–strain relationship. European Journal of Mechanics / A Solids 102 (2023) 105077 7 C. Gatta and D. Addessi Table 1 Main steps of the solution of the damage evolution problem and constitutive law at the Gauss point of the FE. 1. Compute strains 𝐄𝑘 𝑋𝑌 starting from displacements 𝐔𝑒𝑘 𝑋𝑌 : 𝐄𝑘 𝑋𝑌 =𝐋𝑒𝐔𝑒𝑘 𝑋𝑌 2. Project strains 𝐄𝑘 𝑋𝑌 to the material coordinate system: 𝐄𝑘 𝑇 𝑁 =𝜳−1𝐄𝑘 𝑋𝑌 3. Calculate the local and nonlocal damage associated variables by using Eqs. (7) and (15) 4. Evaluate the damage limit functions according to Eqs. (9)–(12) 5. Determine the damage thresholds following the example procedure illustrated in Fig. 4(a) 6. Solve damage evolution problem according to Table 2 and define damage matrix 𝐃𝑘as in Eq. (5) 7. Update the limit surface after damage following the example procedure illustrated in Fig. 4 8. Compute damaged stiffness matrix  𝐂𝑘 𝑇 𝑁 and stresses 𝜮𝑘 𝑇 𝑁 :  𝐂𝑘 𝑇 𝑁 = (𝐈−𝐃𝑘)𝐂𝑇 𝑁 (𝐈−𝐃𝑘)𝑇 𝜮𝑘 𝑇 𝑁 = 𝐂𝑘 𝑇 𝑁 𝐄𝑘 𝑇 𝑁 9. Evaluate  𝐂𝑘 𝑋𝑌 and stresses 𝜮𝑘 𝑋𝑌 :  𝐂𝑘 𝑋𝑌 =𝜱 𝐂𝑘 𝑇 𝑁 𝜳−1 𝜮𝑘 𝑋𝑌 =𝜱𝜮𝑘 𝑇 𝑁 being  𝑌𝑖the nonlocal quantities at point 𝐗, evaluated on the basis of the corresponding local variables 𝑌𝑖at points placed in its neighborhood on the surface 𝛺. The influence on 𝐗of the point 𝐒is weighted by means of the classical Gaussian function 𝜓, which, in turn, depends on the nonlocal radius 𝐿𝑐, as follows: 𝜓(𝐗,𝐒)=𝑒 −(‖𝐗−𝐒‖ 𝐿𝑐)2 .(16) After computating the integral quantities in Eq. (15), these are introduced in the limit functions in Eqs. (8)–(12) and, in case, are used to solve the evolution problem of the damage variables at each integration point of the FE discretized problem, according to Eqs. (14). 3.2. Solution algorithm Table 1 summarizes the main steps involved in the solution of the damage evolution problem at the typical Gauss point of each FE, making explicit the passage from the global to the material coordinate system and vice versa. Reference is made to the generic Newton– Raphson iteration ‘𝑘’ of the current time step ‘𝑡𝑛+1’ within the global solving algorithm. In short words, the FE program provides the global nodal displacement vector 𝐔𝑘 𝑋𝑌 , from which the element displacements 𝐔𝑒𝑘 𝑋𝑌 are extracted. Then, the strains 𝐄𝑘 𝑋𝑌 are computed through the compatibility matrix 𝐋𝑒=𝐁𝐍𝑒, obtained by applying the standard compatibility operator 𝐁to the shape function matrix 𝐍𝑒. Afterwards, the strain vector 𝐄𝑘 𝑇 𝑁 and the nonlocal damage associated variables  𝑌𝑘 𝑖(𝑖=𝑇 , 𝑁, 𝑇 𝑁) are evaluated. On the basis of these, the damage evolution problem is solved by computing the damage limit functions and each damage parameter according to the following general form: 𝐷𝑘 𝑖=𝐷𝑛 𝑖+𝛥𝐷𝑘 𝑖𝑖= (𝑇 𝑡, 𝑇 𝑐, 𝑁𝑡, 𝑁𝑐, 𝑇 𝑁),(17) with 𝐷𝑛 𝑖denoting the damage value at the previous time step 𝑡𝑛and 𝛥𝐷𝑘 𝑖 the damage increment evaluated at the current time 𝑡𝑛+1 and iteration 𝑘 as reported in Table 2, where apex ‘𝑛+1’ is omitted for simplicity. From Table 2, it emerges that each damage increment is computed using the current damage thresholds, 𝑌𝑘 𝑇0,𝑌𝑘 𝑁0and 𝑌𝑘 𝑇 𝑁0, and the material parameters 𝑏𝑖(𝑖=𝑇 𝑡, 𝑁𝑡, 𝑇 𝑐, 𝑁𝑐, 𝑇 𝑁). Finally, the solving algorithm ends with the computation of the stress vector and the effective stiffness matrix referred to the material system and, after, to the global reference system, as detailed in Table 1. To be noted is that the secant stiffness matrix is adopted in the Newton– Raphson procedure, as the material tangent stiffness is cumbersome to obtain within nonlocal integral formulations (Jirásek and Patzák, 2002). 4. Model validation In this section, the presented model is employed to analyze masonry response both at material and structural level. First, simple uni-axial Table 2 Solution algorithm for the evolution laws of the damage variables. IF 𝐹𝐷𝑘 ℎ<0 (ℎ= 1,2,3) 𝛥𝐷𝑘 𝑇 𝑡 =𝛥𝐷𝑘 𝑇 𝑐 =𝛥𝐷𝑘 𝑁𝑡 =𝛥𝐷𝑘 𝑁𝑐 =𝛥𝐷𝑘 𝑇 𝑁 = 0 ELSE 𝛥𝐷𝑘 𝑇 𝑁 =  𝑌𝑘 𝑇 𝑁 −𝑌𝑘 𝑇 𝑁0  𝑌𝑘 𝑇 𝑁 +𝑏𝑇 𝑁 𝑌𝑘 𝑇 𝑁0 →𝐷𝑘 𝑇 𝑁 =𝑚𝑖𝑛 (𝐷𝑛 𝑇 𝑁 +𝛥𝐷𝑘 𝑇 𝑁 ,1) IF  𝑌𝑘 𝑇>0THEN 𝐷𝑘 𝑇 𝑐 =𝐷𝑛 𝑇 𝑐 𝛥𝐷𝑘 𝑇 𝑡 =  𝑌𝑘 𝑇−𝑌𝑘 𝑇0  𝑌𝑘 𝑇+𝑏𝑇 𝑡 𝑌𝑘 𝑇0 →𝐷𝑘 𝑇 𝑡 =𝑚𝑎𝑥 (𝐷𝑘 𝑇 𝑐 , 𝑚𝑖𝑛 (𝐷𝑛 𝑇 𝑡 +𝛥𝐷𝑘 𝑇 𝑡,1)) ELSE 𝛥𝐷𝑘 𝑇 𝑐 =| 𝑌𝑘 𝑇|−|𝑌𝑘 𝑇0| | 𝑌𝑘 𝑇|+𝑏𝑇 𝑐 |𝑌𝑘 𝑇0|→𝐷𝑘 𝑇 𝑐 =𝑚𝑖𝑛 (𝐷𝑛 𝑇 𝑐 +𝛥𝐷𝑘 𝑇 𝑐 ,1) 𝐷𝑘 𝑇 𝑡 =𝑚𝑎𝑥 (𝐷𝑘 𝑇 𝑐 , 𝐷𝑛 𝑇 𝑡) END IF  𝑌𝑘 𝑁>0THEN 𝐷𝑘 𝑁𝑐 =𝐷𝑛 𝑁𝑐 𝛥𝐷𝑘 𝑁𝑡 =  𝑌𝑘 𝑁−𝑌𝑘 𝑁0  𝑌𝑘 𝑁+𝑏𝑁𝑡 𝑌𝑘 𝑁0 →𝐷𝑘 𝑁𝑡 =𝑚𝑎𝑥 (𝐷𝑘 𝑁𝑐 , 𝑚𝑖𝑛 (𝐷𝑛 𝑁𝑡 +𝛥𝐷𝑘 𝑁𝑡,1)) ELSE 𝛥𝐷𝑘 𝑁𝑐 =| 𝑌𝑘 𝑁|−|𝑌𝑘 𝑁0| | 𝑌𝑘 𝑁|+𝑏𝑁𝑐 |𝑌𝑘 𝑁0|→𝐷𝑘 𝑁𝑐 =𝑚𝑖𝑛 (𝐷𝑛 𝑁𝑐 +𝛥𝐷𝑘 𝑁𝑐 ,1) 𝐷𝑘 𝑁𝑡 =𝑚𝑎𝑥 (𝐷𝑘 𝑁𝑐 , 𝐷𝑛 𝑁𝑡) END END stress tests are performed, then, the exploration is moved towards more complex bi-axial loading conditions and, finally, structural applications on shear walls are presented. The numerical results are validated by comparison with experimental outcomes. 4.1. Uni-axial monotonic and cyclic behavior Uni-axial tests are performed to show the reliability of the proposed model in describing masonry orthotropic response under monotonic and cyclic loads. Numerical simulations are performed by adopting material parameters contained in Table 3 and setting 𝛽= 1.5and 𝑏𝑖= 1.5(𝑖=𝑇 𝑡, 𝑁𝑡, 𝑇 𝑐, 𝑁𝑐, 𝑇 𝑁). The first example explores the response of a masonry element subjected to uni-axial horizontal tension (i.e. tensile load acting along the global 𝑋-axis applied by a displacement controlled procedure), considering three values of the orthotropy angle 𝜗, that is 𝜗= 0◦, 𝜗= 45◦and 𝜗= 90◦. Numerical results, in terms of stress–strain relationship, are shown in Fig. 6(a) and prove the model capability European Journal of Mechanics / A Solids 102 (2023) 105077 8 C. Gatta and D. Addessi Table 3 Material parameters for uni-axial tests in Figs. 6,7and 8. Elastic parameters Damage parameters 𝐸𝑇 𝑇 [MPa]𝐸𝑁𝑁 [MPa]𝜈𝑇 𝑁 𝐺𝑇 𝑁 [MPa]𝑌𝑇 𝑡0𝑌𝑁𝑡0𝑌𝑇 𝑐0𝑌𝑁𝑐0𝑌𝑠0 4000 2000 0.1 1500 9.95E−05 9.95E−05 7.46E−04 1.99E−03 3.33E−04 Table 4 Material parameters for bi-axial tests in Fig. 9. Elastic parameters Damage parameters 𝐸𝑇 𝑇 [MPa]𝐸𝑁𝑁 [MPa]𝜈𝑇 𝑁 𝐺𝑇 𝑁 [MPa]𝑌𝑇 𝑡0𝑌𝑁𝑡0𝑌𝑇 𝑐0𝑌𝑁𝑐0𝑌𝑠0 5700 5600 0.19 2350 6.8E−05 4.1E−05 6.8E−04 1.2E−03 1.9E−04 Fig. 6. Uni-axial tensile response for different values of the orthotropy angle 𝜗: (a) stress–strain relationships, (b) variations of the damage variables 𝐷𝑇,𝐷𝑁and 𝐷𝑇 𝑁 with respect to the strain 𝐸𝑋. of accounting for orientation of the applied load with respect to bed joints direction. Indeed, different initial elastic stiffnesses and maximum strengths are obtained for the three values of the 𝜗angle adopted. This is a consequence of the stress and strain fields acting in the material axis system, which cause activation of different damaging mechanisms. With reference to Fig. 6(b), it appears that the damage parameter 𝐷𝑇 𝑁 is activated only in case of 𝜗= 45◦, as this is related to the shear strain 𝛤𝑇 𝑁 . Conversely, when 𝜗= 0◦or 𝜗= 90◦, damage 𝐷𝑇 𝑁 disappears and only damage variables 𝐷𝑇and 𝐷𝑁arise. In particular, 𝐷𝑇starts and evolves when 𝜗= 0◦since the 𝑇-axis coincides with the 𝑋-axis, whereas 𝐷𝑁appears in case of 𝜗= 90◦as a consequence of the 𝑋and 𝑁-axis overlap. As evident from Fig. 6(b), the three damage variables, 𝐷𝑇,𝐷𝑁and 𝐷𝑇 𝑁 , evolve in the same way in case of 𝜗= 45◦. This is a special condition due to the parameters chosen to rule the damage evolution, that is 𝑏𝑇 𝑡 =𝑏𝑁𝑡 =𝑏𝑇 𝑁 . Indeed, different evolutive laws and constitutive responses could be obtained by removing these assumptions. To clarify, Fig. 7(a) shows the stress–strain relationship and the damage variations corresponding to 𝑏𝑇 𝑡 =𝑏𝑁𝑡 = 1.5and 𝑏𝑇 𝑁 = 2.5. It emerges that the evolution of 𝐷𝑇 𝑁 differs from that of 𝐷𝑇and 𝐷𝑁. Finally, completely different damage increments appear in the most general case, as testified in Fig. 7(b), where results obtained assuming 𝑏𝑇 𝑡 = 2.5, 𝑏𝑁𝑡 = 1.5,𝑏𝑇 𝑁 = 3 are reported. In these monotonic tests, damages 𝐷𝑇and 𝐷𝑁represent the material degradation caused by tensile load and, accordingly, always correspond to 𝐷𝑇 𝑡 and 𝐷𝑁𝑡, respectively. Instead, when dealing with the cyclic response, it is useful decompose the damage parameters in their tensile and compressive part, especially if the re-closure crack phenomenon is to be analyzed. As an example, Fig. 8(a) shows the stress–strain relationship obtained by applying the deformation history at the top of Fig. 8(b) to the masonry specimen with horizontal bed joints, i.e. 𝜗= 0◦(it is assumed 𝑏𝑇 𝑡 =𝑏𝑇 𝑐 = 1.5). It can be noticed that the stiffness recovery occurs when passing from tension to compression (A-B phase), as this is related to the re-closure of the tensile cracks under compressive state. The subsequent reloading in tension, which leads to point C, is slightly affected by the accumulated compressive damage, because of the constraint 𝐷𝑇 𝑡 ≥𝐷𝑇 𝑐 . The phenomenon is clearly illustrated in the lower part of Fig. 8(b), where the variations of the damage variables 𝐷𝑇 𝑡,𝐷𝑇 𝑐 and 𝐷𝑇are plotted with respect to the fictitious time variable with red, green and dashed black curve, respectively. Apparently, 𝐷𝑇assumes the same value of 𝐷𝑇 𝑡 for tensile states and, then, when a reversal strain occurs, returns equal to 𝐷𝑇 𝑐 , allowing a proper representation of the unilateral damage recovering upon load reversal. 4.2. Bi-axial response: comparison with experimental data The experimental data provided by Page (1981,1983) are used as reference solutions to study the in-plane anisotropic response of masonry. Masonry panels, made of half-scale solid clay units arranged in running bond, were tested under bi-axial loading conditions by using a proper device to impose uniform stress states. The applied loads were oriented at various angle 𝜗with respect to bed joints and the resulting failure surfaces were obtained in terms of principal stresses and their orientation to the bed joints. Test results proved that mortar joints act as planes of weakness, causing distinct directional properties. European Journal of Mechanics / A Solids 102 (2023) 105077 9 C. Gatta and D. Addessi Fig. 7. Uni-axial tensile stress–strain relationships and variations of the damage variables for 𝜗= 45◦: (a) 𝑏𝑇 𝑡 =𝑏𝑁𝑡 = 1.5and 𝑏𝑇 𝑁 = 2.5, (b) 𝑏𝑇 𝑡 = 2.5,𝑏𝑁𝑡 = 1.5,𝑏𝑇 𝑁 = 3. Fig. 8. Uni-axial cyclic response for 𝜗= 0◦: (a) stress–strain law, (b) evolution of the strain 𝐸𝑋and damage variables 𝐷𝑇,𝐷𝑇 𝑡 and 𝐷𝑇 𝑐 during the loading history. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) This clearly emerges from Fig. 9(a–c), where the experimental failure domains are depicted with dots for cases of 𝜗= 0◦,22.5◦and 45◦. Results referred to 𝜗= 67.5◦and 𝜗= 90◦are implicitly contained in those of 𝜗= 22.5◦and 𝜗= 0◦, respectively. Fig. 9(a–c) also shows the numerical failure surfaces (solid lines) derived with the proposed model adopting the material parameters contained in Table 4, set according to Page (1981,1983) and Page et al. (1985). For the analyses is assumed 𝛽= 1.5and 𝑏𝑖= 2 (𝑖=𝑇 𝑡, 𝑁𝑡, 𝑇 𝑐, 𝑁𝑐, 𝑇 𝑁). On the overall, a pretty good agreement emerges between experimental and numerical outcomes. In fact, according to the experimental data, the numerical surface exhibits a non-symmetric shape in case of 𝜗= 0◦, being the compressive strengths normal and parallel to bed joints significantly different. Then, by increasing 𝜗, the shape of the failure domain varies. Despite some slight discrepancies between experimental and numerical data, the proposed model well describes the symmetric shape characterizing the case of 𝜗= 45◦. These results clearly testify the model capability of accounting for the bed joints orientation, thus describing in a phenomenological way the preferential direction of microcracks evolution due to the spatial arrangement of mortar and bricks. 4.3. Numerical and experimental response of shearing walls To explore the capability of the model of reproducing response of masonry structural elements, the panels experimentally tested by Raijmakers and Vermeltfoort (1992) are numerically studied. The walls were built by assembling 18 courses of solid bricks with dimensions 210 × 52 × 100 mm3and 10 mm thick mortar. Only 16 courses were activated, thus resulting in overall width 𝑊= 990 mm and height 𝐻= 1000 mm (Fig. 10). The experiments involved two phases: first, a vertical pressure 𝑝was applied on the top side, then, a monotonically increasing horizontal displacement 𝑠was imposed through a steel beam, preventing any vertical movement of the upper boundary (Fig. 10). Four specimens, labeled as JD walls, were tested assuming different levels of the compression load, i.e. 𝑝= 0.3 MPa for J4D and J5D, 𝑝= 1.21 MPa for J6D and 𝑝= 2.12 MPa for J7D. To perform the numerical simulations, the effective elastic properties of the material were derived via a homogenization procedure. This is based on the selection of a masonry unit cell (UC) representative of the regular arrangement (UC in Fig. 10), both in terms of geometric characteristics and constitutive properties of the components. The cell