Multiscale modeling of composite materials with DECM approach : shape effect of inclusions

This paper addresses the study of the stress field in composites continua with the multiscale approach of the DECM (Discrete Element modeling with the Cell Method). The analysis focuses on composites consisting of a matrix with inclusions of various shapes, to investigate whether and how the shape of the inclusions changes the stress field. The purpose is to provide a numerical explanation for some of the main failure mechanisms of concrete, which is precisely a composite consisting of a cement-based matrix and aggregates of various shapes. Actually, while extensive experimental campaigns detailed the shape effect of concrete aggregates in the past, so far it has not been possible to model the stress field within the inclusions and on the interfaces accurately. The reason for this lies in the limits of the differential formulation, which is the basis of the most commonly used numerical methods. The Cell Method (CM), on the contrary, is an algebraic method that provides descriptions up to the micro-scale, independently of the presence of rheological discontinuities or concentrated sources. This makes the CM useful for describing the shape effect of the inclusions, on the micro-scale. When used together with a multiscale approach, it also models the macro-scale behavior of periodic composite continua, without losing accuracy on the micro-scale. The DECM uses discrete elements precisely to provide the CM with a multiscale approach. Keywords—Cell Method (CM), Discrete Element Method (DEM), multiscale modeling, periodic composite continua.

of view, since it allows the use of small stiffness matrices, leading to a significant reduction in computation time. The enforcement of both equilibrium and compatibility on the boundaries of the discrete elements then allows the individual workspaces to exchange information, restoring the continuity of the domain.
The main reason for using the CM on the micro-scale, leaving to the DECM the task of managing the discrete elements to capture the behavior on the macro-scale, is that the algebraic approach of the CM overcomes some of the typical disadvantages of the differential formulation. In particular, the CM can easily treat any type of singularity, including concentrated forces and discontinuities in the rheological properties. Therefore-unlike the Finite Element Method (FEM), the most used differential method in numerical analysis-the CM allows the modeling of phenomena up to the scale of single inclusions or interfaces. This makes the DECM particularly useful for modeling periodic composite continua. The purpose of this paper is to exploit precisely this capacity of the DECM, in order to study the stress field in cement-based composites under compression loads, with particular attention to the shape of the coarse aggregates.
Some experimental tests performed by various researchers in the past [11]- [19] already clarified how the strength of cement-based composites depends on the coarse aggregates, as well as on the water to cement ratio and the degree of compaction. For example, it is an established fact that cracks initiate at the interfaces between aggregates and mortar matrix [20]. Furthermore, the stress for which cracks develop [21]- [23] and the crack propagation speed along the interfaces [24] vary with the aggregate properties. In particular, the critical stress level at the interfaces depends on the difference between the elastic moduli of the matrix and the aggregates, since large differences induce higher tangential, radial and/or shear stresses at matrix aggregate interface. When the strength of the aggregates is comparable to the strength of the matrix, the cracks can also propagate through the aggregates [23]. Otherwise, the cracks make their way mainly through the matrix, or the matrix/aggregates interfaces. The shape of the coarse aggregate is decisive in inducing crack initiation. In fact, the smooth gravel leads to cracks at lower stresses than the rough and angular crushed rock [25]. Moreover, in the Multiscale modeling of composite materials with DECM approach: shape effect of inclusions E. Ferretti T absence of an interface treatment that makes the interface stronger, the concrete with spherical aggregates show a wide debonding of the aggregates, whereas the concrete with crushed aggregates show large crushing of the aggregates [26]. In both cases, the strength of the concrete is much lower than the strength of the matrix. Actually, the main function of the aggregates is to make concrete economic by reducing the requirements of cement, the most expensive ingredient, in addition to giving volume stability and increasing durability. The aforementioned experimental observations also deserve a numerical evaluation. The DECM approach now allows us to investigate them in more detail.

