# Establishment and Application of Stochastic Mesoscopic Concrete Model

Authors: Shufang Wang; Xiaoliang Zhai; Yangyang Gao; Jialu Li; Wei Hu
Journal: Advances in Civil Engineering
 (2022)

---

## Abstract

Concrete is a kind of complex multiphase composite material, and the spatial distribution of each component of concrete is random, which has an impact on its physical and mechanical properties. This paper presents algorithms for generating two-dimensional (2D) and three-dimensional (3D) arbitrarily graded stochastic concrete aggregate models, graphical format conversion methods, and finite element applications. Aggregate shapes include circular, ellipsoidal, convex polygon, ellipsoidal, and spherical. These models take into account the randomness of appearance, position distribution, and particle size distribution, including the interface transition zone, and the aggregate content can reach a high level. An efficient interference algorithm of ellipsoidal aggregate is proposed, which can reasonably control the thickness of the mortar layer around the aggregate. The corresponding program is compiled by MATLAB to realize the generation of random mesoscopic concrete aggregate geometric model. A graphic format conversion method is provided, which greatly extends the applicability and generality of the random aggregation model. The usefulness of the method in this paper was verified by the mechanical analysis of random aggregate model specimens with different volume fractions and the chloride transport properties of different aggregate shape specimens. In addition, the results of the two cases show that the random aggregate geometry and aggregate content have a great impact on the physical and mechanical properties of concrete, and the randomness of the spatial distribution of aggregates makes the chloride ion transport path and the elastic modulus and Poisson’s ratio of concrete be random. The results of this paper can provide an important research foundation for the study of mesoscopic physical and mechanical properties, gradation theory, and durability of concrete.


---

## Body

## 1. Introduction

Cement concrete is made of coarse aggregate, cement, sand, water, and admixtures in a certain proportion through mixing, hydration, and solidification. The complexity of the original composition and molding process of concrete leads to its complex physical and mechanical properties. With the in-depth study of concrete materials and the complexity of concrete ratio design, some properties of concrete are difficult to explain from a macroscopic research and need to be explained from mesoscopic or even microscopic perspective. Therefore, researchers have studied the composition of concrete from a mesoscopic perspective and established a mesomechanical model to study the mechanical and physical properties of concrete [1–3]. Among them, concrete aggregates are the most important components of the mesomechanical model, whose shape, spatial distribution, volume rate, and gradation determine the applicability and reliability of the meso-stochastic research.In the study of concrete stochastic models, concrete can be assumed to be a two-phase composite consisting of aggregate and cement mortar [4–7] or a three-phase composite consisting of aggregate, cement mortar, and the interface transition zone (ITZ) [2, 3, 8, 9]. Pebble aggregates are generally assumed to be round or oval, and broken aggregates are usually assumed to be convex polygon or polyhedron [2–5, 10]. Among them, the generation of polygonal aggregate is more complex than that of circular aggregate. There are mainly four methods to generate polygonal aggregate: one is to generate inscribed polygon in the circle or ellipse [4, 7]; the second is to generate convex polygon within the circle and then conduct convex extension [10]; the third is to randomly generate triangular or quadrilateral aggregate base, then randomly grow into polygon [5]; the fourth is to use polar coordinates controlled by radius and direction angle [11]. Wittmann et al. [1] used polar coordinates, expressed the radius as a function of angle to generate the two-dimensional circular and polygonal aggregates based on the average radius. However, this method did not adequately consider the aggregate gradation. Wang et al. [2] proposed a procedure for generating random aggregate structures for rounded and angular aggregates based on the Monte Carlo random sampling principle and developed a method for generating meshes using the forward front approach. However, the polygonal aggregates generated by this method may have concave surfaces and sharp corners, reducing the reliability of the calculation. Leite et al. [6] introduced an algorithm to generate ellipsoids of a given size and shape and implemented the consideration of fiber reinforcement. Pan et al. [10] took the circle as the benchmark, obtained the elliptical aggregate by stretching, and established the convex polygon aggregate by generating the inscribed polygon of the circle. But the generated aggregates were mostly narrow and long with many sharp angles, which did not match the actual situation. Ma et al. [8] demonstrated a convex extension method to construct convex polygons and polyhedrons represented as crushed coarse aggregates in two- or three-dimensional spaces and proposed a occupation and removal method (ORM) to improve the efficiency of generating mesostructure models with a high aggregate content. Sayyafzadeh et al. [12] put forward the ‘‘equivalent space method,” which appears to be more convenient for both low and high volume fraction specimens and leads to more realistic concrete models with less random numbers. Qin et al. [9] used the spatial allocation method to generate on fully-graded polyhedral random aggregate, which gets the randomness of the appearance and distribution of aggregates involved, and the cohesive elements are globally inserted to capture the fracture propagation. But the method does not consider the control of the sharp angle of polygon, and it is not convenient to control the volume ratio of aggregate.The establishment of stochastic aggregate model is a complex mathematical problem, which often requires software programming. But the mathematical analysis software generally does not have the function of finite element calculation, which makes the stochastic aggregate model not convenient to use. Therefore, it is necessary to use appropriate methods to convert the aggregate graphic format constructed by the programming software into the file format that can be read by the finite element software. Hu et al. [11] established two-dimensional pebble and gravel mesoscopic concrete models based on MATLAB and converted format of the graphics file to improve the generality of the stochastic aggregation model. Liang et al. [13] established an aggregate structure that can be used by a variety of finite element software.This paper establishes a graded mesoscopic stochastic concrete model with arbitrary gradation and high volume fraction, which is capable of generating gravel or pebble aggregate in accordance with the actual situation; the aggregate model can be transferred to general finite element software. Firstly, the establishment process and generation algorithm of random aggregate models are presented, including circular aggregate, circular and elliptical random mixed aggregate, convex polygonal aggregate, and spherical-ellipsoidal random mixed aggregate mesoscopic models. Then, a method is proposed to convert random aggregate files in graphical format to a common format so that the generated mesoscopic aggregation models can be easily used by finite element software. At the same time, the mesh generation method after importing stochastic aggregates into the finite element software is introduced. Finally, the main role and significance of the stochastic concrete model are illustrated by two examples, including mechanical analysis of stochastic aggregate model specimens with different volume fractions and chloride ion transport properties of different aggregate shape specimens.


## 2. Principle of Random Aggregate Generation for Concrete

### 2.1. Parametric Space Representation of Aggregate Particles

In a certain space (Ω), the position of any aggregate particle in Ω is randomly distributed. In Cartesian coordinate system, the coordinates of reference pointO(XO, YO, ZO) of any aggregate particle can be expressed as(1)XO=f11x,y,z+η1f12x,y,z−f11x,y,zYO=f21x,y,z+η2f22x,y,z−f21x,y,zZO=f31x,y,z+η3f32x,y,z−f31x,y,z,where f11(x, y, z), f12(x, y, z), f21(x, y, z), f22(x, y, z), f31(x, y, z), and f32(x, y, z) are boundary curve functions of domain Ω, respectively; η1, η2, η3 are independent random numbers between 0 and 1. Specifically for prism specimens, the coordinates of any reference point O(XO, YO, ZO) can be expressed as(2)XO=Xmin+η1Xmax−XminYO=Ymin+η2Ymax−YminZO=Zmin+η3Zmax−Zmin,where Xmax, Xmin, Ymax, Ymin, Zmax, and Zmin are the maximum and minimum values of 3D Cartesian coordinate system on the boundary of domain Ω, respectively.Generally, pebble aggregate particles can be represented by sphere or ellipsoid and can be represented by circle or ellipse in 2D. Gravel aggregate can be represented by convex polyhedron or convex polygon. The circular aggregate particles can be described by three characteristic parameters: center coordinates (reference point)O(XO, YO) and radius r. The elliptical aggregate particles can be determined by the reference point O(XO, YO), the long half axis a, the short half axis b (or the ratio of long axis to short axis e), and the direction angle θ of the long axis. θ represents the angle between the long axis and the horizontal axis x. The convex polygon aggregate can be described by the number of sides n, the polar angle α, and the polar radius rp. Firstly, n is generated randomly:(3)ni=nmin+roundη1inmax−nmin,where η1i is a random number with uniform distribution between 0 and 1, namely η1i∼U0,1; nmax and nmin are the minimum and maximum values of polygon sides, generally within 4∼12; and round () is the rounding function.The polar angle (α) is also a random number. Assuming that α also obeys the uniform distribution, α can be generated by the following:(4)αi=η2i2π−0,where η2i∼U0,1. After obtaining n polar angles α, the sequence of n ordered polar angles α is obtained, which is actually the direction of the vertices of the polygon. The next step is to determine the specific positions of these vertices in these directions. Therefore, a polar radius rp is generated randomly as(5)rp,i=rp,min+η3irp,max−rp,min,where η3i∼U0,1; rp,min and rp,max are the minimum and maximum of the polar radius of the polygon, respectively. It could be noted that η1i, η2i, and η3i are three independent random numbers.The expression method of spherical aggregate is similar to that of circular aggregate, which can be described by the coordinate componentsXO, YO, ZO of reference point and radius r. Ellipsoidal aggregate particles can be determined by nine characteristic parameters: reference point (O(XO, YO, ZO)), semimajor axis of ellipsoid (a, b, c), precession angle (γ), nutation angle (φ), and spin angle (ψ).


