The role of particle morphology on concrete fracture behaviour: A meso-scale modelling approach
颗粒形态对混凝土断裂行为的影响：一种中尺度建模方法

1. Introduction 1. 引言

In the context of particulate composites with second phase particles embedded in a host matrix material, the discrepancy of mechanical properties between the two phases can be used to produce superior overall material performance. For example, hard particulate inclusions can strengthen a non-composite matrix material, such as in metal matrices enhanced by ceramic particles [1,2]. In contrast, soft particulates can be distributed in a brittle matrix to enhance the universal fracture toughness [3]. Concrete is a typical example of a matrix-inclusion composite containing hard particulates. Its deformation, damage and failure can occur across the different length scales: macro, meso, micro and nano [4,5].
在含有第二相颗粒嵌入基体材料的颗粒复合材料的背景下，两相之间机械性能的差异可以用来产生优越的整体材料性能。例如，硬质颗粒夹杂物可以增强非复合基体材料，例如通过陶瓷颗粒增强的金属基体 [1,2]。相反，软质颗粒可以分布在脆性基体中，以增强普遍的断裂韧性 [3]。混凝土是一个典型的包含硬质颗粒的基体-夹杂复合材料的例子。其变形、损伤和失效可以在不同的长度尺度上发生：宏观、介观、微观和纳米 [4,5]。

Macroscopically, there is no doubt that excellent progress has been achieved by representing concrete as a homogenous material. However, it is recognised to a certain extent that some mechanical properties, such as the inelastic deformation and brittle fracture behaviour, of heterogenous materials cannot be effectively captured by local continuum mechanics and linear elastic fracture models [6]. Despite the advancement and development of nonlocal continuum mechanics and nonlinear elastic fracture models [[7], [8], [9]], through which mean macro material responses can be described, the spatial variability of local responses is difficult to capture [10]. The first level of concrete inhomogeneity can be found at the meso-scale, where the cementitious composite is a combination of coarse aggregates, interfacial transition zone (ITZ) and mortar paste for normal strength concrete. Fine aggregates and surrounding hardened cement paste hold the mortar paste. High strength concretes typically have low water to cement ratio (e.g. <0.3 by mass) and may incorporate only mortar paste containing fine sands. Both of them, normal or high strength, are called concrete in some studies [73,[96], [97], [98]]. Up to now, most numerical studies meso-scale concrete using explicitly-generated particle shapes focus on normal strength concrete with coarse aggregates and mortar paste [[24], [25], [26], [27], [28], [29],[42], [43], [44],53,65,66]. Due to the computational limitations, requiring similar mesh sizes for coarse and fine aggregates, the mortar paste must be treated as a homogeneous medium, which hinders investigation of fracture nucleation in finer-scale phases within the paste. In this study, we focus on the class of concretes without coarse aggregates, to reveal the interaction and competition of failure mechanisms between the cement paste matrix, fine aggregates and the ITZ, for compression induced fracture behaviour. The numerical framework can be extended for other classes of concretes with coarse aggregates, though at a significantly increased computational cost.
从宏观上看，毫无疑问，通过将混凝土视为均质材料，取得了显著进展。然而，在某种程度上，人们认识到异质材料的一些机械性能，如非弹性变形和脆性断裂行为，无法通过局部连续介质力学和线性弹性断裂模型有效捕捉[6]。尽管非局部连续介质力学和非线性弹性断裂模型[[7]，[8]，[9]]的发展和进步使得可以描述平均宏观材料响应，但局部响应的空间变异性仍然难以捕捉[10]。混凝土非均质性的第一级可以在中观尺度上找到，其中水泥复合材料是粗骨料、界面过渡区（ITZ）和普通强度混凝土的砂浆浆料的组合。细骨料和周围的硬化水泥浆保持砂浆浆料。高强度混凝土通常具有低水胶比（例如，质量比<0.3），并且可能仅包含含有细砂的砂浆浆料。 在某些研究中，正常强度或高强度的混凝土都被称为混凝土[73,[96], [97], [98]]。到目前为止，大多数数值研究中使用显式生成的颗粒形状的中观混凝土主要集中在粗骨料和砂浆的正常强度混凝土上[[24], [25], [26], [27], [28], [29],[42], [43], [44],53,65,66]。由于计算限制，粗骨料和细骨料需要相似的网格尺寸，因此砂浆必须被视为均匀介质，这阻碍了对砂浆中细尺度相的断裂成核的研究。在本研究中，我们关注不含粗骨料的混凝土类别，以揭示水泥浆基体、细骨料和界面过渡区（ITZ）之间的失效机制的相互作用和竞争，以研究压缩引起的断裂行为。尽管数值框架可以扩展到其他含粗骨料的混凝土类别，但计算成本将显著增加。

In the context of the prediction of mechanical properties for cementitious concrete, a major source of discrepancy with conventional composite theories [11], originates from the complexity of the interfacial transition zone (ITZ, e.g., [12,13]). For high strength concrete, the strengths of ITZ and mortar are widely assumed to be the same as aggregates [14], though actual measurement data remains scarce. Due to the experimental difficulties in controlling other factors of inclusions influencing the overall concrete strength (e.g. aggregate fabric, size, morphology and spatial distribution), the solid fraction of sand is widely adopted as the first primitive approximation for a net characterisation of these factors. Accordingly, contrasting findings appear in the literature. For example, when the aggregate volume fraction increased from 45% to 60%, Amparano et al. [15] indicate that the fracture energy and compressive strength would decrease, whereas Tasdemir and Karihaloo [16] conclude the opposite; Guinea et al. [17] found that rough aggregates can decrease the overall strength due to stress concentration at sharp corners, however in Wu et al.'s [14] experiments, 10–20% higher strength could be achieved using crushed quartzite instead of rounded marble aggregates. A plausible explanation for such an observation in Wu et al.'s experiments is that the stress concentration could induce cracks to propagate across harder aggregates when their morphology is rough, and hence more energy would be consumed when compared with fracturing along the relatively weaker ITZ and cement phase, therefore enhancing the overall compressive strength. In a short conclusion, experimental approaches at this scale are complex and difficult to control, particularly when attempting to isolate a single factor from the many other controllable and uncontrollable features.
在水泥混凝土机械性能预测的背景下，与传统复合理论[11]的主要差异来源于界面过渡区（ITZ，参见[12,13]）的复杂性。对于高强度混凝土，ITZ 和砂浆的强度通常被广泛假设与骨料相同[14]，尽管实际测量数据仍然稀缺。由于在控制影响整体混凝土强度的其他因素（如骨料结构、尺寸、形态和空间分布）方面存在实验困难，砂的固体分数被广泛采用作为这些因素净特征的初步近似。因此，文献中出现了相互对立的发现。例如，当骨料体积分数从 45%增加到 60%时，Amparano 等[15]指出，断裂能和抗压强度会降低，而 Tasdemir 和 Karihaloo[16]得出相反的结论；Guinea 等[17]发现，粗糙骨料由于在尖角处的应力集中可能会降低整体强度，然而在 Wu 等的研究中。在[14]的实验中，使用破碎的石英岩代替圆形大理石骨料可以实现 10-20%的强度提高。对于吴等人的实验中观察到的这种现象，一个合理的解释是，压力集中可能会导致裂缝在形态粗糙的硬骨料中扩展，因此与沿相对较弱的界面过渡区和水泥相断裂相比，会消耗更多的能量，从而增强整体的抗压强度。简而言之，这种规模的实验方法复杂且难以控制，特别是在试图从众多可控和不可控特征中孤立出单一因素时。

A complementary tool to investigate the meso-scale fracture behaviour of concrete can be based on computational mechanics. In this scope, many works have explicitly depicted particle shapes using different approaches: molecular dynamics [18,19], DEM [20,21], FEM [22,23], FEM with various fracture models (e.g., cohesive element methods: [24,25]; phase field: [26,27]; and damage plasticity: [28,29]). Given its ability to describe discrete cracks, DEM is a promising approach for the purpose in this paper. Although the classical parallel bond DEM with rigid circular or spherical elements [30] may not be the most suitable for fracture simulation, due to its inherent limitations in describing elasticity and the fact that it requires the fitting of many local model parameters [31], many attempts have been made with DEM due to its novelty in dealing with the search and determination of contacts between fractured surfaces [[32], [33], [34]], considering that friction dissipates most input energy. To bypass the fractured surface contacting behaviour, most continuum-mechanics-based FEM simulations assume that tensile-induced behaviour and fractured surfaces after debonding would not contact intrinsically as encountered in DEM [35]. For meso-scale concrete simulations, where fracture can initiate and propagate in different phases, contact between fractured surfaces is an important feature of the underlying mechanisms.
一种研究混凝土中尺度断裂行为的补充工具可以基于计算力学。在这个范围内，许多研究明确描绘了粒子形状，采用了不同的方法：分子动力学 [18,19]、离散元法（DEM） [20,21]、有限元法（FEM） [22,23]、以及采用各种断裂模型的有限元法（例如，粘结元方法：[24,25]；相场法：[26,27]；和损伤塑性：[28,29]）。鉴于其描述离散裂缝的能力，离散元法（DEM）是本文目的的一个有前景的方法。尽管经典的平行键离散元法使用刚性圆形或球形元素 [30] 可能不是最适合的断裂模拟方法，因为它在描述弹性方面存在固有的局限性，并且需要拟合许多局部模型参数 [31]，但由于其在处理断裂表面之间接触的搜索和确定方面的新颖性，许多尝试已在离散元法中进行 [[32], [33], [34]]，考虑到摩擦消耗了大部分输入能量。 为了绕过破裂表面接触行为，大多数基于连续介质力学的有限元法（FEM）模拟假设拉伸引起的行为和脱粘后的破裂表面不会像离散元法（DEM）中那样本质上接触[35]。对于中尺度混凝土模拟，破裂可以在不同阶段启动和传播，破裂表面之间的接触是潜在机制的重要特征。

