Next Article in Journal
Optimizing Financial Allocation for Maintenance and Rehabilitation of Munster’s Road Network Using the World Bank’s RONET Model
Previous Article in Journal
Application of a Non-Invasive Technique for the Preservation of a Fortified Masonry Tower
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Discrete Element Bonded-Block Models for Detailed Analysis of Masonry

1
National Laboratory for Civil Engineering (LNEC), Av. do Brasil 101, 1700-066 Lisbon, Portugal
2
School of Civil Engineering, University of Leeds, Woodhouse Lane, Leeds LS2 9JT, UK
*
Author to whom correspondence should be addressed.
Infrastructures 2022, 7(3), 31; https://doi.org/10.3390/infrastructures7030031
Submission received: 10 January 2022 / Revised: 15 February 2022 / Accepted: 23 February 2022 / Published: 25 February 2022
(This article belongs to the Section Infrastructures and Structural Engineering)

Abstract

:
A detailed modelling approach to represent masonry at the meso-scale is proposed, based on the discrete element method, considering the nonlinear behavior of the joints and the units. The fracture of units is represented by the bonded-block concept, in which a random network of potential cracks is created, allowing the progressive development of failure mechanisms. For simplicity, only the 2D case is presented, but the extension to 3D is straightforward. A key component of the proposed model is a framework for a joint or interface constitutive model, including the post-peak softening range, taking into account the experimental fracture energies. In this model, the softening curves in tension or shear are defined by piecewise linear segments, calibrated to reproduce the most common masonry constitutive models. The essential issues involved in the application of bonded-block models to masonry are examined, namely the block shape, either Voronoi polygons or triangles; size; deformability; and the influence of the main constitutive parameters. Uniaxial compression tests are analyzed in detail. The simulation of a well-known experiment of a brick panel under shear shows the good performance of the proposed approach. The investigation results demonstrate the model’s capabilities and provide guidelines for its application.

1. Introduction

Masonry is a heterogeneous material composed of units, such as bricks, stones, or blocks of various materials and shapes, bonded together with or without mortar. The local source of materials, the construction techniques and the age of the constructions makes for a great diversity of structural types and performance. In common, the discontinuous nature of these structures and the low tensile strength of the materials (units and binders) create a complex nonlinear behavior. It is to be expected that different conceptual models can be used to analyze the mechanical response of masonry. The purpose of the analysis is important in the selection of the most appropriate idealization. Models can be used in the design of new structures, the assessment of existing ones, the evaluation of the need for reinforcing historical constructions, and other situations. Whether the aim is the examination of the ultimate capacity or in-service performance, different tools of analysis, possibly with different underlying assumptions, may be chosen for each type of structure. Elaborate models always require more extensive material characterization, so the available experimental information always limits the complexity of the idealization. A wide class of numerical models has been applied successfully to masonry, as discussed in various critical review studies (Lourenço [1]; Roca et al. [2]; D’Altri et al. [3]).
The scale of the analysis is important in the choice of model. The global response of large constructions is usually better addressed with equivalent continuum models (“macro-models”) using the finite element method (FEM) [1]. An explicit representation of the joints, as in the designated “micro-models”, is often possible in the study of panels, walls, or smaller buildings. For large structures, the computational demands increase, but these more detailed models are increasingly applied to more complex problems. Both the finite element method (FEM) and the discrete element method (DEM) can address this type of model, in which the blocky structure of masonry is reproduced. DEM belongs to a broad class of models that view the structure as a set of blocks in mechanical interaction. Its power and versatility has been demonstrated in many applications to masonry (e.g., Sarhosis et al. [4]). It is conceptually linked to other discrete approaches, which have been devised to extend the scope of classical limit equilibrium methods, often tailored to specific problems, namely rigid block analysis, macro-blocks, FEM/DEM and related techniques (e.g., Milani et al. [5]; Caliò et al. [6]; Chiozzi et al. [7]; Baraldi et al. [8]). The DEM class also includes more elaborate numerical idealization, namely addressing block internal deformation and fracture, thus blurring the frontier between FEM and DEM. The point contact representation, a key feature to efficiently address the large displacement range, generally differentiates DEM codes (Lemos [9]). DEM models, as with all detailed numerical models, may require significant computational resources, thus limiting the size of structures that can be analyzed. The separate representation of unit and joint behavior is also more demanding in terms of input material parameters. Various alternative approaches were proposed to provide more a efficient analysis of larger structures. In the FEM field, the combination of micro- and macro-models to simulate different parts of the same structure is shown to be a very effective tool (Funari et al. [10]). The automation of procedures to identify the multiple domains of analysis is an important aspect of these multi-scale modelling techniques (Driesen et al. [11]).
At a finer scale of analysis, models of lab tests and small assemblages allow for a detailed representation of each component. These are models that are intended to investigate fundamental aspects of behavior, providing insights into the outcome of experiments. An active field of research is the study of rock fracture with bonded-particle models, a class of DEM models in which random assemblies of circular or spherical particles are used to simulate uniaxial or triaxial tests of rock specimens (Potyondy and Cundall [12]). An application of particle models to a masonry stone wall was presented by Azevedo et al. [13]. The use of polygonal blocks instead of circular particles, while increasing computational demands, provides a closer representation of the structure of geomaterials. Random networks of potential cracks are used to define the possible paths of fracture propagation through the solid material. An early example of use of random assemblies of polygonal blocks to define the potential fracturing paths in a concrete beam was presented by Lorig and Cundall [14]. Vonk [15] addressed the detailed modelling of concrete compression tests using triangular and polygonal blocks. More recently, these block models, designated as “bonded-block models (BBM)”, have been more extensively used in fracture studies of rock and other geomaterials (Garza-Cruz and Pierce [16]; Lorig et al. [17]). A random network of cracks creates blocks, shaped as triangles/tetrahedra or as Voronoi polygons/polyhedra (Herbst et al. [18]), through which the fracture process progressively develops (Yoon & Hazzard [19]).
In masonry, these detailed analyses at the meso-scale may address the fracture and breakage of the units or the mortar–unit interaction, mainly focusing on the analysis of lab experiments or structural components. Pina-Henriques and Lourenço [20] used FEM models for these detailed analyses. In the DEM context, Sarhosis and Lemos [21] showed that bonded-block models are capable of reproducing the observed behavior of mortared joints in tension and shear. The compressive failure of a brick stack was also addressed, representing both the units and the mortar by bonded-block assemblies. Chen et al. [22] modelled brick triplet tests, using Voronoi block assemblies to examine the fracture propagation in the mortar and mortar–brick interface. Baratucci et al. [23] modelled the same type of test under cyclic loads. Pulatsu et al. [24,25] employed 3D bonded-block models to study the fracture behavior of masonry prisms and wallettes. Bretas et al. [26] applied a similar concept to the safety analysis of masonry dams.
From its original development by Cundall, DEM attempted to represent the physical reality by means of assemblies of particles or blocks using the simplest constitutive assumptions for the mechanical interaction between them (Cundall [27]). In masonry applications, Mohr–Coulomb models are the most common, with the breakage of bonded contacts following a brittle failure model (e.g., [21]). It is widely recognized that a closer representation of fracture processes can be achieved with constitutive models displaying post-peak weakening, in tension or shear, in such a way that the fracture energies are taken into account. The widely used joint model proposed by Lourenço and Rots [28] includes an exponential softening regime. Joint constitutive laws with softening developed by Macorini and Izzudin [29], other authors for FEM models, and Pulatsu et al. [24,25] for DEM models will also be discussed.
In the present paper, the concept of bonded-block models for the detailed analysis of masonry at the meso-scale is examined. The fundamental assumptions for these models are presented and analyzed in a comprehensive manner, followed by a discussion of the various options available, which are tested in a set of numerical simulations. A novel general framework for joint constitutive laws incorporating post-peak weakening and the normal and shear directions is proposed. It is shown that it can reproduce the key aspects of masonry joint behavior, while also being able to represent the fracture processes in the bonded joint network inside the units. The results of a series of numerical simulations of unit failure under compression provide an assessment of the main variables in a bonded-block model, including the geometrical patterns and constitutive parameters. Several random block generations for each case are compared, to evaluate the variability of results. Numerical and computational issues, namely efficiency questions, are also critically examined. Finally, an application to a shear brick panel, previously analyzed by various authors with other numerical models is presented and the performance of the proposed approach is discussed.