### 2.2. Spatial Random Distribution of Aggregate Particles

Reasonable aggregate size distribution can effectively improve the compactness and durability of concrete. Fuller grading curve is recognized by scholars as one of the most ideal grading curve, because the concrete formulated by Fuller grading has theoretically the maximum density and strength. Fuller curve can be illustrated by the following [14]:(6)pd=100ddmaxn,where p(d) is cumulative volume fraction of aggregate within the particle size d, dmax is the maximum aggregate size, and n is the index of the equation, n = 0.45∼0.70. For the graded section [ds, ds+1], the aggregate volume contained is(7)Vagds,ds+1=pds+1−pdspdmax−pdminRagVcon,where Rag is the volume ratio of aggregate (%); Vcon is the volume of concrete; dmax is the maximum size of aggregate; and dmin is the minimum size of aggregate. When the concrete area, the aggregate volume ratio, and the aggregate grading section [ds, ds+1] are determined, the aggregate volume of the graded section can be calculated by (7).


## 2.1. Parametric Space Representation of Aggregate Particles

In a certain space (Ω), the position of any aggregate particle in Ω is randomly distributed. In Cartesian coordinate system, the coordinates of reference pointO(XO, YO, ZO) of any aggregate particle can be expressed as(1)XO=f11x,y,z+η1f12x,y,z−f11x,y,zYO=f21x,y,z+η2f22x,y,z−f21x,y,zZO=f31x,y,z+η3f32x,y,z−f31x,y,z,where f11(x, y, z), f12(x, y, z), f21(x, y, z), f22(x, y, z), f31(x, y, z), and f32(x, y, z) are boundary curve functions of domain Ω, respectively; η1, η2, η3 are independent random numbers between 0 and 1. Specifically for prism specimens, the coordinates of any reference point O(XO, YO, ZO) can be expressed as(2)XO=Xmin+η1Xmax−XminYO=Ymin+η2Ymax−YminZO=Zmin+η3Zmax−Zmin,where Xmax, Xmin, Ymax, Ymin, Zmax, and Zmin are the maximum and minimum values of 3D Cartesian coordinate system on the boundary of domain Ω, respectively.Generally, pebble aggregate particles can be represented by sphere or ellipsoid and can be represented by circle or ellipse in 2D. Gravel aggregate can be represented by convex polyhedron or convex polygon. The circular aggregate particles can be described by three characteristic parameters: center coordinates (reference point)O(XO, YO) and radius r. The elliptical aggregate particles can be determined by the reference point O(XO, YO), the long half axis a, the short half axis b (or the ratio of long axis to short axis e), and the direction angle θ of the long axis. θ represents the angle between the long axis and the horizontal axis x. The convex polygon aggregate can be described by the number of sides n, the polar angle α, and the polar radius rp. Firstly, n is generated randomly:(3)ni=nmin+roundη1inmax−nmin,where η1i is a random number with uniform distribution between 0 and 1, namely η1i∼U0,1; nmax and nmin are the minimum and maximum values of polygon sides, generally within 4∼12; and round () is the rounding function.The polar angle (α) is also a random number. Assuming that α also obeys the uniform distribution, α can be generated by the following:(4)αi=η2i2π−0,where η2i∼U0,1. After obtaining n polar angles α, the sequence of n ordered polar angles α is obtained, which is actually the direction of the vertices of the polygon. The next step is to determine the specific positions of these vertices in these directions. Therefore, a polar radius rp is generated randomly as(5)rp,i=rp,min+η3irp,max−rp,min,where η3i∼U0,1; rp,min and rp,max are the minimum and maximum of the polar radius of the polygon, respectively. It could be noted that η1i, η2i, and η3i are three independent random numbers.The expression method of spherical aggregate is similar to that of circular aggregate, which can be described by the coordinate componentsXO, YO, ZO of reference point and radius r. Ellipsoidal aggregate particles can be determined by nine characteristic parameters: reference point (O(XO, YO, ZO)), semimajor axis of ellipsoid (a, b, c), precession angle (γ), nutation angle (φ), and spin angle (ψ).


## 2.2. Spatial Random Distribution of Aggregate Particles

Reasonable aggregate size distribution can effectively improve the compactness and durability of concrete. Fuller grading curve is recognized by scholars as one of the most ideal grading curve, because the concrete formulated by Fuller grading has theoretically the maximum density and strength. Fuller curve can be illustrated by the following [14]:(6)pd=100ddmaxn,where p(d) is cumulative volume fraction of aggregate within the particle size d, dmax is the maximum aggregate size, and n is the index of the equation, n = 0.45∼0.70. For the graded section [ds, ds+1], the aggregate volume contained is(7)Vagds,ds+1=pds+1−pdspdmax−pdminRagVcon,where Rag is the volume ratio of aggregate (%); Vcon is the volume of concrete; dmax is the maximum size of aggregate; and dmin is the minimum size of aggregate. When the concrete area, the aggregate volume ratio, and the aggregate grading section [ds, ds+1] are determined, the aggregate volume of the graded section can be calculated by (7).


## 3. Interference Conditions between Aggregate Particles

### 3.1. Interference Judgment between Aggregate Particles

