Modeling of Near-Surface Aggregate Size Distributions in Concrete

1. Introduction

Concrete is a composite material fundamentally consisting of cement, mixing water, fine aggregates (sand), and coarse aggregates (gravel). Upon hardening, the hydration products form a solid matrix surrounding the aggregates. On a mesoscopic scale (${10}^{-4}$ m < l < ${10}^{-1}$ m) [1], the heterogeneous mixture can be simplified as a two-phase system, consisting of a cementitious matrix and aggregates. The positions of the polydisperse aggregates are constrained by the size of the aggregates and the physical boundaries of the formwork. In addition, between the aggregates and the cementitious matrix, interfacial transition zones (ITZ) are formed [2,3,4,5].

The generation of realistic virtual concrete mesostructures is a significant area of research, driven by the need for accurate input in numerical simulations to predict concrete behavior [6,7]. These models aim to statistically replicate the configuration, size, and distribution of aggregates within the cementitious matrix [8]. The literature describes a variety of different methods of generating numerical mesoscale models. A comprehensive review was given by Ren [8]. A common approach is described as the take-and-place method. In this method, spherical or ellipsoidal aggregates, sampled from a known grading curve, are randomly placed inside the boundaries of a virtual specimen, ensuring no overlap with previously placed aggregates [9]. Furthermore, methods have been explored to generate more realistic aggregate geometries, moving from simple spheres or ellipsoids to polyhedral shapes, with the ability to consider features like irregular and concave faces [6,10,11]. Other approaches include techniques based on Voronoi tessellation and Delaunay triangulation to control particle shape and size in the system [12,13,14]. Some models also consider the ITZ between aggregates and mortar, which is known to be a critical area for crack initiation and concrete durability [15]. In addition, current research explores methods to consider non-uniform aggregate distributions due to segregation, which can significantly impact concrete properties [7].

Figure 1 (a) exemplarily shows a random packing of polydisperse spheres within a cube. The volume fraction $F_V$ describes the ratio between the total volume of the boundary geometry and the volume of all spheres. A vertical slice through the cube at any position, as depicted in Figure 1 (b), yields the area fraction $F_A$, describing the ratio between the total area of the sliced boundary geometry and the area of all intersected sphere segments. The function of area fraction $F_A\left(x\right)$ describes the progression of $F_A$ along the axis perpendicular to the vertical slices, see Figure 1 (c). It can be seen that $F_A\left(x\right)$ is never constant and exhibits values meandering around $F_V$.

In this study, a novel methodology is presented, where concrete mesostructures are generated by minimizing the difference between the actual and idealized homogeneous cumulative volume fraction functions $F_V\left(x\right)$ of a polydisperse distribution of spherical aggregates. The functions for $F_V\left(x\right)$ have been derived from integrating $F_A\left(x\right)$ over one axis of the geometry, see equation 1.

$$F_V\left(x\right)={\int }_0^x{F}_A\left(x\right)dx$$

(1)

The aggregates in the system are assembled based on the assumption of a constant $F_A\left(x\right)$ and thus linear $F_V\left(x\right)$ within the geometry, which results in smaller particles migrating towards the sample edges as a consequence of the optimization process. It is important to note that the method presented in this study does not determine the exact three-dimensional position and rotation of each individual aggregate in the system. This method rather assumes the depth of the aggregates according to the probability distribution functions of the different grain size groups in the system and their mutual interferences. The focus of this method is on the realistic representation of the border zones of the concrete laboratory specimens (wall effect) [16].

Existing wall-effect modeling approaches are typically based on explicit particle packing simulations or geometric placement algorithms, in which boundary effects arise from direct particle-to-particle or particle-to-wall interactions [11,17,18]. While these approaches provide detailed spatial representations, they require high computational effort. As reported in prior work, packing algorithms become increasingly inefficient at high volume fractions [6]. In contrast, the proposed distribution-optimized mesostructure estimation method formulates the problem as an inverse cumulative-volume matching task. Instead of explicitly simulating contact mechanics, the method enforces a target cumulative volume function along the specimen depth and determines depth-dependent aggregate distributions through optimization. This formulation reduces computational complexity while preserving the geometric characteristics of boundary-induced redistribution.

