Simulating masonry wall behaviour using a simplified micro-model approach

In this paper, a simplified Micro-model approach utilising a combination of plasticity-based constitutive models and the extended finite element method (XFEM) is proposed. The approach is shown to be an efficient means of simulating the three-dimensional non-linear behaviour of masonry under monotonic in-plane, out of plane and cyclic loads. The constitutive models include surface-based cohesive behaviour to capture the elastic and plastic behaviour of masonry joints and a Drucker Prager (DP) plasticity model to simulate crushing of masonry under compression. The novel use of XFEM in simulating crack propagation within masonry units without initial definition of crack location is detailed. Analysis is conducted using standard finite element software (Abaqus 6.13) following a Newton Raphson algorithm solution without employing user-defined subroutines. The capability of the model in terms of capturing non-linear behaviour and failure modes of masonry under vertical and horizontal loads is demonstrated via comparison with a number of published experimental studies.


Introduction:
Masonry is one of the oldest and most widespread structural materials; it has been and is still used for various construction purposes. Masonry consists of units and mortar, these constituents have their own mechanical properties and their geometry and arrangement can vary forming different masonry assemblages. Thus masonry is classified as a heterogeneous anisotropic material, and analysis, understanding and capture of the structural behaviour of masonry is therefore complex. For design of non-standard masonry structures or assessment of existing structures, recourse to numerical modelling is often required to understand the structural behaviour under various loading conditions. Nowadays, numerical models offer a viable alternative to physical experiments. Different numerical methods such as the finite element method (FEM), discrete element method (DEM), limit analysis [1] and [2] and the applied element method (AEM) [3] have been employed to conduct numerical analysis and simulate linear and non-linear behaviour of masonry. The finite element method (FEM) is the focus of this paper. FEM for masonry is based on two main modelling approaches, namely, Micro-modelling and Macro-modelling, the choice depending on the level of accuracy and detail required.
In the Micro-model approach, the simulation can be detailed; the units and mortar are modelled as continuum elements and unit-mortar interfaces are modelled as discontinuum elements. The detailed Micro-model Fig. 1 (a) can provide accurate results, but it is computationally intensive and thus limited to simulating relatively small masonry elements.
Alternatively, a simplified Micro-modelling approach Fig. 1 (b) can be adopted to address the disadvantages of the detailed micro approach. In the simplified approach, the units are expanded by adding the mortar thickness, the expanded units are modelled as a series of continuum elements and the interaction between the expanded units is modelled as series of discontinuum elements.
In the Macro-model approach, Fig. 1 (c), the masonry is considered as a homogenous material with no distinction between units and mortar, the material properties are obtained from average properties of masonry constituents and the masonry is modelled as a series of continuum elements [4]. This approach is adopted where relatively larger and more complex masonry structures are modelled and the global behaviour is of interest, but it cannot capture detailed failure modes. Over the past four decades, finite element techniques have continuously evolved to capture the complex structural behaviour of masonry walls and associated structures. Arya and Hegemier [5] and Page [6] attempted to model masonry using a simplified Micro-modelling approach by taking masonry units as continuum elements and mortar joints as interface elements. This approach was then adopted by Lotfi and Shing [7]  Shing and Cao [9] conducted finite element analysis for partially grouted masonry shear walls, a smeared crack model was adopted to simulate the fracture behaviour of the masonry units and plasticity-based interface elements were used to simulate the mortar joint responses under tensile and shear stress. Although, the model successfully simulated the failure modes of the masonry walls, the lateral resistance of the walls was higher than the resistance obtained from experiments. For example in one of the models reported, the lateral resistance from the numerical analysis was 60% higher than the experimental results. Sutcliffe, et al.
[10] conducted a lower bound limit analysis for masonry shear walls where both tensile and shear failure in the brick units and a compression cap for the interface elements were included, but material softening behaviour was ignored. Citto [11] and Kumar et al. [12] developed an interface model to simulate initiation and propagation of cracks in masonry joints and potential vertical cracks in the middle of masonry units under normal and shear stresses, a compression cap was also included in the model to simulate the plastic response under compression. The proposed model in [11] and [12] was analysed via Abaqus by making use of a user defined subroutine, which defined the constitutive behaviour. In all the aforementioned studies, simplified 2D Micro-models were used to simulate only in-plane behaviour of masonry under normal and shear stresses.
The aforementioned studies have dealt with monotonic in-plane load regimes. With regards to modelling cyclic in-plane loading of masonry, a study conducted by Oliveira and Lourenço [13], proposed a 2D model to simulate the behaviour of masonry under cyclic loadings using interface elements between masonry units. In a more recent study conducted by Miglietta et al. [14], finite/discrete element modelling (DEM) is implemented to simulate the behaviour of Furthermore, 3D models are necessary to simulate reinforced masonry walls because simulation is either not possible or limited in 2D FE analysis. A few of the proposed 3D models available in the literature such as in [15] and [17] have relied on employing a user subroutine and an explicit dynamic analysis procedure. In addition, it can be noted that the crack propagation in the brick units, which plays an important role in non-linear degradation of masonry assemblages, is either ignored or defined via interface elements. In the latter case, the formation of potential cracks has always been assumed to be vertical in the middle of units.
This paper presents a simplified 3D Micro-modelling approach to simulating the 3D

