Composite

A computationally efficient method for generating virtual periodic representative volume element (RVE), capable of handling arbitrary inclusion shapes, is developed. A universal collision/overlap detection and repair method is proposed, where each inclusion shape is represented as a union of n-Spheres (UnS). A constrained optimization problem is formulated and solved to remove inclusion overlaps; a closed-form solution is derived for calculating the degree of inclusions overlap and its gradient vector with respect to inclusion position. RVE generation is illustrated with circular, spherical, four non-circular and four non-spherical inclusion shapes. Computational efficiency is demonstrated using an elaborate RVE generation time study. The generated RVEs are evaluated using various statistical metrics; results confirm the random distribution of inclusions. Effective properties of RVEs, representing unidirectional composites, are determined using homogenization with various fibre cross-section shapes; obtained mechanical properties have shown transverse isotropy.


Introduction
Over the past few decades, composites presence is increasing significantly in various industries due to their superior properties, such as higher specific stiffness and strength.Evaluation of these properties (and damage behaviour) is a non-trivial task.They are controlled by multiple factors such as constituent material properties, the relative proportion of constituents, geometrical arrangement of inclusions and properties of the inclusion-matrix interface.Experimental and theoretical approaches exist for evaluating composite properties, but the former approach is prohibitive due to the enormous amount of human and capital resources involved in the process.In theoretical approaches, the material's representative volume element (RVE) is used to find effective properties.RVE is chosen to contain all the typical heterogeneities of the microstructure, and its constitutive behaviour is the same as that of the whole composite material.Among the different theoretical methods, finite element analysis based numerical homogenization has become very popular due to its capability in capturing microstructural morphology.
Large-scale manufacturing of composites leads to random placement of fibres instead of regular arrangement.It was observed that in the case of long fibre reinforced composites, RVE of regularly arranged inclusions can predict effective properties up to a reasonable accuracy, but it is not a good choice for damage predictions [1][2][3][4].Also, the possible.The jamming limit varies from RVE generation algorithm to algorithm and generating RVEs with volume fraction close to jamming limit would take significantly longer times [7].For example, in the case of RSA, generating RVE of unidirectional composites, containing fibres of circular cross-section, with volume fractions greater than 50% would take very long time as the jamming occurs at about 54.6% inclusion volume fraction.Similarly, with RSA algorithm, jamming occurs at about 38% inclusion volume fraction in the case of particle reinforced composites of spherical shaped inclusions [14][15][16].These limits are less than the inclusion volume fractions of high strength composites used in aerospace applications and theoretical limits(90% for UDCs and 74% for PRCs).
The lower jamming limit of RSA is due to poor space utilization, as the position of the accepted inclusions is fixed.In order to overcome this, one line of RVE generation models [12,17] are focused on perturbing accepted inclusions so that space is made to accommodate extra inclusions, thereby increasing the volume fraction limit.Recently, Wenlong Tian et al. [18] has coupled RSA with molecular dynamics simulations to generate RVE of mono-disperse particles (or spheres) with volume fractions greater than 50%.It is reported that this model [18] requires large computational times for volume fractions ≥50% or number of particles >100.Vaughan et al. [19], and Yang et al. [20] developed growth-based models where inclusions are added about an initial seed point (i.e., inclusion) at distances drawn from experimentally determined distributions.These models are able to reach fibre volume fractions up to about 65%.In another line of RVE generation models [21][22][23], instead of adding inclusions in series, all inclusions are arranged regularly upfront (for example, in rectangular or hexagonal patterns) then each of them is given random perturbation to achieve random distribution of inclusions.It is observed that it is challenging to obtain randomness at higher volume fractions using these methods [22].
In a more recent line of research [13,24,25], the required number of inclusions are randomly initialized, then inclusion overlaps are removed using an iterative procedure.Pathan et al. [13] used L-BFGS-B constrained optimization algorithm to minimize/eliminate the inclusion overlaps.Though it can generate RVEs with higher inclusion volume fractions, very high computational times are reported.Herráez et al. [25] considered inclusions as rigid bodies where overlaps are removed using repulsive force and momentum.This model, [25], generates RVEs of non-circular cross-sections with higher volume fractions in relatively shorter times.In this work, [25], RVE overlap is quantified using the area of inclusions intersection, which may not be easily determined for arbitrary inclusion shapes.Inclusion overlap checking, in general, is non-trivial for non-circular or non-spherical shapes, especially when they are arbitrarily aligned in the space [26,27].Among the various existing methods, Gilbert-Johnson-Keerthi (GJK) algorithm [28] is more efficient for checking overlap of two arbitrary shapes in arbitrary positions.However, it is iterative in nature and intended for convex shapes.It can be extended to concave shapes as well by representing as a combination of convex shapes with increased complexity.
Machine Learning (ML) models [29][30][31][32] were developed to predict the effective properties of the composite materials.ML models, in general, are data-intensive and require a variety of data for better learning.In the case of RVE, this variety in data can be obtained from a broad spectrum of inclusion volume fractions, different inclusion shapes and different degrees of inclusion randomness.So, it is of interest to have an RVE generation algorithm handling these variations while maintaining the computational speed.
In this work, the RVE generation method based on bounds constrained optimization is proposed which works for any 2 and 3 inclusion shapes.An optimization problem is formulated to eliminate the inclusions overlap that occurred due to random initialization.As explained, algorithms like GJK can be used to efficiently detect the overlap between arbitrary shapes but finding the closed form expression for overlap cost and especially its gradient is not easy.Hence, in this work, a new method is developed for this purpose where each inclusion shape is represented as a union of n-Spheres(UnS).Also, the same UnS representation of a shape is used to determine its periodic copies.In Section 2, the proposed RVE generation methodology is described.Then, the method is applied to generate RVEs of five different 2 and five different 3 inclusion shapes in Section 3. Finally, in Sections 4 and 5, statistical and micromechanical validation is performed on the generated RVEs.