Focusing on particles with regular or irregular morphology features, packing (for granular media) or parking (for cementitious materials) models for heterogenous composite materials can be classified into two main categories: real image-based and virtual computer-generated models. The former type of model takes real 2D pixels [27,[36], [37]] or 3D voxels [[38], [39], [40]] as inputs and generates the corresponding computational model. This bottom-up method includes filtering for reducing noise of raw data, a threshold to separate different phases, and mesh generation for numerical methods (e.g. FEM and DEM). Each step of processing contains many parameters and nearly any of them can influence the final computational mechanics model to a certain degree. Further, besides the dependence on image resolution, the meshing of contacting aggregates is still a challenge. This is so because the approximation of the particle geometry would bring unavoidable intersections between contacting particles, which become unrealistic for meso-scale modelling. Some studies have directly implemented voxel data as the FEM mesh. However, after debonding, the contact between non-smooth voxelised surfaces becomes rather inaccurate. Even the approximation for Herztian contact for contacting spheres may not be fully satisfactory [41].
针对具有规则或不规则形态特征的颗粒，异质复合材料的填充（针对颗粒介质）或停放（针对水泥材料）模型可以分为两大类：基于真实图像的模型和虚拟计算机生成的模型。前者模型以真实的二维像素[27,[36], [37]]或三维体素[[38], [39], [40]]作为输入，生成相应的计算模型。这种自下而上的方法包括过滤以减少原始数据的噪声、阈值以分离不同相位，以及用于数值方法（例如有限元法和离散元法）的网格生成。每个处理步骤包含许多参数，几乎任何一个参数都可能在一定程度上影响最终的计算力学模型。此外，除了对图像分辨率的依赖外，接触颗粒的网格化仍然是一个挑战。这是因为颗粒几何形状的近似会导致接触颗粒之间不可避免的交叉，这在中尺度建模中变得不现实。一些研究已直接将体素数据作为有限元网格实现。 然而，在脱粘之后，非光滑体素化表面之间的接触变得相当不准确。即使是对接触球体的赫兹接触的近似也可能并不完全令人满意[41]。

In a top-down fashion, on the contrary, the computer-generated composite model firstly generates single aggregate shapes, and then park or pack them into a predefined domain. As early in 1980s, Zaitsev and Wittmann [42] pioneered the implementation of circular pores and circular or polygonal aggregates in 2D to simulate the fracture nucleation of meso-scale concrete, which attracted a number of subsequent studies [6,10,43,44]. Since then, more types of irregular particle shapes have been covered. Although many irregular shapes have been applied in recent computer-generated models, most are still far from being realistic and are completely convex for ease of contact detection. Considering the complexity of hierarchical particulate geometry, an effective shape parameter is required for describing its universal morphology features in addition to the conventional shape indices (e.g. aspect ratio, sphericity, convexity, roundness and roughness). The fractal dimension (D_f), a cross-scale descriptor that incorporates localised and overall particle morphological features, has the potential for serving as such a parameter [45,46].
相反，计算机生成的复合模型采用自上而下的方式，首先生成单一的聚合物形状，然后将其停放或打包到预定义的区域。早在 1980 年代，Zaitsev 和 Wittmann [42] 首创了在二维中实现圆形孔和圆形或多边形聚合物，以模拟中尺度混凝土的断裂成核，这引发了许多后续研究 [6,10,43,44]。自那时起，更多类型的不规则颗粒形状被涵盖。尽管许多不规则形状已在最近的计算机生成模型中应用，但大多数仍远未逼真，并且为了便于接触检测，完全是凸形的。考虑到分层颗粒几何形状的复杂性，除了传统的形状指标（例如长宽比、球形度、凸度、圆度和粗糙度）外，还需要一个有效的形状参数来描述其普遍的形态特征。分形维数（D）作为一个跨尺度描述符，结合了局部和整体颗粒形态特征，有潜力作为这样的参数 [45,46]。

Based on the above clarification, we are motivated to revisit the meso-scale concrete fracture behaviour with a specific focus on the role of particle morphology. The study is organised as follows. First in Section 2, a combined Spherical Harmonic (SH) analysis and Voronoi tessellation method is proposed to effectively produce realistic particle parking. Then, a combined finite and discrete element method (FDEM) approach is presented, where FEM with contact detection and interaction which is able to simulate not only the continuum behaviour within different phases, but also the initiation and development of internal cracks and the interactions between fractured parts, is adopted to simulate concrete fracture behaviour in Section 3. For the validation of the developed FDEM scheme, in-situ experimental XCT data of a meso-scale concrete sample is imported as a benchmark, including diffraction data for comparing aggregate stress tensors. For better clarity on the meaning and necessity of the simulations with virtual concrete specimen, additional inspirations from the in-situ experiment are also included. In Section 4, results of virtual meso-scale concrete simulations, focusing on various factors, such as aggregate volume fraction, particle morphology and particulate fabric, are discussed. Factors influencing concrete fracture behaviour are studied, highlighting the significant contributions from the particle morphology. Finally, discussions and conclusions are drawn in Section 5.
基于上述阐述，我们有动力重新审视中尺度混凝土的断裂行为，特别关注颗粒形态的作用。研究的组织结构如下。首先，在第二节中，提出了一种结合球谐分析（SH）和 Voronoi 镶嵌方法的方案，以有效生成真实的颗粒停放。然后，在第三节中，提出了一种结合有限元和离散元方法（FDEM）的方案，其中采用具有接触检测和相互作用的有限元方法，能够模拟不同相内的连续行为，以及内部裂缝的产生和发展以及断裂部分之间的相互作用，以模拟混凝土的断裂行为。为了验证所开发的 FDEM 方案，导入了一种中尺度混凝土样本的原位实验 XCT 数据作为基准，包括用于比较骨料应力张量的衍射数据。为了更清晰地理解虚拟混凝土试件模拟的意义和必要性，还包括了来自原位实验的额外启示。 在第 4 节中，讨论了虚拟中尺度混凝土模拟的结果，重点关注各种因素，如骨料体积分数、颗粒形态和颗粒结构。研究了影响混凝土断裂行为的因素，强调了颗粒形态的重要贡献。最后，在第 5 节中进行了讨论和总结。

2. Generation of meso-scale concrete samples
2. 中尺度混凝土样本的生成

In general, for non-overlapping particles via top-down fashion, there are 4 main methods for parking: i) the take-and-place method [24,47]; ii) size scaling method [48]; iii) Tetris method for minimum rectangle or cube particle boundary box [49]; and iv) Graph-based (e.g. Delaunay triangulation in 2D and Voronoi diagram in 3D) shrinking method [50,51]. By keeping the concept of DEM contact search [52] in mind, many dynamic packing methods, where particles are allowed to overlap slightly with one another, have been proposed [22,53]. When such methods are used, a further scaling of particles is needed after packing to eliminate overlapping from meshing. However, both existing parking and packing algorithms are mostly applicable for simplified particle shapes, and with the increase of particle vertices, they become significantly time-consuming [[54], [55], [56]]. On the other hand, due to the limited computational resources, meso-scale simulations must be performed on a representative volume element (RVE). Aggregates across RVE boundaries should be segmented, which is satisfied only by some research works containing simple-shaped inclusions [57,58]. When rough shapes comprising of many fine scale features are cut, a compromise must be made between merging short edges for computational efficiency and the retention of geometrical particle features.
一般来说，对于非重叠粒子采用自上而下的方法，停车主要有四种方法：i) 取放法 [24,47]；ii) 尺寸缩放法 [48]；iii) 最小矩形或立方体粒子边界框的俄罗斯方块法 [49]；iv) 基于图的方法（例如，二维的德劳内三角剖分和三维的沃罗诺伊图）收缩法 [50,51]。在考虑离散元法（DEM）接触搜索的概念 [52] 的基础上，提出了许多动态打包方法，其中粒子可以稍微重叠 [22,53]。使用这些方法时，打包后需要进一步缩放粒子以消除网格中的重叠。然而，现有的停车和打包算法大多适用于简化的粒子形状，随着粒子顶点数量的增加，它们的计算时间显著增加 [[54], [55], [56]]。另一方面，由于计算资源有限，中观尺度的模拟必须在代表性体积单元（RVE）上进行。跨越 RVE 边界的聚集体应被分割，这仅通过一些包含简单形状夹杂物的研究工作得以满足 [57,58]。 当包含许多细小特征的粗糙形状被切割时，必须在合并短边以提高计算效率和保留几何粒子特征之间做出妥协。