In practice, aggregates rarely touch or overlap each other, so an impact zone can be created around each aggregate to keep other aggregate particles out, as shown in Figure1. For round aggregate, the interference can be judged by the center distance of adjacent aggregate particles:(8)xOj−xOi2+yOj−yOi2≥ς3rj+ri,where ς3 is aggregate influence range coefficient [15, 16], which is 1.05 in this paper.Figure 1
Circular random aggregate models with same geometric parameters (rmin = 0.475, rmax = 4, Qs = 35.5%).Although many scholars have pointed out that the round aggregate could simulate the pebble aggregate, the actual aggregate is not a pure round, but a relatively smooth ellipsoid or sphere. Hence, this paper regards two-dimensional pebble aggregates as a mixture of round and ellipsoidal shapes (Figure2). When the lengths of the two axes of the ellipse are equal, the ellipse degenerates into a circle. Therefore, the circle can be regarded as a special form of the ellipse. Based on this idea, we can consider a given parameter ei, which is expressed as the ratio of the long axis to the short axis of the ellipse, to control the random change of the ellipse shape and ensure the appropriate length width ratio, namely:(9)ei=2ai2bi,where a and b are the long and short semiaxes of the ellipse, respectively.Figure 2
Elliptical and circular random mix aggregate with same geometric parameters (rmin = 0.475, rmax = 4, Qs = 35.5%).The method to judge the interference relationship between two ellipses is as follows [17]. An inscribed polygon is generated in one ellipse, and the geometric relationship of two ellipses is determined by judging the geometric relationship between the inscribed polygon and another ellipse. Now that the sum of the distances from any point of an ellipse to its two focal points is twice the length of the long semiaxis. If the length from one vertex of a polygon inscribed by an ellipse to the two focal points of another ellipse is less than or equal to a fixed value, the two ellipses will interfere; otherwise, the two ellipses will separate.Meanwhile, the problem can also be solved by using the quadratic form of the spherical equation. For any two ellipsesAXTAX = 0 and BXTBX = 0, the generalized characteristic polynomial of these two ellipses could be expressed as(10)fλ=detλA+B.If the characteristic polynomial has two distinct positive roots (λ), then the two ellipses are separated, else the two ellipses interfere with each other.In order to guarantee the convexity of polygons, various methods that can generate arbitrarily convex polygons have been proposed, but they still cannot fully satisfy the requirement of randomness or the algorithm is too complex. In order to solve the problem of randomness, the following methods are adopted in this paper: firstly, the range of polygon particle size is given according to gradation requirements of aggregate. Then the reference point and particle size of polygon are randomly generated. After generating a particle size each time, the relationship between the existing polygon and the newly generated polygon is compared. If (11) holds, the requirement is not satisfied; otherwise, the requirement is satisfied.(11)xi−x2+yi−y2<ri+rς4,where ς4 is adjustment coefficient, which is 1.05.The number of polygon sides and angles is randomly generated, and the angles are sorted from small to large anticlockwise. Finally, the coordinates of every vertex of the polygon corresponding to each angle can be obtained by(12)x=lricosang+X0y=lrisinang+Y0,where l is the elongation coefficient, reflecting the length width ratio of convex polygon. Since the actual gravel aggregate is neither narrowly pointed nor standardly round, its elongation should be controlled. The elongation coefficient l can be defined as follows:(13)l=lmaxlmin,where lmax and lmin are the maximum and minimum distances from the polygon vertex to the reference coordinate, respectively. Finally, the random convex polygon can be generated by connecting these vertices. The algorithm is easy to program and meets the requirements of randomness.The interference judgment methods of many scholars could be roughly divided into two types: area method [5] and angle method [11]. Most scholars use area as the measure to determine the interference relationship. For any convex polygon, set its vertex A1, A2, A3, … , Ai−1, Ai, Ai+1, … , An sorted counter clockwise, the coordinate of the vertex Ai is (xi, yi), and p(x, y) is a point in the plane; then the area of the triangle pAiAi+1 is Si, as shown in Figure 3. Assuming that the inner region enclosed by convex polygon A1A2A3 … An is ΩA and the boundary is ΩB, then p∈ΩA,Si>0i=0,1,2,⋯,np∈ΩB,atleastoneSi=0i=0,1,2,⋯,np∉ΩA,atleastoneSi<0i=0,1,2,⋯,n.Figure 3
The relationship betweena point p and a polygon. (a) Inside the polygon. (b) Outside the polygon. (c) Outside the polygon.
(a)(b)(c)The algorithm of area method is cumbersome, and it is easy to miss the special case like Figure3(c), so it is necessary to do edge check again. The angle method use included angles for interference judgment. Assuming that p is any point in the field Ω, the point p is connected with the vertices of the polygon A1A2A3 … An, and θi = ∠AipAi+1, θn = ∠AnpA1. Then,(14)p∈ΩA,∑i=1nθi=360∘p∉ΩA,∑i=1nθi=0,where θi can be obtained by cosine theorem, and the positive and negative of θi can be defined in the following way (taking ∠A1pA2 for example): set(15)W=x1−xpy2−yp−y1−ypx2−xp.IfW < 0, it is clockwise, taking θi as negative, otherwise, taking θi as positive.According to the above method, arbitrary convex polygon can be generated. However, the polygon generated by the above method may have some sharp corners, and the phenomenon can also be seen in other researches. The existence of sharp corners is not consistent with the actual aggregate shape, and the sharp corners are extremely detrimental to the finite element meshing and the convergence of the calculation. For the mechanical model, the sharp corners can also cause stress concentration, which affects the calculation results. Therefore, in order to reduce the probability of sharp corners, the following improvement methods are proposed:(16)angm+1−angm<α1n&angm+1−angm>α2n,where α1 and α2 are defined angles, α1 > α2. The difference between the two adjacent angles is limited to a certain range. Figure 4 shows the aggregate models at different angles α1 and α2, which proves the necessity of the control of the adjacent angles.Figure 4
2D random aggregate model of crushed stone concrete. (a)rmin = 0.475, rmax = 4, Qs = 30%, n = 4–10, α1 = 200°, α2 = 10°; (b) rmin = 0.475, rmax = 4, Qs = 30%, n = 4–10, α1 = 200°, α2 = 20°; (c) rmin = 0.475, rmax = 4, Qs = 30%, n = 4–10, α1 = 200°, α2 = 25°.
(a)(b)(c)The random aggregate of three-dimensional pebble can be expressed as spherical or ellipsoid. The interference relationship between two spherical aggregates can be judged by the center distance of aggregate particles:(17)XOj−XOi2+YOj−YOi2+ZOj−ZOi2≥ς3rj+ri,where ς3 is influence range coefficient of aggregate, which is 1.05 in this paper.Figure5 shows the spherical random aggregate model with particle size d∈(2.36, 40) mm and volume fraction of 40%.Figure 5
Geometric model of single-graded spherical aggregate.According to the grading range of coarse aggregate specified in industry recommended standards of the PRC (JTG/T 3650–2020) [18], standard concrete specimens are generated with the size of 150 mm × 150 mm × 150 mm. The corresponding grading parameters are shown in Table 1. The total volume fraction of coarse aggregate is 45%, as shown in Figure 6.Table 1
Gradation of coarse aggregates.
Nominal particle size (mm)26.5∼31.519∼26.516∼199.5∼164.75∼9.52.36∼4.75Volume content (%)5303020150Volume fractionQs (%)2.2513.513.596.750Figure 6
Geometric model of continuous-graded spherical aggregates.Although many scholars point out that the spherical aggregate could simulate the pebble aggregate, the actual aggregate is closer to the random mixing of sphere and ellipsoid. In polar coordinates, the standard equation of ellipsoid can be expressed as(18)x=asinθcosφy=bsinθsinφz=ccosθ,where a, b, and c are the three semimajor axes of the ellipse; θ, φ is the centrifugal angle of ellipsoid, respectively. When a = b = c, the elliptic equation degenerates into a spherical equation.The quadratic form corresponding to ellipsoidA is(19)XTAX=XT1a200001b200001c20000−1X=0,where A is ellipsoid coefficient matrix, X is homogeneous coordinates of any point on the ellipsoid, and X = (x, y, z, 1)T.The ellipsoid generated above is the standard ellipsoid with the reference point at the origin. In fact, the aggregate is arranged at any position according to any angle in the concrete, so it should be translated and rotated. After translating the ellipsoid to a random position (XO, YO, ZO), the coordinate of any point on the ellipsoid is as follows:(20)X′=BX=100XO010YO001ZO0001X.After rotation, the coordinate of any point on the ellipsoid is(21)X″=C1C2C2X′,or(22)X″=cosψsinψ00−sinψcosψ000010000110000cosϕsinϕ00−sinϕcosϕ00001cosγsinγ00−sinγcosγ0000100001X′,where ψ, ϕ, and γ are three independent random numbers between 0 and π/2, respectively. The algebraic expression of the ellipsoid after rotation is as follows:(23)X″TAX″=C1C2C3BXTAC1C2C3BX=XTC1C2C3BTAC1C2C3BX=0.SetD=C1C2C3BTAC1C2C3B; then, any elliptic equation in space is XTDX=0.In actual concrete, there is no intersection and inclusion between any two ellipsoids, so the minimum distance between them should be controlled. The positional relationship between the two ellipsoids can be judged by solving the characteristic equation of the two ellipsoids [19, 20]:(24)detλA1+A2=0,where A1 and A2 are the coefficient matrix (diagonal matrix) of any two standard ellipsoids, respectively and λ is the generalized eigenvalue. When there are two different positive answers to the generalized characteristic equation, the two ellipsoids are separated. The eigenvalue method can separate ellipsoids, but cannot give the minimum distance between ellipsoids. At the same time, the method [19] does not study the ellipsoidal interference judgment after rotation in space.Suppose that two ellipsoid equations areE1:XTD1X, E2:XTD2X. As shown in Figure 7(a) [21], xk, ykk∈N+ is a point on two ellipsoids, respectively; then, the shortest distance between the two ellipsoids is as follows:(25)dE1,E2=limk⟶∞xk−yk.Figure 7
Diagram for judging the interference relationship between two ellipsoids.
(a)(b)In thek-th iteration, we construct sphere Qc1,k,r1,k inscribed at E1,Qc2,k,r2,k inscribed at E2 with xk, yk as the internal tangent points. If Q(c1,kr2,k) is completely contained in E1∪E2, then the two ellipsoids are tangent or intersecting. Otherwise, the next iteration is carried out with xk+1, yk+1 as the internal tangent points. If ∀k∈N+, Q(c1,kr2,k)⊄E1∪E2, then the two ellipsoids are separated. Each iteration of this algorithm needs to generate an inscribed sphere and calculate the position relationship between Q(c1,kr2,k) and E1∪E2. The calculation process is complicated and the amount of calculation is large.Based on the above problems, this paper proposes a new method to determine the positional relationship between two ellipsoids, as shown in Figure7(b). Firstly, assume that the center of the ellipse E1 is O1. Then, select point p on the surface of the ellipsoid E2, and determine the line segment L by O1 and p points. The position relationship between ellipsoids E1 and E2 could be simplified as the relationship between point g on the line segment L and the ellipsoid E1. For any straight line segment L, if ∃g∉E1 is guaranteed, then point p is separated from E1. The specific algorithm steps are as follows:(1)
By determining the pointpi(xip, ypi, zpi) on the ellipsoid E2, the parameter equation of the straight line segment Li could be determined, so as to determine the point gi on the straight line segment Li.(2)
When gi is substituted into the equation of ellipsoidE1, if ∀gi∈E1 holds, then E1 and E2 will not be separated and the calculation will be terminated. Otherwise, the point pi + 1 will be determined to judge the position relationship between gi+1 and E1.(3)
If∀pi∈E2,∃gi∉E1, the two ellipsoids are separated.The thickness between cement mortar aggregates cannot be quantified due to the separation of the two ellipsoids. Therefore, the amplification factor matrixξ is introduced:(26)ξ=ηηη1,where η is the aggregate size adjustment factor, η = 1/1.052.After enlarging the ellipsoidE1, the virtual ellipsoid E1e is obtained, and its coefficient matrix is(27)D1′=C1C2C3BTξA1C1C2C3B.After ensuring that the virtual ellipsoidE1e and E2 are tangent or separated, the distance between the ellipsoid E1e and E1 is the minimum mortar thickness between the aggregates. Figure 8 shows the single-graded ellipsoidal random aggregate cube with a side length of 150 mm. Figure 9 shows the continuous-graded ellipsoidal random aggregate model, and the aggregate grading parameters are the same as Table 1.Figure 8
Single-graded ellipsoidal random aggregate,d ∈ (4.75, 31.5) mm, Qs = 38.4%.Figure 9
Continuous-graded ellipsoidal random aggregate.


