Discrete mechanical models of concrete fracture

Engineering Fracture Mechanics 257 (2021) 108030 Available online 30 September 2021 0013-7944/© 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech Discrete mechanical models of concrete fracture John E. Bolandera,∗, Jan Eliášb, Gianluca Cusatisc, Kohei Nagaid aDepartment of Civil and Environmental Engineering, University of California, Davis, One Shields Avenue, Davis, CA 95616, USA bInstitute of Structural Mechanics, Faculty of Civil Engineering, Brno University of Technology, Veveří 331/95, Brno, 60200, Czech Republic cDepartment of Civil and Environmental Engineering, Northwestern University, 2145 Sheridan Road, Tech A134, Evanston, IL 60208-3109, USA dInstitute of Industrial Science, University of Tokyo, 4-6-1, Komaba, Meguro-ku, Tokyo, 153-8505, Japan ARTICLE INFO Keywords: Discrete model Lattice model Rigid-body–spring model Elasticity Fracture ABSTRACT Discrete models of solids have been motivated, in large part, by the discontinuous and heterogeneous nature of material structure and its breakdown under loading. The capabilities of discrete models have evolved over the past several decades, offering novel means for investigating material structure–property relationships. However, lack of understanding of both the utilities and disadvantages of discrete models limits their further development and applications. This paper reviews relevant features of discrete approaches applied to modeling the mechanical behavior of geomaterials, concrete materials in particular. The discrete models are classified according to their form and abilities to represent elastic and fracture behaviors in the presence of large-scale material heterogeneity. Discretization of the material domain plays a large role in this respect. Emphasis is placed on particle-based lattice models. The relative merits of various strategies for introducing reinforcing components, which are essential for many applications, are outlined. Recent advances are highlighted, including the use of discrete models for coupled, multi-field analysis. The merits of discrete approaches are summarized in the conclusions. 1. Introduction The mechanical behavior of solid materials has been a long-running subject of investigation. A core assumption of many theories for describing and predicting mechanical behavior is continuity of the displacement field [1]. In recent decades, a class of models that abandons the displacement continuity assumption has come into view. Those models are typically composed of either basic structural elements (e.g., truss or frame elements) that form a lattice network or discrete particles that interact at their mutual points of contact. Concrete is a composite material consisting of aggregates, matrix, interphases and pores that form a complex and random internal structure. This material mesostructure strongly affects macroscopic mechanical behavior under loading. Distributed microcracking often precedes the formation of larger cracks that coalesce within the failure zone(s) [2]. At each scale of observation, these cracks or slip mechanisms are discontinuities in the displacement field. Due to the presence of these features, and other examples of large scale heterogeneity, tensile fracture exhibits post-peak softening [3]. Modeling the localization of distributed microcracking into a finite-size fracture process zone, and the resulting quasi-brittle behavior of the material volume, is essential to many applications. This is especially important when the connectivity and openings of the crack network influence properties of interest, such as residual strength or permeability. Moreover, the heterogeneous nature of these materials brings a characteristic internal length that needs to be reflected by the model [4]. ∗Corresponding author. E-mail address: [email protected] (J.E. Bolander). https://doi.org/10.1016/j.engfracmech.2021.108030 Received 17 July 2021; Received in revised form 16 September 2021; Accepted 22 September 2021 Engineering Fracture Mechanics 257 (2021) 108030 2 J.E. Bolander et al. Table 1 Partial listing of discrete modeling approaches applied to concrete materials and structures. Modeling approach Model variants Classical lattice models Section 2.1 Particle-based lattice models Section 2.2 Rigid-body–spring model Lattice discrete particle model Discrete element method Section 2.3 It is relatively difficult to incorporate such features into models that assume continuity of the displacement field, though a number of such capable models exist, including [5–7] as recent examples. Alternatively, continuum formulations can account for crack kinematics by accommodating discontinuities within the solution space [8–10]. In contrast, the complex cracking mechanisms can be effectively simulated using discrete models. The discontinuity that appears at the particle boundaries provides a simple, natural means for crack representation. Discretization of the region of interest into lattice elements or particles can be of two forms. The discretization is either (i) non-physical or (ii) based on some physical units of the material structure (e.g., stiff aggregate particles embedded in a more compliant matrix for the case of concrete.) The non-physical form represents a pure discretization technique, such that the model response should be sufficiently independent of the discrete element size or shape, and without correspondence to material structure. The effects of material structure can be represented either phenomenologically via internal length parameter(s) within the model formulation (homogeneous models) or by spatially varying the mechanical parameters according to some projected image of material heterogeneity (mesoscale models). In contrast, physical discretizations explicitly represent relevant aspects of the material mesostructure. Such mesoscale models contain an internal length scale that is physically significant. These discretization types are referenced throughout the paper as they align with different motivations, methods and capabilities in discrete modeling. Discrete modeling approaches, of both non-physical and physical discretization types, are the focus of this review paper. Emphasis is placed on modeling mechanical behavior (particularly elasticity and fracture), since it is within that realm that discrete and continuum models greatly differ. The differences between the two model categories are less apparent when simulating scalar field behavior, such as that of heat transfer or mass transport. Although the discussions of discrete models are relevant to the broader field of geomechanics, the main emphasis is on modeling concrete materials and structures. This review paper starts by classifying discrete models according to their form, function and purpose. Some of the prominent modeling approaches are briefly described. The efficacy of these approaches depends on the methods of domain discretization, which are reviewed next. With that background, the elasticity and fracture properties of discrete models are covered, noting their merits and demerits relative to continuum modeling approaches. In particular, efforts to simulate the Poisson effect with discrete models are given special attention. Thereafter, the abilities of discrete models to simulate fragmentation and other phenomena that occur under dynamical loading are discussed. Methods for incorporating reinforcing materials (e.g., reinforcing bars or fibers) are highlighted, since most applications of concrete require reinforcement. The article then discusses the various computational strategies employed for discrete model analyses of fracture, including means for reducing computational cost and complexity. As many current developments and applications of discrete models involve coupling between multiple field quantities, the use of discrete models for multi-field analysis is reviewed. Application examples, involving actual structures, are presented. The advantages and disadvantages of discrete material descriptions are summarized in the conclusions. 2. Classification of discrete models The term ‘‘discrete model’’, which appears widely in the fracture mechanics literature, can take on different meanings [11]. We now provide a basic categorization of discrete models. Common to most modeling strategies, it is assumed that the field variables are determined at nodal points placed within the computational domain. From a general perspective, the models considered herein are based on pair-wise interactions between nodal points. These interactions are typically short-ranged, such that they can be defined by elements that interconnect neighboring nodal points. In some cases, one might want to enrich the pair-wise interaction by additional information from the neighborhood, e.g., in order to fully control the Poisson’s ratio (see Sections 4.1.1 and 4.1.2) or to account for lateral confinement. Discrete models can be classified according to the nature of these nodal/elemental interactions. Table 1 lists discrete models that have been applied to concrete materials and structures. The term lattice model will be used for methods where the nodal connectivity remains constant. This is in contrast to models that rely on contact determination to allow for nodal neighbors to change during the course of the analysis to accommodate material flow or large localized deformations (such as the DEM models discussed in Section 2.3). To maintain focus within this review paper, attention is mainly given to particle-based lattice models. 2.1. Classical lattice models Classical lattice models are composed of ordinary truss, beam, or frame elements for structural problems (or 1-D conduit elements for problems that involve scalar field quantities, such as temperature or electrical potential). Values of the field variable are determined at the nodal points and within a limited nodal-pair neighborhood (e.g., along the element axis), not continuously within the domain. Development of the lattice models originated in the field of theoretical physics, where they have been used to Engineering Fracture Mechanics 257 (2021) 108030 3 J.E. Bolander et al. Fig. 1. (a) Basic mechanical element of twoand three-dimensional rigid-body–spring models with contact dimension 𝑆and length ℎ; and (b) series construct for visco-elastic and viscous creep modeling. model the breakdown of various systems and establish scaling laws associated with such processes [12–14]. These models typically featured regular triangular [15–18] or square [19,20] geometry. In mechanics, truss elements have been found unsatisfactory due, in part, to their inability to transmit shear force; frame elements are needed to produce realistic results [21]. Whereas most frame element formulations rely on Bernoulli–Euler beam theory, geometrically nonlinear shear-deformable beams allow for more general deformation patterns [22]. With the growing power of modern computers, the necessary move to three dimensional versions has been made [23–27]. Recent review papers, devoted mainly to these classical lattice models, cover the subject in more detail [28,29]. Due to their relative simplicity, lattice models have flourished as academic in-house computer codes. One can easily modify open source software (e.g. [30]) for that purpose as well. 2.2. Particle-based lattice models Whereas particle-based lattice models share many features with classical lattice models, they differ in that each node is positioned within a geometric construct, or particle. Various models have been developed, differing mainly with respect to: •shape and other properties of the particle; •form of linkages between particles, •means for identifying linked neighbors and •the physical or non-physical meanings of the particle structure. The earliest particle-based lattice models were applied to simulating concrete fracture. Bažant et al. [31–33] placed nodes at the centers of circular aggregate particles and interconnected neighboring particles by truss elements. These early works demonstrated the simplicity, robustness and promising potential of discrete models, but they also highlighted some limitations. For example, truss connections do not enable control of Poisson’s ratio of the particle assemblage; planar analysis frameworks do not account for three-dimensionality of the material structure and fracture processes. 2.2.1. Rigid-Body-Spring Models For rigid-body–spring models (RBSM), the domain is discretized into ideally rigid units interconnected along their common boundaries. Each connection is composed of a zero-size set of springs and rigid-arm constraints that relate the spring set and nodal degrees of freedom [34], as shown in Fig. 1 where ℎ𝐼𝐽 and 𝑆𝐼𝐽 are the length and contact area associated with element 𝐼𝐽 . Equivalently, the flexibility of the spring set can be distributed over the contact surface. Domain discretization can be done by the Delaunay/Voronoi tessellation of randomly placed generator points. The nodal degrees of freedom are defined at either the generator points [35] or at the centroids of the Voronoi cells [36,37]. The rigid-body–spring elements are similar to the frame elements commonly used in classical lattice models [17]. Indeed, the RBSM and classical lattice models result in the same stiffness formulations for regular square lattices [35]. In contrast to ordinary lattice models, however, these contacts involve geometrical parameters (e.g., area and moment of inertia) derived from the tessellation geometry. Furthermore, since the spring set is typically positioned eccentrically to the element axis, the elemental axial and rotational stiffness terms are coupled. Though the RBSM are in most cases used for quasi-static analyses, there are also dynamical versions [38–40]. To account for creep deformations, the normal and tangential springs on the contact surface are replaced by the series construct shown in Fig. 1b [41,42]. Properties of the spring and rheological components can be assigned, for example, according to solidification theory [43] or code type formulations [44]. Although the original formulation is based on the assumption of small generalized displacements, recent extensions accommodate large displacements and rotations [45,46]. Engineering Fracture Mechanics 257 (2021) 108030 4 J.E. Bolander et al. 2.2.2. Lattice Discrete Particle Model The Lattice Discrete Particle Model (LDPM) utilizes the rigid-body–spring concept in conjunction with a physical discretization of the concrete mesostructure [47–50]. The concrete material is discretized by placing nodes at the volume centroids of each aggregate particle above a specified size threshold (Fig. 2d). The size distribution of aggregate particles follows an actual grading curve, as determined by sieve analysis. Elements are defined by the edges of the Delaunay simplices constructed from the nodal points. The elements represent the combined actions of the aggregates and intermediary matrix. The nature of interaction has been enriched to account for the effects of triaxial stress conditions on element constitutive behavior. In this way, the LDPM has been successful in simulating concrete behavior for the entire range of loading conditions, spanning uniaxial tension/compression to general triaxial stress states [51–54]. LDPM was also used successfully to simulate rock behavior [55,56]. It employs an explicit dynamical solver and has been applied to both quasi-static and dynamical loading cases, including those involving projectile penetration through reinforced concrete slabs [57–59]. Variations in the material properties can be accommodated using a stochastic modeling approach featuring random fields [60,61]. For some basic loading conditions, model parameters can be identified from experiments conducted at the concrete mesoscale [62]. 2.3. Discrete Element Methods The Discrete Element Method (DEM), sometimes also referred to as the Distinct Element Method, is an effective means for simulating the stability and flow of granular assemblies [63]. Unlike the lattice or RBSM with fixed inter-particle connectivity, the DEM updates the connectivity in time and accounts for the change in topology due to large displacements and rotations. Particle motion is tracked by (usually explicit) time integration of the dynamical equations of equilibrium. Contact modeling and stability of the solution scheme are central requirements of DEM [11]. The particle shapes have been mostly spherical, but other shapes (e.g., polyhedrons [64,65]) have also been implemented. Particle clusters [66,67] are also widely applied when complex grain shapes are needed. Recently, the level set DEM [68] has been formulated to accurately capture complex particle shapes. The rigid units in the system typically correspond to some real physical units of the heterogeneous material, i.e., the DEM models are mesoscale models. Applications to cohesive heterogeneous materials often focus on rock formations [69], yet a variety of other materials (e.g., nacre [70]) have also been studied. DEM has been used to simulate fracture of concrete [71–74]. In such analyses, concrete is represented by a collection of rigid particles that interact through point-wise contacts; the contact conditions are modeled by normal and tangential interactions dependent on the particle overlap and relative sliding. Particle type can be associated with the different phase fractions of the material. Parameters of the constitutive models acting at the contacts are calibrated, such that the DEM model represents the macroscopic elastic and fracture properties of concrete. Although the contact properties depend on the nature of the particle discretization scheme, realistic simulations of concrete fracture behavior can be achieved [67,75–77], including for cases of cyclic loading [78]. A noticeable advantage of DEM is the existence of several commercial or open-source software products directly optimized for particle simulations, for example YADE [79] or LIGGGHTS [80]. Due to the use of explicit time integration, the calculations can be easily parallelized. 2.4. Meshless particle methods There are various other methods that use particles as local representatives of a continuum. These particles are used to track movements in the continuum and integrate its governing equations. These meshless methods are therefore continuum methods, unlike the discrete methods considered herein. Such meshless methods are not studied in this review paper, but briefly covered in this section to provide perspective. One such method, which has been applied to concrete materials, is peridynamics. Peridynamic theory [81,82] assumes the material to be composed of infinitesimally small particles whose movements are governed by Newton’s second law. Particles within a finite distance from one another (called the material horizon) interact, in a pair-wise manner, according to a force density function. The particle interactions are (in a sense) nonlocal and tend to be of longer range, relative to the other discrete approaches considered herein. Newton’s second law is then used to track particle motion, in response to prescribed force or displacement boundary conditions. The force density function can include moment components, such that the pairwise links between particles are akin to frame elements. As peridynamic theory makes no assumptions regarding the continuity of displacement, both continuous and discontinuous deformations emerge naturally from the model formulation. Peridynamics has been used to simulate qualitatively microcracking and fracture of structural concrete [83], yet the extent of its applications to concrete has been quite limited. Fundamental concepts and workings of peridynamic formulations are critiqued elsewhere [84]. Another example of meshless particle method is Smoothed Particle Hydrodynamics (SPH), which discretizes the computational domain as a set of particles. Each particle interacts with its neighboring particles within a support domain at each time instant without requiring information about the initial connectivity of the particles. The field variables are tracked along with the particles even under large deformations. Applications to concrete materials include the simulation of fresh properties [85] and failure analyses of structural concrete [86,87]. Wu et al. [88] critique the application of SPH to concrete failure analyses and offer an alternative smoothed particle Galerkin method for the same purposes. Engineering Fracture Mechanics 257 (2021) 108030 5 J.E. Bolander et al. Fig. 2. (a) 2D sketch of artificially generated concrete mesostructure with connectivity provided by Delaunay triangulation. (b)–(d) Various tessellations of the domain into rigid bodies. 3. Domain discretization This section describes domain discretization for mechanical analysis. Different approaches may be used for discretizing scalar fields, as discussed in Section 7.1. 3.1. Nodal grid arrangement and connectivity Many of the early lattice models were based on regular (e.g., triangular) nodal patterns. It has since been recognized [21,33] that the pattern of nodal placement may strongly affect both the elastic and fracture properties of the model. Due to spatial isotropy of the elements framing into a node, regular (structured) lattices are elastically homogeneous under uniform modes of straining. However, this attribute is offset by the tendency for cracks to run preferentially along the principle directions of the structured grid. Irregular (unstructured) grids reduce bias on potential cracking directions but are, in general, not elastically homogeneous. These issues are described further in Section 4. The construction of irregular grids typically involves: (i) placement of generator points in the computational domain; and (ii) rules for defining nodal connectivity. Nodal placement may be done using a process of random sequential addition. Similar results can be achieved by using a Poisson distribution. The Delaunay tessellation of the nodal point set is an effective means for defining nodal connectivity [89]. The placement of points can be modified for physical discretization to account for size of individual particles (Fig. 2a). Furthermore, the dual Voronoi tessellation (Fig. 2b) uniquely tiles the computational domain. The Voronoi cells define the volumes and shapes associated with each nodal point. As explained later, the volume/shape information defines properties of mechanical elements and is also relevant to dynamical or potential field analyses. In principle, any other tessellation can be used, such as the power (Laguerre) tessellation ([90], Fig. 2c), which is convenient for physical discretization where the polyhedral cells represent the larger aggregates and surrounding matrix. LDPM features its own tessellation type, depicted in 2D in Fig. 2d, which provides a more effective representation of aggregate diameter and its relationship with the surrounding matrix. Nodal density can be graded within non-physical random lattices to provide higher resolution in regions of interest. Geometrical features can be explicitly represented through the strategic placement of nodes to sequentially define vertices, edges and faces, followed by random nodal filling of the interior [91]. However, the Delaunay/Voronoi discretization of complex features or domain geometries is challenging. In contrast to finite element methods, where nodes individually connect to define vertices, edges and surfaces, multiple coordinated nodes are needed to define such features with the Voronoi diagram. 3.2. Representation of material heterogeneity Concrete is often viewed as a three-phase material composed of aggregates, matrix, and aggregate–matrix interfaces. A fundamental aspect of domain discretization, and model development in general, is the degree to which such features (at the relevant length scales) are explicitly represented by the mesh design. Fig. 3 shows, in schematic form, various potential means for representing heterogeneity with a discrete approach. Whereas line-element links between nodes are shown in the figure, the concepts pertain to particle-based lattices, as well. Images of material structure can be synthetically generated or obtained via spatial mapping techniques, including computed microtomography [92,93]. Some mapping techniques can be applied in situ, such that 3D crack development can be observed under loading [94]. With respect to synthetically generated images, disks or spheres have often been used to the represent the aggregate particles [17]. A range of disk/sphere sizes can be introduced to account for the particle size distribution of the concrete. Real-shape particles can be generated using geometric models [95,96]. Discrete modeling approaches offer a variety of ways to represent the concrete meso-structure and its influences on mechanical behavior. Based on some form of image data, either physically mapped or synthetically generated, there are several possibilities for representing an aggregate particle embedded within a cement-based matrix (Fig. 3). Engineering Fracture Mechanics 257 (2021) 108030 6 J.E. Bolander et al. Fig. 3. Image data and corresponding levels of mesoscale material representation: (a) three-phase model featuring a refined discretization of the concrete– aggregate interface; (b) three-phase model based on regular or irregular arrangement of lattice elements; (c) two-phase model with nodes placed at the aggregate centroids; (d) mono-phase model in which mesh size relates to coarse aggregate size; and (e) mono-phase model acting as a pure discretization technique of a homogeneous medium (i.e., the effects of individual aggregates are absent). (a) Nodal points are pre-positioned to define the aggregate–matrix interface. This can be done to control interface thickness and the orientation of lattice elements spanning the interface. Control of orientation allows for refined constitutive modeling of interface behavior since, for example, normal and shear components of the elements can correspond to aggregate surface features. Thereafter, the aggregate interior and surrounding matrix region are populated with nodes resulting in a semi-random, three-phase model of the concrete [97–101]. (b) By mapping the lattice network or RBSM onto an image of the material structure, each lattice element can be associated with one of the material phases. Depending on the positions of the element nodes, an element is assigned the properties of the matrix, aggregate, or matrix–aggregate interface [17,26,102]. By mapping this non-physical discretization of the domain volume onto the material structure, the fracture response of these models is to a large extent regularized. Alternatively, the discrete elements can be formulated to accommodate weak and strong discontinuities that correspond to bi-phase elastic properties and cracking, respectively [103]. In this way, the spatial arrangement of the matrix–inclusion interface can be precisely represented. The lower bound of the size range of aggregate inclusions is constrained by available computational resources, since about three elements (at a minimum) should span the smallest significant features in the image [104]. Only small material volumes, relative to the size of inclusions, can be modeled with such mesoscale descriptions. (c) By placing nodes at the volume centroids of aggregate particles, and establishing element connectivity on that nodal point set, the resulting physical discretization captures mesoscale information at a coarser scale. This has several advantages, including reductions in computational expense that enable simulations of much larger volumes of concrete. The load transfer mechanisms between the main features of heterogeneity are explicitly represented. In that sense, the orientation of elements has physical significance [50]. The displacement jump between rigid units can be viewed as a lumping of the deformation of the soft matrix between stiff aggregates. The local stress variations resulting from the random geometry approximate the actual stress oscillations in the material. (d) Average element size is related to coarse aggregate size [105,106]. At this coarsest representation of the meso-scale, full-size structural components can be analyzed. The properties of structural elements effectively represent aggregated information about the mesostructure that spatially belongs to the individual elements. Alternatively, the influence of mesostructure can be represented by material parameters that fluctuate according to a random field [107]. (e) The lattice may be employed as a pure non-physical discretization technique [35], such that element size and network geometry do not correspond to material features. Whereas this allows for coarser discretizations that extend the size range of modeling capabilities, internal length scale(s) need to be included in the constitutive formulations, such that realistic, energy-conserving simulations can be made. Although this is attractive in concept, effective methods have so far been limited to basic loading cases (e.g., tension dominant loading, as described in Section 4.2.1). 3.3. Domain dimensionality The material structure and damage processes are typically three-dimensional, such that 3D models are required for quantitative evaluation of most problems. For example, the planar representations of aggregate inclusions in 2D models, and the associated interfaces, extend through the width of the material volume. Within such models, the area ratio of fractured interface and matrix material is incorrect. Furthermore, any load transfer or toughening mechanisms that might develop in the through-width direction are not present. Such inaccuracies cannot be corrected through adjustments in the fracture energy of the various phases, while keeping the planar representation of those phases. 2D models are also largely inadequate to simulate fracture due to compression loading. Engineering Fracture Mechanics 257 (2021) 108030 7 J.E. Bolander et al. 4. Mechanical behavior 4.1. Elastic behavior For the particle assemblies considered herein, elastic behavior is governed by the geometrical structure of the assembly and two material constants associated with the particle contacts. Let us denote these constants as the normal stiffness, 𝐸0, and tangential/normal stiffness ratio, 𝛼. In three dimensions, each contact between discrete units has some local orthonormal coordinate system with normal direction 𝐧and tangential directions 𝐦and 𝐥(see Fig. 1a). It is important for realistic shear behavior to have the normal vector parallel to the line connecting the nodes [50]. The particles are assumed to be rigid, such that the displacements of a given particle are represented by a single nodal point. The movements of two neighboring particles produce a displacement jump vector at their point of contact. The geometrical equation relating nodal displacements to the strains typically assumes that the strain vector is the displacement jump divided by Euclidean distance between the nodes ℎ. The normal (𝑠𝑁) and tangential (𝐬𝑇) components of stress are calculated from normal (𝑒𝑁) and tangential (𝐞𝑇) components of strain. The elastic constitutive relation is typically expressed as 𝑠𝑁=𝐸0𝑒𝑁𝐬𝑇=𝛼𝐸0𝐞𝑇(1) Note that the first equation involves scalar quantities, while the latter one features vectors in three-dimensional problems. For some modeling approaches, rotational stiffness between the particles is also defined. The stresses, integrated over the contact area, result in forces and moments acting at the nodes of the discrete model. The contact stiffness coefficients take the form 𝑘𝑁=𝑆𝐸0 ℎ𝑘𝑇=𝛼𝑆𝐸0 ℎ(2) where 𝑆and ℎare parameters associated with the cross-section area and length, respectively, of each element (i.e., of each twoparticle subassembly). It has become common to determine the area from the tessellation of the domain into polyhedral cells. Based on the elemental relations, the set of equations for static or dynamical equilibrium is assembled and solved. The equivalent macroscopic elastic behavior of discrete assemblies has been extensively studied for decades. Approximate equations relating the material parameters 𝐸0and 𝛼to elastic constants from continuum mechanics (e.g. elastic modulus and Poisson’s ratio) are known [50,108,109]. They impose a limit on the maximum macroscopic Poisson’s ratio 𝜈of the discrete system: 1/3 and 1/4 for plane stress and plane strain simplifications, respectively, and 1/4 for three dimensional models. These limitations complicate the application of discrete modeling to materials with higher Poisson’s ratio. The derivations were done under the assumption of Voigt’s hypothesis that places kinematic constraint on the system. Therefore, they provide a lower bound estimate of 𝜈(the achievable values of 𝜈are actually higher) and an upper bound estimate of elastic modulus. Liao et al. [110] developed a best fit hypothesis based on static constraint, which delivers an upper bound estimate of the Poisson’s ratio and a lower bound estimate of elastic modulus. These two estimates are similar for 𝛼around 1 for which 𝜈= 0, but significantly differ otherwise. Some remedies enabling arbitrary Poisson’s ratio are described in Section 4.1.2 but they are suitable only for non-physical discretizations. The geometry of the assembly plays an important role. To obtain an isotropic structure, there cannot be any directional bias in the orientations of the contacts [111]. On the other hand if one desires to obtain an anisotropic structure, it can be achieved by modifying the orientation distribution of model geometry [112]. It has been also shown that any deviation of the facet-normal direction 𝐧from the direction connecting individual nodes further decreases the Poisson’s ratio limits [113]. For example, such situations occur when other than Voronoi (or Power) tessellation is used in RBSM or when non-spherical particles are applied in DEM. Discrete models also produce a wall effect. The presence of boundaries affects the directional distribution in the tessellation: mechanical elements tend to be perpendicular to the boundary surface. In effect, the macroscopic elastic properties of the boundary layer differ from those of the interior [111]. With physical discretizations, however, one can argue that this difference corresponds to a real wall effect in concrete. Similarly to the homogenization of mechanical properties, one often needs to estimate the macroscopic stress tensor in a discrete assembly. The discrete system with rotational degrees of freedom asymptotically corresponds to micropolar continua [114,115] and therefore the stress tensor obtained in this way is generally non-symmetric and couple stress tensor emerges. In both elastic and inelastic regimes, the average stress and couple stress tensors in a system of volume 𝑉with 𝑘𝑓external forces 𝐟𝑖and 𝑘𝑚external moments 𝐦𝑖at locations 𝐱𝑖read [115–118] 𝝈=1 𝑉 𝑘𝑓 ∑ 𝑖=1 𝐱𝑖⊗𝐟𝑖(3) 𝝁=1 𝑉 𝑘𝑚 ∑ 𝑖=1 𝐱𝑖⊗𝐦𝑖+1 𝑉 𝑘𝑓 ∑ 𝑖=1 𝐱𝑖⊗(𝐱𝑖×𝐟𝑖)(4) The first equation is independent of the reference system while the second equation requires use of the centroid as the origin. Engineering Fracture Mechanics 257 (2021) 108030 8 J.E. Bolander et al. Fig. 4. Lattice configurations (left); relative error in field quantity (right). Fig. 5. Sub-assembly of a particle system under uniaxial compression: (a) macroscopic Poisson’s ratio fixed to 1/4 in 3D; (b) macroscopic Poisson’s ratio ranging from −1 to 1/4 in 3D depending on 𝛼; (c) all thermodynamically admissible values for Poisson’s ratio. 4.1.1. Local stress oscillations, elastic uniformity The random discrete structures generate local stress oscillations. Models based on physical discretization assume that these oscillations arise from mesoscale internal structure. Local stress peaks enable transversal tensile cracks under compressive loading, a fundamental phenomenon missing in homogeneous continuum models. This is understood as one of the major advantages of discrete modeling with physical discretization. In contrast, models based on non-physical discretization deem these oscillations as nonphysical, not controllable and spuriously affecting fracture development. Much attention has been therefore devoted to removing the oscillations, i.e. to obtain elastically homogeneous models exhibiting uniform strain at all contacts when subjected to uniform external straining. Regular lattices are elastically homogeneous by virtue of symmetry of the elements framing into a given nodal site. However, strong bias in the direction of crack propagation arises from such regularity in geometry [119]. Irregular lattices reduce the directional bias but, in general, they are not elastically uniform. Schlangen and Garboczi [119] developed the first method to remove stress oscillations by iteratively adjusting the properties of the elements, albeit with some elements exhibiting non-physical properties (e.g., negative cross-section areas and second area moments). It is worth mentioning that for the case of 𝛼= 1 (𝜈= 0), the lattice is elastically uniform when geometrical properties of the elements (i.e., area 𝑆and length ℎ) are defined by the Voronoi tessellation, even if the lattice discretization is irregular [35]. Consider a homogeneous material discretized by RBSM, which is a form of irregular lattice, as shown in Fig. 4. For the elasticity results presented in the figure, the lattice elements are defined by the Delaunay edges formed by a set of approximately 1000 nodal points. Under uniform uniaxial loading, applied to the nodes along two opposing faces, uniform strain is produced throughout the domain. This is evidenced by the computed displacements in the loading direction 𝑢𝑖that match (to machine/algorithmic precision) the theoretical values 𝑢𝑖for each node 𝑖. In the same way, the uniform strain condition is obtainable when simulating creep behavior, using the series construction (Fig. 1b) for each lineal spring in the RBSM [41]. The same Voronoi scaling produces homogeneous models for potential field (scalar) problems, as described in Section 7. For the potential flow results presented in Fig. 4, nodal potential was prescribed along opposing faces of the cube and 𝑢𝑖represents the calculated potential at node 𝑖. Furthermore, even when the displacement field is homogeneous the response of discrete models is not necessary equivalent to the corresponding continuous solution. This aspect is illustrated in Fig. 5 and it is intrinsically associated with the issue of Poisson’s ratio discussed in the next section. Fig. 5a shows a triangular subassembly of lattice elements subjected to a remote uniaxial compressive field. The inclined lattice struts are in compression and, if the lattice elements do not carry shear, equilibrium requires the horizontal strut to be in tension. Furthermore, if the constitutive equation is the one reported in Eq. (1) (with 𝛼= 0, which corresponds to no shear stress), both the axial strain and stress in the horizontal lattice element are tensile. This leads to a lateral expansion that corresponds to a fixed macroscopic Poisson’s ratio of 1∕4 in 3D. Obviously this solution, although homogeneous for a regular triangular lattice system, does not correspond to the continuous solution in which lateral deformation occurs at zero transverse Engineering Fracture Mechanics 257 (2021) 108030 9 J.E. Bolander et al. stress. The case in which the lattice elements carry shear (i.e, when 𝛼≠0in Eq. (1)) is portrayed in Fig. 5b, where 𝑔(𝛼)represents the dependence of 𝑠𝑁on 𝛼. By varying the shear stiffness, different values of stress in the horizontal lattice strut can be attained, even compression. Correspondingly, the lateral deformation can be controlled and different values of macroscopic Poisson’s ratio can be simulated within the limits of −1 (for 𝛼→∞) and 0.25 (for 𝛼= 0). However, the stress in the horizontal strut is equal to zero only for 𝛼= 1 which corresponds to a zero Poisson’s ratio. The reason for this behavior is that the constitutive laws in Eq. (1) do not allow for the deviatoric–volumetric split of the elastic energy, which is one of the fundamental assumptions of the theory of elasticity of continua. Indeed, such split can be included in the formulation of discrete models if one has 𝑠𝑁=𝐸𝑉𝜀𝑉+𝐸𝐷(𝑒𝑁−𝜀𝑉)(5) 𝐬𝑇=𝐸𝐷𝐞𝑇(6) where 𝜀𝑉is the volumetric strain calculated in a volume of material adjacent to the lattice element, 𝐸𝑣is the volumetric modulus and 𝐸𝐷is the deviatoric modulus. In this case (Fig. 5c) the exact continuum solution is obtained with no limitation on the value of the Poisson’s ratio [120]. In this case, however, the calculation of the volumetric strain requires information beyond that supplied by each two-node element, such that the fracture formulation becomes more cumbersome. It also leads to complete loss of stress oscillations and transverse tension, which are main advantages of discrete models with physical discretization. Other methods for partially or fully eliminating stress oscillations exist. They are described in Section 4.1.2, because stress uniformity is directly tied with an ability to represent the full range of Poisson’s ratio. 4.1.2. Simulating the Poisson effect As described in the previous section, the Poisson effect is typically simulated by changing the parameter 𝛼in Eq. (1). By reducing 𝛼, a discrete assembly of particles exhibits a larger macroscopic Poisson’s ratio. For lattice models composed of structural elements, this is done by reducing the shear stiffness of the elements. However, the full theoretical range of Poisson’s ratio is unattainable (as noted in Section 4.1). This is not a major limitation in the analysis of concrete materials, which exhibit a relatively low Poisson’s ratio, but it poses a problem for the analysis of some geomaterials. To address this problem and the question of elastic uniformity, remedies have been proposed to achieve the full range of thermodynamically admissible Poisson’s ratio. Although it is not often obvious, these remedies are based on separate treatments of the volumetric and deviatoric parts of the strain and stress tensors as discussed in the previous section. Furthermore, some nonlocal operator acting over several neighboring nodes is involved in obtaining these tensorial quantities. Consequently, the stress oscillations caused by the irregular geometry of the system can be fully or partially eliminated. Several of these remedies are described below. •Provided the starting lattice is elastically homogeneous for the case where 𝜈= 0, a hybrid combination of lattice and finite element technology can be used to introduce the Poisson effect [121]. The Delaunay tessellation can be used to form volumefilling simplex finite elements that provide only the Poisson effect (i.e., the other stiffness contributions of the finite elements have been removed). Whereas this approach is effective in simulating elastic behavior, the complications of continuum elements arise whether modeling fracture. •As discussed in the previous Section 4.1.1, the constitutive relation can be directly cast as a function of the volumetric and deviatoric parts of the strain tensor [120], see Eq. (6). The volumetric part of the strain tensor is computed from changes in volume of the corresponding Delaunay simplices. •The distinct lattice spring model (DLSM) [122] uses multi-body shear springs for each pair of contacting particles to overcome the limitations of ordinary lattice models in representing Poisson’s ratio. Shear strain is determined, in a rotationally invariant manner, from the displacements of neighboring particles using a least squares method. •An auxiliary stress can be introduced to achieve the Poisson effect [123,124]. Setting 𝛼= 1, Eq. (1) is supplemented by an additional vector term, which is determined via an iterative procedure (L. Khazanovich, personal communication, December 24, 2009). This method has been extended to represent transversely isotropic media [125]. Capabilities of the approach are demonstrated in Fig. 6, which involves a circular inclusion that is embedded within a homogeneous matrix under far-field uniform compression. The nodal stress values are calculated using Eq. (4). Since 𝑟=𝐿∕100 ≪ 𝐿, the solution for an infinite domain [126] can be used for comparison. When there is a significant mismatch in the elastic moduli of the different phases, the stress field local to such inclusions depends greatly on the values of Poisson’s ratio. •Similar modifications are done in Ref. [127] and Ref. [128], either by adding a lateral stress term into Eq. (1) or by modifying the strain definition to account for an average deformation of the bodies. In summary, these various remedies are counterproductive for physical discretizations because crucial information about the internal structure of the material would be lost. In contrast, the non-physical discretizations achieve the following benefits: the models are elastically homogeneous and the effects of Poisson’s ratio are precisely represented, both globally and locally. However, the local behavior no longer depends exclusively on the displacements of neighboring particle pairs, which complicates the formulation of fracture behavior. Engineering Fracture Mechanics 257 (2021) 108030 16 J.E. Bolander et al. Fig. 12. Explicit modeling of reinforcing components via frame and bond-link elements: (a) two-dimensional model showing the bond link configuration; and (b) three-dimensional model. Fig. 13. Transfer of force between matrix and embedded fiber crossing a crack. Bridging force effects: (a) lumped at the crack faces (in the form of a zero-size spring aligned with the fiber and having stiffness 𝑘𝑓, which evolves during fiber debonding and pullout); and (b) distributed along the embedded lengths as generalized nodal forces [199]. fracture of fiber-reinforced cement composites [186,189]. Since these reinforcing components possess nodal degrees of freedom, however, their use is cost prohibitive for modeling large volumes of fiber reinforced materials, where the number of fibers can be extremely large. Similar approaches have been developed for explicitly representing reinforcing components within DEM analysis frameworks [190]. 6.3. Line element discretization — embedded approach By embedding the reinforcing components within one or more discrete elements, they are not assigned degrees of freedom. Instead, the preand/or post-cracking actions of the reinforcement depend on the motion of the background mesh representing the matrix [191,192] or concrete materials [193–196]. A similar approach has been developed for simulating fiber reinforced polymers [197]. In this way, large numbers of reinforcing components can be represented enabling, for example, simulations of fracture of realistic volumes of fiber reinforced concrete. In combination with the explicit modeling approach, described in the preceding section, this embedded approach has been used to simulate the tensile fracture behavior of fiber-reinforced concrete containing steel reinforcing bars [198]. In one approach based on the rigid-body–spring concept, a fiber element is formed whenever the fiber crosses a boundary between matrix particles, as shown in Fig. 13a. A zero-size spring, directed along the fiber axis, is placed at the boundary crossing. The axial stretching/contraction of the spring depends (via rigid arms 𝑖−𝑝and 𝑗−𝑝) on the generalized displacements of the associated matrix nodes. Prior to matrix cracking, spring stiffness is determined according to elastic shear lag theory [191,200], such that axial stress in the fiber agrees well with theory. After matrix cracking, the spring forces evolve according to the nonlinear processes of fiber debonding and pullout from the matrix [191,201]. It has recently become evident that, when a fiber intersects multiple boundaries between matrix particles, lumping the fiber pullout force at the entry points from the crack faces is not objective with respect to non-physical mesh refinement [202]. Distributing the pullout force along the embedment lengths of individual fibers, as pictured in Fig. 13b, suppresses the spurious localization of damage and reduces mesh size dependency. Moreover, proper transfer of bridging force is necessary for accurate simulation of crack spacing and opening in fiber reinforced cement composites [199]. 7. Multi-field analyses Concrete behavior, including its susceptibility to cracking, is governed by various chemical and physical phenomena. On the one hand, cement hydration consumes water and produces heat as the cement paste becomes stiffer and stronger. On the other hand, changes in moisture content or temperature cause volumetric changes that, in the presence of restraint, produce stress in the concrete. In addition to applied mechanical loading, hygral or thermal effects may also contribute to crack initiation and Engineering Fracture Mechanics 257 (2021) 108030 17 J.E. Bolander et al. Fig. 14. Discretization strategies for multi-field analyses: (a) structural and flow network nodes are coincident; (b) structural and flow networks are defined by the Delaunay and Voronoi tessellations, respectively; and (c) dual-lattice construction in three dimensions. development. Concrete creep, and its interplay with fracture phenomena, also depends on the moisture and temperature conditions of the concrete. The durability mechanics of structural concrete is thus affected by a combination of mechanical, hygral and/or thermal processes, which typically commence at early ages [203,204] and evolve over the life-cycle of the structure. Therefore, for many applications, concrete fracture should be modeled using a multi-field analysis framework. This article centers on the mechanical modeling of concrete and its fracture. Within this subject area, continuum and discrete models differ greatly. As noted above, however, the strain or stress fields that drive fracture may arise from mechanical, hygral or thermal loadings. Accordingly, this section briefly reviews the capabilities of discrete models for simulating scalar field (e.g., relative humidity or temperature) behaviors, along with their potential advantages in multi-field analyses. 7.1. Discretization of the scalar potential field variable(s) The field quantity is typically transported along 1D conduit elements interconnected at the nodal sites. In the absence of cracking, it is common to use the same nodal sites for representing the mechanical (i.e., displacements) and flow-related quantities (e.g., moisture content and temperature) [41,205], as shown in Fig. 14a. By scaling the cross-section areas of the conduit elements according to the corresponding Voronoi facet areas, the lattice of conduit elements can precisely represent uniform flow conditions (Fig. 4). After crack formation, however, the conduit elements are highly skew to the crack path, which (in a conceptual sense) follows the boundaries of the discrete bodies (Fig. 14b). This misalignment of the flow elements with the crack path does not allow for realistic simulations of crack-assisted flow and its dependence on crack opening. By defining the flow network on the particle boundaries (e.g., on the edges of the Voronoi tessellation or generally on the dual grid [206]), both pre-cracking and postcracking flow can be simulated [37,207]. The cross-section areas of the flow elements scale according to the corresponding facets of the Delaunay tetrahedra, enabling precise simulations of uniform flow through the uncracked material (Fig. 4). Furthermore, as seen in Fig. 14b, flow elements are now aligned with potential cracks. The structural analysis provides values of crack opening, which govern the properties of the flow elements along the crack path [37,207,208]. This one-way coupling between fracture analyses and crack-assisted transport is relevant to a variety of durability problems of structural concrete [209], since most of the deterioration mechanisms are affected by the presence and transport of water or chemical species in solution [210]. For example, Wang and Ueda [211] used this form of one-way coupling to compute profiles of chloride concentration in cracked concrete specimens. Such analysis is difficult in homogeneous continuum models as the information about crack connectivity and opening is not directly available. Using this dual network complex, two-way couplings between fluid pressure and evolving crack networks enable simulations of hydraulic driven fracture in concrete and other geomaterials [206,212–215]. Whereas these representations of post-cracking flow or flow-driven fracture are realistic, this dual-lattice approach introduces several complicating factors. •In three dimensions, the dual lattice defined by the edges of the tessellation is computationally more expensive than its Delaunay counterpart. The numbers of dual lattice vertices and flow lattice elements are dramatically greater than those of the Delaunay vertices and structural lattice elements, respectively. •Coupling of the displacement and flow fields is complicated by the use of different nodal sets for each respective field (Figs. 14b and c). •The occurrence of essentially zero-length conduit elements can ill-condition the system of equations associated with the flow field analyses. Such zero-length edges rarely occur from the random placement of structural nodes, but some regular arrangements of nodes lead to degenerate Delaunay tessellations and zero-length Voronoi edges. The problem can be mitigated by assigning a sufficiently large length to the elements in question. In summary, crack development introduces new pathways for mass transport, which can be readily handled by discrete approaches, albeit with increased computational cost. In contrast, continuum approaches may require remeshing to accommodate crack development and the associated new pathways for mass transport. This is complicated for various reasons, including difficulties in simulating the transition from diffused to localized cracking, the need to transfer history-dependent internal variables onto the new mesh configurations and occurrence of crack-induced flow anisotropy. Engineering Fracture Mechanics 257 (2021) 108030 18 J.E. Bolander et al. 7.2. Flow field analyses By virtue of nodal site symmetry, regular lattices of conduit elements represent the condition of uniform flow through a homogeneous medium. The same capability is achieved with irregular lattices, provided the element areas are based on the corresponding Voronoi facet areas of the nodal point set (or potentially some other means of tiling the domain with the individual element volumes) [208]. This condition of flow uniformity is analogous to the condition of elastic uniformity, described in Section 4.1.1, for the case of a mechanically loaded lattice. Mass conservation calculations for flow entering/leaving each Voronoi cell (analogous to the force balance condition represented by Eq. (4)) provide nodal flux values [41], which can serve as a basis for stream line visualization. The previous discussions of physical and non-physical discretization approaches and their relationship with material structure, given in Section 3.2, also apply to these discrete models of potential flow. 7.3. Coupling schemes As noted, discrete models can represent oneor two-way couplings involving the network structure and field quantities of interest. Several forms of coupling exist. Simultaneous coupling of the field quantities can be introduced within the theoretical formulations (using, e.g., the effective stress concept of Biot). Alternatively a sequential procedure can be used, where one set of field quantities is calculated separately, assuming the other field quantities are in a steady-state condition. In many cases, results of the potential field analyses feed forward into the mechanical analyses [41,205]. Iterative solutions can account for interplay between the differing field variables within each computational cycle, which may better account for the path dependent, irreversible nature of crack propagation [215]. Coupled hydro-mechanical processes are central to several geoscience applications, where many advancements in modeling have been made [216]. 8. Practical applications Several applications of discrete models in realistic scenarios are discussed now. We begin with mechanical models of fracture and continue with multi-field applications. 8.1. Applications in structural mechanics The structural applications of mesoscale models of concrete are rather limited due to excessive computational cost when used for real-sized structural members. Exceptions include the simulations of reinforced concrete beams (up to length 4.6m) [51,217] and adhesive anchors [218] by LDPM. Bhaduri et al. [219] used LDPM for detailed analyses of the interplay between reinforcing bars and short fiber reinforcement during shear and bending failure of ultra high performance concrete beams. The mesoscale models of concrete are advantageously used to validate homogeneous continuum models [107,220], or for detailed study of the fracture processes in heterogeneous materials [106]. Interesting applications of mesoscale discrete models at the structural level are available for various forms of masonry [221–225]. The homogeneous or macroscale models, on the other hand, more readily accommodate the structural dimensions used in practice. RBSM has been used for structural concrete analyses, including high-fidelity simulations of concrete walls under reversed cyclic loading [226] or bridge piers [184]. Rasmussen and de Farias [227] and Kim et al. [228] apply discrete models for analysis of rock tunnels, accounting for elastic and strength anisotropy; Chang et al. [229] employed a lattice model to assess the buildability performance of 3D-printed concrete; Luković et al. [186] study effective of repair systems made of strain hardening cementitious composites. Pullout capacities and cracking patterns are simulated, considering the weakened zone caused by closely spaced reinforcing bars [177]. Eddy and Nagai [105] developed models that explicitly represent the complicated arrangement of reinforcing bars with mechanical anchorages in the beam–column joints. Effects of the arrangement of the reinforcement in the joint region could be analyzed, including the patterns of load transfer and crack propagation (Fig. 15). As another form of discrete model, Deng et al. [230] represent planar RC walls and other RC components as assemblages of nonlinear line (truss) elements. Truss cells, composed of vertical, horizontal and diagonal truss elements, are assembled into macro-elements that permit practical-scale analyses of RC structures. The representation of local behavior with two-node elements effectively captures shear failure and other forms of damage localization. 8.2. Multi-field applications With respect to concrete structures, and geomaterials in general, many application needs involve a coupling between mechanical behavior and other field quantities, such as temperature or fluid pressure. As previously noted, variations in those scalar field quantities often contribute to crack initiation and propagation. In turn, cracking may profoundly affect mass transport through the material domain. Such coupled analyses benefit from the accurate representation of cracking in discrete models. The most common coupling schemes involve mechanical behavior and fluid transport. Two-way couplings realize both the effect of fluid pressure on material structure and the effect of deformation and cracking on mass transport. One obvious application of such coupling schemes is hydraulic fracturing. Various discrete model types have been used to simulate hydraulic fracturing, for example particle-based lattice models [214,231,232], DEM [212,233–235] and peridynamics [236]. Engineering Fracture Mechanics 257 (2021) 108030 19 J.E. Bolander et al. Fig. 15. Discrete modeling of (a) shear failure of ultra high performance concrete beams [219] and (b) beam–column joint [105]: discretization of reinforcement and crack pattern after loading. Fig. 16. Cracking due to corrosion of rebars: (a) 2D simulation from [241] and (b) 3D model from [240] (red color indicates high stress level); (c) thermal spalling simulated by LDPM [245] (red color indicates large crack openings associated with fragments). Another largely studied topic is cracking due to expansion of steel corrosion products within structural concrete. The model is often only mechanical in nature. To represent the effects of expansion of the corrosion products, pressure is applied at the interface between the rebar and concrete [237,238]. Since the reinforcement is discretized in its circumferential direction, the expansion can be introduced differentially in the radial direction according to the corrosion conditions. Pullout tests and simple beam tests of reinforced concrete with corrosion damage were analyzed in [106,239,240] (Fig. 16b). One way coupling between transport and mechanical behavior was considered in simulations of deterioration due to corrosion in [241,242] (Fig. 16a). The cracks induced by externally applied pressure at the interface served as channels for the corrosion products in the transport model. Amalia et al. [243] considered also reduction of the pressure due to the movement of the corrosion products through the cracks. This scheme was expanded to simulate electric current flow during the corrosion process so as to simulate the crack propagation simultaneously [244]. Other coupling schemes have been developed to simulate various practical problems. Liu et al. [246] simulated drying shrinkage coupled with the transport of moisture by mesoscale lattice models. Creep and shrinkage of concrete have been simulated by coupling LDPM with the prevailing thermal and humidity fields [247,248]. A similar setup was applied in [249] to simulate concrete deterioration due to alkali-silica reaction and in [245] to simulate thermal spalling (Fig. 16c). RBSM was used to model deterioration due to frost damage in [250,251] by coupling mechanical behavior with heat and moisture transport. Chloride diffusion affected by cracking was studied in [252,253] using a lattice model with projected concrete mesostructure. Mass transport, heat generation and diffusion, and mechanical behavior were coupled within homogeneous RBSM to simulate the early-age conditions of concrete bridge decks [42]. Wang et al. [254] simulated cracking due to alkali silica reaction and delayed ettringite formation by mesoscale RBSM. 8.3. Identification of model parameters Many practical modeling applications that involve concrete fracture treat the concrete as a homogeneous continuum. For cases of tensile fracture, the parameters of interest include tensile strength, specific fracture energy, and the shape of the relationship that governs tension softening. Due to difficulties in conducting uniaxial tension tests, however, these parameters are often not measured directly, but rather calibrated through inverse analyses of flexural test results [255,256] or some other indirect means. This situation is more complex when other materials or loading configurations are considered [257]. In general, such calibration is needed on a case-by-case basis, due to the relative lack of physical bases of the model formulation. Engineering Fracture Mechanics 257 (2021) 108030 20 J.E. Bolander et al. One intention of mesoscale modeling is to provide increased physical bases to the analysis efforts. The degree to which this goal is met depends on the validity of the model and accurate assignments of the model parameters. Attempts have been made to measure parameter values through testing at the mesoscale [62,258] or at the microscale using nanoindentation techniques [259– 261]. In most cases, however, the parameter values are difficult to measure directly due to limitations in experimental capabilities or couplings between the behaviors represented by the parameters. Furthermore, the potential benefits of mesoscale modeling typically come with an increase in the number of parameters. The problem of parameter identification is further complicated when considering multiple, interacting field quantities in the analyses. For both homogeneous and mesoscale concrete models, parameter identification involves the fitting of one or more aspects of the model response to corresponding sets of experimental data. In most cases, a trial set of parameter values is updated by an iterative procedure. Although this updating process has been done heuristically (i.e., as a trial-and-error process), modern needs strongly favor classical algorithmic [262–264] or machine learning [265–267] approaches for inverse analyses. Janouchová et al. [268] present a probabilistic calibration of the LDPM, which uses Bayesian inference to solve the inverse problem in a robust way, while providing quantitative information regarding the uncertainties of the parameter values. Calibration and validation, at all relevant scales covered by the model formulation, are essential for both scientific and practical applications of the model. 9. Conclusions This paper reviews discrete mechanical and multi-physical models of concrete with particular focus on simulating fracture. These discrete models are described according to their purpose, formulation, and form of domain discretization. In particular, physical and non-physical forms of discretization offer differing, complimentary means for simulating material and structural behavior. In the former case, relevant aspects of material heterogeneity are directly captured by the discrete particles and their modes of interaction. In the latter case, the mechanical behavior of homogeneous media can be simulated, which can be used for structural scale modeling or the explicit representation of material phases at the concrete mesoscale. Both discrete and continuum models are capable of matching a range of experimental results as recently demonstrated, e.g., in Ref. [269]. There are however fundamental differences that should be acknowledged and understood. This paper aims to provide such understanding and clarify the merits and demerits of discrete models relative to continuum models in general, as well as mutually across the various types of discrete models. The following statements can be made. •Most discrete models are composed of two-node elements or pairwise assemblies of neighboring particles. As a consequence, stress and strain are defined as vectorial quantities that (for physical forms of discretization) tend to correspond directly with material structure. The vectorial basis of the constitutive relations provides simple, transparent representations of material behavior. •The pairwise interactions between particles, or lattice nodes, offer effective means for representing discontinuities in the displacement field due to cracking or slippage. Furthermore, since only two nodes are involved in modeling the local process of material separation (i.e., fracture), the stress-locking problems of ordinary continuum formulations are largely avoided. •By capturing the displacement discontinuity of cracks, and the stress variations associated with crack opening, discrete models are effective in representing the anisotropy induced by crack formation. When heterogeneity is present, discrete models are adept at simulating the development of diffuse cracking and its transition to localized fracture. •On the other hand, the vectorial descriptions of material behavior pose challenges in simulating elasticity of the composite material or its constituents. –The discrete structure exhibits stress variations that (for the case of non-physical discretizations) are viewed as spurious heterogeneity. When using physical discretizations, however, these variations are considered to be beneficial, as they reflect the presence of actual heterogeneity. Under uniform compressive load, for example, patches of transverse tension appear throughout the domain, which can induce splitting cracks that align with the direction of compression. –There are upper limits on the achievable Poisson’s ratio of discrete assemblies. Several remedies exist for allowing the full range of thermodynamically admissible Poisson’s ratio (see Section 4.1.2). These remedies also remove the aforementioned stress oscillations caused by the discrete structure. –Through systematic scaling of the cross-sectional areas of elements, or the contact areas between particles, discrete models can be rendered elastically homogeneous under uniform modes of straining. This is advantageous when employing non-physical forms of discretization to model homogeneous behavior, either at the structural or material component levels. –The vectorial constitutive law needs to be equipped with volumetric deformation dependence to simulate triaxial compression. This requires use of information from the element neighborhood, beyond what is available from the individual two-node elements. •It is widely accepted that irregular, or quasi-random, nodal arrangements should be used since the resulting models exert less bias on potential crack directions. As a consequence, the model response includes random features, as well. When physical discretizations are employed, these features may be associated with randomness of the natural material; on the other hand, the features appear as spurious artifacts within non-physical discretizations of the material domain. •The discretization technique unavoidably introduces a boundary layer with different material properties than the bulk material [111]. It may be interpreted as realistic for physical discretizations, but may introduce unwanted bias within non-physical ones. Engineering Fracture Mechanics 257 (2021) 108030 21 J.E. Bolander et al. •Discrete models are often the method of choice for solving coupled multi-physical problems relevant to concrete materials. Several points can be made. –Discrete models are applicable to scalar field analyses, in which the field quantities are transported via lineal conduit elements. The observations made regarding physical and non-physical forms of discretization hold for scalar field analyses, as well. Through appropriate scaling of the cross-section areas of the conduit elements, the model can be rendered homogeneous, such that uniform potential flow can be simulated without mesh bias. –The explicit description of cracking makes discrete models highly suitable for simulating transport through cracked media. –In the absence of cracking, the same set of nodal points can be used for the mechanical and scalar field analyses. In the presence of cracking, however, the dual-lattice approach (Section 7) offers a direct and robust technique for such types of coupling. The conduit elements reside in the potential planes of cracking, allowing for transport both along individual cracks and within the uncracked bulk material, as the fracture network evolves under loading. •Discrete modeling approaches accommodate the inclusion of reinforcing components, which are essential for many applications of concrete materials. Several possible approaches, differing in computational demands and accuracy, are available (Section 6). Both steel bar and short fiber reinforcement can be introduced irrespective of the background discretization of the material volume. •There are many practical applications of discrete models at the structural level (Section 8). Computational complexity is of paramount importance in applied numerical modeling. Models with physical discretization are arguably the most effective kinematic representation of mesoscale material structure. Yet, practical applications use mostly scalable non-physical discretization as the mesoscale models are still expensive, even though several reduction techniques are available (Section 5.2). The discrete models reviewed herein, and concrete models in general, rely on input parameters that are often difficult to measure and quantify. Multiscale modeling addresses some of the issues associated with parameter identification, but difficulties may persist. This calls for continued collaboration with experimentalists and recognition of the potential capabilities of inverse analysis techniques. As a related matter, verification and validation exercises, extending across the length and time scales of the model, should be a focal point of future research and model development. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgment Jan Eliáš acknowledges financial support by the project International Mobility of Researchers of Brno University of Technology, Czechia (EF16_027/0008371). References [1] Malvern LE. Introduction to the mechanics of a continuous medium. Englewood Cliffs, N.J., USA: Prentice Hall; 1969. [2] Cedolin L, Dei Poli S, Iori I. Experimental determination of the fracture process zone in concrete. Cem Concr Res 1983;13(4):557–67. http://dx.doi.org/ 10.1016/0008-8846(83)90015-7. [3] Hillerborg A, Modéer M, Petersson P-E. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 1976;6(6):773–81. http://dx.doi.org/10.1016/0008-8846(76)90007-7. [4] Bažant ZP, Pijaudier-Cabot G. Measurement of characteristic length of nonlocal continuum. J Eng Mech 1989;115(4):755–67. http://dx.doi.org/10.1061/ (ASCE)0733-9399(1989)115:4(755). [5] Negi A, Kumar S. Localizing gradient damage model with smoothed stress based anisotropic nonlocal interactions. Eng Fract Mech 2019;214:21–39. http://dx.doi.org/10.1016/j.engfracmech.2019.04.011. [6] Xue J, Kirane K. Strength size effect and post-peak softening in textile composites analyzed by cohesive zone and crack band models. Eng Fract Mech 2019;212:106–22. http://dx.doi.org/10.1016/j.engfracmech.2019.03.025. [7] Yang L, Yang Y, Zheng H. A phase field numerical manifold method for crack propagation in quasi-brittle materials. Eng Fract Mech 2021;241:107427. http://dx.doi.org/10.1016/j.engfracmech.2020.107427. [8] Oliver J, Huespe AE, Pulido MDG, Chaves E. From continuum mechanics to fracture mechanics: the strong discontinuity approach. Eng Fract Mech 2002;69(2):113–36. http://dx.doi.org/10.1016/SS013-7944(01)00060-1. [9] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A geometrically non-linear three-dimensional cohesive crack method for reinforced concrete structures. Eng Fract Mech 2008;75:4740–58. http://dx.doi.org/10.1016/j.engfracmech.2008.06.019. [10] Sukumar N, Dolbow JE, Moës N. Extended finite element method in computational fracture mechanics: a retrospective examination. Int J Fract 2015;196:189–206. http://dx.doi.org/10.1007/s10704-015-0064-8. [11] Bićanić N. Discrete element methods. In: Encyclopedia of computational mechanics, vol. 1, John Wiley & Sons, Ltd.; 2007, p. 1–33. [12] Herrmann HJ, Roux S. Statistical Models for the Fracture of Disordered Media. Elsevier; 1990, p. 353. http://dx.doi.org/10.1016/C2009-0-14278-2. [13] Herrmann HJ, Hansen A, Roux S. Fracture of disordered, elastic lattices in two dimensions. Phys Rev B 1989;39(1):637–48. http://dx.doi.org/10.1103/ PhysRevB.39.637. [14] Roux S, Guyon EG. Mechanical percolation: a small beam lattice study. J Phys Lett 1985;46:999–1004. http://dx.doi.org/10.1051/jphyslet: 019850046021099900. Engineering Fracture Mechanics 257 (2021) 108030 22 J.E. Bolander et al. [15] Bolander JE, Shiraishi T, Isogawa Y. An adaptive procedure for fracture simulation in extensive lattice networks. Eng Fract Mech 1996;53(3):325–34. http://dx.doi.org/10.1016/0013-7944(95)00200-6. [16] Moukarzel C, Herrmann HJ. A vectorizable random lattice. J Stat Phys 1992;68:911–23. http://dx.doi.org/10.1007/BF01048880. [17] Schlangen E, van Mier JGM. Experimental and numerical analysis of micromechanisms of fracture of cement-based composites. Cem Concr Comp 1992;14(2):105–18. http://dx.doi.org/10.1016/0958-9465(92)90004-F. [18] van Mier JGM, van Vliet MRA, Wang TK. Fracture mechanisms in particle composites: statistical aspects in lattice type analysis. Mech Mater 2002;34(11):705–24. http://dx.doi.org/10.1016/S0167-6636(02)00170-9. [19] Ince R, Arslan A, Karihaloo BL. Lattice modelling of size effect in concrete strength. Eng Fract Mech 2003;70(16):2307–20. http://dx.doi.org/10.1016/ S0013-7944(02)00219-9. [20] Ostoja-Starzewski M, Sheng PY, Jasiuk IM. Damage patterns and constitutive response of random matrix-inclusion composites. Eng Fract Mech 1997;58(5):581–606. http://dx.doi.org/10.1016/S0013-7944(97)00046-5. [21] Schlangen E, Garboczi EJ. Fracture simulations of concrete using lattice models: Computational aspects. Eng Fract Mech 1997;57(2–3):319–32. http: //dx.doi.org/10.1016/S0013-7944(97)00010-6. [22] Ibrahimbegovic A, Delaplace A. Microscale and mesoscale discrete models for dynamic fracture of structures built of brittle material. Comput Struct 2003;81(12):1255–65. http://dx.doi.org/10.1016/S0045-7949(03)00040-3. [23] Jivkov AP, Engelberg DL, Stein R, Petkovski M. Pore space and brittle damage evolution in concrete. Eng Fract Mech 2013;110:378–95. http: //dx.doi.org/10.1016/j.engfracmech.2013.05.007. [24] Luković M, Šavija B, Schlangen E, Ye G, van Breugel K. A 3D lattice modelling study of drying shrinkage damage in concrete repair systems. Materials 2016;9:575. http://dx.doi.org/10.3390/ma9070575. [25] Man H-K, van Mier JGM. Influence of particle density on 3D size effects in the fracture of (numerical) concrete. Mech Mater 2008;40(6):470–86. http://dx.doi.org/10.1016/j.mechmat.2007.11.003. [26] Man H-K, van Mier JGM. Damage distribution and size effect in numerical concrete from lattice analyses. Cem Concr Comp 2011;33(9):867–80. http://dx.doi.org/10.1016/j.cemconcomp.2011.01.008. [27] Mungule M, Raghuprasad BK, Muralidhara S. Fracture studies on 3D geometrically similar beams. Eng Fract Mech 2013;98:407–22. http://dx.doi.org/ 10.1016/j.engfracmech.2012.11.012. [28] Nikolić M, Karavelić E, Ibrahimbegovic A, Miščević P. Lattice element models and their peculiarities. Arch Comput Methods Eng 2018;25:757–84. http://dx.doi.org/10.1007/s11831-017-9210-y. [29] Pan Z, Ma R, Wang D, Chen A. A review of lattice type model in fracture mechanics: theory, applications, and perspectives. Eng Fract Mech 2018;190:382–409. http://dx.doi.org/10.1016/j.engfracmech.2017.12.037. [30] Patzák B, Rypl D. Object-oriented, parallel finite element framework with dynamic load balancing. Adv Eng Softw 2012;47(1):35–50. http://dx.doi.org/ 10.1016/j.advengsoft.2011.12.008. [31] Bažant ZP, Tabbara MR, Kazemi MT, Pijaudier-Cabot G. Random particle model for fracture of aggregate or fiber composites. J Eng Mech 1990;116(8):1686–705. http://dx.doi.org/10.1061/(ASCE)0733-9399(1990)116:8(1686). [32] Jirásek M, Bažant ZP. Macroscopic fracture characteristics of random particle systems. Int J Fract 1994;69:201–28. http://dx.doi.org/10.1007/BF00034763. [33] Jirásek M, Bažant ZP. Particle model for quasibrittle fracture and application to sea ice. J Eng Mech 1995;121(9):1016–25. http://dx.doi.org/10.1061/ (ASCE)0733-9399(1995)121:9(1016). [34] Kawai T. New discrete models and their application to seismic response analysis of structures. Nucl Eng Des 1978;48(1):207–29. http://dx.doi.org/10. 1016/0029-5493(78)90217-0. [35] Bolander JE, Saito S. Fracture analyses using spring networks with random geometry. Eng Fract Mech 1998;61(5–6):569–91. http://dx.doi.org/10.1016/ S0013-7944(98)00069-1. [36] Nagai K, Sato Y, Ueda T. Mesoscopic simulation of failure of mortar and concrete by 3D RBSM. J Adv Concr Technol 2005;3(3):385–402. http: //dx.doi.org/10.3151/jact.3.385. [37] Nakamura H, Worapong S, Yashiro R, Kunieda M. Time-dependent structural analysis considering mass transfer to evaluate deterioration processes of RC structures. J Adv Concr Technol 2006;4(1):147–58. http://dx.doi.org/10.3151/jact.4.147. [38] Hwang YK, Bolander JE, Lim YM. Simulation of concrete tensile failure under high loading rates using three-dimensional irregular lattice models. Mech Mater 2016;101:136–46. http://dx.doi.org/10.1016/j.mechmat.2016.08.002. [39] Hwang YK, Bolander JE, Lim YM. Evaluation of dynamic tensile strength of concrete using lattice-based simulations of spalling tests. Int J Fract 2020. http://dx.doi.org/10.1007/s10704-020-00422-w. [40] Le J-L, Eliáš J, Gorgogianni A, Vievering J, Květoň J. Rate-dependent scaling of dynamic tensile strength of quasibrittle structures. J Appl Mech 2018;85(2):021003. http://dx.doi.org/10.1115/1.4038496. [41] Bolander JE, Berton S. Simulation of shrinkage induced cracking in cement composite overlays. Cem Concr Comp 2004;26(7):861–71. http://dx.doi.org/ 10.1016/j.cemconcomp.2003.04.001. [42] Pan Y, Prado A, Porras R, Hafez OM, Bolander JE. Lattice modeling of early-age behavior of structural concrete. Materials 2017;10:231. http: //dx.doi.org/10.3390/ma10030231. [43] Di Luzio G, Cusatis G. Solidification-microprestress-microplane (SMM) theory for concrete at early age: theory, validation and application. Int J Solids Struct 2013;50:957–75. http://dx.doi.org/10.1016/j.ijsolstr.2012.11.022. [44] Bažant ZP, Panula L, Kim J-K, Xi Y. Improved prediction model for time-dependent deformations of concrete: Part 6 - Simplified code-type formulations. Mater Struct 1992;25:219–23. http://dx.doi.org/10.1007/BF02473066. [45] Hwang YK, Bolander JE, Hong J-W, Lim YM. Irregular lattice model for geometrically nonlinear dynamics of structures. Mech Res Commun 2020;107:103554. http://dx.doi.org/10.1016/j.mechrescom.2020.103554. [46] Yamamoto Y, Isaji Y, Nakamura H, Miura T. Collapse simulation of reinforced concrete including localized failure and large rotation using extended RBSM. In: G. Pijaudier-Cabot, P. Grassl, C. La Borderie (Eds.) 10th international conference on fracture mechanics of concrete and concrete structures; 2019. p. 1–10. https://doi.org/10.21012/FC10.235632. [47] Cusatis G, Bažant ZP, Cedolin L. Confinement-shear lattice model for concrete damage in tension and compression: I. Theory. J Eng Mech 2003;129(12):1439–48. http://dx.doi.org/10.1061/(ASCE)0733-9399(2003)129:12(1439). [48] Cusatis G, Bažant ZP, Cedolin L. Confinement-shear lattice model for concrete damage in tension and compression: II. Computation and validation. J Eng Mech 2003;129(12):1449–58. http://dx.doi.org/10.1061/(ASCE)0733-9399(2003)129:12(1449). [49] Cusatis G, Bažant ZP, Cedolin L. Confinement-shear lattice CSL model for fracture propagation in concrete. Comput Methods Appl Mech Engrg 2006;195(52):7154–71. http://dx.doi.org/10.1016/j.cma.2005.04.019. [50] Cusatis G, Pelessone D, Mencarelli A. Lattice Discrete Particle Model (LDPM) for failure behavior of concrete. I: Theory. Cem Concr Comp 2011;33(9):881–90. http://dx.doi.org/10.1016/j.cemconcomp.2011.02.011. [51] Alnaggar M, Pelessone D, Cusatis G. Lattice Discrete Particle Modeling of reinforced concrete flexural behavior. J Struct Eng 2019;145(1):04018231. http://dx.doi.org/10.1061/(ASCE)ST.1943-541X.0002230. Engineering Fracture Mechanics 257 (2021) 108030 23 J.E. Bolander et al. [52] Cusatis G, Mencarelli A, Pelessone D, Baylot J. Lattice Discrete Particle Model (LDPM) for failure behavior of concrete. II: Calibration and validation. Cem Concr Comp 2011;33(9):891–905. http://dx.doi.org/10.1016/j.cemconcomp.2011.02.010. [53] Han L, Pathirage M, Akono A, Cusatis G. Lattice discrete particle modeling of size effect in slab scratch tests. J Appl Mech 2021;88:021009. http://dx.doi.org/10.1115/1.4048989. [54] Wan L, Wendner R, Liang B, Cusatis G. Analysis of the behavior of ultra high performance concrete at early age. Cem Concr Comp 2016;74:120–35. http://dx.doi.org/10.1016/j.cemconcomp.2016.08.005. [55] Esna Ashari S, Buscarnera G, Cusatis G. A lattice discrete particle model for pressure-dependent inelasticity in granular rocks. Int J Rock Mech Min 2017;91:49–58. http://dx.doi.org/10.1016/j.ijrmms.2016.10.007. [56] Li W, Jin Z, Cusatis G. Size effect analysis for the characterization of marcellus shale quasi-brittle fracture properties. Rock Mech Rock Eng 2019;52:1–18. http://dx.doi.org/10.1007/s00603-018-1570-6. [57] Feng J, Sun W-W, Li B-M. Numerical study of size effect in concrete penetration with LDPM. Defence Technology 2018;14(5):560–9. http://dx.doi.org/ 10.1016/j.dt.2018.07.006. [58] Smith J, Cusatis G. Numerical analysis of projectile penetration and perforation of plain and fiber reinforced concrete slabs. Int J Numer Anal Methods Geomech 2017;41(3):315–37. http://dx.doi.org/10.1002/nag.2555. [59] Smith J, Cusatis G, Pelessone D, Landis E, O’Daniel J, Baylot J. Discrete modeling of ultra-high-performance concrete with application to projectile penetration. Int J Impact Eng 2014;65:13–32. http://dx.doi.org/10.1016/j.ijimpeng.2013.10.008. [60] Eliáš J, Vořechovský M. Fracture in random quasibrittle media: I. Discrete mesoscale simulations of load capacity and fracture process zone. Eng Fract Mech 2020;235:107160. http://dx.doi.org/10.1016/j.engfracmech.2020.107160. [61] Eliáš J, Vořechovský M, Skoček J, Bažant ZP. Stochastic discrete meso-scale simulations of concrete fracture: comparison to experimental data. Eng Fract Mech 2015;135:1–16. http://dx.doi.org/10.1016/j.engfracmech.2015.01.004. [62] Fascetti A, Bolander JE, Nisticó N. Lattice discrete particle modeling of concrete under compressive loading: Multiscale experimental approach for parameter determination. J Eng Mech 2018;144(8):04018058. http://dx.doi.org/10.1061/(ASCE)EM.1943-7889.0001480. [63] Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Géotechnique 1979;29(1):47–65. http://dx.doi.org/10.1680/geot.1979.29. 1.47. [64] Boon CW, Houlsby GT, Utili S. A new algorithm for contact detection between convex polygonal and polyhedral particles in the Discrete Element Method. Comput Geotech 2012;44:73–82. http://dx.doi.org/10.1016/j.compgeo.2012.03.012. [65] Nezami EG, Hashash YMA, Zhao D, Ghaboussi J. A fast contact detection algorithm for 3-D Discrete Element Method. Comput Geotech 2004;31(7):575–87. http://dx.doi.org/10.1016/j.compgeo.2004.08.002. [66] Favier JF, Abbaspour-Fard MH, Kremmer M, Raji AO. Shape representation of axi-symmetrical, non-spherical particles in discrete element simulation using multi-element model particles. Eng Comput 1999;16(4):467–80. http://dx.doi.org/10.1108/02644409910271894. [67] Wang P, Gao N, Ji K, Stewart L, Arson C. DEM analysis on the role of aggregates on concrete strength. Comput Geotech 2020;119:103290. http://dx.doi.org/10.1016/j.compgeo.2019.103290. [68] Kawamoto R, Andò E, Viggiani G, Andrade JE. Level set discrete element method for three-dimensional computations with triaxial case study. J Mech Phys Solids 2016;91:1–13. http://dx.doi.org/10.1016/j.jmps.2016.02.021. [69] Lisjak A, Grasselli G. A review of discrete modeling techniques for fracturing processes in discontinuous rock masses. J Rock Mech Geotech Eng 2014;6(4):301–14. http://dx.doi.org/10.1016/j.jrmge.2013.12.007. [70] Abid N, Pro JW, Barthelat F. Fracture mechanics of nacre-like materials using discrete-element models: Effects of microstructure, interfaces and randomness. J Mech Phys Solids 2019;124:350–65. http://dx.doi.org/10.1016/j.jmps.2018.10.012. [71] Donzé FV, Magnier S-A, Daudeville L, Mariotti C, Davenne L. Numerical study of compressive behavior of concrete at high strain rates. J Eng Mech 1999;125(10):1154–63. http://dx.doi.org/10.1061/(ASCE)0733-9399(1999)125:10(1154). [72] Hentz S, Daudeville L, Donzé FV. Identification and validation of a discrete element model for concrete. J Eng Mech 2004;130(6):709–19. http: //dx.doi.org/10.1061/(ASCE)0733-9399(2004)130:6(709). [73] Nguyen TT, Bui HH, Ngo TD, Nguyen GD, Kreher MU, Darve F. A micromechanical investigation for the effects of pore size and its distribution on geopolymer foam concrete under uniaxial compression. Eng Fract Mech 2019;209:228–44. http://dx.doi.org/10.1016/j.engfracmech.2019.01.033. [74] Nitka M, Tejchman J. Modelling of concrete behaviour in uniaxial compression and tension with DEM. Granul Matter 2015;17(1):145–64. http: //dx.doi.org/10.1007/s10035-015-0546-4. [75] Beckmann B, Schicktanz K, Reischl D, Curbach M. DEM simulation of concrete fracture and crack evolution. Struct Concr 2012;13(4):213–20. http://dx.doi.org/10.1002/suco.201100036. [76] Nitka M, Tejchman J. A three-dimensional meso-scale approach to concrete fracture based on combined DEM with 𝜇CT images. Cem Concr Res 2018;107:11–29. http://dx.doi.org/10.1016/j.cemconres.2018.02.006. [77] Sinaie S, Heidarpour A, Zhao X-L. A micro-mechanical parametric study on the strength degradation of concrete due to temperature exposure using the discrete element method. Int J Solids Struct 2016;88–89:165–77. http://dx.doi.org/10.1016/j.ijsolstr.2016.03.009. [78] Sinaie S, Ngo TD, Nguyen VP. A discrete element model of concrete for cyclic loading. Computer Struct 2018;196:173–85. http://dx.doi.org/10.1016/j. compstruc.2017.11.014. [79] Šmilauer et al. V. Yade documentation. 2nd ed.. The Yade project; 2015, http://dx.doi.org/10.5281/zenodo.34073,http://yade-dem.org/doc/. [80] Kloss C, Goniva C, Hager A, Amberger S, Pirker S. Models, algorithms and validation for opensource DEM and CFD–DEM. Progr Comput Fluid Dynamics Int J 2012;12(2–3):140–52. http://dx.doi.org/10.1504/PCFD.2012.047457. [81] Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. J Mech Phys Solids 2000;48(1):175–209. http://dx.doi.org/10. 1016/S0022-5096(99)00029-0. [82] Silling SA, Lehoucq RB. Peridynamic theory of solid mechanics. Adv Appl Mech 2010;44:73–168. http://dx.doi.org/10.1016/S0065-2156(10)44002-8. [83] Gerstle W, Sau N, Silling S. Peridynamic modeling of concrete structures. Nucl Eng Des 2007;237(12):1250–8. http://dx.doi.org/10.1016/j.nucengdes. 2006.10.002. [84] Bažant ZP, Luo W, Chau VT, Bessa MA. Wave dispersion and basic concepts of peridynamics compared to classical nonlocal damage models. J Appl Mech 2016;83. http://dx.doi.org/10.1115/1.4034319, 111004–1. [85] Abo Dhaheer MS, Kulasegaram S, Karihaloo BL. Simulation of self-compacting concrete flow in the J-ring test using smoothed particle hydrodynamics (SPH). Cem Concr Res 2016;89:27–34. http://dx.doi.org/10.1016/j.cemconres.2016.07.016. [86] Mokhatar SN, Sonoda Y, Kueh ABH, Jaini ZM. Quantitative impact response analysis of reinforced concrete beam using the Smoothed Particle Hydrodynamics (SPH) method. Struct Eng Mech 2015;56(6):917–38. http://dx.doi.org/10.12989/sem.2015.56.6.917. [87] Rabczuk T, Eibl J. Simulation of high velocity concrete fragmentation using SPH/MLSPH. Int J Numer Methods Eng 2003;56:1421–44. http://dx.doi.org/ 10.1002/nme.617. [88] Wu CT, Wu Y, Crawford JE, Magallanes JM. Three-dimensional concrete impact and penetration simulations using the smoothed particle Galerkin method. Int J Impact Eng 2017;106:1–17. http://dx.doi.org/10.1016/j.ijimpeng.2017.03.005. [89] Okabe A, Boots B, Sugihara K. Spatial tessellations: Concepts and applications of Voronoi diagrams. Chichester, England: John Wiley & Sons; 1992. [90] Aurenhammer F. Power diagrams: Properties, algorithms and applications. SIAM J Comput 1987;16(1):78–96. http://dx.doi.org/10.1137/0216006. Engineering Fracture Mechanics 257 (2021) 108030 24 J.E. Bolander et al. [91] Asahina D, Bolander JE. Voronoi-based discretizations for fracture analysis of particulate materials. Powder Technol 2011;213(1):92–9. http://dx.doi.org/ 10.1016/j.powtec.2011.07.010. [92] Landis EN, Keane DT. X-ray microtomography. Mater Charact 2010;61(12):1305–16. http://dx.doi.org/10.1016/j.matchar.2010.09.012. [93] Pan Z, Chen A, Ma R, Wang D, Tian H. Three-dimensional lattice modeling of concrete carbonation at meso-scale based on reconstructed coarse aggregates. Constr Build Mater 2018;192:253–71. http://dx.doi.org/10.1016/j.conbuildmat.2018.10.052. [94] Landis EN, Nagy EN, Keane DT. Microstructure and fracture in three dimensions. Eng Fract Mech 2003;70(7–8):911–25. http://dx.doi.org/10.1016/S00137944(02)00157-1. [95] Holla V, Vu G, Timothy JJ, Diewald F, Gehlen C, Meschke G. Computational generation of virtual concrete mesostructures. Materials 2021;14:3782. http://dx.doi.org/10.3390/ma14143782. [96] Qian Z, Garboczi EJ, Ye G, Schlangen E. Anm: a geometrical model for the composite structure of mortar and concrete using real-shape particles. Mater Struct 2016;49:149–58. http://dx.doi.org/10.1617/s11527-014-0482-5. [97] Asahina D, Landis EN, Bolander JE. Modeling of phase interfaces during pre-critical crack growth in concrete. Cem Concr Comp 2011;33(9):966–77. http://dx.doi.org/10.1016/j.cemconcomp.2011.01.007. [98] Grassl P, Rempling R. A damage-plasticity interface approach to the meso-scale modelling of concrete subjected to cyclic compressive loading. Eng Fract Mech 2008;75(16):4804–18. http://dx.doi.org/10.1016/j.engfracmech.2008.06.005. [99] Grassl P, Grégoire D, Rojas Solano L, Pijaudier-Cabot G. Meso-scale modelling of the size effect on the fracture process zone of concrete. Int J Solids Struct 2012;49(13):1818–27. http://dx.doi.org/10.1016/j.ijsolstr.2012.03.023. [100] Grégoire D, Verdon L, Lefort V, Grassl P, Saliba J, Regoin J-P, Loukili A, Pijaudier-Cabot G. Mesoscale analysis of failure in quasi-brittle materials: comparison between lattice model and acoustic emission data. Int J Numer Anal Methods Geomech 2015;39(15):1639–64. http://dx.doi.org/10.1002/ nag.2363. [101] Gu X, Hong L, Wang Z, Lin F. A modified rigid-body-spring concrete model for prediction of initial defects and aggregates distribution effect on behavior of concrete. Comput Mater Sci 2013;77:355–65. http://dx.doi.org/10.1016/j.commatsci.2013.04.050, URL http://www.sciencedirect.com/science/article/ pii/S092702561300222X. [102] Eliáš J, Stang H. Lattice modeling of aggregate interlocking in concrete. Int J Fract 2012;175:1–11. http://dx.doi.org/10.1007/s10704-012-9677-3. [103] Benkemoun N, Hautefeuille M, Colliat J-B, Ibrahimbegovic A. Failure of heterogeneous materials: 3D meso-scale FE models with embedded discontinuities. Int J Numer Methods Eng 2010;82(13):1671–88. http://dx.doi.org/10.1002/nme.2816. [104] van Mier JGM, Chiaia BM, Vervuurt A. Numerical simulation of chaotic and self-organizing damage in brittle disordered materials. Comput Methods Appl Mech Engrg 1997;142:189–201. http://dx.doi.org/10.1016/S0045-7825(96)01128-0. [105] Eddy L, Nagai K. Numerical simulation of beam-column knee joints with mechanical anchorages by 3D rigid body spring model. Eng Struct 2016;126:547–58. http://dx.doi.org/10.1016/j.engstruct.2016.07.054. [106] Jiradilok P, Wang Y, Nagai K, Matsumoto K. Development of discrete meso-scale bond model for corrosion damage at steel-concrete interface based on tests with/without concrete damage. Constr Build Mater 2020;236:117615. http://dx.doi.org/10.1016/j.conbuildmat.2019.117615. [107] Xenos D, Grégoire D, Morel S, Grassl P. Calibration of nonlocal models for tensile fracture in quasi-brittle heterogeneous materials. J Mech Phys Solids 2015;82:48–60. http://dx.doi.org/10.1016/j.jmps.2015.05.019. [108] Bathurst RJ, Rothenburg L. Micromechanical aspects of isotropic granular assemblies with linear contact interactions. J Appl Mech 1988;55(1):17–23. http://dx.doi.org/10.1115/1.3173626. [109] Kuhl E, D’Addetta GA, Herrmann HJ, Ramm E. A comparison of discrete granular material models with continuous microplane formulations. Granul Matter 2000;2(3):113–21. http://dx.doi.org/10.1007/s100350050003. [110] Liao C-L, Chang T-P, Young D-H, Chang CS. Stress-strain relationship for granular materials based on the hypothesis of best fit. Int J Solids Struct 1997;34(31):4087–100. http://dx.doi.org/10.1016/S0020-7683(97)00015-2. [111] Eliáš J. Boundary layer effect on behavior of discrete models. Materials 2017;10:157. http://dx.doi.org/10.3390/ma10020157. [112] Yao C, Jiang Q-H, Shao J-F, Zhou C-B. A discrete approach for modeling damage and failure in anisotropic cohesive brittle materials. Eng Fract Mech 2016;155:102–18. http://dx.doi.org/10.1016/j.engfracmech.2016.01.012. [113] Eliáš J. Elastic properties of isotropic discrete systems: Connections between geometric structure and Poisson’s ratio. Int J Solids Struct 2020. http://dx.doi.org/10.1016/j.ijsolstr.2019.12.012, in press. [114] Alassi HT, Holt R. Relating discrete element method parameters to rock properties using classical and micropolar elasticity theories. Int J Numer Anal Methods Geomech 2012;36(10):1350–67. http://dx.doi.org/10.1002/nag.1056. [115] Rezakhani R, Cusatis G. Asymptotic expansion homogenization of discrete fine-scale models with rotational degrees of freedom for the simulation of quasi-brittle materials. J Mech Phys Solids 2016;88:320–45. http://dx.doi.org/10.1016/j.jmps.2016.01.001. [116] Bardet JP, Vardoulakis I. The asymmetry of stress in granular media. Int J Solids Struct 2001;38(2):353–67. http://dx.doi.org/10.1016/S0020-7683(00) 00021-4. [117] Christoffersen J, Mehrabadi MM, Nemat-Nasser S. A micromechanical description of granular material behavior. J Appl Mech 1981;48(2):339–44. http://dx.doi.org/10.1115/1.3157619. [118] Weber J. Recherches concernant les contraintes intergranulaires dans les milieux pulvérulents. Bull de Liaison Des Ponts-Et-ChaussÉEs 1966;20:1–20. [119] Schlangen E, Garboczi EJ. New method for simulating fracture using an elastically uniform random geometry lattice. Int J Eng Sci 1996;34(10):1131–44. http://dx.doi.org/10.1016/0020-7225(96)00019-5. [120] Cusatis G, Rezakhani R, Schauffert EA. Discontinuous Cell Method (DCM) for the simulation of cohesive fracture and fragmentation of continuous media. Eng Fract Mech 2017;170:1–22. http://dx.doi.org/10.1016/j.engfracmech.2016.11.026. [121] Bolander JE, Yip M, Moriizumi K, Kunieda M. Rigid-body-spring network modeling of cement-based composites. In: de Borst R, Mazars J, PijaudierCabot Gilles, van Mier J G M, editors. In: Fracture mechanics of concrete structures, 2, Lisse, The Netherlands: A. A. Balkema Publishers; 2001, p. 773–80. [122] Zhao G-F, Fang J, Zhao J. A 3D distinct lattice spring model for elasticity and dynamic failure. Int J Numer Anal Methods Geomech 2011;35(8):859–85. http://dx.doi.org/10.1002/nag.930. [123] Asahina D, Ito K, Houseworth JE, Birkholzer JT, Bolander JE. Simulating the Poisson effect in lattice models of elastic continua. Comput Geotech 2015;70:60–7. http://dx.doi.org/10.1016/j.compgeo.2015.07.013. [124] Asahina D, Aoyagi K, Kim K, Birkholzer JT, Bolander JE. Elastically-homogeneous lattice models of damage in geomaterials. Comput Geotech 2017;81:195–206. http://dx.doi.org/10.1016/j.compgeo.2016.08.015. [125] Rasmussen LL, de Assis AP. Elastically-homogeneous lattice modelling of transversely isotropic rocks. Comput Geotech 2018;104:96–108. http://dx.doi. org/10.1016/j.compgeo.2018.08.016. [126] Deryugin YY, Lasko GV. Field of stresses in an isotropic plane with circular inclusion under tensile stress. Engineering 2012;4(9):583–9. http: //dx.doi.org/10.4236/eng.2012.49074. [127] Celigueta MA, Latorre S, Arrufat F, Oñate E. Accurate modelling of the elastic behavior of a continuum with the Discrete Element Method. Comput Mech 2017;60(6):997–1010. http://dx.doi.org/10.1007/s00466-017-1453-9. Engineering Fracture Mechanics 257 (2021) 108030 25 J.E. Bolander et al. [128] Rojek J, Zubelewicz A, Madan N, Nosewicz S. The discrete element method with deformable particles. Int J Numer Methods Eng 2018;114(8):828–60. http://dx.doi.org/10.1002/nme.5767. [129] Rots JG, Belletti B, Invernizzi S. Robust modeling of RC structures with an ‘‘event-by-event’’ strategy. Eng Fract Mech 2008;75(3):590–614. http: //dx.doi.org/10.1016/j.engfracmech.2007.03.027. [130] Cusatis G, Cedolin L. Two-scale study of concrete fracturing behavior. Eng Fract Mech 2007;74(1–2):3–17. http://dx.doi.org/10.1016/j.engfracmech.2006. 01.021. [131] Grassl P, Bažant ZP. Random lattice-particle simulation of statistical size effect in quasi-brittle structures failing at crack initiation. J Eng Mech 2009;135:85–92. http://dx.doi.org/10.1061/(ASCE)0733-9399(2009)135:2(85). [132] Cusatis G. Strain-rate effects on concrete behavior. Int J Impact Eng 2011;38(4):162–70. http://dx.doi.org/10.1016/j.ijimpeng.2010.10.030. [133] Frantík P, Veselý V, Keršner Z. Parallelization of lattice modelling for estimation of fracture process zone extent in cementitious composites. Adv Eng Softw 2013;60–61:48–57. http://dx.doi.org/10.1016/j.advengsoft.2012.11.020. [134] Bishop JE. Simulating the pervasive fracture of materials and structures using randomly close packed voronoi tessellations. Comput Mech 2009;44(4):455–71. http://dx.doi.org/10.1007/s00466-009-0383-6. [135] Leon SE, Spring DW, Paulino GH. Reduction in mesh bias for dynamic fracture using adaptive splitting of polygonal finite elements. Int J Numer Methods Eng 2014;100:555–76. http://dx.doi.org/10.1002/nme.4744. [136] Manzoli OL, Maedo MA, Bitencourt Jr LAG, Rodrigues EA. On the use of finite elements with a high aspect ratio for modeling cracks in quasi-brittle materials. Eng Fract Mech 2016;153:151–70. http://dx.doi.org/10.1016/j.engfracmech.2015.12.026. [137] Bažant ZP, Oh BH. Crack band theory for fracture of concrete. Matériaux et Construction 1983;16(3):155–77. http://dx.doi.org/10.1007/BF02486267. [138] Berton S, Bolander JE. Crack band model of fracture in irregular lattices. Comput Methods Appl Mech Engrg 2006;195(52):7172–81. http://dx.doi.org/ 10.1016/j.cma.2005.04.020. [139] Bolander JE, Sukumar N. Irregular lattice model for quasistatic crack propagation. Phys Rev B 2005;71(9):094106. http://dx.doi.org/10.1103/PhysRevB. 71.094106. [140] Yamamoto Y, Nakamura H, Kuroda I, Furuya N. Simulation of crack propagation in RC shear wall using a 3D Rigid-Body-Spring Model with random geometry. In: Proceedings of the 8th international conference on fracture mechanics of concrete and concrete structures; 2013. [141] Luković M, Yang Y, Schlangen E, Hordijk D. On the potential of lattice type model for predicting shear capacity of reinforced concrete and SHCC structures. In: Proceedings of the 2017 fib symposium, maastricht, the Netherlands; 2017. p. 804–13. ISBN 978-3-319-59470-5. https://doi.org/10.1007/978-3-31959471-2_94. [142] Lale E, Rezakhani R, Alnaggar M, Cusatis G. Homogenization coarse graining (HCG) of the lattice discrete particle model (LDPM) for the analysis of reinforced concrete structures. Eng Fract Mech 2018;197:259–77. http://dx.doi.org/10.1016/j.engfracmech.2018.04.043. [143] Baktheer A, Aguilar M, Chudoba R. Microplane fatigue model MS1 for plain concrete under compression with damage evolution driven by cumulative inelastic shear strain. Int J Plast 2021;102950. http://dx.doi.org/10.1016/j.ijplas.2021.102950. [144] Park K, Paulino GH. Cohesive zone models: a critical review of traction-separation relationships across fracture surfaces. Appl Mech Rev 2011;64(6). http://dx.doi.org/10.1115/1.4023110, 060802–1. [145] Kupfer H, Hilsdorf HK, Rusch H. Behavior of concrete under biaxial stresses. ACI Struct J 1969;66(8):656–66. [146] Nguyen H, Pathirage M, Rezaei M, Issa M, Cusatis G, Bažant ZP. New perspective of fracture mechanics inspired by gap test with crack-parallel compression. Proc Natl Acad Sci 2020;117(25):14015–20. http://dx.doi.org/10.1073/pnas.2005646117. [147] Chang Z, Zhang H, Schlangen E, Šavija B. Lattice fracture model for concrete fracture revisited: Calibration and validation. Appl Sci-Basel 2020;10(14). http://dx.doi.org/10.3390/app10144822. [148] Javidan F, Shahbeyk S, Safarnejad M. Lattice discrete particle modeling of compressive failure in hollow concrete blocks. Comput Concr 2014;13:437–56. http://dx.doi.org/10.12989/cac.2014.13.4.437. [149] Mendoza Jr. R, Yamamoto Y, Nakamura H, Miura T. Numerical evaluation of localization and softening behavior of concrete confined by steel tubes. Struct Concr 2018;19(6):1956–70. http://dx.doi.org/10.1002/suco.201700266. [150] Landis EN. Microplanes and microstructure: Connecting abstractions and reality. Eng Fract Mech 2018;200:42–9. http://dx.doi.org/10.1016/j.engfracmech. 2018.07.015. [151] Wittel FK, Carmona HA, Kun F, Herrmann HJ. Mechanisms in impact fragmentation. Int J Fract 2008;154:105–17. http://dx.doi.org/10.1007/s10704008-9267-6. [152] Cusatis G, Mencarelli A, Pelessone D, Baylot J. Lattice discrete particle model (LDPM) for fracture dynamics and rate effect in concrete. ASCE 2008;315(42):1–11. http://dx.doi.org/10.1016/j.cemconcomp.2011.02.011. [153] Hentz S, Donzé FV, Daudeville L. Discrete element modelling of concrete submitted to dynamic loading at high strain rates. Comput Struct 2004;82:2509–24. http://dx.doi.org/10.1016/j.compstruc.2004.05.016. [154] Lilliu G, van Mier JGM. 3D lattice type fracture model for concrete. Eng Fract Mech 2003;70(7):927–41. http://dx.doi.org/10.1016/S0013-7944(02)001583. [155] Yip M, Li Z, Liao B-S, Bolander JE. Irregular lattice models of fracture of multiphase particulate materials. Int J Fract 2006;140:113–24. http: //dx.doi.org/10.1007/s10704-006-7636-6. [156] Pari M, Hendriks MAN, Rots JG. Non-proportional loading in sequentially linear analysis for 3D stress states. Int J Numer Methods Eng 2019;119(6):506–31. http://dx.doi.org/10.1002/nme.6060. [157] Pari M, Van de Graaf AV, Hendriks MAN, Rots JG. A multi-surface interface model for sequentially linear methods to analyse masonry structures. Eng Struct 2021;238:112123. http://dx.doi.org/10.1016/j.engstruct.2021.112123. [158] DeJong MJ, Hendriks MAN, Rots JG. Sequentially linear analysis of fracture under non-proportional loading. Eng Fract Mech 2008;75(18):5042–56. http://dx.doi.org/10.1016/j.engfracmech.2008.07.003. [159] Eliáš J, Frantík P, Vořechovský M. Improved sequentially linear solution procedure. Eng Fract Mech 2010;77(12):2263–76. http://dx.doi.org/10.1016/j. engfracmech.2010.05.018. [160] Graça-E-Costa R, Alfaiate J, Dias-da Costa D, Sluys LJ. A non-iterative approach for the modelling of quasi-brittle materials. Int J Fract 2012;178(1–2):281–98. http://dx.doi.org/10.1007/s10704-012-9768-1. [161] Liu J, Liang N. Algorithm for simulating fracture processes in concrete by lattice modeling. Theor Appl Fract Mech 2009;52(1):26–39. http://dx.doi.org/ 10.1016/j.tafmec.2009.06.004. [162] Yu C, Hoogenboom PCJ, Rots JG. Incremental sequentially linear analysis to control failure for quasi-brittle materials and structures including non-proportional loading. Eng Fract Mech 2018;202:332–49. http://dx.doi.org/10.1016/j.engfracmech.2018.07.036. [163] Pari M, Hendriks MAN, Rots JG. Non-proportional loading in sequentially linear solution procedures for quasi-brittle fracture: A comparison and perspective on the mechanism of stress redistribution. Eng Fract Mech 2020;230:106960. http://dx.doi.org/10.1016/j.engfracmech.2020.106960. [164] Eliáš J. Generalization of load-unload and force-release sequentially linear methods. Int J Damage Mech 2015;24(2):279–93. http://dx.doi.org/10.1177/ 1056789514531001. [165] Liu J, El Sayed TS. On the load-unload (L-U) and force-release (F-R) algorithms for simulating brittle fracture processes via lattice models. Int J Damage Mech 2012;21(7):960–88. http://dx.doi.org/10.1177/1056789511424585.