RVE generation methodology
RVE generation is accomplished majorly in two steps.First, inclusions are initialized randomly within the bounds of RVE, which may lead to inclusion-overlapping.Second, eliminating overlaps by inclusion repositioning while keeping them in contact with RVE.Fig. 1 shows all the steps involved in generating RVE of arbitrary inclusion shapes and a sample RVE generated using this algorithm.

Initialization
The inclusions are randomly placed, one after the other allowing overlaps, within the RVE boundary until the desired volume fraction is reached.In the present work, inclusion shapes are categorized as circular, spherical, non-circular and non-spherical shapes and their dimensions can be constant or drawn from a distribution of choice.Inclusions location and alignment are drawn from uniform distribution where each inclusion takes a position in the RVE domain with equal probability.For example, Eq. ( 1) shows optimization variables and their initialization formulae for rectangular or cubical shaped RVE.The set of optimization variables applicable for an inclusion depends on its geometry.Circular and spherical shapes need (, ) and (, , ) respectively, as their orientations are immaterial so that the remaining terms can be set to zero or some arbitrary constant.Similarly, noncircular and non-spherical shapes need (, , ) and (, , , , ) respectively.
where , , , ,  ∈ [0, 1] are drawn from uniform distribution,  is azimuthal angle,  is polar angle and (  ,   ,   ,   ,   ,   ) are bounds of the RVE in (, , ) directions.The choice of , in Eq. ( 1), is to ensure uniform distribution of points on the unit sphere, otherwise choosing  =  would cluster points near the poles [33].Random initialization, as per Eq. ( 1), may lead to the overlapping of inclusions, which must be removed.For this purpose, an optimizationbased iterative procedure for inclusions repositioning is presented in the next section.