II. SIMILARITIES AND DIFFERENCES BETWEEN THE DECM
AND PREVIOUS DEM APPROACHES The basic idea of the DECM is to perform the CM analysis on a set of subdivisions of the modeling domain, forming stiffness matrices for each subdomain. Contrary to a Finite Element (FE) approach, the DECM does not assemble the local matrices into a global matrix, but enforces equilibrium and compatibility through an iterative analysis on the subdomains, by means of stabilization cycles (Section III.A). This reduces the computation time significantly and allows the modeling of large structures.
The distinct elements approach connotes the multiscale extension of the CM as a DE approach, in a broad sense. In fact, according to the original definition of discrete element method (DEM) proposed by Cundall and Hart, a DEM is any numerical technique that allows finite displacements and rotations of discrete bodies-including complete detachment-and automatically recognizes new contacts, as the simulation progresses. Well, the CM is actually able to manage finite displacements and rotations of deformable bodies, identify the points of enucleation of the cracks and model the propagation of the cracks, also in the cases of bifurcation of the cracks, automatically recognizing the positions in which the two edges of a crack come into direct contact [27]. Now the DECM inherits all these features from the CM.
To tell the truth, there are many differences between the DECM and the previous DEM approaches. The first difference is the generation of the geometry of the model: unlike the other DEM approaches, the DECM does not reproduce the shape of the modeling continuum by filling a geometric shape of particles. In fact, the discrete elements of the DECM are continua, not particles, and the set of the discrete elements covers the entire modeling area. This makes it possible not to use the filling procedure, which is one of the most delicate and time-consuming phases of the DEM approaches.
Even the subsequent calibrations and scaling processes of the DEM are time-consuming. Since they are a direct consequence of using a filling process, the DECM no longer needs them. In particular, unlike the DEM [28]- [31], the DECM does not need particular contact constitutive relations to establish interactions between the aggregates not in direct contact, when they are within a predetermined interaction area. In fact, the intrinsic nonlocal nature of the DECM variables-which include the lengths-scales--makes it possible to automatically take into account medium-and longrange interactions.
It is worth noting that the discrete elements of the DEM can be both rigid and deformable. In particular, they are rigid [32], [33] or deformable [34] particles in the modeling of geomaterials and particulate matter, such as powders or granules, and deformable blocks in the modeling of continua [29], [35]- [39]. The second case, however, involves some additional computational problems with respect to the first one, due to the poor performance of the internal mesh used for deformable blocks [40]. Therefore, the DEM is not useful for studying the stress field within deformable blocks. Also for this reason, some DEM approaches use rigid blocks even for the modeling of continua, in particular for the modeling of historical masonry structures [41], [42]. In this latter case, however, the DEM analysis could underestimate the strength capabilities of a real building, predicting collapse conditions not confirmed by the experimental results [43].
The DECM overcomes the drawbacks associated with the use of deformable and rigid blocks for the modeling of continua, as it has no problem generating efficient meshes for the sub-domains, even when they consist of more than one material [44]. Consequently, the DECM can provide descriptions of the stress field within single inclusions and on the interfaces between different materials [7].
Another difference between the two approaches concerns the solving equations. All the DEM approaches developed so far perform the modeling in the time domain, even for static or quasi-static problems. In fact, the static DEM solution is actually a solution to the steady state of a dynamic problem, since the discrete elements move according to Newton's second law of motion. This involves some disadvantages, from the numerical point of view, because the dynamic approach requires the calibration of the stable time step to allow the convergence of the numerical solution. Moreover, the computational cost of the dynamic relaxation technique also requires a preliminary evaluation of the minimum number of contact points to obtain the correct solution, when the discrete elements of the DEM are deformable.
Unlike the DEM, the DECM solves static problems in the space domain. Therefore, the DECM does not require any calibration of the stable time step and has no limits in the number of contact points. This makes it easier to find the static solution with the DECM than with the DEM.
Lastly, the time steps are often very small in dynamic analyses with the DEM, due to numerical stability requirements [45]. Moreover, determining the critical time step requires some approximations [46]. Consequently, the explicit time stepping DEM algorithms are quite effective just for quasi-static analyses. Performing dynamic analyses with the CM, on the contrary, does not require time-consuming calibrations of the critical time step [47].
As already mentioned in Section I, the DECM approach is particularly useful for modeling periodic continua and, specifically, periodic composite continua. In these cases, the DECM considers each basic unit of the periodic model as a distinct discrete element. Then, the DECM arranges the discrete elements into rows and columns, as if they were elements of an array, identifying the row and column indexes based on the coordinates of the boundary nodes. Therefore, the dimension of the array depends on the dimension of the problem: it is 2 for two-dimensional continua and 3 for threedimensional continua. For the sake of simplicity, the following Sections will deal with applications of the DECM to two-dimensional periodic composite continua, consisting of a matrix and inclusions of different shapes.