Based on above discussions, this section presents the three-main parts of the procedure adopted for the generation of meso-scale concrete samples. This procedure includes, along with the up-down algorithm, three main parts: Voronoi tessellation, Spherical Harmonic (SH) analysis and high-quality surficial mesh generation. First, the concrete domain is tessellated with a Voronoi tessellation, with the resulting cells becoming FEM elements. Compared with relevant research by Mollon and Zhao [55], the faces of the resulting FEM elements have a similar area, thereby avoiding distortions in SH expansion, and every Voronoi cell vertex is FEM node, well preserving cell geometry. For better reconstruction of cell surfaces, more vertices via sub-dividing FEM meshes are covered. Then, an SH analysis is conducted for every meshed cell, of which SH coefficients are subsequently altered to conform to the desired D_f and aspect ratio. The generated single particles are also scaled via their SH coefficients to be totally contained within the associated cell. At last, the SH coefficients are scaled by the given aggregate fraction. High-quality meshes are obtained by reconstructing particle geometry using spherical coordinates (θ, φ) from icosahedron-based geodesic sphere and unit octahedron-based geodesic sphere for complete and cut aggregates, of which the surface is composed of nearly uniform triangles.
基于上述讨论，本节介绍了生成中尺度混凝土样本所采用的程序的三个主要部分。该程序包括上下算法、Voronoi 剖分、球面谐波（SH）分析和高质量表面网格生成。首先，混凝土域通过 Voronoi 剖分进行剖分，生成的单元成为有限元（FEM）元素。与 Mollon 和 Zhao 的相关研究[55]相比，生成的 FEM 元素的面具有相似的面积，从而避免了 SH 展开中的扭曲，并且每个 Voronoi 单元的顶点都是 FEM 节点，良好地保持了单元几何形状。为了更好地重建单元表面，通过细分 FEM 网格覆盖了更多的顶点。然后，对每个网格单元进行 SH 分析，随后调整 SH 系数以符合所需的 D 和长宽比。生成的单个粒子也通过其 SH 系数进行缩放，以完全包含在相关单元内。最后，SH 系数根据给定的骨料比例进行缩放。 高质量网格是通过使用来自基于二十面体的测地球和基于单位八面体的测地球的球坐标（θ，φ）重建粒子几何形状获得的，适用于完整和切割的骨料，其表面由几乎均匀的三角形组成。

2.1. Combined SH analysis and Voronoi tessellation
2.1. 组合 SH 分析与 Voronoi 镶嵌

The distance r_I (x_I(θ, φ), y_I (θ, φ), z_I (θ, φ)) between surface points on a star-shaped particle and the particle centroid in polar coordinate system, can be represented by the orthogonal Spherical Harmonic (SH) function:(1)$r_I\left(\theta ,\phi \right)={\sum }_{n=0}^{\infty }{\sum }_{m=-n}^n{c}_n^m{Y}_n^m\left(\theta ,\phi \right),$(2)$r_I\left(\theta ,\phi \right)=\sqrt{{\left(x_I-x_0\right)}^2+{\left(y_I-y_0\right)}^2+{\left(z_I-z_0\right)}^2},$where I denoted the I-th selected node on the sub-divided cell mesh, (x_I, y_I, z_I) and (x₀, y₀, z₀) are the Cartesian coordinates of I-th node and the mass centre of the intact particle; θ∈[0,π] and φ∈[0,2π] are the latitudinal and longitudinal coordinates respectively, and c_n^m are the SH coefficients of degree n and order m. Y_n^m(θ, φ) (n ∈ ℕ, −n ≤ m ≤ n, ℕ denotes natural number) is the so-called SH function defined on the surface of a sphere as:(3)${\mathrm{Y}}_n^m\left(\theta ,\phi \right)=\sqrt{\frac{\left(2n+1\left)\right(n-\left|m\right|\right)!}{4\pi \left(n+\left|m\right|\right)!}}{P}_n^m\left(cos\theta \right)e^{\mathit{im\phi }},$(4)$Y_n^{-m}\left(\theta ,\phi \right)={\left(-1\right)}^m{\left[\sqrt{\frac{\left(2n+1\left)\right(n-\left|m\right|\right)!}{4\pi \left(n+\left|m\right|\right)!}}{P}_n^m\left(cos\theta \right)e^{\mathit{im\phi }}\right]}^{\ast },$where [.]^⁎ denotes the complex conjugate and P_n^m(x) are called associated Legendre functions, which can be expressed by Rodrigues' formula [59]:(5)$P_n^m\left(x\right)={\left(1-x^2\right)}^{\frac{\left|m\right|}{2}}·\frac{d^{\left|m\right|}}{d{x}^{\left|m\right|}}\left[\frac{1}{2^{n}n!}·\frac{d^n}{d{x}^n}{\left(x^2-1\right)}^n\right].$
在极坐标系中，星形粒子表面点与粒子质心之间的距离 r (x(θ, φ), y(θ, φ), z(θ, φ)) 可以用正交球谐函数 (SH) 表示： (1)$r_I\left(\theta ,\phi \right)={\sum }_{n=0}^{\infty }{\sum }_{m=-n}^n{c}_n^m{Y}_n^m\left(\theta ,\phi \right),$ (2)$r_I\left(\theta ,\phi \right)=\sqrt{{\left(x_I-x_0\right)}^2+{\left(y_I-y_0\right)}^2+{\left(z_I-z_0\right)}^2},$ ，其中 I 表示在细分单元网格上选择的第 I 个节点，(x, y, z) 和 (x0, y0, z0) 分别是第 I 个节点和完整粒子的质心的笛卡尔坐标；θ∈[0,π] 和 φ∈[0,2π] 分别是纬度和经度坐标，cnm 是度 n 和阶 m 的 SH 系数。Ynm(θ,φ) (n ∈ ℕ, −n ≤ m ≤ n, ℕ 表示自然数) 是定义在球面上的所谓 SH 函数，如下所示： (3)${\mathrm{Y}}_n^m\left(\theta ,\phi \right)=\sqrt{\frac{\left(2n+1\left)\right(n-\left|m\right|\right)!}{4\pi \left(n+\left|m\right|\right)!}}{P}_n^m\left(cos\theta \right)e^{\mathit{im\phi }},$ (4)$Y_n^{-m}\left(\theta ,\phi \right)={\left(-1\right)}^m{\left[\sqrt{\frac{\left(2n+1\left)\right(n-\left|m\right|\right)!}{4\pi \left(n+\left|m\right|\right)!}}{P}_n^m\left(cos\theta \right)e^{\mathit{im\phi }}\right]}^{\ast },$ ，其中 [.]⁎ 表示复共轭，Pnm(x) 被称为关联勒让德函数，可以通过罗德里格斯公式 [59] 表示： (5)$P_n^m\left(x\right)={\left(1-x^2\right)}^{\frac{\left|m\right|}{2}}·\frac{d^{\left|m\right|}}{d{x}^{\left|m\right|}}\left[\frac{1}{2^{n}n!}·\frac{d^n}{d{x}^n}{\left(x^2-1\right)}^n\right].$

The derivations and definitions of SH-based fractal dimension (D_f) are included in Appendix A. The relation between SH descriptors and degree n can be bridged by fractal dimension (D_f) as:(6)$D_n=D_2·{\left(\frac{n}{2}\right)}^{2D_f-6},$where D_n is the SH descriptor defined by the normalized amplitude (L_n) of SH coefficient c_n by L₀ in Appendix A, where six kinds of commonly used concrete particulate materials demonstrate the novelty of SH D_f in denoting realistic aggregate morphology. The relative roughness (R_r) quantifying how globally different the irregular aggregate shape is from its c₀-determined sphere can be expressed as (see the derivation in Appendix A)(7)$R_r=\sqrt{\sum _{n=2}^{n_{max}}{\left(D_2·{\left(\frac{n}{2}\right)}^{2D_f-6}\right)}^2},$from which it is clear that such aggregate roughness is dependent on both D₂ and D_f.
SH 基于分形维数 (D) 的推导和定义包含在附录 A 中。SH 描述符与阶数 n 之间的关系可以通过分形维数 (D) 连接，如下所示： (6)$D_n=D_2·{\left(\frac{n}{2}\right)}^{2D_f-6},$ 其中 Dn 是由附录 A 中的 L0 归一化幅度 (Ln) 定义的 SH 描述符 cn，六种常用的混凝土颗粒材料展示了 SH D 在表示真实骨料形态方面的新颖性。相对粗糙度 (Rr) 量化了不规则骨料形状与其 c0 确定的球体之间的全球差异，可以表示为（见附录 A 中的推导） (7)$R_r=\sqrt{\sum _{n=2}^{n_{max}}{\left(D_2·{\left(\frac{n}{2}\right)}^{2D_f-6}\right)}^2},$ ，由此可以清楚地看出，这种骨料粗糙度依赖于 D2 和 D。