Overlap elimination
In this section, all the steps involved in inclusion overlap elimination are presented in detail.Let  and  be any two inclusions with centres at   and   in a RVE domain  as shown in Fig. 2. The domain  is chosen as periodic to avoid gaps and overlaps when it is repeated in the space.Overlap between a pair of inclusions can be quantified either by using the intersection area or penetration distance along the line connecting centres.Due to its simplicity, we chose the distancebased approach.Let   and   , respectively, denote centre-to-centre Euclidean distance and distance of the closest approach for  and .Here, the distance of the closest approach is defined as the centreto-centre Euclidean distance when inclusions are just touching each other externally.Then, the total magnitude of inclusion overlap,  , in the RVE is defined as the cumulative sum of the overlap between all inclusion pairs as shown in Eq. ( 2) where () is Heaviside step function with unit value for all  > 0.
An optimization problem is formulated to minimize the total overlap,  , subject to the constraints that centre points,   , of inclusions does not cross the RVE domain .
In Eq. (3), the RVE domain () can be in any shape of interest and it decides the nature of the constraint equations.In the present study, the RVE domain  is chosen as rectangular in 2 and cuboid in 3 as it ensures simple bounds on the spatial variables (, , ).Non-monotone spectral projected gradient (NMSPG) is used, in the present work, to solve the Eq. ( 3) as explained in Section 3 while any constrained optimization solver can be used.

Shape representation as a union of n-Spheres (UnS)
The overlap evaluation is a non-trivial task for non-circular and nonspherical shapes.In this section, we describe a novel, time-efficient/ scalable approach to identify overlaps, thereby allowing our method to handle non-circular and non-spherical inclusion shapes.The idea is, overlap detection between a pair of arbitrary shapes can be converted to overlap detection between two groups of n-Spheres by representing inclusion shape as a union of n-Spheres (UnS).Here, the term ''n-Sphere'' implies a circle in 2D and a sphere in 3D, as used in mathematics.Fig. 3 shows the UnS representation of inclusion shapes considered in the present work.For example, checking the overlap between an ellipse and a rectangle requires finding the overlap between circles representing these two shapes.Overlap computation time is directly related to the number of n-Spheres used for each shape.Hence, using the minimum possible number of n-Spheres is desirable.
We proceed to find the centres and radii of n-Spheres to represent a shape.Let C  and C  represent the equations of  ℎ shape boundary and a concentric outer curve around C  with a thickness of   , minimum distance between boundaries of any two shapes.Here, C  represents the region within which the boundary of the other shapes is barred from entering.Assuming that we know the position  0 and radii  0 of the starting n-Sphere, then   and   of the successive n-Spheres are chosen such that the concentric n-Spheres of radius   +   would intersect on C  .This would ensure that the successive n-Spheres are separated by an optimum distance  *  .Overlap of inclusions may not be captured if n-Spheres are placed at a distance greater than  *  because such a UnS may not represent the shape accurately.Also, n-Spheres closer than  *  would unnecessarily increase the number of n-Spheres there by increasing the computation time.In Appendix A, formulae for finding the n-Spheres centre and radius are given for representing the ellipse/spheroid shape, the region between a pair of lines/planes separated by an angle 2.

Periodic copies addition
In a micromechanical analysis of the RVE, applying periodic boundary conditions (PBC) is more advantageous over homogeneous boundary conditions [34,35].In order to apply PBC, RVE must be continuous across its boundary when it is repeated in the space.To achieve this, if some part of an inclusion leaves the RVE boundary, then an appropriate number of its copies must be added to RVE at appropriate locations.Here, the number of periodic copies of an inclusion equals the number of RVE faces &∕ edges &∕ a vertex intersecting the inclusion.These periodic copies must be placed on the corresponding opposite face &∕ edge &∕ a vertex.For example, in 2D RVE, if an inclusion intersects two edges and a vertex, then three copies are added on two opposite edges and the opposite vertex.Finding intersections of RVE faces/edges/vertices with circular (in 2D) and spherical (in 3D) inclusion shapes is relatively simple.So, in this work, periodic copies are determined using the inclusions Union of n-Spheres (explained in the previous section) representation.
For each inclusion, the intersection of each of its n-Sphere with the RVE boundary is evaluated.Instances of inclusion's intersection with RVE is the same as the set of all n-Sphere's intersections with RVE faces &∕ edges &∕ a vertex.Accordingly, periodic copies of inclusions are added to RVE.As it involves running over each n-Sphere of each inclusion, filtering out boundary inclusions saves the computational effort.This can be done by checking the overlap of inclusion's bounding box with that of the RVE.