### 3.2. Delivery of Aggregate

After determining the reference point and orientation of aggregate particles, the aggregate particles are placed into the concrete domain Ω. According to the spatial relationship between aggregate and cement mortar matrix in actual concrete, each aggregate particle should meet the following requirements:(1)
The aggregate particles must be all located in the domain Ω of concrete.(2)
The aggregate particles put in later cannot overlap with those already put in.(3)
The distance between the edge of aggregate particles and the edge of concrete area shall not be less than the minimum thickness of cement mortar coating.(4)
The minimum distance between aggregate particles shall be maintained.The input requirement (1) should meet the requirements of (1), and the input requirement (3) should meet the requirements of (28) for the two-dimensional aggregate model and (29) for the three-dimensional aggregate model. The release conditions (2) and (4) can be controlled by the interference judgment discussed above.(28)XO−ς1ri≥f11&XO+ς2ri≤f12YO−ς1ri≥f21&YO+ς2ri≤f22,(29)XO−ς1ri≥f11&XO+ς2ri≤f12,YO−ς1ri≥f21&YO+ς2ri≤f22,ZO−ς1ri≥f31&ZO+ς2ri≤f32,where ς1 and ς2 are the adjustment factors used to make the aggregate particles not tangent to the boundary of the domain Ω or pass through the outside of the domain or too close to the boundary to meet the mortar wrapping thickness. In this paper, ς1 = ς2 = 1.05.The spatial distribution of aggregate particles can be divided into two steps: “taking” and “releasing.” “Taking” is to take the aggregate particles that meet the (shape and size) requirements from the given grading curve, while “release” refers to placing the collected aggregate particles randomly in the concrete area according to certain rules to form a random combination sequence of aggregate particles. In the process of “releasing,” attention should be paid to the boundary relationship between the aggregate particle generated each time and the concrete area Ω, and the interference relationship between the aggregates generated successively. Therefore, the release process is more complex than the generation process. Generally speaking, “take” and “release” are two continuous processes, that is, generation ⟶ release ⟶ generation, which is conducive to the judgment and control of the generated quantity. Based on all the analysis above, the general process of concrete aggregate generation and delivery can be summarized as Figure10.Figure 10
The basic process of stochastic aggregate particles generation and delivery.


### 3.3. Interfacial Transition Zone

It is generally believed that the meso concrete is composed of mortar, coarse aggregate, and their ITZ, and the thickness of ITZ is generally 15∼120μm [22–26]. Because of the different properties between cement paste and aggregate, ITZ is the first area of concrete to damage and fracture. These defects not only have a great impact on the mechanical properties of concrete, but also have certain impact on the freezing resistance, permeability resistance, and corrosion resistance of concrete. Therefore, ITZ should be involved in the mesoscale model generally.The ITZs of mesoscopic stochastic concrete could be obtained by appropriately enlarging the aggregate. Another method is to endow ITZs with corresponding property around the aggregate when meshing [10], but the first method is more versatile. Now that the thickness of the interfacial transition zone around each aggregate is not uniform, the influence of randomness on the formation of interfacial transition zone should be considered. For convex polygonal and circular aggregates, a random number can be generated as the amplification factor, and the original aggregate can be amplified by the coordinate. For ellipsoidal or spherical aggregates, it can be achieved by generating a magnification factor matrix M:(30)M=ηaηbηc1,where ηk is a random number, and k = a, b, c, ηk∈(1.003, 1.032). In this paper, the calculation method of ηk is as follows:(31)ηk=m+hm,h=0.003∼0.032m,m=a,b,c.


## 3.1. Interference Judgment between Aggregate Particles