2. Bonded-Block Models for Meso-Scale Analysis of Masonry

The designation “bonded-block model” has been employed in rock mechanics [16] through comparisons with the bonded-particle models, based on spherical particles [12]. In the masonry literature, “detailed micro-model” has also been applied to similar approaches, as in Sarhosis and Lemos [21], where the essential concepts of this type of meso-scale representation are presented.

2.1. Block Geometry and Mechanics

The DEM bonded-block model is composed of discrete blocks, separated by discontinuity surfaces of two types: (i) the real masonry joints and interfaces; (ii) the random network of discontinuities that define the potential cracks inside the units and the mortar. Figure 1a shows the case in which both the unit and the mortar are divided into Voronoi-shaped polygonal blocks, hereafter named “inner-blocks”, which are initially bonded. The real discontinuities in this case are the mortar–unit interfaces. The random cracks are assigned the strength of either the intact material of the unit or the mortar. A simplified version is shown in Figure 1b, where only the unit is divided into bonded inner-blocks. The brick-to-brick joint, which could be dry mortared, is represented by the corresponding joint constitutive model. This option is mainly intended to reduce the computational effort, in cases where the masonry joints have a well-known behavior.
Voronoi patterns resemble the grain structure of many rocks, as used by Lan et al. [30], in which inner-blocks correspond to different minerals; other block shapes have also been used. Figure 1c shows a triangular inner-block structure generated by a Delaunay algorithm, which corresponds to connecting the centers of adjacent Voronoi polygons. An alternative procedure is to subdivide each Voronoi polygon into radial triangles. If the Voronoi polygons are intended to mimic a grain structure of a rock material, the radial joints can be assigned different properties from the peripheral joints, allowing the distinction between inter-grain and intra-grain cracking (Gao and Stead [31]). More often, the inner-block shape is simply intended to provide a network of potential cracks with random orientations. Voronoi shapes typically lead to cracking patterns that resemble the observed tensile cracking in masonry units. Tetrahedral inner-blocks, generated by a Delaunay triangulation algorithm, have also been used [16]. The triangular or tetrahedral patterns provide a network of potential cracks with more continuity, as many planes meet at the nodes, so the shear strength is generally lower, as discussed in the following calibration examples.
The inner-blocks can be assumed to behave as rigid blocks, thus concentrating all the system deformation in the joints. A more versatile representation allows the blocks to deform, which can be achieved with an internal division into elements (or zones), which are standard uniform strain finite elements, triangular or tetrahedral. The interaction between inner-blocks follows the standard DEM representation using point contacts. Typically, 2-point contacts are present at the ends of each edge, as shown in Figure 2. However, if a finer internal element mesh is used, the edge will have more contacts, and thus a more refined representation of the distribution of the contact stresses. The contact mechanics will be addressed in the next section.

2.2. Micro- and Macro-Properties

The input material parameters that govern the mechanical behavior of inner-blocks and their interfaces will be designated as “micro-properties”. These have to be selected in such a way that the global response matches the experimental behavior of the masonry materials, which will be referred to as “macro-properties”. The deformability of the unit depends on the inner-block moduli and the stiffnesses of the joints between them. Globally, the random-shaped block system should match the experimental moduli of the unit material. Various options can be adopted for the split of the deformation between blocks and joints. If the inner-blocks are rigid, all of the deformability is given by the joints. At the other extreme, if the inner-blocks are assigned the unit Young’s modulus, the joints will require a very high stiffness, which may damage the computational efficiency. Intermediate solutions are possible, as will be analyzed below (ref. Section 4).
Similar considerations apply to the mortar material if it is also discretized into inner-blocks (Figure 1a). A simplified representation of the mortar joint as a zero-thickness interface (Figure 1b), the type of joint models used in the standard simplified micro-modelling approach, can be applied with the joint stiffnesses assigned to the experimental values for mortared or dry joints [1]. Rigid or elastic inner-blocks are the common choices. In principle, non-elastic (e.g., elasto-plastic) inner-blocks can be used, but the extra number of parameters required makes this choice too complex, losing the main advantages of the bonded-block concept.
The discretization of the unit into inner-blocks provides a potential fracture network of fictitious random fractures. The mechanical response of the unit is, to some extent, dependent on the pattern and size of the inner-blocks. In general, a calibration procedure is required to provide the properties of the various model components (inner-blocks and fictitious joints) that correctly reproduce the experimental behavior of each material or interface.
The first step is to reproduce the elastic response. For the units, the unit Young’s modulus needs to be matched, by splitting its deformability into inner-block moduli and joints. For a rigid block system, a first approximation can be obtained using a joint normal stiffness kn given by:
kn = E/d
where E is the unit Young’s modulus and d is the average spacing. For random patterns, this approximation has to be verified by a numerical test, and some calibration is usually required to meet the target value, E. The Poisson’s ratio (ν) shows a large dependence on the ratio between joint shear and normal stiffness, ks/kn, increasing with this ratio. A first approximation of the joint shear stiffness, for rigid block models, assumes a ratio equal to the ratio G/E, between shear modulus and Young’s modulus of the material:
ks = G/d
For a given block system, a numerical test can provide an approximation closer to the desired Poisson’s ratio of the unit material. When deformable blocks are employed, the total unit deformability will depend on both the joints and the inner-blocks. Different ratios can be adopted, as exemplified in Section 4.
Bonded-block models typically assume that blocks are rigid or elastic, thus the material strength is a function of the joint properties. In bonded-particle models, the calibration of strength parameters is a key step, as the sample strength is dependent on the particle sizes. For a given particle size, the contact micro-properties are chosen to match the known material macro-properties, using numerical tests following established procedures [12]. In the present bonded-block models with joint softening laws that incorporate a given fracture energy, the dependence on inner-block size can be reduced, but not totally eliminated, as the block size and contact discretization are typically not fine enough, given the computational constraints, to provide an accurate fracture propagation analysis. The size effect issue is analyzed in Section 4.
The calibration procedure can be performed by numerical simulation of elementary loading tests on a sample of the BBM (bonded-block model) with the block pattern and size that will be used in the global model. Uniaxial tension, uniaxial compression and, potentially, biaxial compression tests can be undertaken. The tensile strength can be calibrated by simulating a uniaxial tensile test. If the fracture energy is respected, the micro- and macro-properties are typically not very different, so the experimental tensile strength is often assigned to the joints between inner-blocks. The uniaxial compression test is the critical test for masonry structures, and it will be examined in detail in Section 4. For a Mohr–Coulomb type of model, this test will provide the cohesion and the friction angle of the BBM joints which reproduce the experimental uniaxial compression strength.

2.3. DEM Solution and Numerical Issues

The present bonded-block model was implemented in the discrete element code UDEC (Itasca [32]). This code follows Cundall’s DEM approach, employing a dynamic time-stepping algorithm to solve either quasi-static or dynamic problems. By introducing artificially high damping, the dynamic algorithm can be made to converge with the static solution, either an equilibrium state or a failure process, a procedure known as dynamic relaxation (Cundall [32]). The progressive updating of block and contact locations allow the analysis to continue into the large displacement range. When deformable blocks are used, the algorithm performs the time integration of the equations of motion of the nodes (or grid-points) that define the block geometry in Figure 2, which can be expressed as:
m   u ¨ + F d ( u ˙ ) + F e ( u ) + F c ( u ) = F a
where m is the nodal mass; u is the nodal displacement vector; u ˙ and u ¨ are the nodal velocity and acceleration vectors; F d ( u ˙ ) is the damping force vector, a function of the nodal velocities; F e ( u ) is the nodal forces equivalent to the internal element stresses, a function of the nodal displacements; F c ( u ) is the forces from the contacts on the adjacent edges; and F a is the external applied loads. The element and contact forces applied to the node are calculated according to standard interpolation procedures. The damping force may have several forms. Viscous damping can be adopted, as is the usual case in dynamic analysis. For quasi-static processes, however, viscous damping displays a high degree of dependence on system elastic properties, requiring adaptive algorithms to control the parameters for an efficient solution. A better performance and a smoother convergence to the solution can be achieved with Cundall’s “local damping” [33], a non-viscous energy dissipation procedure in which the magnitude of the damping force is proportional to the nodal out-of-balance force, while its direction always opposes the nodal velocity:
F d = α   | F u |   s i g n ( u ˙ )
where F u stands for the nodal unbalanced force vector and s i g n ( u ˙ ) is a vector of unit components with the sign of the nodal velocity components. The parameter α is typically assigned a value in the vicinity of 0.8, which is shown to produce an efficient convergence to the solution in most cases [32].
The equations of motion are solved by an explicit finite-difference algorithm, which is very appropriate for simulating a large displacement range, as the contact types and locations are updated as the system evolves. The main disadvantage is that it is only conditionally stable, requiring small time steps. To obtain a better performance, it is advisable to avoid large values of joints stiffness or element moduli, as will be discussed in Section 4. A very stiff block may be better represented as a rigid block, while joint stiffness should be limited to values that are physically reasonable.