Cost function and gradient evaluation
Let inclusions  and , as shown in Fig. 2, have   and   number of n-Spheres in their UnS form.The case of circular/spherical inclusion shape can be obtained by choosing   = 1 and   = 1.The cost of overlap between  ℎ n-Sphere of  ℎ inclusion and  ℎ n-Sphere of  ℎ inclusion is defined as with, where,  is the Heaviside step function,   is the minimum surface-tosurface distance between any two inclusions,   and   are centres of  ℎ and  ℎ n-Spheres of  ℎ and  ℎ inclusions.The total cost of overlap,  , for all  inclusions is given by Gradient, the rate of change of the cost function,  , due to change in  ℎ inclusion variable   ∈ {  ,   ,   ,   ,   } is evaluated as (refer to Appendix B for detailed derivation) where, The cost function and its gradient vector, as given by Eqs. ( 6) and (7), can be used with any inclusion geometry by considering the appropriate set of optimization variables.For example, in the case of circular shapes (, , ) are constant hence gradients vanish corresponding to these variables.

Generating RVE of different inclusion shapes
This section demonstrates RVE generation for circular, spherical, four non-circular and four non-spherical inclusion shapes.Fig. 3 shows the UnS form of non-circular and non-spherical inclusion shapes.Generating an RVE requires solving the optimization problem given in Eq. (3).Constraints of Eq. ( 3) become simple bounds on the spatial variables for box shaped RVE (i.e., rectangle in 2D and cuboid in 3D).Gradient projection methods, an extension of the steepest descent method to the constrained optimization, are proved to be effective choices [36] for solving optimization problems with such simple bounds.Hence, Non-monotone spectral projected gradient (NMSPG) algorithm [37,38], a variant of gradient projection methods, is used to solve the Eq. ( 3).The particular choice of the algorithm is based on two reasons.One, the speed of convergence due to the choice of spectral step length.Two, guaranteed global convergence due to the non-monotone nature of the line-search.In NMSPG implementation, parameters  = 50,   = 10 −6 ,   = 10 6 ,  1 = 0.1,  2 = 0.9,  = 10 −4 are used with the same naming convention as in [37].RVE generation procedure is designed to restart if the solution does not converge in a predefined number of iterations.Table 1, lists three different sets of RVE parameters used in the present work.In this table,   is the inclusion volume fraction, (  ,   ) are mean and standard deviation of reference circle/sphere radii,  = ∕  is the RVE size, L is RVE side length and   is inclusions minimum surface-to-surface distance.Dimensions of different inclusion shapes, see Fig. C.11, are chosen to give the same area/volume as that of the reference circle/sphere so that their RVE generation times can be compared.Although RVEs can be generated for variable inclusion sizes, we have used   = 0 to study RVE generation for monodisperse inclusions.This is because the monodisperse case requires more computational effort and have lower jamming limits than those with polydisperse inclusions [13].Two and three-dimensional RVEs generated using the set 1 and set 3, respectively, are shown in Fig. 4 3.1.RVE generation time.
The RVE generation module is written in the open-source computational language Julia [39] and used on a computer with 1.19 GHz, four-core, 8 GB RAM specifications.Lower RVE generation times are desirable for producing large data sets which may be used for studying composites with data-hungry approaches like deep learning.To assess the computational performance, RVE generation time study (  ) is performed with various inclusion shapes having different aspect ratios and RVE sizes ().For this purpose, RVE parameters of sets 2 and 3 of Table 1, are used respectively, with  ∈ {25, 50, 75, 100} in 2 and  ∈ {15, 20, 30, 35} in 3.It is observed that, a large part of the   is due to overlap cost function ( ) and its gradient (∇ ) evaluations, and it gets enhanced with number of n-Spheres used in UnS form of inclusion(see Section 2.2.1).The  and ∇ are evaluated between a pair of inclusions only if their bounding boxes overlap to avoid unnecessary overlap check with far-off inclusions.This reduces computation costs significantly while dealing with non-convex shapes.
The minimum (  ), mean (  ), standard deviation (  ) and maximum (  ) values of   are reported for 2D and 3D RVEs of different inclusion shapes in Tables C. 5 and C.6, respectively, obtained from 20 realizations in each case of 3 RVEs and 50 realizations in each case of 2 RVEs.Of all four indicators, meantime is more reliable, as the other metrics are sensitive to outliers.If the random initialization is far from the solution, convergence may not be achieved in the predefined number of iterations, so the RVE generation restarts with new initialization.This leads to higher   and   in some cases.Thus, the parameters   ,   ,  max can be used to understand the uncertainty involved in that particular case.It is found that   is directly proportional to the RVE size and aspect ratio of the inclusion because the increase in RVE size leads to more number of inclusions, and an increase in aspect ratio leads to more chances of inclusion overlap.For n-tip stars,   is decreasing from a 3-tip star to a 5-tip star, as the 5-tip star has fewer chances of entanglement/overlap with other inclusions.This is because more tips would result in a smaller tip height for a given star area.
For circle-shaped inclusions, the present algorithm generated RVE of size  = 50 with 65% inclusion volume fraction in about 0.2-0.3s, while [13] has reported 107.02 min on a computer with similar specifications.This drastic reduction   may be majorly due to explicit gradient evaluation, which otherwise has to be calculated using computationally expensive numerical methods.Recently, [25] has reported RVE generation times for non-circular inclusion shapes.In comparison with this study, for generating the RVE of approximately 2070 elliptical and rectangular inclusions, the present algorithm takes 35 and 20 s, respectively, while [25] has reported 69 and 86 s.Higher   of 3D RVEs, compared to 2D RVEs, is due to an increased number of inclusions for a given RVE size and volume fraction.Note that, we have not reported   for RVEs of cylinders with  ∈ {30, 35} as the total number of spheres in the system (representing all cylinders in RVE) has reached 200,000 to 300,000, thus leading to very high   .RVE generation of such shapes can be made faster by developing cost-effective UnS representation in future work.