III. PERIODIC ROUND INCLUSIONS
Let the two-dimensional modeling domain be a 3 3  array of square discrete elements, with sides equal to 10 mm , and consider the two cases of presence ( The load condition for both specimens in Fig. 1

A. How to Enforce Equilibrium and Compatibility on the Discrete Elements
The boundary conditions on the sides of the discrete elements are unknown at the beginning of the computation, in terms of both forces and displacements, except for the sides of the outer boundary. Therefore, the DECM code must identify the boundary conditions for the inner sides of the domain, before proceeding with the computation. To this aim, the DECM code generates a series of twin nodes along the common sides of adjacent discrete elements. The twin nodes are pairs of nodes with the same coordinates, one on a discrete element and the other on the adjacent discrete element.
At the beginning of the computation, both the nodal displacements and the nodal forces on the twin nodes are set equal to zero and each discrete element moves independently of the presence of adjacent elements. The DECM code performs an iterative analysis on the twin nodes, in order to find the pairs of nodal displacements and nodal forces that restore continuity between the discrete elements.
The iterative analysis consists of two phases [7]: in the first phase, the DECM code processes the discrete elements with the same row index and, in the second phase, the DECM code processes the rows. Both phases, in turn, consist of two iterations, since the DECM code processes the discrete elements and the rows twice, in each phase. In particular, in the first phase (iteration on the column index), the DECM code processes the discrete elements, the first time, from left to right (with increasing values of the column index) and, the second time, from right to left (with decreasing values of the column index). Similarly, in the second phase (iteration on the rows), the DECM code processes the rows, the first time, from top to bottom (with increasing values of the row index) and, the second time, from bottom to top (with decreasing values of the row index).
The iterative analysis begins by executing the first phase on the first row. In particular, in the iteration with increasing value of the column index, the DECM code imposes the forces on the left sides of the c n elements and calculates the forces on the right sides (where the displacements are known). In the iteration with decreasing value of the column index, the code imposes the displacements on the right sides and calculates the displacements on the left sides (where the forces are known).
The DECM code applies the left-to-right-to-left procedure to the first row iteratively, updating both displacements and forces through a bisection technique. The procedure ends when the minimum relative displacement between the twin nodes of the first row is less than a predetermined value. This restores compatibility on the inner vertical sides of the first row (see the first row in Fig. 3, for 1 i k  ). Then, the DECM code starts processing the rows (second phase). In particular, in the iteration with increasing value of the row index, the code imposes the forces on the c n upper sides of the rows and calculates the forces on the c n lower ) sides (where the displacements are known). In the iteration with decreasing value of the row index, the code imposes the displacements on the c n lower sides and calculates the displacements on the c n upper sides (where the forces are known). Whenever the code updates the row index, it also performs a new left-to-right-to-left procedure on the c n elements of the new row (second and third rows in Fig. 3, for 1 i k  ). The top-to-bottom-to-top procedure updates the loads and the displacements of the r n rows iteratively, using a bisection technique. The procedure ends when the minimum relative displacement between all twin nodes is less than the same predetermined value of the first phase. This restores compatibility on the inner horizontal sides of the rows. Fig. 3 shows the deformations of the discrete elements after the first 12 top-to-bottom-to-top iterations, for the 3 3  array of Fig. 1b.