3. A Framework for Joint Constitutive Models with Post-Peak Softening

3.1. Constitutive Models for Masonry Joints and Interfaces

Constitutive models for joints in masonry structures are generally based on the concept of a zero-thickness interface, in which the deformation is characterized by the difference in displacements across the interface, designated as interface or joint displacements. The interface displacements are typically expressed in terms of normal and shear components, and the constitutive model provides the corresponding normal and shear stresses, as a function of the displacements and possibly other state parameters. Various authors proposed constitutive relations to reproduce experimental data, using elasto-plastic (e.g., Lourenço and Rots [28]) or damage concepts (e.g., Gambarotta and Lagomarsino [34]). A post-peak weakening curve, typically with an exponential shape, is introduced in the normal and shear behavior, the area under the curve representing the fracture energies in the tension and shear modes. The softening rate can be scaled to match the experimental values of fracture energies for masonry materials (Lourenço [35]). Some models also introduce a nonlinear curve in the compressive branch in order to simulate the failure of the material under high compression.
These constitutive models with post-peak softening (or displacement weakening) have been implemented in FEM codes by various authors using interface elements (e.g., Macorini and Izzuddin [29]). In DEM codes, they have been used less often, as the Mohr–Coulomb model with a brittle failure after peak is the most common idealization. Resende et al. [36] implemented a joint model with linear softening in tension and shear in the code 3DEC, applied to concrete and rock fracture problems. More recently, several authors began applying softening laws in DEM models. Pulatsu et al. [24] proposed a model with linear or polynomial tension-softening curves. This model was extended by Pulatsu et al. [25] by adopting a shear softening behavior following linear and exponential laws and validated against tests on masonry brick wallets, and by Pulatsu et al. [37] to include yielding in compression, according to a linear softening law. Yuen et al. [38] proposed a damage–plasticity model for masonry joints, implemented in a discrete finite element context. Bisoffi-Sauve et al. [39] developed a constitutive model for the analysis of mixed mode fracture of mortar joints, also for a DEM code.
In the present BBM model, a general framework is proposed, which is intended to encompass various types of softening curves, extending the model developed by Resende et al. [36]. A piecewise linear approximation is employed for the softening in tension and shear, which has the ability to approximate the exponential curve or other shapes, even with a limited number of segments. It is shown in the examples below that a softening curve defined by two segments, i.e., a bilinear representation, provides an acceptable approximation for many problems. BBM models represent compressive failure through the progressive cracking and slip of the random joint network; therefore, the present constitutive model assumes elastic behavior in compression. The inclusion of compressive failure in a joint constitutive model provides a very appealing solution in a small displacement analysis. When collapse mechanisms are followed into the large displacement range, however, the progressive overlap of the block edges tends to pose numerical difficulties.

3.2. Proposed Constitutive Framework Based on a Piecewise Linear Weakening

The proposed constitutive framework is intended to cover the various types of joints present in the BBM model, including real joints represented as zero-thickness interfaces, as well as the potential cracks between inner-blocks in the units or in the mortar. The model comprises an elastic range, followed by a post-peak weakening curve, defined by a set of linear segments until the residual state is reached. Post-peak curves, either in tension or in shear, are defined by pairs of non-dimensional parameters (Figure 3): the joint displacement ratio, defined as the ratio of the displacement (in the normal or shear direction) over the peak elastic displacement and the joint strength ratio, defined as the ratio of the strength above residual level over the difference between the peak and residual strengths. The fracture energy in the post-peak stage is measured by the area under this curve. Once the peak strength is defined, the post-peak points can be scaled to achieve the desired fracture energy.

3.3. Tensile Behavior

Under normal stresses, the behavior is assumed to be linear elastic in compression, and for tensile stresses, below the peak tensile stress. The normal stress trial increment, assuming elastic response, is calculated as:
σ n = k n   u n
where u n is the normal displacement increment and k n is the joint normal stiffness. The stress increment is added to the existing normal stress and corrected if it exceeds the current tensile strength.
In tension, the peak normal elastic displacement is:
u np = t p k n
where t p is the peak tensile strength.
The normal displacement ratio is defined as a function of the peak elastic displacement as:
u ¯ n = u n u np
where u n is the total joint normal displacement.
Once the tensile stress exceeds the peak tensile strength, the softening law of Figure 3 is applied as a function of the joint normal displacement until the residual tensile strength, typically zero, is reached. The resulting evolution of the tensile normal stress as a function of the joint normal displacement ratio, shown in Figure 4, is given by:
t ( u ¯ n ) = t r + t ¯ ( u ¯ n )   ( t p t r )
The strength input parameters are the peak tensile strength t p , the residual tensile strength t r , and the non-dimensional post-peak curve t ¯ ( u ¯ n ) of Figure 3. The normalized softening curve is defined by a table of pairs ( u ¯ n , t ¯ ), where u ¯ n starts at 1 at peak, and t ¯ drops from 1 at peak to zero at the residual state. Alternatively, the fracture energy in mode I, herein designated as Gn, may be prescribed and a typical shape of post-peak curve selected, as exemplified in the following section.
An important component of the constitutive law is the unloading and reloading behavior, namely for cyclic loading (e.g., Oliveira and Lourenço [40]). However, even for monotonous loading, in a numerical simulation, there are always regions undergoing unloading, and occasional minor numerical unloading–reloading events, which have to be dealt with in a consistent manner. In the present model, after the peak unloading and reloading, the curves are assumed to go through the origin without permanent deformations, as generally assumed.

3.4. Shear Behavior

In shear, a Mohr–Coulomb failure criterion is used, defined by the peak and residual values of cohesion and friction angle. A dilation angle can also be prescribed. The behavior is elastic in shear until the peak envelope is reached. The shear stress trial increment is calculated as:
σ s = k s   u s
where u s is the shear displacement increment and k s the joint shear stiffness. The shear stress increment is added to the existing shear stress, and then corrected according to the failure criterion.
The peak shear elastic displacement for a cohesive-only joint is:
u sp = c p k s
where c p is the peak cohesive strength. The nonlinear behavior starts as the shear stress σ s reaches the peak Mohr–Coulomb envelope:
| σ s | = c p σ n   tan ( φ )
where σ n is the normal stress, assumed negative in compression, and φ is the friction angle.
The shear displacement ratio is defined as:
u ¯ s = u s u sp
where u s is the total joint shear displacement. The softening law is applied to the cohesion as a function of the joint shear displacement until the residual state is reached. The shear stress softening curve, shown in Figure 5, is given by:
c ( u ¯ s ) = c r + c ¯ ( u ¯ s )   ( c p c r )
where c r is the residual cohesion, and the non-dimensional softening law of Figure 3 is applied to the cohesion c ¯ ( u ¯ s ) .
The friction angle in masonry joints is typically assumed not to change in the post-peak range. However, if different peak and residual values of friction are used, the same softening ratio may be applied for friction and cohesion, as assumed, for example, by Lourenço and Rots [28].
The strength input parameters are the peak and residual cohesive strengths, c p and c r , and the friction angle. The normalized softening curve is defined by a table of pairs ( u ¯ s , c ¯ ), where u ¯ s takes the value of 1 at peak and c ¯ goes from 1 at peak to zero at the residual state. Alternatively, the fracture energy in mode II, herein designated as Gs, may be prescribed, and a typical shape of post-peak curve selected, as exemplified in the following section.
In shear, unloading and reloading are assumed to follow the elastic joint stiffness, with permanent deformations, as shown in Figure 5. It is possible to generalize this behavior by considering a reduction in the stiffness in the unloading and reloading branches (e.g., [34]).
An issue that still requires further research is the coupling between normal and shear weakening processes. In the present model, following Lourenço and Rots [28], a coupling between damage occurring in tension and shear is assumed; therefore, the post-peak strength ratio is taken at every step as the minimum of the current tensile and cohesive strength ratios.