Given a domain Ω ∈ R^d, d is the dimension, and it had N Voronoi cell points or generators, {p_i}_i=1^N, with its corresponding Voronoi tessellations or diagrams, {V_i}_i=1^N. The so-called Voronoi tessellation within specific domain can be generated [60]:(8)$V_i=\left\{x|x∊\Omega \wedge ‖x-p_i‖\le ‖x-p_j‖,i\ne j,i,j∊\left[1,N\right]\wedge \mathbb{N}\right\}.$where i and j are i-th and j-th points of N; only if x belongs to the Voronoi cell boundary totally inside the domain boundary, ‖x − p_i‖ = ‖x − p_j‖. Consequently, the whole domain is entirely discretised:(9)$\left\{\begin{array}{c}\bigcup _{i=1}^N{V}_i=\Omega \\ V_i\cap V_j=\varnothing \end{array}\right..$
给定一个域 Ω ∈ Rd，d 是维度，并且它有 N 个 Voronoi 单元点或生成点 {p}i=1N，以及相应的 Voronoi 划分或图 {V}i=1N。所谓的 Voronoi 划分可以在特定域内生成 [60]： (8)$V_i=\left\{x|x∊\Omega \wedge \left(\left.x-p_i\right)\le \left(\left.x-p_j\right)\right.\right.,i\ne j,i,j∊\left[1,N\right]\wedge \mathbb{N}\right\}.$ 其中 i 和 j 是 N 的第 i 和第 j 个点；仅当 x 完全属于域边界内的 Voronoi 单元边界时，‖x − p‖ = ‖x − p‖。因此，整个域被完全离散化： (9)$\left\{\begin{array}{c}\bigcup _{i=1}^N{V}_i=\Omega \\ V_i\cap V_j=\varnothing \end{array}\right..$

If the mass density of the domain or Voronoi tessellations is ρ(x), the concurrent mass centre {m_i}_i=1^N of single cell with its volume equal to v_i is(10)$m_i=\frac{1}{v_i}{\int }_{v_i}x·\rho \left(x\right)\mathit{dx},i∊\left[1,N\right]\wedge \mathbb{N}.$
如果区域或 Voronoi 镶嵌的质量密度为ρ(x)，则单个体积为 v 的单元的并发质心{m}i=1N 为 (10)$m_i=\frac{1}{v_i}{\int }_{v_i}x·\rho \left(x\right)\mathit{dx},i∊\left[1,N\right]\wedge \mathbb{N}.$

To introduce the differences between cells of a domain, a geometry energy function is defined as:(11)$\mathbb{E}\left({\left\{p_i\right\}}_{i=1}^N,{\left\{V_i\right\}}_{i=1}^N\right)=\sum _{i=1}^N{\int }_{v_i}‖x-p_i‖·\rho \left(x\right)\mathit{dx}.$
为了引入一个领域中细胞之间的差异，定义了一个几何能量函数为： (11)$\mathbb{E}\left({\left\{p_i\right\}}_{i=1}^N,{\left\{V_i\right\}}_{i=1}^N\right)=\sum _{i=1}^N{\int }_{v_i}\left(x-p_i\right)·\rho \left(x\right)\mathit{dx}.$

Evidently, when {p_i}_i=1^N coincides with {m_i}_i=1^N the energy is the least and the cell geometry is the most stable. However, these two points are usually far from each other, and a perfect coincidence is rather difficult to be obtained. Iterations are applied to the easy case where {p_i}_i=1^N is randomly distributed in the domain via moving generators of next loop to the mass centre of current cells:(12)${\left\{{\left\{p_i\right\}}_{i=1}^N\right\}}_{\mathrm{i}+1}={\left\{{\left\{m_i\right\}}_{i=1}^N\right\}}_{\mathrm{i}}.$where ⅈ means ⅈ-th loop or iterations. In this study, aggregate shapes are generated by the shrinking of Voronoi cells, via SH analysis, and are statistically similar, though not yet identical, to each other. With the given D₂ and D_f, the process to obtain temporary SH coefficients ${\left\{{c^{\prime \prime }}_n\right\}}_{i=1}^N$ of realistic aggregate shapes from those ({c′_{i, n}}_i=1^N) of Voronoi cells is given in Appendix B. To guarantee the aggregate to be completely contained in the cell, one can have:(13)$\forall {\left\{r_I\left(\theta ,\phi \right)\right\}}_{i=1}^N\le {\left\{d_I\left(\theta ,\phi \right)\right\}}_{i=1}^N$where d_I(θ, φ) is the distance between cell surface vertices and its mass centre.
显然，当 {p}i=1N 与 {m}i=1N 重合时，能量最低，细胞几何形状最稳定。然而，这两个点通常相距较远，完美重合相当困难。迭代应用于简单情况，其中 {p}i=1N 在域内随机分布，通过将下一个循环的生成器移动到当前细胞的质心： (12)${\left\{{\left\{p_i\right\}}_{i=1}^N\right\}}_{\mathrm{i}+1}={\left\{{\left\{m_i\right\}}_{i=1}^N\right\}}_{\mathrm{i}}.$ ，其中 ⅈ 表示 ⅈ-th 循环或迭代。在本研究中，通过 SH 分析生成的聚合物形状是通过 Voronoi 细胞的收缩而形成的，虽然在统计上相似，但尚未完全相同。给定 D2 和 D，从 Voronoi 细胞的 ({c′i, n}i=1N) 中获得现实聚合物形状的临时 SH 系数 ${\left\{{c^{\prime \prime }}_n\right\}}_{i=1}^N$ 的过程在附录 B 中给出。为了确保聚合物完全包含在细胞内，可以得到： (13)$\forall {\left\{r_I\left(\theta ,\phi \right)\right\}}_{i=1}^N\le {\left\{d_I\left(\theta ,\phi \right)\right\}}_{i=1}^N$ ，其中 d(θ,φ) 是细胞表面顶点与其质心之间的距离。

The scaling from mass centre for SH coefficients {c_{i, n}}_i=1^N of final aggregate shapes must be done:(14)${\left\{c_{i,n}\right\}}_{i=1}^N\le \frac{{\left\{{c^{\prime \prime }}_n\right\}}_{i=1}^N}{max\left({\left\{\frac{r_I\left(\theta ,\phi \right)}{d_I\left(\theta ,\phi \right)}\right\}}_{i=1}^N\right)}.$
最终聚合物形状的 SH 系数{ci, n}i=1N 必须从质心进行缩放： (14)${\left\{c_{i,n}\right\}}_{i=1}^N\le \frac{{\left\{{c^{\prime \prime }}_n\right\}}_{i=1}^N}{max\left({\left\{\frac{r_I\left(\theta ,\phi \right)}{d_I\left(\theta ,\phi \right)}\right\}}_{i=1}^N\right)}.$

For further improving the efficiency of proposed framework in producing realistic particle shapes, five shape parameters are defined as:(15)$\mathcal{F}=\frac{\mathit{III}}{\mathit{II}};$(16)$\mathcal{E}=\frac{\mathit{II}}{I};$(17)$\mathcal{S}=\frac{\sqrt[3]{36\pi {{\mathbb{V}}_a}^2}}{\mathbb{S}};$(18)$\mathcal{C}=\frac{{\mathbb{V}}_a}{{\mathbb{V}}_c};$(19)$\mathcal{R}=\frac{\sum \left({\mathbb{S}}_l·\frac{k_{\mathit{in}}}{k_{M,l}}\right)}{\sum {\mathbb{S}}_l},\forall k_{M,l}\le k_{\mathit{in}},$where$\mathcal{F}$,$\mathcal{E}$,$\mathcal{S}$, $\mathcal{C}$ and $\mathcal{R}$ are shape indices of flatness, elongation, sphericity, convexity and roundness; I, II and III (I > II > III) denote the lengths of the longest, intermediate and shortest dimensions of the box exactly containing the aggregate using the Principal Component Analysis (PCA) in Matlab environment [94], respectively; $\mathbb{S}$ and ${\mathbb{S}}_l$ are the surface area of the aggregate and l-th triangular mesh to compose its surface; ${\mathbb{V}}_a$ and ${\mathbb{V}}_c$ are the volumes of the aggregate and its perfect convex hull; k_in is curvature value of the maximum inscribed sphere, while k_{M, l} is the average of median curvature values of l-th triangle's three vertices [61].
为了进一步提高所提框架在生成真实粒子形状方面的效率，定义了五个形状参数： (15)$\mathcal{F}=\frac{\mathit{III}}{\mathit{II}};$ (16)$\mathcal{E}=\frac{\mathit{II}}{I};$ (17)$\mathcal{S}=\frac{\sqrt[3]{36\pi {{\mathbb{V}}_a}^2}}{\mathbb{S}};$ (18)$\mathcal{C}=\frac{{\mathbb{V}}_a}{{\mathbb{V}}_c};$ (19)$\mathcal{R}=\frac{\sum \left({\mathbb{S}}_l·\frac{k_{\mathit{in}}}{k_{M,l}}\right)}{\sum {\mathbb{S}}_l},\forall k_{M,l}\le k_{\mathit{in}},$ ，其中 $\mathcal{F}$ 、 $\mathcal{E}$ 、 $\mathcal{S}$ 、 $\mathcal{C}$ 和 $\mathcal{R}$ 分别是平坦度、延伸度、球形度、凸度和圆度的形状指标；I、II 和 III（I > II > III）表示使用 Matlab 环境中的主成分分析（PCA）精确包含聚合物的盒子的最长、中间和最短维度的长度； $\mathbb{S}$ 和 ${\mathbb{S}}_l$ 分别是聚合物的表面积和第 l 个三角网格的表面积； ${\mathbb{V}}_a$ 和 ${\mathbb{V}}_c$ 是聚合物及其完美凸包的体积；kin 是最大内切球的曲率值，而 kM，l 是第 l 个三角形三个顶点的中位曲率值的平均值[61]。

A domain composed of 200 Voronoi cells is generated, of which the parental particles are two kinds of real sands, one coarse sand (MA-1117) and one fine sand (MA-106A) in Fig. A.1 of Appendix A. Fig. 1 illustrates the morphological features of the real and virtual particles, of which shape indices of virtual particulates are compared with their corresponding virtual samples in Fig. 2. The shape parameters of real aggregates in Fig. 2 are from an open-source software, Virtual Cement and Concrete Testing Library (VCCTL, Bullard, 2014). We use the Spherical Harmonic (SH) coefficients of various kinds of aggregates in VCCTL to conduct SH expansion in Eq. (1) and icosahedral geodesic structures with 1280 faces to reconstruct their shapes. Then, the parameters of the virtual aggregates are calculated according to Eqs. (15), (16), (17), (18), (19). D_f and D₂ conform to the fitted normal distributions of those from parental ones. Since the packing and generation of realistic-shaped aggregates are obtained from the shrinking process of Voronoi cells, the aggregate shapes are somewhat dependent on their corresponding cells, although irregular morphology features can be imposed via Spherical Harmonic (SH) analysis. According to properties of centroidal Voronoi tessellations [60], increasing iterations (moving generators of the next loop to the mass centre of current cells) makes the shapes and sizes of Voronoi cells more uniform. For example, if a large enough number of iterations are made for a cubic domain containing L³ generators, L³ congruent cubic cells would be composing the whole region. Hence, the generated aggregates are more realistic than after 100 iterations in their spatial distributions, size distributions and shape parameters. One may expect the number of iterations for original Voronoi generators would much influence the shape parameters of the final particle shapes. However, from our quantitative results of shape indices of 200 particles in Fig. 2 and snapshots of virtual particle morphology in Fig. 1(b) to (c) and (e) to (f), such an influence is not too evident and a convergence tendency of particle shape indices appears with the increase of iterations of the initial Voronoi tessellations. Similar results, where Voronoi shapes have limited influence on shape parameters of their correspondingly realistic particles, have also been reported by [55], who first produce realistic-shaped particles and then put them into Voronoi cells. For coarse sands, clear angular faces can be seen from Fig. 1(d) while in Fig. 1(a) the fine sands have more curved surfaces. These key features also occurred in virtual samples from the high normalized curvature value ($\frac{k_{M,I}}{k_{\mathit{in}}}$) in Fig. 1(e) and (f). Notably, to facilitate the study on the effects of D₂ and D_f in concrete fracture behaviour, each aggregate had the same D₂ and D_f values in the following parts for the virtual concrete sample.
生成了一个由 200 个 Voronoi 单元组成的区域，其中母粒子为两种真实沙子，一种粗砂（MA-1117）和一种细砂（MA-106A），见附录 A 的图 A.1。图 1 展示了真实粒子和虚拟粒子的形态特征，虚拟颗粒的形状指数与其对应的虚拟样本在图 2 中进行了比较。图 2 中真实集料的形状参数来自开源软件虚拟水泥和混凝土测试库（VCCTL，Bullard，2014）。我们使用 VCCTL 中各种集料的球面谐波（SH）系数在公式（1）中进行 SH 展开，并使用 1280 面体的二十面体地理结构重建它们的形状。然后，根据公式（15）、（16）、（17）、（18）、（19）计算虚拟集料的参数。D 和 D2 符合来自母粒子的拟合正态分布。 由于真实形状的聚合物的打包和生成是通过 Voronoi 单元的收缩过程获得的，因此聚合物的形状在一定程度上依赖于其对应的单元，尽管可以通过球面谐波（SH）分析施加不规则的形态特征。根据质心 Voronoi 镶嵌的特性[60]，增加迭代次数（将下一个循环的生成器移动到当前单元的质心）使 Voronoi 单元的形状和大小更加均匀。例如，如果对包含 L3 个生成器的立方体域进行足够多的迭代，L3 个全等的立方体单元将组成整个区域。因此，生成的聚合物在空间分布、大小分布和形状参数上比经过 100 次迭代后的聚合物更为真实。人们可能会预期原始 Voronoi 生成器的迭代次数会对最终颗粒形状的形状参数产生很大影响。然而，从我们在图 2 中对 200 个颗粒的形状指数的定量结果以及图中的虚拟颗粒形态快照来看， 1(b) 到 (c) 和 (e) 到 (f) 的影响并不明显，随着初始 Voronoi 划分迭代次数的增加，颗粒形状指标出现了趋同的趋势。类似的结果也在 [55] 中报道，该研究首先生成了真实形状的颗粒，然后将其放入 Voronoi 单元中。对于粗砂，从图 1(d) 可以清晰地看到棱角分明的面，而在图 1(a) 中，细砂则具有更为弯曲的表面。这些关键特征也出现在图 1(e) 和 (f) 中高归一化曲率值 ( $\frac{k_{M,I}}{k_{\mathit{in}}}$ ) 的虚拟样本中。值得注意的是，为了便于研究 D2 和 D 对混凝土断裂行为的影响，以下部分的虚拟混凝土样本中每种骨料的 D2 和 D 值均相同。

1. Download: Download high-res image (664KB)
   下载：下载高分辨率图像（664KB）
2. Download: Download full-size image
   下载：下载完整尺寸图像

Fig. 1. (a) and (d): Snapshots of real fine sand (MA-106A) and coarse sand (MA-1117) after [62], of which mean size is about 1 mm and 3 cm; (b) and (c): virtual fine sand packings from Voronoi cells after 0 and 100 iterations; (e) and (f): virtual coarse sand packings from Voronoi cells after 0 and 100 iterations. The colour bar denotes mean curvature value of vertices normalized by that of its correspondingly maximum sphere ($\frac{k_{M,I}}{k_{\mathit{in}}}$). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
图 1. (a)和(d)：真实细沙（MA-106A）和粗沙（MA-1117）的快照，来源于[62]，其平均粒径约为 1 毫米和 3 厘米；(b)和(c)：经过 0 次和 100 次迭代的 Voronoi 单元的虚拟细沙堆积；(e)和(f)：经过 0 次和 100 次迭代的 Voronoi 单元的虚拟粗沙堆积。颜色条表示顶点的平均曲率值，经过其对应最大球体的归一化处理（ $\frac{k_{M,I}}{k_{\mathit{in}}}$ ）。(关于本图例中颜色的解释，读者可参考本文的网络版本。)

1. Download: Download high-res image (286KB)
   下载：下载高分辨率图像（286KB）
2. Download: Download full-size image
   下载：下载完整尺寸图像

Fig. 2. Cumulative distribution functions (CDF) of various shape parameters for real and virtual fine ((a)–(e)) and coarse aggregates ((f)–(j)). $\mathcal{F}$,$\mathcal{E}$,$\mathcal{S}$, $\mathcal{C}$ and $\mathcal{R}$ denote flatness, elongation, sphericity, convexity and roundness.
图 2. 实际和虚拟细骨料((a)–(e))及粗骨料((f)–(j))的各种形状参数的累积分布函数(CDF)。 $\mathcal{F}$ 、 $\mathcal{E}$ 、 $\mathcal{S}$ 、 $\mathcal{C}$ 和 $\mathcal{R}$ 分别表示平坦度、延伸度、球形度、凸度和圆度。

Voronoi tessellation has boundary effects since the cells connecting the domain boundary must have the boundary shapes of the domain [63]. To this end, as shown in Fig. 3(a), a larger periodic domain with size equal to 3 L × 3 L × 3 L is defined for its uniformly distributed sub-boxes, each of size L × L × L. Ten iterations are first conducted in a sub-box having 64 random generators, which is more than the 24 aggregates for the RVE in composite theory [57]. The distributions of the Voronoi cell points are then copied into other 26 sub-domains. Finally, Voronoi cells are generated in the larger domain. Consequently, to avoid lengthy computation only cells including complete or cut aggregates are selected in Fig. 3(b), of which the entire concrete boundary is shown in Fig. 3(c).
Voronoi 划分存在边界效应，因为连接域边界的单元必须具有域的边界形状 [63]。为此，如图 3(a) 所示，定义了一个大小为 3 L × 3 L × 3 L 的更大周期性域，以便其均匀分布的子盒，每个子盒的大小为 L × L × L。首先在一个具有 64 个随机生成器的子盒中进行十次迭代，这个数量超过了复合材料理论中 RVE 的 24 个聚集体 [57]。然后，Voronoi 单元点的分布被复制到其他 26 个子域中。最后，在更大的域中生成 Voronoi 单元。因此，为了避免冗长的计算，仅选择图 3(b) 中包含完整或切割聚集体的单元，其整个混凝土边界如图 3(c) 所示。

1. Download: Download high-res image (310KB)
   下载：下载高分辨率图像（310KB）
2. Download: Download full-size image
   下载：下载完整尺寸图像

Fig. 3. Periodic Voronoi domain with 27 sub-domains (a), the enlargement of interest cells containing complete and cut aggregates (b) and cut Voronoi cells in the domain of concrete (c).
图 3. 具有 27 个子域的周期性 Voronoi 域(a)，包含完整和切割聚集体的感兴趣单元的放大(b)以及混凝土域中的切割 Voronoi 单元(c)。

2.2. Surface mesh generation
2.2. 表面网格生成

High-quality FEM mesh generation is still a big challenge to successfully simulate meso-scale concrete fracture behaviour, especially for realistically shaped aggregates, for which four neighboured vertices do not fall on the same flat plane. This may be the reason why voxel data, after being re-scaled and resolution sacrificed, is implemented as a FEM cubic mesh for most XCT-based simulations [29,64,65]. Some researches [66,67] also generate a tetrahedral aggregate volume mesh via a mapping into an unstructured FEM mesh of the whole concrete. However, the resulting aggregate surface tomography is dependent on the pre-defined mesh of concrete, which is not sufficient for incorporating the influence of the aggregate morphology features. For example, even though diverse morphological details have been given to aggregates in terms of initial surficial meshes, the morphology would be transformed after mapping to that of pre-defined volume meshes. In this study, the generated aggregate mesh is first generated before the determination of cement paste matrix by cutting the aggregate mesh from the entire concrete domain.
高质量的有限元网格生成仍然是成功模拟中尺度混凝土断裂行为的一大挑战，特别是对于形状真实的骨料，其四个相邻顶点并不位于同一平面上。这可能是为什么在大多数基于 XCT 的模拟中，体素数据在重新缩放和牺牲分辨率后被用作有限元立方体网格的原因 [29,64,65]。一些研究 [66,67] 还通过映射到整个混凝土的非结构化有限元网格生成四面体骨料体积网格。然而，生成的骨料表面断层成像依赖于预定义的混凝土网格，这不足以考虑骨料形态特征的影响。例如，即使在初始表面网格方面给骨料赋予了多样的形态细节，映射到预定义体积网格后，形态也会发生变化。在本研究中，生成的骨料网格首先在确定水泥浆基体之前生成，通过从整个混凝土域中切割骨料网格。

One main advantage of SH-generated surfaces is the ability to reconstruct particle shapes based on the implicit form of aggregate boundary, as long as its SH coefficient is known in Eq. (1). On the other hand, to bypass the computational burden of FEM, an appropriate number of high-quality surficial triangular meshes is needed. To mesh implicit boundary with SH function, many researches have implemented icosahedron-based geodesic spheres for this purpose [55,68], of which nodes are projected to those of implicit boundary with the same spherical angles. However, the gap between the number (20×4^ℵ, ℵ ∊ ℕ) of meshes for this type of spherical angles is very long. Here, we approximate the arbitrary number of uniform triangular tessellating the unit sphere surface via searching for the minimum electrostatic energy [69]:(20)$\mathrm{Є}=\sum _{j=1}^{\mathrm{\aleph }}\sum _{k=1}^{\mathrm{\aleph }}\frac{q_j·q_k}{\sqrt{d_{j,k}}},j\ne k,$where ℵ is the number of vertices on the sphere; q is the product of length and area distortion or the so-called particle charge; j and k denote j- and k-th vertex; d is the geodesic distance between two vertices. Appendix C Details of surficial mesh generation, Appendix D Approximation of introduce surficial mesh generation of complete aggregates, as well as cut ones using octahedron-based geodesic structure. Fig. 4 shows the final surficial mesh of aggregates in the concrete sample. Clearly there are no short edges due to the cutting process, and the mesh sizes of uncut aggregates are nearly identical to each other.
SH 生成表面的一个主要优点是能够根据聚合边界的隐式形式重建粒子形状，只要其在公式(1)中的 SH 系数已知。另一方面，为了绕过有限元法的计算负担，需要适当数量的高质量表面三角网格。为了用 SH 函数对隐式边界进行网格划分，许多研究已实现基于二十面体的测地球面[55,68]，其节点被投影到具有相同球面角度的隐式边界上。然而，这种球面角度的网格数量（20×4ℵ，ℵ ∊ ℕ）之间的差距非常大。在这里，我们通过寻找最小静电能量来近似均匀三角形镶嵌单位球面上的任意数量[69]： (20)$\mathrm{Є}=\sum _{j=1}^{\mathrm{\aleph }}\sum _{k=1}^{\mathrm{\aleph }}\frac{q_j·q_k}{\sqrt{d_{j,k}}},j\ne k,$ ，其中ℵ是球面上的顶点数量；q 是长度和面积畸变的乘积或所谓的粒子电荷；j 和 k 分别表示第 j 和第 k 个顶点；d 是两个顶点之间的测地距离。 附录 C 表面网格生成的详细信息，附录 D 完整骨料的表面网格生成近似，以及使用八面体基础的测地结构对切割骨料的表面网格生成。图 4 显示了混凝土样本中骨料的最终表面网格。显然，由于切割过程，没有短边，并且未切割骨料的网格大小几乎相同。

1. Download: Download high-res image (176KB)
   下载：下载高分辨率图像（176KB）
2. Download: Download full-size image
   下载：下载完整尺寸图像

Fig. 4. Final surficial meshes of aggregates in concrete sample with size equal to 1 mm × 1 mm × 1 mm.
图 4. 混凝土样本中尺寸为 1 mm × 1 mm × 1 mm 的骨料最终表面网格。

3. FDEM simulation and its validation
3. FDEM 模拟及其验证

In this section the combined finite and discrete element method (FDEM) is briefly reviewed [70,71] and the key details are summarised. For more information please refer to a previous paper for single sand crushing behaviour [72]. The combined XCT and diffraction data from a uniaxial compression experiment on concrete [73] is also used as a benchmark for the simulations.
在本节中，简要回顾了组合有限和离散元方法（FDEM）[70,71]，并总结了关键细节。有关更多信息，请参阅之前关于单砂破碎行为的论文[72]。来自混凝土单轴压缩实验的组合 XCT 和衍射数据[73]也被用作模拟的基准。

3.1. FDEM and its damage model
3.1. FDEM 及其损伤模型

Within FDEM, solid parts are discretised using 4-noded tetrahedral FEM elements and the deformation and stress in them can be captured as in FEM, as shown in Fig. 5(a). Simultaneously, Fig. 5(b) shows the capability of DEM to deal with contact, thus FDEM can efficiently model the transition of continuum to discrete states, including the presence of fracture or even breakage of concrete. Although most aggregates undergo fracture, they can still bear a significant portion of the macro load due to their high strength, see Fig. 5(c). An explicit solver via central difference integration framework is adopted, due to its higher efficiency in dealing with a large number of contacts and deformation when compared with implicit solvers. A penalty function is implemented to calculate contact forces and prevent the overlapping of contacting elements. Viscous damping, together with normal elastic response, is also introduced for numerical stability, while a Coulomb-type friction law is employed for friction forces. Fracture is simulated via the debonding of 6-node cohesive interfacial elements (CIEs) inserted between every two conjunct 4-node solid elements. In addition to solid aggregate elements used to simulate fracture of three-phase composed of concrete, another two kinds of CIEs and two types of solid elements are needed, see the experimental sample in Fig. 6 and corresponding solid aggregate elements in Fig. 7(d). After the deletion of the CIEs that have failed, completely fractured surfaces can still develop contact via the penalty function. The material failure process had two stages: fracture initiation and propagation. The maximum stress criterion, including tensile and shear strengths, is applied for crack initiation. Bilinear traction-separation damage laws are employed for the stiffness degradation during crack evolution. More details about the damage models are provided in Appendix E. For mixed damage-mode opening, the Benzeggagh-Kenane criterion [74] is employed:(21)$\left\{\begin{array}{c}{G}_C=G_{n,C}+\left(G_{s,C}-G_{n,C}\right){\left(\frac{G_{\text{shear}}}{G_T}\right)}^{\eta }\\ G_{\text{shear}}=G_s+G_t\\ G_T=G_n+G_s+G_t\end{array}\right.,$where G_n,C and G_shear are the respective critical energies for pure mode I (tension) and combined modes II and III (shear), η is a material parameter, and G_T is the mixed fracture energy for the CIEs. For simplicity, we set G_{s, C} = G_{n, C} and thus only G_{n, C} influenced the energy-based crack evolution criterion.
在 FDEM 中，固体部分使用 4 节点四面体有限元（FEM）单元进行离散化，其变形和应力可以像在 FEM 中一样被捕捉，如图 5(a)所示。同时，图 5(b)展示了 DEM 处理接触的能力，因此 FDEM 可以有效地模拟连续体到离散状态的过渡，包括混凝土的裂缝甚至破裂的存在。尽管大多数骨料会发生断裂，但由于其高强度，它们仍然能够承受相当一部分的宏观载荷，见图 5(c)。采用显式求解器通过中心差分积分框架，由于其在处理大量接触和变形时相较于隐式求解器具有更高的效率。实施了惩罚函数以计算接触力并防止接触元素的重叠。引入了粘性阻尼以及正常弹性响应以确保数值稳定性，同时采用了库仑型摩擦定律来处理摩擦力。通过在每两个相邻的 4 节点固体元素之间插入 6 节点粘结界面元素（CIEs）的脱粘来模拟断裂。 除了用于模拟由混凝土组成的三相断裂的固体骨料元素外，还需要另外两种 CIE 和两种类型的固体元素，参见图 6 中的实验样本和图 7(d)中的相应固体骨料元素。在删除已失效的 CIE 后，完全断裂的表面仍然可以通过惩罚函数发展接触。材料失效过程分为两个阶段：断裂起始和扩展。最大应力准则，包括拉伸和剪切强度，用于裂纹的起始。双线性拉伸-分离损伤法则用于裂纹演化过程中的刚度退化。关于损伤模型的更多细节见附录 E。对于混合损伤模式开启，采用 Benzeggagh-Kenane 准则 (21)$\left\{\begin{array}{c}{G}_C=G_{n,C}+\left(G_{s,C}-G_{n,C}\right){\left(\frac{G_{\text{shear}}}{G_T}\right)}^{\eta }\\ G_{\text{shear}}=G_s+G_t\\ G_T=G_n+G_s+G_t\end{array}\right.,$ ，其中 Gn,C 和 Gshear 分别是纯模式 I（拉伸）和组合模式 II 和 III（剪切）的临界能量，η是材料参数，GT 是 CIE 的混合断裂能量。为简化起见，我们设定 Gs,C = Gn,C，因此只有 Gn,C 影响基于能量的裂纹演化准则。

1. Download: Download high-res image (415KB)
   下载：下载高分辨率图像（415KB）
2. Download: Download full-size image
   下载：下载完整尺寸图像

Fig. 5. (a) FEM meshes with displacement boundary conditions at top and bottom layer nodes, (b) fracture pattern and (c) stress contour of one typical virtual meso-scale concrete.
图 5. (a) 顶层和底层节点的位移边界条件的有限元网格，(b) 断裂模式，(c) 一种典型虚拟中尺度混凝土的应力轮廓。

1. Download: Download high-res image (323KB)
   下载：下载高分辨率图像（323KB）
2. Download: Download full-size image
   下载：下载完整尺寸图像

Fig. 6. Different FDEM elements of the experimental sample of which aggregate elements are in Fig. 7(d).
图 6. 实验样本的不同 FDEM 元素，其聚合元素见图 7(d)。

1. Download: Download high-res image (550KB)
   下载：下载高分辨率图像（550KB）
2. Download: Download full-size image
   下载：下载完整尺寸图像

Fig. 7. Image processing from XCT images to FEM mesh.
图 7. 从 XCT 图像到 FEM 网格的图像处理。

Table 1 summarise the mechanical parameters adopted in this study. Note that the solid properties are taken from experimental researches by [75,76]; the stiffness values and mesh size are determined following [25,77]; other mechanical properties of CIEs are taken directly from relevant research by [77]; and the friction coefficient was 0.5.
表 1 总结了本研究中采用的机械参数。请注意，固体属性取自[75,76]的实验研究；刚度值和网格大小根据[25,77]确定；CIE 的其他机械属性直接取自[77]的相关研究；摩擦系数为 0.5。

Table 1. FDEM material parameters of elements in concrete.
表 1. 混凝土中元素的 FDEM 材料参数。

Empty Cell	Parameter 参数	Value 价值
Solid elements 固体元素	Density, ρ (kg/m3) 密度，ρ (kg/m³)	Agg: 2650
		CP: 2200
	Elastic modulus, E (GPa) 弹性模量，E (GPa)	Agg: 96 聚合物：96
		CP: 16
	Poisson's ratio, υ 泊松比，υ	Agg: 0.08
		CP: 0.23
Cohesive interface elements (CIEs) 凝聚性界面元素 (CIEs)	Normal stiffness, kn (N/m3) 正常刚度，kn (N/m³)	Agg: 9.6 × 1013 聚合物：9.6 × 10^13
		CP: 2.5 × 1013 CP: 2.5 × 10^13
		ITZ: 1.36 × 1013 ITZ: 1.36 × 10¹³
	Shear stiffness, ks (N/m3) 剪切刚度，ks (N/m³)	Agg: 4.44 × 1013 聚合物：4.44 × 10^13
		CP: 1.04 × 1013 CP: 1.04 × 10¹³
		ITZ: 5.53 × 1013 ITZ: 5.53 × 10^13
	Tensile strength, Nmax (MPa) 抗拉强度，Nmax (MPa)	Agg: 10 聚合物：10
		CP: 4
		ITZ: 2.4 ITZ：2.4
	Shear strengths, Smax and Tmax (MPa) 剪切强度，Smax 和 Tmax (MPa)	Agg: 5 聚合物：5
		CP: 2
		ITZ: 1.2 ITZ：1.2
	Fracture energy, GC (N/m) 断裂能，GC (N/m)	Agg: 60 聚合度：60
		CP: 50
		ITZ: 30 ITZ：30
Contact law 接触法	Friction coefficient, μ 摩擦系数，μ	0.5

Quasi-static loading conditions are approximated in the simulations by keeping kinetic energy below 10% of the associated internal energy. In this study, no loading platens are explicitly constructed, while the boundary conditions are set on the nodes at the top and bottom surfaces of the hexahedron concrete along loading direction, as in Fig. 5(a). The bottom surface is fixed in the loading direction, while the nodes at the top surface are applied a downward linearly ramp displacement velocity up to 0.04 m/s, which is then kept constant until an overall displacement of 0.1 mm was met. Nodes at both bottom and top surfaces are allowed to expand in other directions than loading direction.
在模拟中，通过将动能保持在相关内能的 10%以下来近似准静态加载条件。在本研究中，没有明确构建加载平板，而是在六面体混凝土的顶部和底部表面节点上设置边界条件，如图 5(a)所示。底部表面在加载方向上固定，而顶部表面的节点施加向下的线性 ramp 位移速度，最高可达 0.04 m/s，随后保持不变，直到整体位移达到 0.1 mm。底部和顶部表面的节点允许在加载方向以外的其他方向扩展。

3.2. In-situ experiment and its simulation
3.2. 原位实验及其模拟

The experimental uniaxial compressive test of a concrete sample described in our recently published paper [73] is selected as the benchmark for the FDEM simulation. The size of the concrete sample was about 1 mm × 1 mm × 1 mm and it contains 52 small aggregates. For more information about the experimental procedures and sample preparation, please refer to [73]. At the experimental concrete boundary where the concrete sample contacted with the steel loading platens, a voxel-mixed zone with fuzzy 2D images occurred. Since the XCT images became unclear near the loading boundaries due to imaging artefacts, we heuristically cut out an interest volume that is about 85% of the full height of the sample, which is shown in Fig. 7(a). Given that [73] used image processing to separate the different aggregate phases, in Fig. 7(b) the grey value intensities could be directly taken for binarization. However, there is still image noise that could be identified incorrectly as small ‘aggregates’, for example inclusions in the red ellipse of Fig. 7(b). These are removed after the labelling of various particulate phases in Fig. 7(c) by measuring and deleting all those small ‘aggregates’. The volume ratio of the largest deleted small ‘aggregate’ to the smallest real aggregate was <1%, which confirms that such deletions do not significantly alter the underlying composition. Finally, in Fig. 7(d) surficial meshes of different aggregates are generated with given numbers to retain a nearly linear relation between aggregate surface mesh number and surface area, using an open-source mesh tool, Iso2mesh [78]. Care is taken to avoid mesh intersection between touching aggregates, which are connected by ITZ CIEs, the same as the connection between aggregate and cement paste (CP) phases. The volume mesh of aggregates is produced using the open-source mesh tool Gmsh [79] with controlled mesh quality and number. Each simulation of the experiment and other virtual concrete samples involved about 350,000 solid and 700,000 cohesive elements. According to relevant studies of Wang et al. [24] having nearly identical number of aggregates to our model, 236,260 solid elements are found to be sufficient for the simulation of fracture behaviour of a meso-scale cubic concrete sample, from which the number of elements herein adopted are deemed adequate and no further mesh sensitivity studies are performed.
我们最近发表的论文[73]中描述的混凝土样本的实验单轴压缩试验被选为 FDEM 模拟的基准。混凝土样本的尺寸约为 1 mm × 1 mm × 1 mm，包含 52 个小骨料。有关实验程序和样本制备的更多信息，请参阅[73]。在混凝土样本与钢加载板接触的实验混凝土边界处，出现了模糊的二维图像的体素混合区。由于成像伪影，XCT 图像在加载边界附近变得不清晰，因此我们启发性地切割出一个约占样本全高 85%的感兴趣体积，如图 7(a)所示。鉴于[73]使用图像处理分离不同的骨料相，在图 7(b)中，灰度值强度可以直接用于二值化。然而，仍然存在图像噪声可能被错误识别为小“骨料”，例如图 7(b)中红色椭圆内的夹杂物。这些在图 7(c)中通过测量和删除所有这些小“骨料”后被标记的各种颗粒相移除。 最大删除的小“骨料”与最小真实骨料的体积比小于 1%，这确认了此类删除不会显著改变基础组成。最后，在图 7(d)中，使用开源网格工具 Iso2mesh [78]生成了不同骨料的表面网格，给定数量以保持骨料表面网格数量与表面积之间的近线性关系。注意避免接触骨料之间的网格交叉，这些骨料通过 ITZ CIEs 连接，类似于骨料与水泥浆(CP)相之间的连接。骨料的体积网格是使用开源网格工具 Gmsh [79]生成的，控制了网格质量和数量。每次实验模拟及其他虚拟混凝土样本涉及约 350,000 个固体元素和 700,000 个粘结元素。根据王等人的相关研究。 [24] 与我们的模型几乎相同数量的聚合物，发现 236,260 个固体单元足以模拟中尺度立方体混凝土样本的断裂行为，因此此处采用的单元数量被认为是足够的，且未进行进一步的网格敏感性研究。

Notably, with the aid of XCT at a resolution of 1.48 μm, some initial voids were found in the experimental concrete sample, though at a low volume ratio of about 1%. Hence, voids are not considered in this study since the emphasis is on pure influences of particle morphology features. Fig. 8(a) illustrates the force-displacement curves of experiments and simulations. In experiments, the 1 mm³ concrete sample, surrounded by aluminium cylinder to prevent fragments from being lost, was inserted into rotational and axial motion system (RAMS) loading grips at the Cornell High Energy Synchrotron Source (CHESS). The sample was compressed by lowering the stainless-steel platen in 2 μm increments until desired loads were reached. With the strain held constant, the concrete sample was rotated first 180° and then 360° for XCT and diffraction tests, respectively. The forces were obtained from the load cell; the sample stress was calculated by dividing the load cell reading during measurements by the initial sample area, 1mm².
值得注意的是，在 1.48 微米分辨率的 XCT 帮助下，实验混凝土样本中发现了一些初始空隙，尽管其体积分数约为 1%。因此，本研究不考虑空隙，因为重点在于颗粒形态特征的纯影响。图 8(a)展示了实验和模拟的力-位移曲线。在实验中，1 立方毫米的混凝土样本被铝圆筒包围，以防止碎片丢失，并被插入康奈尔高能同步辐射源（CHESS）的旋转和轴向运动系统（RAMS）加载夹具中。通过以 2 微米的增量降低不锈钢平板，样本被压缩，直到达到所需的载荷。在应变保持不变的情况下，混凝土样本首先旋转 180°，然后旋转 360°，分别进行 XCT 和衍射测试。力是通过负载传感器获得的；样本应力通过在测量过程中将负载传感器读数除以初始样本面积 1 平方毫米来计算。

1. Download: Download high-res image (915KB)
   下载：下载高分辨率图像（915KB）
2. Download: Download full-size image
   下载：下载完整尺寸图像

Fig. 8. Comparisons between simulations and experiments of the in-situ concrete fracture behaviour: (a) force-displacement curves; (b) evolutions of stress tensors with σ and σ′ denoting aggregate stress tensors of simulations and experiments; (c) and (e) raw CT images without any image processing of the concrete before and after fracture and (d) and (f) associated FEM results before and after fracture with the same normalized forces by their maximum compression forces. The slice of XCT image in Fig. 8(c) is from the bottom of the cubic volume of interest and is highlighted by bold rectangular edges in Fig. 7(a).
图 8. 原位混凝土断裂行为的模拟与实验比较：(a) 力-位移曲线；(b) 应力张量的演变，其中σ和σ′分别表示模拟和实验的骨料应力张量；(c) 和 (e) 断裂前后混凝土的原始 CT 图像，未经过任何图像处理；(d) 和 (f) 断裂前后与其最大压缩力相同归一化力的有限元分析结果。图 8(c)中的 XCT 图像切片来自感兴趣立方体体积的底部，并在图 7(a)中以粗体矩形边缘突出显示。

For our FDEM numerical simulations of the experiments, the loading platens are not explicitly modelled and the boundary conditions are applied to bottom and up layers of nodes of the concrete sample, which are fixed only on loading directions. We are therefore confident that our calculated sample strain in all loading steps matched those in steps 4 to 7 in experiments.
在我们对实验的 FDEM 数值模拟中，加载板没有被明确建模，边界条件应用于混凝土样本的底层和上层节点，仅在加载方向上固定。因此，我们有信心在所有加载步骤中计算的样本应变与实验中的第 4 到第 7 步相匹配。

It is plausible that relaxation occurs especially for the step with maximum displacement and the much longer scanning time due to the time required for experimental measurements, during which the upper loading platen is held in a fixed position to enable rotations of the sample aiming to conduct XCT and diffraction scan. The creep strains after this loading step are also calculated with the guidance of AS 3600:2018 [80] and then removed from Fig. 8(a). There are however still differences between the resulting curve and the numerical simulation. One reason for this phenomenon may be attributed to the roughness of sample surface planes in contact with the loading platens, on which newly observable fractures appear. From the raw XCT in Fig. 8(c) and (e) nearly every aggregate undergoes fracture and even breakage, indicating the necessity of their inclusion in the simulations of meso-scale concrete fracturing process. There are four different phases in them: voids or fractures, aggregates, cement paste and other high-density constituents. Black pixels with low greyscale values are deemed fractures or voids; white pixels with high greyscale values are for high-density constituents, such as calcium silicates; pixels with intermediate greyscale values denote cement paste or aggregates. Thresholds of grey values are carefully selected for segmenting phases, as described in [73]. A standard deviation filter is also used to threshold the aggregates from the cement paste, both of which have similar greyscale values (see [73] for details). Based on the fracture mechanics parameters for spherical aggregates in the meso-scale concrete simulation by [25], where no particle fracture is obtained, our simulations of in-situ experiments do exhibit aggregate fractures, as shown in Fig. 8(d) and (f). Although there is no further calibration of parameters for the simulation of in-situ experiments, the fracture patterns between numerical and experimental results are comparable, as well as the force-displacement curves. It seems that the simulations underestimate intra-aggregate fracturing, when comparing Fig. 8(e) and (f). Two main reasons may explain it: First, as we have described in experimental procedures, during every load step of the experiment, the sample is kept at constant stress and strain during rotation for XCT and diffraction measurements. Considering that these operations take at least half an hour, additional fractures due to the time-dependent behaviour may occur. Meanwhile, in our computational simulations, no strain rate influence on fractures have been considered. Second, aggregates in our simulations are homogenous and defect free, while in experiments the aggregates are heterogenous and include natural flaws along which fractures can initiate and evolute. It is anticipated that fracture behaviour between simulations and experiments would be closer if these two problems are solved by using more advanced computational resources.
放松现象的发生是合理的，特别是在最大位移的步骤中，以及由于实验测量所需的时间而导致的更长扫描时间。在此期间，上部加载板保持固定位置，以便样品进行旋转，以便进行 XCT 和衍射扫描。此加载步骤后的蠕变应变也根据 AS 3600:2018 [80] 的指导进行计算，并从图 8(a)中去除。然而，结果曲线与数值模拟之间仍然存在差异。这一现象的一个原因可能归因于与加载板接触的样品表面平面的粗糙度，在这些表面上出现了新可观察到的裂缝。从图 8(c)和(e)中的原始 XCT 来看，几乎每个骨料都经历了裂缝甚至破裂，这表明在中尺度混凝土破裂过程的模拟中包含它们的必要性。它们中有四个不同的相：孔隙或裂缝、骨料、水泥浆和其他高密度成分。 黑色像素具有低灰度值，被视为裂缝或空隙；白色像素具有高灰度值，代表高密度成分，如硅酸钙；中间灰度值的像素表示水泥浆或骨料。灰度值的阈值经过精心选择以进行相位分割，如[73]所述。还使用标准差滤波器将骨料与水泥浆进行阈值分离，两者具有相似的灰度值（详细信息见[73]）。基于[25]中对中尺度混凝土模拟中球形骨料的裂纹力学参数，未获得颗粒破裂，我们的原位实验模拟确实显示了骨料破裂，如图 8(d)和(f)所示。尽管在原位实验的模拟中没有进一步的参数校准，但数值结果与实验结果之间的破裂模式是可比的，力-位移曲线也是如此。比较图 8(e)和(f)时，模拟似乎低估了骨料内部的破裂。 两个主要原因可以解释这一点：首先，正如我们在实验程序中所描述的，在实验的每个加载步骤中，样品在进行 XCT 和衍射测量时保持恒定的应力和应变。考虑到这些操作至少需要半个小时，可能会由于时间依赖行为而发生额外的断裂。同时，在我们的计算模拟中，没有考虑应变速率对断裂的影响。其次，我们模拟中的聚合物是均匀且无缺陷的，而在实验中，聚合物是非均匀的，并包含自然缺陷，沿着这些缺陷可以引发和演化断裂。如果通过使用更先进的计算资源解决这两个问题，预计模拟与实验之间的断裂行为会更接近。

To compare stress of aggregates between simulations and experiments, volume-weighted stress tensors are defined as:(22)${\sigma }_{\mathit{ij}}=\frac{1}{\sum _{k=1}^{N_{\mathrm{agg}}}{V}_{\mathrm{agg}}^k}·\sum _{k=1}^{N_{\mathrm{agg}}}{\sigma }_{\mathit{ij}}^k{V}_{\mathrm{agg}}^k,$where N_agg is the number of aggregates; k meant k-th aggregate; V and σ_ij are volume and stress tensors. Considering the differences between the numerical and experimental force-displacement curves, σ_ij is also normalized by F_max/A, where F_max was the peak force and A was the approximate horizontal cross-section area of the concrete (1 mm²). Load steps denoted by the same number in the simulation and the experiment are of the same normalized force with respect to their associated F_max. The main stress component σ_zz, towards the compression direction, nearly coincide with each other during loading process as shown in Fig. 8(b), which indicates the accuracy of the proposed FDEM schematic. The other two diagonal stress tensor components, σ_xx and σ_yy, have a similar evolution tendency in the numerical and experimental results.
为了比较模拟和实验中集料的应力，定义了体积加权应力张量为： (22)${\sigma }_{\mathit{ij}}=\frac{1}{\sum _{k=1}^{N_{\mathrm{agg}}}{V}_{\mathrm{agg}}^k}·\sum _{k=1}^{N_{\mathrm{agg}}}{\sigma }_{\mathit{ij}}^k{V}_{\mathrm{agg}}^k,$ ，其中 Nagg 是集料的数量；k 表示第 k 个集料；V 和σij 分别是体积和应力张量。考虑到数值与实验力-位移曲线之间的差异，σij 也通过 Fmax/A 进行归一化，其中 Fmax 是峰值力，A 是混凝土的近似水平截面积（1 mm²）。在模拟和实验中用相同数字表示的加载步骤，其归一化力与相关的 Fmax 相同。主应力分量σzz 在压缩方向上，在加载过程中几乎重合，如图 8(b)所示，这表明所提出的 FDEM 示意图的准确性。其他两个对角应力张量分量σxx 和σyy 在数值和实验结果中具有类似的演变趋势。

4. The role of aggregate morphology
4. 聚合物形态的作用

In this part, focus is given to the statistical analyses and discussions of the obtained results correlating particle morphology with concrete fracture behaviour. Three factors are discussed: D_f, aspect ratio in terms of D₂, and aggregate fraction. In real concrete, fine aggregate fraction in mortar is typically higher than in this study [95]. However, obtaining such high packing density of irregular-shaped particles is difficult. With the purpose of isolating the influence of the particle morphology on concrete fracture behaviour, we tried to eliminate other factors associated with particles, such as the particle size distribution, particle spatial position in the sample and its longest axis orientation. Generating densely-packed mono-sized aggregate fractions while ensuring that particle position and fabric (the longest axis orientation) remain unchanged in concrete containing different fine aggregate shapes is challenging. For this reason, and because the existing experimental data we used for validation is from a sample with an aggregate fraction of 20%, we keep the fine aggregate fraction of 20% for most simulations and focus on the relative contributions from the aggregate morphology on material behaviour.
在这一部分，重点讨论了统计分析和所获得结果的讨论，关联颗粒形态与混凝土断裂行为。讨论了三个因素：D、以 D²表示的长宽比和骨料分数。在实际混凝土中，砂浆中的细骨料分数通常高于本研究中的比例[95]。然而，获得如此高的非规则颗粒堆积密度是困难的。为了隔离颗粒形态对混凝土断裂行为的影响，我们尝试消除与颗粒相关的其他因素，如颗粒大小分布、样本中颗粒的空间位置及其最长轴的方向。在确保颗粒位置和结构（最长轴方向）在含有不同细骨料形状的混凝土中保持不变的情况下，生成密集堆积的单一尺寸骨料分数是具有挑战性的。 因此，由于我们用于验证的现有实验数据来自一个集料比例为 20%的样本，我们在大多数模拟中保持细集料比例为 20%，并专注于集料形态对材料行为的相对贡献。