In practice, aggregates rarely touch or overlap each other, so an impact zone can be created around each aggregate to keep other aggregate particles out, as shown in Figure1. For round aggregate, the interference can be judged by the center distance of adjacent aggregate particles:(8)xOj−xOi2+yOj−yOi2≥ς3rj+ri,where ς3 is aggregate influence range coefficient [15, 16], which is 1.05 in this paper.Figure 1
Circular random aggregate models with same geometric parameters (rmin = 0.475, rmax = 4, Qs = 35.5%).Although many scholars have pointed out that the round aggregate could simulate the pebble aggregate, the actual aggregate is not a pure round, but a relatively smooth ellipsoid or sphere. Hence, this paper regards two-dimensional pebble aggregates as a mixture of round and ellipsoidal shapes (Figure2). When the lengths of the two axes of the ellipse are equal, the ellipse degenerates into a circle. Therefore, the circle can be regarded as a special form of the ellipse. Based on this idea, we can consider a given parameter ei, which is expressed as the ratio of the long axis to the short axis of the ellipse, to control the random change of the ellipse shape and ensure the appropriate length width ratio, namely:(9)ei=2ai2bi,where a and b are the long and short semiaxes of the ellipse, respectively.Figure 2
Elliptical and circular random mix aggregate with same geometric parameters (rmin = 0.475, rmax = 4, Qs = 35.5%).The method to judge the interference relationship between two ellipses is as follows [17]. An inscribed polygon is generated in one ellipse, and the geometric relationship of two ellipses is determined by judging the geometric relationship between the inscribed polygon and another ellipse. Now that the sum of the distances from any point of an ellipse to its two focal points is twice the length of the long semiaxis. If the length from one vertex of a polygon inscribed by an ellipse to the two focal points of another ellipse is less than or equal to a fixed value, the two ellipses will interfere; otherwise, the two ellipses will separate.Meanwhile, the problem can also be solved by using the quadratic form of the spherical equation. For any two ellipsesAXTAX = 0 and BXTBX = 0, the generalized characteristic polynomial of these two ellipses could be expressed as(10)fλ=detλA+B.If the characteristic polynomial has two distinct positive roots (λ), then the two ellipses are separated, else the two ellipses interfere with each other.In order to guarantee the convexity of polygons, various methods that can generate arbitrarily convex polygons have been proposed, but they still cannot fully satisfy the requirement of randomness or the algorithm is too complex. In order to solve the problem of randomness, the following methods are adopted in this paper: firstly, the range of polygon particle size is given according to gradation requirements of aggregate. Then the reference point and particle size of polygon are randomly generated. After generating a particle size each time, the relationship between the existing polygon and the newly generated polygon is compared. If (11) holds, the requirement is not satisfied; otherwise, the requirement is satisfied.(11)xi−x2+yi−y2<ri+rς4,where ς4 is adjustment coefficient, which is 1.05.The number of polygon sides and angles is randomly generated, and the angles are sorted from small to large anticlockwise. Finally, the coordinates of every vertex of the polygon corresponding to each angle can be obtained by(12)x=lricosang+X0y=lrisinang+Y0,where l is the elongation coefficient, reflecting the length width ratio of convex polygon. Since the actual gravel aggregate is neither narrowly pointed nor standardly round, its elongation should be controlled. The elongation coefficient l can be defined as follows:(13)l=lmaxlmin,where lmax and lmin are the maximum and minimum distances from the polygon vertex to the reference coordinate, respectively. Finally, the random convex polygon can be generated by connecting these vertices. The algorithm is easy to program and meets the requirements of randomness.The interference judgment methods of many scholars could be roughly divided into two types: area method [5] and angle method [11]. Most scholars use area as the measure to determine the interference relationship. For any convex polygon, set its vertex A1, A2, A3, … , Ai−1, Ai, Ai+1, … , An sorted counter clockwise, the coordinate of the vertex Ai is (xi, yi), and p(x, y) is a point in the plane; then the area of the triangle pAiAi+1 is Si, as shown in Figure 3. Assuming that the inner region enclosed by convex polygon A1A2A3 … An is ΩA and the boundary is ΩB, then p∈ΩA,Si>0i=0,1,2,⋯,np∈ΩB,atleastoneSi=0i=0,1,2,⋯,np∉ΩA,atleastoneSi<0i=0,1,2,⋯,n.Figure 3
The relationship betweena point p and a polygon. (a) Inside the polygon. (b) Outside the polygon. (c) Outside the polygon.
(a)(b)(c)The algorithm of area method is cumbersome, and it is easy to miss the special case like Figure3(c), so it is necessary to do edge check again. The angle method use included angles for interference judgment. Assuming that p is any point in the field Ω, the point p is connected with the vertices of the polygon A1A2A3 … An, and θi = ∠AipAi+1, θn = ∠AnpA1. Then,(14)p∈ΩA,∑i=1nθi=360∘p∉ΩA,∑i=1nθi=0,where θi can be obtained by cosine theorem, and the positive and negative of θi can be defined in the following way (taking ∠A1pA2 for example): set(15)W=x1−xpy2−yp−y1−ypx2−xp.IfW < 0, it is clockwise, taking θi as negative, otherwise, taking θi as positive.According to the above method, arbitrary convex polygon can be generated. However, the polygon generated by the above method may have some sharp corners, and the phenomenon can also be seen in other researches. The existence of sharp corners is not consistent with the actual aggregate shape, and the sharp corners are extremely detrimental to the finite element meshing and the convergence of the calculation. For the mechanical model, the sharp corners can also cause stress concentration, which affects the calculation results. Therefore, in order to reduce the probability of sharp corners, the following improvement methods are proposed:(16)angm+1−angm<α1n&angm+1−angm>α2n,where α1 and α2 are defined angles, α1 > α2. The difference between the two adjacent angles is limited to a certain range. Figure 4 shows the aggregate models at different angles α1 and α2, which proves the necessity of the control of the adjacent angles.Figure 4
2D random aggregate model of crushed stone concrete. (a)rmin = 0.475, rmax = 4, Qs = 30%, n = 4–10, α1 = 200°, α2 = 10°; (b) rmin = 0.475, rmax = 4, Qs = 30%, n = 4–10, α1 = 200°, α2 = 20°; (c) rmin = 0.475, rmax = 4, Qs = 30%, n = 4–10, α1 = 200°, α2 = 25°.
(a)(b)(c)The random aggregate of three-dimensional pebble can be expressed as spherical or ellipsoid. The interference relationship between two spherical aggregates can be judged by the center distance of aggregate particles:(17)XOj−XOi2+YOj−YOi2+ZOj−ZOi2≥ς3rj+ri,where ς3 is influence range coefficient of aggregate, which is 1.05 in this paper.Figure5 shows the spherical random aggregate model with particle size d∈(2.36, 40) mm and volume fraction of 40%.Figure 5
Geometric model of single-graded spherical aggregate.According to the grading range of coarse aggregate specified in industry recommended standards of the PRC (JTG/T 3650–2020) [18], standard concrete specimens are generated with the size of 150 mm × 150 mm × 150 mm. The corresponding grading parameters are shown in Table 1. The total volume fraction of coarse aggregate is 45%, as shown in Figure 6.Table 1
Gradation of coarse aggregates.
Nominal particle size (mm)26.5∼31.519∼26.516∼199.5∼164.75∼9.52.36∼4.75Volume content (%)5303020150Volume fractionQs (%)2.2513.513.596.750Figure 6
Geometric model of continuous-graded spherical aggregates.Although many scholars point out that the spherical aggregate could simulate the pebble aggregate, the actual aggregate is closer to the random mixing of sphere and ellipsoid. In polar coordinates, the standard equation of ellipsoid can be expressed as(18)x=asinθcosφy=bsinθsinφz=ccosθ,where a, b, and c are the three semimajor axes of the ellipse; θ, φ is the centrifugal angle of ellipsoid, respectively. When a = b = c, the elliptic equation degenerates into a spherical equation.The quadratic form corresponding to ellipsoidA is(19)XTAX=XT1a200001b200001c20000−1X=0,where A is ellipsoid coefficient matrix, X is homogeneous coordinates of any point on the ellipsoid, and X = (x, y, z, 1)T.The ellipsoid generated above is the standard ellipsoid with the reference point at the origin. In fact, the aggregate is arranged at any position according to any angle in the concrete, so it should be translated and rotated. After translating the ellipsoid to a random position (XO, YO, ZO), the coordinate of any point on the ellipsoid is as follows:(20)X′=BX=100XO010YO001ZO0001X.After rotation, the coordinate of any point on the ellipsoid is(21)X″=C1C2C2X′,or(22)X″=cosψsinψ00−sinψcosψ000010000110000cosϕsinϕ00−sinϕcosϕ00001cosγsinγ00−sinγcosγ0000100001X′,where ψ, ϕ, and γ are three independent random numbers between 0 and π/2, respectively. The algebraic expression of the ellipsoid after rotation is as follows:(23)X″TAX″=C1C2C3BXTAC1C2C3BX=XTC1C2C3BTAC1C2C3BX=0.SetD=C1C2C3BTAC1C2C3B; then, any elliptic equation in space is XTDX=0.In actual concrete, there is no intersection and inclusion between any two ellipsoids, so the minimum distance between them should be controlled. The positional relationship between the two ellipsoids can be judged by solving the characteristic equation of the two ellipsoids [19, 20]:(24)detλA1+A2=0,where A1 and A2 are the coefficient matrix (diagonal matrix) of any two standard ellipsoids, respectively and λ is the generalized eigenvalue. When there are two different positive answers to the generalized characteristic equation, the two ellipsoids are separated. The eigenvalue method can separate ellipsoids, but cannot give the minimum distance between ellipsoids. At the same time, the method [19] does not study the ellipsoidal interference judgment after rotation in space.Suppose that two ellipsoid equations areE1:XTD1X, E2:XTD2X. As shown in Figure 7(a) [21], xk, ykk∈N+ is a point on two ellipsoids, respectively; then, the shortest distance between the two ellipsoids is as follows:(25)dE1,E2=limk⟶∞xk−yk.Figure 7
Diagram for judging the interference relationship between two ellipsoids.
(a)(b)In thek-th iteration, we construct sphere Qc1,k,r1,k inscribed at E1,Qc2,k,r2,k inscribed at E2 with xk, yk as the internal tangent points. If Q(c1,kr2,k) is completely contained in E1∪E2, then the two ellipsoids are tangent or intersecting. Otherwise, the next iteration is carried out with xk+1, yk+1 as the internal tangent points. If ∀k∈N+, Q(c1,kr2,k)⊄E1∪E2, then the two ellipsoids are separated. Each iteration of this algorithm needs to generate an inscribed sphere and calculate the position relationship between Q(c1,kr2,k) and E1∪E2. The calculation process is complicated and the amount of calculation is large.Based on the above problems, this paper proposes a new method to determine the positional relationship between two ellipsoids, as shown in Figure7(b). Firstly, assume that the center of the ellipse E1 is O1. Then, select point p on the surface of the ellipsoid E2, and determine the line segment L by O1 and p points. The position relationship between ellipsoids E1 and E2 could be simplified as the relationship between point g on the line segment L and the ellipsoid E1. For any straight line segment L, if ∃g∉E1 is guaranteed, then point p is separated from E1. The specific algorithm steps are as follows:(1)
By determining the pointpi(xip, ypi, zpi) on the ellipsoid E2, the parameter equation of the straight line segment Li could be determined, so as to determine the point gi on the straight line segment Li.(2)
When gi is substituted into the equation of ellipsoidE1, if ∀gi∈E1 holds, then E1 and E2 will not be separated and the calculation will be terminated. Otherwise, the point pi + 1 will be determined to judge the position relationship between gi+1 and E1.(3)
If∀pi∈E2,∃gi∉E1, the two ellipsoids are separated.The thickness between cement mortar aggregates cannot be quantified due to the separation of the two ellipsoids. Therefore, the amplification factor matrixξ is introduced:(26)ξ=ηηη1,where η is the aggregate size adjustment factor, η = 1/1.052.After enlarging the ellipsoidE1, the virtual ellipsoid E1e is obtained, and its coefficient matrix is(27)D1′=C1C2C3BTξA1C1C2C3B.After ensuring that the virtual ellipsoidE1e and E2 are tangent or separated, the distance between the ellipsoid E1e and E1 is the minimum mortar thickness between the aggregates. Figure 8 shows the single-graded ellipsoidal random aggregate cube with a side length of 150 mm. Figure 9 shows the continuous-graded ellipsoidal random aggregate model, and the aggregate grading parameters are the same as Table 1.Figure 8
Single-graded ellipsoidal random aggregate,d ∈ (4.75, 31.5) mm, Qs = 38.4%.Figure 9
Continuous-graded ellipsoidal random aggregate.


## 3.2. Delivery of Aggregate