3.5. Post-Peak Curves

Exponential post-peak weakening curves are often assumed to be a good approximation of the experiments in masonry joints [25,28,41]. This curve shape can be approximated by the piecewise linear framework in various manners. Four types of curves are shown in Figure 6, all with the same fracture energy as the exponential curve, which is also shown: a simple linear weakening law, two bilinear curves, and a trilinear curve. They lead to a better progressive match for the target. The two bilinear curves differ in the value where they reach the residual state. The parameters that provide the energy equality for the four cases are given in Table 1 for the normal direction. The displacement ratio values are a function of the parameter (Bn) and are equal to:
B n = 2   G n t p t r u np
which governs the displacement ratio value at which the residual state is reached for the linear softening case, where G n is the fracture energy in mode I. In the shear direction, the definitions are analogous, using cohesive strength instead of tensile strength, and mode II fracture energy. A comparison of the effects of these approximations on the response of the uniaxial compression test is presented below.
A simple comparison of the bilinear A curve with the experimental results of a shear test under 3 levels of normal stress (Van der Pluijm [42]) is shown in Figure 7, where the fracture energies follow Lourenço and Rots [28].

4. Analysis of the Influence of the Governing Parameters in Bonded-Block Models

4.1. Compressive Failure Simulations

In a bonded-block model, as discussed above, the parameters that govern the random joints between inner-blocks have to be calibrated in order to achieve the correct behavior of the material, either the masonry unit, or the mortar. Calibration tests are performed on a small sample of the numerical representation of the material, created with the inner-block patterns and a size that will be used in the complete model. An elementary stress field is applied to the numerical sample, and the micro-properties are adjusted until the desired macro-properties are obtained. This section will focus on the failure under compressive stresses. The parameters that define the bonded-block model will be examined in turn, in order to clarify their influence on the response of the numerical model. In the next section, the base model for the uniaxial compression test is presented. Then, a series of parametric studies can highlight the role of the governing parameters. First, the inner-block deformability, size and shape, followed by the micro-properties of the joints between them are discussed. In these simulations, several random generations of block geometries with similar parameters are performed; typically, four runs are used to obtain a mean value, which allows for an assessment of the variability of results. Mechanical properties, however, are assumed to be uniform in the sample. Random distributions of properties may be equally introduced, obeying the experimentally measured variations, to increase the realism of the representation (e.g., Sarhosis et al. [43]).

4.2. Base Model

The base model for the calibration tests under uniaxial compression is shown in Figure 8. The numerical sample has dimensions 130 mm × 280 mm. A Voronoi pattern of inner-blocks with an average dimension of 10 mm is generated, creating a system of about 400 blocks. Two rigid blocks are added to provide the boundary conditions.
In the base model, rigid blocks were used, a comparison with deformable block models being presented in the next section. The joint properties are listed in Table 2. The joint stiffnesses were calibrated to provide a global Young’s modulus of 10 GPa and a Poisson’s ratio of 0.2. This was verified by measuring the sample deformations at a vertical strain of 2.5 × 10−3, for which the model is essentially in the linear range.
The tensile strength was assumed to be 3.5 MPa, and the cohesion and friction angles were estimated to obtain a uniaxial compression strength in the order of 35 MPa. The cohesion was set at 10 MPa, with a null residual value, and the friction angle at 25 degrees, assuming no dilation. The post-peak curves were based on the bilinear A model of Figure 6, scaled to obtain fracture energies of 0.09 and 0.50 N/mm, in mode I and mode II. With these valuesm, the residual state is reached with joint normal displacement of 6.8 × 10−5 m in tension and a joint shear displacement of 1.3 × 10−4 m for shearing under zero normal stress.
The tests were performed by applying a vertical velocity to the loading block at the top, while keeping its horizontal displacement and rotation fixed. The base block was also fixed. The joints between the platen blocks and the inner-blocks were assigned a purely frictional behavior with a very low friction angle of 2 degrees. The applied velocity was 5 × 10−3 m/s, which provided a slow loading of the sample in a quasi-static manner. The load velocity was determined by a convergence study, being progressively reduced until no change in peak stress was observed; local damping was used. The applied stress was calculated by measuring the vertical reaction at the loading block. The reaction at the base block was also monitored and it was found to be very similar, showing that the simulation proceeded slowly, with the fracturing starting from a fairly uniform stress state in the sample.
For the model shown in Figure 8, the peak strength was 35.3 MPa at a vertical strain of 3.6 × 10−3. Figure 9 shows the applied stress versus the vertical strain. A brittle post-peak response typical of masonry units is obtained. Figure 10 shows the failure mode, showing the dominant vertical cracks, which were formed first, followed by the progressive separation of the blocks and lateral expansion.
An insight into the progression of the failure process can be gained by monitoring the increase in the damage in the contacts. If we define a damage indicator as the percentage of the loss of strength, starting at 0 at peak and reaching 1 at the residual state, Figure 11 shows the number of contacts with damage greater than or equal to 0.1, 0.5 and 0.90. We can see that small amounts of damage start at low vertical strains and increase slowly. Contacts with damage above 0.5 only start to grow at 3 × 10−3, near the peak. The curve for damage above 0.90 is only noticeable very close to the peak. At peak stress, many contacts reach the residual state (damage of 1), and the number increases rapidly after the peak, denoting a fairly brittle failure mode. Figure 12 shows the contacts with damage levels above 0.9 for two stages of loading, corresponding to vertical strains of 3.5 × 10−3, just before the peak, and 3.7 × 10−3, after the peak. The figures illustrate the process of coalescence of vertical cracks into longer fractures, leading to the axial splitting failure mode shown in Figure 10.
In order to measure the variability of the results due to the random nature of the block system, three other models were created using different seeds to generate the Voronoi system with the same block size. These three runs led to peak strengths of 34.0, 36.0 and 36.6 MPa, a mean value of 35.5 MPa, with an average deviation from the mean of 0.83 MPa (2.3%).

4.3. Block Deformability

The use of rigid blocks, as in the base model, has the advantage of simplicity and computational efficiency. However, deformable block models provide a more satisfactory representation of the strain field, also allowing the determination of internal stresses. In UDEC, deformable blocks are simulated by an internal mesh of triangular finite elements (or zones). Coarse meshes were used in this study, simply triangulating the Voronoi polygons.
The objective is to reproduce the unit deformability, so the Young’s modulus of the zones and the joint stiffnesses have to be calculated to give the desired global Young’s modulus. Three cases were considered in which the inner-block deformation accounts for 30%, 50% and 90% of the total elastic deformation. Of course, it would be possible to assign the total deformability to the inner-blocks, but this would imply a very high joint stiffness, which is inefficient in explicit algorithm codes such as UDEC. The 90% case is considered close enough. As in the rigid block case, trial joint stiffnesses was estimated and then corrected by a full simulation. Poisson’s ratios close to 0.20 were obtained for all cases. The properties used and the macro-properties obtained are listed in Table 3. The elastic properties are well-matched. Peak strengths are given for the base model and the mean value of four random generations. The compressive strength shows a small increase with the increase in block deformability, about 8% from rigid blocks to a case of 90% deformability from the blocks. However, the results show that it is possible to obtain a target value of sample strength with different splits of deformability between inner-blocks and joints.

4.4. Block Size

In the standard bonded-particle analysis, typically using circular particle and brittle failure contact models, the calibration of the contact strength for the chosen particle size is a key step, following well-documented procedures [12]. Potyondy and Ivars [44] indicate that 25 to 50 particles along the cross-section of the numerical sample are required to minimize the size effect on strength. In the present case of polygonal blocks with joint constitutive models allowing the specification of fracture energies, the dependence on block size is reduced but not eliminated. An accurate fracture propagation analysis would require a fine contact discretization along the joints, which is not possible for models with many blocks. In this study, only two contacts are used along each joint between inner-blocks.
The results of runs with 5 and 20 mm Voronoi blocks, in addition to the 10 mm case of the previous section, are summarized in Table 4, for the case of rigid blocks and 50% deformability in the blocks. The number of blocks in the three models is: 110, 400 and 1500. Therefore, there is a significant increase in model resolution and computer run times.
The results show a decrease in strength as the block size is increased, about 15% between the extreme cases, with a block size increased by a factor of 4. It may be concluded that the calibration of micro-properties is a recommended step, using the block size selected for the complete model.