B. Stress Field Analysis within the Discrete Elements
Figs. [4][5][6] show the three stress components for the specimen without inclusions, while Figs. [7][8][9] show the stresses of the specimen with round inclusions. To allow a direct comparison, the upper and lower limits of the color scales are the same for each stress component, even for the following Sections. The stress values that exceed the upper limits are in dark red and the stress values below the lower limits are in dark blue.
From the comparison between Figs. 4 and 7, it follows that the normal stress y  tends to concentrate within the round inclusions, which causes y  to decrease in the matrix around the inclusions. As shown in [7], a similar behavior also occurs for shear loads. In particular, both for shear and axial loads, the normal stresses concentrate near the boundaries of the inclusions but y  reaches its maximum value (in absolute value) at the ends of the constraint. Moreover, the maximum normal stress y  (in absolute value) is lower in the case with inclusions than in the case without inclusions. This is a consequence of the combined action of the stiffness difference (between matrix and inclusions) and the Poisson effect. In fact, due to the greater macroscopic stiffness of the specimen with inclusions, the horizontal displacements caused by the Poisson effect are smaller for the specimen in Fig. 1b than for the specimen in Fig. 1a. Therefore, also the horizontal components of constraining reaction that cancel the horizontal displacements are smaller for the specimen with inclusions than for the specimen without inclusions. Due to the Poisson effect, the cancelation of smaller horizontal displacements causes less extra-strains in the longitudinal direction and, ultimately, a lower concentration of y  at the ends of the constraint.
Furthermore, the high stiffness of the inclusions leads the matrix enclosed between inclusions of the same column to compress more than the remaining material of the matrix. Consequently, the matrix between inclusions of the same column bears normal stresses y y p   (in absolute value), greater than the average axial stress of the matrix. The high stiffness of the inclusions also leads the matrix to the right and left of the inclusions to compress less than the remaining areas of the matrix. Therefore, the matrix enclosed between inclusions of the same row bears normal stresses y y p   (in absolute value), which are lower than the average axial stress of the matrix and change sign near the inclusions. This effect is similar to that observed in [48] for radiant heat floors. The reason for the complex stress field x  in Fig. 8 is the Poisson effect on the bi-material cross-sections. In fact, a single column of three elements with round inclusions and subjected to a uniformly distributed compressive axial load deforms as shown in Fig. 10 (for the uniformly distributed tensile axial load, see [44]). That is, the transverse strain x  is positive and reaches its maximum value where the axial strain y  is maximum (in absolute value), due to the lower stiffness on the cross-section. Therefore, the deformation of the crosssections is greater for the higher values of y  , which is maximum on the cross-sections that do not intersect the inclusions and are at a distance from the constraint at least equal to half the length of the constraint.
In the 3 3  array, the material continuity along the common vertical sides of adjacent elements of the same row prevents the cross-sections from deforming as in Fig. 10 (Fig.  11). Consequently, adjacent discrete elements interact with each other through mutual actions, which modify the stress field along the inner vertical sides of the discrete elements and give rise to the positive and negative normal stresses x  in Fig. 8.
In brittle materials, which generally have a low tensile strength, the positive values of x  on the vertical inner sides can trigger the propagation of vertical cracks between the discrete elements.
Furthermore, the high difference in stiffness between the matrix and the inclusions causes some negative x  to concentrate on the upper and lower portions of the interfaces between matrix and inclusions, in both the 3 1  (Fig. 10) and the 3 3  arrays (Fig. 8).  array (Fig. 11). It is worth noting that a downward curvature also characterizes the arrays without inclusions, but, in the latter case, the curvature is very weak (see Fig. 12 for the 3 3  array of Fig. 1a). Consequently, the normal stresses The difference in stiffness is also responsible for the shear stresses xy  that arise along the boundaries of the inclusions ( Fig. 9), significantly modifying the typical behavior of xy  for specimens without inclusions (Fig. 6). The presence of shear stresses along the boundaries of the inclusions is particularly harmful for the interfaces, as it can cause the detachment of the inclusions from the matrix. Furthermore, also the matrix in the upper part of the first row can damage due to the high values of xy  caused by the inclusions. This is a further effect of the mutual actions between adjacent elements of the same row. Lastly, the effect of the round inclusions on the maximum value of xy  at the ends of the constraint is similar to that already observed for y  , since the maximum shear stress on the constrained corners is lower for the specimen with inclusions ( Fig. 9) than for the specimen without inclusions (Fig. 6). In contrast, the maximum absolute value of x  is greater for the specimen with inclusions (Fig. 8) than for the specimen without inclusions (Fig. 5). However, the points where x  reaches its maximum value in the presence of inclusions are no longer the ends of the constraint-as instead is the case in Fig. 5-since x  concentrates on the points of the constraint that lie below the inclusions (Fig. 8). Therefore, the lower corners of the specimen benefit from the presence of the inclusions in any case.
In conclusion, in the case of uniaxial loading, the round inclusions decrease the concentration of the three components of stress at the ends of the constraint, delaying damaging effects on the lower corners of the specimen. On the other hand, however, the stress concentration along the boundaries of the inclusions can cause detachments between matrix and inclusions. This explains the experimental evidence that cracks initiate on the interfaces between mortar matrix and aggregates. Lastly, some vertical cracks can enucleate on the upper sides or between the inclusions of the same row.
It is worth noting that the stress concentrations at the boundaries of the inclusions, the upper sides, and the inner vertical sides cause cracks to start propagating at lower axial loads, compared to the case without inclusions. This confirms the experimental results on concrete with untreated matrixaggregate interfaces [26], according to which the strength of a composite consisting of matrix and inclusions is much lower than the strength of a specimen consisting only of the matrix.