Statistical validation
It is important to assess the distribution of inclusions in the RVE, as it influences the mechanical response and damage initiation.In the following sections, the spatial distribution of inclusions is evaluated in the neighbourhood and at several radial distances about each inclusion.For this purpose, 20 different realizations of RVEs are generated for each inclusion shape, using the parameters given in set 2 (for 2D) and set 3 (for 3D) of Table 1.

Voronoi regions areas and volumes
The RVE domain is discretized using Voronoi tessellation [40], with inclusions centres as seed points, to demarcate a unique region (hence neighbours) for each inclusion.Voronoi tessellation of a RVE with circular inclusions of 65% volume fraction is shown in Fig. 5(a).The regular arrangement of inclusions generates Voronoi regions with the uniform area (or volume) and equidistant neighbours.So, these metrics are used to determine the randomness in the immediate neighbourhood of the inclusion.Coefficient of variation, mean normalized standard deviation(  = ∕), of Voronoi region areas (  , for 2 RVE) and Voronoi region volumes (  , for 3 RVE) is shown with box plots in Fig. 5.For regular arrangement,   = 0 due to vanishing standard deviation of Voronoi region areas or volumes and higher   imply more randomness.Table 2 shows good agreement of   for RVE of circular inclusions when compared with that reported in literature.For the same set of RVE parameters,   and   of all inclusions, except stars, are in the range of 0.1 to 0.15, see Fig. 5, in line with the observations of [22,25].It is observed that   is increasing with aspect ratio and decreasing with the number of tips in the case of the n-tip star.As explained previously, a 3-tip star has more chances of entanglement than a 5-tip star.