4.5. Block Shape

An alternative to the Voronoi inner-block shape is the use of triangular blocks. A triangular block pattern can be obtained by splitting each Voronoi polygon into triangles, formed by adding a central node and connecting to the polygon nodes (Gao and Stead [31]). An alternative procedure is to create a Delaunay triangulation, which consists essentially in connecting the centers of adjacent Voronoi polygons. This option has been used by Garza-Cruz and Pierce [16] and Lorig et al. [17]. It was also applied in the present study.
Figure 13 shows the triangular pattern of inner-blocks and the final collapse of the model, showing two well-defined inclined fractures. The triangular pattern favors the development of more continuous fractures, where shearing is more likely to take place, unlike the Voronoi pattern, in which the fractures are composed by step segments, therefore acquiring a geometrical “roughness”, as shown in Figure 10 above. As a result, for equal joint properties, the triangular block pattern is weaker. Table 5 compares the compressive strength for Voronoi pattern with the triangular pattern. The cases of rigid blocks and blocks with 50% of the global deformability are shown, for the base model and the mean of four random generations. It can be seen that to obtain strength values in a similar range, the joint friction angle of the triangular block models must be increased, in this case from 25° to 35°.

4.6. Shear Strength Properties

Figure 14 displays the results of a few parametric studies assessing the influence of the shear strength parameters of the joints for equal values of tensile strength and fracture energies (as in the base model). It can be seen that the peak cohesive strength and friction angle change substantially the overall strength. Note that changing the cohesive strength by maintaining the same fracture energy is not realistic in practice, as it implies that the low cohesion case displays a longer post-peak curve. This is somehow reflected in the curve for the lower cohesion value in the figure, which approaches the other curves for higher friction angles.

4.7. Fracture Energy

The influence of the fracture energies on the compressive strength of the base model sample is depicted in Figure 15. Bilinear post-peak softening was used, as in the base model. Keeping the other parameters unchanged, the fracture energies in mode I (Gn) and mode II (Gs) were multiplied by factors of 2, 5 and 10, both independently and jointly. It may be noticed that the shear behavior has a dominant effect. The tensile failure of the joints is the key to the initiation of vertical cracks, but the path to total collapse is then governed by the shear properties.

4.8. Shape of Post-Peak Curve

The four types of post-peak curves in Figure 6 were trialed for the base model properties, i.e., keeping the fracture energies unchanged. Table 6 shows the resulting compressive strengths. It can be seen that the effect is relatively small, given the variability of the results already discussed. Evidently, the trilinear curve is the closest approximation to the exponential softening, so it may be considered the best choice. The bilinear curve results are, however, fairly close to those of trilinear curve. The linear softening option gives a strength 6% higher than the trilinear case, showing a larger deviation from the other cases.

4.9. Computational Parameters

The numerical simulations were performed assuming a quasi-static mechanical behavior. UDEC solves this type of problem by employing a dynamic algorithm in which damping provides the energy dissipation that prevents vibratory motion and the onset of dynamic response. Local damping and the non-viscous damping for explicit codes, was employed given its ability to provide a smooth response, often superior to adaptive viscous damping (Cundall [33]). In the compressive failure runs, the base block was fixed and the top rigid block was given a uniform vertical velocity, while preventing horizontal and rotational motion. The value of the applied velocity is an important parameter, which has to be determined by numerical trial and convergence studies. For computational efficiency, it is desirable that it be as high as possible, but it should not induce dynamical effects during loading. In the base model, a velocity of 0.005 m/s was used, while the time step, calculated by the code, was 1.3 × 10−6 s. As the time scale in these quasi-static runs is fictitious, the practical result showed that about 44,000 steps were required to apply a vertical strain increment of 0.001. The selected applied velocity was determined by running a series of trial simulations with progressively lower velocities until no change in peak strength was noticed. The values of applied stress reported herein were obtained by monitoring the vertical reactions in the top block. It is also advisable to monitor the reaction in the fixed base block and to compare the two curves. If the two reactions are different, an average value is likely a better indicator of the stress state within the model. In the runs performed, top and base reactions were practically similar.

5. Application: Brick Panel in Shear

The shear tests of brick panels performed by Vermeltfoort and Raijmakers [45], shown in Figure 16, have been used by various researchers to validate their numerical models. The simulation of these experiments was also trialed with the present Voronoi bonded-block model. The brick specimens had a length of 990 mm and height of 1000 mm and were composed of 18 brick layers, the top and bottom ones being fixed to steel beams. In the tests, the panels were first preloaded with a vertical top pressure, pv = 0.3 MPa for J4D and J5D and pv = 2.12 MPa for J7D, and a horizontal load Fh was then applied in the plane of the walls at the top edge under displacement control up to collapse.
Lourenço and Rots [28] presented the first numerical simulation of the tests employing the simplified micro-modelling approach. Each brick was modelled by elastic finite elements with a vertical joint in the middle to simulate possible tensile splitting. The masonry joints were modelled with the constitutive models proposed in their paper. Compressive failure was simulated at the joints with an elasto-plastic cap model, which was important for representing the crushing of the highly overstressed brick at the left lower corner. Macorini and Izzuddin [29] applied a 3D finite element model in which nonlinear behavior was assigned to interface elements obeying a constitutive model based on work-softening multi-surface plasticity. Different properties were assigned to the mortar joints and to the vertical interfaces that simulated brick cracking. D’Altri et al. [46] addressed the same experiments adopting a detailed micro-modeling methodology in which both units and mortar were discretized into solid finite elements. Brick–mortar interface elements were incorporated. The materials and the interfaces displayed non-elastic behavior according to plastic-damage constitutive laws. Therefore, representations of increasing complexity have been applied to this problem, providing a closer approximation of the experimental curves, but also increasing the computational demands and need for meaningful input material parameters. Pulatsu et al. [37] simulated the same tests with a DEM model, using a simplified micro-modelling approach, in which only the mortared joint and a vertical crack in the middle of the bricks were included. A joint model with softening laws was used, which also simulated the compressive failure of the units.
The present bonded-block model employed a Voronoi discretization of the units, but not the mortar. The mortar joints were represented without thickness as in the simplified micro-modeling approach. In this problem, this option was considered sufficient to account for the dominant shearing along the horizontal joint planes. The additional complexity in this model is the division of the bricks into inner-blocks in order to simulate the tensile failure in the region around the main diagonal of the panel, as well the cracking induced by the high compressive stresses in the toe block for the high vertical stress case. The inner-blocks were assumed elastic. The nonlinear behavior of the units was represented by the simple joint constitutive model with softening, as already described. The UDEC model is shown in Figure 17, composed of about 3600 Voronoi-shaped blocks, with an average size of 20 mm.
A model with deformable inner-blocks was selected, in which 90% of the unit deformation was accounted for by the blocks and the remainder by the fictitious joints inside the units. The block Young’s modulus was 18.6 GPa, and the Poisson’s ratio was 0.15. The joint material properties, listed in Table 7, follow the data provided by Vermeltfoort and Raijmakers [45], as well as additional data used in the numerical simulations by the researchers cited above [28,29]. In particular, the tensile and cohesive strengths used for the mortar joints in the tests with different vertical stresses, as well as the fracture energies, follow Macorini and Izzuddin [29]. The post-peak weakening assumed a bilinear shape, as in the bilinear A option in section above.
The horizontal load was applied by a velocity at the top blocks; while the bottom row of blocks was fixed. The curves of applied force against the horizontal displacement are compared with the experimental data in Figure 18. The test with a lower vertical stress (J4D and J5D) involves a larger rotational motion, with the opening of a diagonal crack from top right to bottom left. The numerical simulation accurately matches the experimental curve. For the test with a vertical stress of 2.12 MPa (J7D), more cracking of the units is observed, with the global force curve dropping significantly from the peak, which is reasonably attained in the numerical model. In Figure 19, the displacement contours highlight the path of cracked joints in the numerical simulation, involving the opening of vertical mortar joints and many cracks through the units, in the vicinity of the main failure path, which agrees with the behavior observed in the test.
It may be concluded that for the simulation of these experiments, the mortar joints can be adequately represented by the zero-thickness hypothesis, saving the extra computational cost of the discretization of the mortar into Voronoi blocks. This alternative is however available, and it is shown to be important in other tests, such as the compression of the brick stack analyzed in [21].