Modelling approach:
Herein the constitutive models used to simulate 3D masonry under the simplified modelling approach are described in detail. In addition, the failure modes associated with the models are also presented.

Surface-based cohesive behaviour model
This surface-based cohesive model is employed to obtain the structural response of masonry along bed and head joints. In other words, linear and fracture behaviour of joints which is based on traction separation behaviour between masonry units, is captured. Thus, failure modes of masonry joints, namely, tensile cracking of the joints Fig. 2 (a) and shear sliding of the joints Fig. 2 (b) and (c) are simulated. So, the equivalent stiffness for joint interfaces is expressed as a function of the mortar's and unit's moduli of elasticity, and the thickness of the mortar eq. 2 and eq. 3 [4].

Plastic response of the joint interfaces:
The initial linear response of the joints is followed by crack propagation. When the damage initiation criterion is achieved based on the user defined tractions between the masonry interfaces i.e. shear and tensile strength of the joints, cracking propagates in the masonry joints. The quadratic stress criterion is used to define damage initiation; this criterion is met when the quadratic stress ratios of masonry interfaces are equal to one. This criterion is adopted as it effectively predicts the damage initiation of joints subjected to mixed-mode loadings [18], which is the case in masonry joint interfaces (the masonry joint interfaces are subjected to tensile stress in the normal direction and shear stress in the two shear directions).
The criterion is expressed as in eq. 4.
The Macaulay bracket in eq. 4 indicates the exclusion of the compressive stresses on the fracture behaviour of the joints in the normal direction. Tensile cracking of masonry joints is governed by the defined tensile strength of masonry joints. Critical shear stress of joints prior to failure is described by Mohr-Coulomb failure eq. 5.

= +
The shear strength of masonry joints is calculated based on eq. 5, in which the cohesion, coefficient of friction and normal compressive stress are taken into consideration, thus is used to define the shear strength of masonry joints ( and ). Correspondingly, the possible pre-failure enhancement in shear behaviour due to frictional resistance is considered in the crack initiation criterion of masonry joints in the surface-based cohesive model.
In addition, the coefficient of friction of the masonry joints is defined to simulate the postfailure shear sliding behaviour (Tangential behaviour). The critical sliding shear stress ( ) is obtained based on the friction law eq. 6, which is governed by a linear relationship between the coefficient of friction and normal compressive stress.
The above friction formulation indicates the sliding of masonry units when the shear stress in the failed masonry joints is more than the critical sliding shear stress ( ).
Once the damage initiation criterion is reached, the propagation of cracks in the masonry joints causes stiffness degradation at a defined rate which leads to total strength loss and failure of joints. Thus, eq. 1 is rewritten as eq. 7: is the damage evolution variable, the value increases from 0 to 1 as per continuity of traction stresses after the damage initiation criterion met. In this study, a linear damage evolution variable is assumed by specifying the energy dissipated as a result of the damage process Fig. 3. The damage variable is expressed as: The effective separation , is given by [19]: The effective separation at complete failure , is also expressed as: is obtained from the Benzeggagh-Kenane (BK) law [20] since it is the most suitable when the critical fracture energies of both shear directions (mode ІІ and mode ІІІ) are the same, which is the case in masonry joints. The exponent, , in the BK law is set as 2 assuming brittle behaviour [20] in masonry joints. The critical fracture energy , under the mixed-mode in the BK law is expressed as:

Extended finite element method:
The extended FEM approach (XFEM) was originally developed by Belytschko and Black [21] to simulate crack propagation in an element based on the nodal displacements of the element around the crack tip and without the requirement for re-meshing. In the XFEM, discontinuous enrichment functions are added to the classical FEM based on the partition of unity concept proposed by Melenk and Babuska [22]. The enrichment functions, which make the crack independent of the mesh, are expressed as the approximation for a displacement vector function, , and is written as below [23]: In the above enrichment functions, N I (x) is associated with nodal shape functions, u I is nodal displacement vector, H(x) is associated with discontinuous jump functions to form the crack path, a I is vector of the nodal enriched degree of freedom, F ∝ (x) is associated with the cracktip functions to develop cracks at the tip and b I ∝ is the vector of the nodal enriched degree of freedom.
In this paper XFEM is used to simulate crack initiation and propagation in masonry units under tension and shear Fig. 4 without a priori-specification of a crack path and crack location. XFEM cracks are simulated based on the cohesive segment method [24]. The method can follow the traction separation law described above for the surface based cohesive behaviour to initiate and propagate cracks in the elements of masonry units. Where: The hydrostatic pressure stress, p, is given by: And the von Mises equivalent stress, q, is given by: = √