A Python implementation of the methodology is available online: https://github.com/ahaynack/DOME (accessed on 11 March 2026).

2. Mesostructure Generation

The first fundamental concept of the presented method is the representation of the volume of a spherical aggregate as the cumulative sum of its intersected circle segment areas along an axis, i.e. the integral of its area. Figure 2 (a) depicts a spherical aggregate, which is intersected by vertical slices. At the locations of the slices, a circle segment area with radius $r_{s_i}$ is formed. The function of circle segment areas $A_{s_i}\left(x\right)$ is shown in Figure 2 (b). Integrating $A_{s_i}\left(x\right)$ ultimately results in the function of sphere volume $V_{s_i}\left(x\right)$, as depicted in Figure 2 (c).

The shape of the volume function of each aggregate i is governed by the radius $r_i$ and the x-position of the center of the aggregate $x_{c_i}$. Figure 3 (a) shows four individual volume functions $V_i\left(x\right)$ of spheres with different radii and random positions along the x-axis. In Figure 3 (b), the total volume function $V_{tot}\left(x\right)$ is formed from the individual volume functions.

The second fundamental concept of the presented method is the comparison of the total volume function $V_{tot}\left(x\right)$ with an idealized homogeneous volume function $V_{lin}\left(x\right)$. This function assumes an even distribution of sphere volume along the x-axis. Hence, it represents the integral of constant area fractions across the x-axis, resulting in a linear volume function. Minimizing the differences between those two functions by rearranging the positions of the individual sphere centers $x_{c_i}$ results in an even distribution of the actual spheres in the system, constrained by the physical boundaries of the sample geometry. Figure 4 shows the total volume functions $V_{tot}\left(x\right)$ for three different spatial configurations of spheres within the system. Each function is compared to their linear volume function $V_{lin}\left(x\right)$ by calculating the coefficient of determination $R^2$, defined in equation 2.

$$R^2=1-{\displaystyle \frac{S{S}_{\mathrm{res}}}{S{S}_{\mathrm{tot}}}}$$

(2)

The term $S{S}_{\mathrm{res}}$ represents the residual sum of squares and is calculated as shown in equation 3:

$$S{S}_{\mathrm{res}}=\sum _{i=1}^n{(y_i-{\widehat{y}}_i)}^2$$

(3)

The term $S{S}_{\mathrm{tot}}$ represents the total sum of squares, proportional to the variance of the observed data, defined in equation 4:

$$S{S}_{\mathrm{tot}}=\sum _{i=1}^n{(y_i-\overline{y})}^2$$

(4)

where:

• $R^2$: The coefficient of determination.
• $y_i$: The y-value of the i-th point in the observed data.
• ${\widehat{y}}_i$: The corresponding y-value on the fitted linear line (the predicted value).
• $\overline{y}$: The mean of the observed y-values, calculated as ${\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{y}_i$.
• n: The total number of data points.

While the explicit rearrangement of individual sphere positions works well for a small number of spheres, simulations based on actual grading curves for concrete or mortar samples, potentially containing millions of aggregates, can become computationally prohibitive. To address this, the polydisperse grading curve is divided into distinct grain size groups. For each group, the spatial arrangement (depth) of the aggregates is modeled using a symmetrical Beta distribution, defined in equation 5. This approach relies on the assumption that the boundary conditions are identical at both ends of the sample. By setting the two standard shape parameters equal to each other ($\beta =\alpha$), the distribution becomes symmetric around $x=0.5$, controlled by a single parameter $\alpha$. Since the standard Beta distribution is defined on the interval $[0,1]$, the physical positions are scaled by the total sample length.

$$B(x,\alpha )={\displaystyle \frac{\Gamma \left(2\alpha \right)x^{\alpha -1}{(1-x)}^{\alpha -1}}{\Gamma {\left(\alpha \right)}^2}}$$

(5)

where:

• $B(x,\alpha )$: The probability density indicating the likelihood of an aggregate from a specific grain size group being located at normalized depth x.
• x: The normalized depth of the aggregate, defined as $x=d/L$, where d is the physical depth and L is the total length of the sample ($0\le d\le L$).
• $\alpha$: The shape parameter for the specific grain size group ($\alpha >0$). Due to the symmetry assumption ($\beta =\alpha$), this parameter controls the concentration of aggregates relative to the boundaries:

  –

  If $\alpha >1$, aggregates concentrate toward the center of the sample (bell-shaped distribution).

  –

  If $\alpha <1$, aggregates concentrate toward the border zones (U-shaped distribution).

  –

  If $\alpha =1$, aggregates are uniformly distributed along the length L.

• $\Gamma (·)$: The Gamma function.

Analogously to the individual aggregates, the distributions of aggregates result in the total volume function $V_{tot}\left(x\right)$, which can be compared to the linear volume function $V_{lin}\left(x\right)$. Figure 5 (a) exemplarily shows the distribution of n = 805 spheres of the same radius. The sphere positions are segmented into $n_{\mathrm{bins}}$ = 30 bins along the x-axis. These exemplary values were chosen to account for an adequate visual clarity. The magnitudes of sphere bins correspond to a symmetrical Beta distribution, which can be described solely with parameter $\alpha$. The $R^2$ score of the spheres’ total volume function $V_{tot}\left(x\right)$, generated from the arrangement of sphere bins, and the linear volume function $V_{lin}\left(x\right)$ is depicted in Figure 5 (b). Afterwards, both volume functions can be derived, resulting in the area functions along the x-axis, see Figure 5 (c).

Figure 6 displays a flowchart of the main optimization algorithm for the generation of the mesostructure. At first, an initial set of $\alpha$ parameters is set for the initial guess. These parameters are then passed on to the minimization module. The optimization problem is solved by minimizing the loss function, which consists of five sub-modules. For each aggregate group i, a histogram of spheres along the x-axis is generated based on the parameter of the symmetrical Beta distribution ${\alpha }_i$. Subsequently, the respective sphere volume function $V_i\left(x\right)$ is generated. Afterwards, the sum of each aggregate group’s sphere total volume function $V_{tot}\left(x\right)$ is calculated according to equation 6, which is being compared to the linear volume function $V_{lin}\left(x\right)$.

$$V_{tot}\left(x\right)=\sum _{i=1}^n{V}_i\left(x\right)$$

(6)

Where n denotes to total number of grain size groups. The resulting $R^2$ score is being minimized according to the loss function $L\left(\mathsf{\alpha }\right)$, as depicted in equation 7.

$$L\left(\mathsf{\alpha }\right)=1-R^2\left(\mathsf{\alpha }\right)$$

(7)

Where $R^2\left(\mathsf{\alpha }\right)$ denotes the coefficient of determination between the target $V_{lin}\left(x\right)$ and generated $V_{tot}\left(x\right)$ values. The formulation of $L\left(\mathsf{\alpha }\right)$ ensures that the loss decreases as the fit improves, i.e., as $R^2\left(\mathsf{\alpha }\right)\to 1$, the loss $L\left(\mathsf{\alpha }\right)\to 0$.

Calculating the aggregate histograms involves rounding continuous distributions to integer sphere counts. This creates a step-wise function where small changes in the shape parameters $\mathsf{\alpha }$ may not change the sphere count, causing standard gradient-based optimization methods to fail (local gradient is zero) [19]. To overcome this, the optimization is performed using Powell’s method [20]. This is a derivative-free optimization algorithm, which, instead of calculating a slope, minimizes the loss by searching along specific directions one by one (sequential line searches). The shape parameters $\mathsf{\alpha }$ are updated iteratively according to equation 8 [21].

$${\mathsf{\alpha }}_{\mathrm{new}}={\mathsf{\alpha }}_{\mathrm{old}}+\lambda \mathbf{d}$$

(8)

where:

• ${\mathsf{\alpha }}_{\mathrm{old}}/{\mathsf{\alpha }}_{\mathrm{new}}$: The vector containing the shape parameters for all grain size groups (i.e., $\mathsf{\alpha }=[{\alpha }_1,{\alpha }_2,\dots ,{\alpha }_n]$).
• $\mathbf{d}$: The search direction vector. In the initial steps, this vector isolates individual grain size groups (i.e., varying one parameter ${\alpha }_i$ while keeping the others fixed). As the algorithm progresses, $\mathbf{d}$ evolves into a combined direction, adjusting multiple shape parameters simultaneously.
• $\lambda$: The scalar step size that determines how far the algorithm moves the parameters along the direction $\mathbf{d}$ to reach the minimum loss for that specific search step.

By cycling through these line searches, the algorithm converges to the optimal aggregate distribution without requiring an analytical gradient.

Once the minimization is successful, the algorithm yields the optimized aggregate distribution according to the shape parameters $\mathsf{\alpha }$ of the final step of the optimization process. Ultimately, the aggregate distribution can be used to either calculate the depth-dependent mean bulk density or calculate the depth-dependent area fraction of the system. It has to be noted that the depth-dependent mean bulk density can only be determined if the individual density values of the cementitious matrix and aggregates are known.

The proposed mesostructure generation approach minimizes the deviation between the actual and idealized cumulative aggregate volume functions. In principle, the coefficient of determination ($R^2$) could be evaluated over the entire specimen depth, resulting in a global optimization objective. However, the primary structural effect investigated in this study is the near-surface accumulation of aggregates and thus, densification of the material caused by geometric boundary effects. When a purely global objective is used, deviations near the specimen boundaries may only have a minor influence on the overall loss function because the interior region dominates the cumulative volume profile.

To increase the sensitivity of the optimization to boundary effects, the loss function can be evaluated over a limited portion of the specimen depth. This region is defined by the so-called edge percentage $p_{\mathrm{edge}}$, which specifies the fraction of the specimen depth included in the objective function. For example, an edge percentage of 100 % corresponds to a global optimization, while a value of 20 % restricts the evaluation of the $R^2$ score to the first 20 % of the specimen depth. This formulation allows the optimization process to place stronger emphasis on accurately reproducing the near-boundary aggregate distribution.

In Appendix A, a sensitivity analysis of the optimization parameters was conducted. Based on the obtained results, a parameter combination of $n_{\mathrm{bins}}$ = 801 and $p_{\mathrm{edge}}$ = 20 % was selected for subsequent simulations, representing a compromise between numerical accuracy, near-surface agreement, physical plausibility, and computational efficiency.

3. Experimental Validation

3.1. General

For the experimental validation of the developed method, the depth-dependent density change of a concrete sample has been determined. With the aid of a high-resolution 3D laser scanner, the volume of a concrete sample was detected. The corresponding mass was measured with a precision scale. Afterwards, the experimental results are compared with the density values derived from the generated mesostructure.

3.2. Materials

The mix design investigated in this study consists of an Ordinary Portland Cement (CEM I 42.5 N), according to DIN EN 197-1 [22], with a w/c ratio of 0.55. The maximum aggregate size of the mixture was 8 mm (quarzitic, naturally rounded gravel). The grain size distribution can be taken from Table 1.

Table 2 shows the concrete mixture proportion used in this study. A volume fraction of aggregates $F_V$ of 0.67 was chosen. The density of aggregates ${\rho }_{aggregates}$ was determined as 2622 kg/m³.

The experimental determination of the depth-dependent density change of concrete was performed exemplarily on one cube with a length of 15 cm. After casting, the sample has been kept in the formwork for 1 day. Afterwards, the sample was placed in underwater storage for 6 days, followed by a storage at 20 °C and 65 % relative humidity for 21 days. The experiments were conducted at a sample age of 28 days. Prior to testing, the total sample density ${\rho }_{\mathrm{total}}$ was determined as 2252 kg/m³.

3.3. Methods

3.3.1. Experimental Determination of Depth-Dependent Density Change of Concrete

The depth-dependent density change of the sample has been determined with the aid of a high-resolution 3D laser scanner, analogously to the methodology presented in [23]. One vertical side of the concrete cube has been incrementally abraded with a grinding machine. After each grinding step, the mass and volume of the residual sample were determined. The mass was measured with a scale, the volume was measured by the 3D laser scanner. The experimental setup is displayed in Figure 7.

The initial mass and volume serve as the reference values for the determination of the depth-dependent density change. For each measurement point, the change in mass and volume in respect to the reference has been determined. Afterwards, the respective change in mean bulk density was calculated by division of change in mass by change in volume, see Figure 8. The depth of each grinding step was determined by calculating the mean height loss of the abraded sample side.

3.3.2. Generation of Mesostructure

In preparation for the generation of the mesostructure, the aggregates in the system were generated according to their grain size distribution related to $F_V$, see Figure 9.

The stepwise generation of the aggregates results in the following numbers of spherical aggregates for each mesh size, see Table 3. It is important to note that no aggregates smaller than 0.125 mm or greater than 8 mm have been considered in the generation of the aggregates. These fractions were attributed as part of the volume of the cementitious matrix. As a result, the target volume fraction $F_V$ of 0.67 was reduced to 0.657. In addition, the mesh sizes have been chosen analogously to the aggregates’ grain size distribution. Choosing a smaller mesh distance can result in the generation of a more continuous distribution of aggregates.

Figure 10 shows the $R^2$ scores of both the initial (non-optimized) and final (optimized) aggregate distributions compared to the linear target distribution. It can be seen that the total $R_{\mathrm{total}}^2$ score for both distributions yields values close to 1.0 (0.99956 and 0.99994). In contrast, with the focus on optimizing the distibution of aggregates near the edge of the sample, the coefficient of determination was also determined for the first 10 % of the sample length. This $R_{\mathrm{edge}}^2$ of the initial distribution results in a value of 0.89446, which is a difference of approximately 11 % compared to the respective $R_{\mathrm{total}}^2$. The final distribution yields an $R_{\mathrm{edge}}^2$ of 0.99720. Ultimately, the percentage difference of $R_{\mathrm{total}}^2$ and $R_{\mathrm{edge}}^2$ of the final distribution is approximately 0.27 %.

Figure 11 shows a detailed view of the edge regions previously presented in the zoomed-in sections of Figure 10. Although the earlier magnifications yielded high $R^2$ scores (close to 1), suggesting an almost linear cumulative distribution, local deviations near the boundaries were not clearly visible at that scale. To better illustrate these effects, Figure 11 directly compares the initial and final distributions in the edge region together with the linear target distribution. This representation reveals the nonlinearity at the sample boundaries more clearly. A perfectly linear cumulative volume distribution cannot be achieved due to the boundary conditions of the mesostructure generation algorithm. The cumulative distributions represent the integrated aggregate volume along the x-axis. Due to the aggregates’ spherical shape, their contribution to the total cumulative distribution follows a characteristic sigmoidal profile, as displayed in Figure 2 (c), rather than a linear increase. This geometric constraint inherently causes deviations from the linear target distribution. The deviation is strongest in the initial configuration. After optimization, where smaller aggregates move towards the boundary regions, the edge deviation is reduced, leading to a closer approximation to linear behavior in the final distribution.

Figure 12 depicts the depth-dependent mean bulk density values of the generated mesostructures. The green line shows the values derived from the initial $\alpha$ parameters. The purple line displays the density change derived from the optimized $\alpha$ parameters. It can be observed that both mesostructures form a pronounced depth-dependent density profile in the near-surface areas. However, due to the increased $R_{\mathrm{edge}}^2$ after the optimization process, the optimized mesostructure exhibits higher density values in the border zone compared to the non-optimized mesostructure.

The displayed values have been determined by equation 9.

$$\rho \left(x\right)={\displaystyle \frac{{\rho }_{\mathrm{aggregates}}{V}_{\mathrm{aggregates}}\left(x\right)+{\rho }_{\mathrm{cement}}{V}_{\mathrm{cement}}\left(x\right)}{V_{\mathrm{formwork}}\left(x\right)}}$$

(9)

Where $V_{\mathrm{aggregates}}\left(x\right)$, $V_{\mathrm{cement}}\left(x\right)$, and $V_{\mathrm{formwork}}\left(x\right)$ denote the cumulative volume functions of the aggregates, cementitious matrix, and formwork of the system. Function $V_{\mathrm{aggregates}}\left(x\right)$ was determined from the optimization algorithm. Function $V_{\mathrm{cement}}\left(x\right)$ was obtained by equation 10.

$$V_{\mathrm{cement}}\left(x\right)=V_{\mathrm{formwork}}\left(x\right)-V_{\mathrm{aggregates}}\left(x\right)$$

(10)

With known density of aggregates ${\rho }_{\mathrm{aggregates}}$ (2622 kg/m³), total sample density ${\rho }_{\mathrm{total}}$ (2252 kg/m³), and volume fraction of aggregates $F_V$ (0.67), the density of the cementitious matrix ${\rho }_{\mathrm{cement}}$ was derived as 1502 kg/m³.

Figure 13 shows the depth-dependent area fractions derived from the generated mesostructures with the initial (non-optimized) and final (optimized) $\alpha$ parameters. Analogously to the mean bulk density change, the area fraction of the final distribution shows peaks close to the border zones of the sample, indicating a higher accumulation of (finer) aggregates compared to the initial uniform distribution.

3.4. Comparison of Mean Bulk Densities

Figure 14 shows the comparison of the depth-dependent mean bulk density changes of the measurement data and the generated mesostructures (initial and optimized $\alpha$ parameters). It can be seen that the initial distribution contains lower density values towards the edge of the specimen compared to the measurement data ($R^2=0.71$). After the optimization process, which causes a migration of smaller particles towards the edges, the density values of the generated mesostructure increases and overlaps with the measurement data ($R^2=0.95$). With increasing depth, the mean bulk density values of both the optimized mesostructure and measurement data converge towards the total density of the sample (2252 kg/m³).

4. Comparison of Volume Fractions

The grading curve of Section 3.2 has been used to generate aggregates for different volume fractions $F_V$. Afterwards, the resulting aggregates were used to generate different mesostructures and ultimately derive the depth-dependent area fractions $F_A$. Figure 15 shows the resulting area fractions $F_A$ for volume fractions $F_V$ ranging from 0.10 to 0.70. These values have been chosen purely to visualize the effect of increased aggregate volume fractions on the resulting area fractions. The volume fraction of $F_V$ = 0.67 of the concrete sample used in the experimental validation is visualized as a red line. It can be observed that, with increasing $F_V$, the migration of smaller particles also rises and causes an increasing $F_A$ in the edges of the specimen. It is important to note that the peak values of $F_V$ and thus $F_A$ are limited by the maximum packing density of the specific grading curve. It is subject of future research to implement the functionality of limiting the generation local $F_V$ values based on the maximum packing density. For comparison, random close packing fractions of polydisperse circles in a two-dimensional space have been reported between 0.84 and 0.94, depending on the investigated grading curves [24,25]. These values provide a physical plausibility constraint for the generated mesostructures.

5. Discussion

5.1. Geometric Origin of the Near-Surface Densification

The proposed distribution-optimized mesostructure estimation method enforces a linear cumulative aggregate volume function along the specimen depth, which corresponds to a constant area fraction. Due to the spherical shape of the aggregates, their individual cumulative volume contributions follow sigmoidal curves. Near the specimen boundaries, only a portion of the spherical volume contributes to the cumulative volume within the considered depth interval. As a result, the effective volume contribution per unit depth is reduced close to the boundaries. To compensate for this geometric effect, the optimization algorithm shifts smaller aggregates toward the boundary regions. This leads to an increased local area fraction and therefore to higher mean bulk density values in the near-surface zone. Importantly, this behavior results from geometric constraints and is not imposed by predefined segregation rules. The method therefore captures the geometric component of the wall effect, i.e. the modification of aggregate distributions caused by boundary-induced geometric constraints on particle contributions [11,17,18]. It does not explicitly model rheological processes during casting, gravitational segregation, bleeding, or vibration effects. Instead, it isolates the structural consequences of geometric confinement under the assumption of statistical homogeneity in the bulk material.

5.2. Interpretation of the Experimental Agreement

The optimized mesostructure shows strong agreement with the experimentally determined depth-dependent mean bulk density ($R^2=0.95$). In contrast, the non-optimized configuration underestimates the near-surface density ($R^2=0.71$), indicating that a uniform aggregate depth distribution cannot reproduce the measured profile. The improved agreement suggests that geometric boundary effects contribute to the observed density increase in the near-surface region of the investigated specimen [26,27]. However, the experimental validation was performed for only one concrete mixture and one specimen geometry. Although the numerical investigations for different volume fractions show a consistent trend of increasing near-surface fine aggregate accumulation with increasing aggregate content, further experimental studies are required to confirm the general applicability of the method for different grading curves, maximum aggregate sizes, and casting conditions.

5.3. Symmetry Assumption of Boundary Conditions

The use of symmetric Beta distributions ($\beta =\alpha$) assumes identical boundary conditions at both specimen faces [11,17,18]. This assumption is appropriate for cubic specimens with similar formwork interfaces on both sides. In practical applications, however, casting direction, bleeding, sedimentation, or surface finishing may introduce asymmetric boundary effects [28,29]. The presented framework can be extended by allowing $\beta \ne \alpha$, which would enable asymmetric depth distributions when required. Such an extension would allow the method to account for directional casting effects or asymmetrical phenomena.

5.4. Limitation Regarding Maximum Packing Density

The current implementation does not explicitly include a local maximum packing constraint. Although the global aggregate volume fraction remains within realistic limits, local area fractions near the boundaries may approach values that should be verified against established packing density models [30,31]. Future work should therefore include local packing constraints to ensure the physical legitimacy of generated mesostructures, particularly for high aggregate volume fractions. Including such constraints would further improve the physical consistency of the model.

6. Conclusions

• A novel methodology has been presented, which generates concrete mesostructures by minimizing the $R^2$ score between the actual and idealized homogeneous cumulative volume fraction functions $F_V\left(x\right)$ of a polydisperse distribution of spherical aggregates.
• Ensuring the assumption of a constant $F_A\left(x\right)$ and therefore linear $F_V\left(x\right)$ within the geometry, the aggregates in the system are assembled in respect to the physical boundaries of the formwork. This results in smaller particles migrating towards the sample edges as a consequence of the optimization process.
• The presented method isolates and models the geometric contribution of boundary-induced aggregate redistribution. While it does not explicitly simulate casting-induced segregation mechanisms, it provides a framework for incorporating wall effects into mesoscale concrete simulations.
• For systems with high amounts of individual aggregates, the presented method assumes groups of equal aggregate sizes (based on grading curve). The distributions of the aggregate groups are described by a symmetrical Beta function. As a result, each aggregate group can be described with one parameter $\alpha$.
• The current iteration of the method only considers spherical aggregates. For a more realistic representation of the aggregates, a different method for the generation of the aggregates has to be chosen. The cumulative-volume-based optimization principle itself is not limited to spherical particles. However, during the generation of the mesostructure, the segmentation of grain size groups and representation of the aggregates with Beta distributions is not possible for unique, non-spherical aggregates. In this case, the individual positions of the aggregates need to be minimized. In addition, a rotation angle can be added to the optimization process to account for different aggregate rotations.
• Generated mesostructures were compared with experimental data. For this, the mean bulk density of a concrete sample was determined both experimentally and virtually. The mesostructure generated with the optimization method yielded good agreement with the measurement data ($R^2=0.95$). This can be attributed to the migration of smaller paricles towards the edges of the sample, which increased the density of the border zone. The density of a uniform, non-minimized aggregate distribution resulted in lower values compared to the measurement data ($R^2=0.71$).
• The method assumes the densities of both aggregates and cement matrix as constant values. In future iterations, the density of the cement matrix should also be viewed as a depth-dependent function. Analogously, different density values for the individual aggregate groups should be considered.
• The generation of mesostructures with different volume fractions $F_V$ shows an increasing accumulation of smaller aggregates towards the edges of the sample.
• Future research will focus on implementing local packing constraints, asymmetric boundary conditions, and extensions toward non-spherical aggregate geometries.