IV. PERIODIC IRREGULAR INCLUSIONS
In order to verify how the crushed rocks modify the stress field with respect to rounded river gravels, the DECM code presented in this paper uses a new inclusions tool, which starts from a round inclusion to generate inclusions with random shapes. In particular, the inclusions tool generates a round inclusion that is actually a regular polygonal inclusion with a large number of sides (the number of sides in Fig. 2 is equal to  40). Then, the inclusions tool eliminates a random number of nodes located at random positions. This results in randomly shaped irregular polygons that can also include some rounded sides (Fig. 13).    The irregularity of the inclusions has an effect also on x  and xy  . In fact, as already observed for y  , Figs. 18 and 19 show that a greater degree of irregularity of the inclusions increases the maximum values of both x  and xy  , which tend to concentrate at the corners of the inclusions, within the inclusions. If compared to the stress fields for round inclusions (Figs. [7][8][9], this shifts the vulnerability from the interfaces to the inclusions and explains the two failure mechanisms observed in [26] for concrete: wide debonding of the aggregates for spherical aggregates and large crushing of the aggregates for crushed aggregates. Since the strength of the inclusions is greater than the strength of the interfaces, the concentration of the stresses within angular inclusions also explains why the smooth gravel leads to cracks at lower stresses than the rough and angular crushed rock [25]. The inclusions tool presented in Section IV can also diversify the random shapes of the polygonal inclusions into each discrete element. This allows a more realistic modeling of the stress fields in cementitious matrices with aggregates of natural origin. Moreover, diversifying the random shapes allows us to compare the stress fields inside and around inclusions of different shapes, also highlighting the local effect on the stress field deriving from particularly sharp inclusions.

V. IRREGULAR INCLUSIONS OF RANDOM SHAPE
Figs. [20][21][22] show the deformations and the stress fields for the same 3 3  array with three different sets of random inclusions (the distributed load and constraint conditions are the same for the three cases). As we already observed for periodic irregular inclusions, the three components of stress concentrate within the more pointed inclusions. The exact position of stress concentration depends on the component of stress and the orientation of the inclusion with respect to the load direction. In particular, y  concentrates on the sides that are almost parallel to the load, while x  concentrates on the sides that are almost orthogonal to the load. The stress xy  has a more complex behavior and can both concentrate on the boundary and spread into the inclusion. In this latter case, in fact, the direction of the sides is not as decisive as, rather, the relative positions of the pointed corners.
In conclusion, the three components of stress are highly variable in the arrays with irregular inclusions of random shape, as also results from the great variability of the maximum stress values for the three cases (Figs. [20][21][22]. Therefore, in addition to the shape effect that concerns the geometric characteristics of the specimen [49], the stress fields also suffer from a shape effect related to the geometric characteristics of the inclusions. A secondary effect of this is that the deformed configurations lose symmetry, despite the symmetry of the axial load. In particular, the asymmetry of the deformed configurations for the three sets of inclusions decreases from Fig. 20 (plotted for the set of inclusions with the highest degree of irregularity) to Fig. 22 (plotted for the set of inclusions with the lowest degree of irregularity).
The concentration of stresses within the inclusions makes the angular inclusions more vulnerable than the rounded inclusions, as observed also experimentally [26]. In the case of brittle inclusions, the most serious dangers for the angular inclusions derives from the tensile stresses x  and the shear stresses xy  , which can both split the inclusions. On the contrary, the occurrence of crushing of the inclusions due to high values of negative y  is more rare.  array with irregular inclusions of random shape does not deform symmetrically, while the 6 4  array with round inclusions has a symmetric deformation. However, since each column has twice the number of elements compared to the previous examples, the maximum difference in the average degree of irregularity calculated along the columns (that is, in the direction of the axial load) is lower for the 6 4  array than for the 3 3  array. Therefore, on average, the columns subjected to axial load deform more similarly in the 6 4  array than in the 3 3  array. In other In the 6 4  array with round inclusions the constraint modifies the stress field up to a distance almost equal to half the length of the constraint (Figs. [24][25][26]. Beyond this distance and excluding the first row, the stresses have a periodic behavior, symmetric for y  (Fig. 24) and x  (Fig. 25) and skew-symmetric for xy  (Fig. 26).
For inclusions of random shape, on the contrary, it is not possible to identify an extinction distance from the constraint, beyond which the three components of stress become periodic. In fact, the local perturbations due to the irregularity of the inclusions do not allow the stress fields to assume a recurrent behavior. For the same reason, it is not even possible that the stress fields of x  and y  are symmetric or that the stress field of xy  is skew-symmetric.
The multiscale approach with separate workspaces to store the data of each discrete element makes it possible to focus attention on individual discrete elements, just as it is possible to extract the entry in the i-th row and j-th column of a twodimensional array. By way of example, Figs. [27][28][29] show the discrete plots of the three stress components for the   It is worth noting that the stresses y  in Fig. 27 change sign crossing the upper-right side of the inclusion, a side that is almost parallel to the axial load, due to the high difference in stiffness between matrix and inclusion. A high gradient of y  also characterizes the lower-right side of the inclusion, but to a lesser extent as this second sub-vertical side forms a wider angle with the load direction. Lastly, the y  stress field within the inclusion finds its privileged direction in the most sub-vertical line among those that connect the corners of the inclusion boundary.  VI. CONCLUSION This paper shows the results of a DECM analysis for periodic composite continua consisting of a matrix with dispersed inclusions, such as concrete solids with aggregates of various shapes. The analysis allows the modeling of both micro-and macro-behaviors of periodic composite continua, thanks to the multiscale approach with discrete elements (DE) and the micro-scale analysis of the Cell Method (CM).
The simulations carried out for uniform uniaxial loads with various shapes of the inclusions confirm many experimental results-known from uniaxial tests on concrete-that so far did not find a numerical explanation. In particular, the paper highlights how the roundness and roughness of the inclusions affect the stress field:  Round inclusions modify the stress field significantly, in particular when the difference in stiffness between matrix and inclusions is high. The alterations in the stress field consist of stress concentrations on the matrix-inclusions interfaces, between the inclusions aligned along the direction orthogonal to the uniaxial load, and on the upper sides of the specimens. These stress concentrations are detrimental for composite continua, particularly if they consist of brittle materials, and can both cause detachments of the inclusions and trigger the propagation of cracks for lower load values than in specimens without inclusions.
Since the matrix-inclusions interfaces are generally weaker than the matrix and the inclusions, the cracks initiate on the interfaces. Furthermore, the crack initiation for relatively low load values means that the strength of composite continua is lower than the strength of the matrix.  Irregular inclusions modify the stress field more significantly than round inclusions, increasing stress concentrations on the matrix-inclusions interfaces, in particular near the corners of the inclusions. This decreases the load of crack initiation on the interfaces further. Therefore, as with round inclusions, the strength of the specimens with angular inclusions is less than the strength of the matrix.  Irregularities in the shape of the boundary amplify the shape effect of angular inclusions. In fact, for each assigned number of sides of the angular inclusion, the stresses concentrate more at the corners when the irregularity degree of the polygonal boundary is greater.  With round inclusions, the stresses tend to concentrate on the matrix-inclusions interfaces and in the matrix in the immediate vicinity, while, with angular inclusions, the stresses tend to concentrate within the inclusions. This means that the failure mechanism depends on the shape of the inclusions, leading to the debonding of the inclusions when they have a round shape and to the crushing or splitting of the inclusions when they have a polygonal shape. Since the matrix-inclusions interfaces are weaker than the inclusions, the first failure mechanism occurs for lower load values compared to the second failure mechanism. Therefore, the strength of specimens with round inclusions is less than the strength of specimens with angular inclusions. The DECM numerical analysis has also clarified some peculiar aspects of the positions of concentration of x  and y  along the boundaries of the inclusions, never highlighted by experimental tests. That is, the normal stress in the direction of the axial load concentrates on the sides that are almost parallel to the load and the normal stress in the direction orthogonal to the axial load concentrates on the sides that are almost orthogonal to the load.
The latter and the previous numerical results significantly improve the knowledge of the stress field in composite continua, compared to the FE (Finite Element) analyses. The DECM code also offers a more performing numerical analysis with discrete elements since, unlike other DEM (Discrete Element Method) codes, it performs the static analysis directly in the space domain. This allows us to avoid many typical drawbacks of the DEM analysis.  array with irregular random inclusions in Fig. 26