Elastic behaviour of expanded units:
The elastic modulus of the expanded masonry units must be adjusted and made to have an equivalent elastic response to the original masonry assemblage (unit and mortar). It is to be determined by taking the original masonry unit and mortar moduli of elasticity and geometry of the masonry assemblage into account. For this purpose, eq. 20 is proposed based on the assumption of a stack bond between masonry units and uniform stress distribution in masonry constituents. It is presented as:

Finite element modelling and analysis:
In The steps adopted to conduct the numerical analysis impose the actions (load or displacement) to the model either based on load control or displacement control. In both cases, the actions are incrementally imposed. Large displacement non-linear geometry effects were considered in all models. A general non-linear static procedure was adopted to follow a Newton-Raphson algorithm solution which iteratively solves for equilibrium in each increment. In addition, stiffness degradation and softening behaviour of masonry joints induce numerical instability; therefore viscous regularisation is required as a damage stabiliser in the surface-based cohesive behaviour model to simulate the full failure of masonry joints without numerical convergence difficulties. To obtain an appropriate viscosity parameter for specifying the stabilisation in masonry joints, parametric studies were conducted. In this parametric study, the viscosity parameter was adopted by considering the computational time, the effect of the parameter on the overall response of the model and the convergence of numerical analysis. For this purpose, three different viscosity parameters were tested by simulating the wall used as the validation example under in-plane load; while mesh size was kept constant (each masonry unit was modelled with 7 x 2 x 3 elements Fig. 7 (a)). The detailed comparison between numerical results using different viscosity parameters are shown in Fig. 9. Based on this study and by considering a balance between the accuracy of results and computational time, the viscosity parameter was adopted as 0.002 for the rest of the numerical analyses.

Response of masonry under in-plane loading:
The data and results of experimentally tested masonry shear walls undertaken and reported in After applying the compressive load, the vertical movement of the top beam was restrained, and then a monotonic load was horizontally applied to the walls via the top beam Fig. 10. corresponding experiments [28]. The mechanical properties and data used for defining the numerical models were obtained from the experimental results reported in [28], the data used in [4] and using eq. 2 and eq. 20. The mortar elastic modulus ( ) was calculated by considering the given values of and and using eq. 2. Then, the adjusted elastic modulus ( ) of the masonry wall was calculated using eq. 20. All parameters used in this validation study are summarised in tables 1, 2 and 3.    Axial strain (mm/mm) softening parts of axial compressive stresses versus plastic strain of the masonry assemblages are required. In addition, other required parameters such as dilation angle (ψ), friction angles (β) and flow stress ratio (R) are defined as follows; ψ was selected as a lower bound as 11.3 degrees based on [29], β is set at 36 degrees based on the coefficient of friction for masonry units and K was set as the default value equal to 1 [27]. The compressive stress-strain curves required to define the compressive behaviour of the numerical models were obtained by considering the study conducted by Kaushik et al. [30], making use of the ultimate compressive strength values (σ c ) given in the experiments and the adjusted elastic modulus value ( ) calculated from eq. 20. For detailed information, the reader is referred to [30].
The compressive stress strain curves are shown in Fig. 11 (a) and (b).  The term STATUSXFEM in Fig. 14