After determining the reference point and orientation of aggregate particles, the aggregate particles are placed into the concrete domain Ω. According to the spatial relationship between aggregate and cement mortar matrix in actual concrete, each aggregate particle should meet the following requirements:(1)
The aggregate particles must be all located in the domain Ω of concrete.(2)
The aggregate particles put in later cannot overlap with those already put in.(3)
The distance between the edge of aggregate particles and the edge of concrete area shall not be less than the minimum thickness of cement mortar coating.(4)
The minimum distance between aggregate particles shall be maintained.The input requirement (1) should meet the requirements of (1), and the input requirement (3) should meet the requirements of (28) for the two-dimensional aggregate model and (29) for the three-dimensional aggregate model. The release conditions (2) and (4) can be controlled by the interference judgment discussed above.(28)XO−ς1ri≥f11&XO+ς2ri≤f12YO−ς1ri≥f21&YO+ς2ri≤f22,(29)XO−ς1ri≥f11&XO+ς2ri≤f12,YO−ς1ri≥f21&YO+ς2ri≤f22,ZO−ς1ri≥f31&ZO+ς2ri≤f32,where ς1 and ς2 are the adjustment factors used to make the aggregate particles not tangent to the boundary of the domain Ω or pass through the outside of the domain or too close to the boundary to meet the mortar wrapping thickness. In this paper, ς1 = ς2 = 1.05.The spatial distribution of aggregate particles can be divided into two steps: “taking” and “releasing.” “Taking” is to take the aggregate particles that meet the (shape and size) requirements from the given grading curve, while “release” refers to placing the collected aggregate particles randomly in the concrete area according to certain rules to form a random combination sequence of aggregate particles. In the process of “releasing,” attention should be paid to the boundary relationship between the aggregate particle generated each time and the concrete area Ω, and the interference relationship between the aggregates generated successively. Therefore, the release process is more complex than the generation process. Generally speaking, “take” and “release” are two continuous processes, that is, generation ⟶ release ⟶ generation, which is conducive to the judgment and control of the generated quantity. Based on all the analysis above, the general process of concrete aggregate generation and delivery can be summarized as Figure10.Figure 10
The basic process of stochastic aggregate particles generation and delivery.


## 3.3. Interfacial Transition Zone

It is generally believed that the meso concrete is composed of mortar, coarse aggregate, and their ITZ, and the thickness of ITZ is generally 15∼120μm [22–26]. Because of the different properties between cement paste and aggregate, ITZ is the first area of concrete to damage and fracture. These defects not only have a great impact on the mechanical properties of concrete, but also have certain impact on the freezing resistance, permeability resistance, and corrosion resistance of concrete. Therefore, ITZ should be involved in the mesoscale model generally.The ITZs of mesoscopic stochastic concrete could be obtained by appropriately enlarging the aggregate. Another method is to endow ITZs with corresponding property around the aggregate when meshing [10], but the first method is more versatile. Now that the thickness of the interfacial transition zone around each aggregate is not uniform, the influence of randomness on the formation of interfacial transition zone should be considered. For convex polygonal and circular aggregates, a random number can be generated as the amplification factor, and the original aggregate can be amplified by the coordinate. For ellipsoidal or spherical aggregates, it can be achieved by generating a magnification factor matrix M:(30)M=ηaηbηc1,where ηk is a random number, and k = a, b, c, ηk∈(1.003, 1.032). In this paper, the calculation method of ηk is as follows:(31)ηk=m+hm,h=0.003∼0.032m,m=a,b,c.


## 4. Random Mesoscopic Concrete Aggregate Geometric Model

### 4.1. Graphic Format Conversion of Aggregate Model

It is very important that mesoscopic stochastic concrete structures can be opened by finite element software or other drawing software for other computational and analytical activities. In this paper, random concrete model is generated by MATLAB software. MATLAB has powerful numerical calculation and drawing functions, but its output graphics are not versatile and cannot be directly opened and used by finite element software. Autodesk Computer Aided Design (AutoCAD) has powerful interface function, which can save image files in more abundant formats, so it can connect with most of the finite element softwares such as ANSYS. Using AutoCAD to read the SCR script file generated by MATLAB to realize the transformation of graphic file format has strong universality and high efficiency [27–29].Surface objects in MATLAB can be generated by different functions, such as the surf() type function and the mesh() type function. The color of the face image can be determined by the CData matrix. Surf type function can fill every cell in the surface with color. The surface generated by mesh type function only has color on the line (color changes with height). Because of lack of corresponding color filling data in AutoCAD script file, the color information will be lost when SCR file is generated from MATLAB graphics file. This conversion process is as follows: fopen() type function in MATLAB is used to generate an empty SCR file, and then fprintf() type function is used to read the graphic data generated by MATLAB. Finally, AutoCAD is used to read SCR script files and draw graphics files. Figure11 shows the examples of the transformation.Figure 11
Format conversion method.


### 4.2. Finite Element Model of Stochastic Concrete

After the establishment of stochastic concrete model, it can be imported into the finite element software for mesh and related analysis. The commonly used meshing methods include mapping mesh method, Delaunay triangulation method, improved Delaunay triangulation method, and improved advancing front technique (AFT) method. Due to the complex shape of meso concrete model, especially the convex polygon aggregate model, mesh generation is prone to be malformed. In the finite element software COMSOL, the mesh quality can be improved by controlling the maximum and minimum element size, maximum element growth rate, curvature factor, resolution of narrow region, etc. The malformed elements in the mesh can be removed until the shape passed the checking. The bandwidth of the stiffness matrix can be reduced by using the method of node rearrangement, so that the finite element numerical simulation of complex mesostructure of concrete can be realized. The steps of meshing are as follows:(1)
The aggregate model is generated by MATLAB according to the algorithm proposed, the SCR script file is used to realize the data format conversion, and the converted data are imported into the finite element software to generate the geometric model. For the model with ITZs, the default repair tolerance should be changed to 10−6∼10−7.(2)
As shown in Figure12, the imported graphics include two kinds of surface (body): one is the surface (volume) of aggregate particles, mortar, and ITZs; the other is the whole area (including aggregate particles, mortar, and ITZs). Because the aggregate and mortar need to be meshed separately, the mortar and ITZ should be separated from the whole area.(3)
The material properties of mortar, ITZs, and aggregates are assigned, respectively, and the generated mortar, ITZ, and aggregates are divided into elements, and the nodes are rearranged to reduce the bandwidth of stiffness matrix in finite element calculation. Finally, the data of nodes and elements are saved.Figure 12
Mesh generation of random concrete models.


## 4.1. Graphic Format Conversion of Aggregate Model

It is very important that mesoscopic stochastic concrete structures can be opened by finite element software or other drawing software for other computational and analytical activities. In this paper, random concrete model is generated by MATLAB software. MATLAB has powerful numerical calculation and drawing functions, but its output graphics are not versatile and cannot be directly opened and used by finite element software. Autodesk Computer Aided Design (AutoCAD) has powerful interface function, which can save image files in more abundant formats, so it can connect with most of the finite element softwares such as ANSYS. Using AutoCAD to read the SCR script file generated by MATLAB to realize the transformation of graphic file format has strong universality and high efficiency [27–29].Surface objects in MATLAB can be generated by different functions, such as the surf() type function and the mesh() type function. The color of the face image can be determined by the CData matrix. Surf type function can fill every cell in the surface with color. The surface generated by mesh type function only has color on the line (color changes with height). Because of lack of corresponding color filling data in AutoCAD script file, the color information will be lost when SCR file is generated from MATLAB graphics file. This conversion process is as follows: fopen() type function in MATLAB is used to generate an empty SCR file, and then fprintf() type function is used to read the graphic data generated by MATLAB. Finally, AutoCAD is used to read SCR script files and draw graphics files. Figure11 shows the examples of the transformation.Figure 11
Format conversion method.


## 4.2. Finite Element Model of Stochastic Concrete

After the establishment of stochastic concrete model, it can be imported into the finite element software for mesh and related analysis. The commonly used meshing methods include mapping mesh method, Delaunay triangulation method, improved Delaunay triangulation method, and improved advancing front technique (AFT) method. Due to the complex shape of meso concrete model, especially the convex polygon aggregate model, mesh generation is prone to be malformed. In the finite element software COMSOL, the mesh quality can be improved by controlling the maximum and minimum element size, maximum element growth rate, curvature factor, resolution of narrow region, etc. The malformed elements in the mesh can be removed until the shape passed the checking. The bandwidth of the stiffness matrix can be reduced by using the method of node rearrangement, so that the finite element numerical simulation of complex mesostructure of concrete can be realized. The steps of meshing are as follows:(1)
The aggregate model is generated by MATLAB according to the algorithm proposed, the SCR script file is used to realize the data format conversion, and the converted data are imported into the finite element software to generate the geometric model. For the model with ITZs, the default repair tolerance should be changed to 10−6∼10−7.(2)
As shown in Figure12, the imported graphics include two kinds of surface (body): one is the surface (volume) of aggregate particles, mortar, and ITZs; the other is the whole area (including aggregate particles, mortar, and ITZs). Because the aggregate and mortar need to be meshed separately, the mortar and ITZ should be separated from the whole area.(3)
The material properties of mortar, ITZs, and aggregates are assigned, respectively, and the generated mortar, ITZ, and aggregates are divided into elements, and the nodes are rearranged to reduce the bandwidth of stiffness matrix in finite element calculation. Finally, the data of nodes and elements are saved.Figure 12
Mesh generation of random concrete models.


