Engineering Fracture Mechanics xxx (2017) xxx–xxx Contents lists available at ScienceDirect Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images W. Trawiński, J. Tejchman ⇑, J. Bobiński Gdańsk University of Technology, Faculty of Civil and Environmental Engineering, 80-233 Gdańsk, Poland a r t i c l e i n f o Article history: Received 8 December 2016 Received in revised form 21 September 2017 Accepted 9 October 2017 Available online xxxx Keywords: Cohesive elements Concrete beam Fracture Interfacial transitional zones 3D meso-scale approach FEM X-ray lCT a b s t r a c t The paper shows three-dimensional numerical results regarding concrete fracture at aggregate level in a notched concrete beam under bending using cohesive interface elements. Concrete heterogeneity was taken into account by considering a 3-phase material description: aggregate, cement matrix and interfacial transitional zones (ITZs). The real shape and location of aggregate particles in concrete were assumed in numerical analyses based on Xray lCT scans. The 3D meso-scale finite element method (FEM) with cohesive elements embedded in the cement matrix and ITZs was used in simulations. The effect of different properties of ITZs and cement matrix was investigated in parametric studies. A satisfactory agreement in terms of the vertical force versus crack mouth opening displacement evolution and crack geometry was achieved between analyses and laboratory tests. Ó 2017 Elsevier Ltd. All rights reserved. 1. Introduction Fracture process in quasi-brittle materials such as concrete is very complex due to branching, coalescence, kinking, tortuousness and interlocking of cracks [1–3]. Experiments and corresponding numerical analyses with concrete under different loading conditions apparently indicate that the concrete fracture behaviour depends on its internal structure (e.g. aggregate volume, aggregate size, aggregate roughness, aggregate stiffness, particle size distribution curve, mortar volume and macro-porosity). Therefore, a realistic description of fracture in concrete can be achieved only when the abovementioned factors are taken into account. In particular, aggregate shape and placement should be considered as they are the major ingredients of the concrete mixture (approximately 70–80% of its volume, including sand particles) [4]. Since the onset of micro-cracks takes place in interfacial transitional zones (ITZs) between aggregate and cement matrix, thus it is crucial to properly determine their properties. The width of ITZs, which are a highly porous regions of the cement paste around aggregates, is between 0 and 50 mm depending upon the aggregate roughness [5–7]. Meso-scale simulations of concrete fracture at aggregate level may be performed using either continuous [8–14] or discontinuous models [15–18] within continuum mechanics. It is also possible to model the concrete behaviour with the use of discrete models which include lattice models [19–22], discrete element method (DEM) [23–26], interface element models with constitutive laws based on non-linear fracture mechanics [27,28] and rigid body-spring models [29–31]. Classical continuous models are sensitive to a finite element mesh discretization and thus they need to be enriched with a characteristic length of micro-structure in ⇑ Corresponding author. E-mail addresses: [email protected] (W. Trawiński), [email protected] (J. Tejchman), [email protected] (J. Bobiński). https://doi.org/10.1016/j.engfracmech.2017.10.003 0013-7944/Ó 2017 Elsevier Ltd. All rights reserved. Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 2 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx Nomenclature B beam width C bulk stiffness matrix d50 mean particle diameter D grain diameter D scalar damage variable E modulus of elasticity fc uniaxial compressive strength ft uniaxial tensile strength tensile strength under bending ft,flex F vertical force (load) GF tensile fracture energy hc crack height H beam height kn, ks, kt current normal and tangential diagonal components of stiffness matrix K kn0, ks0, kt0 initial normal and tangential diagonal components of stiffness matrix K K stiffness matrix of cohesive elements lc crack length L beam length T traction vector tn, ts, tt current normal and tangential components of traction vector tn0, ts0, tt0 critical normal and tangential components of traction vector U beam deflection a softening curve parameter d relative displacements vector dm effective relative displacement d0m effective relative displacement at damage initiation dmf effective relative displacement at de-cohesion maximum effective relative displacement dmax m dn, ds, dt normal and tangential components of relative displacement vector e strain vector v Poisson’s ratio r standard deviation r stress vector 2D two-dimensional 3D three-dimensional CMOD crack mouth opening displacement COH3D6 6-node three-dimensional cohesive finite element C3D4 4-node tetrahedral solid finite element DEM discrete element method FE finite element FEM finite element method ITZ interfacial transitional zone V volume lCT micro-computed tomography ‘a’ aggregate ‘cm’ cement matrix ‘max’ maximum ‘mid’ medium ‘min’ minimum order to provide mesh-independent results [32–40]. In contrast, discontinuous models using cohesive elements or the eXtended Finite Element Method for modelling concrete fracture do not need any extension, since a crack or shear zone is created in a discrete way [3,18,41,42]. The main aim of this paper was to carefully examine a complex fracture process in a concrete beam subjected to threepoint bending using a 3D meso-scale approach with zero thickness cohesive elements embedded in the initial mesh of solid (bulk) elements. Interface elements inserted into the cement matrix and along aggregate/cement matrix faces were equipped with a crack initiation criterion and traction-separation law [15–18]. A heterogeneous nature of concrete was taken into consideration by modelling 3 different phases: aggregate (with the diameter da 2 mm), cement matrix and ITZs Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx 3 (around aggregate particles) at the meso-scale level and assigning different material properties to each of them. First, aggregate particles (da 2 mm) were modelled as randomly distributed spheres of a varying diameter based on the sieve curve. In the second step, aggregate real shape and placement were considered in concrete based on X-ray lCT images of internal structure [4,13,18]. The segmentation process of aggregate was carried out using the free open-source software 3D Slicer [43,44] (the software package for visualization and medical image computing). The 3D numerical simulations were carried out with the commercial FE package Abaqus [45]. The effect of the mesh size and key fracture parameters on forcedisplacement curves and crack patterns was investigated. The presented results are a continuation of our 2D simulation results described in [18]. The developing of mesoscopic models for concrete is important because they can replace costly laboratory experiments. The innovative points of the present research works are: (1) description of a method generating the FE mesh for 3D aggregate particles of a real shape by means of generally available tools, (2) concrete fracture simulations with the use of a 3D mesoscopic continuum approach with real angularly-shaped aggregate particles (which is significantly more complex than a 2D approach), (3) a direct comparison of numerical outcomes regarding fracture with experimental images from the microcomputed tomography [13] and (4) a direct comparison between 2D and 3D numerical mesoscopic results. The paper is organized as follows. Laboratory tests on concrete beams under bending are described in Section 2. Details of fracture modelling in concrete beams are presented in Section 3. Section 4 provides insights into a 3D FE model generation. Results of 3D numerical simulations as compared to experimental ones are described in Section 5. Parametric studies are depicted in Section 6 and main conclusions are listed in Section 7. 2. Own experiments on concrete beams under loading The detailed description of laboratory tests using X-ray micro-computed tomography with two simply-supported notched concrete beams subjected to three-point bending may be found in [13]. The concrete mixture was composed of aggregate and sand grains (dmax = 16 mm, d50 = 2 mm), Portland cement (CEM I 32.5 R) and water. The water/cement, sand/cement and aggregate/cement ratios were 0.42, 0.80 and 1.43, respectively. Standard curing conditions with the relative humidity greater than 90% and temperature of 20 °C were applied to concrete specimens. The total aggregate and sand volumetric content was V = 75%.The rectangular beam cross-section had the height of H = 80 mm and width of B = 40 mm. The beam’s total length was L = 320 mm (4 H) and the span between the supports was 3 H = 240 mm. A macro-crack was triggered in the beam’s mid-part due to a notch with the height of H/10 = 8 mm and width of 3 mm (Fig. 1). Three months after the concrete mixing, the loading machine Instron 5569 was used for quasi-static bending tests on notched concrete beams. The tests were performed with a constant-controlled notch opening displacement rate of 0.002 mm/min (crack mouth opening displacement (CMOD)). Each test ended for CMOD = 0.10–0.15 mm prior to the unstable failure in order to provide the images of internal micro-structure by means of lCT. After the tests, six concrete specimens 10 10 10 cm3 were cut out from the concrete beams in order to determine the average uniaxial compressive strength fc (the maximum vertical force divided by the cross-section area). It was equal to fc = 51.81 MPa with the standard deviation r = 1.12 MPa. Next, other six cylindrical concrete specimens 1530 cm2 were used to measure the Young’s modulus E and Poisson’s ratio v, which were equal to E = 36.1 GPa (r = 0.46 GPa) and v = 0.22 (r = 0.01), respectively. The flexural tensile strength, calculated on the basis of four-point bending test of concrete beams 4 4 16 cm3, was equal to ft,flex = 3.7–4.3 MPa Fig. 1. Geometry and boundary conditions of concrete beams subjected to three-point bending. Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 4 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx (the maximum bending moment divided by the sectional modulus). The accompanying laboratory tests were carried out according to the standard regulations. Two typical experimental diagrams of the vertical force F versus CMOD for the concrete beams (with the displacement CMOD up to 0.15 mm for the beam ‘1’ and up to 0.10 mm for the beam ‘2’) are shown in Fig. 2A. The ultimate vertical force F was equal to 2.15–2.25 kN (the flexural tensile strength was 3.73–3.91 MPa). The distribution of aggregate grains and macro-voids in concrete beams was determined based on lCT-images of the concrete hexahedrons (80 50 40 mm3) extracted from the mid-part of each beam after tests [13]. Fig. 2B shows 3D images of concrete internal structure obtained with the aid of SkyScan 1173 which represents a new generation in high-resolution desktop X-ray micro-computed tomography systems [13]. Fig. 2. Experiments with concrete beam under bending: (A) measured vertical force F versus CMOD for two concrete beams, (B) X-ray lCT images of cracked hexahedral specimens 80 50 40 mm3 cut out from beam ‘1’ (a) and beam ‘2’ (b) and C) scanning electron micrographs of ITZs around aggregate particles (magnification factor 1000) [13] Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 5 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx The main discrete macro-crack was strongly curved along the beam height and width due to a random presence of aggregate grains (Fig. 2B) and it propagated mainly through the weakest phase in concrete which were interfacial transitional zones (ITZs). Thus micro-cracking occurred first in ITZs and when two interfacial cracks occurred around adjacent aggregates, a crack inside the cement matrix initiated to bridge the interfacial cracks so that a connected crack path was formed. Sometimes the crack propagated through macro-voids and very rarely through a single weak aggregate grain. A very nonuniform porous structure of ITZs as well as the presence of separated small sand grains was possible to be observed with the maximum magnification factor 30,000 using the scanning electro-microscope (SEM) Hitachi TM3030. The width of porous ITZs did not depend on aggregate particles diameter and changed between 30 and 50 lm [13]. 3. Concrete fracture modelling In order to simulate the fracture process, a cohesive crack model was applied in the beam’s mid-region. This approach was first introduced by Hillerborg et al. [46] based on the Dugdale model for steel [47]. The proposed model of cohesive elements allowed for describing cracks as jumps in a displacement field [48]. Although this approach does not require an enhancement with a characteristic length of micro-structure to describe cracks, it possesses several drawbacks. Crack trajectories are dominated by preferred mesh orientations [49,50]. Cracks are limited to bulk element boundaries and a diffuse cracking pattern is strongly mesh-dependent [39]. The initial stiffness of interface elements is a penalty parameter which may lead to traction oscillations in the elastic regime [51]. A stochastic distribution of the critical traction in cohesive elements and sufficient mesh-refinement may improve the aforementioned disadvantages [50–52]. A random distribution of aggregates contributes to mesh-independent results in terms of force-deflection curves and crack patterns [17,18]. The assumption of ITZs around aggregate particles improves the calculation convergence, since they act as attractors for cracks [17,18]. The meso-scale finite element method for concrete described as a three-four phase-material with cohesive elements turned out to be a very effective approach for modelling concrete fracture under uniaxial tension in both 2D [15,16] and 3D conditions [17] and under bending in 3D conditions [18] with randomly distributed spherical and angular aggregates. The FE mesh of the concrete beam (Section 2) was composed of the tetrahedral solid elements C3D4 and the interface elements COH3D6. Since the FE mesh included two element types, it was necessary to define the constitutive behaviour in each element type. 3.1. Solid elements A linearly elastic model was assumed for tetrahedral solid elements throughout the entire analysis. Thus, the only required material constants for continuum elements were: Young’s modulus E and Poisson’s ratio v. The constitutive matrix C linked the stress r with the strain e as: r ¼ C e: ð1Þ The material properties for all solid elements were the same as in our 2D numerical analyses [18]. The Young’s modulus of aggregate grains, Ea = 47.2 GPa, was calculated as the weighted average of the moduli of individual components [4,18]. For the cement matrix bulk elements the value of Ecm = 29.2 GPa was chosen and Eh = 36.1 GPa was taken for the homogeneous beam’s part wherein a 1-phase description was applied (see Fig. 9a). The Poisson’s ratio for each material phase was taken as v = 0.22. Nowadays, several multi-scale models for estimating the elastic stiffness, compressive strength and creep behaviour of the cement paste are available based on the amount of material constituents and a hydration degree [53–56]. These models are in particular attractive for complex analyses when the time-dependence is considered. In our calculations we used simply a weighted average value to compute the Young’s modulus of the cement matrix, since the preliminary parametric studies revealed its small influence on the structural beam response [18]. When using the formula established e.g. from the method by Mori-Tanaka and generalized self-consistent method [57], the Young’s modulus of the homogeneous concrete beam part, Eh, was Eh ¼ Ecm þ V a ðEa Ecm Þ 0:375ð47:2 29:2Þ ¼ 29:2 þ ¼ 34:6 GPa 47:229:2 1 þ ð1 0:375Þ 29:2þ411:97=3 1 þ ð1 V a Þ Ea E4lcmcm Ecm þ ð2Þ 3 with (mcm = m = 0.22) lcm ¼ Ecm 29:2 ¼ 11:97: ¼ 2ð1 þ mcm Þ 2ð1 þ 0:22Þ ð3Þ The value of 34.6 GPa in Eq. (2) was solely smaller by 4% than the value of 36.1 GPa that was assumed in simulations. 3.2. Cohesive elements The behaviour of cohesive elements was governed by the following equation: t ¼ Kd ð4Þ Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 6 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx Table 1 Basic material properties used in 3D FE analyses of concrete fracture at aggregate level. Solid elements Constant Aggregate (da 2 mm) Cement matrix Homog. region E [GPa] v [–] 47.2 0.2 29.2 0.2 36.1 0.2 Cohesive elements Constant Cement matrix ITZs around aggregate k0 [MPa/mm] ft [MPa] GF [N/m] a [–] 106 4.4 40 7.5 0.071 106 1.6 20 7.5 0.098 f dm [mm] Fig. 3. Laboratory aggregate size distribution curve [13] discretized into 3 regions: da,min = 6 mm (aggregate with diameter da = 4–8 mm), da,mid = 10 mm (aggregate with diameter da = 8–12 mm) and da,max = 4 mm (aggregate with diameter da = 12–16 mm). Fig. 4. Exemplary FE model of spherical aggregate particles with different volumetric content Va: (a) Va = 25% and (b) Va = 35%. Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx 7 or 2 tn 3 2 kn 6 7 6 4 ts 5 ¼ 4 0 tt 0 0 ks 0 0 32 dn 3 76 7 0 54 ds 5; dt kt ð5Þ where the tractions across the discontinuity t were related to the relative displacements d via the stiffness matrix K. The symbols tn, ts and tt are the tractions, kn, ks and kt denote the interface stiffnesses, dn, ds and dt are the relative displacements Fig. 5. Real concrete meso-structure displayed in three orthogonal sections (a–c) and in 3D view (d) for segmentation process by using 3D Slicer [43,44]. Fig. 6. 3D aggregate model created in program Slicer [43,44]. Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 8 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx for a given load increment in a normal and tangential directions. The stiffnesses kn, ks and kt were taken before cracking as the initial stiffnesses kn0, ks0 and kt0, respectively. The issue of the appropriate value of the initial stiffness of cohesive elements was discussed in [17,18,58,59]. In our numerical simulations the initial stiffness in a normal and tangential directions was chosen as kn0 = ks0 = kt0 = 106 MPa/mm for all interface types, similarly to 2D computations by Wang et al. [16], Trawiński et al. [18] and Lopez et al. [60] and in 3D studies by Wang et al. [17] and Caballero et al. [28]. The assumed penalty stiffness ensured ITZs to be intact (almost without relative displacements) in the elastic deformation range. After cracking the normal kn and shear ks, kt components of the stiffness matrix K were influenced by the scalar variable D (defined by Eq. (10) in Section 3.3) according to: kn ¼ ð1 DÞkn0 ; ks ¼ ð1 DÞks0 ; kt ¼ ð1 DÞkt0 : ð6Þ 3.3. Crack initiation and propagation A crack nucleated when the quadric nominal stress criterion was fulfilled: 2 2 2 ht n i ts tt þ þ ¼ 1: t n0 t s0 tt0 ð7Þ The critical traction for the fracture mode I tn0, fracture mode II ts0 and fracture mode III tt0 was different for each interface type. Based on the 2D numerical simulations at the aggregate level with cohesive elements by Trawiński et al. [18], the following tensile strength constants were assumed in most of analyses: tn0,ITZ (ft,ITZ) = 1.6 MPa and tn0,cm (ft,cm) = 4.4 MPa (ft,ITZ/ft, cm = 0.35) for the inter-phase (aggregate/cement matrix) and intra-phase (cement matrix/cement matrix) interface elements. In addition, the FE calculations were carried out with ft,ITZ/ft,cm = 0.5, ft,ITZ/ft,cm = 0.7 and ft,ITZ/ft,cm = 1 (ft,cm = 4.4 MPa). The tensile strength of the cement matrix was taken from own laboratory bending tests using the formula that related the uniaxial tensile strength with the flexural one [4]. In general it is possible to experimentally determine the strength properties of aggregate and cement matrix. However the determination of ITZ-properties is extremely difficult in experiments. Therefore some parametric studies were needed to estimate the tensile strength of ITZs with respect to the tensile strength of the Fig. 7. FE surface model of real angular aggregate grains: (A) description of single aggregate particle: (a) initial shape, (b) fine mesh description and (c) coarse mesh description and (B) description of assembly of aggregate particles using: (a) fine mesh description and (b) coarse mesh description. Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx 9 Fig. 8. FE model generation: (a) views on initial surface mesh of single particle, (b) conversion of surface mesh into solid mesh and (c) single aggregate particle embedded in cement matrix. Fig. 9. Concrete beam divided into homogeneous part and heterogeneous one described as 3-phase mid-region (a) and prismatic cohesive elements between tetrahedral solid elements (b). Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 10 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx cement matrix. Due to the lack of experimental results, the shear critical tractions ts0 and tt0 for all cohesive elements were assumed to be the same as the normal ones (similarly as in numerical analyses by Ren et al. [15], Wang et al. [16,17], Trawiński et al. [18], Su et al. [61] and Yang et al. [62]). Two-dimensional parametric studies by Trawiński et al. [18] revealed that the specimen strength slightly increased with growing shear strength ts0 while the macro-crack path shape remained unaffected. In order to describe damage under a combination of normal and shear deformations across the interface, the effective relative displacement dm was introduced: dm ¼ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ hdn i2 þ d2s þ d2t ; ð8Þ where h i is the Macaulay bracket: hxi ¼ 0; x < 0 xP0 x; : ð9Þ The exponential softening curve was chosen to describe the behaviour of the interface elements in a post-peak regime. The scalar damage variable D, dependent on the effective relative displacement, was calculated as [18]: ( D¼1 d0m dmax m 9 8 max 0 > )> > 1 exp a dmf d0m > < = dm dm 6 dmf ; 1 for dmax m > > 1 expðaÞ > > : ; ð10Þ > dmf ; D ¼ 1:0 for dmax m where d0m and dmf are the effective relative displacements at the crack nucleation and complete damage, respectively. The slope of the softening curve was affected by the parameter a whereas the maximum effective relative displacement obtained Fig. 10. Vertical force-deflection curves obtained from numerical analyses with spherical aggregate particles for 3 realizations with aggregate volumetric content Va: (a) Va = 25% and (b) Va = 35%. Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 11 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx during the loading history was dmax m . The fracture energy was assumed as GF,ITZ = 20 N/m (ITZs) and GF,cm = 40 N/m (cement matrix) [15–18]. The parameters dmf and a in Eq. (10) were the same as in our previous 2D fracture simulations [18]. Summarized, the constitutive model has the following constants: Ea, va, Ecm, vcm, kcm0, kITZ0, ft,cm, GF,cm, ft,ITZ, GF,ITZ, Eh and vh (’a’ - aggregate, ’cm’ - cement matrix, ’ITZ’ - interfacial transitional zone and ’h’ - homogeneous beam part). The calculations were mainly carried out with the constants (called standard) listed in Table 1. The weakest phase were obviously ITZs (ft,ITZ/ft, cm = 0.35). Their strength was rather low due to the fact that their thickness was assumed to be zero in order to reduce the computational effort, although their measured value was between 30 and 50 lm [13]. The three-point bending test was simulated by enforcing the vertical displacement of the loading plate with the width of 3 mm up to its final value of 0.1 mm. a) a) b) b) c) c) left crack face right crack face A) left crack face right crack face B) Fig. 11. Crack patterns from 3D numerical simulations with spherical aggregate for aggregate volume fraction of: (A) Va = 25% and (B) Va = 35% for 3 realisations: ‘1’ (a), ‘2’ (b) and ‘3’ (c) of Fig. 10 (green colour - cement matrix, blue colour - aggregate). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 12 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx 4. FE mesh generation In mesoscopic 3D FE calculations for concrete we considered two different aggregate descriptions with: (1) randomly located spherical aggregate grains and (2) angularly-shaped particles based on lCT scans of the real internal concrete structure (Fig. 2B). The aggregate particles were located in the mid-region of the concrete beam only. In the first case (initial analyses), the laboratory aggregate size distribution curve was discretized (Fig. 3) and 3 different aggregate diameters were distinguished: da,min = 6 mm (da = 4–8 mm), da,mid = 10 mm (da = 8–12 mm) and da,max = 14 mm (da = 12–16 mm).The data regarding the random meso-structure (co-ordinates of particle centres and their diameters) were calculated using the object-oriented application (see Box 1) and provided in the form of a Python script that was next used in Abaqus [45] to generate meso-structural mesh data. An exemplary FE model with spherical aggregate particles with a different volumetric content is shown in Fig. 4. In the second case, the real shape and placement of aggregates based on lCT images of the internal structure were taken into account. The segmentation process of aggregate grains was carried out according to the procedure listed in Box 2. At this stage, there were around 300 particles described by a very dense triangular surface mesh. Next, the free software Autodesk Meshmixer [63] working with triangle meshes was used in order to smoothen aggregate shapes and to reduce the number of nodes and triangular faces. Each 3D model of a single aggregate particle was reduced to a reasonable mesh accuracy level in order to decrease the computational time of FE analyses (Fig. 7A). Two mesh accuracy levels were assumed for each particle: a coarse one - described by 100 triangles and a fine one - described by 300 triangles. All surface models of aggregates were merged together into a single ⁄.stl file (Fig. 7B) that served as the input data for creating a 3D concrete model for computations by Abaqus [45]. The entire generation process of the full 3D FE model consisted of two main stages (presented in Box 3). In the first stage, the initial FE-mesh without cohesive elements was defined (within Abaqus CAE). The surface mesh of the beam’s mid-part Box 1 Algorithm for generating the random spherical structure. discretize the laboratory sieve curve (da,min, da,mid, da,max), calculate the number of spheres for each diameter (nmin, nmid, nmax), initialize the pseudo-random number generator based on the current time, place the particles inside the beam’s mid-region by starting from the spheres with the maximum diameter, a. get the trial centre’s coordinates using the pseudo-random number generator and check if the new particle does not overlap the existing ones, b. if the new sphere overlaps at least one particle, apply the translate method (see below) until the non-overlapping condition is satisfied, c. add the particle to the located particle list and go back to the point ‘4a’ until all spheres are located in the beam’s mid-region, (5) prepare the Python script for further pre-processing in Abaqus [45]. (1) (2) (3) (4) translate method uses the following procedure: move the particle by changing only one coordinate while keeping the remaining ones equal to the trial ones (e.g. keep ‘x’ and ‘z ‘as the trial coordinates and move the particle vertically only); decision on the choice of the axis direction is made at random, if the allowable position is still not found, change simultaneously all centre coordinates; decision on the choice of both the coordinate axis and coordinate increment direction is made at random. Box 2 Segmentation process procedure. reduce the resolution of lCT scans by using the re-size option in CT-Analyser, convert the images into the DICOM (Digital Imaging and Communications in Medicine) format using SkyScan’s Data Viewer, import the DICOM files to the free open-source software 3D Slicer [43,44] and reconstruct the hexahedral concrete volume (Fig. 5), extract manually the aggregate grains from the cement matrix and create the 3D surface models of particles using the Editor and Model Maker module in 3D Slicer [43,44] (Fig. 6), save the aggregate grain models in the file format *.stl (stereolithography - the form of the 3D printing technology). Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx 13 Box 3 Generation of the finite element model. Stage 1: definition of the initial mesh (Abaqus CAE) import the aggregate surface model and prepare the surface mesh of the beam’s mid-part, convert the surface meshes into the solid ones (Fig. 8a and b), merge two solid meshes (aggregate and cement matrix) in the Assembly Module of Abaqus [45] (Fig. 8c), create the homogeneous beam’s part and attach it to the heterogeneous region (Fig. 9a). Stage 2: addition of the cohesive elements. save the input file containing information regarding the initial mesh topology, based on the input file, prepare the following text files: – nodes_old (data on node coordinates creating the initial mesh), – agg_old, cm_old, h_old (the solid element nodal connectivity for aggregates, cement matrix and homogeneous parts, respectively), – mid_part (3-phase beam’s region: top and bottom face coordinates and left and right face abscissas). run the object-oriented application with the above files as the initial data: – define the tables/lists of the nodes and elements, – find the elements connected to each node in the heterogeneous beam’s region, – find the adjacent elements for all cement matrix tetrahedrons and renumber the nodes, – generate the cohesive elements between the cement matrix/cement matrix and cement matrix/aggregate, – save the new mesh topology, prepare the new input file for Abaqus simulations. was prepared with the aid of the aggregate surface model. First, the surface hexahedron (50 80 40 mm3) was defined. Next two surface meshes (hexahedron and aggregates) were artificially merged (since there were no common nodes to be combined) to form one region - a surface hexahedron with holes at aggregate locations. As a result after converting the surface mesh of the beam’s mid-region into the solid one, there were empty spots for further merging with the solid aggregate model. In the second stage, cohesive elements were added to the model. Cohesive elements were placed between bulk elements in the beam’s mid-part (volume 50 80 40 mm3) only. The following interfaces were formed: cement matrix/cement matrix and aggregate/cement matrix (Fig. 9b). A crack was allowed to nucleate in each interface between the cement matrix elements or in ITZ-elements between the aggregate and cement matrix. The crack propagation path was only constrained by solid elements’ boundaries and thus it could not appear in aggregates. 5. Meso-scale 3D numerical simulations The most frequent method to generate a FE model in Abaqus CAE [45] is to create a geometrical object (a collection of cells, faces and edges) and next to mesh it. If a meso-scale approach is used, two concrete phases have to be at least prepared: Fig. 12. Effect of ratio between tensile strength of ITZs and cement matrix on vertical force-deflection curve evolution in 3D computations (spherical aggregate particles with Va = 35%, ft,cm = 4.4 MPa). Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 14 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx aggregate and cement matrix. Solid parts of these two phases need to be created and assembled using a Boolean operation of merging the geometry with retaining intersecting boundaries. As a result there are not overlapping regions since each area belongs to aggregates or a cement matrix. This approach was solely used in our simulations with spherical particles since aggregate grains might be generated in Abaqus based on input data (particle centre coordinates and particle diameters). However, the objects with complex shapes are usually imported as the discretized ones with the aid of the other software. Fig. 13. Macro-crack propagation from 3D FE simulations for different ft,ITZ/ft,cm ratio (spherical aggregate, Va = 35%, ft,cm = 4.4 MPa): (A) 100%, (B) 50% and (C) 35% in three representative cross-sections: (a) at depth of 10 mm, (b) at depth of 20 mm and (c) at depth of 30 mm (red colour - crack, green colour cement matrix, blue colour - aggregate). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx 15 In our analyses with true aggregate, the surface model of angular particles was prepared by Autodesk Meshmixer [63]. In the Abaqus [45] nomenclature, a mesh of an imported discretized part is called an orphan one and both Boolean operations and partial remeshing are not possible to be performed due to the lack of information concerning the associated geometrical object since the mesh is regarded as a collection of nodes and elements. In order to obtain a model that is equivalent to the one generated based on geometrical objects, the procedure outlined in the Box 3 was applied; the model was assembled by merging common nodes of aggregate grains and a cement matrix instead of performing Boolean operations. In contrast to our 2D numerical analyses, macro-voids were neglected in present calculations. The generation of macro-voids around aggregate grains requires namely two Boolean operations: (1) producing holes (cutting geometry operation) in the cement matrix where macro-voids existed and (2) embedding aggregate grains in the cement matrix (merging the geometry operation with retaining intersecting boundaries). A possible way to include macro-voids is to assume their presence at a small A) B) a) b) left crack face right crack face Fig. 14. 3D numerical results of concrete beam with true angular aggregate: (A) vertical force-CMOD curves and (B) crack paths for coarse (a) and fine mesh description of particles (b) (green colour - cement matrix, blue colour - aggregate). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 16 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx distance from aggregate grains (by creating very thin regions of the cement matrix between macro-voids). However this operation contributes to a very dense FE mesh and higher computation time. Contrary to our 2D simulations [18], an explicit solver was used in our 3D numerical simulations in order to reduce the computation time, similarly as in the papers [15–17,61,62]. Dynamic effects were negligible since a very long loading time was assumed. It was chosen by monitoring the ratio between the kinetic and internal energy and shape of force-deflection curves. The time increment Dt = 106 s was assumed. The default values of the numerical linear and quadric bulk viscosity were chosen [45]: 0.06 and 1.2. The linear bulk viscosity was responsible for dumping oscillations occurring when a shock wave propagated through a mesh since it generated the bulk viscosity pressure that was linear with the volumetric strain Fig. 15. Calculated macro-crack propagation in beam’s mid-region above notch for coarse (A) and fine 3D aggregate discretization (B) (red colour - crack, green colour - cement matrix, blue colour - aggregate) for 5 different cross-sections as compared with 2D results [18] (C) and experimental results (D): (a) at depth of 5 mm, (b) at depth of 10 mm, (c) at depth of 20 mm, (d) at depth of 30 mm and (e) at depth of 35 mm. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx 17 rate. The quadric bulk viscosity smeared the shock over several elements and it was introduced to prevent elements from collapsing under extremely high velocity gradients. This bulk viscosity was quadratic with the volumetric strain rate. 5.1. Initial calculations with spherical aggregates First, the calculations were performed with random spherical aggregate particles of a varying diameter according to the discretized sieve curve of Fig. 3. The force-deflection curves for the different aggregate volume fraction (Fig. 10a and b) show that a global beam response was almost insensitive to a stochastic particles’ location. The beam strength decreased with increasing aggregate content due to a larger number of weak ITZs (Fig. 10). The macro-crack propagated mainly through the interfaces between aggregate grains and cement matrix (Fig. 11). The beam strength and brittleness decreased with decreasing tensile strength of ITZs (Fig. 12). The main crack was more curved and touched more aggregates with diminishing tensile strength of ITZs (Fig. 13). Fig. 16. Vertical force-CMOD curves from FE analyses of concrete beam as compared to experiment: (A) 3D result with standard parameters, (B) 2D results standard parameters in 3 different cross-sections [18] for: (a) depth of 5 mm, (b) depth of 20 mm and (c) depth of 35 mm and (C) 3D result with improved material parameters. Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 18 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx 5.2. Calculations with true aggregate Next, the calculations with a real grain size distribution were performed. Initially, in order to check the influence of the mesh density on numerical outcomes, the aggregate particles in the range 2 mm da 16 mm were described with a coarse and fine mesh of Fig. 8. Although the force-CMOD curves were mesh-independent (Fig. 14A), the calculated crack paths slightly depended on the mesh accuracy used for a particle description (Fig. 14B). For a fine mesh description, the macrocrack path was slightly more curved. Since a better agreement between the FE results and experimental ones with respect to the crack geometry was achieved for a fine discretization (Fig. 15), it was used in further parametric studies (300 triangles per one aggregate particle). In order to complete a single 3D analysis with the model consisting of 1,770,984 nodes and 1,831,363 elements (the approximate finite element size was 1.15 mm in the heterogeneous part), we needed 3 days using the parallel processing with 8 CPUs (3.2 GHz) or 20 days when using PC 3.27 GHz with one processor. Fig. 16a shows that the calculated force-CMOD curve was in a satisfactory agreement with the laboratory test result for the same cohesive element properties (called standard parameters) as in the corresponding 2D analyses [18]. The slope of the softening curve obtained from the numerical analysis was the same as in the experiment. The maximum vertical force Fmax = 2.45 kN was higher by 10% than the experimental one (Fmax = 2.25 kN). The calculated crack height hc varied between 56 mm and 58 mm for CMOD = 0.10 mm (Fig. 15) (in the experiment hc = 45–56 mm). The crack tortuosity, expressed by the crack length lc, was between 66 and 71 mm (in the experiment lc = 57–63 mm). The fractured zone (wherein cohesive element stiffness was decreased by more than 95%) included 29,022 cohesive elements while the number of interface elements forming the traction-free crack (D = 1.0) was 6390. In general, very good agreement was achieved with the experiment (Fig. 15). In addition, a better agreement was achieved for 3D numerical analyses than for corresponding 2D simulations [18] (Figs. 15 and 16). For 3D conditions the slope of the force-deflection curve in a post-peak regime was more realistic as compared to the 2D results. The macro-crack paths obtained from 3D numerical analyses were also closer to the experimental results than outcomes from 2D meso-scale simulations (wherein the calculated crack pattern was identical to the experiment only in a cross-section at the depth of 5 mm). The 3D numerical analyses provided very good agreement regarding the crack geometry in several cross-sections (at the depth of 5, 20, 30 and 35 mm). Note that in 2D slices obtained from 3D FE analyses some aggregate grains were missing due to particle mesh simplification prior to the numerical simulations. Finally it may be seen from Fig. 17 that the calculated crack’s width diminished along the beam height (from 0.1 mm at the bottom down to zero at the top) due to bending. In order to obtain a better agreement in terms of the force-CMOD curve, the material parameters for cohesive elements (called improved parameters) were updated: ft,ITZ = 1.6 MPa, ft,cm = 3.6 MPa, GF,ITZ = 20 N/m and GF,cm = 25 N/m (Fig. 16C). 6. Parametric study The meso-scale modelling might be used for a practical design of concrete with the improved desired performance in terms of both the higher tensile strength and ductility (instead of costly laboratory tests) by investigating the quantitative effect of all 4 concrete meso-phases on the concrete behaviour. Below we investigated the influence of fracture properties of the cement matrix and ITZs (ft,cm, GF,cm, ft,ITZ, GF,ITZ) on the global (force-deflection curve) and local (crack pattern) beam response for the true aggregate, since numerical outcomes revealed the importance of the accurate determination of these two individual phases’ properties in our meso-scale model. Due to the lack of the experimental data, the tensile strength of ITZs was estimated with the aid of parametric studies by comparing crack paths and force-deflection curves with laboratory test outcomes for the given tensile strength of the cement matrix. H Fig. 17. Crack width change above notch along beam height H from FE analyses (with displacement scale attached in [mm]). Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx 19 Fig. 18. 3D FE analyses: influence of tensile strength of ITZs (ft,cm = 4.4 MPa) on: (A) vertical force-deflection curve and (B) macro-crack propagation in beam’s mid-region for different ft,ITZ/ft,cm-ratio of 100–35% (from top to bottom) in 5 different cross-sections for: (a) depth of 5 mm, (b) depth of 10 mm, (c) depth of 20 mm, (d) depth of 30 mm and (e) depth of 35 mm (red colour - crack, green colour - cement matrix, blue colour - aggregate). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 20 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx 6.1. Effect of ITZ properties The beam’s strength and brittleness increased with the higher critical traction of ITZs (Fig. 18A). The macro-crack changed its shape from an almost straight to a strongly curved one with decreasing tensile strength of ITZs from 1 down to 0.35–0.70 with respect to the cement matrix strength. (Fig. 18B). The calculated geometry of the macro-crack was almost identical with ft,ITZ/ft,cm = 0.35 and ft,ITZ/ft,cm = 0.50 and very similar with ft,ITZ/ft,cm = 0.70. Thus, the calculated crack shape was in agreement with the experiment for the ratio of ft,ITZ/ft,cm = 0.35–0.70. The higher the fracture energy of ITZs, the larger was the maximum vertical force (Fig. 19A). In contrast to our 2D simulations [18], the dominant crack pattern was affected by the change of the ITZ-fracture energy (Fig. 19B). In 2D simulation outcomes the crack propagated however on the other side of several aggregate grains in cross-sections at the depth of 20 and 35 mm. The best agreement between the FE 3D results and experiment was obtained for ft,ITZ = 1.6 MPa (ft,cm = 4.4 MPa) and GF,ITZ = 20 N/m (GF,cm = 40 N/m). 6.2. Effect of cement matrix properties The strength growth of the cement matrix contributed to the increase of both the beam strength and brittleness (Fig. 20A). The macro-crack was more curved with increasing strength of the cement matrix (Fig. 20B). The increase of Fig. 19. 3D FE analyses: influence of fracture energy of ITZs (GF,cm = 40 N/m) on: (A) vertical force-deflection curve and (B) macro-crack propagation in beam’s mid-region: (a) GF,ITZ/GF,cm = 100% (top row) and GF,ITZ/GF,cm = 50% (bottom row) in 5 different cross-sections for: (a) depth of 5 mm, (b) depth of 10 mm, (c) depth of 20 mm, (d) depth of 30 mm and (e) depth of 35 mm (red colour - crack, green colour - cement matrix, blue colour - aggregate). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx 21 Fig. 20. 3D FE analyses: influence of tensile strength of cement matrix (ft,ITZ = 1.6 MPa) on: (A) force-deflection curve and (B) macro-crack propagation in beam’s mid-region: a) ft,cm = 1.6 MPa (top row), ft,cm = 2.2 MPa (mid-row) and ft,cm = 4.4 MPa (bottom row) in 5 different cross-sections for: (a) depth of 5 mm, (b) depth of 10 mm, (c) depth of 20 mm, (d) depth of 30 mm and (e) depth of 35 mm (red colour - crack, green colour - cement matrix, blue colour aggregate). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) the fracture energy of the cement matrix resulted in the higher beam strength (Fig. 21) without the influence on the macrocrack’s geometry. For spherical particles, the maximum vertical force was slightly higher than for the real angular aggregate (Fig. 22) due to the local stress concentrations at sharp edges of angular grains. Similar conclusions were derived by Trawiński et al. [18], He et al. [64], He [65], Kim and Abu Al-Rub [11] and Du et al. [66] using a mesoscopic continuous FE approach. While manufacturing concrete, its tensile strength may be improved by reducing the porosity in the cement matrix (mainly in the macro-voids and ITZs). This can be done by adding very fine silica particles to the concrete mixture or reducing the water/cement ratio by applying plasticizer. However the concrete brittleness grows at the same time. The concrete Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 22 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx Fig. 21. Influence of fracture energy of cement matrix with respect to GF,ITZ = 20 N/m on calculated vertical force-deflection curve (3D simulations, Va = 35%). Fig. 22. Vertical force-deflection curve from 3D FE simulations for real angular aggregates as compared to averaged response for randomly spaced spherical particles with Va = 35%. tensile strength may be also increased by a reduction of both the thickness of porous ITZs (through applying round smooth aggregate) and number of ITZs (by decreasing the aggregate volumetric content). With respect to the tensile strength, round aggregates are more advantageous than the angular ones and a lower difference between the strength of the cement matrix and ITZs is also more favourable. The concrete brittleness may be diminished by increasing porosity, ITZ-thickness, aggregate volumetric content, aggregate angularity and difference between the strength of the cement matrix and ITZs in order to lengthen the crack propagation way. Thus it is possible by means of comprehensive mesoscopic calculations (by changing several factors at the same time) to elaborate quantitative recommendations for the optimum aggregate volume content, aggregate shape and mortar macro-porosity (in order to simultaneously obtain both the relatively high strength and ductility of concrete). In the next calculation step, both macro-voids and ITZs with a defined width will be taken into account in 3D simulations. The real shape of cement matrix particles will be also considered [53]. The simulation results will be compared with the experimental ones for the different w/c-ratio, porosity, aggregate and fine particle volume and aggregate shape. 7. Conclusions Three-dimensional numerical analyses of concrete fracture at aggregate level in beams under bending brought the following conclusions: – A complex fracture process in notched concrete beams subjected to three-point bending was successfully captured using the 3D meso-scale finite element method (FEM) with cohesive elements embedded in the initial mesh of solid elements. The beam strength was practically insensitive to the accuracy of the mesh describing the aggregate shape. However, the aggregate mesh size slightly affected the geometry of a propagating macro-crack. With a finer aggregate mesh, the macrocrack was slightly more curved. – Considering the real shape and placement of aggregate particles allowed for obtaining a satisfactory agreement between FE outcomes and experimental results in terms of the global (force-displacement) and local (crack pattern) beam Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx 23 response. Three-dimensional numerical results remained in a better agreement with laboratory tests than our previous 2D analyses [18]. – The parametric studies revealed the strong influence of fracture properties of the cement matrix and ITZs. The increase of the tensile strength of both phases contributed to a higher beam’s strength and beam’s brittleness and a smaller crack’s curvature. The effect was more pronounced for intra-phase cohesive elements, since under 3D conditions a crack could not form such a diffuse pattern as in 2D analyses and propagated through a larger number of cement matrix/cement matrix interfaces. The maximum vertical force was higher when increasing fracture energy of the cement matrix and ITZs. – The force-deflection curves for spherical aggregate grains were very similar for the same aggregate volumetric content. However, the aggregate volume fraction influenced the maximum vertical force (the smaller the aggregate volume fraction, the higher was the vertical force). The averaged beam strength was higher for spherical aggregate than for angular one and the stiffness degradation initiated faster for angular aggregate. Fracture simulations with substitute spherical aggregate are able to predict the beam strength but not the crack geometry. The real shape and placement of aggregate have be taken into account in order to obtain a macro-crack path comparable with the experiment. – The simulations of the concrete behaviour at the meso-scale lead to a better design of a concrete mix. The concrete behaviour can be improved by e.g. increasing the mortar porosity (to make the tensile strength higher) and increasing both aggregate volumetric content and aggregate angularity (to make the brittleness lower through a growth of the number of ITZs and their thickness). Acknowledgements The research work has been carried out within the project ‘‘Experimental and numerical analysis of coupled deterministicstatistical size effect in brittle materials” financed by the Polish National Science Centre (NCN) (UMO-2013/09/B/ST8/03598). The FE calculations were carried out at the Academic Computer Centre in Gdańsk (Poland). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] Bažant Z, Planas J. Fracture and size effect in concrete and other quasi-brittle materials. Boca Raton: CRC Press LLC; 1997. Lilliu G, van Mier JGM. 3D lattice type fracture model for concrete. Eng Fract Mech 2003;70:927–41. Tejchman J, Bobiński J. Continuous and discontinuous modelling of fracture in concrete using FEM. Berlin-Heidelberg: Springer; 2013. _ Skarzyński Ł, Nitka M, Tejchman J. Modelling of concrete fracture at aggregate level using FEM and DEM based on X-ray lCT images of internal structure. Eng Fract Mech 2015;147:13–35. Scrivener KL, Crumbie AK, Laugesen P. The interfacial transition zone (ITZ) between cement paste and aggregate in concrete. Interface Sci 2004;12:411–21. Mondal P, Shah SP, Marks LD. Nanomechanical properties of interfacial transition zone in concrete. Nanotechnol Constr 2009;3:315–20. Königsberger M, Pichler B, Hellmich CH. Micromechanics of ITZ-aggregate interaction in concrete Part II: Strength upscaling. J Am Ceram Soc 2014;97:543–51. Gitman IM, Askes H, Sluys LJ. Coupled-volume multi-scale modelling of quasi-brittle material. Eur J Mech A/Solids 2008;27:302–27. _ Ł, Tejchman J. Calculations of fracture process zones on meso-scale in notched concrete beams subjected to three-point bending. Eur J Mech Skarzyński A/Solids 2010;29:746–60. _ Skarzyński Ł, Tejchman J. Modelling the effect of composition on the tensile properties of concrete. Understanding the tensile properties of concrete. Woodhead Publishing Limited; 2013, p. 52–97. Kim SM, Abu Al-Rub RK. Meso-scale computational modelling of the plastic-damage response of cementitious composites. Cem Concr Res 2011;41:339–58. Shahbeyk S, Hosseini M, Yaghoobi M. Meso-scale finite element prediction of concrete failure. Comput Mater Sci 2011;50:1973–90. _ Skarzyński Ł, Tejchman J. Experimental investigations of fracture process in concrete by means of X-ray micro-computed tomography. Strain 2016;52:26–45. Li YY, Schmitt DR. Drilling induced core fractures and in situ stress. J Geophys Res 1998;103:5225–39. Ren W, Yang Z, Sharma R, Zhang Ch, Withers PJ. Two-dimensional X-ray CT image based meso-scale fracture modelling of concrete. Eng Fract Mech 2015;133:24–39. Wang X, Yang Z, Jivkov AP. Monte Carlo simulations of mesoscale fracture of concrete with random aggregates and pores: a size effect study. Constr Build Mater 2015;80:262–72. Wang X, Zhang M, Jivkov AP. Computational technology for analysis of 3D meso-structure effects on damage and failure of concrete. Int J Solids Struct 2016;80:310–33. Trawiński W, Bobiński J, Tejchman J. Two-dimensional simulations of concrete fracture at aggregate level with cohesive elements based on X-ray lCT. Eng Fract Mech 2016;168:204–26. Herrmann HJ, Hansen A, Roux S. Fracture of disordered elastic lattices in two dimensions. Phys Rev B 1989;39:637–47. Schlangen E, Garboczi EJ. Fracture simulations of concrete using lattice models: computational aspects. Eng Fract Mech 1997;57:319–32. Bolander JE, Sukumar N. Irregular lattice model for quasi-static crack propagation. Phys Rev B 2005;71:94–106. Kozicki J, Tejchman J. Effect of aggregate structure on fracture process in concrete using 2D lattice model. Arch Mech 2007;59:1–20. Sakaguchi H, Muhlhaus HB. Mesh free modelling of failure and localization in brittle materials. Deform Prog Fail Geomech 1997:15–21. D’Addetta GA, Kun F, Ramm E. On the application of a discrete model to the fracture process of cohesive granular materials. Granular Matter 2002;4:77–90. Hentz S, Donze FV, Daudeville L. Discrete element modelling of concrete submitted to dynamic loading at high strain rates. Comput Struct 2004;82:2509–24. Nitka M, Tejchman J. Modelling of concrete behaviour in uniaxial compression and tension with DEM. Granular Matter 2015;17:145–64. Carol I, López CM, Roa O. Micromechanical analysis of quasi-brittle materials using fracture-based interface elements. Int J Numer Meth Eng 2001;52:193–215. Caballero A, Carol I, López CM. New results in 3D meso-mechanical analysis of concrete specimens using interface elements. Computational modelling of concrete structures, EURO-C. Taylor and Francis. Group 2006:43–52. Kikuchi A, Kawai T, Suzuki N. The rigid bodies-spring models and their applications to three-dimensional crack problems. Comput Struct 1992;44:469–80. Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003 24 W. Trawiński et al. / Engineering Fracture Mechanics xxx (2017) xxx–xxx [30] Cusatis G, Bažant Z, Cedolin I. Confinement-shear lattice CSL model for fracture propagation in concrete. Comput Methods Appl Mech Eng 2006;195:7154–71. [31] Mariotti C. Lamb’s problem with the lattice model. Geophys J Int 2007;171:857–64. [32] Pamin J. Gradient-dependent plasticity in numerical simulation of localization phenomena. PhD thesis. TU Delft; 1994. [33] de Borst R, Pamin J. Some novel developments in finite element procedures for gradient-dependent plasticity. Int J Numer Meth Eng 1996;39:2477–505. [34] Askes H, Sluys LJ. A classification of higher-order strain gradient models in damage mechanics. Arch Appl Mech 2003;53:448–65. [35] Pamin J. Gradient-enhanced continuum models: formulation, discretization and applications. Habilitation. Cracow University of Technology; 2004. [36] Pijaudier-Cabot G, Bažant ZP. Nonlocal damage theory. ASCE J Eng Mech 1987;113:1512–33. [37] Bažant ZP, Lin F. Non-local yield limit degradation. Int J Numer Meth Eng 1988;35:1805–23. [38] Brinkgreve RBJ. Geomaterial models and numerical analysis of softening. PhD thesis. TU Delft; 1994. [39] Bažant ZP, Jirásek M. Non-local integral formulations of plasticity and damage: survey of progress. J Eng Mech 2002;128:1119–49. [40] Jirásek M, Rolshoven S. Comparison of integral-type nonlocal plasticity models for strain-softening materials. Int J Eng Sci 2003;41:1553–602. [41] Bobiński J, Tejchman J. Comparison of continuous and discontinuous constitutive models to simulate concrete behaviour under mixed mode failure conditions. Int J Numer Anal Meth Geomech 2016;40:406–35. [42] Bobiński J, Tejchman J. A coupled constitutive model for fracture in plain concrete based on continuum theory with non-local softening and eXtended Finite Element Method. Finite Elem Anal Des 2016;114:1–21. [43] http://www.slicer.org. [44] Fedorov A, Beichel R, Kalpathy-Cramer J, Finet J, Fillion-Robin J-C, Pujol S, et al. 3D slicer as an image computing platform for the quantitative imaging network. Magn Reson Imaging 2012;30(9):1323–41. [45] Dessault Systemes. Abaqus documentation, Version 6.13; 2013. [46] Hillerborg A, Modeer M, Petersson PE. Analysis of crack propagation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 1976;6:773–82. [47] Dugdale DS. Yielding of steel sheers containing slits. J Mech Phys Solids 1960;8:100–8. [48] Camacho GT, Ortiz M. Computational modelling of impact damage in brittle materials. Int J Solids Struct 1996;33:2899–938. [49] Tijssens MGA, Sluys BLJ, van der Giessen E. Numerical simulation of quasi-brittle fracture using damaging cohesive surfaces. Eur J Mech A/Solids 2000;19:761–79. [50] Zhou F, Molinari JF. Dynamic crack propagation with cohesive elements: a methodology to address mesh dependency. Int J Numer Meth Eng 2004;59:1–24. [51] de Borst R, Remmers JJC. Computational modelling of delamination. Compos Sci Technol 2006;66:713–22. [52] Yang Z, Xu XF. A heterogeneous cohesive model for quasi-brittle materials considering spatially varying random fracture properties. Comput Methods Appl Mech Eng 2008;197:4027–39. [53] Pichler B, Hellmich Ch. Upscaling quasi-brittle strength of cement paste and mortar: a multi-scale engineering mechanics model. Cem Concr Res 2011;41:467–76. [54] Pichler B, Hellmich Ch, Eberhardsteiner J, Wasserbauer J, Termkhajornkit P, Barbarulo R, et al. Effect of gel-space ratio and microstructure strength of hydrating cementitious materials: an engineering micromechanics approach. Cem Concr Res 2013;45:55–68. [55] Sanahuja J, Dormieux L, Chanvillard G. Modelling elasticity of hydrating cement paste. Cem Concr Res 2007;37:1427–39. [56] Pichler B, Hellmich Ch, Eberhardsteiner J. Spherical and acicular representation of hydrates in a micromechanical model for cement paste: prediction of early-age elasticity and strength. Acta Mech 2009;203:137–62. [57] Tan H, Huang Y, Liu C, Geubelle PH. The Mori-Tanaka method for composite materials with nonlinear interface debonding. Int J Plast 2005;21:1890–918. [58] Qiang TC, Kou XD, Zhou WY. Three-dimensional FEM interface coupled method and its application to arch dam fracture analysis. Chin J Rock Mech Eng 2000;19:562–6. [59] Willam K, Rhee I, Shing B. Interface damage model for thermomechanical degradation of heterogeneous materials. Special WCCMV issue on computational failure mechanics in CMAME; 2003. [60] López CM, Carol I, Aguado A. Meso-structural study of concrete fracture using interface elements. I: Numerical model and tensile behavior. Mater Struct 2008;41:583–99. [61] Su XT, Yang ZJ, Liu GH. Monte Carlo simulation of complex cohesive fracture in random heterogeneous quasi-brittle materials: a 3D study. Int J Solids Struct 2010;47:2336–45. [62] Yang ZJ, Su XT, Chen JF, Liu GH. Monte Carlo simulation of complex cohesive fracture in random heterogeneous quasi-brittle materials. Int J Solids Struct 2009;46:3222–34. [63] http://www.meshmixer.com. [64] He H, Guo Z, Stroeven P, Stroeven M, Sluys LJ. Influence of particle packing on elastic properties of concrete. Proc first international conference on computational technologies in concrete structures (CTCS’09); 2009, p. 1177–97. [65] He H. Computational modeling of particle packing in concrete. PhD thesis. TU Delft; 2010. [66] Du C, Sun L, Jiang S, Ying Z. Numerical simulation of aggregate shapes of three-dimensional concrete and its applications. J Aerosp Eng 2013:515–27. Please cite this article in press as: Trawiński W et al. A three-dimensional meso-scale modelling of concrete fracture, based on cohesive elements and X-ray lCT images. Engng Fract Mech (2017), https://doi.org/10.1016/j.engfracmech.2017.10.003

1/--страниц