6. Conclusions

An enhanced understanding of the behavior of masonry structures requires insights into its essential governing mechanisms. The interaction between experimental work and detailed numerical modelling is the key to advancing this understanding, ultimately aiming to improve the predictive capabilities of engineering models. The bonded-block model approach is a promising tool for the analysis of masonry at the meso-scale. It shares the discrete element underlying proposition of modelling systems of geomaterials by means of complex assemblies of blocks and contacts which obey simple constitutive laws. In the masonry models presented, non-elastic behavior is confined to the discontinuities, either in real joints or potential crack networks through units and mortar, thus mirroring the discrete nature of fracture. This approach allows the progressive damage and failure processes in masonry structures to be studied. The proposed constitutive framework, based on the Mohr–Coulomb failure criterion and representing the post-peak range by a piecewise linear approximation, has the flexibility to encompass experimental data, taking into account the energy criteria of fracture mechanics.
The main options available to build bonded-block models have been examined, focusing on the simulation of the uniaxial compressive failure of masonry units. Block geometrical patterns and deformability assumptions were compared, and the influence of main joint material parameters assessed. Voronoi-shaped polygonal blocks have shown a good performance in the simulation of realistic cracking patterns for uniaxial compression tests. It was shown that while deformable blocks are more versatile and realistic, rigid blocks can also be used with reasonably good results. The comparison with the shear panel experiments shows the ability of the bonded-block approach to represent the fracture and breakage of the units without requiring excessive model refinement, and using simple constitutive laws, easily calibrated to match experimental data. In assemblages with mortared joints, a decision has to be made on whether the mortar needs to be discretized into inner-blocks, which substantially increases the computational costs. In many cases, as in the shear panel analyses presented, the Voronoi inner-blocks may be employed only in the units, while the mortared joints follow the standard interface constitutive models.
The comparison of various randomly generated assemblies using the same geometrical and mechanical properties showed that the dispersion of results is not excessive. It is always advisable to run several instances of a model for more representative results, but a small number seems acceptable for practical purposes. Some dependence of the results on the inner-block size was found, possibly linked to coarse contact discretization, which is inevitable in systems with many blocks. Therefore, there is a need to calibrate the micro-properties for the selected block pattern and size in order to obtain the desired experimental macro-properties of the masonry materials.
The block-bonded model was presented here in 2D for simplicity, but the extension to 3D is straightforward. The computational costs, however, increase substantially, thus limiting the size of the systems that can be addressed in practice. An issue deserving further investigation is the effect of block shape, namely the performance of triangular or tetrahedral inner-blocks, for which various generation procedures are available. In addition to the random geometries, the use of random distributions of material properties to achieve a more realistic simulation of specimen sets is also worthy of attention [43].

Author Contributions

Conceptualization, J.V.L. and V.S.; methodology, J.V.L. and V.S.; validation, J.V.L. and V.S.; formal analysis, J.V.L.; investigation, J.V.L. and V.S.; writing—original draft preparation, J.V.L.; writing—review and editing, V.S. All authors have read and agreed to the published version of the manuscript.

Funding