## 5. Application Examples

### 5.1. Meso Analysis of Chloride Ions Transport

In order to study the influence of aggregate shape on chloride ion transport, six random aggregate models were randomly generated with the volume fraction of 38.5% and size of 50 mm × 50 mm. The chloride transports under the same input conditions were calculated. The calculation conditions are as follows.(1) Chloride diffusion coefficient of cement mortar:Dc = 1.45e − 11 × (28/(28 + t/1[d]))0.35 m2/s, and chloride diffusion coefficient of aggregate: Da = 0. (2) Transport direction is one dimensional transmission, from left to right. (3) Boundary condition: Cs (x = 0, t) = 0.21%. (4) Initial condition: C (x, t = 0) = 0. (5) Diffusion time of chloride ions for computing is 8 a.The calculation results are shown in Figure13. It can be seen that there is a significant difference between the calculation results of the meso model and the homogeneous mortar model, and the calculation results of the different aggregate models also differ significantly. Chloride ions will flow around the aggregate in the process of propagation, so the streamline is not straight, but around the aggregate. Due to the impermeability of the aggregate, the aggregate has a significant blocking effect on the chloride ion transport, which complicates the direction of the chloride ion concentration gradient and the transport path. For models with the same aggregate shape and volume fraction, the chloride ion concentration distribution is similar.Figure 13
Chloride ions concentration distribution and streamline. (a) Mortar; (b) circle; (c) convex polygon; (d) mixed circle and ellipse; (e) mixed circle and ellipse; (f) mixed circle and ellipse.
(a)(b)(c)(d)(e)(f)To more clearly demonstrate the differences between these models, the distribution of chloride ion concentrations at the same depth for the models is given in Figure14. On the concrete surface, the chloride ion concentration increases due to the blocking effect of aggregates, and the calculation results of the fine view model are significantly higher than the macroscopic model. The concentration distribution of chloride ion has multiple peaks due to aggregation. As the concrete depth increases, the aggregation effect diminishes, and the results of mesoscopic model are close to those of the macroscopic model, but there are still many peaks and valleys. By comparing the calculated results of these three mesoscopic aggregate models, it was found that the convex polygonal aggregates showed the greatest variation, followed by the circular and elliptical mixed aggregates. The results show that the more complex the shape of aggregate, the greater the impact on chloride ions transport.Figure 14
Chloride concentration distribution of different models at the same depth. (a)x = 5 mm; (b) x = 10 mm; (c) x = 20 mm; (d) C-Qs curve.
(a)(b)(c)(d)To verify the effect of volume fraction of aggregates on chloride ion transport, six random aggregate models with different aggregate volume fractions were randomly generated and numerical experiments on chloride ion transport were conducted, as shown in Figure15. The relevant parameters of aggregates are shown in Table 2.Figure 15
Chloride ions concentration distribution and streamline. (a) Model 1. (b) Model 2. (c) Model 3. (d) Model 4. (e) Model 5. (f) Model 9.
(a)(b)(c)(d)(e)(f)Table 2
Parameters of random concrete aggregate models.
Modelrmin (mm)rmax (mm)Qs (%)Number of particlesJ (mol/m)·10−4Mortar00003.7659Aggregate0.4754100800Model 10.475419.1803.0037Model 20.475423.41502.7525Model 30.475428.52002.5188Model 40.475430.52502.4569Model 50.475434.43002.3449Model 60.475438.55002.1529Model 70.475438.55002.0977Model 80.475438.55002.1334Model 90.475441.15502.0628J is defined as the integral of chloride concentration (C) in space; that is(32)J=∫Cdxdy.With the increase of aggregate volume fraction, the total amount of material transport decreases. Figure14(d) is the curve of the relationship between aggregate volume fraction (Qs) and total chloride content (J), which shows a linear relationship. The curve is obtained by fitting the calculated results of the test pieces:(33)J=−0.0423Qs+3.7658,where Qs is the aggregate volume fraction (%). When Qs = 0, J = 3.7658, the calculated results of (33) are consistent with the calculated values of the homogeneous mortar model in Figure 13(a). Therefore, a linear model of total chloride ion and aggregate volume fraction at a certain time is presented as the following:(34)J=aQs+b,where a and b are constants at a certain time, which are only related to materials, and a < 0, b > 0.


### 5.2. Mesoanalysis of Mechanical Properties

In order to illustrate the influence of aggregate on the stress distribution of concrete structure, the concrete aggregate structure model in Section4.1 is still used for finite element analysis. For the aggregate, the elastic modulus of aggregate (E0) is 60 GPa, and Poisson’s ratio (v) is 0.16. For cement mortar, E0 = 20 GPa, and v = 0.2. It is assumed that the vertical deformation of the concrete cube is 10 μm. The calculated stress distribution is shown in Figure 16. The results show that the local tensile stress comes into being due to the influence of heterogeneity caused by aggregate. The distribution of compressive stress in concrete is also uneven, which shows obvious differences from macroscopic theory. The content of aggregate has great influence on the local tensile stress, the tensile stress shows an upward trend with the increase of volume fraction of aggregate.Figure 16
Stress distribution of concrete with different aggregate volume fraction. (a)σx of model 1 to model 3; (b) σy of model 1 to model 3; (c) σx of models 4, 5, and 9; (d) σy of models 4, 5, and 9; (e) σx of model 6 to model 8; (f) σy model 6 to model 8.
(a)(b)(c)(d)(e)(f)In order to quantitatively explain the influence of aggregate on the mechanical properties of concrete, the elastic modulus (E) and Poisson’s ratio (v) are taken as indicators for analysis, and the definitions are as follows:(35)σs¯=∫σsdxdy∫dxdyεs¯=∫εsdxdy∫dxdy,where σs¯ is the average stress in s direction, εs¯ is the average strain in s direction, and s = x, y. According to the generalized Hooke’s Law:(36)εx¯=1E¯σx¯−v¯σy¯+σz¯εy¯=1E¯σy¯−v¯σx¯+σz¯,where E¯ and v¯ can be regarded as macroscopic elastic modulus and macroscopic Poisson’s ratio, respectively.The calculated results shown in Table3 indicate that the macro elastic modulus of numerical concrete changes significantly with the change of aggregate. With the increase of aggregate volume fraction, the macro elastic modulus increases and macro Poisson’s ratio decreases. The elastic moduli of the same volume fraction are basically equal.Table 3
Macroscopic elastic modulus and Poisson’s ratio of stochastic mesoscopic concrete model.
Modelσx¯ (MPa)σy¯ (MPa)εx¯εy¯E¯ (GPa)v¯Mortar————200.2Aggregate————600.16Model 1−0.1355−4.72933.40E − 05−2.00E − 0423.64650.1702Model 2−0.1370−4.90253.36E − 05−2.00E − 0424.51250.1679Model 3−0.1429−5.13603.37E − 05−2.00E − 0425.68000.1685Model 4−0.1416−5.23413.35E − 05−2.00E − 0426.17050.1676Model 5−0.1550−5.40023.28E − 05−2.00E − 0427.00100.1639Model 6−0.1587−5.62133.28E − 05−2.00E − 0428.10650.1639Model 7−0.1543−5.63583.27E − 05−2.00E − 0428.17900.1635Model 8−0.1539−5.61953.25E − 05−2.00E − 0428.09750.1626Model 9−0.1645−5.77483.24E − 05−2.00E − 0428.87400.1622In order to study the macroscopic characteristics of concrete, it is necessary to analyze the random characteristics of aggregate size, shape, distribution density, and distribution form. In addition, according to the required physical and mechanical properties, a mesoscopic stochastic concrete model can be used instead of manual tests for mix proportion design to achieve refined and intelligent design.


## 5.1. Meso Analysis of Chloride Ions Transport