Response of masonry under out-of-plane loading:
The focus of this validation study is the on-plan C-shaped masonry wall tested by Griffith and Vaculik [31] under out of plane loading, see   [32]. The tensile bond strength of joints ( ) was considered as one third of the flexural strength ( ) of the masonry system as recommended in [33]. The flexural strength of the masonry system was calculated based on the virtual work method by considering the experimental results reported in [31] i.e. the ultimate horizontal load carrying capacity (3.04 kPa), the collapse mechanism (yield line) was idealised from the experimental result as shown in Fig. 17. In addition, the horizontal ( ℎ ) and diagonal ( ) bending moment capacities for the case shown in Fig. 17 were determined by using eqs. 21 and 22 reported in [34] and [35], respectively.
In the above equations, k b is a numerical factor and taken as 0.214 for stretcher bond masonry walls, and ∅ is the angle of the diagonal crack line. Fig. 17 The idealised failure mechanism considered in the virtual work method The value of tensile fracture energy is not dependant on the tensile bond strength of joints, the average value of 0.012 N/mm is recommended in the absence of detailed information by Angelillo et al. [36]. The cohesion value of the joint interfaces was considered as 1.4 of the tensile bond strength as implemented in [4]. The shear fracture energy of masonry joints was set at 0.04 N/mm. According to [4], the shear fracture energy of masonry joints with a cohesion strength (c) ranging from 0.1 to 1.8 MPa ranges from 0.01 to 0.25 N/mm. In this study, 0.04 N/mm was adopted with a view of giving best agreement between experimental and numerical results. In    Table 6 Properties for the masonry units Similarly above, the dilation angle (ψ), friction angles (β) and flow stress ratio (R) are defined as 11.3 degrees, 36 degrees and 1, respectively, for the Drucker Prager model. The compressive stress strain curve was obtained Fig. 18 based on [30], making use of the ultimate compressive strength value ( ) and masonry elastic modulus ( ) given in the experiment [31]. The numerical surface pressure versus displacement response of the wall showed a good agreement with the experimental response. It was reported in the experiment that, the post peak strength response of the wall was relatively constant from the beginning of non-linear behaviour until reaching excessive displacement. This response was reported as being due to redistribution; initially mobilising the diagonal bending moment resistance via stepped diagonal cracks, then the horizontal bending moment resistance provided by the vertical edges of the return walls. The same response was captured in the numerical model Fig. 19. The experimental mode of failure in [31] is reported as the formation of horizontal cracks at mid height of the inner wall and stepped diagonal cracks extending from the middle to the corners at the inner surface of the main wall Fig. 20(a). This scenario was also well captured in the numerical model Fig. 20(b). At this stage in the experiment, the wall was unloaded and no further damage in the wall was reported in [31],  was considered as one third of as recommended in [33]. Similar to the table 5, the tensile fracture energy of the joint interfaces was set as 0.012 N/mm as recommended in [36]. The cohesion value of the joint interface was set as 0.64 MPa as reported in the experiment [37]. The shear fracture energy was considered as 0.04 N/mm, which was within the range of values recommended in [4].   [36] c (MPa) [37] µ, [28] (N/mm) (  c (MPa) [37] µ, [37] (N/mm) (Calibrated) 0.028 0.0010 0.038 0.425 0.0038 Table 9 Non-linear material properties for the joint interface where the DPC layer presents i.e. in the bed joint between the first and second courses In addition, the experimental results showed that the failure mode was due to formation of horizontal cracking along the first bed joint containing the DPC layer and subsequent sliding of the upper part of the wall over the first bed joint. Furthermore, vertical cracks occurred in the first course of the wall as a result of the cyclic sliding movements of the upper part of the wall over the first bed joint Fig. 27 (a). Similar failure modes were observed in the numerical results Fig. 27 (b) and (c). It is also worth mentioning that the numerical compressive stresses attained were below the compressive strength of the wall during all the cycles Fig. 28, similarly in the experiment no crushing under compression was observed.

Conclusions
In this paper a combination of constitutive models has been employed together with the extended finite element method (XFEM) to simulate 3D masonry structures using a simplified Micro-modelling approach. In the new approach progressive cracking, and nonlinear post-failure behaviour between the masonry joint interfaces were well-captured by using a cohesive, surface-based approach with a traction separation law. In addition, crack plane, out-of-plane and in-plane cyclic loads, which were able to reproduce experimentally observed behaviour with accuracy and without numerical convergence difficulties.
The approach described has a variety of novel and beneficial features: for the first time masonry structural systems can be modelled in 3D when subject to arbitrary combinations of in-plane and out-of-plane loads without resort to user generated numerical code. This is a significant step forward because it allows masonry to be modelled by practitioners without the resources to develop, and validate, their own code. Moreover, the approach uses a quasistatic solution procedure rather than relying on dynamic-explicit procedures that are onerous to use and interpret. The approach is also the first to capture the cracking of masonry units without the need to specify crack locations in advance, thus removing the need to make significant assumptions, different for each possible load, that are inherent in the approaches available previously. However, the limitation of the proposed model in capturing the failure mode under out-of-plane loads is demonstrated by the fact that multiple distributed cracks rather than single cracks were predicted. This outcome may be due in part to the assumption of uniform properties for unit-mortar interfaces. Further study on the effect of the interface parameters is required. Such further study could include random assignment of varying interface properties across the wall to simulate the variability that occurs in practice.
The approach opens the possibility of analysing complete masonry structures under complex load combinations. It also offers the possibility of examining the efficacy of strengthening systems when applied to masonry structures. Such systems typically require consideration of 3D and out-of-plane behaviour for their effects to be fully captured, something that was previously not so readily possible to obtain numerically. This ability will allow the development and optimisation of strengthening systems without the need for extensive and expensive experimental programmes. This line of research is being pursued by the authors.