Part of this work was funded by the EPSRC project “Exploiting the resilience of masonry arch bridge infrastructure: a 3D multi-level modelling framework” (ref. EP/T001348/1). The financial contribution is very much appreciated.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Lourenço, P.B. Computations of historical masonry constructions. Prog. Struct. Eng. Mater. 2002, 4, 301–319. [Google Scholar] [CrossRef]
  2. Roca, P.; Cervera, P.M.; Gariup, G.; Pela, L. Structural Analysis of Masonry Historical Constructions—Classical and Advanced Approaches. Arch. Comput. Methods Eng. 2010, 17, 299–325. [Google Scholar] [CrossRef] [Green Version]
  3. D’Altri, A.M.; Sarhosis, V.; Milani, G.; Rots, J.; Cattari, S.; Lagomarsino, S.; Sacco, E.; Tralli, A.; Castellazzi, G.; Miranda, S. Modeling Strategies for the Computational Analysis of Unreinforced Masonry Structures: Review and Classification. Arch. Comput. Methods Eng. 2020, 27, 1153–1185. [Google Scholar] [CrossRef]
  4. Sarhosis, V.; Lemos, J.V.; Bagi, K. Discrete element modeling. In Numerical Modeling of Masonry and Historical Structures—From Theory to Application; Ghiassi, B., Milani, G., Eds.; Woodhead Publishing: Cambridge, UK, 2019; pp. 469–501. [Google Scholar]
  5. Milani, G.; Casolo, S.; Naliato, A.; Tralli, A. Seismic assessment of a medieval masonry tower in northern Italy by limit, nonlinear static, and full dynamic analyses. Int. J. Archit. Herit. 2012, 6, 489–524. [Google Scholar] [CrossRef]
  6. Caliò, I.; Marletta, M.; Pantò, B. A new discrete element model for the evaluation of the seismic behaviour of unreinforced masonry buildings. Eng. Struct. 2012, 40, 327–338. [Google Scholar] [CrossRef]
  7. Chiozzi, A.; Milani, G.; Grillanda, N.; Tralli, A. A fast and general upper-bound limit analysis approach for out-of-plane loaded masonry walls. Meccanica 2018, 53, 1875–1898. [Google Scholar] [CrossRef]
  8. Baraldi, D.; Reccia, E.; Cecchi, A. In plane loaded masonry walls: DEM and FEM/DEM models. A critical review. Meccanica 2018, 53, 1613–1628. [Google Scholar] [CrossRef] [Green Version]
  9. Lemos, J.V. Discrete element modeling of masonry structures. Int. J. Archit. Herit. 2007, 1, 190–213. [Google Scholar] [CrossRef]
  10. Funari, M.F.; Silva, L.C.; Savalle, N.; Lourenço, P.B. A concurrent micro/macro FE-model optimized with a limit analysis tool for the assessment of dry-joint masonry structures. In International Journal for Multiscale Computational Engineering; Begell House Inc.: New York, NY, USA, 2022. [Google Scholar] [CrossRef]
  11. Driesen, C.; Degée, H.; Vandoren, B. Efficient modeling of masonry failure using a multiscale domain activation approach. Comput. Struct. 2021, 251, 106543. [Google Scholar] [CrossRef]
  12. Potyondy, D.O.; Cundall, P.A. A bonded-particle model for rock. Int. J. Rock Mech. Min. Sci. 2004, 41, 1329–1364. [Google Scholar] [CrossRef]
  13. Azevedo, N.M.; Lemos, J.V.; Almeida, J.R. Discrete Element Particle Modelling of Stone Masonry. In Computational Modeling of Masonry Structures Using the Discrete Element Method; Sarhosis, V., Bagi, K., Lemos, J., Milani, G., Eds.; IGI Global: Hershey, PA, USA, 2016; pp. 146–169. [Google Scholar]
  14. Lorig, L.J.; Cundall, P.A. Modeling of Reinforced Concrete Using the Distinct Element Method. In Fracture of Concrete and Rock; Shah, S.P., Swartz, S.E., Eds.; SEM: Bethel, ME, USA, 1987; pp. 459–471. [Google Scholar]
  15. Vonk, R.A. Softening of Concrete Loaded in Compression. Ph.D. Thesis, Technische Universiteit Eindhoven, Eindhoven, The Netherlands, 1992. [Google Scholar] [CrossRef]
  16. Garza-Cruz, T.V.; Pierce, M. A 3DEC Model for Heavily Veined Massive Rock Masses. In Proceedings of the 48th US Rock Mechanics/Geomechanics Symposium, Minneapolis, MN, USA, 1–4 June 2014; American Rock Mechanics Association: Alexandria, VA, USA, 2014; pp. 14–7660. [Google Scholar]
  17. Lorig, L.; Potyondy, D.; Varun. Quantifying excavation-induced rock mass damage in large open pits. In Slope Stability 2020; Dight, P.M., Ed.; Australian Centre for Geomechanics: Perth, Australia, 2020. [Google Scholar]
  18. Herbst, M.; Konietzky, H.; Walter, K. 3D microstructural modeling. In Continuum and Distinct Element Numerical Modeling in Geo-Engineering; Hart, R., Detournay, C., Cundall, P., Eds.; Itasca: Minneapolis, MN, USA, 2008; pp. 435–441. [Google Scholar]
  19. Yoon, J.S.; Hazzard, J. Coupled Fracture Modelling with Distinct Element Methods. In Modelling Rock Fracturing Processes Theories, Methods and Applications; Shen, B., Stephansson, O., Rinne, M., Eds.; Springer: Berlin/Heidelberg, Germany, 2020; pp. 227–254. [Google Scholar] [CrossRef]
  20. Pina-Henriques, J.; Lourenço, P.B. Masonry compression: A numerical investigation at the meso-level. Eng. Comput. 2006, 23, 382–407. [Google Scholar] [CrossRef] [Green Version]
  21. Sarhosis, V.; Lemos, J.V. A detailed micro-modelling approach for the structural analysis of masonry assemblages. Comput. Struct. 2018, 206, 66–81. [Google Scholar] [CrossRef]
  22. Chen, W.; Konietzky, H.; Liu, C.; Fu, H.; Zhang, J. Prediction of Brickwork Failure Using Discrete-Element Method. J. Mater. Civ. Eng. ASCE 2018, 30, 06018012. [Google Scholar] [CrossRef]
  23. Barattucci, S.; Sarhosis, V.; Bruno, A.G.; D’Altri, A.M.; Miranda, S.; Castellazzi, G. An experimental and numerical study on masonry triplets subjected to monotonic and cyclic shear loadings. Constr. Build. Mater. 2020, 254, 119313. [Google Scholar] [CrossRef]
  24. Pulatsu, B.; Erdogmus, E.; Lourenço, P.B.; Quey, R. Simulation of uniaxial tensile behavior of quasi-brittle materials using softening contact models in DEM. Int. J. Fract. 2019, 217, 105–125. [Google Scholar] [CrossRef]
  25. Pulatsu, B.; Erdogmus, E.; Lourenço, P.B.; Lemos, J.V.; Hazzard, J. Discontinuum analysis of the fracture mechanism in masonry prisms and wallettes via discrete element method. Meccanica 2020, 55, 505–523. [Google Scholar] [CrossRef]
  26. Bretas, E.M.; Lemos, J.V.; Lourenço, P.B. A DEM based tool for the safety analysis of masonry gravity dams. Eng. Struct. 2014, 59, 248–260. [Google Scholar] [CrossRef] [Green Version]
  27. Cundall, P.A. A discontinuous future for numerical modelling in geomechanics? Proc. Inst. Civ. Eng.-Geotech. Eng. 2001, 149, 41–47. [Google Scholar] [CrossRef]
  28. Lourenço, P.B.; Rots, J.G. Multisurface interface model for analysis of masonry structures. J. Eng. Mech. 1997, 123, 660–668. [Google Scholar] [CrossRef]
  29. Macorini, L.; Izzuddin, B.A. A non-linear interface element for 3D mesoscale analysis of brick-masonry structures. Int. J. Numer. Methods Eng. 2011, 85, 1584–1608. [Google Scholar] [CrossRef] [Green Version]
  30. Lan, H.; Martin, C.D.; Hu, B. Effect of heterogeneity of brittle rock on micromechanical extensile behavior during compression loading. J. Geophys. Res. 2010, 115, B01202. [Google Scholar] [CrossRef]
  31. Gao, F.Q.; Stead, D. The application of a modified Voronoi logic to brittle fracture modelling at the laboratory and field scale. Int. J. Rock Mech. Min. Sci. 2014, 68, 1–14. [Google Scholar] [CrossRef]
  32. Itasca. UDEC—Universal Distinct Element Code: Theory and Background; Version 6.0; Itasca Consulting Group: Minneapolis, MN, USA, 2014. [Google Scholar]
  33. Cundall, P.A. Distinct element models of rock and soil structure. In Analytical and Computational Methods in Engineering Rock Mechanics; Brown, E.T., Ed.; George Allen & Unwin: Sydney, Australia, 1987; pp. 129–163. [Google Scholar]
  34. Gambarotta, L.; Lagomarsino, S. Damage models for the seismic response of brick masonry shear walls. Part I: The mortar joint model and its applications. Earthq. Eng. Struct Dyn. 1997, 26, 423–439. [Google Scholar] [CrossRef]
  35. Lourenço, P.B. Recent advances in masonry modelling: Micromodelling and homogenization. In Multiscale Modeling in Solid Mechanics—Computational Approaches; Galvanetto, U., Ferri Aliabadi, M.H., Eds.; Imperial College Press: London, UK, 2009; pp. 251–294. [Google Scholar] [CrossRef] [Green Version]
  36. Resende, R.; Lemos, J.V.; Dinis, P.B. Application of a discontinuity model with softening to the analysis of dam foundations using the discrete element method. In Numerical Modelling of Discrete Materials in Geotechnical Engineering, Civil Engineering and Earth Sciences; Konietzky, H., Ed.; CRC Press: London, UK, 2004; pp. 249–255. [Google Scholar] [CrossRef]
  37. Pulatsu, B.; Erdogmus, E.; Lourenço, P.B.; Lemos, J.V.; Tuncay, K. Simulation of the in-plane structural behavior of unreinforced masonry walls and buildings using DEM. Structures 2020, 27, 2274–2287. [Google Scholar] [CrossRef]
  38. Yuen, T.Y.P.; Deb, T.; Zhang, H.; Liu, Y. A fracture energy based damage-plasticity interfacial constitutive law for discrete finite element modelling of masonry structures. Comput. Struct. 2019, 220, 92–113. [Google Scholar] [CrossRef]
  39. Bisoffi-Sauve, M.; Morel, S.; Dubois, F. Modelling mixed mode fracture of mortar joints in masonry buildings. Eng. Struct. 2019, 182, 316–330. [Google Scholar] [CrossRef] [Green Version]
  40. Oliveira, D.V.; Lourenço, P.B. Implementation and validation of a constitutive model for the cyclic behaviour of interface elements. Comput. Struct. 2004, 82, 1451–1461. [Google Scholar] [CrossRef] [Green Version]
  41. van der Pluijm, R. Material properties of masonry and its components under tension and shear. In Proceedings of the 6th Canadian Masonry Symposium, Saskatoon, SK, Canada, 15–17 June 1992. [Google Scholar]
  42. van der Pluijm, R. Shear behaviour of bed joints. In Proceedings of the 6th North American Masonry Conference, Philadelphia, PA, USA, 6–9 June 1993. [Google Scholar]
  43. Sarhosis, V.; Forgács, T.; Lemos, J.V. Stochastic strength prediction of masonry structures: A methodological approach or a way forward? RILEM Tech. Lett. 2019, 4, 122–129. [Google Scholar] [CrossRef]
  44. Potyondy, D.O.; Ivars, D.M. Simulating spalling with a flat-jointed material. In Applied Numerical Modeling in Geomechanics—2020; Billaux, D., Hazzard, J., Nelson, M., Schöpfer, M., Eds.; Itasca: Minneapolis, MN, USA, 2020; Paper 03-01. [Google Scholar]
  45. Vermeltfoort, A.T.; Raijmakers, T.M.J. Deformation Controlled Meso Shear Tests on Masonry Piers—Part 2; Draft Report; Department of BKO, TU Eindhoven: Eindhoven, The Netherlands, 1993. [Google Scholar]
  46. D’Altri, A.M.; Miranda, S.; Castellazzi, G.; Sarhosis, V. A 3D detailed micro-model for the in-plane and out-of-plane numerical analysis of masonry panels. Comput. Struct. 2018, 206, 18–30. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Representation of masonry units and mortar joints by means of polygonal elements (a). Representation of units by Voronoi inner-blocks (b); or triangular inner-blocks (c).