In order to study the influence of aggregate shape on chloride ion transport, six random aggregate models were randomly generated with the volume fraction of 38.5% and size of 50 mm × 50 mm. The chloride transports under the same input conditions were calculated. The calculation conditions are as follows.(1) Chloride diffusion coefficient of cement mortar:Dc = 1.45e − 11 × (28/(28 + t/1[d]))0.35 m2/s, and chloride diffusion coefficient of aggregate: Da = 0. (2) Transport direction is one dimensional transmission, from left to right. (3) Boundary condition: Cs (x = 0, t) = 0.21%. (4) Initial condition: C (x, t = 0) = 0. (5) Diffusion time of chloride ions for computing is 8 a.The calculation results are shown in Figure13. It can be seen that there is a significant difference between the calculation results of the meso model and the homogeneous mortar model, and the calculation results of the different aggregate models also differ significantly. Chloride ions will flow around the aggregate in the process of propagation, so the streamline is not straight, but around the aggregate. Due to the impermeability of the aggregate, the aggregate has a significant blocking effect on the chloride ion transport, which complicates the direction of the chloride ion concentration gradient and the transport path. For models with the same aggregate shape and volume fraction, the chloride ion concentration distribution is similar.Figure 13
Chloride ions concentration distribution and streamline. (a) Mortar; (b) circle; (c) convex polygon; (d) mixed circle and ellipse; (e) mixed circle and ellipse; (f) mixed circle and ellipse.
(a)(b)(c)(d)(e)(f)To more clearly demonstrate the differences between these models, the distribution of chloride ion concentrations at the same depth for the models is given in Figure14. On the concrete surface, the chloride ion concentration increases due to the blocking effect of aggregates, and the calculation results of the fine view model are significantly higher than the macroscopic model. The concentration distribution of chloride ion has multiple peaks due to aggregation. As the concrete depth increases, the aggregation effect diminishes, and the results of mesoscopic model are close to those of the macroscopic model, but there are still many peaks and valleys. By comparing the calculated results of these three mesoscopic aggregate models, it was found that the convex polygonal aggregates showed the greatest variation, followed by the circular and elliptical mixed aggregates. The results show that the more complex the shape of aggregate, the greater the impact on chloride ions transport.Figure 14
Chloride concentration distribution of different models at the same depth. (a)x = 5 mm; (b) x = 10 mm; (c) x = 20 mm; (d) C-Qs curve.
(a)(b)(c)(d)To verify the effect of volume fraction of aggregates on chloride ion transport, six random aggregate models with different aggregate volume fractions were randomly generated and numerical experiments on chloride ion transport were conducted, as shown in Figure15. The relevant parameters of aggregates are shown in Table 2.Figure 15
Chloride ions concentration distribution and streamline. (a) Model 1. (b) Model 2. (c) Model 3. (d) Model 4. (e) Model 5. (f) Model 9.
(a)(b)(c)(d)(e)(f)Table 2
Parameters of random concrete aggregate models.
Modelrmin (mm)rmax (mm)Qs (%)Number of particlesJ (mol/m)·10−4Mortar00003.7659Aggregate0.4754100800Model 10.475419.1803.0037Model 20.475423.41502.7525Model 30.475428.52002.5188Model 40.475430.52502.4569Model 50.475434.43002.3449Model 60.475438.55002.1529Model 70.475438.55002.0977Model 80.475438.55002.1334Model 90.475441.15502.0628J is defined as the integral of chloride concentration (C) in space; that is(32)J=∫Cdxdy.With the increase of aggregate volume fraction, the total amount of material transport decreases. Figure14(d) is the curve of the relationship between aggregate volume fraction (Qs) and total chloride content (J), which shows a linear relationship. The curve is obtained by fitting the calculated results of the test pieces:(33)J=−0.0423Qs+3.7658,where Qs is the aggregate volume fraction (%). When Qs = 0, J = 3.7658, the calculated results of (33) are consistent with the calculated values of the homogeneous mortar model in Figure 13(a). Therefore, a linear model of total chloride ion and aggregate volume fraction at a certain time is presented as the following:(34)J=aQs+b,where a and b are constants at a certain time, which are only related to materials, and a < 0, b > 0.


## 5.2. Mesoanalysis of Mechanical Properties

In order to illustrate the influence of aggregate on the stress distribution of concrete structure, the concrete aggregate structure model in Section4.1 is still used for finite element analysis. For the aggregate, the elastic modulus of aggregate (E0) is 60 GPa, and Poisson’s ratio (v) is 0.16. For cement mortar, E0 = 20 GPa, and v = 0.2. It is assumed that the vertical deformation of the concrete cube is 10 μm. The calculated stress distribution is shown in Figure 16. The results show that the local tensile stress comes into being due to the influence of heterogeneity caused by aggregate. The distribution of compressive stress in concrete is also uneven, which shows obvious differences from macroscopic theory. The content of aggregate has great influence on the local tensile stress, the tensile stress shows an upward trend with the increase of volume fraction of aggregate.Figure 16
Stress distribution of concrete with different aggregate volume fraction. (a)σx of model 1 to model 3; (b) σy of model 1 to model 3; (c) σx of models 4, 5, and 9; (d) σy of models 4, 5, and 9; (e) σx of model 6 to model 8; (f) σy model 6 to model 8.
(a)(b)(c)(d)(e)(f)In order to quantitatively explain the influence of aggregate on the mechanical properties of concrete, the elastic modulus (E) and Poisson’s ratio (v) are taken as indicators for analysis, and the definitions are as follows:(35)σs¯=∫σsdxdy∫dxdyεs¯=∫εsdxdy∫dxdy,where σs¯ is the average stress in s direction, εs¯ is the average strain in s direction, and s = x, y. According to the generalized Hooke’s Law:(36)εx¯=1E¯σx¯−v¯σy¯+σz¯εy¯=1E¯σy¯−v¯σx¯+σz¯,where E¯ and v¯ can be regarded as macroscopic elastic modulus and macroscopic Poisson’s ratio, respectively.The calculated results shown in Table3 indicate that the macro elastic modulus of numerical concrete changes significantly with the change of aggregate. With the increase of aggregate volume fraction, the macro elastic modulus increases and macro Poisson’s ratio decreases. The elastic moduli of the same volume fraction are basically equal.Table 3
Macroscopic elastic modulus and Poisson’s ratio of stochastic mesoscopic concrete model.
Modelσx¯ (MPa)σy¯ (MPa)εx¯εy¯E¯ (GPa)v¯Mortar————200.2Aggregate————600.16Model 1−0.1355−4.72933.40E − 05−2.00E − 0423.64650.1702Model 2−0.1370−4.90253.36E − 05−2.00E − 0424.51250.1679Model 3−0.1429−5.13603.37E − 05−2.00E − 0425.68000.1685Model 4−0.1416−5.23413.35E − 05−2.00E − 0426.17050.1676Model 5−0.1550−5.40023.28E − 05−2.00E − 0427.00100.1639Model 6−0.1587−5.62133.28E − 05−2.00E − 0428.10650.1639Model 7−0.1543−5.63583.27E − 05−2.00E − 0428.17900.1635Model 8−0.1539−5.61953.25E − 05−2.00E − 0428.09750.1626Model 9−0.1645−5.77483.24E − 05−2.00E − 0428.87400.1622In order to study the macroscopic characteristics of concrete, it is necessary to analyze the random characteristics of aggregate size, shape, distribution density, and distribution form. In addition, according to the required physical and mechanical properties, a mesoscopic stochastic concrete model can be used instead of manual tests for mix proportion design to achieve refined and intelligent design.


## 6. Conclusions

This study summarizes the construction of mesoscopic concrete models and proposes a systematic method for efficient construction of 2D and 3D concrete stochastic mesoscopic structures. The approach can take into account the requirements of randomness in shape, spatial distribution, and particle size distribution, and the level of aggregate content can be controlled independent of single, continuous, or interstitial grading. The significance and application of the stochastic aggregate model are illustrated by two examples, and the feasibility of the developed stochastic aggregation model is demonstrated. The conclusion is as follows:(1)
A more efficient and realistic generation mechanism is proposed to generate arbitrary convex polygon not based on circle or ellipse, which can effectively control the aspect ratio and sharp angle of polygon.(2)
A new algorithm has been proposed for interference judgment of ellipsoidal aggregates, which can reasonably control the thickness of cement mortar around the aggregate. By controlling the axes ratio of ellipsoidal aggregate particles, this method can also be used to simulate the pores in concrete.(3)
Aggregate increases the tortuosity and distortion of chloride ion transport path and reduces the global speed. This effect is more notable in convex polygonal aggregate than in elliptical aggregate and circular aggregate. The influence of aggregate makes the concentration distribution of chloride ion present multiple peaks in local tensile stress.(4)
Numerical concrete simulating results show that the macro elastic modulus and Poisson’s ratio change obviously with the aggregate fraction; the local tensile stress in concrete increases with the increase of aggregate volume fraction. There is a good linear relationship between the total chloride content transported in concrete and the volume fraction of aggregate.The mesoscale random concrete model is an approach to realize the well-directed and intelligent design of mix proportion according to the physical and mechanical properties required substituting for manual experiments.



---

*Source: 5918272-2022-08-08.xml*