Nearest neighbour distributions
In this section, distributions of the first two nearest neighbour distances and nearest neighbour orientation are studied.Probability density functions(PDF) of distances from every inclusion to its closest and second-closest inclusion are plotted, for all the considered inclusion shapes, in Fig. 6.PDFs are determined using kernel density estimation with Gaussian kernel and Scott method based bandwidth.PDF becomes a sharp peak at the neighbour distance if all neighbours are equidistant, as in the case of a regular arrangement.PDF gets flat as the neighbour distance distribution changes from regular to random.Fig. 6 shows such PDF curves flattening with increasing aspect ratio and decreasing number of tips on the star.This further supports previous observations with   .
Cumulative distribution function (CDF) of nearest neighbour orientations is evaluated, for different inclusion shapes, using kernel density estimation and plotted in Fig. 7.When inclusions are in complete state of randomness (CSR), nearest neighbour angles are expected to take all possible values with equal probability.These angles are evaluated, using Eqs.(1d) and (1e), and plotted on the same Fig. 7.It is observed that, in all three plots of Fig. 7, CDF of different inclusion shapes is closely following that of the CSR.This indicates retained randomness in the nearest neighbour orientations, after the inclusions re-positioning during overlap removal process.

Ripley's K function, K(r)
Ripley's K function, (), determines the expected number of points inside a search circle (or sphere) of radius , centred at any point , using the Eq. ( 8) where,  is point density,  is the number of points and  is the indicator function with value of 1 if the point  is inside the search circle (or sphere), otherwise zero.Edge correction term,   (), factors the absence of points in the exterior part of the search circle (or sphere) if it crosses the bounds of the domain.The edge correction term can be avoided by replicating RVEs around the RVE due to the virtue of its periodicity.Residual Ripley's K function, (), given in Eq. ( 9) measures the deviation of a given distribution from that of the CSR.In the case of CSR,   () have  2 and 4 3  3 for 2 and 3 RVEs respectively.
In a search circle (or sphere) of a given radius , the number of expected points higher than   () indicates clustering, while the lower number implies dispersion.Hence, () takes the value of zero for CSR distribution and oscillates about zero for a regular arrangement.In Figs.8(a) and (b), () is plotted for 2 and 3 cases, with standard deviation error bars, evaluated from 20 realizations.In the case of 2 inclusion shapes, () is stabilized at about a radial distance ratio of 7 and remained constant but slightly above zero, indicating minor clustering for longer distances.In the case of 3 inclusion shapes, L(r) is not fluctuating but monotonically reduced to zero, indicating the absence of regular arrangement.

Radial distribution function, G(r)
Ripley's K function evaluates point patterns up to radius  in a cumulative manner.Hence, it becomes difficult to identify the radial distance where the given distribution is deviating from the CSR.For this purpose, the radial distribution function, (), given in Eq. (10), is derived from ().It evaluates the probability of finding points in a circular strip (or spherical shell) of inner radius  and small thickness , centred at any arbitrary point.
where,   is the number of points in an annulus of inner radius  and outer radius  +  centred at ℎ inclusion centre.In the case of CSR, () approaches unity at sufficiently longer distances from the inclusion centre.Radial distribution functions for 2 and 3 inclusion shapes, evaluated from 20 realizations, is shown in Figs.8(c) and (d) with standard deviation as error bars.In the case of both 2 and 3 inclusion shapes, () approaches unity at radial distance ratios of about 7, indicating the randomness at sufficiently longer radial distances.

Micro-mechanical validation
In this section, 2D RVEs generated using the proposed algorithm are assessed for transverse isotropy, using finite element based micromechanical analysis.Isotropy is expected on the cross-section normal to fibre direction due to the random distribution of fibres.In order to compare with the literature results, the same RVE parameters and material set E-glass/MY750/HY917/DY063 are considered [41].Four different inclusion shapes, circle, ellipse, capsule and rectangle, are considered in this analysis.The material set contains linear elastic isotropic matrix and fibre with Young's modulus and Poisson's ratio   = 3.35 GPa,   = 0.35,   = 74 GPa,   = 0.2 respectively.RVEs are generated with constant fibre radius   = 2.6 μm, RVE window size  = 50.0,  = 0.07  and inclusion volume fraction  ∈ {0.2, 0.3, 0.4, 0.5, 0.6}.For each  and inclusion shape, 20 different realizations are generated.
Finite element models of RVE are generated in ABAQUS using python scripts.As we are interested only in transverse properties of the composite, four-node plane strain elements (CPE4), along with a small proportion of three-node triangular elements (CPE3), are used in finite element modelling.The number of elements in each model is chosen as approximately 50,000 elements following a mesh convergence study.Periodic boundary conditions (PBC) are chosen over uniform boundary conditions.This is because, when PBC are applied, equivalent properties of composite converge faster with RVE size [42][43][44].PBC are applied on the four edges using dummy nodes (known as reference points in ABAQUS terminology) using the procedure explained in [45].
The effective material properties evaluated from the micromechanical analysis of RVEs, containing 60% volume fraction of circular inclusions, match closely with the values reported in the literature (see Table 3).Also, as shown in Table 4, transverse isotropy is observed in the RVEs of different fibre shapes.In Fig. 9, effective transverse elastic     Hence, RVEs generated using the proposed algorithm are suitable for modelling the microstructure of the unidirectional composite materials.