Figure 1. Representation of masonry units and mortar joints by means of polygonal elements (a). Representation of units by Voronoi inner-blocks (b); or triangular inner-blocks (c).
Infrastructures 07 00031 g001
Figure 2. (a) Inner-blocks discretized into internal elements (zones); (b) contact point between blocks.
Figure 2. (a) Inner-blocks discretized into internal elements (zones); (b) contact point between blocks.
Infrastructures 07 00031 g002
Figure 3. Non-dimensional post-peak weakening behavior defined by linear segments: joint strength ratio vs. joint displacement ratio.
Figure 3. Non-dimensional post-peak weakening behavior defined by linear segments: joint strength ratio vs. joint displacement ratio.
Infrastructures 07 00031 g003
Figure 4. Tensile behavior. Joint normal strength ratio vs. joint normal displacement ratio.
Figure 4. Tensile behavior. Joint normal strength ratio vs. joint normal displacement ratio.
Infrastructures 07 00031 g004
Figure 5. Shear behavior. Joint cohesive strength ratio vs. joint shear displacement ratio.
Figure 5. Shear behavior. Joint cohesive strength ratio vs. joint shear displacement ratio.
Infrastructures 07 00031 g005
Figure 6. Approximation of exponential curve by 4 segmented curves with equal fracture energy.
Figure 6. Approximation of exponential curve by 4 segmented curves with equal fracture energy.
Infrastructures 07 00031 g006
Figure 7. Comparison of experimental shear softening curves (Van der Pluijm [42]) with the approximation by a bilinear model.
Figure 7. Comparison of experimental shear softening curves (Van der Pluijm [42]) with the approximation by a bilinear model.
Infrastructures 07 00031 g007
Figure 8. Base model for the compression test.
Figure 8. Base model for the compression test.
Infrastructures 07 00031 g008
Figure 9. Applied stress vs. vertical strain.
Figure 9. Applied stress vs. vertical strain.
Infrastructures 07 00031 g009
Figure 10. Failure model under vertical compression.
Figure 10. Failure model under vertical compression.
Infrastructures 07 00031 g010
Figure 11. Fraction of contacts with damage levels equal or above 0.1, 0.5, and 0.9 vs. vertical strain.
Figure 11. Fraction of contacts with damage levels equal or above 0.1, 0.5, and 0.9 vs. vertical strain.
Infrastructures 07 00031 g011
Figure 12. Contacts with damage above 0.9: (a) before peak, at vertical strain of 3.5 × 10−3; (b) after peak, at vertical strain of 3.7 × 10−3.
Figure 12. Contacts with damage above 0.9: (a) before peak, at vertical strain of 3.5 × 10−3; (b) after peak, at vertical strain of 3.7 × 10−3.
Infrastructures 07 00031 g012
Figure 13. (a) Model with triangular inner-blocks; (b) Failure mode.
Figure 13. (a) Model with triangular inner-blocks; (b) Failure mode.
Infrastructures 07 00031 g013
Figure 14. Compressive strength for different values of joint friction angle and joint peak cohesion (in MPa).
Figure 14. Compressive strength for different values of joint friction angle and joint peak cohesion (in MPa).
Infrastructures 07 00031 g014
Figure 15. Compressive strength for different fracture energies. Multiplier applied only to Gn, only to Gs, and to both.
Figure 15. Compressive strength for different fracture energies. Multiplier applied only to Gn, only to Gs, and to both.
Infrastructures 07 00031 g015
Figure 16. Shear test of brick panel (Vermeltfoort and Raijmakers [45]).
Figure 16. Shear test of brick panel (Vermeltfoort and Raijmakers [45]).
Infrastructures 07 00031 g016
Figure 17. Bonded-block model of shear test of brick panel.
Figure 17. Bonded-block model of shear test of brick panel.
Infrastructures 07 00031 g017
Figure 18. Shear panel model. Horizontal force vs. horizontal displacement.
Figure 18. Shear panel model. Horizontal force vs. horizontal displacement.
Infrastructures 07 00031 g018
Figure 19. Failure mode of shear panel test with vertical stress of 2.12 MPa. (a) Contours of horizontal displacement (in mm); (b) Distribution of cracked joints and inner-brick cracks.
Figure 19. Failure mode of shear panel test with vertical stress of 2.12 MPa. (a) Contours of horizontal displacement (in mm); (b) Distribution of cracked joints and inner-brick cracks.
Infrastructures 07 00031 g019
Table 1. Definition of the piecewise linear post-peak weakening curves in Figure 6.
Table 1. Definition of the piecewise linear post-peak weakening curves in Figure 6.
LinearBilinear ABilinear BTrilinear
Points u ¯ n t ¯ u ¯ n t ¯ u ¯ n t ¯ u ¯ n t ¯
1 (peak)11111111
21 + Bn01 + Bn/31/21 + Bn/21/31 + Bn/31/2
3--1 + 4Bn/301 + 3Bn/201 + 2Bn/31/4
4------1 + 5Bn/30
Table 2. Joint properties of base model.
Table 2. Joint properties of base model.
Joint Property Value
Normal stiffnesskn1326 GPa/m
Shear stiffnessks552 GPa/m
Tensile strengthtp3.5 MPa
Residual tensile strengthtr0
Peak cohesioncp10.0 MPa
Residual cohesioncr0
Friction angleϕ25°
Fracture energy mode IGn0.09 N/mm
Fracture energy mode IIGs0.50 N/mm
Table 3. Joint stiffness used in 3 cases and results sample for deformability and strength.
Table 3. Joint stiffness used in 3 cases and results sample for deformability and strength.
Input Micro-PropertiesMacro-Property Results
Normal Stiffness
(GPa/m)
Shear Stiffness
(GPa/m)
Inner-Block Young’s Modulus (GPa)Global Young’s Modulus
(GPa)
Peak Strength
(MPa)
(Base Model)
Peak Strength
(MPa)
(Average of 4 Models)
Rigid1326552-9.9735.335.5
30%189478933.39.9836.437.2
50%2652110520.09.9736.837.7
90%13,259552511.110.0338.339.7
Table 4. Compressive strength of specimen for 3 Voronoi block sizes.
Table 4. Compressive strength of specimen for 3 Voronoi block sizes.
Mean Strength
(MPa)
Average Deviation
(MPa)
Mean Strength
(MPa)
Average Deviation
(MPa)
Mean Strength
(MPa)
Average Deviation
(MPa)
Block size5 mm10 mm20 mm
Rigid38.40.8635.50.8333.21.0
50% block deformability40.01.237.71.534.71.1
Table 5. Comparison of compressive strength (MPa) for Voronoi and triangular block models.
Table 5. Comparison of compressive strength (MPa) for Voronoi and triangular block models.
Block ShapeFriction AngleRigid Blocks50% Block Deformability
Base ModelMean 4 ModelsBase ModelMean 4 Models
Voronoi25°35.335.536.837.7
Triangular25°27.427.528.528.7
Triangular35°35.334.134.535.1
Table 6. Compressive strength (MPa) for different post-peak curve shapes.
Table 6. Compressive strength (MPa) for different post-peak curve shapes.
LinearBilinear ABilinear BTrilinear
36.935.335.334.9
Table 7. Input joint properties for shear panel model.
Table 7. Input joint properties for shear panel model.
Mortar Joints
pv = 0.3 MPa
Mortar Joints
pv = 2.12 MPa
Intra-Unit Joints
Normal stiffness (GPa/m)82828350
Shear stiffness (GPa/m)36363630
Tensile strength (MPa)0.160.252.0
Peak cohesion (MPa)0.2240.3750.8
Residual cohesion (MPa)000.2
Friction angle (degrees)36.9°36.9°22°
G I (N/mm)0.0180.0180.080
G II (N/mm)0.1250.0500.030
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Lemos, J.V.; Sarhosis, V. Discrete Element Bonded-Block Models for Detailed Analysis of Masonry. Infrastructures 2022, 7, 31. https://doi.org/10.3390/infrastructures7030031

AMA Style

Lemos JV, Sarhosis V. Discrete Element Bonded-Block Models for Detailed Analysis of Masonry. Infrastructures. 2022; 7(3):31. https://doi.org/10.3390/infrastructures7030031

Chicago/Turabian Style

Lemos, José V., and Vasilis Sarhosis. 2022. "Discrete Element Bonded-Block Models for Detailed Analysis of Masonry" Infrastructures 7, no. 3: 31. https://doi.org/10.3390/infrastructures7030031

APA Style

Lemos, J. V., & Sarhosis, V. (2022). Discrete Element Bonded-Block Models for Detailed Analysis of Masonry. Infrastructures, 7(3), 31. https://doi.org/10.3390/infrastructures7030031

Article Metrics

Back to TopTop