Conclusions
A computationally efficient method is developed for generating periodic RVEs of arbitrary inclusion shapes.Inclusion overlaps that occurred during random initialization are removed by solving a constrained optimization problem.Overlap detection between arbitrary inclusion shapes is accomplished by representing each of them as a union of n-Spheres(UnS).The following conclusions are drawn from a detailed RVE generation time study, statistical and micromechanical analysis of RVEs.In these studies, ten different inclusion shapes are used and the influence of varying aspect ratios is also considered.
• Repositioning of the inclusions using the proposed overlap cost function gradient had enabled faster RVE generation; gradient yields a coordinated movement of inclusions, by considering the position of their neighbours, to reduce the overall inclusion overlap in the RVE.• Increasing aspect ratio of the inclusions has increased computational time, lowered the maximum reachable volume fraction due to the increased chances of overlap.• Statistical analysis of RVEs has shown randomness, both, in the neighbourhood and at longer distances from the inclusion centre.• Micromechanical analysis of RVEs has shown transverse isotropy for circular, elliptical and rectangular inclusion shapes, indicating the random distribution of inclusion shapes.
Inclusions overlap detection using UnS representation can make the existing RVE generation algorithms handle arbitrary inclusion shapes.Also, representation could be used for object/shape collision tection in other fields like robotics and computer graphics.Due to computational efficiency and capability to work with different shapes, large and rich/varied data sets of RVEs can be generated for data-driven studies of composite materials.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix B. Gradient evaluation
The gradient of cost function,  , is evaluated, with respect to an independent variable   of the  ℎ inclusion, as follows For the overlap cost,  , between  ℎ and  ℎ n-Sphere of  ℎ and  ℎ inclusion, •  , =  , , as the cost overlap between  and  should be the same as that between  and .
•  , = 0 for  =  = , as the cost of overlap is defined for a pair of two different inclusions

Fig. 1 .
Fig. 1.RVE generation flow chart and a sample RVE of different inclusion shapes.

Fig. 3 .
Fig. 3. Non-circular and non-spherical shapes representation as a union of n-Spheres.

Fig. 5 .
Fig. 5.A typical Voronoi tessellation and its area and volume metrics for different inclusions.

Fig. 6 .
Fig. 6.Probability density function (PDF) of first two nearest neighbour distances of various inclusion shapes.

R
.Nakka et al.

Table 1
RVE parameter sets used in the present study.

Table 2
Coefficient of variation of area,   for RVE of circular inclusions.

Table 3
Transverse effective properties for RVE of circles.Transverse elastic and shear moduli for different inclusion shapes.

Table 4
Transverse isotropy for different inclusion shapes, with 60% volume fraction.
In the first term, replacing  with  and using  , =  , gives, = 0 for  ≠  and  ≠ , as the overlap magnitude between a pair of inclusions is independent of other inclusions position using the above properties of  , , Eq. (B.1), is simplified as,

Table C .
5 RVE generation times (in seconds) for different 2D inclusion shapes, with different RVE size  and number of inclusions   .RVE generation times (in seconds) for 3D inclusion shapes, with different RVE size  and number of inclusions   .